Mechanism of basal-plane antiferromagnetism in the spin-orbit driven iridate Ba2IrO4

By ab initio many-body quantum chemistry calculations, we determine the strength of the symmetric anisotropy in the 5$d^5$ j $\approx$ 1/2 layered material Ba$_2$IrO$_4$. While the calculated anisotropic couplings come out in the range of a few meV, orders of magnitude stronger than in analogous 3d transition-metal compounds, the Heisenberg superexchange still defines the largest energy scale. The ab initio results reveal that individual layers of Ba$_2$IrO$_4$ provide a close realization of the quantum spin-1/2 Heisenberg-compass model on the square lattice. We show that the experimentally observed basal-plane antiferromagnetism can be accounted for by including additional interlayer interactions and the associated order-by-disorder quantum-mechanical effects, in analogy to undoped layered cuprates.


I. INTRODUCTION
The few varieties of square-lattice effective spin models are emblematic in modern quantum magnetism and extensively investigated in relation to layered superconducting materials such as the copper oxides [1] and the iron pnictides/chalcogenides [2]. While the dominant magnetic energy scale is set in these systems by the isotropic Heisenberg exchange between nearest-neighbor (NN) [3] and possibly next-NN sites [4], there are many examples where the smaller, anisotropic terms become important too, e.g., for correctly describing the antiferromagnetic (AF) ordering pattern in La 2 CuO 4 [5] or in the cuprate oxychlorides [6,7]. This topic, the role of anisotropic interactions in transition-metal compounds has lately received a new impetus with recent insights into the basic electronic structure of 5d systems such as the 5d 5 iridium oxides. Here, a subtle interplay between spinorbit interactions and sizable electron correlations gives rise to insulating ground states and well protected magnetic moments [8][9][10][11][12][13]. Due to the strong spin-orbit couplings, however, these magnetic moments are best described as effective j ≈ 1/2 entities [9,10,14] and the effective anisotropic exchange parameters are orders of magnitude larger than in 3d transitionmetal compounds. For the square-lattice system Sr 2 IrO 4 , for instance, Dzyaloshinskii-Moriya (DM) interactions as large as one quarter of the NN AF superexchange have been predicted [15,16] while in honeycomb iridates the symmetric Kitaev exchange is believed to be even larger than the Heisenberg interaction [17][18][19][20].
Valuable insights into the role of different superexchange processes in correlated d-metal oxides come from the detailed analysis of extended multiorbital Hubbard-type models. The foundations of superexchange theory were laid as early as the 50's with the work of Anderson, Goodenough, and Kanamori [21]. Standard approaches within this theoretical framework proved to be extremely useful in, e.g., better understanding the origin and relative strength of the anisotropic couplings in layered cuprates [22,23]. In two-dimensional (2D) iridates, on the other hand, much less information is presently available on the magnitude of various electronic-structure parameters that enter the superexchange models. While estimates for these effective electronic-structure parameters are normally based on either density-functional band-structure calculations [15,16,18,24,25] or experiments [10,11,13,17,20], we here rely on many-body quantum chemistry methods to directly obtain an ab initio assessment of both the NN Heisenberg exchange and the anisotropic couplings on the square lattice of Ba 2 IrO 4 . Our study reveals uniaxial symmetric anisotropy that is bond dependent, thus giving rise to quantum compass interaction terms [26] superimposed onto the much stronger (due to the 180 • bond geometry) isotropic Heisenberg exchange. We also show that the resulting Heisenbergcompass model for individual layers of Ba 2 IrO 4 is not sufficient to explain the AF ground-state ordering pattern inferred from recent resonant magnetic scattering measurements, with spins ordered along the [110] direction [12]. To rationalize the latter, we carry out a detailed analysis of the role of interlayer couplings and the associated order-by-disorder phenomena. An extended three-dimensional (3D) spin Hamiltonian based on NN exchange terms as found in the ab initio quantum chemistry calculations and additional farther-neighbor interlayer exchange integrals turns out to provide a realistic starting point to explain the magnetism of Ba 2 IrO 4 .

II. GENERAL CONSIDERATIONS
The magnetically active sites, the Ir 4+ ions, have a 5d 5 valence electron configuration in Ba 2 IrO 4 , which under strong octahedral crystal-field and spin-orbit interactions yields an effective j ≈ 1/2 Kramers-doublet ground state, see Refs. [10,14] and [12,27]. The exchange interactions between such pseudospin entities involve both isotropic Heisenberg and anisotropic terms. For a pair of NN pseudospinsS i andS j , the most general bilinear spin Hamiltonian can be cast in the form where J ij is the isotropic Heisenberg exchange, the vector D ij defines the DM anisotropy, and Γ ij is a symmetric traceless second-rank tensor that describes the symmetric portion of the exchange anisotropy. Depending on various geometrical details and the choice of the reference frame, some elements of the DM vector and/or of the Γ αβ ij tensor may be zero. For the square lattice of corner-sharing IrO 6 octahedra in Ba 2 IrO 4 , the symmetry of each block of two NN octahedra is D 2h , with inversion symmetry at the bridging oxygen site [28]. Given the inversion center, the DM anisotropy vanishes. The remaining symmetries require that in the {xyz} frame, with x along the Ir-Ir link and z orthogonal to the IrO 2 layers, Γ ij is diagonal. The two-site effective spin Hamiltonian for an Ir-Ir link along the x axis can then be written as with Γ zz = −(Γ + Γ ⊥ ) since Γ is traceless. Due to the four-fold z-axis symmetry, we analogously have For D 2h symmetry of the two-octahedra unit, the four lowlying (spin-orbit) states, |Ψ S , |Ψ 1 , |Ψ 2 , and |Ψ 3 , transform according to the A 1g , B 2u , B 1u , and A 1u irreducible representations, respectively [19]. As discussed in the following, this symmetry analysis is useful in determining the nature of each of the low-lying many-body states in the quantum chemistry calculations.

A. Computational details
The effective magnetic coupling constants are obtained on the basis of multireference configuration-interaction (MRCI) calculations [29] on units of two corner-sharing IrO 6 octahedra. Since it is important to accurately describe the charge distribution at sites in the immediate neighborhood [30][31][32], we also include in the actual cluster the closest 16 Ba ions and the six adjacent IrO 6 octahedra around the reference [Ir 2 O 11 ] fragment, see Fig. 1 and also Refs. 19,[33][34][35] To make the whole analysis tractable, we however replaced the six Ir 4+ d 5 NN's by closed-shell Pt 4+ d 6 ions, a usual procedure in quantum chemistry investigations on d-metal systems [19,[33][34][35][36][37]. The extended solid-state surroundings were modeled as a large array of point charges fitted to reproduce the crystal Madelung field in the cluster region. We used the crystal structure reported by Okabe et al. [28] All calculations were performed with the MOLPRO quantum chemistry software [38]. Energy-consistent relativistic pseudopotentials from the standard MOLPRO library were used for Ir [39] and Ba [40]. The valence orbitals at the central Ir sites were described by basis sets of quadruple-zeta quality supplemented with two f polarization functions [39] while for the ligand bridging the two magnetically active Ir ions we applied quintuple-zeta valence basis sets and four d polarization functions [41]. The other O's at the two central octahedra were modeled by triple-zeta valence basis sets [41]. For the additional ligands coordinating the six adjacent 5d sites we used minimal atomic-natural-orbital basis functions [42]. At those adjacent 5d sites we applied triple-zeta valence basis sets [39].
Multiconfiguration reference wave functions were first generated by complete-active-space self-consistent-field (CASSCF) calculations [29]. The active space is here given by five electrons and three (t 2g ) orbitals at each of the two magnetically active Ir sites. The orbitals were optimized for an average of the lowest nine singlet and the nine triplet states arising from such an active space. All these states entered the spin-orbit calculations, both at the CASSCF and MRCI levels. In the MRCI treatment, single and double excitations from the six Ir t 2g orbitals and the 2p shell of the bridging ligand site are taken into account. Similar strategies of explicitly dealing only with selected groups of localized ligand orbitals were adopted in earlier studies on both 3d [43][44][45][46] and 5d [19,[33][34][35] compounds, with results in good agreement with the experiment [33,34,[44][45][46]. To separate the metal 5d and O 2p valence orbitals into different groups, we used the orbital localization module available in MOLPRO. The MRCI was performed for each spin multiplicity, singlet or triplet, as a nine-root calculation.
To obtain information on the magnitude of the direct ex-change, we additionally carried out single-configuration restricted open-shell Hartree-Fock (ROHF) calculations [29]. The latter were performed as frozen-orbital calculations, i.e., we used the orbitals obtained by CASSCF (see above), without further optimization. The spin-orbit treatment was carried out according to the procedure described in Ref. [47]. To determine the nature of each spin-orbit state we explicitly compute with MOLPRO the dipole and quadrupole transition matrix elements among those four low-lying states describing the magnetic spectrum of two corner-sharing octahedra, see Table I and the next subsection. Standard selection rules and the nonzero dipole and quadrupole matrix elements in the quantum chemistry outputs then clearly indicate which state is which, see also the analysis and discussion in Ref. [19].

B. Ab initio results
Of the 36 spin-orbit states that are obtained in the ab initio calculations, the low-lying four are listed in Table I [48]. These four states are further mapped onto the eigenvalues of the effective spin Hamiltonian in (2). Energy splittings and the associated effective magnetic couplings are provided at three levels of approximation: single-configuration ROHF (HF+SOC), CASSCF (CAS+SOC), and MRCI (CI+SOC). It is seen that at all levels of theory two of the triplet components, Ψ 1 and Ψ 2 , are degenerate [49]. Given the tetragonal distortions in Ba 2 IrO 4 , with out-of-plane (z-axis) Ir-O bonds significantly stretched as compared to the in-plane (x/y) bonds [28], this degeneracy is somewhat surprising. Using Eqs. (4), this means that two of the diagonal couplings of Γ are equal, Γ zz = Γ ⊥ , which further implies Γ = −2Γ ⊥ . The interaction terms in (2) and (3) can then be rewritten as whereJ ≡ J+Γ ⊥ andΓ ≡ −3Γ ⊥ . Quantum chemistry results forJ andΓ are provided on the lowest line in Table I. The value computed for the HeisenbergJ within the ROHF approximation, −12 meV (see Table I), is sizable and close to the results computed in square-lattice 3d 9 Cu oxides (see, e.g., Ref. 44). It accounts only for direct exchange, since no (intersite) excitations are allowed. In contrast to the ROHFJ, the anisotropicΓ is AF by ROHF.
With correlated wave functions, CASSCF and MRCI, the singlet Ψ S becomes the ground state, well below the triplet components Ψ 1 , Ψ 2 , and Ψ 3 . This shows that the largest energy scale is here defined by the isotropic Heisenberg ex-changeJ (J > 0). In the CASSCF approximation, only intersite d-d excitationsà la Anderson [21] are accounted for, i.e., polar t 6 2g -t 4 2g configurations. Again, the CAS+SOC J, 37.5 meV, is very similar to the CASSCF J's in layered 3d 9 cuprates [44,50,51]. It is seen in Table I that the configuration-interaction treatment, which now includes as well t 5 2g e 1 g -t 4 2g and O 2p to Ir 5d charge-transfer virtual states, enhancesJ by about 70% as compared to the TABLE I. Energy splittings for the four lowest spin-orbit states of two NN IrO6 octahedra and the corresponding effective coupling constants in Ba214, at different levels of approximation (all in meV).
CAS+SOC value, somewhat less spectacular than the ratio between the configuration-interaction and CASSCF J's in layered cuprates. In the latter compounds, this ratio is 3 to 4 [44,52]. If we include in the MRCI treatment only the six Ir t 2g orbitals,J is 49.1 meV (not shown in Table I). The difference between the latter number and the CAS+SOC value given in Table I is indicative of the role of excitation processes via the Ir 5d e g levels. The further increase from 49.1 to 65 meV is due to excitations that additionaly involve the bridging O 2p orbitals. The data in Table I also show that the correlation treatment very much enlarges the symmetric anisotropic couplingΓ , from 0.4 by ROHF to 3.4 meV by MRCI.

IV. COMPARISON TO EFFECTIVE SUPEREXCHANGE MODELS
For the Mott-like insulating regime occurring in the iridates [8][9][10], an effective superexchange model can be in a first approximation set up by considering the leading excited configurations with two holes at the same Ir site. With corner-sharing octahedra and straight Ir-O-Ir bonds along the x axis, the intersite d-d hopping takes place via both in-plane p y and outof-plane p z π-type O orbitals. The relevant effective hopping integrals are t 1 = (t π pd ) 2 /|ǫ xy d − ǫ y p | for the in-plane, xy pair of NN Ir t 2g functions and t 2 = (t π pd ) 2 /|ǫ xz d − ǫ z p | for the outof-plane, xz t 2g functions. ǫ y/z p and ǫ xy/xz d = ǫ 1/2 are here crystal-field split energy levels while the p-d π-type hopping amplitude t π pd is assumed to be the same for both chanels. For tetragonal distortions, ǫ 1 = ǫ 2 , ǫ y p = ǫ z p and therefore t 1 and t 2 may acquire quite different values. A hole hopping between NN Ir ions is then described by the Hamiltonian where d † imσ (d imσ ) is the creation (annihilation) operator of a hole with spin σ in the orbital d xy for m = 1 and d xz for m = 2 at site i. For a bond along the y axis, p y is replaced by p x , d xz by d yz , ǫ y p = ǫ x p , ǫ 3 = ǫ yz d = ǫ xz d = ǫ 2 , and the hopping Hamiltonian in (6) has the same form.
The interaction of two holes in the t 2g subshell is described by Hund's coupling J H and the Coulomb repulsion integrals U mm ′ ≃ U − 2J H , if m = m ′ , and U mm = U . While the isotropic exchange is related to second-order processes that concern transitions between the lowest spin-orbit Kramers doublets, i.e., J ∼ t 2 1/2 /U , the symmetric anisotropy is entirely determined by third-order processes that involve excited Kramers doublets, i.e., is dependent on t 2 1/2 J H /U 2 . The lowest Kramers doublet wave functions as well as those for the excited Kramers doublets are here parametrized as in Ref. [10], with the angle θ given by By collecting the second-and third-order processes in this effective superexchange model, we arrive at the pseudospin Hamiltonian in (2), with Here η = J H /U and γ = − η U cos 2 θ[(t 1 − t 2 ) 2 sin 2 θ + t 2 1 sin 2 θ + 1 2 t 2 2 cos 2 θ]. Now, for Γ zz = Γ ⊥ , the model described by (5) displays uniaxial compass-like anisotropy [26]. That is obviously the case for perfect, cubic octahedra with ∆ = 0, t 1 = t 2 = t, and cos θ c = √ 2 sin θ c = 2/3. In the cubic limit we further have from Eqs. (8): J c = (16/9)t 2 /U + γ c , γ c = −(4η/9)t 2 /U , Γ c = −γ c , and Γ c ⊥ = Γ c zz = (−2η/3)t 2 /U − γ c . For tetragonal distortions as found in Ba 2 IrO 4 [28], Γ ⊥ = Γ zz implies that (t 2 /t 1 ) 2 = 2 tan 2 θ. As measure of how large the departure from the cubic limit is we can take the ratio between the tetragonal t 2g splitting ∆ and the strength of the spin-orbit coupling λ. The quantum chemistry calculations yield ∆ = 65 meV in Ba 2 IrO 4 , see the discussion in Ref. [53], in agreement with estimates based on experiment [27]. A direct estimate of the spin-orbit coupling can be also obtained from the splitting of the j = 1/2 and j = 3/2 t 5 2g states for idealized cubic octahedra. It turns out that for perfect octahedra λ = 0.47 eV [34,53], close to values of 0.4-0.5 eV earlier derived from electron spin resonance and optical measurements on 5d 5 ions [54][55][56][57]. The ratio ∆/λ is therefore rather small, ≈0. 15.
Estimates for the parameters that enter the effective superexchange model can be most easily obtained in the cubic limit. Using Eqs. (8) we find thatΓ /J ≈ (3/8)η. The CI+SOC values of Table I,Γ = 3.4 andJ = 65 meV, then lead to η ≈ 0.14 and 4t 2 /U ≈ 149 meV. Interestingly, estimates of the hopping integral t from calculations based on density-functional theory (DFT) are t DFT ≈ 260 meV, while the on-site Coulomb repulsion comes out from constrained calculations in the random phase approximation (RPA) as U RPA ≈ 1.65 eV [25]. The ratio 4t 2 DFT /U RPA is therefore ≈164 meV, close to the result derived on the basis of the CI+SOC effective couplings listed in Table I. On the other hand, the η parameter extracted from the periodic DFT calculations [25] is η DFT ≈ 0.08, much smaller than the above value of 0.14. Using the latter value for η, η DFT ≈ 0.08, an estimate for the symmetric anisotropic couplingΓ = 3 8 ηJ would be significantly smaller than the quantum chemistry result.

V. GROUND STATE PHASE DIAGRAM
Having established the strength of the dominant in-plane exchange interactions and anisotropies, we now turn to the nature of the magnetic ground state of Ba 2 IrO 4 , focusing first on a single square-lattice IrO 2 layer. In the classical limit, the compass-Heisenberg model defined by Eqs. (5) has an accidental SO(2) ground-state degeneracy, with spins pointing along any direction in the basal xy-plane [26,58,59]. This degeneracy is eventually lifted via thermal [60][61][62] or quantum [23,[62][63][64] order-by-disorder effects, whereby harmonic spin wave fluctuations select the states with spins pointing either along the x or y axis. This is however in sharp contrast to experiments which below ∼240 K show basal-plane AF order with magnetic moments along the [110] direction [12]. It indicates additional anisotropies in the system, large enough to overcome the energy gain from the order-by-disorder mechanism.
The situation is actually analogous to several 3d 9 Cu oxides with the same "214" crystal structure as Ba 2 IrO 4 . It has been shown that in cuprates that particular type of AF order is selected by a subtle interplay between in-plane and interlayer interactions, as discussed in detail in Ref. [23]. Assuming that qualitatively the same 3D mechanism is applicable to Ba 2 IrO 4 , we analyze below the main contributions to the ex- pression of the 3D ground-state energy and derive a generic phase diagram. This exercise provides useful insights into the dependence of the ground-state spin configuration on various interaction parameters in 214 iridates.
It turns out that the most important effects competing with the in-plane NN interactions concern (i) the frustrating nature of the isotropic interlayer exchange and (ii) the symmetric part of the anisotropic exchange between layers. To show this we proceed by parametrizing the global spin direction in each basal plane by an angle φ n , where n is the layer index, and by writing down all relevant energy contributions.
The first contribution is the zero-point energy (per spin) coming from the order-by-disorder mechanism in each individual layer, and ω ± (q) are the two spin wave branches, for which explicit expressions are provided in Appendix A. A numerical analysis of Eq. (9), using the ab initio quantum chemistry values for the in-plane NN effective couplings (see Sec. III), shows that E ZP,2D (φ) is almost identical to the expression with K = 0.86 µeV and E 0 = 56.55 meV. We now turn to the second contribution to the energy, which stems from the interlayer isotropic exchange J out . Despite being the dominant portion of the interlayer interactions, its total contribution to the energy vanishes in the mean-field sense due to geometric frustration in the 214 structure, see Fig. 1. Yet quantum fluctuations driven by J out still give rise to a zeropoint energy contribution where B ≃ 0.032J 2 out /(2J av ) and J av = J + (Γ + Γ ⊥ )/2 [22]. Since B is positive for any sign of J out , this contribution favors collinearity of the staggered magnetization in adjacent layers.
The third contribution to the energy comes from the anisotropic portion of the interlayer couplings. We first note that the antisymmetric DM component vanishes by symmetry since the midpoint of each of these out-of-plane NN Ir-Ir links is an inversion center. The remaining, symmetric portion, can be described by a traceless second-rank tensor Γ out . The structure of the latter is simplified by using the fact that the out-of-plane NN Ir-Ir links are C 2 axes, additionally perpendicular to reflection planes. Adding up the four tensors (related to each other by symmetry) from all four NN bonds above/below the reference layer gives [23] E aniso,3D = −A n sin(φ n + φ n+1 ), where the constant A is fixed by the elements of Γ out . The total energy now reads It can be minimized analytically as described in Appendix B by working it out for a bilayer of Ba214. The resulting phase diagram in the (A/K, B/K) plane is shown in Fig. 2 for positive A (the phase diagram for A < 0 is identical, see Appendix B), and hosts three different phases, two collinear (phases I and II) and one noncollinear (phase III).
In phase I, the staggered magnetizations point along one of the 110 axes and the relative directions between adjacent planes are regularly collinear or anticollinear. In phase III, the AF magnetization prefers one of the 100 axes and the relative directions in adjacent planes are now perpendicular to each other. Finally, in phase II, the relative directions between adjacent planes are again either collinear or anticollinear but the staggered magnetizations in each layer rotate in the basal plane as a function of A/K, see Appendix B. Importantly, the degeneracy is not completely lifted by the above couplings. As explained in Appendix B, all phases have an Ising degree of freedom per layer, which comes from the fact that the energy remains the same if we flip all spins within a given layer. This remaining macroscopic degeneracy may eventually be lifted via higher-order processes or fartherneighbor couplings, see for example the discussion in [22]. The collinear AF structure observed experimentally [12] in Ba 2 IrO 4 can now be naturally explained provided that A and B fall into the broad region of phase I in the phase diagram of Fig. 2 and by taking into account removal of the macroscopic Ising degeneracy by the mechanism mentioned above.
As pointed out by Boseggia et al. [12], the AF component of the ordered momenta in the 214 iridates Sr 2 IrO 4 and Ba 2 IrO 4 is essentially identical -in Sr 2 IrO 4 , the canted AF state is characterized by an AF vector aligned along the 110 direction and a residual FM moment confined to the same basal plane. Staggered rotation of the IrO 6 octahedra as realized in Sr 2 IrO 4 requires the more general single-layer Hamiltonian of Eq. (1) [10,16], with a DM vector along the z axis and a biaxial easy-plane symmetric anisotropy described in our notation by two independent diagonal components Γ > 0 and Γ zz > 0. This model correctly explains the canting angle of the basal-plane AF order [16] but fails in predicting the AF vector alignment along one of the 110 axes. The reason is that the two additional anisotropies, D||z and Γ zz > 0, do not remove the SO(2) basal-plane ground-state degeneracy, at least not in the classical limit. This accidental degeneracy can however again be lifted via the 3D mechanism discussed above, to arrive to an AF ordering pattern similar to that of Ba 2 IrO 4 [12].

VI. CONCLUSIONS
While ab initio quantum chemistry techniques have been earlier used to derive the sign and strength of the symmetric anisotropic (Kitaev) interactions in 5d 5 iridates with edge-sharing IrO 6 octahedra [19], we here employ the same methodology to clarify the signs and magnitude of the symmetric anisotropic couplings for corner-sharing octahedra in the square-lattice compound Ba 2 IrO 4 . The ab initio results reveal effective uniaxial anisotropy, although the actual sym-metry of each of the in-plane Ir-Ir links is lower than D 4h . The anisotropic effective coupling constants are as large as 3.5 meV, comparable in strength with the anisotropic Kitaev exchange in honeycomb Na 2 IrO 3 [19]. However, in contrast to Na 2 IrO 3 , the largest energy scale is here defined by the Heisenberg J, with J ≈ 65 meV. The latter value agrees with estimates based on resonant inelastic x-ray scattering measurements on 214 iridates [11]. Given the uniaxial structure of the exchange coupling tensor, the relevant in-plane (pseudo)spin model is a Heisenberg-compass type of model. Yet to understand the experimentally determined AF ordering pattern, with spins along the [110] direction [12], interlayer interactions must be included in the effective Hamiltonian.
Further investigations are now carried out in our group for quantifying the strength of Dzyaloshinskii-Moriya couplings for the closely related 214 compound Sr 2 IrO 4 , displaying bent Ir-O-Ir links. Another interesting issue is the dependence of the in-plane anisotropic couplings, their signs in particular, on pressure [65] and strain [66], in both Ba 2 IrO 4 and Sr 2 IrO 4 .

VII. ACKNOWLEDGEMENTS
L. H. acknowledges financial support from the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG).

Appendix B: Energy minimization for a bilayer
The ground-state magnetic energy of the layered system can be written as a sum over bilayer contributions (per spin and per layer): where the angles φ 1 and φ 2 define orientations (say, with respect to the x axis) in two adjacent planes and φ ± = φ 1 ± φ 2 . We note that both K and B are positive. In the subsequent discussion, the coupling A is chosen positive as well by taking into account the fact that for A < 0 the simultaneous change of signs, φ 1 → −φ 1 and φ 2 → −φ 2 , retains the expression for E(φ 1 , φ 2 ) invariant.
Minimizing E(φ 1 , φ 2 ) we find four possible extrema solutions for φ 1 and φ 2 and the respective energies (n and m are integers): which is possible if B > K; with the requirement A < 4K; which is possible if B < K; sin φ which may occur in the parameter region B < K, A < K 1+B/K 2 . Comparison of the energies of the four possible ground-state configurations shows that three of them, from (B1) to (B3), occur in different domains of the A-B parameter space. In the region B > K and A > 4K, the most stable is the configuration (B1) with φ (1) 1 = π 4 + n π 2 and φ 1 − mπ, which means that the spins (staggered magnetizations) are along one of 110 axes and in two adjacent planes the spin alignment is either collinear or anticollinear. Next, in the region with B > K, 0 < A < 4K, the second configuration (B2) with φ 1 − mπ is realized. Here, the collinear/anticollinear alignment in successive layers still persists. However, the preferred direction is specified by A/4K. In the region with B < K, A > 4K, the third configuration (B3) with φ 1 − π 2 − mπ is the most stable, which corresponds to having the magnetization along one of the 100 axes with two directions in successive layers being rotated by 90 • . Finally, for B < K and A < 4K, the fourth solution (B4) has the highest energy and two of the other configurations, i.e., (B2) and (B3), compete to give the phase boundary depicted in Fig. 2.
Ising degrees of freedom -It is clear that the above classical minima of a Ba214 bilayer are also the minima of the infinite system. In all phases, however, there is still an Ising degree of freedom per layer, which is not fixed by the couplings considered here. In phase I for example, we may flip the directions of all spins in any plane, since the energy is the same for both collinear and anticollinear relative orientations between adjacent planes. The eventual removal of this remaining macroscopic Ising degree of freedom must originate from higher-order processes or farther-neighbor couplings [22].