
\section{DFT calculation details} 

We carried out an {\it ab initio} band structure calculations within the frameworks of density functional theory (DFT) as implemented in VASP~\cite{kresse_efficient_1996, kresse_efficiency_1996, kresse_ab_1994} and a full-potential local-orbital code (FPLO) \cite{eschrig2004relativistic}. The generalized gradient approximation (GGA)~\cite{perdew_generalized_1996} for the exchange-correlation functional was used. 


{\it VASP.} The projector-augmented wave (PAW) technique was applied~\cite{blochl_projector_1994}. Strong electronic correlations and spin-orbit coupling (SOC) were taken into account via so-called GGA+U+SOC approximation (LDAUTYPE = 1)~\cite{Liechtenstein1995}. The on-site Coulomb repulsion parameter $U$ for Co ions was varied in a range from $5$ to $7$ eV~\cite{Zvereva2016}. The Hund's intra-atomic exchange parameter was taken to be $J_H$ = 0.9 eV. The electronic population numbers for Co were obtained by integration within atomic sphere with radius of $1.302$~\AA. The cutoff energy for plane waves was chosen to be $500$ eV. The integration over the Brillouin zone was performed using the Monkhorst-Pack scheme~\cite{monkhorst_special_1976}. The $7\times 7\times 7$ $k$-mesh was used for conventional DFT calculations (including the Wannier function projection  ~\cite{Mostofi2014} and construction of projected localized orbital - PLO~\cite{Schuler2018}), the $5\times 2\times 3$ $k$-mesh was used for the optimization of the crystal structure, while $3\times 3\times 3$ $k$-grid for calculation of the exchange constants. In two last cases supercells consisting of 4 and 5 formula units (f.u.) were used. 

{\it FPLO.} Hopping integrals were estimated starting from scalar relativistic calculations at the GGA (PBE) level with 12 $\times$ 12 $\times$ 12 $k$-grid. Wannier orbitals were constructed via projection \cite{koepernik2021symmetry}, with the local cubic projection axes chosen such that the $x+y+z$ direction was parallel to the crystallographic $c$-axis, and the angle between each local axis and the corresponding Co-O bond vector was minimized.  


{\it Wien2k.} To analyze formation of the bonding and antibonding $e_g$ states due to large hopping between third nearest neighbors, we used the methods described in Ref.~\onlinecite{Foyevtsova2013} based on the eigenstates from the full-potential linearized augmented plane-wave (LAPW) calculations as realized in Wien2k code~\cite{Wien2k,Blaha2020}. We chose the basis-size controlling parameter RK$_{\rm max}$ = 7 and a mesh of 500 {\bf k} points in the first Brillouin zone (FBZ) of the primitive unit cell. The density of states (DOS) was computed with 1000 {\bf k} points in the full Brillouin zone. 

The crystal structure of BaCo$_2$(AsO$_4$)$_2$ was taken from Ref.~\onlinecite{Dordevic2008}. In order to estimate the exchange interaction parameters we consider the anisotropic-exchange model
\begin{eqnarray}
H=\sum_{i \ne j}\vec S_i \mathbb{J}_{ij}  \vec S_j
\end{eqnarray}
where the summation runs twice over each pair of $S=1/2$, and $\mathbb{J}_{ij}$ is a $3\times 3$ matrix. We used the total energy method, described in Ref.~\onlinecite{Xiang2011} in a supercell with 10 Co ions ordered in a honeycomb plane.

\section{Crystal-field splitting and hopping parameters\label{hoppings}} 

All parameters are given in eV.

\begin{itemize}
    \item FPLO results (basis: $d_{xy},d_{yz},d_{xz},d_{z^2},d_{x^2-y^2}$)

On-site Hamiltonian 
\begin{align}
\mathbb{H}_{\rm CFS} = \left(\begin{array}{rrrrr}
 -0.906&  0.044&  0.044& -0.005&  0.034\\
  0.044& -0.906&  0.044& -0.027& -0.021\\
  0.044&  0.044& -0.906&  0.032& -0.013\\
 -0.005& -0.027&  0.032&  0.144&  0.000\\
  0.034& -0.021& -0.013&  0.000&  0.144
 \end{array} \right)
 \nonumber
\end{align}

Nearest neighbor hoppings
\begin{align}
\mathbb{T}_{1} = \left(\begin{array}{rrrrr}
-0.296&  0.044&  0.021&  0.049& -0.003\\
 0.044&  0.062& -0.019&  0.004& -0.019\\
 0.021& -0.019&  0.067&  0.012&  0.038\\
 0.049&  0.004&  0.012& -0.038& -0.004\\
-0.003& -0.019&  0.038& -0.004& -0.041
 \end{array} \right)
 \nonumber
\end{align}

3rd nearest neighbor hoppings
\begin{align}
\mathbb{T}_{3} = \left(\begin{array}{rrrrr}
 -0.038& 0.000& -0.004& 0.021& 0.012\\
 0.000& 0.000& -0.004& -0.001& 0.001\\
 -0.004& -0.004& 0.001& 0.005& 0.009\\
 0.021& -0.001& 0.005& -0.039& -0.001\\
 0.012& 0.001& 0.009& -0.001& 0.124
 \end{array} \right)
 \nonumber
\end{align}



    \item VASP Wannier90 results (basis: $d_{xy},d_{yz},d_{z^2},d_{xz},d_{x^2-y^2}$)
    

On-site Hamiltonian 
\begin{align}
\mathbb{H}_{\rm CFS} = \left(\begin{array}{rrrrr}
 3.953&  0.040& -0.005&  0.040& -0.027\\
 0.040&  3.954&  0.024&  0.040&  0.009\\
-0.005&  0.024&  4.851& -0.022& -0.000\\
 0.040&  0.040& -0.022&  3.954&  0.017\\
-0.027&  0.009& -0.000&  0.017&  4.851
 \end{array} \right)
 \nonumber
\end{align}

Nearest neighbor hoppings
\begin{align}
\mathbb{T}_{1} = \left(\begin{array}{rrrrr}
-0.296&   0.033&   0.049&  0.028&  0.017\\
 0.033&   0.067&  -0.002& -0.020& -0.019\\
 0.049&  -0.002&  -0.037&  0.013&  0.001\\
 0.028&  -0.020&   0.013&  0.067&  0.035\\
-0.017&  -0.019&   0.001&  0.035& -0.042
 \end{array} \right)
 \nonumber
\end{align}

3rd nearest neighbor hoppings
\begin{align}
\mathbb{T}_{3} = \left(\begin{array}{rrrrr}
 -0.039&  0.000&  0.021& -0.004&  0.002\\
  0.000&  0.001& -0.002& -0.003&  0.003\\
  0.021& -0.002& -0.036&  0.007&  0.001\\
 -0.004& -0.003&  0.007&  0.001&  0.013\\
  0.002&  0.003&  0.001&  0.013&  0.125 
 \end{array} \right)
 \nonumber
\end{align}
    
\end{itemize}

We would like to note that on-site energies and hopping parameters obtained using PAW are very close to what one gets using Wannier function projection technique. 

\section{DFT: Bonding-antibonding $t_{2g}$ bands due to a direct hopping}

In order to clarify importance of the direct hopping between $t_{2g}$ orbitals centered at the neighboring sites we performed Wannier function projection~\cite{Mostofi2014} of DFT Hamiltonian as obtained in VASP. The real space Hamiltonian was transformed to the local coordinate system with axis directed to ligands. Then for each pair nearest Co atoms  $xy/xy$ hopping was put to zero and resulting Hamiltonian rotated back to the global coordinate system and Fourier transformed to the reciprocal space. In Fig.~\ref{t2gbands} resulting band structure (red) is shown together with the one obtained by the tight-binding method using Wannier function projected Hamiltonian (black). One can clearly see that this is the $xy/xy$ hopping, which leads to formation of the gap between two branches of $t_{2g}$ bands.
\begin{figure}[t]
\includegraphics[width=0.45\textwidth]{SM_t2g_bands.pdf}
\caption{\label{t2gbands} Black: band-structure calculated within the tight-binding method using parameters obtained by the Wannier function projection in VASP. Red: the same, but with the $xy/xy$ hopping between nearest neighbors put to zero.}
\end{figure}


\section{DFT: Bonding-antibonding $e_{g}$ bands due to third neighbor hopping}
There is a substantial hopping between ${x^2-y^2}$ orbitals centered on third nearest neighbors. We therefore project the density of states onto bonding and antibonding $x^2-y^2$ orbitals. This is technically easier to perform in Wien2k code~\cite{Foyevtsova2013} and results of the projection are shown in Fig.~\ref{qmonn3}. 

\begin{figure}[ht]
\includegraphics[width=0.45\textwidth]{qmodosNN3.pdf}
\caption{\label{qmonn3} Nonrelativistic GGA density of states projected onto bonding and antibonding $d_{x^2-y^2}$ orbitals.}
\end{figure}

\section{DFT+U+SOC electronic structure}
The electronic structure obtained in GGA+U+SOC calculations is presented in Fig.~\ref{DOS}. 
\begin{figure}[t!]
\includegraphics[width=0.45\textwidth]{BaCoAsO-ggau-soc.png}
\caption{\label{DOS} Total density of states as obtained in GGA+U+SOC calculations. Fermi energy is in zero.}
\end{figure}


\section{DFT+U+SOC: Atomic position relaxation}
Relaxation of the atomic positions was performed keeping the volume and unit cell shape constant while variation of the total energy between ionic iterations was more than 10$^{-6}$ eV per unit cell. The hopping parameters obtained by the subsequent Wannier function projection of non-magnetic DFT Hamiltonian are given below. 


On-site Hamiltonian 

\begin{align}
\mathbb{H}_{\rm CFS} = \left(\begin{array}{ccccc}
 3.770& 0.038&-0.003& 0.039&-0.013\\
 0.038& 3.770& 0.011& 0.039& 0.002\\
-0.003& 0.011& 5.024&-0.009& 0.000\\
 0.039& 0.039&-0.009& 3.771& 0.008\\
-0.013& 0.002& 0.000& 0.008& 5.022
 \end{array} \right)
 \nonumber
\end{align}

Nearest neighbor hoppings
\begin{align}
\mathbb{T}_{1} = \left(\begin{array}{ccccc}
-0.284& 0.032& 0.058& 0.020&-0.014\\
 0.032& 0.066& 0.001&-0.018&-0.022\\
 0.058& 0.001&-0.048& 0.017& 0.000\\
 0.020&-0.018& 0.017& 0.067& 0.039\\
-0.014&-0.022& 0.000& 0.039&-0.035
 \end{array} \right)
 \nonumber
\end{align}

3rd nearest neighbor hoppings
\begin{align}
\mathbb{T}_{3} = \left(\begin{array}{ccccc}
-0.043& 0.002& 0.026&-0.003& 0.005\\
 0.002& 0.000&-0.004&-0.003& 0.007\\
 0.026&-0.004&-0.039& 0.004& 0.000\\
-0.003&-0.003& 0.004& 0.000& 0.009\\
 0.005& 0.007& 0.000& 0.009& 0.129
 \end{array} \right)
 \nonumber
\end{align}
  

We did not perform total energy calculations of the exchange parameters in the relaxed structure because of the too large unit cell. Instead the nearest neighbor exchange matrix was calculated using exact diagonalization method. Results are shown in Table \ref{JED}. As can be seen, there are only small modifications compared to the results for the experimental structure presented in the main text, thus highlighting the strong sensitivity of the ground state to the couplings. 



\begin{table}[t!]
\centering
\caption{Exchange interaction parameters (in units of K) computed from exact diagonalization of effective $d$-orbital model. Values in brackets include estimated corrections for omitted ligand exchange processes.}
\begin{tabular}{l | c | c c c}
\hline
\hline
$J_{H,t_{2g}}$&0.7~eV&\multicolumn{3}{c}{0.9~eV}\\
 $U$&3.25 eV & 5 eV  & 6 eV & 7 eV \\
\hline 
 $\text{J}_1$ (K)                & -106 (-126)& -36 (-56)  & -16 (-36)  &  -7.3 (-27) \\
  $K_1$ (K)               & 31 &  12   &5.9 & 2.9 \\ 
$\Gamma_1$ (K)            & 28 (35)&  14 (21) & 8.1 (15) &  4.9 (12)\\ 
  $\Gamma_1^{\prime}$ (K) & 8.8 (16) & 6 (13) & 3.5 (11) & 1.9 (9) \\ 
\hline
$\text{J}_3$ (K)       & 48 & 34   &30  & 27
\\
$K_3$ (K)              &-1.3 & -0.7 &-0.6 & -0.5
\\
$\Gamma_3$ (K)         &-24 & -14  & -12& -10
\\
$\Gamma_3^\prime$ (K)  & -24 & -14  & -12 & -10
\\
\hline
\hline
\end{tabular}
\label{JED}
\end{table}



