Order from disorder phenomena in BaCoS$_2$

TL;DR

This paper reveals that the magnetic transition in BaCoS$_2$ is driven by an order-from-disorder mechanism involving orbital degrees of freedom, explaining its finite Néel temperature despite magnetic frustration.

Contribution

It demonstrates that the magnetic transition in BaCoS$_2$ is an order-from-disorder phenomenon involving orbital effects, a novel insight into frustrated layered magnets.

Findings

Magnetic transition driven by order-from-disorder mechanism.

Orbital degrees of freedom contribute actively to the transition.

Presence of abundant low-energy excitations around $T_N$.

Abstract

At $T_N\simeq 305~\text{K}$ the layered insulator BaCoS$_2$ transitions to a columnar antiferromagnet that signals non-negligible magnetic frustration despite the relatively high $T_N$, all the more surprising given its quasi two-dimensional structure. Here, we show by combining ab initio and model calculations that the magnetic transition is an order-from-disorder phenomenon, which not only drives the columnar $C_4\to C_2$ symmetry breaking, but also, and more importantly, the inter-layer coherence responsible for the finite N\'eel transition temperature. This uncommon ordering mechanism, actively contributed by orbital degrees of freedom, hints at an abundance of low energy excitations below and, especially, above $T_N$, not in disagreement with experimental evidences, and might as well emerge in other layered correlated compounds showing frustrated magnetism at low temperature.

Tables

Table 1. Table 1: DFT+U, with U = 2.8 ​ eV 𝑈 2.8 eV U=2.8~{}\text{eV} , energies in Kelvin and per formula unit of the low-lying spin and orbital configurations in the tetragonal structure with an 8-site unit cell, assuming a magnetic order AFS- a 𝑎 a , being degenerate with AFS- b 𝑏 b . The lowest energy state sets the zero of energy. Note that some states are doubly degenerate, for instance C ​ ( x ​ y ​ y ​ x ) 𝐶 𝑥 𝑦 𝑦 𝑥 C(xyyx) is degenerate with C ​ ( y ​ x ​ x ​ y ) 𝐶 𝑦 𝑥 𝑥 𝑦 C(yxxy) as well as G ​ ( y ​ y ​ y ​ y ) 𝐺 𝑦 𝑦 𝑦 𝑦 G(yyyy) is degenerate with G ​ ( x ​ x ​ x ​ x ) 𝐺 𝑥 𝑥 𝑥 𝑥 G(xxxx) , and thus we just indicate one of them. Moreover, the table includes also configurations not allowed by the C ​ m ​ m ​ e 𝐶 𝑚 𝑚 𝑒 Cmme orthorhombic space group, which, nonetheless, represent alternative symmetry-breaking paths from the tetragonal structure. Each state is labelled by T n 𝑛 n , T referring to the tetragonal phase and n 𝑛 n being the ascending order in energy.

spin and orbital configurations	$E$(Kelvin)	#
C-AFS-a-C($xyyx$)	0	T0
G-AFS-a-C($xyyx$)	2	T1
C-AFS-a-G($xyyx$)	14	T2
G-AFS-a-G($xyyx$)	22	T3
G-AFS-a-G($xyxy$)	50	T4
C-AFS-a-G($xyxy$)	52	T5
C-AFS-a-C($xyxy$)	52	T6
C-AFS-a-C($yyyy$)	57	T7
G-AFS-a-C($xyxy$)	64	T8
G-AFS-a-C($yyyy$)	73	T9
C-AFS-a-G($xxyy$)	79	T10
C-AFS-a-G($yyyy$)	86	T11
G-AFS-a-G($yyyy$)	89	T12
G-AFS-a-C($xxyy$)	89	T13
C-AFS-a-C($xxyy$)	93	T14
G-AFS-a-G($xxyy$)	95	T15
G-AFS-a-C($xxxx$)	171	T16
C-AFS-a-C($xxxx$)	176	T17

Table 2. Table 2: Upper table: the energies of the different orbital configurations within the assumed nearest neighbour antiferromagnetic Ising model with exchange constants Γ 1 + σ ​ δ ​ Γ 1 subscript Γ 1 𝜎 𝛿 subscript Γ 1 \Gamma_{1}+\sigma\,\delta\Gamma_{1} along a 𝑎 a , Γ 1 − σ ​ δ ​ Γ 1 subscript Γ 1 𝜎 𝛿 subscript Γ 1 \Gamma_{1}-\sigma\,\delta\Gamma_{1} along b 𝑏 b , and uniform pseudo-magnetic field σ ​ B τ 𝜎 subscript 𝐵 𝜏 \sigma\,B_{\tau} . Bottom table: the values of those parameters extracted through Table  1 . We just consider the C-type orbital stacked configurations, since the G-type ones do not allow fixing B τ subscript 𝐵 𝜏 B_{\tau} .

orbital configuration	$E$	$\Delta E$
$xyyx$	$-2{\Gamma }_1$	0
$xyxy$	$-2\sigma \delta {\Gamma }_1$	$2{\Gamma }_1-2\sigma \delta {\Gamma }_1$
$xxyy$	$2\sigma \delta {\Gamma }_1$	$2{\Gamma }_1+2\sigma \delta {\Gamma }_1$
$xxxx$	$2{\Gamma }_1-\sigma B_{\tau }$	$4{\Gamma }_1-\sigma B_{\tau }$
$yyyy$	$2{\Gamma }_1+\sigma B_{\tau }$	$4{\Gamma }_1+\sigma B_{\tau }$

Table 3. Table 3: Same as in Table  1 but for the orthorhombic structure with b > a 𝑏 𝑎 b>a , b / a = 1.008 𝑏 𝑎 1.008 b/a=1.008 . In this case, AFS-a and AFS-b are not degenerate, and thus both have been studied. Only the orbital configurations allowed by symmetry are shown. The states are labelled by O n 𝑛 n , where O refers to the orthorhombic phase and n 𝑛 n is the order.

spin and orbital configurations	$E$(Kelvin)	#
C-AFS-a-C($xyyx$)	0	O0
G-AFS-a-C($xyyx$)	2	O1
C-AFS-a-G($xyyx$)	14	O2
C-AFS-b-C($xyyx$)	20	O3
G-AFS-a-G($xyyx$)	22	O4
G-AFS-b-C($xyyx$)	22	O5
C-AFS-b-G($xyyx$)	34	O6
G-AFS-b-G($xyyx$)	42	O7
C-AFS-b-C($xxxx$)	50	O8
G-AFS-b-C($xxxx$)	65	O9
C-AFS-b-G($xxxx$)	79	O10
G-AFS-b-G($xxxx$)	82	O11
C-AFS-a-C($yyyy$)	85	O12
C-AFS-a-G($xxxx$)	93	O13
G-AFS-a-G($xxxx$)	96	O14
G-AFS-a-C($yyyy$)	101	O15
G-AFS-a-C($xxxx$)	160	O16
C-AFS-a-C($xxxx$)	165	O17
G-AFS-b-C($yyyy$)	203	O18
C-AFS-b-C($yyyy$)	209	O19

Table 4. Table 4: Leading hopping processes T ( n x , n y , n z ) subscript 𝑇 subscript 𝑛 𝑥 subscript 𝑛 𝑦 subscript 𝑛 𝑧 T_{(n_{x},n_{y},n_{z})} , where 𝐫 = ( n x , n y , n z ) 𝐫 subscript 𝑛 𝑥 subscript 𝑛 𝑦 subscript 𝑛 𝑧 \mathbf{r}=(n_{x},n_{y},n_{z}) identifies the bond connecting Co(1), see Fig.  2 , to another cobalt at distance 𝐫 𝐫 \mathbf{r} . The bonds emanating from Co(2) are obtained by the non-symmorphic symmetry, which, in particular, implies n z → − n z → subscript 𝑛 𝑧 subscript 𝑛 𝑧 n_{z}\to-n_{z} . All hopping processes are written as matrices in the subspace ( d x ​ z , d y ​ z ) subscript 𝑑 𝑥 𝑧 subscript 𝑑 𝑦 𝑧 \big{(}d_{xz},d_{yz}\big{)} . The values, in meV, are obtained by the 5-orbital model restricted to the ( d x ​ z , d y ​ z ) subscript 𝑑 𝑥 𝑧 subscript 𝑑 𝑦 𝑧 \big{(}d_{xz},d_{yz}\big{)} subspace.

Equations

$$J_{ax}$$

$${\mathbf{S}}_{\ell ,n}(\mathbf{q})$$

$${\mathbf{S}}_{\ell }(\mathbf{q},q_z)$$

$$\mathbf{S}(\mathbf{q},q_z)$$

$$H$$

$$\displaystyle\qquad\qquad\qquad+J_{3}\,\Big{(}\gamma(\mathbf{q},\delta_{3})\,\mathbf{S}_{1,n}(\mathbf{q})\cdot\mathbf{S}_{2,n+1}(-\mathbf{q})+H.c.\Big{)}\bigg{\}}\equiv H_{2}+H_{1}+H_{3}$$

$$=\frac{1}{NL}\sum _{\mathbf{q},q_z}{\mathbf{S}}^{†}(\mathbf{q},q_z)\left(\begin{array}{cc}{J}_2(\mathbf{q})& J_1\gamma (\mathbf{q},{\delta }_1{)}^{\ast }+J_3\gamma (\mathbf{q},{\delta }_3{)}^{\ast }{\text{e}}^{-i{q}_z}\\ J_1\gamma (\mathbf{q},{\delta }_1)+J_3\gamma (\mathbf{q},{\delta }_3){\text{e}}^{i{q}_z}& J_2(\mathbf{q})\end{array}\right)\mathbf{S}(\mathbf{q},q_z)$$

$$\equiv \frac{1}{NL}\sum _{\mathbf{q},q_z}{\mathbf{S}}^{†}(\mathbf{q},q_z)\hat{J}(\mathbf{q},q_z)\mathbf{S}(\mathbf{q},q_z),$$

$$J_2(\mathbf{q})$$

$$\gamma (\mathbf{q},\delta )$$

$${\lambda }_{Q_z}=-2J_2-2\sqrt{J_1^2{\delta }_1^2+J_3^2{\delta }_3^2+2J_1{J}_3{\delta }_1{\delta }_3\cos Q_z}.$$

$${\mathbf{S}}_{i,n}(\mathbf{q})$$

$$\mathbf{Q}$$

$${\mathbf{S}}_{i,n}(\mathbf{q})\cdot {\mathbf{n}}_{3,i,n}$$

$${\mathbf{S}}_{i,n}(\mathbf{q})\cdot {\mathbf{n}}_{1,i,n}$$

$${\mathbf{S}}_{i,n}(\mathbf{q})\cdot {\mathbf{n}}_{2,i,n}$$

$$\displaystyle\Big{[}\,x^{\phantom{\dagger}}_{i,n}(\mathbf{q})\,,\,p_{j,m}^{\dagger}(\mathbf{q}^{\prime})\,\Big{]}$$

$${\Pi }_{i,n}(\mathbf{q}-\mathbf{Q})$$

$$\displaystyle\qquad\qquad+p_{i,n}^{\dagger}(\mathbf{k})\,p^{\phantom{\dagger}}_{i,n}(\mathbf{k}+\mathbf{q}-\mathbf{Q})-\delta_{\mathbf{q},\mathbf{Q}}\Big{)}\,.$$

$$H_2$$

$$\displaystyle\qquad\qquad\qquad\qquad\quad+p_{i,n}^{\dagger}(\mathbf{q}-\mathbf{Q})\,p^{\phantom{\dagger}}_{i,n}(\mathbf{q}-\mathbf{Q})\Big{)}\,,$$

$$H_1$$

$$H_3$$

$${\mathbf{X}}_{i,n}(\mathbf{q})={\mathbf{n}}_{1,i,n}{x}_{i,n}(\mathbf{q})+{\mathbf{n}}_{2,i,n}{p}_{i,n}(\mathbf{q}-\mathbf{Q}).$$

$$x_{i,n}^{}(\mathbf{q})$$

$$p_{i,n}^{}(\mathbf{q})$$

$$K(\mathbf{q}{)}^2$$

$$H_2$$

$${\omega }_2(\mathbf{q})$$

$$S(\beta )$$

$$F_2$$

$$\displaystyle=-\frac{\displaystyle\Xi_{2}(T)}{\displaystyle J_{2}}\,\sum_{n}\bigg{[}J_{1}^{2}\,\big{(}{\boldsymbol{n}}_{3,1,n}\cdot{\boldsymbol{n}}_{3,2,n}\big{)}^{2}$$

$$\displaystyle\qquad\qquad\qquad\qquad+J_{3}^{2}\,\big{(}{\boldsymbol{n}}_{3,1,n}\cdot{\boldsymbol{n}}_{3,2,n+1}\big{)}^{2}\bigg{]}\,,$$

$${\Xi }_2(T)$$

$${\left(\frac{{\omega }_2(\mathbf{q})}{{\omega }_{\lambda }^2+{\omega }_2(\mathbf{q}{)}^2}\right)}^2>0,$$

$$F_4$$

$$\displaystyle\qquad\Big{[}\big{(}{\boldsymbol{n}}_{3,1,n}\cdot{\boldsymbol{n}}_{3,2,n+1}\big{)}+\big{(}{\boldsymbol{n}}_{3,1,n-1}\cdot{\boldsymbol{n}}_{3,2,n}\big{)}\Big{]}\,,$$

$${\Xi }_4(T)$$

Full text

There is growing interest in the BaCo1-xNixS2, $x\in[0,1]$, series for it displays a rich physics comprising renormalized Dirac states [1, 2], non-Fermi liquid behavior [3, 4], and insulator-to-metal transitions as function of both doping and pressure [5, 6, 7, 8, 9, 10], typical manifestations of underlying strong electronic correlations. Here, we focus on the left-hand side of that series, i.e., $\text{Ba}\text{Co}\text{S}_{2}$, which undergoes at the Néel temperature $T_{N}\simeq 305~{}\text{K}$ [11], more recently estimated to be $T_{N}\simeq 290~{}\text{K}$ [12], a transition from a correlated paramagnetic insulator to an antiferromagnetic one, characterised by columnar spin-ordered planes, which we hereafter refer to as antiferromagnetic striped (AFS) order. The planes are stacked ferromagnetically along the $c$-axis, so called C-type stacking as opposed to the antiferromagnetic G-type one.

Inelastic neutron scattering experiments show that magnetic excitations below $T_{N}$ have pronounced two-dimensional (2D) nature [13, 14]. For instance, a direct estimate of spin exchange constants by neutron diffraction has been recently attempted in doped tetragonal BaCo0.9Ni0.1S1.9 subject to an uniaxial strain [14]. This compound undergoes a Néel transition to the C-type AFS phase at 280 K [14], not far from $T_{N}$ of undoped BaCoS2. The spin model used to fit neutron data was a conventional $J_{1}-J_{2}$ Heisenberg model [15] on each plane, supplemented by an inter-plane ferromagnetic exchange $J_{c}$ estimated to be $\lesssim 0.04~{}\text{meV}$. The latter is much too small to explain through the estimated $J_{2}\sim 9.3~{}\text{meV}$ and $J_{1}\sim-2.3~{}\text{meV}$ [14] the large value of $T_{N}\simeq 290~{}\text{K}$, which is therefore puzzling and highly surprising.

Another startling property is the anomalously broad peak of the magnetic susceptibility at $T_{N}$ [11, 12], which hints at a transition in the Ising universality class rather than the expected Heisenberg one [11]. A possible reason of this behavior might be the spin-orbit coupling [11]. Indeed, a Rashba-like spin-orbit coupling due to the layered structure and the staggered sulfur pyramid orientation, see Fig. 1, has been found to yield sizeable band splittings at specific points within the Brillouin zone, at least in metallic BaNiS2 [16]. The Rashba-like spin-orbit coupling strength may barely differ in BaCoS2, or be weakened by strong correlations [17]. In either case, its main effect is to introduce an easy plane anisotropy, as indeed observed experimentally [11], which, at most, drives the transition towards the $XY$ universality class. It is well possible that the weak orthorhombic distortion may turn the easy plane into an easy axis, but the resulting magnetic anisotropy should be negligibly small, being a higher-order effect, and thus unable to convincingly explain the experimental observations.

Our aim here is to shed light into those puzzles, which we achieve by a combination of ab initio calculations and model analyses.

Before entering into the details of our work, it is worth anticipating the results and placing them in a more general framework. Frustrated magnets, as BaCoS2 certainly is, often display a continuous accidental degeneracy of the classical ground state that leads to the existence of pseudo-Goldstone modes within the harmonic spin-wave approximation [18]. Since those modes are not protected by symmetry, they may acquire a mass once anharmonic terms are included in the spin-wave Hamiltonian. This mass, in turn, cuts off the singularities brought about by the pseudo-Goldstone modes, in that way stabilizing ordered phases otherwise thwarted by fluctuations. Put differently, let us imagine that the classical potential has a manifold of degenerate minima generally not invariant under the symmetry group of the Hamiltonian. It follows that the eigenvalues of the Hessian of the potential change from minimum to minimum. Allowing for quantum or thermal fluctuations is therefore expected to favour the minima with lowest eigenvalue, although the two kinds of fluctuations not necessarily select the same ones [19]. Moreover, it is reasonable to assume that the minima with lowest eigenvalues are those that form subsets invariant under a symmetry transformation of the Hamiltonian, so that choosing any of them corresponds to a spontaneous symmetry breaking. Such a phenomenon, also known as order from disorder [20], emerges in many different contexts [21], from particle physics [22, 23] to condensed matter physics [24, 25], even though frustrated magnets still provide the largest variety of physical realisations [26, 20, 27, 28, 29, 15, 30, 31, 18, 19].

Our results indicate that the antiferromagnetic insulator BaCoS2 is another representative of the class of frustrated magnets where order from disorder effects play a very active role and may explain at once the puzzling phenomenology of this material, including very recent experimental findings, presented and discussed in a parallel publication [12].

I Phase diagram of $\text{Ba}\text{Co}\text{S}_{2}$

$\text{Ba}\text{Co}\text{S}_{2}$ is a metastable layered compound that, quenched from high temperature, crystallises in an orthorhombic structure with space group $Cmme$, no. 67 [32], characterised by in-plane primitive lattice vectors $\mathbf{a}\not=\mathbf{b}$. However, we believe physically more significant to consider as reference structure the higher-symmetry non-symmorphic $P4/nmm$ tetragonal one, thus $\mathbf{a}=\mathbf{b}$, of the opposite end member, BaNiS2, and regard the orthorhombic distortion, which may occur through either $b>a$ or vice versa, as an instability driven by the substitution of Ni with the more correlated Co. The hypothetical tetragonal phase of $\text{Ba}\text{Co}\text{S}_{2}$ is shown in Fig. 1(A). Each CoS $a-b$ plane has two inequivalent cobalt atoms, Co(1) and Co(2), see Fig. 1(B), which are related to each other by the non-symmorphic symmetry.

Below $T_{N}\simeq 290~{}\text{K}$ [12], a striped antiferromagnetic (AFS) phase sets in. In the $a-b$ plane it consists of ferromagnetic chains, either along $a$ (AFS-a) or $b$ (AFS-b), coupled antiferromagnetically, see Fig. 1(C). The stacking between the planes is C-type, i.e., ferromagnetic, thus the labels C-AFS-a and C-AFS-b that we shall use, as well as G-AFS-a and G-AFS-b whenever we discuss the G-type configurations with antiferromagnetic stacking. We mention that the orthorhombic distortion with $b>a$ is associated with C-AFS-a, i.e., ferromagnetic bonds along $a$, and vice versa for $b<a$ [11], at odds with the expectation that ferromagnetic bonds are longer than antiferromagnetic ones. This counterintuitive behaviour represents a key test for the ab initio calculations that we later present.

Let us now briefly discuss how the electronic properties change across the magnetic transition. Above $T_{N}$, BaCoS2 is a very bad conductor with no evidence of a Drude peak and with an optical conductivity that grows linearly in frequency up to around 1 eV [33]. The dc resistivity shows an activated behaviour with activation energy $\simeq 91~{}\text{meV}$ that persists above $T_{N}$ up to $400~{}\text{K}$ [33]. At the Néel transition, a reduction of low-energy spectral weight starts, which is compensated by an increase at an energy of about $40~{}\text{meV}$ [33]. Such $40~{}\text{meV}$-gap opening occurs gradually below $T_{N}$. On the contrary, the supposed Mott or charge-transfer peak at $\sim 1~{}\text{eV}$ [34, 33] is rather insensitive to the transition. These evidences suggest that, irrespective of the system being a poor metal or a weak semiconductor above $T_{N}$, there is an abundance of low energy excitations both above the Néel transition as well as below in the insulating phase.

Neutron scattering refinement and magnetic structure modelling in the low-temperature phase point to an ordered moment of $\mu_{Co}\sim 2.63-2.9\mu_{B}$ [11, 35], suggesting that each Co2+ is in a $S=3/2$ spin configuration, in agreement with the high-temperature magnetic susceptibility [11]. Moreover, the form factor analysis of the neutron diffraction data [35] indicates that the three 1/2-spins lie one in the $d_{3z^{2}-r^{2}}$, the other in the $d_{x^{2}-y^{2}}$, and the third either in the $d_{xz}$ or $d_{yz}$ $3d$-orbitals of Co. Since $d_{xz}$ and $d_{yz}$, which we hereafter denote shortly as $x$ and $y$ orbitals, form in the $P4/nmm$ tetragonal structure a degenerate $E_{g}$ doublet occupied by a single hole, such degeneracy is going to be resolved at low-temperature. That hints at the existence of some kind of orbital order, besides the spin one, in the magnetic orthorhombic phase. Let us try to anticipate by symmetry arguments which kind of order can be stabilised.

We observe that in the $Cmme$ orthorhombic structure the cobalt atoms occupy the Wyckoff positions $4g$, which, for convenience, we denote as $\text{Co}(1)\equiv(0,0,z)$, $\text{Co}(2)\equiv(1/2,0,-z)$, $\text{Co}(3)\equiv(0,1/2,-z)$, $\text{Co}(4)\equiv(1/2,1/2,z)$, and have symmetry $mm2$. As a consequence, the hole must occupy either the $x$ orbital or the $y$ one, but not a linear combination, and the chosen orbital must be the same for Co(1) and Co(4), as well as for Co(2) and Co(3). Therefore, if we denote as $d_{n}$, $d=x,y$, the orbital occupied by the hole on Co($n$), $n=1,\dots,4$, and as $d_{1}\,d_{2}\,d_{3}\,d_{4}$ a generic orbital configuration, then there are only four of them that are symmetry-allowed: $xxxx$, $yyyy$, $xyyx$ and $yxxy$, see Fig. 2.

We remark that the density matrix of a unit cell in the lattice model can have finite diagonal components for all four configurations, while off-diagonal terms are prohibited by symmetry. We further note that $xxxx$ is degenerate with $yyyy$ in the tetragonal phase. The choice of either of them is associated with the same $C_{4}\to C_{2}$ symmetry breaking that characterises both the AFS-a or AFS-b spin order and the orthorhombic distortion, $b>a$ or $a>b$. All these three choices can be associated with three Ising variables $\tau$, $\sigma$ and $X$ such that $\tau=+1$ corresponds to $xxxx$, $\sigma=+1$ to AFS-a, $X=+1$ to $b>a$, and vice versa. Since they all have the same symmetry, odd under $C_{4}$, they would be coupled to each other should we describe the transition by a Landau-Ginzburg functional. We shall hereafter denote as $Z_{2}(C_{4})$ the Ising sector that describes the $C_{4}\to C_{2}$ symmetry breaking.

The other two allowed orbital configurations $xyyx$ and $yxxy$, see Fig. 2, are instead degenerate both in the tetragonal and orthorhombic phases, but break the non-symmorphic symmetry (NS) that connects, e.g., Co(1) with Co(2) and Co(3). We can therefore associate to those configurations a new Ising sector $Z_{2}(\text{NS})$.

We emphasise that the above conclusions rely on the assumption of a $Cmme$ space group. A mixing between $x$ and $y$ orbitals is instead allowed by the $Pba2$ space group proposed in Ref. [36] as an alternative scenario for BaCoS2 at room temperature. As a matter of fact, the two symmetry-lowering routes, $P4/nmm\to Cmme$ and $P4/nmm\to Pba2$, correspond to different Jahn-Teller-like distortions involving the $d_{xz}$-$d_{yz}$ doublet and the $E_{g}$ phonon mode of the $P4/nmm$ structure at the $\mathbf{M}$ point, which is found to have imaginary frequency by ab initio calculations [36]. However, the latest high-accuracy X-ray diffraction data of Ref. [12] confirm the $Cmme$ orthorhombic structure even at room temperature, thus supporting our assumption.

II Ab initio analysis

We carried out ab initio DFT and DFT+U calculations using the Quantum ESPRESSO package [37, 38]. The density functional is of GGA type, namely the Perdew-Burke-Ernzerhof (PBE) functional [39], on which local Hubbard interactions were added in case of the GGA+U fully rotational invariant framework [40, 41]. If not stated otherwise, the geometry of the unit cell and the internal coordinates of the atomic positions in the orthorhombic structure were those determined experimentally, taken from Ref. 32. For paramagnetic calculations, the relative atomic positions were kept fixed and the in-plane lattice constants $a=b$ chosen such that the unit cell volume matched the one of the orthorhombic structure. Co and S atoms are described by norm-conserving pseudopotentials (PP) with non-linear core corrections, Ba atoms are described by ultrasoft pseudopotentials. The Co PP contains 13 valence electrons (3s2,3p6,3d7), the Ba PP 10 electrons (5s2,5p6,6s2), and S PPs are in a (3s2,3p3) configuration. The plane-waves cutoff has been set to $120$Ry and we used a Gaussian smearing of $0.01$Ry. The $k$-point sampling of the electron-momentum grid was at least $8\times 8\times 8$ points.

To determine the band structure and derive an effective low-energy model, we performed a Wannier interpolation with maximally localized Wannier functions (MLWF) [42, 43] using the Wannier90 package [44].

In the first place, we checked if the tetragonal phase is unstable towards magnetism, considering both a conventional Néel order (AFM) compatible with the bipartite lattice and the observed AFS. We found, using $U=2.8~{}\text{eV}$ as motivated by a cRPA analysis [4], that the lowest energy state is indeed the AFS, the AFM and paramagnetic phases lying above by about 0.5 eV and 2.3 eV, respectively.

Let us therefore restrict our analysis to AFS and stick to $U=2.8~{}\text{eV}$. We use an 8-site unit cell that includes two planes, which allows us to compare C-AFS with G-AFS. In addition, we consider both the tetragonal structure with AFS-a, since AFS-b is degenerate, and the orthorhombic structure with $b>a$, in which case we analyse both AFS-a and AFS-b. For all cases, we consider all four symmetry-allowed orbital configurations, $xxxx$, $yyyy$, $xyyx$ and $yxxy$, assuming either a C-type or G-type orbital stacking between the two planes of the unit cell, so that, for instance, $\text{G}(xxxx)$ means that one plane is in the $xxxx$ configuration and the other in the $yyyy$ one.

In Table 1 we report the energies per formula unit of several possible configurations in the tetragonal structure, including those that would be forbidden in the orthorhombic one. All energies are measured with respect to the lowest energy state and are expressed in Kelvin. In agreement with experiments, the lowest energy state T0 has spin order C-AFS, a or b being degenerate. In addition, it has C-type antiferro-orbital order, C$(xyyx)$. Remarkably, its G-type spin counterpart T1 is only 2 K above, supporting our observation that $T_{N}=290~{}\text{K}$ is anomalously large if compared to these magnetic excitations.

The energy differences between C-type orbital stacked configurations and their G-type counterparts are too small and variable to allow obtaining a reliable modelling of the inter-plane orbital coupling.

On the contrary, the energy differences between in-plane orbital configurations can be accurately reproduced by a rather simple modelling. We assume on each Co-site an Ising variable $\tau_{3}$ equal to the difference between the hole occupations of orbital $x$ and of orbital $y$. The Ising variable on a given site is coupled only to those of the four nearest neighbour sites in the $a-b$ plane, with exchange constants $\Gamma_{1a}=\Gamma_{1}+\sigma\,\delta\Gamma_{1}$ and $\Gamma_{1b}=\Gamma_{1}-\sigma\,\delta\Gamma_{1}$ along $a$ and $b$, respectively. In addition, the Ising variables feel a uniform field $B_{\tau}\,\sigma$. Here, $\sigma$ is the Ising $Z_{2}(C_{4})$ order parameter that distinguishes AFS-a, $\sigma=+1$ from AFS-b, $\sigma=-1$.

We find that the spectrum is well reproduced by the parameters in Table 2.

We now move to the physical orthorhombic structure, assuming $b>a$ with $b/a=1.008$, and recalculate all above energies but considering only the orbital configurations allowed by the $Cmme$ space group. In this case, we have to distinguish between AFS-a and AFS-b, which are not anymore degenerate. The results are shown in Table 3.

We remark that the ab initio calculation correctly predicts that the lowest energy state O0 has ferromagnetic bonds along $a$ despite $b>a$, which, as we mentioned, is an important test for the theory. The energy difference between AFS-a and AFS-b, e.g., O0 and O3, is about 20 K, and gives a measure of the spin-exchange spatial anisotropy in the $a$ and $b$ directions due to the orthorhombic distortion. This small value implies that the Néel transition temperature $T_{N}\simeq 290~{}\text{K}$ is largely insensitive to the orthorhombic distortion that exists also above $T_{N}$ [11, 12]. In particular, the energy difference between C-type and G-type stacking, O0 and O1 in Table 3, remains the same tiny value found in the tetragonal phase. Table 3 thus suggests that the spin configurations C-AFS-a, C-AFS-b, G-AFS-a and G-AFS-b are almost equally probable at the Néel transition, and that despite the orthorhombic structure.

Concerning the effective Ising model that controls the in-plane orbital configuration, we find that the parameter that is most affected by the orthorhombic distortion is $B_{\tau}$, which decreases in the AFS-a configurations and increases in AFS-b ones.

The calculated magnetic moment per Co atom in the lowest energy state, O0 in Table 3, is $\mu_{AFS}\sim 2.65~{}\mu_{B}$, in quite good agreement with experiments [11, 35]. We mention that the actual magnetic moment of Co is lower, $\mu_{AFS}^{Co}\sim 2.29~{}\mu_{B}$. This difference is due to the strong Co-S hybridization that yields a non-negligible magnetic moment on the sulfur sites.

II.1 Orthorhombic distortion.

BaCoS2 has a small orthorhombic distortion quantified by $d=\frac{a-b}{\bar{a}}\sim 0.8\%$ [32] with $\bar{a}=\frac{a+b}{2}$, which is further reduced to $d\sim 0.4\%$ under high pressure synthesis [12]. Therefore, the stability of the different phases with respect to such an orthorhombic distortion can be regarded as a further criterion to identify the correct electronic ground state of the system. Whereas the PM phase favours a strong orthorhombic deformation of $d\sim 8\%$, the magnetic stripe phases have minimal energy variations for small distortions. The precise values differ between the different orbital configurations as shown in Fig. 3, and the lowest energy is found for the orbital ordered configuration $xyyx$ at $d^{xyyx}_{\mathrm{min}}\sim 0.5\%$. For the nematic $xxxx$ configuration, the distortion is in good agreement with the experimental values ($d^{xxxx}_{\mathrm{min}}\sim-1.0\%$) albeit being much higher in energy. The energy minimum of the nematic $yyyy$ phase is competing with $xyyx$, but located at positive instead of negative distortion ($d^{yyyy}_{\mathrm{min}}\sim 0.8\%$).

II.2 Insulator-to-metal transition.

The insulator-metal transition can be controlled either by chemical doping, for instance with Ni atoms [5], or by applying pressure [35] . We focus here on the latter case, where a critical pressure of $p_{cr}\sim 1.3$ GPa [6, 8, 7, 10] was found in experiments to be sufficient to render the system metallic. Close to $p_{cr}$ the transition from AFM insulator to PM metal occurs via the formation of an intermediate antiferromagnetic metal phase.

Based on the structural changes under pressure reported in Ref. 7, we performed DFT+U calculations for the low-energy solutions identified earlier 111In order to be consistent with the rest of the paper, we used the relative structural changes of the lattice constants under pressure reported in Ref. 7 and applied them to the ambient pressure data of Ref. 32.. In Ref. 7 the structural parameters show a rather continuous evolution upon pressure within the metallic and insulating regimes. At the phase transition, however, the structural changes are large. Nevertheless, one can parametrize the distinct low- and high-pressure regimes separately whose structural parameters vary as a function of pressure. Here, we focus on the low-pressure parametrization and compare the energy gap of the two C-AFS-a competing solutions C(xyyx) and C(yyyy) as a function of pressure. In Fig. 4 we plot the density of states (DOS) of both configurations for several pressures. We note that, while C(xyyx) remains insulating at all pressures, the gap of the C(yyyy) configuration closes around $p\simeq 8~{}\text{kbar}$. Such behaviour is robust either upon choosing the atomic positions of Ref. 7 or those of Ref. 32. On the contrary, the relative energies of C(xyyx) and C(yyyy), which are extremely sensitive to the distance between Co and apical S, do depend on that choice. Indeed, Ref. 7 and Ref. 32 report significantly different Co-apical S distances. The insulating C(xyyx) remains the lowest energy state at all investigated pressures using the atomic positions of Ref. 32, see top panel of Fig. 4(B). Conversely, the reduced Co-apical S distance of the crystal structure reported in Ref. 7 predicts a first order transition at $p\sim 12$ kbar between an insulating C(xyyx) and a metallic C(yyyy), consistent with the metal-insulator transition reported experimentally, see bottom panel. Should the latter scenario be representative of BaCoS2 under pressure, it would imply that the metal-insulator transition is also accompanied by the recovery of the non-symmorphic symmetry.

II.3 Wannierization

In order to gain further insight into the electronic structure of BaCoS2 in general and, in particular, of the bands with pronounced $d_{xz}$ and $d_{yz}$ character, we generate two theoretical model Hamiltonians with maximally-localized Wannier functions for Co-$d$-like and Co-$d_{xz/yz}$-like orbitals respectively. To this aim, we construct Wannier fits using the Wannier90 tool [44] based on the paramagnetic DFT+U calculation using a $4\times 4\times 4$ k-grid with a doubled in-plane unit cell comprising 4 Co atoms.

For both models, the resulting Wannier Hamiltonian reproduces overall well the DFT band structure, see Fig. 5. Whereas the fit of the 5-orbital model is nearly perfect, the 2-orbital model shows small deviations along the $\mathbf{M}-\boldsymbol{\Gamma}$ direction due to missing hybridization with the other Co-$d$ orbitals.

Table 4 shows the leading hopping processes of the 5-band model restricted to the \big{(}d_{xz},d_{yz}\big{)} subspace.

We note that, because of the staggered shift of the Co atoms out of the sulfur basal plane, the largest intra-layer hopping is between next-nearest neighbour (NNN) cobalt atoms instead of nearest-neighbour (NN) ones. Moreover, the stacking of the sulfur pyramids and the position of the intercalated Ba atoms makes the inter-layer NN hopping negligible, contrary to the NNN one that is actually larger than the in-plane NN hopping, but still smaller than the in-plane NNN one.

We finally remark that the orthorhombic distortion has a very weak effect on the inter-layer hopping, which is consistent with the tiny energy difference between C-AFS and G-AFS being insensitive to the distortion, compare, e.g., the energies of T1 and O1 in Tables 1 and 3, respectively.

III Effective Heisenberg model

We already mentioned that the largest energy scales that emerge from the ab initio calculations are the magnetic ones separating the lowest-energy AFS configuration from the Néel and paramagnetic states. Therefore, even though BaCoS2 seems not to lie deep inside a Mott insulating regime, we think it is worth discussing qualitatively the spin dynamics in terms of an effective $S=3/2$ Heisenberg model. If we assume that the leading contribution to the exchange constants derives from the hopping processes within the $d_{xz}-d_{yz}$ subspace, then Table 4 suggests the Heisenberg model shown in Fig. 6. According to this figure, the exchange constants $J_{1x/y},J_{2},$ and $J_{3x/y}$ are related to the hopping terms $T_{(\pm 1,0,0)/(0,\pm 1,0)},T_{\pm(1,-1,0)/\pm(1,1,0)},$ and $T_{(\pm 1,0,1)/(0,\pm 1,1)}$, reported in Table 4. This model consists of frustrated $J_{1}-J_{2}$ planes [15, 46, 47, 48, 49] coupled to each other by a still frustrating $J_{3}$. In order to be consistent with the observed columnar magnetic order, the exchange constants have to satisfy the inequality $2J_{2}>\left|J_{3}\right|+\left|J_{1}\right|$. Moreover, $J_{3}$ forces to deal with a two sites unit cell, highlighted in yellow in the left panel of Fig. 6 where the non-equivalent cobalt sites are referred to as Co(1), in blue, and Co(2), in red, respectively. The reason is that Co(1) on a plane is only coupled to Co(2) on the plane above but not below, and vice versa for Co(2).

We write, for $a=1,3$,

[TABLE]

where, in analogy with the single-layer $J_{1}-J_{2}$ model [15], $\delta_{a}\not=0$ are Ising order parameters associated with the $C_{4}\to C_{2}$ symmetry breaking. Moreover, we define the in-plane Fourier transforms

[TABLE]

where $\mathbf{R}$ labels the $N$ unit cells in the $a-b$ plane, $\ell=1,2$ the two sites (sublattices) within each unit cell, and $n=1,\dots,L$ the layer index. Similarly, we introduce the three-dimensional Fourier transforms

[TABLE]

as well as the spinor

[TABLE]

With those definitions, the Hamiltonian reads

[TABLE]

where $H_{a}$ is proportional to $J_{a}$, $a=1,2,3$,

[TABLE]

The classical ground state corresponds to the modulation wavevector $\mathbf{Q}=(\pi,0,Q_{z})\equiv(0,\pi,Q_{z})$, related to an antiferromagnetic order within each sublattice on each layer, yielding the lowest eigenvalue of $\hat{J}(\mathbf{q},q_{z})$, i.e.,

[TABLE]

The expression of $\lambda_{Q_{z}}$ shows that inter-layer magnetic coherence sets in only when the two Ising-like order parameters, $\delta_{1}$ and $\delta_{3}$, lock together. Specifically, $J_{1}J_{3}\delta_{1}\delta_{3}>0$ stabilises C-AFS, $Q_{z}=0$, otherwise G-AFS, $Q_{z}=\pi$. We already know that the former is lower in energy, though by only few Kelvins, see Table 3. Moreover, an orthorhombic distortion $b>a$ favours AFS-a, which implies $J_{1}\,\delta_{1}+J_{3}\,\delta_{3}>0$, even though AFS-b is higher by only $20~{}\text{K}$ according to DFT+U, see O3 in Table 3.

If we assume that the magnetic excitations in the AFS state of strained tetragonal BaCo0.9Ni0.1S1.9 are representative of those of orthorhombic BaCoS2, then the inelastic neutron scattering data of Ref. [14] at $200~{}\text{K}$ and for in-plane transfer momentum are consistent with $J_{2}\simeq 9.3~{}\text{meV}$, $J_{1}+J_{3}\simeq-2.34~{}\text{meV}$, and $J_{1}\delta_{1}+J_{3}\delta_{3}\simeq 0.53~{}\text{meV}$. However, the scattering data for out-of-plane transfer momentum only yield an upper bound to the propagation velocity of the Goldstone mode along $z$, which actually corresponds to $4J_{1}J_{3}\delta_{1}\delta_{3}<0.08~{}\text{meV}^{2}$ for $J_{3}\simeq J_{1}$. Such small bound, around a quarter of $\left(J_{1}\delta_{1}+J_{3}\delta_{3}\right)^{2}$, suggests that the two order parameters $\delta_{1}$ and $\delta_{3}$ are already formed at $80~{}\text{K}$ below $T_{N}\simeq 280~{}\text{K}$, whereas their mutual locking is still suffering from fluctuations. We finally observe that the ferromagnetic sign of $J_{1}$ and $J_{3}$ is consistent with the diagonal hopping matrices in the corresponding directions, see Table 4, and the antiferro-orbital order.

To better understand the interplay between the $Z_{2}\left(C_{4}\right)$ Ising degrees of freedom and the magnetic order at $T_{N}$, we investigate in more detail the Hamiltonian (5) with $\delta_{1}=\delta_{3}=0$, thus $J_{1x}=J_{1y}=J_{1}$ and $J_{3x}=J_{3y}=J_{3}$. Since $J_{2}>0$ is the dominant exchange process, the classical ground state corresponds to the spin configuration

[TABLE]

where $S=3/2$ is the spin magnitude, ${\boldsymbol{n}}_{3,i,n}$ is a unit vector, and

[TABLE]

the equivalence holding since ${\boldsymbol{G}}=(\pi,\pi,0)$ is a primitive lattice vector for the two-site unit cell. In other words, each sublattice on each plane is Néel ordered, and its staggered magnetisation ${\boldsymbol{n}}_{3,i,n}$ is arbitrary. We therefore expect that quantum and thermal fluctuations may yield a standard order-from-disorder phenomenon [20].

Within spin-wave approximation, the spin operators can be written as

[TABLE]

where ${\boldsymbol{n}}_{1,i,n}$, ${\boldsymbol{n}}_{2,i,n}$ and ${\boldsymbol{n}}_{3,i,n}$ are orthogonal unit vectors, $x_{i,n}^{\dagger}(\mathbf{q})=x^{\phantom{\dagger}}_{i,n}(-\mathbf{q})$ and $p_{i,n}^{\dagger}(\mathbf{q})=p^{\phantom{\dagger}}_{i,n}(-\mathbf{q})$ are conjugate variables, i.e.,

[TABLE]

and

[TABLE]

The three terms of the Hamiltonian (5) thus read, at leading order in quantum fluctuations, i.e., in the harmonic approximation,

[TABLE]

where $E_{0}=2NL\,S(S+1)\,J_{2}(\mathbf{Q})$, $\gamma(\mathbf{q})=\gamma(\mathbf{q},\delta=0)$, and

[TABLE]

We note that $H_{2}$ does not depend on the choice of ${\boldsymbol{n}}_{3,i,n}$, reflecting the classical accidental degeneracy, unlike $H_{1}+H_{3}$. We start treating $H_{1}$ and $H_{3}$ within perturbation theory. The unperturbed Hamiltonian $H_{2}$ is diagonalised through the canonical transformation

[TABLE]

where

[TABLE]

In the transformed basis,

[TABLE]

is diagonal, and the spin-wave energy dispersion is

[TABLE]

The free energy in perturbation theory can be written as $F=\sum_{\ell}\,F_{\ell}$, where $F_{\ell}$ is of $\ell$-th order in $H_{1}+H_{3}$, and $F_{0}$ is the unperturbed free energy of the Hamiltonian $H_{2}$. Notice that only even-order terms are non vanishing, thus $\ell=0,2,4,\dots$ . Given the evolution operator in imaginary time,

[TABLE]

where $H_{a}(\tau)$, $a=1,3$, evolves with the Hamiltonian $H_{2}$, the second order correction to the free energy is readily found to be

[TABLE]

where

[TABLE]

with $\omega_{\lambda}=2\pi\lambda T$, $\lambda\in\mathbb{Z}$, bosonic Matsubara frequencies. Even without explicitly evaluating $\Xi_{2}$, we can conclude that the free-energy gain at second order in $H_{1}+H_{3}$ is maximised by ${\boldsymbol{n}}_{3,1,n}\cdot{\boldsymbol{n}}_{3,2,m}=\pm 1$, with $m=n,n+1$, which reduces the classical degeneracy to $4^{L}$ configurations, where $L$ is the total number of layers in the system. Such residual degeneracy is split by a fourth order correction to the free energy proportional to $J_{1}^{2}\,J_{3}^{2}$ that reads

[TABLE]

where

[TABLE]

We remark that, despite $\omega_{2}(\mathbf{q})$ vanishes linearly at $\mathbf{q}=\mathbf{0},\mathbf{Q}$, both $\Xi_{2}(T)$ and $\Xi_{4}(T)$ are non-singular.

The fourth order correction $F_{4}$ in Eq. (22) has a twofold effect: it forces ${\boldsymbol{n}}_{3,1,n}\cdot{\boldsymbol{n}}_{3,2,n}$ to be the same on all layers and, in addition, stabilises a ferromagnetic inter-layer stacking. Therefore, the ground state manifold at fourth order in $H_{1}+H_{3}$ is spanned by ${\boldsymbol{n}}_{3,1,n}={\boldsymbol{n}}_{3}$ and ${\boldsymbol{n}}_{3,2,n}=\sigma\,{\boldsymbol{n}}_{3}$, where ${\boldsymbol{n}}_{3}$ is an arbitrary unit vector reflecting the spin $SU(2)$ symmetry, and $\sigma=\pm 1$ is associated with the global $C_{4}\to C_{2}$ symmetry breaking.

Similarly to the single-plane $J_{1}-J_{2}$ model [15], the above results imply that an additional term must be added to the semiclassical spin action. Specifically, if we introduce the Ising-like fields $\sigma_{n}(\mathbf{R})={\boldsymbol{n}}_{3,1,n}(\mathbf{R})\cdot{\boldsymbol{n}}_{3,2,n}(\mathbf{R})$ and $\sigma_{n+1/2}(\mathbf{R})={\boldsymbol{n}}_{3,1,n}(\mathbf{R})\cdot{\boldsymbol{n}}_{3,2,n+1}(\mathbf{R})$, Eqs (20) and (22) imply that, at the leading orders in $J_{1}$ and $J_{3}$, the effective action in the continuum limit includes the quadrupolar coupling term [15]

[TABLE]

We expect [15] an Ising transition to occur at a critical temperature $T_{c}$, below which $\langle\sigma_{n}(\mathbf{R})\rangle=m_{1}$, $\langle\sigma_{n+1/2}(\mathbf{R})\rangle=m_{3}$, with $m_{1}\,m_{3}>0$. However, this transition is now three-dimensional and brings along the AFS magnetic ordering, thus a finite Néel temperature $T_{N}\simeq T_{c}$.

To get a rough estimate of $T_{N}$, we assume that, upon integrating out the spin degrees of freedom, the classical action describes an anisotropic three-dimensional ferromagnetic Ising model with exchange constants $I_{1}$ on layers $n$, $I_{3}$ on layers $n+1/2$, and $I_{\perp}<I_{1}$ between layers. Hereafter, we take for simplicity $J_{1}=J_{3}$, thus $I_{1}=I_{3}\equiv I_{\parallel}$. We then note that, if we use the above estimates of $J_{2}\simeq 9.3~{}\text{meV}$ and $J_{1}=J_{3}\simeq-1.17~{}\text{meV}$, the 2D Ising critical temperature with $S=3/2$ of each layer $n$ and $n+1/2$ is about $0.4\,(S+1/2)^{2}J_{2}\simeq 173~{}\text{K}$ [47]. That corresponds to $I_{\parallel}\simeq 6.6~{}\text{meV}$ [50]. The 3D Ising critical temperature $T_{c}\simeq T_{N}$ depends on the anisotropy ratio $I_{\perp}/I_{\parallel}$ [50], being at most $345~{}\text{K}$ for $I_{\perp}/I_{1}\simeq 1$, and reaching the observed $T_{N}\simeq 290$ [12] at $I_{\perp}/I_{1}\simeq 0.75$, which we find not unrealistic.

However, BaCoS2 remains orthorhombic above $T_{N}$, which implies that the structural $C_{4}\to C_{2}$ symmetry breaking occurs earlier than magnetic ordering. Therefore, even though the effects of the orthorhombic distortion on the electronic structure are rather small, see Table 4, it is worth repeating the above discussion assuming from the start that $\sigma_{n}(\mathbf{R})={\boldsymbol{n}}_{3,1,n}(\mathbf{R})\cdot{\boldsymbol{n}}_{3,2,n}(\mathbf{R})=1$ (AFS-a), so that

[TABLE]

It follows that the quadrupolar term (24) becomes

[TABLE]

with $K(T)\gg h(T)>0$, and provides a direct coupling between adjacent planes able to stabilise the 3D magnetic order. For $h(T)=0$, C-AFS, ${\boldsymbol{n}}_{3,1,n}\cdot{\boldsymbol{n}}_{3,1,n+1}=1$, and G-AFS, ${\boldsymbol{n}}_{3,1,n}\cdot{\boldsymbol{n}}_{3,1,n+1}=-1$, are degenerate. In this case, the Néel transition occurs simultaneously with the Ising transition that makes the system choose either C-AFS or G-AFS. Moreover, since the two types of stacking are separated by a macroscopic energy barrier, $T_{N}$ can be large despite a weak spin-wave dispersion along the $c$-axis, not in disagreement with neutron scattering experiments [13, 14]. The finite $h(T)$, which DFT+U predicts to be only few kelvins at $T=0$, eventually favours C-type stacking. Nonetheless, we expect that the Néel transition should still maintain a strong Ising-like character.

IV Conclusions

We showed by ab initio calculations that the low-temperature phase of BaCoS2 critically depends on the specific Co $d$-orbitals where the three unpaired electrons lie. Indeed, the directionality of those orbitals and the significant hybridisation with ligand sulfur atoms are ultimately responsible of the magnetic frustration implied by the columnar magnetic order despite the lattice being bipartite. Moreover, the $d_{xz}-d_{yz}$ doublet hosting a single hole entails orbital degeneracy in tetragonal symmetry that has to be lifted at low temperature. However, contrary to the expectation that the orthorhombic distortion should make the hole choose either $d_{xz}$ or $d_{yz}$, we find that the lowest-energy state has an antiferro-orbital order, thus both orbitals almost equally occupied on average, which breaks the non-symmorphic symmetry and conspires with magnetism to stabilise an insulating phase. For the same reason, the orthorhombic distortion that persists above the Néel transition temperature [11, 12] does not produce the large nematicity that could efficiently remove magnetic frustration should the occupancies of $d_{xz}$ and $d_{yz}$ be very different.

Magnetic frustration and orbital degeneracy are also responsible of the multitude of states that we find within 20 meV above the lowest energy one, which corresponds to planes with columnar spin-order and antiferro-orbital order stacked ferromagnetically, both in spin and orbital, along the $c$-axis. This abundance of low-energy excitations revealed by GGA+U already hints at the order-from-disorder origin of the magnetic transition that we uncover by studying the effective Heisenberg model built from the Wannierization of the band structure, and which represents a three-dimensional generalisation of the $J_{1}-J_{2}$ Heisenberg model on a square lattice. Similarly to the latter [15], the model in Fig. 6 displays a finite-temperature order-from-disorder Ising-like transition, which is now three-dimensional and thus able to drive also a finite temperature Néel transition. This result rationalises in simple terms some of the puzzling features of BaCoS2, specifically, the strong Ising-like character of the Néel transition [11, 12], the large $T_{N}$ despite the weak $c$-axis dispersion of the magnetic excitations [13, 14], and the anomalous electronic properties above $T_{N}$ [33] hinting at the presence of strong low-energy fluctuations.

We finally draw a parallel to another class of materials where orbital, magnetic and structural order parameters are intertwined and play an important role, namely iron pnictide superconductors (FeSC) [51] and FeSe  [52, 53, 54, 55]. In many (mainly electron-doped) iron-based superconductors the disordered phase at high temperature first spontaneously breaks $C_{4}$ symmetry when cooling below a critical temperature $T_{\mathrm{nem}}$, thus entering a nematic phase [51]. Only at a lower temperature $T_{N}<T_{\mathrm{nem}}$, also spin $SU(2)$ is broken and a stripe-ordered magnetic long-range order emerges [51, 54]. The Néel temperature can be quite large as in BaCoS2  or even vanishing within experimental accuracy as in FeSe, where $T_{N}\not=0$ is observed only under pressure [56].

Also from a model point of view, our description of BaCoS2 fits into the modelling of these materials. Indeed, the key role played by several order parameters in BaCoS2 has parallels to the phenomenological Landau free energy description of FeSC [51]. Besides itinerant multi-orbital Hubbard models [57, 58, 59], also spin-1 Heisenberg models conceptually similar to ours have been proposed for iron-based superconductors and FeSe [60, 52]. There, however, instead of achieving the $C_{4}$ symmetry breaking via an order-from-disorder phenomenon, that is often assumed from the start [61] or by explicitly adding bi-quadratic spin exchanges [60] that mimic the quadrupolar terms (24).

Moreover, our ab initio simulations for BaCoS2 predict that the lowest-energy phase has not ferro-orbital order, as often discussed in the context of FeSC, but rather an anti-ferro orbital ordering that breaks the non-symmorphic symmetry instead of $C_{4}$. Therefore, BaCoS2 seems to realise a situation where collinear magnetism and orthorhombicity do not imply orbital nematicity, unless for specific crystal structures under pressure.

BaCoS2 could therefore add yet another piece to the puzzle of understanding the intricate interplay of structural, magnetic and electronic degrees of freedom in strongly correlated materials.

Acknowledgements.

We are thankful to H. Abushammala, A. Gauzzi, and Y. Klein for fruitful discussions. We acknowledge the allocation for computer resources by the French Grand Équipement National de Calcul Intensif (GENCI) under the project numbers A0010906493 and A0110912043. M.F. acknowledges financial support from the European Research Council (ERC), under the European Union’s Horizon 2020 research and innovation programme, Grant agreement No. 692670 ”FIRSTORM”.