Order from disorder phenomena in BaCoS$_2$

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.

There is growing interest in the BaCo 1−x Ni x S 2 , x ∈ [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., BaCoS 2 , which undergoes at the Néel temperature T N 305 K [11], more recently estimated to be T N 290 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 twodimensional (2D) nature [13,14]. For instance, a direct estimate of spin exchange constants by neutron diffraction has been recently attempted in doped tetragonal BaCo 0.9 Ni 0.1 S 1.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 BaCoS 2 . 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 0.04 meV. The latter is much too small to explain through the estimated J 2 ∼ 9.3 meV and J 1 ∼ −2.3 meV [14] the large value of T N 290 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 BaNiS 2 [16]. The Rashbalike spin-orbit coupling strength may barely differ in BaCoS 2 , 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 BaCoS 2 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 as-arXiv:2302.12179v1 [cond-mat.str-el] 23 Feb 2023 sume 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 [15,[18][19][20][26][27][28][29][30][31].
Our results indicate that the antiferromagnetic insulator BaCoS 2 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 BaCoS2
BaCoS 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 a = b. However, we believe physically more significant to consider as reference structure the higher-symmetry non-symmorphic P 4/nmm tetragonal one, thus a = b, of the opposite end member, BaNiS 2 , 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 BaCoS 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 290 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 , BaCoS 2 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 91 meV that persists above T N up to 400 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 meV [33]. Such 40 meV-gap opening occurs gradually below T N . On the contrary, the supposed Mott or charge-transfer peak at ∼ 1 eV [33,34] 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 µ Co ∼ 2.63 − 2.9µ B [11,35], suggesting that each Co 2+ 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/2spins 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 P 4/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 Co(1) ≡ (0, 0, z), Co(2) ≡ (1/2, 0, −z), Co(3) ≡ (0, 1/2, −z), Co(4) ≡ (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, . . . , 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 o↵-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 ! C 2 symmetry breaking that characterise both the AFS-a or AFS-b spin order and the orthorhombic distortion, b > a or a > b. All those three choices can be associated with three Ising variables ⌧ , and X such that ⌧ = +1 corresponds to xxxx, = +1 to AFS-a, X = +1 to b > a, and viceversa. 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 ! 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 nonsymmorphic symmetry (NS) that connects, e.g., Co(1) with Co(2) and Co (3). We therefore associate to those configurations a new Z 2 (NS) symmetry, and denote as m the symmetry-breaking order parameter, assuming that m > 0 means prevailing xyyx and m < 0 the opposite. 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, . . . , 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 → 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 τ , σ and X such that τ = +1 corresponds to xxxx, σ = +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 → 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 nonsymmorphic 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 (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 P ba2 space group proposed in Ref. [36] as an alternative scenario for BaCoS 2 at room temperature. As a matter of fact, the two symmetry-lowering routes, P 4/nmm → Cmme and P 4/nmm → P ba2, correspond to different Jahn-Teller-like distortions involving the d xz -d yz doublet and the E g phonon mode of the P 4/nmm structure at the 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 (3s 2 ,3p 6 ,3d 7 ), the Ba PP 10 electrons (5s 2 ,5p 6 ,6s 2 ), and S PPs are in a (3s 2 ,3p 3 ) configuration. The plane-waves cutoff has been set to 120Ry and we used a Gaussian smearing of 0.01Ry. The k-point sampling of the electron-momentum grid was at least 8 × 8 × 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 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 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, G(xxxx) means that one spin and orbital configurations E(Kelvin) T17   TABLE I. DFT+U, with U = 2.8 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, being degenerate with AFSb. The lowest energy state sets the zero of energy. Note that some states are doubly degenerate, for instance C(xyyx) is degenerate with C(yxxy) as well as G(yyyy) is degenerate with G(xxxx), and thus we just indicate one of them. Moreover, the table includes also configurations not allowed by the Cmme orthorhombic space group, which, nonetheless, represent alternative symmetry-breaking paths from the tetragonal structure. Each state is labelled by Tn, T referring to the tetragonal phase and n being the ascending order in energy.
plane is in the xxxx configuration and the other in the yyyy one. In Table I 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 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 τ 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 Γ 1a = Γ 1 + σ δΓ 1 and Γ 1b = Γ 1 − σ δΓ 1 along a and b, respectively. In addition, the Ising variables feel a uniform field B τ σ. Here, σ is the Ising Z 2 (C 4 ) order parameter that distinguishes AFS-a, σ = +1 from AFS-b, σ = −1. We find that the spectrum is well reproduced by the parameters in Table II.
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 AFSb, which are not anymore degenerate. The results are shown in Table III. 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 290 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 III, remains the same tiny value found in the tetragonal phase. Table III thus suggests that the spin configurations C-AFS-a, C-AFS-b, G-AFS-a and G-AFSb are almost equally probable at the Néel transition, and that despite the orthorhombic structure.
Concerning the effective Ising model that controls the inplane orbital configuration, we find that the parameter that is most affected by the orthorhombic distortion is B τ , which decreases in the AFS-a configurations and increases in AFS-b ones. spin and orbital configurations E(Kelvin) The calculated magnetic moment per Co atom in the lowest energy state, O0 in Table III, is µ AF S ∼ 2.65 µ B , in quite good agreement with experiments [11,35]. We mention that the actual magnetic moment of Co is lower, µ Co AF S ∼ 2.29 µ B . This difference is due to the strong Co-S hybridization that yields a non-negligible magnetic moment on the sulfur sites.
BaCoS 2 has a small orthorhombic distortion quantified by d = a−b a ∼ 0.8% [32] withā = a+b 2 , which is further reduced to d ∼ 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 ∼ 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 min ∼ 0.5%. For the nematic xxxx configuration, the distortion is in good agreement with the experimental values  Table I). We compare C(xyyx) with the two nematic configurations C(xxxx) and C(yyyy) assuming a magnetic order C-AFS-a. Experimental data are taken from Refs. 7, 11, and 32. The phases T0, T7, and T17 refer to (d xxxx min ∼ −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 min ∼ 0.8%).
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 ∼ 1.3 GPa [6-8, 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. 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 8 kbar. Such behaviour is robust either upon choosing the atomic positions 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 ∼ 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 BaCoS 2 under pressure, it would imply that the metal-insulator transition is also accompanied by the recovery of the non-symmorphic symmetry.

C. Wannierization
In order to gain further insight into the electronic structure of BaCoS 2 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 × 4 × 4 k-grid with a doubled in-plane unit cell comprising 4 Co atoms. For both models, the     IV shows the leading hopping processes of the 5-band model restricted to the d xz , d yz 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 I and III, respectively.

III. EFFECTIVE HEISENBERG MODEL
7 smaller than the in-plane NNN one. We finally remark that the orthorhombic distortion has a very weak e↵ect on the inter-layer hopping, which is consistent with the tiny energy di↵erence between C-AFS and G-AFS being insensitive to the distortion, compare, e.g., the energies of T1 and O1 in Tables I and III, respectively. We already mentioned that the largest energy scales that emerge from the ab initio calculations are the magnetic ones separating the lowest-energy AFS configura-tion from the Néel and paramagnetic states. Therefore, even though BaCoS 2 seems not to lie deep inside a Mott insulating regime, we think it is worth discussing qualitatively the spin dynamics in terms of an e↵ective 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 IV suggests the Heisenberg model shown in Fig. 6. It consists of frustrated J 1 J 2 planes [15,[46][47][48][49] coupled to each other by a still frustrating J 3 . The exchange constants satisfy the inequality 2J 2 > |J 3 | + |J 1 |, consistent with the observed columnar magnetic order. Moreover, J 3 forces to deal with a two sites unit cell, highlighted in yellow in the left panel of Fig. 6 where`= 1 and`= 2 refer to 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,

III. EFFECTIVE HEISENBERG MODEL
where, in analogy with the single-layer J 1 J 2 model [15], a 6 = 0 are Ising order parameters associated with the C 4 ! C 2 symmetry breaking. Moreover, we define the in-plane Fourier transforms where R labels the N unit cells in the a b plane,`= 1, 2 the two sites (sublattices) within each unit cell, and n = 1, . . . , L the layer index. Similarly, we introduce the three-dimensional Fourier transforms S`(q, q z ) = X n e iq z n S`, n (q) , as well as the spinor With those definitions, the Hamiltonian reads 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 BaCoS 2 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 IV 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 (±1,0,0)/(0,±1,0) , T ±(1,−1,0)/± (1,1,0) , and T (±1,0,1)/(0,±1,1) , reported in Table IV. 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 > |J 3 |+|J 1 |. 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, where, in analogy with the single-layer J 1 −J 2 model [15], δ a = 0 are Ising order parameters associated with the C 4 → C 2 symmetry breaking. Moreover, we define the in-plane Fourier transforms where R labels the N unit cells in the a − b plane, = 1, 2 the two sites (sublattices) within each unit cell, and n = 1, . . . , L the layer index. Similarly, we introduce the three-dimensional Fourier transforms S (q, q z ) = n e −iqzn S ,n (q) , as well as the spinor With those definitions, the Hamiltonian reads H = 1 N nq J 2 (q) S 1,n (q) · S 1,n (−q) + S 2,n (q) · S 2,n (−q) + J 1 γ(q, δ 1 ) S 1,n (q) · S 2,n (−q) + H.c.
The classical ground state corresponds to the modulation wavevector Q = (π, 0, Q z ) ≡ (0, π, Q z ), related to an antiferromagnetic order within each sublattice on each layer, yielding the lowest eigenvalue ofĴ(q, q z ), i.e., The expression of λ Qz shows that inter-layer magnetic coherence sets in only when the two Ising-like order parameters, δ 1 and δ 3 , lock together. Specifically, J 1 J 3 δ 1 δ 3 > 0 stabilises C-AFS, Q z = 0, otherwise G-AFS, Q z = π. We already know that the former is lower in energy, though by only few Kelvins, see Table III. Moreover, an orthorhombic distortion b > a favours AFS-a, which implies J 1 δ 1 + J 3 δ 3 > 0, even though AFS-b is higher by only 20 K according to DFT+U, see O3 in Table III.
If we assume that the magnetic excitations in the AFS state of strained tetragonal BaCo 0.9 Ni 0.1 S 1.9 are representative of those of orthorhombic BaCoS 2 , then the inelastic neutron scattering data of Ref. [14] at 200 K and for in-plane transfer momentum are consistent with J 2 9.3 meV, J 1 + J 3 −2.34 meV, and J 1 δ 1 + J 3 δ 3 0.53 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 δ 1 δ 3 < 0.08 meV 2 for J 3 J 1 . Such small bound, around a quarter of (J 1 δ 1 + J 3 δ 3 ) 2 , suggests that the two order parameters δ 1 and δ 3 are already formed at 80 K below T N 280 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 IV, and the antiferro-orbital order.
To better understand the interplay between the Z 2 (C 4 ) Ising degrees of freedom and the magnetic order at T N , we investigate in more detail the Hamiltonian (3) with δ 1 = δ 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 S i,n (q) = N S n 3,i,n δ q,Q , i = 1, 2 , n = 1, . . . , L , where S = 3/2 is the spin magnitude, n 3,i,n is a unit vector, and the equivalence holding since G = (π, π, 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 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 where n 1,i,n , n 2,i,n and n 3,i,n are orthogonal unit vectors, x † i,n (q) = x i,n (−q) and p † i,n (q) = p i,n (−q) are conjugate variables, i.e., The three terms of the Hamiltonian (3) thus read, at leading order in quantum fluctuations, i.e., in the harmonic approximation, n,q γ(q) X 1,n (q) · X 2,n (−q) + H.c. , where E 0 = 2N L S(S + 1) J 2 (Q), γ(q) = γ(q, δ = 0), and X i,n (q) = n 1,i,n x i,n (q) + n 2,i,n p i,n (q − Q) .
We note that H 2 does not depend on the choice of 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 In the transformed basis, is diagonal, and the spin-wave energy dispersion is The free energy in perturbation theory can be written as F = F , where F is of -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 = 0, 2, 4, . . . . Given the evolution operator in imaginary time, where H a (τ ), a = 1, 3, evolves with the Hamiltonian H 2 , the second order correction to the free energy is readily found to be n J 2 1 n 3,1,n · n 3,2,n 2 + J 2 3 n 3,1,n · n 3,2,n+1 2 , where with ω λ = 2πλT , λ ∈ Z, bosonic Matsubara frequencies.
Even without explicitly evaluating Ξ 2 , we can conclude that the free-energy gain at second order in H 1 + H 3 is maximised by n 3,1,n · n 3,2,m = ±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 2 1 J 2 3 that reads n n 3,1,n · n 3,2,n n 3,1,n · n 3,2,n+1 + n 3,1,n−1 · n 3,2,n , where We remark that, despite ω 2 (q) vanishes linearly at q = 0, Q, both Ξ 2 (T ) and Ξ 4 (T ) are non-singular. The fourth order correction F 4 in Eq. (7) has a twofold effect: it forces n 3,1,n · 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 n 3,1,n = n 3 and n 3,2,n = σ n 3 , where n 3 is an arbitrary unit vector reflecting the spin SU (2) symmetry, and σ = ±1 is associated with the global C 4 → 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 σ n (R) = n 3,1,n (R) · n 3,2,n (R) and σ n+1/2 (R) = n 3,1,n (R)·n 3,2,n+1 (R), Eqs (6) and (7) 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] We expect [15] an Ising transition to occur at a critical temperature T c , below which σ n (R) = m 1 , σ n+1/2 (R) = 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 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 ⊥ < I 1 between layers. Hereafter, we take for simplicity J 1 = J 3 , thus I 1 = I 3 ≡ I . We then note that, if we use the above estimates of J 2 9.3 meV and J 1 = J 3 −1.17 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 173 K [47]. That corresponds to I 6.6 meV [50]. The 3D Ising critical temperature T c T N depends on the anisotropy ratio I ⊥ /I [50], being at most 345 K for I ⊥ /I 1 1, and reaching the observed T N 290 [12] at I ⊥ /I 1 0.75, which we find not unrealistic. However, BaCoS 2 remains orthorhombic above T N , which implies that the structural C 4 → 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 IV, it is worth repeating the above discussion assuming from the start that σ n (R) = n 3,1,n (R)·n 3,2,n (R) = 1 (AFS-a), so that σ n+1/2 (R) = n 3,1,n (R) · n 3,2,n+1 (R) = n 3,1,n (R) · n 3,1,n+1 (R) .

It follows that the quadrupolar term (8) becomes
A Q − n dR K(T ) n 3,1,n (R) · n 3,1,n+1 (R) 2 +h(T ) n 3,1,n (R) · n 3,1,n+1 (R) , (9) with K(T ) h(T ) > 0, and provides a direct coupling between adjacent planes able to stabilise the 3D magnetic order. For h(T ) = 0, C-AFS, n 3,1,n · n 3,1,n+1 = 1, and G-AFS, n 3,1,n · 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 caxis, 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 lowtemperature phase of BaCoS 2 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 nonsymmorphic 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 BaCoS 2 , 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 supercon-ductors the disordered phase at high temperature first spontaneously breaks C 4 symmetry when cooling below a critical temperature T nem , thus entering a nematic phase [51]. Only at a lower temperature T N < T 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 BaCoS 2 or even vanishing within experimental accuracy as in FeSe, where T N = 0 is observed only under pressure [56]. Also from a model point of view, our description of BaCoS 2 fits into the modelling of these materials. Indeed, the key role played by several order parameters in BaCoS 2 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 [52,60]. 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 (8). Moreover, our ab initio simulations for BaCoS 2 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 nonsymmorphic symmetry instead of C 4 . Therefore, BaCoS 2 seems to realise a situation where collinear magnetism and orthorhombicity do not imply orbital nematicity, unless for specific crystal structures under pressure. BaCoS 2 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.