Ordered ground states of kagome magnets with generic exchange anisotropy

There is a growing family of rare-earth kagome materials with dominant nearest-neighbor interactions and strong spin orbit coupling. The low symmetry of these materials makes theoretical description complicated, with six distinct nearest-neighbor coupling parameters allowed. In this Article, we ask what kinds of classical, ordered, ground states can be expected to occur in these materials, assuming generic (i.e. non-fine-tuned) sets of exchange parameters. We use symmetry analysis to show that there are only five distinct classical ground state phases occurring for generic parameters. The five phases are: (i) a coplanar, 2-fold degenerate, state with vanishing magnetization (${\sf A_1}$), (ii) a noncoplanar, 2-fold degenerate, state with magnetization perpendicular to the kagome plane (${\sf A_2}$), (iii) a coplanar, 6-fold degenerate, state with magnetization lying within the kagome plane (${\sf E}$-coplanar), (iv) a noncoplanar, 6-fold degenerate, state with magnetization lying within a mirror plane of the lattice (${\sf E}$-noncoplanar$_{6}$), (v) a noncoplanar, 12-fold degenerate, state with magnetization in an arbitrary direction (${\sf E}$-noncoplanar$_{12}$). All five are translation invariant (${\bf q}=0$) states. Having found the set of possible ground states, the ground state phase diagram is obtained by comparing numerically optimized energies for each possibility as a function of the coupling parameters. The state ${\sf E}$ noncoplanar$_{12}$ is extremely rare, occupying $<1\%$ of the full phase diagram, so for practical purposes there are four main ordered states likely to occur in anisotropic kagome magnets with dominant nearest neighbor interactions. These results can aid in interpreting recent experiments on ``tripod kagome'' systems R$_3$A$_2$Sb$_3$O$_{14}$, as well as materials closer to the isotropic limit such as Cr- and Fe- jarosites.


I. INTRODUCTION
Frustration can come from various sources. This is certainly true of the frustration exhibited by many magnetic materials, which may be generated by the geometry of the lattice 1,2 , by competition between interactions of different kinds 3,4 or by bond-dependent anisotropies 5,6 . Sometimes, all of these sources of frustration are present at once, making the problem of determining a ground state both more challenging and more rich [7][8][9] .
Kagome lattice rare-earth materials 10-20 provide a realization of this scenario. The kagome lattice [ Fig. 1(a)] is paradigmatic of geometrical frustration while the strong spin-orbit coupling inherent to many rare-earth ions produces complicated anisotropic exchange interactions with distinct, competing, contributions and bond-dependence.
In this Article we study a model of anisotropic exchange on the kagome lattice, including all possible nearest neighbor interactions consistent with the lattice symmetries 8 . This model has six independent coupling parameters, once one allows for the absence of reflection symmetry in the kagome plane, as is appropriate for many materials.
Several previous works have investigated different types of allowed anisotropic nearest-neighbor interaction on the kagome lattice 8,[21][22][23][24][25][26][27][28][29][30] , but none has treated all possible interactions at once, in the absence of reflection symmetry in the plane. Thus, in some sense, these previous works can be viewed as higher-symmetry limits of the generic case studied here. Our goal in this work is to identify the ordered, classical, ground states which are stable over a finite fraction of the six dimensional parameter space of the full model. We will not address the physics at the phase boundaries between different states or limits featuring high symmetry beyond time reversal and lattice symmetries, or cases of accidental degeneracy, although these can be of interest. In this sense, we are study-ing those ground states stable in the presence of "generic" exchange anisotropy.
We find that in the full six-dimensional parameter space there are only five such distinct ground states. They are all translationally invariant, and may be classified by how they transform under the C 3v point group symmetries of the kagome lattice. Example spin configurations for each are shown in Figs. 2-6. In addition to materials with strong exchange anisotropy, our approach is also useful for understanding materials where anisotropy is weak but nevertheless plays a key role in selecting the ground state due to the frustrated nature of Heisenberg interactions on the kagome lattice. Our results can be viewed as illuminating the spectrum of possible ground states which can be obtained by perturbing an isotropic kagome magnet with various allowed forms of nearest-neighbor exchange anisotropy. This may be of use in understanding the ordered ground states of materials including the Cr-and Fejarosites [31][32][33][34][35] and Cd-kapellasite 36 .
The remainder of this Article is organised as follows: • In Section II we review the most general symmetry allowed nearest neighbor exchange Hamiltonian for the kagome lattice 8,35 . We then analyse it in terms of the irreducible representations of the point group C 3v .
• Building on this symmetry analysis, in Section III, we demonstrate the five forms of magnetic order which may arise from the generic Hamiltonian.
• In Section IV we use numerical calculations to calculate the ground state phase diagram of the generic Hamiltonian, delineating the regions of parameter space covered by each of the five ordered phases.
• In Section V we discuss experimental results on kagome materials in the light of our calculations.  1. (a) The kagome lattice, a network of corner sharing triangles. The labels 0, 1, 2 indicate the convention used to label the three sublattices of the kagome lattice in this work. t1 and t2 are the basic translations under which the lattice is symmetric (b) The C3v point group, composed of three reflection symmetries and a threefold rotation axis through the center of the triangle.
• In Section VI we close with a brief summary and discussion of open directions for future work.

II. HAMILTONIAN AND SYMMETRY ANALYSIS
We consider generalized bilinear anisotropic exchange interactions on a kagome lattice [ Fig. 1(a)], We require that the interactions respect the following symmetries: • Time reversal • Lattice translations, t 1 , t 2 , indicated in Fig. 1(a) • Spatial inversion through lattice sites • C 3 rotations around the center of each kagome triangle [ Fig. 1

(b)]
• Reflections in the mirror planes indicated in Fig. 1

(b)
We do not assume any reflection symmetry in the plane of the lattice itself.
We assume that the spins S i transform like magnetic moments, i.e. as axial vectors, odd under time reversal symmetry. This will apply not only when S i is a true magnetic moment but also when it is a pseudospin-1/2 degree of freedom describing the 2-fold degenerate crystal electric field (CEF) ground states of a Kramers ion. In this case the actual magnetic moment is related to the pseudospin via the g-tensor An alternative case is possible in which S i is a pseudospin describing the low energy CEF states of a non-Kramers ion, which will generally be non-degenerate due to a lack of protection from time reversal symmetry. In the non-Kramers case, the pseudospin operators S i will transform differently under time-reversal and the discussion in this section will not apply 37,38 .
We now proceed to constrain the form of the exchange matrices J ij using the symmetries listed above. Time reversal symmetry T is guaranteed by the bilinear form of Eq. (1). There are three spins in the unit cell, which we label S 0 , S 1 , S 2 according to the convention in Fig. 1(b). Translational symmetry imposes that the coupling matrices J ij may only depend on which sublattices i and j belong to and whether the bond ij is on an 'up' or 'down' triangle (red or blue triangles in Fig. 1(a)). Inversion symmetry then guarantees that 'up' and 'down' triangles have the same coupling matrices.
There are thus three different coupling matrices J 01 , J 12 , J 20 entering Eq. (1) which define the interactions between nearest neighbour spins on each pair of sublattices.
The form of the matrices J ij is constrained by the C 3v point group symmetry at the center of each triangle [ Fig. 1(b)], and was given in Refs. 8 and 35: There are six independent parameters in these exchange matrices: three diagonal exchanges J x , J y , J z , two Dzyaloshinskii-Moriya (DM) interactions D y , D z and one symmetric off-diagonal exchange K.
An additional symmetry which could, in principle, be present is reflection symmetry in the plane of the kagome lattice. The presence of such a symmetry would reduce the set of allowed exchange parameters to four, by setting D y = K = 0. This case was discussed in detail in Ref. 8. In this work, we will continue to assume that there is no reflection symmetry in the kagome plane, as is appropriate for many rare-earth kagome materials 17 . Therefore, we shall take both D y and K to be nonzero.
To begin in determining the phase diagram it is helpful to rewrite the Hamiltonian in terms of objects m γ,k transforming according to the irreducible representations (irreps) of the point group. m γ,k are defined for each triangle of the lattice, which we index using k. This approach is discussed for the kagome lattice in Ref. 8  These objects can function as local order parameters for the different kinds of 3-sublattice order which we will encounter on the phase diagram of the anisotropic exchange model. They also aid in the determination of the phase diagram itself. The appropriate objects are defined in Ref. 8 but we reintroduce them here since they are essential to our discussion.
Firstly, there is one scalar object transforming according to the trivial A 1 representation of C 3v . A nonzero average value of this field breaks none of the point-group symmetries, only breaking time-reversal symmetry.
Here S α i,k is the α component of the spin belonging to sublattice i and triangle k.
There are then two linearly independent scalars, which transform according to the A 2 representation. A nonzero average value of these fields breaks time reversal symmetry and all three mirror symmetries of C 3v but preserves the 3-fold rotational symmetry.
Finally, there are three two-component vectors, transforming according to the two dimensional E-irrep of C 3v In terms of these objects the Hamiltonian may be written where k indexes the triangles of the lattice and the final term in Eq. (13) should be interpreted as a sum of 9 scalar products between the vectors m Ei,k . The coefficients λ γ are: It is then useful to write Eq. (13) in a new basis chosen to diagonalize the matrices Λ A2 and Λ E .
Here ω A2i (i = 0, 1) are the eigenvalues of Λ A2 and m A2i are linear combinations of m A2a and m A2b corresponding to the associated eigenvector of Λ A2 [(Eq. 13)]. Similarly, ω Ei (i = 0, 1, 2) are the eigenvalues of Λ E and m Ei are linear combinations of m Ea , m Eb and m Ec corresponding to the associated eigenvector of Λ E . We define, without loss of generality, In this work we will treat the spins as classical vectors of fixed length |S i | = 1. Due to this condition, the following FIG. 2. A1 ordered state, occurring as the ground state of Eq. (1) when λ A 1 < ω A 2 0 , ω E0 [Eq. (24)]. All spins lie in the kagome plane at an angle of 2π 3 to one another and perpendicular to the line joining the spin to the centers of the two neighboring triangles. The spin configuration has vanishing total magnetization, twofold degeneracy, and preserves the lattice symmetries of Eq. (1) while breaking time reversal.
constraint applies to fields m γ defined in Eqs. (7)-(12): Eq. (26) is a necessary but not sufficient condition for the proper normalisation of the spins. It should be emphasised that the reformulation of the problem in terms of variables m γ,k does not require any further assumptions beyond the nearest-neighbor bilinear, nature of the interactions and the symmetries enumerated at the beginning of this section.
In what follows we will seek to find the classical ground states of Eq. (1).

III. WHAT KINDS OF CLASSICAL GROUND STATE ARE POSSIBLE?
In this section we seek to establish the possible classical ordered phases which may occur on the ground state phase diagram of Eq. (1). Our focus is on classical ground states which are stable over finite regions of the full 6-dimensional parameter space. So, although there may be additional ground states which become relevant in particular high symmetry limits of Eq. (1), these are not the subject of our present discussion as they rely on fine-tuning of parameters.
Our conclusions may be summarized as follows: 1. A translation invariant (q = 0) ground state exists for all values of exchange parameters.

If
the ground state is the antiferromagnetic state shown in Fig. 2 and discussed in Section III B.

3.
If There is one spin lying in the kagome plane and two canted out of it in such a way that the three spins remain coplanar, with the plane of coplanarity being tilted with respect to the kagome plane. The plane of coplanarity is indicated by the translucent red planes. There is a net magnetization within the kagome plane. This state breaks time reversal and all of the point group symmetries of the Hamiltonian, apart from a single reflection symmetry which is preserved. It is sixfold degenerate.
the ground state is the noncoplanar state, with magnetization perpendicular to the plane, shown in Fig. 3 and discussed in Section III C.

4.
If the ground state may be one of three states (E-coplanar, E-noncoplanar 6 , E-noncoplanar 12 ) shown in Figs. 4-6 and discussed in Section III D.
A summary of the five phases in terms of the values of local order parameters m γ [Eqs. (7)- (12)] is given in Table I.
In what follows we will demonstrate these results.
TABLE I. Description of the five possible classical ground states in terms of the local order parameters defined for a triangle in Eqs. (7)- (12). E order parameters m Eα are expressed in polar form, with m Eα and ψα as defined in Eqs. (39)- (41). n is an integer, with different choices of n corresponding to different degenerate ground states in E-coplanar and E-noncoplanar6 phases.  A. Existence of q = 0 classical ground state for all parameter sets Here, for completeness, we give the proof that Eq. (1) possesses a q = 0 classical ground state for all values of ex-change parameters, following arguments previously given in Refs. 7 and 8. We follow a strategy of building up the global ground state from the ground states of corner sharing units, as is frequently done for models on lattices with a corner-sharing structure 7,8,29,30,39 .
As we have shown above, the nearest-neighbor exchange Hamiltonian Eq. (1) can be rewritten as a sum over triangles: with H being the same on every triangle of the lattice, as a consequence of inversion and translation symmetries. This formulation makes it clear that any configuration which minimizes the energy of each individual triangle, also minimizes the energy of the system as a whole. Such a configuration may readily be obtained by minimizing the energy on a single "up-pointing" triangle (red triangles in Fig. 1(a)) and then tiling the solution over all "uppointing" triangles of the lattice. The "up-pointing" triangles will then all be in a ground state by construction, and the "down-pointing" triangles will be too, because they have the same exchange matrices as "up-pointing" triangles and the same spin orientation on each sublattice.
This naturally results in a translation invariant (q = 0) state, which is guaranteed to be a ground state. Moreover, it means that the ground state problem on the whole lattice can be reduced to finding the ground state of three spins on a triangle.
In Sections III B-III D we examine the various possible solutions to this problem, that occur in different regions of parameter space.
The argument above does not rule out the existence of additional, q = 0, ground states, degenerate with the q = 0 ones. We regard it, however, as unlikely that such accidental degeneracies are present over finite regions of the 6-dimensional parameter space. Such a robust accidental degeneracy, would require a pair of states not related by any symmetry, to be degenerate with respect to each of the six independent terms of the Hamiltonian individually, which would seem to require a rather large coincidence. The Heisenberg-Kitaev model on the kagome lattice 29,30 exhibits an extended, accidental degeneracy, in the classical limit, but since that model only has two symmetry distinct terms the required coincidence is not so large.
From now on, we assume translationally invariant ground states built by tiling the ground states of a single triangle, and therefore drop the triangle index k from the fields m γ,k and spins S i,k .
We can use the solutions of the single triangle problem to check the validity of the assumption that there are only q = 0 ground states. We do this by checking whether two distinct single triangle ground states can be placed on neighboring triangles without causing an inconsistency at the shared site. If they cannot, then only q = 0 ground states are possible. This is explicitly checked for each single triangle ground state below, and in each case we find that q = 0 states are only possible with fine tuning.

B. A1 order
We first consider the parameter regime defined by inequality (27) where λ A1 is the lowest coefficient in Eq. (24).
We can use Eq. (26) to write: and so eliminate m A1 from the Hamiltonian [Eq. (24)]: All the remaining fields m A2i , m Ei now appear as quadratic forms with positive coefficients, due to inequalities (25) and (27).
Therefore any spin configuration where all these fields vanish is necessarily a ground state, for all parameter sets fulfilling the inequality (27). There are exactly two such configurations, related to each other by time reversal symmetry: These are the ground state spin configurations of the A 1 phase. The only remaining degree of freedom on a triangle is the choice of the + or − sign in Eq. (33). Once this sign is chosen for one triangle, consistency at the shared spin requires that the same sign is chosen on the neighboring triangles. Propagating this throughout the lattice we see that only q = 0 tilings are possible. This phase preserves all lattice symmetries of the original Hamiltonian but breaks time reversal symmetry. One of the ground states is illustrated in Fig. 2.
C. A2 order Next we consider parameter sets falling in the regime described by inequality (28), such that ω A20 is the lowest coefficient in Eq. (24).
Under these conditions we can use Eq. (26) to remove m A20 from the Hamiltonian [Eq. (24)] in a similar manner to the analysis in Section III B. By this means one can show that the ground states for parameter sets obeying the inequality (28) are of the form With the out-of-plane canting angle η being determined by the content of the lowest eigenvector of Λ A2 [Eq. (13)]. In terms of the coupling parameters, η obeys the relation: With η fixed by Eq. (35), the only remaining degree of freedom on a single triangle is the choice of sign in Eq. (34). Once this sign is chosen for one triangle, consistency at the shared spin requires that the same sign is chosen on the neighboring triangles. Propagating this throughout the lattice we see that only q = 0 tilings are possible.
The A 2 configurations have nonzero scalar chirality on the triangle: This phase breaks the reflection and time reversal symmetry of H but preserves the C 3 rotational symmetry. An example ground state in this phase is illustrated in Fig. 3.

D. E orders
We then come to the case Applying the same type of arguments as in Sections III B-III C, we might expect to find a ground state with m A1 = m A2a = m A2b = 0 and with the values of m Ea,b,c being determined by the lowest eigenvector of Λ E . However, for typical eigenvectors of Λ E this is incompatible with the spin length constraints The resolution of this is that the system must mix small values of m A1 , m A2a , m A2b into the ground state, so as to respect the spin length constraints while retaining a large value of |m E0 | as favoured by the Hamiltonian.
We can distinguish the different ways that this can happen by further consideration of the symmetries of the problem. Specifically, we can ask what symmetries of the Hamiltonian can be preserved in the presence of nonzero values of m Ea , m Eb , m Ec .
There are three possibilities consistent with nonzero values of m Eα .
1. One of the reflection symmetries of C 3v is preserved.
This corresponds to the E-coplanar phase discussed below in Section III D 1.
2. The combination of one of the reflection symmetries of C 3v with time reversal is preserved. This corresponds to the E-noncoplanar 6 phase discussed below in Section III D 2.
3. None of the point group symmetries, nor any of their combinations with time reversal symmetry are preserved. This corresponds to the E-noncoplanar 12 phase discussed below in Section III D 3.

E-coplanar
In the E-coplanar phase one of the reflection symmetries of C 3v is preserved. For concreteness, let us suppose that the preserved symmetry is reflection in the yz plane, i.e. the mirror plane that runs through site 2 in Fig. 1 defining the angles ψ Ei to be lie in the interval [0, π), and allowing the scalars m Ei to take either sign ±. Imposing preservation of reflection symmetry in the yz plane constrains ψ Eα More generally, if we had chosen one of the other mirror planes [ Fig. 1(b)] to be preserved, we would have ψ Eα = The symmetry further implies that but a nonzero value of m A1 is allowed and will be mixed into the ground state in such a way as to satisfy the spin length constraints. The magnitudes and relative signs of m Ei , m A1 are fixed by minimizing the energy.
An example spin configuration on the three sublattices in this phase has the form (taking n = 0) where φ and θ are functions of the exchange parameters, which must be determined by minimizing the energy. Degenerate spin configurations can be obtained by applying time reversal and lattice symmetries to Eq. (45) and there is a total degeneracy of six. The spins are in a common plane, which is generally not the plane of the kagome lattice. The total magnetization of the configuration is normal to the unbroken mirror plane. An example configuration is shown in Fig. 4.
Minimizing the energy with respect to θ and φ gives a pair of equations which relate the ground state canting angles to the coupling parameters.
If the angles θ and φ are measured for a given material (e.g. from refinement of Bragg peaks) then Eqs. (46)-(47) can be used to give constraints on the coupling parameters, at least at the level of a classical description. Unless the angles φ, θ are fine tuned to special values (which requires fine tuning of exchange parameters), there is no way to place different members of the set of 6 singletriangle ground states on neighboring triangles. This implies that only q = 0 configurations are possible within this phase, for generic parameters.

E-noncoplanar6
In the E-noncoplanar 6 phase the combination of time reversal with one of the reflection symmetries of C 3v is preserved. For concreteness, let us suppose the preserved symmetry is the combination of time reversal with the mirror plane that runs through site 2 in Fig. 1(b).
This symmetry constrains the angles ψ Eα [Eqs. (39)- (41)], remembering that ψ Eα is defined to lie in the interval [0, π): More generally, if we had chosen one of the other mirror planes [ Fig. 1(b)] to be preserved when in combination with T , we would have ψ Eα = (2n+1)π 6 , n ∈ {0, 1, 2}. If the mirror plane preserved in combination with T runs through site 2 of the unit cell [see Fig. 1(a)] then n = 1, if through site 0 then n = 2, if through site 1 then n = 0.
The symmetry implies that but nonzero values of m A2a and m A2b appear in the ground state as a way to satisfy the spin length constraints An example spin configuration for this phase is S 0 = (cos(ν) sin(µ), sin(ν) sin(µ), cos(µ)) S 1 = (− cos(ν) sin(µ), sin(ν) sin(µ), cos(µ)) S 2 = (0, cos(κ), sin(κ)) The parameters ν, µ and κ are functions of the exchange parameters and must be determined by minimizing the energy. The E-noncoplanar 6 configurations have nonzero scalar chirality on the triangle: The magnetization of the configuration lies within the mirror plane which is unbroken when combined with time reversal. Degenerate spin configurations can be obtained by applying time reversal and lattice symmetries to Eq. (51) and there is a total degeneracy of six.
Different members of the set of 6 ground states cannot be placed on neighboring triangles without causing an inconsistency, unless the angles µ, ν, κ are fine tuned to special values, via fine tuning of exchange parameters. This confirms that only q = 0 configurations are possible within this phase, for generic parameter sets.

E-noncoplanar12
Finally, there is the possibility that time reversal, all point group symmetries and all combinations of the two are broken in the ground state, leaving only translation and inversion symmetries intact.
As shall be shown using numerics in Section IV, this low symmetry state does appear on the ground state phase diagram, but only in a very small region of parameter space.
Minimizing the energy with respect to ζ i , υ i (i = 0, 1, 2) gives a total of six equations relating the canting angles to the coupling parameters.
Thus, if for a system in the E-noncoplanar 12 phase, all six angles are known it should be possible to use Eqs. (59)-(60) to uniquely determine the six exchange parameters.

IV. PHASE DIAGRAM
In this section we calculate the ground state phase diagram of Eq. (1) numerically, by comparing optimized energies for the five phases described in Section III. The numerical optimization of the energy was done by a combination of random search, simulated annealing and iterative minimisation 40  The boundaries of the A 1 and A 2 phases can also be calculated analytically using conditions (27) and (28). These analytic boundaries are shown as white lines in Figs. 7-10, and agree with the results of the numerics. The boundaries between the different E phases are only calculated numerically.
One notable feature of Figs. 7-10 is that the E-coplanar phase is generally found bordering the A 1 phase, whereas the E-noncoplanar 6 phase is generally found bordering the A 2 phase. This is natural since the E-coplanar phase mixes in a finite value of the A 1 order parameter and likewise the Enoncoplanar 6 includes a finite A 2 order parameter.
Another striking feature of the phase diagram is the rarity of the E-noncoplanar 12 phase. This low-symmetry configuration occupies only small portions of the phase diagrams in Figs. 7-10, with its stability generally being increased by a strong negative value of D z .
To investigate the relative frequency of the different phases in the overall parameter space we have calculated the ground state for 100000 different parameter sets, randomly chosen from a uniform distribution on the surface of the 6dimensional hypersphere defined by The pie chart in Fig. 11(a) shows the relative frequency of The phase diagram is obtained by comparing numerically optimized energies for the five phases described in Section III. The numerical optimization procedure is described in Appendix A. The white lines show analytic calculations of the boundaries of the A1 and A2 phases, using conditions (27) and (28). each of the five phases obtained from this procedure. It confirms that E-noncoplanar 12 is indeed a rare phase, found as the ground state for only ∼ 0.5% of randomly generated parameter sets. The four other phases are comparatively common.
This leads us to conclude although the E-noncoplanar 12 state does not require perfect fine tuning to be realized in a kagome material (i.e. it occupies a finite fraction of parameter space), it is unlikely to be realized serendipitously. The other four phases should constitute the classical ground states for the vast majority of kagome materials to which the theory in this paper can be applied (i.e. those with nearest-neighbour, anisotropic interactions).
The above assumes a probability distribution of parameter sets which is isotropic in the 6-dimensional space (J x , J y , J z , D y , D z , K). This may not be the case physically, and indeed it is frequently assumed that the off-diagonal components of the exchange tensor D y , D z , K should be smaller than the diagonal ones J x , J y , J z . We have investigated the distribution of ground states under this assumption, by generating 100000 random parameter sets by choosing J x , J y , J z from a uniform distribution on the surface of the unit sphere: and indepently choosing D y , D z , K from a uniform distribution on the surface of a smaller sphere: The phase diagram is obtained by comparing numerically optimized energies for the five phases described in Section III. The numerical optimization procedure is described in Appendix A. The white lines show analytic calculations of the boundaries of the A1 and A2 phases, using conditions (27) and (28).
The resulting distribution of ground states is shown in Fig.  11(b). The relative frequency of different phases is very similar to that with an isotropic distribution of parameters, although the prevalence of the E-noncoplanar 12 phase increases from ∼ 0.5/% to ∼ 2/%.
A. Phase diagram in the vicinity of the Antiferromagnetic Heisenberg limit The limit J x = J y = J z = J > 0, D y = K = D z = 0, gives the well studied nearest neighbor antiferromagnetic Heisenberg model, which is known to have a highly degenerate ground state 39 . Generic perturbations away from this limit lift the degeneracy, stabilizing a ground state which is unique up to global symmetry operations. Fig. 12 shows the effect of perturbing the Heisenberg model with finite off-diagonal couplings D y , D z , K. D z > 0 strongly favours A 2 order, while D z < 0 favours ordering into the E-coplanar or E-noncoplanar 6 phases depending on which of D y or K is the more dominant perturbation. Our results are in agreement with those of Elhajal et al 21 , who considered the case of perturbing the Heisenberg model with Dzyaloshinskii-Moriya interactions D y , D z , fixing K = 0.
When comparing the results here with those of [21] one should note that the ground state configurations of the E − The phase diagram is obtained by comparing numerically optimized energies for the five phases described in Section III. The numerical optimization procedure is described in Appendix A. The white lines show analytic calculations of the boundaries of the A1 and A2 phases, using conditions (27) and (28). noncoplanar 6 phase become coplanar in the limit of strong positive J and K = 0. This agrees with the labelling of the same phase as coplanar in [21]. Once all symmetry allowed couplings (particularly K) are present, this phase becomes non-coplanar, as identified here.
It is notable that the A 1 phase does not appear at all in Fig.  12. This can be readily understood from the couplings in Eqs. (14)(15)(16)(17)(18)(19)(20). When J x = J y , λ A1 = λ A2,bb . This then implies that ω A20 ≤ λ A1 [cf. Eqs. (25), (27)] with the equality only applying when λ A2,ab = 2( √ 3D y + K) = 0. Thus, when J x = J y the A 2 phase will quite generally have a lower energy than the A 1 phase. A necessary (but not sufficient) condition for the A 1 configurations to be the sole ground states is that λ A1 < λ A2bb =⇒ J x < J y .
The effect of allowing small anisotropy in the transverse exchanges J x , J y is illustrated in Fig. 13. Here we set and vary D y /J and K/J. As implied by the discussion above, δJ ⊥ < 0 favours A 1 order, becoming unstable to the Ecoplanar phase on increasing K. Conversely, when δJ ⊥ > 0 favours A 2 order, which gives way to the E-noncoplanar 6 phase for strong K.  63)]. This models the effect of the assumption that the scale of off-diagonal couplings is lower. The effect on the distribution of phases is minor overall, although assuming weaker off-diagonal exchange expands the size of the rare E-noncoplanar12 from ∼ 0.5% to ∼ 2% In each case, frequencies are determined by numerically finding the ground state for 100000 random parameter sets generated according to the stated distributions.

V. RELEVANCE TO KAGOME MATERIALS
In this section we discuss the application of our results to real kagome materials. We divide our discussion into two areas: firstly, rare-earth magnets belonging to the family 13-20 (sometimes referred to as "tripod kagome" materials 15,17 ), and secondly, Cu, Fe and Cr based magnets where exchange anisotropy should be weaker but nevertheless plays a role in ground state selection. Aside from the systems mentioned below, we anticipate that ongoing work in synthesizing frustrated magnets with strong spin-orbit coupling will reveal new kagome systems to which our results can be applied in the coming years.

A. R3A2Sb3O14 family
In the last few years several rare-earth kagome materials with the general formula R 3 A 2 Sb 3 O 14 have been synthesized. This includes materials with A=Mg, Zn and R=Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb 13,14,17,19 .
Where R is a non-Kramers ion (Pr, Eu, Tb, Ho, Tm), the crystal electric field (CEF) will generally have a non-magnetic singlet ground state, due to the low symmetry of the rare-earth environment. If the gap between this singlet and higher CEF states is smaller than or comparable to the energy scale of interactions, interesting physics may ensue. If the CEF gap is large, the overall ground state of the system will be a trivial singlet driven by the onsite physics. Either way, Eq. (1) cannot describe such physics without being augmented by additional terms, so we will not discuss non-Kramers materials further here.
Where R is a Kramers ion, the CEF will split the 2J+1 multiplet into a series of doublets. At energy and temperature scales below the gap between the lowest and first excited doublet, the magnetism may be represented by pseudospin-1/2 operators S i . S i does not correspond precisely to the magnetic moment, but relates to it via the g-tensor [Eq. 2]. The important thing for our purposes is that S i transforms like a magnetic moment with respect to time-reversal and lattice symmetries, in which case Eqs. (1)-(6) describe the exchange interactions. Below we briefly discuss the various members of the R 3 A 2 Sb 3 O 14 family, with Kramers ions R, in the light of the predictions made in this Article.
The scalar chiral order observed in Nd 3 Mg 2 Sb 3 O 14 16,20 corresponds precisely to the A 2 phase predicted in this work. The magnetic order of the sister compound Nd 3 Zn 2 Sb 3 O 14 has not yet been characterized, but given its essentially similar thermodynamic properties 17 and crystal field environment 18 it seems likely to fall in the same phase as Nd 3 Mg 2 Sb 3 O 14 . Er 3 Mg 2 Sb 3 O 14 was reported in Ref. 17 to avoid long range order down to very low temperatures. It thus appears to be a candidate spin liquid material. The regions near the phase boundaries of the classical phase diagram presented here are likely to be particularly fertile ground for the formation of spin liquid states, and this will be an interesting direction for future research. Er 3 Zn 2 Sb 3 O 14 exhibits strong structural disorder and associated glassy behavior of the magnetic properties 17 , which is beyond the scope of our present discussion.
Yb 3 Mg 2 Sb 3 O 14 exhibits long range order at T N ≈ 0.88K 17 . The form of this magnetic order has yet to be reported in the literature. Based on the expectation that, as a rare earth magnet with moderate magnetic moment, the theory in this manuscript should be applicable to Yb 3 Mg 2 Sb 3 O 14 , we expect that the order will be one of the states discussed in this work. Like Er 3 Zn 2 Sb 3 O 14 , Yb 3 Zn 2 Sb 3 O 14 has strong structural disorder, although unlike the Er compound it does not show clear signs of spin freezing 17 .
Sm 3 Mg 2 Sb 3 O 14 13 and Sm 3 Zn 2 Sb 3 O 14 14 have both been synthesized but their low temperature magnetism has yet to be characterized in detail. This may be challenging due to the small magnetic moment of the Sm 3+ ion, but recent experiments on the pyrochlores Sm 2 Ti 2 O 7 and Sm 2 Sn 2 O 7 indicate that this is possible 42 . There is some evidence of hysteresis in the low temperature magnetization curve for Sm 3 Zn 2 Sb 3 O 14 14 but not for Sm 3 Mg 2 Sb 3 O 14 13 , which may provide some clue as to the low temperature state.
Materials with R=Gd present a somewhat different case, because Hund's rules imply vanishing orbital angular momentum L = 0 for the Gd 3+ ion. The magnetism on the Gd sites thus comes from a pure S = 7/2 spin and anisotropies in the interactions should be much weaker. Some understanding of this case can be gained from considering a model with nearest neighbor Heisenberg exchange and the nearest-neighbor part of the dipolar interaction: In terms of the symmetry-allowed interaction matrices [Eqs. with the full long ranged dipolar interaction D, and found the A 1 configuration as the ground state for weak to moderate D and antiferromagnetic J. It also agrees with previous predictions about the ground state of Gd 3 Mg 2 Sb 3 O 14 15 , and with the observed antiferromagnetic transition at T N ≈ 1.7K 15,44 , although differences between the field cooled and zero-field cooled susceptibility 44 remain to be understood.
For R=Dy the ionic magnetic moment is very large and the long range component of the dipolar interaction cannot be ignored. Dy 3 Mg 2 Sb 3 O 14 exhibits an unusual "fragmented" 45 phase where there is an ordering of emergent "charge" degrees of freedom while spins remain partially disordered 46 . The long-range dipole-dipole interaction plays a crucial role in this phenomenon 47,48 and thus it is beyond the scope of the theory presented in this Article.

B. Nearly isotropic systems
While the most obvious application of the results in this Article is found in systems where exchange anisotropy is strong, our results can also be applied to understand cases where isotropic Heisenberg exchange is weakly perturbed by short ranged anisotropic interactions. This is the case in the Fe-and Cr-jarosites AM 3 (OH) 6 (SO 4 ) 2 where M= {Fe, Cr} and A={K, Rb, NH 4 , Na} [31][32][33][34][35] . These are found to order in the A 2 phasethe most prevalent of our phase diagram. This is generally understood to be a consequence of antiferromagnetic Heisenberg exchange perturbed by a weak D y . This interpretation fully agrees with the results presented here: it can readily be checked that inserting into Eqs. (14)- (23) gives an outcome obeying condition (28) and hence a ground state in the A 2 phase [cf. Fig. 12]. What this work adds to the discussion is a simple and systematic approach to finding the preferred ground state for general kinds of anisotropic nearest neighbor perturbation. An example where weak anisotropic perturbations away from a Heisenberg model lead to something other than A 2 order is given by Cd-kapellasite 36 . The weak ferromagnetic moment confined within the kagome planes in that material is only consistent with the E-coplanar phase, out of the phases in this Article.

VI. SUMMARY AND DISCUSSION
In this Article we have developed a theory of the magnetic orders induced by nearest-neighbor exchange anisotropy in kagome magnets. Our theory reveals that five distinct magnetic orders can be expected from such interactions, all retaining the translational symmetry of the lattice, but being distinguished from one another by their transformations under timereversal and point group symmetries. The five phases are: A 1 [ Fig. 2], A 2 [Fig. 3], E-coplanar [ Fig. 4], E-noncoplanar 6 [ Fig. 5], E-noncoplanar 12 [ Fig. 6]. They are labelled according to the irreducible representation of the point group C 3v with which the primary order parameter transforms, their coplanar or noncoplanar nature and their degeneracy. Eqs. (27)- (28) give exact conditions for the A 1 and A 2 configurations to be classical ground states.
We have used numerical calculations to determine the full zero temperature phase diagram of the most general anisotropic nearest-neighbor exchange model, showing the extent of these five phases . One of the five phases (E-noncoplanar 12 ) is found to be exceedingly rare in the parameter space [ Fig. 11].
We have discussed how this theory relates to various real kagome materials [Section V], with both strong and weak exchange anisotropy.
The dominance of noncollinear (A 1 , E-coplanar) and noncoplanar (A 2 , E-noncoplanar 6,12 ) states on the phase diagram suggests a high possibility of spin excitations with topological band structures in many kagome materials [49][50][51] . It is likely that the five phases identified here from analysis of broken symmetries can be subdivided further by the topology of the excitation bands. Relatedly, the possibility of coupling to itinerant electrons is an interesting area for future research with a view to investigating topological transport phenomena.
The approach used in this work relies on the ability to decompose the Hamiltonian into a sum over blocks, such that the ground state is obtained by finding the ground state on each block and tiling it over the lattice. This would seem to limit the usefulness of the approach for systems with further neighbor interactions, since such a decomposition may either not be possible or may require such large blocks that the decomposition is no longer a useful simplification. Applying the method from this work to quantum systems will also not be possible in general -even for nearest neighbor interactions -because the Hamiltonians on neighboring blocks will usually not commute. There are, however, some specific, fine-tuned, cases where the exact ground state of a quantum system can be built up by such a block-by-block approach 28,52 .
While we have restricted ourselves here to phases which are stable over finite regions of the classical phase diagram, a study of the phase boundaries may also be interesting. As has been studied elsewhere 7,8 phase boundaries between competing classical phases can host non-trivial enlarged manifolds of zero-energy states, which in some cases are associated with new forms of spin liquid 53 . In general, the greater the degree of degeneracy around the phase boundary, the more more favorable the situation becomes towards the formation of spin liquids. Different phase boundaries will have different amounts of additional degeneracy and so some will be more favorable for spin liquid formation than others. Boundaries where 3 (rather than just 2) phases meet may host particularly interesting physics as seen in (e.g.) [53]. An analysis of each possible phase boundary would be an interesting undertaking, which we leave open for future work.
Here we describe the numerical optimization used to obtain the phase diagrams in Figs. 7-10 and the estimates of the relative frequency of phases in Fig. 11.
For a given parameter set, the energy is optimized separately for each of the five phases described in Section III and then the optimized energies are compared to determine which is the lowest.
Due to the argument in Section III A, we need only optimize the configuration on a single triangle, since we know that a ground state on the full lattice can be obtained by tiling the ground state of a single triangle everywhere.
The optimization for each phase is done by either random search or simulated annealing combined with iterative minimization 40 , apart from the A 1 phase where the spin configuration is fixed [Eq. 33] and thus the corresponding energy can directly be calculated without any optimization being necessary: For the other four phases (A 2 , E-coplanar, E-noncoplanar 6 , E-noncoplanar 12 ), the optimization procedure is as described below.

Optimizing A2 configuration
The form for the A 2 configurations is given in Eq. (34). This can be written as with (s a , s b ) on the unit circle Initially, we calculate the energy for 10 5 randomly generated values of (s a , s b ) on the unit circle. The lowest energy configuration obtained from this random search is then used as input for the iterative minimization step.
In the iterative minimization step (s a , s b ) are updated as For sufficiently small, positive, c this update is guaranteed to reduce the energy, unless the system is already in a locally optimal configuration before the update. The parameter c is initially set to 0.1. If the update (A6) does not reduce the energy then c is reduced by a factor of 2 and the update is attempted again. This procedure is repeated until the configuration converges.

Optimizing E-coplanar configuration
The form for an E-coplanar configuration is given in Eq. (45). This can be rewritten as with (σ x , σ y , σ z ) on the unit sphere Initially, we calculate the energy for 10 5 randomly generated values of (σ x , σ y , σ z ) on the unit sphere. The lowest energy configuration obtained from this random search is then used as input for the iterative minimization step.
In the iterative minimization step (σ x , σ y , σ z ) are updated as The parameter c is initially set to 0.1. If the update (A11) does not reduce the energy then c is reduced by a factor of 2 and the update is attempted again. This procedure is repeated until the configuration converges.
The set of configurations covered by the E-coplanar ansatz (45) includes the A 1 configurations (when φ = 4π 3 , θ = π 2 ). Because of this, if the E-coplanar optimization is found to give the lowest energy of the five possibilities we must check that the obtained configuration has a nonzero value of at least one of the order parameters m Eα . In practice we check that If the E-coplanar optimization obtains the lowest energy but the inequality (A12) is not fulfilled, the ground state is assigned to the A 1 phase.
or where J 2 x + J 2 y + J 2 z = 1 (for Fig. 11(b)). The triangle is swept 10 5 times at a given temperature, and the temperature is then reduced by a factor of 0.9. This procedure is repeated 200 times. There are than 10 5 sweeps of the triangle with T = 0, i.e. only accepting energy reducing updates.
The whole annealing procedure is performed from the start 3 times for each parameter set with the final output being the lowest energy configuration obtained over all three sweeps.
To converge the configuration further, there is then an iterative minimisation step where each spin component is updated as: The parameter c is initially set to 0.1. If the update (A19) does not reduce the energy then c is reduced by a factor of 2 and the update is attempted again. This procedure is repeated until the configuration converges. If the energy produced from this procedure is lower than the energy produced from optimizing within the A 1 , A 2 , Ecoplanar or E-noncoplanar 6 phases, then the ground state may be within the E-noncoplanar 12 phase. Because the configuration on the triangle is completely general, to confirm that the configuration has not converged to one of the other phases we check that the inequality (A12) is satisfied, and also check that: If inequalities (A12), (A20), (A21) are not satisfied, the ground state is assigned to one of the other phases depending on the values of the various m γ [ Table I].
Supplemental Material: Ordered ground states of kagome magnets with generic exchange anisotropy In this Supplemental Material we present phase diagrams for a series of additional values of D y /|J z |, D z /|J z |, K/|J z | from {−0.75, −0.25, 0.25, 0.75}, with both signs of J z .