Frustration relief and reorientation transition in the kagome-like dolerophanite Cu 2 OSO 4

We present a theoretical study of dolerophanite Cu 2 OSO 4 , a layered kagome-like spin-12 magnetic insulator that can be described either as a system of chains coupled through dimers or as a kagome lattice where every third spin is replaced by a ferromagnetic spin dimer. Building on insights from ab initio calculations, classical numerical minimizations, and semiclassical expansions, we arrive at a minimal microscopic description that accounts for the experimental data reported so far, including the nature of the magnetic order, the reported spin length, and the observed anisotropy. The latter arises by a peculiar competition between the antisymmetric (Dzyaloshinskii–Moriya) and the symmetric part of the exchange anisotropy, which gives rise to a two-step re-orientation process involving two successive continuous phase transitions. Our work uncovers mechanisms stabilizing canted ferrimagnetic order in kagome systems, and highlights strong magnetic anisotropy in the presence of dissimilar magnetic orbitals on crystallographically nonequivalent Cu sites. We also show how these anisotropy terms affect the spin-wave spectrum and how they can be tracked experimentally.


I. Introduction
The search for quantum spin liquids remains at the forefront of condensed matter and quantum magnetism research for many decades [1][2][3][4][5][6][7][8][9], since P. W. Anderson's first proposal of the resonating valence bond idea [10][11][12][13].One of the main classes of candidate materials that have been explored intensively over the years are transition-metal compounds with dominant isotropic interactions and strong geometric frustration, most notably the 3D pyrochlores and quasi 2D kagome magnets [1,2,14,15].In such materials, spins can evade magnetic ordering down to very low temperatures, despite their strong exchange interactions, due to the proliferation of infinite competing states at low energy scales.In conjunction with low dimensionality and low spin quantum number S , this competition magnifies the effects of quantum-mechanical fluctuations, and opens the possibility for unconventional classical and nonclassical states, including spin nematics, valence bond crystals, and gapped and gapless spin liquids [1][2][3][4][6][7][8][9].
While the above ingredients -frustration, low dimensionality and low spin -are present in many of the candidate materials that have been identified and characterized over the years, the broader consensus is that these same ingredients are also responsible for the high sensitivity to perturbations that are inevitably present in these compounds, such as longer-range interactions, magnetic anisotropy and various forms of randomness.Understanding the role of these perturbations is therefore crucial for explaining experimental data and for developing a useful phenomenology for the search of unconventional states.
Here we present a theoretical study of dolerophanite Cu 2 OSO 4 , a highly-frustrated, quasi-2D, kagome-like magnet, where magnetic anisotropy plays a key role.Cu 2 OSO 4 has been known as a mineral since the early 1960s [16].Its structural characterization has been presented in Refs.[17][18][19], but its magnetic behaviour was explored only recently [20][21][22].The geometry of a given layer in Cu 2 OSO 4 can be thought of as a system of chains joined through dimers or, alternatively, as a kagome lattice where every third spin is replaced by a dimer (Fig. 1).The material consists of two symmetry-nonequivalent Cu 2+ ions denoted by Cu1 and Cu2; the ions form the kagome lattice where every third vertex is occupied by a Cu2 dimer, and Cu1 ions form chains along the b-axis.As it turns out, the two spin- 1  2 moments that make up a given dimer are so strongly ferromagnetically coupled that it would also be possible to consider an effective mixed-spin kagome model where one third of the sites is represented by S = 1 moments.
At first glance, experimental data for Cu 2 OSO 4 suggest only a weak magnetic anisotropy revealed by the Curie-Weiss temperatures (Θ CW ) of −71 K, −75 K and −70 K measured along the a, b, and c * directions, respectively [22].While the negative Θ CW values would suggest predominantly antiferromagnetic couplings, the magnetic ground state of Cu 2 OSO 4 determined by neutron diffraction below T N = 20 K is best described as a canted ferrimagnet [22].Spins adopt a uniform coplanar configuration with a 120 • -like structure in the ab-plane.The ferromagnetically coupled Cu2 spins point along the baxis, resulting in an uncompensated total magnetic moment of ∼ 0.23(3)µ B /Cu along this direction [21,22].This value is remarkably close to the classical value of µ B /4 expected for the coplanar 120 • spin arrangement in a mixed-spin kagome magnet.Moreover, the ferromagnetic alignment of adjacent kagome layers implies that the total moments from different layers do not compensate each other.
Reflecting on these preliminary reports raises some important questions, which we seek to answer in our work.First and foremost, we wish to construct a minimal microscopic model that accounts for the experimental data [20][21][22], and elucidate to what extent frustration of the kagome layers can account for the physics of Cu 2 OSO 4 .We also wish to examine the origin of the uncompensated magnetic moment, understand why it is pinned to a specific direction in the crystal structure, and what determines the preferred plane for the coplanar spins.
Ultimately, in addressing these questions we also aim to understand what sets Cu 2 OSO 4 apart from the conventional kagome physics.At first glance, such a deviation is not unexpected, due to the structural peculiarities of Cu 2 OSO 4 , in particular, the presence of four Cu sites in the unit cell, instead of three, and the absence of a three-fold symmetry and the inequivalence between the nearest-neighbour (NN) exchange paths.However, a closer analysis reveals that, despite these differences, the Heisenberg model of Cu 2 OSO 4 remains highly frustrated, with an infinite ground state manifold that is similar, e.g., to that of the kagome francisites [23], and undistorted kagome magnets.Yet, the observation of a robust lowtemperature ordered state with an almost classical value of the spin lengths [21,22] indicates minimal quantum fluctuations despite the frustration.
This brings us to the next main focus of this work, which is to uncover the key role played by the magnetic anisotropy in this compound.As we show below, the lifting of the infinite ground state degeneracy of the isotropic model, and the selection of the observed, almost classical order, arises from a non-trivial interplay of the antisymmetric part of the exchange anisotropy, i.e., the Dzyaloshinskii-Moriya (DM) interactions [24,25], and the symmetric (traceless) part, which we shall denote by T in this work.According to our Density Functional Theory (DFT) calculations, the relevant DM vectors feature significant components along the crystallographic c * axis, and much weaker components in the ab plane.This hierarchy explains the selection of the observed uniform coplanar state and the pinning of the spins on the ab plane, but fails to reproduce the direction of the uncompensated moment in that plane, as the in-plane DM components select the a and not the (observed) b axis.The inclusion of the symmetric part of the exchange anisotropy resolves this final puzzle and leads to a reorientation of the total moment along b, provided this anisotropy has large enough strength to override the effect of the in-plane DM components, and that it has the appropriate sign.Our DFT calculations for one of the most relevant bonds confirm explicitly that such an anisotropy is present in Cu 2 OSO 4 .We further show that the re-orientation of the total moment proceeds via a two-step process involving an intermediate phase, and that the spin gap depends very sensitively on the interplay between the in-plane DM components and the symmetric part of the exchange anisotropy.
Having identified the right minimal model, our final aim is to carry out self-consistent checks but also make predictions for further experiments, which, in turn, will allow for more quantitative estimates of the microscopic coupling parameters in this complex, multi-sublattice system.To that end, we shall present a study of the response of Cu 2 OSO 4 under a magnetic field along the three crystallographic directions, and uncover the presence of appreciable transverse magnetization which can be probed by magnetic torque measurements.We shall also study the effect of quantum fluctuations and confirm that the reduction of the spin length is indeed very small, in agreement with experiment.Finally, we shall present the magnon band structure from linear spin-wave calculations.The remainder of the paper is organized as follows.In Sec.II we highlight the main aspects of the crystal structure of Cu 2 OSO 4 , including the global symmetries, the different coordination of the two types of Cu ions, and their lattice topology.In Sec.III we discuss details of DFT calculations (Sec.III A) and the resulting electronic structure (Sec.III B), which sheds light into the nature of the magnetic orbitals of the two types of Cu sites in Cu 2 OSO 4 .In Sec.IV we present our study of the effect of isotropic Heisenberg interactions, the predicted hierarchy of energy scales from DFT (Sec.IV A), the resulting picture for the infinite classical ground state manifold and its close similarity to the 2D KHAFM (Sec.IV B), and the effect of the interlayer couplings (Sec.IV C).In Sec.V we present what happens when we include the microscopic DM interactions.This includes an analysis of the constraints imposed by symmetry on the DM vectors (Sec.V A), their numerical values from DFT (Sec.V B), the resulting picture for the structure of the selected ground states (Sec.V C), and a comparison to experiments.In Sec.VI we incorporate the symmetric part of the exchange anisotropy tensors, discuss their constraints from symmetry (Sec.VI A), identify particular elements of these tensors that can turn the total moment along the observed direction (Sec.VI B), and present DFT results on one of the relevant bonds that corroborate this picture (Sec.VI C).In Sec.VI D we show, more generally, that the rotation of the total moment proceeds via a two-step reorientation transition, and discuss the evolution of the spin gap through this transition.In Sec.VII we study some additional topics and make predictions for further experiments: The response of Cu 2 OSO 4 under a magnetic field along the three crystallographic directions (Sec.VII A), the effect of quantum fluctuations at the quadratic level (Sec, VIII), the associated reduction of the spin length (Sec.VIII 1) and the magnon band structure (Sec.VIII 2).In Sec.IX we give our conclusions and a broader perspective of our work.Finally, we include two appendices (App.A and B) which contain technical details and auxiliary information.Thus, each layer can be thought of as either a system of chains joined through dimers or as a kagome lattice with every third spin replaced by a dimer.The overall spin lattice can be described in terms of the underlying monoclinic Bravais lattice, plus a basis of four atoms, two Cu1 and two Cu2, see Fig. 1.

II. Main structural and symmetry aspects
Cu 2 OSO 4 crystallizes in a monoclinic crystal system belonging to the space group C2/m [16], where the monoclinic c-axis forms an angle of β = 122.34• with the a-axis.The vector c * is perpendicular to the ab plane, as shown below Apart from primitive translations, the C2/m symmetry group includes: i) inversion center in the middle of the Cu2 dimers; ii) inversion center in the middle of the J ⊥,2 bond; iii) twofold rotation axis C 2b passing through the middle of the Cu2 dimers; iv) screw axis along the Cu1 chain corresponding to a non-primitive translation by b/2 followed by a two-fold rotation C 2b around the b axis; and v) ac * mirror plane that goes through the midpoints of NN Cu1 (or J 3 ) bonds and includes the Cu2 dimers.The experimental magnetic structure does not break any of these symmetries (it belongs to the irrep Γ 1 of the C2/m space group at zero momentum [22]).
Another important aspect of the Cu 2 OSO 4 structure concerns the local coordination of the two types of Cu 2+ ions, which affects the choice of the magnetic orbital.The Cu1 atoms are surrounded by four oxygen atoms with the Cu-O distances of 2 × 1.882 Å (Cu1-O3) and 2 × 2.070 Å (Cu1-O4), see Fig. 2(d).Two further oxygen atoms are located 2.526 Å away from copper and belong to the SO 4 groups.Therefore, Cu1 has the 4+2 oxygen coordination with the square plaquette of four nearest oxygen atoms obtained by an axial elongation of the CuO 6 octahedron.By contrast, Cu2 features two short Cu2-O2 distances of 1.907 Å followed by three further oxygens at 2.000 Å (Cu2-O3) and 2 × 2.155 Å (Cu2-O4), resulting in the 2+3 trigonal-bipyramidal coordination.This structural configuration gives rise to a rather unusual coupling regime that we discuss in the following.

A. Methods
Density-functional (DFT) band-structure calculations were performed in the FPLO [26] and VASP codes [27,28] using the Perdew-Burke-Ernzerhof (PBE) flavor of the exchangecorrelation potential [29].Correlation effects in the Cu 3d shell were taken into account on the mean-field DFT+U level with the on-site Coulomb repulsion U d = 8.5 eV (FPLO) or 9.5 eV (VASP), Hund's coupling J d =1 eV, and atomic limit for the double-counting correction [30,31].Note that U d is applied to the Cu 3d orbitals, hence there is usually a difference in the optimal U d value depending on the DFT code and its basis set [32].Experimental structural parameters from Ref. [17] were used in all calculations.The various microscopic spin interactions, including the Heisenberg exchange parameters presented below in Sec.IV), the DM vectors (Sec.V) and the symmetric exchange anisotropy T 3 (Sec.VI), have been obtained by a mapping procedure [33,34] from total energies of magnetically ordered states.These energies were converged on a k mesh with 64 points in the first Brillouin zone for the supercells with 64 atoms (double the crystallographic unit cell of Cu 2 OSO 4 ).
Each interaction parameter is derived from total energies of four spin configurations as follows, where S = 1 2 and ↑↑, ↑↓, ↓↑, ↓↓ denote spin arrangements on the interacting atoms 1 and 2 while spins on all other atoms are kept fixed to a globally ferromagnetic configuration.The αβ component of the interaction tensor J is obtained in a similar way by setting the spin on the atom 1 along α and the spin on the atom 2 along β, whereas the spins on all other atoms are kept along γ perpendicular to both α and β.When α β, the resulting component is a combination of the DM interaction (D γ ) and off-diagonal symmetric anisotropy (T αβ ) that are separated as (3) where ϵ αβγ is the Levi-Civita symbol.The FPLO values are for the Heisenberg exchange parameters only, because noncollinear spin configurations are not implemented in this DFT code.Additional data related to the DFT results of this work can be found in Ref. [35].

B. Electronic structure
The uncorrelated (PBE) band structure of Cu 2 OSO 4 shows predominantly oxygen states below −3 eV followed by Cu 3d bands that extend to slightly above the Fermi level, see Fig. 2(a).The band structure is metallic because correlations in the Cu 3d shell have not been included in the calculation.The states near the Fermi level show contributions from both the d x 2 −y 2 and d 3z 2 −r 2 orbitals, although only one of them should be eventually half-filled and magnetic in the 3d 9 Cu 2+ ion.Choosing the x 1 and y 1 local coordinate axes along the shorter Cu1-O3 and Cu1-O4 bonds within the square plaquette results in a separation of the two orbitals, with the d x 2 −y 2 states lying at higher energies than the d 3z 2 −r 2 ones, see Fig. 2(b).Therefore, d x 2 −y 2 is the likely magnetic orbital for Cu1.A similar separation is possible in the case of Cu2, but here the local z 2 axis has to be directed along the two shortest Cu2-O2 bonds.The d 3z 2 −r 2 states are then higher in energy than the d x 2 −y 2 ones, and d 3z 2 −r 2 is the likely magnetic orbital for Cu2, see Fig. 2(c).This analysis is verified by DFT+U calculations that produce an insulating solution.Correlations stabilize the orbital state with the unpaired electron on the d x 2 −y 2 orbital for Cu1 and on the d 3z 2 −r 2 orbital for Cu2.This is similar to the hightemperature phase of another kagome-like mineral, volborthite [36,37], but different from the majority of other Cu-based kagome magnets, including francisite [23] where both Cu sites feature d x 2 −y 2 as the magnetic orbital.I.Both sets of energies are given relative to the reference energy E 0 that corresponds to the nonmagnetic state.

IV. Isotropic model
We begin our analysis with an investigation of the isotropic Heisenberg spin Hamiltonian where J i j denotes the exchange parameter between spin-1/2 operators S i and S j at Cu sites i and j, respectively.Exchange couplings for all Cu-Cu pairs with the distances of less than 7 Å were calculated using supercells doubled either in the ab plane or along the c direction.Two complementary DFT codes produced consistent results except for the coupling J 1 that is twice larger in VASP compared to FPLO.Both sets of parameters will be considered in the following.The leading Heisenberg exchange couplings are listed in Table I and are also indicated in Fig. 1.Further exchange couplings can be neglected in the minimal model because they are below 10 K. Fig. 3 shows that the six leading couplings considered in this work allow a good description of total energies for different spin configurations in Cu 2 OSO 4 .
The coupling on the Cu2 dimers, which is denoted by J d , sets the dominant energy scale and is ferromagnetic.The NN coupling on the Cu1 chains, which is denoted by J 3 , is antiferromagnetic and gives the second strongest energy scale.The couplings J 1 and J 2 between Cu1 and Cu2 ions are also antiferromagnetic and compete with J 3 .Finally, there are two types of interlayer interactions, J ⊥ and J ⊥2 , which couple two Cu1 atoms to a single Cu2 atom, and two Cu2 atoms to a single Cu2 atom, respectively.The predominance of these couplings over other possible superexchange pathways between the kagome layers can be explained by the presence of bridging SO 4 tetrahedra, as shown in panel (d) of Fig. 2. The interlayer couplings are weaker than J 2 and J 3 , but overall significant.Here, J ⊥2 is non-frustrated and at first glance would couple the layers antiferromagnetically.On the other hand, J ⊥ forms triangular loops with J 3 in the same way as J 1 does.The coupling J ⊥ would not be satisfied in the scenario of simple antiferromagnetic order imposed by J ⊥2 .In the experimental magnetic structure, the J ⊥2 − J ⊥2 − J 3 triangles end up having the same type of canted order as the J 1 − J 1 − J 3 triangles, so J ⊥2 may act as an additional driving force of the canted order, see detailed analysis in Sec.IV B.
Let us now compare the experimentally measured values of the Curie-Weiss temperature Θ CW to the values obtained from the DFT Heisenberg couplings of Table I and the expression (see discussion in App.B) where the sign convention of Θ CW corresponds to the Curie-Weiss approximation for the susceptibility χ = C T −Θ CW , and C is the Curie constant.The resulting values are provided in the first four rows of the second column of Table II, where, for comparison, we give the values corresponding to VASP and FPLO, as well as the values without (2D) and with (3D) the interlayer couplings (J ⊥ and J ⊥,2 ) included.The experimentally reported values are −68 K for a powder sample [21], and −71 K, −75 K and −70 K along the a, b, and c * axes, respectively, for a single crystal [22]).The comparison shows that: i) including the interlayer couplings gives a more satisfactory agreement, and ii) the VASP values give a better agreement compared to FPLO.
B. The isotropic model of a single layer: Highly-frustrated magnetism & similarity to the KHAFM The isotropic limit of a single layer of Cu 2 OSO 4 features an infinite number of classical ground states, which are closely related to the ones of the ideal KHAFM [38][39][40][41].To see this, we follow a cluster minimization method (see, e.g., Ref. [42]) and examine the four-site unit cell of the structure.

The building block of classical ground states
Let us first rewrite the total Heisenberg energy of a single layer, H iso,1 , as a sum over contributions from unit cells, (7) where J d is replaced by J d /2 because this bond is shared by two cells.We have, where As is turns out, in the minimum energy configuration, the spins S 3 and S 4 point along the same direction and the energy reduces to that of an 'isosceles Heisenberg triangle' with competing couplings J 3 and J 1 + J 2 , where S is the classical spin length.The minimum energy state takes the schematic form (11) or, more explicitly, where e 1 and e 2 are any pair of orthogonal axes, and the angle θ can be found by minimizing the energy The values of θ predicted from FPLO and VASP are 54.5 • and 43.8 • , respectively.From these, the former is closer to the value of 60 • realized in the 'equilateral triangle' limit of J 1 + J 2 = J 3 , which seems to be the relevant region for Cu 2 OSO 4 [22].The total magnetic moment per Cu can be obtained by Taking S = 1/2 and assuming g = 2, we find m ≃ 0.21µ B (FPLO) and 0.14µ B (VASP) per Cu.The former is only slightly lower than the value 0.23(3)µ B found experimentally [21].So, at the level of the isotropic model of a single layer, the FPLO values give a better agreement for the the length of the total moment.This comparison is premature however, since, as we show below (see also fourth column of Table II), these numbers show significant variation when including the interlayer couplings and/or the anisotropies.

Tilling the 'ABC' state onto the lattice
Having found the ground states of one building block of the Hamiltonian, we can now proceed to tile these on the full lattice.The resulting ground state manifold consists of infinite ground states, including an infinite subset of coplanar states and an infinite subset of non-coplanar states.Some representative members of this manifold are shown in Fig. 4.
Coplanar States -The uniform coplanar state of Fig. 4 (a) is the simplest way to tile the ABC state of Eq. ( 12) to the full lattice, and is also the one that has been observed experimentally [22].This state has an SO(3) order parameter space, as there are two angles needed to fix the plane of A, B and C, and an additional angle to fix the direction of (say) the A sublattice within that plane.
However, this is not the only coplanar state.Indeed, starting from the uniform state of Fig. 4 (a), one can generate the coplanar state shown, e.g., in Fig. 4 (b) by choosing one of the L horizontal rows of the lattice and swapping B↔C along that row.Such an operation does not affect the energetics on the particular row of triangles, or on any other triangles of the lattice, and is therefore a ground state.Counting all possible sequences of swaps on horizontal chains gives 2 L such coplanar states.
Similar swaps can be performed on non-horizontal chains as well.However, one cannot perform a sequence of swaps along chains of different direction, so the total degeneracy of coplanar states is 3 × 2 L (and not 2 3L ).The uniform state is a special member of this sub-extensive family of coplanar states.
Non-coplanar states -The above operation of swapping B↔C along a given row corresponds to a π-rotation of the spins along that row around the direction of the A sublattice.This procedure can be generalized to rotations by any angle ϕ.The spins along the given row will rotate around A from B and C to some new directions B' and C', which are out-of-the original ABC plane, see Fig. 4 (c).Given that all neighbours of that row of sites point along A, a rotation around A preserves the local structure and the relative angles necessitated by Eqs.(12)(13) and does not alter the energy.In total, we can perform such a rotation by an angle ϕ r on the r-th row of the lattice, and these angles are in general different, leading to an SO(2) degeneracy for each row (i.e., SO(2) L in total).Similar rotations can be performed along the non-horizontal chains as well.
Besides the above non-coplanar ground states, one can envisage other ground states which could arise by performing a sequence of rotations along open or closed paths that are different from chains, provided we do not affect the energy of any triangle.While our numerical minimizations show evidence for such more complex ground states, it is not clear whether they belong to the above manifold of non-coplanar states that are generated by performing a sequence of rotations along chains.
Summarizing the physics of the isotropic model on a single layer of Cu 2 OSO 4 , we have found a classical ground state manifold with infinite coplanar and infinite non-coplanar ground states, which can all be built from the local 3-sublattice building block of Eq. ( 12).The uniform coplanar state found experimentally is a special member of this manifold.While this state could potentially be the one that is selected by quantummechanical fluctuations, we anticipate (based on what happens, e.g., in the case of the KHAFM [39,41,[43][44][45][46]) that the respective order-by-disorder energy scale is much smaller compared to the interlayer couplings or the anisotropy.And given that  given by FPLO code (see Table I), and H iso + H DM (VASP, 3D) corresponds to the actual multilayer system, with isotropic and DM couplings (including J ⊥ , J ⊥,2 , D ⊥ and D ′ ⊥ ), and numerical values given by the VASP code.The second column gives the Curie-Weiss temperature Θ CW , which is given by Eq. ( 6) for the 2nd, 4th and 6th rows, and by the same equation but with J ⊥ and J ⊥,2 set to zero for the 1st, 3rd and 5th rows.For the last two rows, the three values given correspond to Θ aa CW , Θ bb CW and Θ c * c * CW , which are given by Eq. (B17) in the appendix.The vector m denotes the total moment per Cu site and m its length.For the meaning of the angles θ, α, ψ and ξ see Eq. ( 29)- (31).the latter may lift the classical ground state degeneracy already at the mean-field level, we shall disregard the order by disorder physics of the pure isotropic model, and proceed to investigate the effect of the interlayer couplings J ⊥ and J ⊥,2 and the anisotropy.
C. Effect of interlayer couplings Referring back to the right panel of Fig. 1, one sees that J ⊥,2 connects one Cu2 dimer from one layer to a Cu2 dimer on the next layer.So the main effect of this interaction is to fix the relative orientation of Cu2 dimers on successive layers, leading to an alternating A, -A, A, • • • structure of the Cu2 dimers along the c * axis.However, J ⊥,2 does not lift the infinite degeneracy discussed above, and the system still features infinite coplanar and infinite non-coplanar ground states.For example, we can start from one of the coplanar states in a single layer and then 'copy' its time-reversed version on the next layer, and continue in an alternating fashion with the next layers, so that we satisfy J ⊥,2 .However, one can still swap B↔C (or -B↔-C) along any Cu1 chain of any layer, without affecting the energy, leading again to a sub-extensive number of coplanar ground states.The situation for non-coplanar states can be explained in a similar way.
Importantly, none of the members of the classical ground state manifold matches the one found experimentally, because even when the in-plane configuration is uniform, the relative orientation is AFM across the layers due to J ⊥,2 .Hence, this interaction does not seem to be relevant for the explanation of the uncompensated moment in Cu 2 OSO 4 .

Effect of J ⊥ alone
Let us now discuss the effect of J ⊥ .According to the right panel of Fig. 1, J ⊥ connects a NN pair of Cu1 sites on one layer to a Cu2 site on the next layer.Effectively then, this coupling plays a role similar to that of J 1 and J 2 , competing with J 3 .At the mean-field level, this coupling gives rise to a parallel arrangement of the Cu2 dimers across different layers, in contrast to the antiparallel arrangement favoured by J ⊥,2 .Additionally, J ⊥ preserves the ABC structure of Eq. ( 12) of the building block of the lattice, but renormalizes the angle θ of Eq. ( 13) to which arises from Eq. ( 13) by replacing reflecting the fact that J ⊥ plays a similar role with J 1 and J 2 , as mentioned above.The ab initio parameters of Table I give θ ≃ 39.6 • (FPLO) and 25.5 • (VASP), which are further away from the 60 • of the equilateral triangle case, compared to the ones without J ⊥ , see comparison in fourth column of Table II.In turn, these numbers give a total moment of m ≃ 0.11µ B (FPLO) and 0.05µ B (VASP) based on Eq. ( 14), which are significantly reduced compared to the experimental values.Finally, using the same arguments as in the case of J ⊥,2 above, one finds that J ⊥ too fails to select a unique ground state, and the classical ground state manifold contains again infinite coplanar and infinite non-coplanar states.The only differences to the case of J ⊥,2 then are: i) the renormalization of the angle θ, and ii) the fact that now the Cu2 dimers align ferromagnetically across the layers.
Importantly, the experimentally reported configuration is one of the members of the ground state manifold, although there is an appreciable deviation in the value of θ.This suggests that either DFT overestimates the value of J ⊥ , and/or that the anisotropic couplings can remedy this problem (see below).

Combined effect of J ⊥ and J ⊥,2
Based on the above, the two interlayer couplings compete with each other, since J ⊥ favours FM alignment of Cu2 dimers across the layers, whereas J ⊥,2 favours AFM alignment.When both couplings are taken into account, the numerical minimization of the classical energy on large finite-size clusters, based on the ab initio values of Table I for J ⊥ and J ⊥,2 , delivers a unique ground state.This is a coplanar state with eight spin sublattices in total, and with successive layers arranged antiparallel to each other.The configuration of each layer is uniform with four sublattices.Hence the combined effect of J ⊥ and J ⊥,2 is to lift the infinite degeneracy completely, but the selected state is not the one found experimentally.
We therefore conclude that the isotropic model of Cu 2 OSO 4 is not enough to account for the experimentally reported state and that further couplings seem to play a qualitative role.In what follows we then turn to the investigation of the anisotropic interactions in Cu 2 OSO 4 , starting from Dzyaloshinskii-Moriya anisotropy [24,25].

V. Dzyaloshinskii-Moriya Interactions
We now turn to the effect of Dzyaloshinskii-Moriya (DM) interactions, which enter the spin Hamiltonian in the general form where D i j is the microscopic DM vector on the bond (i j).
The presence of DM interactions in Cu 2 OSO 4 has been discussed by Takahashi et al, who have also estimated an effective (coarse-grained) DM parameter |D| ≃ 7 K, from the value of the magnon gap measured by antiferromagnetic resonance [20].Below, we show that this value underestimates significantly the actual size of the DM vectors because, i) not all of the DM components contribute to the magnon gap, and ii) the latter is affected by the symmetric part of the exchange anisotropy.

A. Constraints from symmetry
Let us first obtain the structure of the DM interactions using symmetry arguments.First of all, the inversion centers on the middle of the J d and J ⊥,2 bonds necessitate that the DM vectors on these bonds vanish identically.This leaves the DM vectors on the remaining bonds, J 1 , J 2 , J 3 and J ⊥ .
Let us consider the local geometry of two neighbouring building blocks shown in Fig. 5 (a), which shows three consecutive Cu1 sites, labeled as Cu11, Cu14, Cu14s, and four Cu2 sites (residing on the two adjacent Cu2 dimers), labeled as Cu22, Cu22', Cu22s and Cu22's.Now, the ac * mirror plane that crosses through the midpoint of J 3 bonds and contains the Cu2 dimers, maps (S a , S b , S c * ) → (−S a , S b , −S c * ) in spin space.Hence, the DM vectors connected by this operation are related to each other by a π-rotation.This implies, in particular, that the DM vectors on the J 3 bonds lie on the ac * plane.
Next, let us consider the constraints imposed by the screw axis along b, which consists of a C 2b rotation in spinorbit space followed by a translation by b/2.This operation maps Cu11→Cu14, Cu14→Cu14s, Cu22→Cu22s, and Cu22'→Cu22's in real space, and (S a , S b , S c * ) → (−S a , S b , −S c * ) in spin space.Hence, the DM vectors on any two bonds that are connected by this operation must also be related by a π-rotation around b.
Altogether, the symmetries impose the following interrelations between the various DM vectors in Fig. 5 The above relations lead to a sign structure of the DM vectors shown in Fig. 5 (b).
Similarly, the DM vectors, D ⊥ and D ′ ⊥ , connecting two NN Cu1 sites with a Cu2 site on the next layer [see Fig. 1) (b)] are related to each other by the ac * mirror plane, and therefore

Insights from cluster decomposition
Similarly to what we did in Eq. ( 8), we can re-write the total DM energy as a sum over contributions from the building blocks of Fig. 5 (a), namely with As we shall see in Sec.V C, the classical ground state of the lattice model (with Heisenberg and DM couplings included) is a uniform, 3-sublattice state, similar to the state of Eq. ( 11), with all Cu2 sites parallel to each other, due to the strong J d coupling.We can use this result and replace S 4 → S 3 in ( 22) to obtain the simpler expression where One can also incorporate the contribution from the interlayer DM couplings by simply replacing these expressions with The above tell us that, due to the strong FM coupling

B. Insights from DFT calculations
The DM vectors calculated using the VASP code are provided in the last column of Table I.The results show large DM vectors for all of the J 1 , J 2 , and J 3 bonds.The fact that these vectors are comparable in size to the Heisenberg terms seems quite unusual and can be traced back to the close competition between the d x 2 −y 2 and d 3z 2 −r 2 magnetic orbitals.Microscopically, DM interactions can be understood as the effect of electron hoppings between the nonmagnetic (empty/filled) and magnetic (half-filled) d-orbitals in the presence of spin-orbit coupling [25,47].By fitting the energy bands between −1 eV and +0.5 eV with a tight-binding model containing two orbitals (d x 2 −y 2 and d 3z 2 −r 2 ) on each of the Cu sites, we find that such hoppings are comparable in size to hoppings between the magnetic orbitals.For example, we find t 2 = −0.136eV for the hopping between d x 2 −y 2 of Cu1 and d 3z 2 −r 2 of Cu2 (magneticto-magnetic) vs. t ′ 2 = 0.131 eV for the hopping between d x 2 −y 2 of Cu1 and d x 2 −y 2 of Cu2 (magnetic-to-nonmagnetic). This would still lead to |D|/J < 1, but the Heisenberg term can be further lowered by a ferromagnetic contribution in the closeto-90 • coupling geometry [48], and DM interactions can be eventually favored over isotropic exchange.
Finally, besides the DM vectors on the J 1 , J 2 and J 3 bonds, there also exists an interlayer DM interaction that lives on the J ⊥ bonds, which is however much weaker (see Table I) than the remaining anisotropies.
There are two main insights from the DFT results.Either way, these vectors have a very large component along c * (±68 K), like D 3 .Hence, we anticipate that the main effect of the DM interactions is to align the spins almost entirely on the ab plane, which is consistent with experiment.Another effect is to increase the angle 2θ between NN Cu1 spins beyond the value given by Eqs. ( 13) and ( 15), as the DM favours a 90 • orientation.
Turning to the much weaker in-plane DM components, their role is to fix the orientation of the spin structure within the ab plane and additionally give rise to a weak canting of the spins out of that plane.The specifics of these effects will be examined in Sec.V C.

Results from cluster decompositions and unconstrained minimizations
Let us return to the cluster decomposition of Eq. ( 21), and examine the ground state of H DM, .A numerical minimization of this four-site cluster delivers a ground state which is similar to that of Eq. ( 11) but with one crucial difference: the spins 3 and 4 of the Cu2 dimer show a weak misalignment.Due to the latter, this state cannot be tiled on the lattice and therefore the decomposition of Eq. ( 21) does not help to identify the exact ground state of the full lattice model.We then turn to the second smallest cluster, with two tetrahedral blocks, (28) where we have also indicated the direction of the DM vectors on all bonds.The minimum energy configuration of this cluster shows no misalignment between the two spins of the Cu2 dimer (i.e., S 4 = S 3 ), and in addition S 5 = S 1 and S 6 = S 2 .This state can be tiled on the lattice, leading to a uniform, threesublattice state, as shown in Fig. 5 (b).This state coincides with the state we find independently via unconstrained numerical minimizations on large finite-size clusters with periodic boundary conditions.The 3-sublattices A, B and C of the above state can be depicted as (29) and their directions can be written as follows where This configuration is qualitatively similar to the one of Eq. ( 12).The angle 2θ is the angle between the two Cu1 spins [as in Eq. ( 12)], the angle α specifies the plane of the Cu1 spins (S 1 and S 2 ), the angle ψ fixes the overall orientation within that plane, and the angle ξ accounts for a small canting of the Cu2 spins (S 3 and S 4 ) out of that plane.The state of Eq. ( 30) reduces to that of (12) when ξ = ψ = α = 0, v 1 → e 1 and v 2 → e 2 .
Let us discuss some further qualitative features of this state: i) The state is almost coplanar.The plane of the Cu1 spins lies very close to the ab-plane, and the Cu2 sites are tilted slightly away from that plane.Namely, the angles α and ξ are both very small, see 6th and 8th column in Table II.This can be traced back to the dominance of the c * components of D 3 , D eff 1,3 and D eff 2,3 .Furthermore, our unconstrained numerical minimizations of large clusters shows that this result applies to both the single-layer and multiple-layer cases, with only small quantitative changes when the interlayer interactions are included (energy is lowered slightly, and angles between spins are renormalized).
ii) The angle 2θ between the two Cu1 spins is larger than the value given by Eqs. ( 13) and ( 15), as the DM favours a 90 • orientation, see fifth column in Table II.
iii) The total moment per site is given by (assuming an isotropic g tensor, with g = 2) which reduces to Eq. ( 14) for ξ = 0 and v ′ 2 → e 2 .The VASP values of Table I give a total moment per Cu site m = 0.22µ B , for both the single layer and multilayer cases (see 4th column in Table II), which is very close to the experimental value [21,22].
iv) The direction of the total moment lies in the ac * plane, and predominantly along the a-axis (see 3rd column in Table II), and the angle ψ = π/2.The selection of this direction stems from the in-plane components of the DM vectors.This can be seen more directly by replacing the ansatz of Eq. ( 30) into the expression for the energy of Eq. ( 23), which reveals that the terms involving ψ are all proportional to sin ψ, and and therefore the minimum energy corresponds to ψ = ±π/2.This result disagrees with experimental data, which show that the total moment is along the b axis and not the a axis [21,22].Additionally, this disagrees with the refinement of neutron diffraction data, which delivers a state that belongs to the identity irrep Γ 1 C2/m at zero momentum [22], while the state favored by the in-plane DM components clearly breaks the ac * mirror plane and belongs to the Γ 3 irrep.So the model with Heisenberg and DM couplings included reproduces all experimental data except one, the direction of the uncompensated moment.According to Eq. ( 33), this deficiency cannot be remedied within the above model, unless we abandon the major insight from DFT on the separation of energy scales between the out-of-plane and the in-plane components of D 3 , D eff 1,3 and D eff 2,3 .As this ingredient is crucial for getting the remaining aspects of the observed state right, the origin of the for the direction of the total moment must be sought in some other type of anisotropy, with the symmetric interactions T being the natural candidate.
The energy scale of this additional anisotropy should be at least as strong as the in-plane components of D 3 , D eff 1,3 and D eff 2,3 .To get an idea of this energy scale we have fixed the parameters α, ξ and θ to the values corresponding to the minimum energy configuration, and then tracked the variation of the total energy as we rotate ψ away from ±π/2.This calculation delivers a variational bound to the energy cost, which varies within a bandwidth of only ∼ 8.77 K as we vary ψ between 0 and 2π.Hence, only a small amount of the additional anisotropy is necessary to rotate the moment along the direction found experimentally.

VI. Symmetric anisotropy T
The symmetric (and traceless) part of the anisotropy T enters the spin Hamiltonian in the form where the 3 × 3 matrices T i j have the properties (where α and β are Cartesian components) and Tr(T i j ) = 0.

A. Constraints from symmetry
The screw axis along b dictates that T 3 is uniform along the chains (unlike the DM vector D 3 which is staggered along the chains).Furthermore, the ac * mirror plane dictates that the only nonzero off-diagonal matrix elements of T 3 are T ac * 3 .The same is true for the matrix T d .By contrast, the matrix elements of T 1 and T 2 in the (abc * ) frame are all nonzero, in principle.Furthermore, the screw axis along b dictates that the signs of the matrix elements T ab 1,2 = T ba 1,2 and T c * b 1,2 = T bc * 1,2 alternate from one building block of the structure to the next, as we travel along the Cu1 chains.For the bonds of Fig. 5 (a), we have, for example, and and similarly for T 11,22 ′ ≡ T 2 and T 14,22 ′ ≡ T ′ 2 .

B. Identifying the relevant element of the symmetric exchange tensor
Similarly to what we did in Eqs. ( 8) and ( 21), we can re-write the total energy from T interactions as a sum over contributions from the building blocks of Fig. 5 (a), namely with Anticipating that the strong FM coupling J d will again enforce parallel (or almost parallel) alignment of the Cu2 spins, we can replace S 4 → S 3 in ( 38) to obtain the simpler expression where In principle, one can also add the contribution from the interlayer couplings (T ⊥ , T ⊥,2 , and their mirror images T ′ ⊥ and T ′ ⊥,2 ), but we shall disregard them as they are likely much weaker compared to T 3 , T d , T eff 1,3 and T eff 2,3 .In the following we shall set out to identify the components of the T matrices which favour a rotation of the total moment along the b axis, and as such can compete with the in-plane DM components.We shall do this for the matrices T d and T 3 that correspond to the bonds with the strongest Heisenberg exchange paths, which are also the ones that lead to simple analytical insights.

Effect of T d
To understand the effect of T d , we examine its spectrum.The eigenvalues are given by and the corresponding eigenvectors by It follows that, if T bb d < 0 then the T d interaction favours alignment along the b axis, whereas if T bb d > 0 it favours alignment in the ac * plane.So, the crucial ingredient is the sign of T bb d .

Effect of T 3
By re-writing where T , we see that one can understand the effect of T 3 by examining the spectrum of the 6 × 6 matrix . Now, T 3 has the same form as T d , so its eigenvalues and eigenstates are, respectively, and The desired eigenvalues and eigenvectors of  can then be expressed in terms of λ (3)  j and u (3) j (with j = 1 − 3) as respectively.From these results, it follows that the crucial ingredient is the sign of the matrix element T bb 3 : If T bb 3 > 0, the minimum eigenvalue is λ (3)  2 and the form of the corresponding eigenvector (u (3)  2 , u (3) 2 ) T tells us that the total moment of the two Cu1 spins, S 1 + S 2 , must lie on the ab plane.By contrast, if T bb 3 < 0, the minimum eigenvalue is λ (3)  1 and the form of the corresponding eigenvector (u (3)  1 , u (3) 1 ) T tells us that S 1 + S 2 will be aligned with the b axis, and will therefore compete with the in-plane components of the DM vectors.So, similarly to the case of T d , the most relevant aspect of T 3 is the sign of T bb d .
C. T 3 matrix from DFT & effect on the ground state Having identified the relevant matrix elements of the symmetric anisotropy, we have carried out DFT calculations based on the VASP code, in order to check if such matrix elements are in principle present in Cu 2 OSO 4 , ii) have the right sign, and iii) are strong enough to compete with the in-plane DM components.The results for the J 3 bond, (in unit of K), show that T bb 3 is negative and of comparable size with the energy scale required to turn the moment along b, as found experimentally [21,22].
To check this explicitly, we have performed unconstrained numerical minimizations of a single layer of Cu 2 OSO 4 using large clusters with periodic boundary conditions, and including the Heisenberg, DM and T 3 couplings delivered by VASP.The results confirm that a T bb 3 coupling of the order of −5 K is enough to rotate the total moment along the b axis, as found experimentally [21,22].We have also found that the angle ξ which describes the canting of the Cu2 spins out of the plane of the Cu1 spins, reduces to zero and the configuration becomes fully coplanar.Additionally, the angle α that describes the tilt of the spin plane away from the ab plane decreases to 3-5 • .Such a small misalignment of the spin plane from the ab-plane can be easily missed experimentally.Finally, the total moment per Cu site is 0.21µ B and 0.9µ B for one and multiple layers, respectively, see Table II.It is likely that the remaining T couplings will renormalize these numbers slightly, but the agreement to experiment [21,22] is already very satisfactory.

D. Nature of the re-orientation process & the spin gap
To uncover the nature of the re-orientation transition we examine the evolution of the classical ground state of the Hamiltonian as we vary the rescaling parameter λ from zero to one.The results are summarized in Fig. 6.We find that the re-orientation proceeds via two continuous phase transitions, involving an intermediate phase (indicated by shading), where the total moment rotates from predominantly along the a-axis to b-axis.
The boundary between the intermediate phase and the largeλ phase involves the spontaneous breaking of the ac * mirror Spin gap -The main curves of Fig. 6 show the evolution of the spin gap through the two-step re-orientation process, as obtained from linear spin-wave theory (see further discussion on magnon dispersions in Sec.VIII below).First of all, the value of the spin gap depends very sensitively on λ across the two-step re-orientation transition.In particular, the spin gap drops to zero at the two boundaries of the intermediate phase.The closing of the gap at these boundaries is related to the spontaneous breaking of the symmetries discussed above.
Away from the intermediate phase, the spin gap reaches about 19 K and 33 K at λ = 0 and λ = 1, respectively.These values are set by the energy scales of the in-plane DM components of D 3 , D eff 1,3 and D eff 1,3 , and by that of T 3 respectively.The out-of-plane DM components alone do not contribute to the spin gap because, in the absence of the in-plane DM components and T 3 , the system has a continuous U(1) symmetry and the states of Fig. 6 break this symmetry spontaneously.This is shown explicitly in Fig. 8 (a) that we discuss in Sec.VIII below.
Comparison to experiment -Let us now compare the calculated values of the spin gap to the one extracted experimentally from antiferromagnetic resonance data, which is about 1.14 K [20].This value is much smaller compared to the calculated gap at λ = 1.However, one should keep in mind that this value has been measured for a very weak mode that becomes visible below 5.3 K only, i.e., at temperatures much lower than T N .The nature of this mode must be clarified in future experiments.It is also worth noting that the above Hamiltonian includes the symmetric anisotropy on the J 3 bond only, and that the T couplings on the remaining bonds may compete with the effect of T bb 3 , thereby effectively reducing the spin gap and/or bringing Cu 2 OSO 4 closer to the boundary with the intermediate phase, where the spin gap vanishes.
A related comment is in order here, which may connect the puzzle of the smallness of the observed spin gap with another puzzle that concerns the reported magnetization data for H ∥ a, see Fig. 5(b) of Ref. [21].According to these data, m a shows a small jump at zero field, from about −0.0125 to +0.0125 Bohr magnetons per Cu site.One plausible explanation can be the presence of an accidental misalignment of the field towards the b axis.Comparing the zero-field jump in m a for H ∥ a to the corresponding jump in m b for H ∥ b, this misalignment would amount to only 3 • , which is a plausible explanation.However, an alternative scenario, which could also resolve the puzzle of the smallness of the spin gap, is that the system may actually be slightly inside the intermediate phase.A definite answer to this requires further dedicated experimental studies.

VII. Other aspects and predictions for further experiments A. Magnetization process
To compare with reported magnetization curves we have also studied the effect of an applied external field H along the three directions a, b and c * .The Zeeman coupling to the field takes the usual form, where we have assumed an isotropic g tensor with g = 2.
To find the evolution of the classical ground state of Cu 2 OSO 4 for different field directions and field strengths, we have first performed unconstrained minimizations on large finite-size clusters with periodic boundary conditions.As it turns out, the state remains uniform for all field strengths.This allows us to decompose H Z in terms of contributions from the clusters shown in (28), namely with where we use the site-labelling of ( 28), and the prefactor of 1/2 in the first term inside the square bracket takes care of double-counting.A numerical minimization of this 6-site cluster delivers, as before, S 4 = S 3 , S 5 = S 1 and S 6 = S 2 , which simplifies Eq. ( 49) to Before we present the numerical results, let us discuss some symmetry aspects.In the absence of a field, the system is invariant under the ac * mirror symmetry M ac * , and also under the time reversal operation T .Now, when the field is along the a or along the c * axis, the system is invariant under the combined operation M ac * × T , which maps The classical ground states are found to break this symmetry spontaneously at a finite field H * , below which the system shows the following nonzero order parameters: S a 1 −S a 2 , S b 1 +S b 2 , S c * 1 − S c * 2 and S b 3 .By contrast, in the presence of the field along the b axis, the ac * mirror plane survives, which means that the states related by the transformation have the same energy.In principle, the classical ground states could break this symmetry spontaneously, but this is not supported by our numerics (the associated order parameters, which, in this case, are S a 1 + S a 2 , S b 1 − S b 2 , S c * 1 + S c * 2 , S a 3 and S c * 3 , all vanish irrespecive of the field strength).
Our numerical minimization data shown in Fig. 7 give insights for the evolution of the classical ground state of a single layer of Cu 2 OSO 4 as a function of an external magnetic field along a (panels a-c), c * (d)-(f), and b (g)-(i).The results for the multilayer case are qualitatively the same.The first row of panels shows the evolution of the three components of the total magnetic moment per Cu, m.The second row shows the combinations of Eq. ( 52) and (53), which, as discussed above, play the role of order parameters.Finally, the last row shows the combination of spin components which are driven explicitly by the field and their evolution in a wider field range.Let us discuss the main features of Fig. 7: 1) Unlike the case of a field along b (panel g), a field along a (panel a) and/or c * (panel d), can induce appreciable transverse magnetization components, which can be detected experimentally by torque measurements.Such measurements can help to extract a more quantitative estimate of the relevant microscopic couplings.
2) Fields along a (panel b) and c * (panel e) induce a phase transition at a characteristic field strength, H * a and H * c , respectively.Below these fields, the symmetry M ac * × T is broken spontaneously, and the order parameters of Eq. ( 52) become nonzero.By contrast, for fields along b (panel h), there is no such phase transition, and the mirror symmetry M ac * remains intact, irrespective of the field strength.
3) At the fields H * a and H * c , the system has almost the same magnetic structure as in zero field, but with the total moment turned along the field direction.Here, the applied field works against the effect of T 3 (which favours alignment along the b axis, as discussed above), and also partly against the in-plane components of the DM vectors (which favour alignment of the total moment in the ac * plane, and predominantly along the a axis).Along with the presence of multiple bonds, this can explain the large values obtained for H * a and H * c .4) The combinations of spin components that are driven explicitly by the field (data shown in the last row of panels in Fig. 7) show that: i) the saturation field H sat is beyond 200 T for the microscopic parameters obtained from VASP, and ii) that the way to saturation in panels (c) and (i) involve an abrupt re-orientation of the Cu1 spins at some very high fields H * *

VIII. Impact of quantum fluctuations
The picture presented so far has been entirely classical.The natural next step is to consider the effect of quantum fluctuations and confirm that their impact is minimal in Cu 2 OSO 4 , as observed experimentally [21,22].As explained earlier, this aspect cannot be accounted for at the level of the isotropic model description of Cu 2 OSO 4 , which remains highly frustrated despite the structural peculiarities of this material.The presence of symmetric and antisymmetric exchange anisotropy is crucial for resolving this puzzle as well, as these interactions open an appreciable spin gap (as shown already in Fig. 6), which in turn alleviates the impact of quantum fluctuations.To show this explicitly we have performed a standard linear spin-wave expansion [50,51] around the ground state of a single layer of Cu 2 OSO 4 and for various anisotropic terms included in the Hamiltonian.The multiple-layer problem does not feature any qualitative changes for the in-plane dispersion of the magnons and the onsite spin length reduction.

Spin length reduction
Let us take the site labelling convention of Fig. 1, and denote by w j (j=1-4) the direction of the four spins of the magnetic unit cell in the classical ground state; Incidentally, recall that these directions are different for T 3 = 0 and for T 3 0, as the corresponding ground states are different.The spin length of the j-th spin in the renormalized ground state is then given by ⟨S w j ⟩ = S − ⟨n j ⟩, where S = 1/2 and n j is the number operator for the spin flips (magnons).
For the model without T 3 (and all DM components included) we find DM only : while for the model with both DM and T 3 included we find DM + T 3 : As expected by symmetry, the two Cu2 sites (j=3 and 4) have equal spin length and the same is true for the two Cu1 sites (j=1 and 2).Moreover, the effect of quantum fluctuations on the Cu2 sites is less severe due to the strong FM J d exchange.
Most importantly, the reduction of the spin length is minimal (at most 12% of the classical length of 1/2).This is consistent with experiment [21,22], and corroborates our earlier assertion that the impact of quantum fluctuations in Cu 2 OSO 4 is heavily mitigated due to the appreciable anisotropy gap, giving further confidence to the above classical description.

Magnon excitations & dynamical structure factor
We will now examine the spin-wave expansion more closely and discuss the magnon excitation spectrum.The uniform states discussed above have a four-site unit cell, and we therefore expect four magnon bands, emerging from the hybridization and delocalization of the four local spin flips on the unit cell.
Figure 8 shows the resulting magnon band structure of a single layer of Cu 2 OSO 4 .The dispersions shown correspond to special cuts on the ab plane, inside the first Brillouin zone (BZ), which is an almost symmetric hexagon (see inset).In panel (a), the spin Hamiltonian includes Heisenberg plus the c * -components of the DM vectors, whereas in panel (b) it includes all DM components as well as the interaction matrix T 3 .In panel (a), the lowest mode is gapless at the Γ point because the Hamitlonian with only the c * -components of the DM vectors included has a U(1) symmetry around the c * axis, and this symmetry is broken spontaneously by the classical ground state around which we expand.In panel (b), this mode is gapped out due to the in-plane DM components and T 3 , as discussed above.Apart from this qualitative difference, the remaining shape and overall positioning of the magnon bands is similar in the two panels, which reflects the fact that the intermediate-and high-energy parts of the excitation spectrum is not affected much by the anisotropies.
On top of the magnon bands of Fig. 8 we also show by symbols (filled circles) the zero-temperature dynamical structure factor intensity (more specifically, its trace), where M(Q, t) is the time-evolved Fourier transform of the magnetization evaluated at wavevector Q and time t.The size of the symbols (and colour) is scaled relatively to the maximum intensity in each separate panel.The majority of the scattering weight is carried by the three lowest magnon branches.The fourth magnon band carries almost no weight at all, and features a much weaker dispersion with the bandwidth of the order of a few K.

Further analysis of the magnon bands
To gain some understanding of the overall magnon band structure we proceed to a more in-depth analysis of the magnon modes.To that end, we take again the labelling of the unit cell of ( 7), and use, for each separate spin, the local quantization axis pertaining to the direction of that spin in the reference classical state around which we wish to expand.With this choice, the reference state is the tensor product of over all unit cells, whre S 34 is the total spin quantum number of the Cu2 dimer and M 34 is its projection along the local quantization axis.At the mean-field level, the four magnon states and their corresponding energies ∆E (measured from the reference state) are where h 3 = h 4 is the magnitude of the local mean field exerted on a Cu2 site from its Cu1 neighbours, and, similarly, h 1 = h 2 is the magnitude of the local mean field exerted on a Cu1 site from all its neighbours.Specifically, where ⟨S j ⟩ (j=1-4) are the spin vectors in the classical reference state around which we expand.For the isotropic model, the magnitudes of these fields are given by whereas, in the presence of anisotropies, their values can be obtained numerically.
We have presented a comprehensive study of dolerophanite Cu 2 OSO 4 , a layered kagome-like magnet where every third site of the kagome structure is occupied by a strongly ferromagnetic spin dimer.We have arrived at a minimal microscopic model that uncovers the origin of most experimental data reported so far [21,22], including the intralayer and interlayer magnetic structure, the reported total magnetic moment, and the observed anisotropy in the direction of the uncompensated moment.
Our results show that the origin of the canted order in Cu 2 OSO 4 comes from antiferromagnetic NN interactions.This fact sets this material apart from francisite and related compounds that feature canted order arising from ferromagnetic NN interactions frustrated by an antiferromagnetic coupling between second neighbors [23,52].In Cu 2+ compounds, the sign of short-range exchange interactions is usually determined by the Cu-O-Cu angles in the crystal structure.The relevant values for Cu 2 OSO 4 are 90.5/104.7 • (J 1 ), 117.6 • (J 2 ), and 114.2 • (J 3 ), all being similar to francisite and closer to 90 • than, e.g., in herbertsmithite (∼ 119 • ) with its antiferromagnetic NN couplings.The antiferromagnetic nature of J 1 , J 2 , and J 3 is probably reinforced by the peculiar orbital configuration with the different symmetries of the magnetic orbital on the Cu1 and Cu2 sites, see Fig. 2(d).
Given the antiferromagnetic couplings in the kagome plane, the canted order can be broadly understood as a deformation of the 120 • state expected in KHAFM in the presence of perturbations, such as anisotropy and weak lattice distortions.This order may [53] or may not [54] feature an uncompensated total moment depending on the exact symmetry of the material and the nature of perturbations therein.Such canted order in Cu 2 OSO 4 is necessarily ferrimagnetic because one out of three sites is occupied by effective spin-1 moments of the Cu2 dimer.An interesting finding of our work is that the θ-angle close to 60 • in Cu 2 OSO 4 , i.e., an almost undistorted 120 • configuration, is a combined effect of geometrical frustration and anisotropy.Heisenberg exchange alone should result in much lower values of θ as a result of the sizable deformation of the kagome lattice.However, large anisotropy terms restore a close-to-120 • spin configuration that allows a larger energy gain from the out-of-plane DM components.
The large DM anisotropy in Cu 2 OSO 4 is quite unusual and striking.It may be traced back to the close competition between d x 2 −y 2 and d 3z 2 −r 2 as possible magnetic orbitals and the eventual stabilization of different orbitals on different Cu sites.A similar situation is realized in the skyrmion material Cu 2 OSeO 3 where the sizable ratio of |D/J| ≃ 0.6 was found for one of the couplings [31].The simultaneous presence of d x 2 −y 2 and d 3z 2 −r 2 magnetic orbitals may thus be a promising route for realizing large-DM interaction regime experimentally.At any rate, a complete understanding of the precise mechanism behind the large DM interactions is still lacking, and warrants further dedicated investigations, e.g., via superexchange expansions that incorporate the precise low-symmetry crystal field of the two Cu sites, the various hopping amplitudes and the effect of the spin-orbit coupling.
The presence of strong anisotropy is consistent with the fact that the observed magnetic moment is almost classical [21,22].This agrees with our linear spin-wave calculations, which demonstrate that the reduction of the spin length is very small in the presence of anisotropy.The latter ingredient is crucial, because the isotropic model description of Cu 2 OSO 4 delivers an infinite ground state degeneracy at the classical level, which in turn is expected to give rise to large quantum fluctuations.
Another key result is that the observed direction of the uncompensated moment results from a peculiar competition between the antisymmetric (Dzyaloshinskii-Moriya) and the symmetric part of the exchange anisotropy, which gives rise to a two-step re-orientation process involving two continuous phase transitions.
One of the remaining open issues is the size of the magnon gap.Previous studies reported a small value of 1.14 K using antiferromagnetic resonance [20], but the exact nature of the mode probed in this experiment remains to be understood.To facilitate this understanding, we have presented predictions for further experiments, which can allow for a more quantitative refining of the microscopic coupling parameters.These include: i) the behaviour of Cu 2 OSO 4 under magnetic fields along different crystallographic directions and the corresponding prediction of appreciable transverse magnetization that can be measured by torque experiments, ii) the prediction of the magnon excitation spectrum (and the associated dynamical spin structure factor) which can be measured by inelastic neutron scattering experiments, and iii) the prediction for the presence of a flat two-particle mode involving the |S 34 = 1, M 34 = −1⟩ member of the triplet on the Cu2 dimers, with excitation energy ∆E = 2h 3 , which can be measured by Raman scattering.
We hope these predictions and the minimal model uncovered in this study will provide the basic framework and necessary guidelines for further theoretical and experimental works in Cu 2 OSO 4 and related materials.the direction of the (thermal averaged) field-induced moments ⟨S i ⟩ are not necessarily the same.The total field exerted at site i, including the external field, is given by where the minus sign in the first term on the right hand side follows from our convention for the Zeeman energy of site i, which is −gµ B S i • H.For spin-1/2, the thermal average ⟨S i ⟩ in the mean-field decoupled problem, at temperature T , equals where h i = |h i | is the magnitude of the local field and ĥi is its direction.In the high-T limit we can replace Using Eq. (B7) then gives where we have defined r ≡ 1/(4k B T ).We note in passing that, for general spin S , this expression must be replaced with r = S (S +1)/(3k B T ).With varying i and µ, Eq. (B9) is a system of coupled equations which can be written in the compact form where, using the abc * frame, and Λ iµ, jν = δ i j δ µν + rK µν i j , or, more compactly, Λ = 1 + rK.The solution of Eq. (B10) can be expanded in powers of r (or, equivalently, in powers of 1/T ): Suppose now that we wish to obtain the diagonal elements of the susceptibility tensor.We apply the field in a fixed direction, say µ, and calculate the total spin (per site) along the same direction: ⟨S µ tot ⟩/N = gµ B r − r 2 where in the second step we used Eq.(B2).Note that the DM interactions do not contribute to the Curie-Weiss temperature because of the antisymmetry of the tensor A i j .Applying Eq. (B16) to the full model that we consider in the main text (which includes Heisenberg exchange, DM vectors and the symmetric anisotropy T 3 ) gives

FIG. 1 .
FIG. 1.(a) (a, b) view of a single layer of Cu 2 OSO 4 .Cu1 and Cu2 atoms are shown by blue and green spheres, respectively.The exchange coupling J 3 (orange) couples NN Cu1 atoms along chains on the b axis, J d (thicker black lines) couples NN Cu2 atoms, and J 1 and J 2 couple neighbouring Cu1 and Cu2 atoms (blue and green lines, respectively).The vectors t 1 = b, t 2 = (a + b)/2 and t 3 = c are primitive translation vectors, and the numbers 1-4 (in red) label the four atoms of the basis.(b) Two layers of Cu 2 OSO 4 , showing the topology of the couplings J ⊥ and J ⊥2 (red and grey lines, respectively).

Figure 1 (
Figure 1 (a) shows the kagome-like structure of each ab layer of Cu 2 OSO 4 , and Fig. 1 (b) shows the way such layers arrange along the c * axis.There are two symmetry nonequivalent Cu 2+ sites in Cu 2 OSO 4 , denoted by Cu1 and Cu2.The nearestneighbour (NN) Cu2 sites form dimers, whereas Cu1 sites form 1D chains along the b axis.Thus, each layer can be thought of as either a system of chains joined through dimers or as a kagome lattice with every third spin replaced by a dimer.The overall spin lattice can be described in terms of the underlying monoclinic Bravais lattice, plus a basis of four atoms, two Cu1 and two Cu2, see Fig.1.Cu 2 OSO 4 crystallizes in a monoclinic crystal system belonging to the space group C2/m[16], where the monoclinic c-axis forms an angle of β = 122.34• with the a-axis.The vector c * is perpendicular to the ab plane, as shown below

FIG. 2 .
FIG. 2. (a) Electronic density of states (DOS) for Cu 2 OSO 4 calculated on the PBE level shows predominant Cu 3d states near the Fermi level, which is denoted by the dashed line at zero energy.(b,c) Orbital-resolved DOS indicates different symmetry of the magnetic orbitals for the Cu1 and Cu2 atoms.(d) Magnetic orbitals visualized using Wannier functions and superexchange pathways for the interlayer couplings J ⊥ and J ⊥2 .Local coordinate frames for the Cu1 and Cu2 atoms are also shown.The notation of oxygen positions follows Ref. [17].

FIG. 3 .
FIG.3.Comparison of total energies for different spin configurations as evaluated by DFT+U and obtained from the minimal model with six exchange couplings listed in TableI.Both sets of energies are given relative to the reference energy E 0 that corresponds to the nonmagnetic state.

FIG. 4 .
FIG. 4. Representative members of the classical ground state manifold of the isotropic model of a single layer of Cu 2 OSO 4 , where Cu2 dimers are represented by one of the three types of vertices of the kagome (A sites).(a) The uniform coplanar ABC state, where A, B and C denote the spin directions of Eq. (12) of the building block of the lattice.(b) One of the 2 L coplanar ground states of a single layer (where L is the number of horizontal rows of the layer).This state results from the one in (a) by swapping B↔C in one of the rows of the layer (shaded).(c) One of the non-coplanar ground states of a single layer.Here, this state results from the one in (a), by rotating the spins belonging to a given row around the direction of A by some arbitrary angle ϕ, which effectively rotates B and C to some new directions B' and C', respectively, out of the plane of A, B and C.

FIG. 5 .
FIG. 5. (a) Local geometry used for the symmetry analysis of the DM vectors, see discussion in the main text.(b) Structure of DM vectors on the layers of Cu 2 OSO 4 , where J 1 and J 2 bonds are projected onto each other.Arrows define the sequence of sites (i, j) in the interaction term D i j • (S i × S j ).The letters A, B and C denote the uniform, 3sublattice state found numerically (see text).

3 = T c * a 3 and T bc * 3 =
T c * b

FIG. 6 .
FIG. 6.(a) Evolution of the spin gap with the rescaling parameter λ of Eq. (47), which characterizes the strength of the symmetric anisotropy T 3 .(b) Same as in (a), but zooming in the intermediate, re-orientation region.(c) Corresponding evolution of the angle ψ, see inset diagrams.

a
and H * * b below H sat .Such a transition is absent for fields along c * .5) According to the results shown in panels (a), (d) and (g) of Fig. 7, the slopes of the longitudinal magnetizations (per Cu) at low fields are approximately 0.015 µ B /T for H ∥ a, 0.008 µ B /T for H ∥ c * , and 0.0015 µ B /T for H ∥ b.The corresponding slopes estimated from the experimental data of Zhao et al at 2 K [Figure 5(b) of Ref. [21]] are 0.015-0.016µ B /T for H ∥ a, 0.012 µ B /T for H ∥ c * , and 0.0028 µ B /T for H ∥ b.The agreement is satisfactory [49].

7 .
FIG. 7. Evolution of the classical ground state of a single layer of Cu 2 OSO 4 as a function of an external magnetic field applied along the a axis (a-c), c * axis (d-f), and b axis (g-i).The results are obtained by numerical minimizations on a single layer of Cu 2 OSO 4 , using the Hamiltonian H iso + H DM + H T 3 + H Z , see text.The first row of panels shows the three components of the total magnetic moment per Cu site, m, in units of µ B .Panels (b) and (e) in the second row shows the combinations of spin components that play the role of order parameters, which set in at the fields H * a and H * c .There is no such phase transition for fields along b, see panel (h).The last row shows the combination of spin components which are driven by the field, in a wider field range.The fields H * * a and H * * b correspond to abrupt re-orientation of the Cu1 spins, before we reach the saturation fields H sat .There is no such transition for fields along c * .

FIG. 8 .
FIG. 8. Linear spin-wave spectra (white dashed lines) and dynamical structure factor intensity S(Q, ω) (symbol size, rescaled to the respective maximum intensity in each separate panel) for a single layer of Cu 2 OSO 4 , along certain cuts of the Brillouin zone, see inset.In (a), the Hamiltonian includes Heisenberg couplings and the c * components of the DM couplings, whereas in (b) it includes Heisenberg, all components of the DM vectors, as well as the interaction matirx T 3 .
i j /N H + O(r 3 ) , (B14) from which we can read off the Curie-Weiss temperature by comparing to the high-T expansion of the Curie-Weiss lawχ µµ = C T − Θ µµ CW = C T + CΘ µµ CW T 2 + O(1/T 3 ) .(B15)We find, for general S , k B Θ µµ CW = − S (S +1) j (J i j + T µµ i j ) ,

TABLE I .
Microscopic magnetic parameters of Cu 2 OSO 4 : leading exchange couplings J i (in K) obtained from FPLO (U d =8.5 eV) and VASP (U d =9.5 eV), and corresponding DM vectors D i j (VASP), calculated for each bond type.Here, the labelling of the two sites ('site i' and 'site j') of each bond follows the notation of Table III in App. A. Note that inversion symmetry forbids the DM interactions on the J d and J ⊥2 bonds.

TABLE II .
Numerical values of various parameters based on DFT values of microscopic couplings.Different rows correspond to different levels of model descriptions.E.g., H iso (FPLO, 2D) corresponds to a single layer of Cu 2 OSO 4 with isotropic Heisenberg couplings only, and values