Spiral order from orientationally correlated random bonds in classical XY models

We discuss the stability of ferromagnetic long-range order in three-dimensional classical XY ferromagnets upon substitution of a small subset of equally oriented bonds by impurity bonds, on which the ferromagnetic exchange J_perp>0 is replaced by a strong antiferromagnetic coupling J_imp<0. In the presence of a single impurity bond, once the absolute value of the frustrating coupling J_imp<0 exceeds a threshold J_c>0, the ground state becomes two-fold degenerate, corresponding to either clockwise or anticlockwise canting of the spins in the vicinity of the impurity bond. In the presence of a small concentration of impurity bonds, the effective low-energy Hamiltonian is that of Ising variables encoding the sense of rotation of the local canting around the impurities. Those degrees of freedom interact through a dipolar interaction mediated by spin waves. A ferromagnetic Ising ground state indicates the instability of the XY ferromagnet towards a spiral state with a wave vector proportional to the concentration of impurity bonds. To analyze under which circumstances such a ground state arises, we study first impurities forming superlattices. For a subclass of those, we can rigorously establish the existence of spiral order. For another class of superlattices, the Ising variables order ferromagnetically in planes perpendicular to the orientation of impurity bonds, but antiferromagnetically parallel to it, which results in a fan-like XY ground state. Second, we consider the case when the impurity bonds are randomly distributed on the three-dimensional host lattice according to a Poisson process. We show the phenomenon of spiral order by disorder with an ordering wave vector proportional to the impurity concentration. The analytical predictions are confirmed by Monte Carlo simulations and are relevant for magnetic materials such as YBaCuFeO_5.


I. INTRODUCTION
Insulating magnets supporting long-range magnetic spiral order are of technological interest as they can display "magnetically" induced ferroelectricity 1-4 . In prototypical spin-spiral multiferroics, e.g., RMnO 3 (R=Tb 3+ , Dy 3+ , etc.) 5,6 , a magnetic spiral phase can be stabilized by the competition between nearest-neighbor and further-neighbor magnetic exchange interactions with opposite signs 7,8 . However, the resulting frustration only induces spiral states if further-neighbor couplings are sufficiently strong as compared to nearest-neighbor couplings. The latter are typically much bigger in magnitude, except under special circumstances that lead to their suppression. In such exceptional cases, the characteristic exchange scale is set by the further-neighbor interactions and is thus very weak, entailing a low spiral ordering temperature.
In order to engineer magnetic insulators with magnetic spiral order establishing at high temperatures, it is of fundamental interest to investigate analogous mechanisms. An interesting route was suggested by the study of Ivanov et al. 9 who considered a Heisenberg antiferromagnet on a square lattice, in which every other horizontal nearestneighbor bond in a staggered pattern was replaced by a ferromagnetic coupling. Sufficiently strongly frustrating bonds were shown to induce a magnetic spiral order. From the experimental side, there are interesting hints that a similar mechanism might be tied to the presence of disorder. Indeed, certain insulating compounds containing some degree of chemical disorder were reported to stabilize magnetic spiral order [10][11][12][13][14] at high temperatures. For example, the transition temperatures to the magnetic spiral phase were found to range from 180 K to 310 K 12,14-17 in YBaCuFeO 5 , whereby several further characteristics of the spiral depend on the degree of disorder. This empiric observation suggests the possibility that, for some materials, a magnetic spiral order might be induced by some "impurity bonds" formed by nearest-neighbor magnetic ions whose exchange coupling frustrates the order that would establish in their absence. Recent Monte Carlo simulations have confirmed this conjecture in a model describing YBaCuFeO 5 with disorder in the spatial location of the magnetic Cu and Fe ions 18 . The latter was assumed to result in a small concentration of locally frustrating bonds along the c-direction, which indeed was shown to induce magnetic spiral order in an experimentally relevant window of parameters.
In this paper, we describe and study the general mechanism that renders the ferromagnetic long-range order of classical XY spins unstable towards spiral order, when a finite fraction of the ferromagnetic interactions is replaced by sufficiently strong antiferromagnetic exchange couplings. In the end, we will confront the theory with experimental data, as shown in Fig. 1, with good quantitative agreement.
The general physical mechanism at work is the fol- 13 for the compound YBaCuFeO5, obtained by different annealing rates that control the fraction of bipyramids hosting magnetically frustrating Fe-Fe pairs instead of low energy Fe-Cu pairs. The transition temperature to the spiral state is plotted versus the spiral wavevector measured at low temperature. Our theory explains the linear relationship between the two and predicts the slope as a function of the three exchange couplings. We plot the predicted slope for different values of the exchange on the driving frustrating bonds, fixing the other two exchanges to the values estimated in Ref. 18. The larger the frustrating coupling, and the higher the fraction of frustrating Fe-Fe bonds in the bipyramids of YBaCuFeO5, the larger the transition temperature and the spiral wavevector. Each triplet of triangular symbols joined by a dashed line corresponds to the same concentration of frustrating bonds in our theory.
lowing. We consider a geometrically unfrustrated lattice Λ in d > 2 dimensions, hosting isotropic spins with a continuous symmetry. The symmetry is broken spontaneously at low temperatures, which implies the existence of Goldstone modes. Dilute but strong impurity bonds embedded in this lattice can induce local cantings which behave as "dipole type" defects with an Ising degree of freedom associated to them. The Goldstone modes mediate an interaction between the defects, decaying as r −d for large separation. Correlations in the distribution of such impurity bonds (e.g., a restriction to bonds that point in a single direction) may ensure a sufficiently nonfrustrated pairwise interaction between these defects so as to favor long-range ferromagnetic order in the orientation of the local cantings. Such long-range Ising order entails a global twist of the ferromagnetic order parameter density, and thus a magnetic spiral, as the local magnetization twists in the same sense across every impurity bond. The wave vector of the resulting magnetic spiral is proportional to the magnetization density of the Ising degrees of freedom, and thus, to the density of impu-rity bonds. We will show that such a spiral state is the ground state of the XY system for a rather wide range of parameters of the impurity bond distribution. For simplicity, we consider a cubic host lattice Λ embedded in three-dimensional Euclidean space with the Cartesian coordinates x, y, and z. We impose a tetragonal symmetry by choosing the ferromagnetic nearestneighbor exchange to be J > 0 for couplings in the x-y plane and J ⊥ > 0 for bonds oriented along the z-axis. We further consider a set of impurity bonds, which form a dilute subset of the nearest-neighbor bonds that are directed along the c-direction of the cubic host lattice Λ. For each impurity bond, the ferromagnetic J ⊥ > 0 is replaced by the antiferromagnetic exchange coupling J imp < 0.
A single impurity bond does not destroy the ferromagnetic long-range order of the ground state. However, it does result in a canting of the classical XY spins in the vicinity of the impurity bond, provided that the local frustration is sufficiently strong, i.e., |J imp | ≥ J c for some threshold value J c > 0. Under these conditions (and with fixed boundary conditions at infinity) the ground state is two-fold degenerate, exhibiting either a clockwise or counter-clockwise sense of the local canting. At low concentration we can thus associate a corresponding low-energy Ising degree of freedom to every impurity bond. Apart from these discrete soft degrees of freedom, the background ferromagnet hosts low-energy spin wave excitations. They mediate an effective interaction between the Ising degrees of freedom, which results in an effective classical Ising model with effective two-body interactions of dipolar type. Their algebraic decay at large distance is a direct consequence of the gaplessness of the spin waves. A similar effective interaction results in any system that spontaneously breaks a continuous symmetry, and thus hosts gapless Goldstone modes mediating algebraic interactions between impurity degrees of freedom. The specific case where the impurity bonds form a Bravais superlattice, with a unit cell that is large compared to that of the cubic host lattice Λ, is analytically tractable, and for certain classes of superlattices we are able to rigorously establish the presence of spiral order.
Although the analysis in this paper assumes ferromagnetic interactions for the host lattice Λ, we note that our results can be readily extended to any unfrustrated XY magnet. For example, if the lattice is bipartite and J imp has sign opposite to J ⊥ , the system can be mapped to the above described ferromagnet as follows. For every spin, a reference frame is chosen such that the unfrustrated ground state of the impurity-free system corresponds to a ferromagnetic configuration. By virtue of this mapping, the low energy effective theory presented below extends to this larger class of magnetic insulators.
We emphasize that for the establishment of ferromagnetic order it is central that the impurity bonds be not randomly oriented. Otherwise, the pair-wise interactions between the associated Ising degrees of freedom would be strongly random in sign, which would most likely lead to spin glass order, as observed in models of dilute, randomly oriented Ising dipoles 19,20 . Since an Ising glass state generally carries no net magnetization, it would not induce a spiral state of the original XY spins. Also, in the limit of a high density of randomly oriented impurity bonds, one expects long-range spin-glass order (observed directly at the level of the XY spins), since the model becomes that of a random-bond XY gauge glass, as was studied by Villain [21][22][23] .
The remainder of this paper is structured as follows. In Section II, we define the spin lattice model. Section III begins with the case of a single impurity bond. We then consider a small concentration of impurity bonds and derive a mapping to an effective Ising model for low energies. Section IV describes how to find the ground state of the effective Ising model when the impurity bonds realize a superlattice. The effective Ising model is solved by analytical and numerical means. Its solution is then compared to Monte Carlo simulations of a model with the same network of exchange interactions, but in which the classical XY spins are replaced by classical Heisenberg spins with an additional easy-plane anisotropy. The latter allows for close contact with experimentally realized magnets, such as YBaCuFeO 5 , which are believed to embody the physical ingredients and mechanisms discussed above. In Section V, we study the analytically tractable case of dilute, randomly placed impurity bonds, and conclude that random samples undergo a ferromagnetic ordering of cantings, and thus form a spiral phase. Section VI addresses the onset temperature of the spiral phase and predicts it to be proportional to the spiral wavevector. We compare our quantitative predictions with experimental observations. Section VII summarizes our findings and discusses how the general mechanism identified here applies to other systems.

II. LATTICE HAMILTONIAN FOR CLASSICAL XY SPINS
A. Definition of the XY −model We consider a magnet of classical spins, described by two-dimensional unit vectors S r with S 2 r = 1. They are located at the sites r = x x + y y + z z (x, y, z ∈ Z) of a cubic lattice Λ made of |Λ| sites spanned by the orthonormal unit vectors x, y, and z of R 3 . In most cases we will impose periodic boundary conditions on these classical spins. However, as usual, the choice of boundary conditions does not affect the bulk properties.
We consider a classical Hamiltonian containing only nearest-neigbor interactions between spins. Here, L denotes the set of impurity bonds, as we will describe below.
The exchange Hamiltonian in Eq. (2.1a) possesses the translation symmetries of the cubic lattice, since its nearest-neighbor ferromagnetic Heisenberg exchange couplings depend only on the relative position of the spins, (2.1c) The in-plane (J ) and out-of-plane (J ⊥ ) couplings are ferromagnetic but can be different, 0 < J = J ⊥ , in which case the cubic point-group symmetry is reduced to the tetragonal one.
The contribution from the disorder in Eq. (2.1a) describes the presence of antiferromagnetic impurity bonds. We label the bonds by the coordinate of the end point with the smaller z-coordinate. These end points form a subset L of the points of the cubic host lattice Λ. This term breaks the lattice translation symmetry. On all impurity bonds the ferromagnetic J ⊥ > 0 is replaced by the antiferromagnetic coupling J imp < 0, inducing local frustration. Hamiltonian (2.1a) is invariant under any rotation of all spins by the same orthogonal 2 × 2 matrix, i.e., H L has a global O(2) symmetry.

B. Impurity-free case
Here we consider an impurity-free system, i.e., an empty set L, The ground state is ferromagnetic with all spins parallel. We choose the polar parametrization with the orthonormal basis x and y of R 2 . In this polar representation, At low temperatures, T J ⊥ , J , we can use the spin-wave approximation, which assumes that the deviations from the ferromagnetic ground state (2.5) are small. In that case, the Hamiltonian (2.4) can be expanded to quadratic order in the angle differences, is the energy of the ferromagnetic ground state, and is the symmetric spin-wave kernel. It only depends on the difference r − r, which we henceforth use as the only subscript. We close Sec. II by establishing a few important properties obeyed by the spin-wave kernel D (0) . We observe that D This is a consequence of spin rotational symmetry, which implies that any global orthogonal transformation (2.6) leaves the bilinear form (2.7a) invariant. Moreover, if we impose that the angles θ r obey periodic boundary conditions, we then have the Fourier transform for any k belonging to the Brillouin zone of the host cubic lattice Λ. We shall denote this Brillouin zone by BZ(Λ). Finally, it is convenient to introduce the inverse of the spin-wave kernel D (0) as the Green's function G (0) , which satisfies r −r = δ r,r , ∀r, r .
(2.10a) Due to the zero mode (2.8), G r−r is defined up to a constant, which we fix by requiring that r∈Λ G (0) r = 0, (2.10b) such that both G (0) and D (0) annihilate constant functions. As the inverse of a symmetric kernel, G r−r is symmetric, (2.10c) Imposing periodic boundary conditions, we have (2.10d) The asymptotic large distance behavior of the Green's function is (2.11) On the right-hand side of Eq. (2.11), we recognize the three-dimensional Coulomb potential for the rescaled coordinates We will see in the next section that impurities couple to each other through the combination of Green's functions with the Fourier transform of Γ .
k=0 does not enter the sum. It is therefore convenient to defineΓ (2.14) where we use the notation of Eq. (2.12).

A. Periodic boundary conditions
We are ultimately interested in describing states with a spiralling magnetic order, where the angles of the local magnetization grow linearly with distance. However, before dealing with this possibility, we first analyze a situation anticipating no spiralling. In this case we can impose periodic boundary conditions on the spin angles.
We recall that the set of all directed impurity bonds defines the set L ⊂ Λ consisting of all the sitesr ∈ Λ such that r,r +z is a directed impurity bond. We anticipate that the twist angles across bonds are relatively small, except on the impurity bonds. We therefore make the approximation where ∆θr is the canting angle across the impurity bond, This approximation is certainly valid on the majority of bonds in the limit of a dilute concentration of impurity bonds, which we will assume from now on. We are going to establish under what conditions there is a configuration of angles other than the ferromagnetic one that minimize the energy (3.1a) in the presence of impurity bonds anchored on the set L.
First, we fix all the angles θr and θr +z withr ∈ L and integrate out the angles on all other sites, i.e., on all sites from Λ \ (L ∪ L + z). Within our spin wave approximation, and treating the angles as non-compact variables, this can be done exactly (at any temperature where the spin wave approximation is justified) since those angular variables enter the Hamiltonian quadratically. Here, we carry out the calculation at T = 0 by solving the saddlepoint equations for all angles on the sites Λ \ (L ∪ L + z). In this way, we shall find the angles θ r , as well as an effective Hamiltonian, expressed solely in terms of the angles {θr,r ∈ L ∪ L + z}.
The minimization over all r ∈ Λ \ (L ∪ L + z) requires wherebyr ∈ L ∪ L + z runs over all sites in Λ that are end points of an impurity bond (two sites per impurity bond). Here, χr has to be chosen in such a way that the saddle-point equations (3.3) are satisfied as θr and θr +z take their prescribed values. Inverting Eqs. (3.3a) and (3.3b) for any r ∈ Λ, we obtain where θ 0 ∈ [0, 2π[ is the angle of the magnetization far away from impurities. Restricting ourselves to the subset r ∈ L ∪ L + z, Eq. (3.4) can be inverted to yield where G (0) is the 2|L| × 2|L| matrix obtained by restricting G The algebraic decay of G (0) r ensures that the distortions in the angular pattern also have algebraic tails.
Second, we evaluate the energy (3.1a) at the saddle point (3.6). We thereby obtain the effective energy H eff θr, θr +z ,r ∈ L . . = 1 2 r,r ∈L∪L+z θr G where we have assumed that the magnetization far away from the impurities are oriented along θ 0 = 0, making use of the fact that the constant θ 0 can be chosen freely, since rotating all spins by the same angle does not affect the energy. In what follows, we will drop such additive constants. Note that the derivation of the effective action (3.7), starting from the Hamiltonian (3.1a), is exact at any temperature, if we use the Gaussian approximation and treat the domain of the angles θ r as non-compact, ignoring that the energy is in fact 2π-periodic in the angles.
For the third and last step we assume T = 0. We minimize the effective action (3.7) with respect to the impurity bond angles θr and θr +z withr ∈ L. This can be done exactly for a single impurity bond, and approximately in the case of a dilute set of impurities.

The case of a single impurity bond
In the case of a single impurity bond r,r + z represented by the single site L = {r}, the inverse of the 2 × 2 symmetric matrix with elements G where r, r ∈ {r,r + z}.
If we use Eq.
Minimization over θr +θr +z imposes the condition that the two angular distortions away from the asymptotic θ 0 on either side of the impurity bond are opposite, The remaining degree of freedom is the canting angle across the impurity bond, ∆θr ≡ θr − θr +z , (3. 12) with −π < ∆θr ≤ π, in terms of which the effective restricted (eff/res) Hamiltonian becomes (3.13a) where we have introduced the short-hand notation The coupling J c depends parametrically on J and J ⊥ . The rationale for the subscript in J c is the following. For eff/res (∆θr) has a single minimum at ∆θr = 0. (3.14) eff/res (∆θr) develops a double well with two degenerate minima. The two degenerate minima occur at the relative canting angles ∆θr = σr ∆θ, (3.15a) with σr = ±1 and ∆θ being the positive solution of When |J imp | > J c and the temperature T is sufficiently small, namely, 16) thermal fluctuations around the minima of the double well are small. We will thus call |∆θr| a hard degree of freedom, while we refer to the Ising variable sgn(∆θr) = σr as a soft degree of freedom. This terminology is motivated by expanding about the minimum of the double potential well, σr ∆θ, which is closest to ∆θr When the temperature is small compared to the curvature at the two minima, thermal fluctuations of the hard degree of freedom |∆θr| are much smaller than fluctuations due to the soft degree of freedom σr. We will thus ignore the latter in the range (3.16) of temperatures. If we use the values of the relative angle (3.1b) at the pair of minima (3.15) in combination with Eqs. (3.11) and (3.6) with θ 0 = 0, we obtain the two canting patterns shown in Fig. 2.
r=0 entering J c in Eq. (3.13b) can be expressed with the help of Eq. (2.13a) as . (3.18) For isotropic couplings J = J ⊥ = J, this evaluates to which yields the critical coupling For general couplings with ratio we plot the threshold value of J c /J ⊥ in Fig. 3 together with asymptotic expressions that become valid in the limit of strong anisotropy.

B. The case of a dilute set of impurity bonds
When the density of impurity bonds is finite, we cannot rely on the explicit representation (3.8). However, at a small impurity concentration, the interaction between the angles on the same impurity bond is much stronger than the coupling between angles on different impurity bonds. It thus makes sense to split G (0) into a bond-local term [with inverse on every bond given by Eq.
a bond-off-diagonal term according to (3. 22) and to approximate its inverse as This yields the Hamiltonian where H (1) eff (θr, θr +z ) is given by Eq. (3.9). Since this term is dominant at low impurity concentration, it is again reasonable to restrict the angular configurations to the subspace given by θr = −θr +z = ∆θr 2 =⇒ θr + θr +z = 0, θr − θr +z = ∆θr, (3. 25) as in the single impurity-bond problem of Eq. (3.9). This leads to the effective restricted Hamiltonian [compare with Eq. (3.13)] where the interaction Γ (0) r−r is seen to be mediated by the combination of Green's functions introduced in Eq. (2.13a) that scales like an anti-dipolar interaction at long distances.
It remains to minimize H eff/res ({∆θr}) with respect to the canting angles ∆θr on the impurity bonds. The corresponding saddle point equation reads, As before, the ferromagnetic state (3.14) is a solution of Eqs. (3.27a) and (3.27b), but it becomes unstable for sufficiently large |J imp |.
The term Ξr is expected to be dominated by the closest neighbor bonds, since the sum is over a set of decreasing terms, which, for a valid saddle point, contribute with alternating signs. Using Eqs. (2.14) and (3.27b), one expects that Ξr scales as n imp . In that case, the condition will be met for sufficiently large values of J imp and sufficiently small values of n imp . Under the condition (3.28), local minima of the Hamiltonian (3.1a), i.e., solutions to Eq. (3.27), exist, which locally look like the solution for a single antiferromagnetic impurity bond.
In the dilute impurity regime, a low energy state will have a canting of angles across the impurities bonds close to the single impurity case, i.e., ∆θr ≈ σr ∆θ. (3.29a) Inserting this Ansatz into the effective Hamiltonian (3.26), we obtain the effective Ising model r−r σr . r between a pair of Ising variables a distancer apart that enters on the right-hand side of Eq. (3.29b) was derived in Eq. (2.14). Its asymptotic behavior is that of Ising dipoles (σr oriented along the direction z), albeit with the opposite sign as compared to the usual dipolar interaction. In the dilute impurity limit, the sign of the two-body interaction Γ (0) r depends on the relative positionr (see Fig. 4). It is ferromagnetic whenx (2.14) between cantings, in theρ 2 −z 2 -plane (with ρ 2 ≡x 2 +ỹ 2 ). Along the z-axis the interactions are antiferromagnetic, while for separation along the xy−plane they are ferromagnetic. The sign change occurs on the cone described by J ⊥ρ 2 = 2J z 2 . In the quasi-one-dimensional limit vanishing on the conical surfacẽ and antiferromagnetic wheñ (3.30c) C. Boundary conditions along the z axis So far, we have considered the ferromagnetic state, θ r = 0 for all r ∈ Λ, of H 0 defined in Eq. (2.4) and performed a spin-wave expansion about it, imposing periodic boundary conditions. However, this precludes the possibility that the reference state for the spin wave expansion is non-collinear, and moreover, the periodic boundary conditions on the angles θ r preclude a spiral state. To overcome these limitations, we allow for a linear growth of along the z direction. Here, the global degree of freedom Q ∈ [−π, π[ describes a constant twist rate, while the local degrees of freedom φ r obey periodic boundary conditions. In finite systems, the twist rate Q along the z direction should be an integer multiple of 2π/L z if we impose periodic boundary conditions on the original spins S r = cos θ r x + sin θ r y, but this discrete constraint is irrelevant in the thermodynamic limit. With the change of variables (3.31), the spin-wave approximation (3.1a) becomes where we recall that |Λ| is the number of sites in the host cubic lattice Λ, and we again denote by ∆φr .
the twist of φ across the impurity bond labelled byr. For a low impurity concentration 0 < n imp 1, we assume and will verify a posteriori, that where ∆θ≥ 0 is the modulus of the canting angle across an isolated impurity bond. The leading effect of the emerging spiral order will appear at order O(n 2 imp ) in the energy per spin. We therefore expand the impurity bond terms in Eq. (3.32a), i.e., the second line on the right-hand side of Eq. (3.32a), up to linear order in Q. In this approximation, the saddle-point values for ∆φr are therefore again given by ∆φr ≈ σr ∆θ (3. 34) up to corrections of order O(Q). We can neglect those since they lead to corrections to the energy per spin of order O(n 3 imp ). The angles of spins that do not belong to impurity bonds are again given by Eq. (3.27), with ∆φr replacing ∆θr. However, minimizing over the twist rate Q after linearization with respect to Q of the second line on the right-hand side of Eq. (3.32a), yields the nontrivial saddle point value We can now verify a posteriori the validity of the assumption (3.33). The maximal value of |Q| is given by Thus, for our assumption is certainly self-consistent. In the opposite regime, as n imp approaches 1 from below, the interaction between the planes starts to be dominated by the impurity bonds, which may induce an entirely different ground state with no spiral order. Injecting the saddle point value of Q, Eq. (3.35), into Eq. (3.32) and expressing the energy as a function of the Ising variables σr leads to the same effective Hamiltonian as in Eq. (3.29b), except for an additional term − |Λ|, which expresses the lowering of the total energy due to the coupling of the canting pattern to the spiral order, where we have introduced the constant and the effective Ising interaction Note that the last term in (3.38a) is non-extensive and thus irrelevant in the thermodynamic limit. There are two additive contributions to the Ising exchange coupling J (I) r−r . The contribution Γ (0) r−r represents an (anti)-dipolar two-body interaction between the effective Ising degrees of freedom σr associated with the dilute antiferromagnetic bonds. This interaction is long ranged and frustrated, owing to the indefinite sign of the kernel Γ (0) r−r . The coupling of the canting pattern to the spiral order instead favors a net (ferromagnetic) bias of the cantings σr and contributes an all-to-all interaction of strength 1 J ⊥ |Λ| , proportional to the inverse volume. According to the saddle-point equation (3.27a), a uniform magnetization of the Ising spins σr favors a spiral state with a nonvanishing Q for the original O(2) spin degrees of freedom. Hence, if the ground state of Eq. (3.38) supports a non-vanishing magnetization σr, the frustration induced by the dilute impurity bonds turns the pristine ferromagnetic order of the impurity-free ground state into spiral order. If instead the ground state supports no net bias of the Ising spins σr, a state with Q = 0 is favored, and thus no net winding of the O(2) spins is induced. An example of a possible resulting ground state is a fan-like magnetic state where the Ising spins order in a layered antiferromagnetic pattern. This translates into a pattern of the original O(2) spins ordering essentially ferromagnetically, but with orientations that alternate slightly between successive layers.
The Ising Hamiltonian with anti-dipolar coupling (3.38) is structurally very similar to Ising systems with standard dipolar couplings, as are realized, e.g., in rare earth compounds with strongly localized magnetic moments, such as LiHo x Y 1−x F 4 at moderate to low dilution x 24 . It has been theoretically and experimentally well established that such random dipolar Ising systems exhibit a glass transition towards an amorphous magnetic order at low temperature 19,25 . It may thus come as a surprise that in our case we will find that a change of sign of the dipolar term, in conjunction with the additional term arising from the coupling to the spiral, suffices to induce ferromagnetic Ising order, in spite of the positional randomness of the Ising spins. This difference is presumably largely due to the additional mean-field like interaction mediated by the formation of the spiral, which stabilizes the ferromagnetic phase. That type of interaction is absent in systems of elementary magnetic dipoles, which therefore fall much more easily into a glassy phase.

IV. SUPERLATTICES OF IMPURITY BONDS
As Eq. (3.38) involves long-range two-body interactions whose sign depends on the relative positions of the impurity bonds, the ground state cannot be found explicitly for an arbitrary choice of L, so that in general one has to resort to numerical methods, or to approximate treatments.
However, if the set of impurity bonds L realizes certain Bravais superlattices, it is possible to establish a sufficient condition for the ground state of the effective Ising Hamiltonian (3.38) to be ferromagnetic and, thus, for the ground-state spin configuration of the Hamiltonian (3.1a) to sustain spiral order.

A. Analytical considerations
We consider the case when the subset L of the cubic host lattice Λ forms a Bravais lattice with the basis vectors A, B, and C given by three independent linear combinations with integer-valued coefficients of a ≡ (1, 0, 0) T , b ≡ (0, 1, 0) T , and c ≡ (0, 0, 1) T . The concentration of the impurity bonds is (4.1) In reciprocal space, the superlattice L defines a small first Brillouin zone BZ(L), which is 1/n imp times smaller than the first Brillouin zone BZ(Λ) of the cubic host lattice Λ. For the |L| Ising degrees of freedom σr associated with impurity bonds anchored atr ∈ L, we use the Fourier representation where BZ(L) is the first Brillouin zone associated with the lattice L. For any p ∈ R 3 , we shall make use of the identity where L denotes the reciprocal lattice of L. The effective Ising Hamiltonian (3.38) can now be written in reciprocal space as where γ is defined by Eq. (3.38b) and Using Eq. (4.3), we obtain For a generic choice of the Bravais lattice L and of the couplings J ⊥ , J , and J imp , the ground state of the Ising Hamiltonian (4.4) cannot be found in closed form. However, an analytical solution is available in certain cases. For instance, if over the reduced Brillouin zone BZ(L), the kernel Υ q assumes its global minimum at a unique momentum q min , such that for allr ∈ L 1 − e −ik z e ik·r 2J 2 − cos k x − cos k y + 2J ⊥ (1 − cos k z ) , (4.6b) with γ defined by Eq. (3.38b) and Examples of q min for which Eq. (4.5) holds are q min = 0, (4.7a) and where A , B , and C are the basis vectors of the reciprocal lattice L . In the former case (4.7a), Q min = 0 and the ground state is an O(2) magnetic spiral, while in the latter case (4.7b), Q min = 0 and thus there is no spiral.

B. Comparison between analytical and numerical results for superlattices
To illustrate that the effective Ising Hamiltonian (3.38) captures the low-energy physics of the microscopic Hamiltonian (2.1), we consider several superlattices L of impurity bonds and compare their microscopic ground state to the ground state of the effective Ising Hamiltonian (4.4).
Instead of directly studying the ground state of the microscopic Hamiltonian (2.1), we actually study the microscopic Hamiltonian To obtain the corresponding effective model in terms of cantings we take the following steps. First, we verify that the exchange couplings allow for a non-trivial solution ∆θ = 0 of Eq. (3.15). Second, we solve for the ground state of the effective Ising Hamiltonian (4.4). To this end we minimize the function Υ q in Eq. (4.4c) with respect to q. If the absolute minimum in the Brillouin zone BZ(L) occurs at q = 0, we predict a spiral ground state (Q = 0). If the absolute minimum occurs at q min = C /2, a ground state with Q = 0 is predicted. If instead q min does not satisfy Eq. (4.5), the ground state has more than one Fourier component and we would need to solve the effective Ising model numerically. Third, from the Ising ground state of Hamiltonian (4.4), the microscopic pattern (2.3) of the O(2) spins is finally obtained from Eq. (4.6), using the value of ∆θ obtained from solving Eq. (3.15).
The effective Ising Hamiltonian (4.4) is found to be very accurate once the impurity coupling |J imp | sufficiently exceeds the critical value J c . This is illustrated by Fig. 5. Its four panels compare the approximate ground state obtained via the effective Ising Hamiltonian (4.4a) (shown on the left) with the ground state of the Hamiltonian (4.8) obtained via MC simulation (shown on the right). This is done for two strengths of impurity couplings and two different superlattices. We choose parameters such that both methods yield a spiral state, with q min = 0 minimizing the kernel Υ q . No coupling anisotropy (J /J ⊥ = 1) was assumed in all these cases.  (4.10) When |J imp |/J c is large, this renormalization has little effect on the solution of the saddle point equation, ∆θ, which will be close to π in any case. However, when |J imp |/J c is close to the threshold of 1, an effective reduction of J c (which is expected for superlattice that favor ferromagnetic Ising order) leads to an increase of ∆θ. For example, in Fig. 5(a), the value of the canting angle would ∆θ = 1.02 according to Eq. (3.15), but increaes to ∆θ = 1.15 after correcting it by Eq. (4.9).
Similar results are found for other superlattices. For instance, panels (c) and (d) show results for a denser superlattice, but with the same exchange couplings as in panels (a) and (b), respectively. In all panels of Fig.  5, the deviations from the local ferromagnetic order at non-impurity bonds are small, justifying a posteriori the spin-wave approximation used to derive the effective Ising Hamiltonian (4.4).
To quantify the quality of the approximations incurred when trading the microscopic Hamiltonian (4.8) for the effective Ising Hamiltonian (4.4), we compare the quantity (4.11) obtained from both Hamiltonians for several superlattices and various ratios |J imp |/J ⊥ in Fig. 6(a). Here, L x , L y , and L z in Eq. (4.11) are the linear dimensions of the lattice, while P is an order parameter for the magnetic spiral phase. On the right-hand side, the sine of the relative angle between S r and S r+z is summed over all sites of the cubic host lattice Λ.
Figure6(a) shows how the value of |P |, evaluated on the minimal energy configuration, increases with increasing |J imp |/J ⊥ ≥ J c /J ⊥ for three superlattices of impurity bonds in an isotropic cubic lattice (J /J ⊥ = 1), whereby all superlattices were chosen so that they induce a spiral state. At relatively large |J imp |/J ⊥ , the results for |P | from the effective Ising Hamiltonian (4.4) (dots) are close to those obtained from the microscopic simulation of Hamiltonian (4.8) (squares), up to corrections of order n imp , as anticipated in the discussion around Eq. (3.28). However, as |J imp | approaches J c from above, deviations become stronger, as we discussed after Eq. (4.9). In this regime the double-well potential defining the Ising degrees of freedom associated with the canting pattern becomes very shallow. Thus, even relatively weak contributions from neighboring impurities, Ξr, can stabilize and enhance the local canting ∆θ and strengthen the spiral wave vector beyond the approximations we used to derive the effective model. For the same reason of mutual stabilization, we still find a finite spiral order, P = 0, even when J imp /J ⊥ J c /J ⊥ = 2.

C. Spiral phase in a realistic model for YBaCuFeO5
It was argued in Ref. 18 that the magnetic degrees of freedom in the insulator YBaCuFeO 5 realize a close cousin of Hamiltonian (4.8), in that J ⊥ > 0 is to be replaced by two distinct values J ⊥ > 0 and J ⊥ > 0 depending on the parity of the z component of the coordinate r of the cubic lattice. From the estimates for J > 0, J ⊥ > 0, J ⊥ > 0, and J imp < 0 made in Ref. 18 we have borrowed the values J = 28.9meV, (4.12) Figure 6(b) compares the dependence of the magnitude |P | of the spiral order parameter P defined in Eq. (4.11) on J imp for the microscopic Hamiltonian (4.8) (squares) with that for the effective Ising Hamiltonian (4.4) (dots) for the case when the impurity bonds form a superlattice that stabilizes a long-range spiral order. Again, good agreement is found once the the impurity bond strength is sufficiently stronger than the threshold, in which case the two possible canting patterns form robust local minima of the Hamiltonian. This is indeed the case in YBaCuFeO 5 , where |J imp |/J ⊥ ≈ 23 > J c /J ⊥ ≈ 9.4.

D. Dependence of the ground state on the superlattice of impurity bonds
We now illustrate how the choice made for the superlattice of impurity bonds affects the ground state. We use again the values (4.12) corresponding to ideal YBaCuFeO 5 when defining the effective Ising Hamiltonian (4.4) and the microscopic Hamiltonian (4.8).
We consider two superlattices of impurity bonds. They are chosen such that Υ q has a global minimum at q min = 0 for one, and at q min = C /2 for the other, cf. Fig. 7). Furthermore, both superlattices are chosen such that they share with YBaCuFeO 5 the additional property that impurity bonds only occur between every other plane.
The essential difference between the two superlattices lies in the relative position of nearest-neighbor impurity bonds. In the first lattice, the majority of nearestneighbor impurity bonds is ferromagnetically coupled. This favors a ferromagnetic Ising phase in the effective Ising Hamiltonian (4.4), i.e., a spiral magnetic phase. In the second lattice, the majority of nearest-neighbor impurity bonds is antiferromagnetically coupled. This favors a layered antiferromagnetic Ising ground state of the effective Ising Hamiltonian (4.4), and thus, a fan-like magnetic order with no net winding of the spins, whereby the orientation of the magnetization of the layers alternate between even and odd pairs of planes.
We have verified using MC simulations of the microscopic Hamiltonian (4.8) that the long-range spiral order, which is present when the impurity bonds are arranged in certain Bravais superlattices, is robust to weak distortions of that superlattice, as expected on theoretical grounds.

E.
Limit of dilute impurities Next we analyze the limits of very dilute tetragonal, face-centered-cubic, and body-centered-cubic superlattices of impurity bonds. These are tractable analytically. We shall show that ferromagnetic order prevails at the level of the Ising degrees of freedom associated with local cantings for dilute face-centered-cubic and bodycentered-cubic superlattices of impurity bonds. This is to say that spiral order for the underlying XY spin degrees of freedom prevails for these diluted superlattices of impurity bonds. The results obtained here will also be useful for the study of random impurity bonds in Sec. V.

Cubic superlattices
We assume that the impurity bonds occupy a cubic sublattice L of the cubic host lattice Λ. As it turns out, this case supports antiferromagnetic order in the Ising model. Some results of this calculation will later help us to establish that the disordered case, in contrast, orders ferromagnetically.
If the cubic host lattice Λ and the superlattice L are finite and not too large, it is possible to calculate the energy (3.38a) for all Ising spin configurations by exact evaluation of the Ising kernel J (I) r−r defined in Eq. (3.38c). In the thermodynamic limit, |Λ| → ∞, with n imp held fixed, this approach is not possible anymore. Instead, we shall restrict ourselves to a few long range-ordered Ising configurations that are likely candidates for the ground state, and compare their energies.
The ferromagnetic Ising configuration is described by which describes a sequence of stacks of m ≥ 1 layers, whose magnetization alternates. Here, where the spin autocorrelation function f C (r) ≡ σr σr +r r = f C (z) (4.14b) only depends on the difference in thez coordinate, owing to Eq. (4.13). Here, . . . r denotes the average over the sitesr of the superlattice L. For configurations F and AF(1), it is given by (4.14c) In the dilute limit n imp → 0, the typical distance between a pair of nearest-neighbor impurities is large. Hence, the typical pair-wise interaction Γ (0) r tends to the dipolar form (2.14) and can be safely used to evaluate ε AF(m) L up to corrections which are subleading in the limit n imp → 0. The case of the ferromagnetic configuration is more subtle, however. Indeed, a naive use of Eq. (2.14) would suggest that the first term in the right-hand side of Eq. (4.14a) vanishes due to the sum over symmetry related directions, while in fact it does not. This is due to corrections to the dipolar interaction (2.14) that scale as the inverse of the volume, but add up to a finite contribution when summed with equal signs over the whole superlattice. In the case of an isotropically shaped, cubic sample with L x = L y = L z and isotropic interactions J = J ⊥ ≡ J, the computation can be done exactly, using the fact that upon averaging over all the permutations k x → k y → k z → k x the kernelΓ (0) k (2.13b) reduces to 1/3J. Using f F (z) = 1 this allows us to evaluate the lattice sum exactly for any impurity density as This finite, negative contribution disfavors the ferromagnet, in analogy to demagnetizing factors known from standard magnetic dipolar systems. Restricting ourselves to the isotropic case and inserting Eq. (4.15) into Eq. (4.14a), we obtain the energy per impurity (4. 16) More generally, it is useful to cast the energy (4.14a) for the ferromagnetic configuration (4.13a) of the Ising variables in a different form, namely, where the reciprocal lattice L of L enters through the identity (4.3),Γ k =0 has been defined in Eq. (2.13b) and we recall thatΓ (0) k=0 = 0. The sum over k thus contains |Λ|/|L| − 1 = 1/n imp − 1 terms. Note that the self-interaction Γ (0) r=0 [cf. Eq. 3. 18)], which appears also in the single impurity energy, is subtracted on the righthand side of Eq. (4.17).
We point out an important difference between the present effective dipolar problem and genuine magnetic dipoles. Genuine dipolar interactions are mediated by magnetic fields which extend everywhere in space, beyond the boundaries of the sample. Therefore they only depend on the relative position of two spins, irrespective of where the spins are deep in the bulk, or close to a surface of a finite sample. However, this is not so in our case where the dipolar interactions arise through the mediation of spin waves, which are confined to the sample. Accordingly, the interactions involving Ising spins at the periphery of the sample are not exactly the same as those for bulk Ising spins with the same relative position. More importantly there are no magnetic stray fields beyond the sample. In real dipolar magnets those store a lot of magnetic energy, which is avoided in the ground state by domain formation. The unavoidable presence of domains complicates the computation of the energy density. In particular, the evaluation for a homogeneously magnetized sample yields a shape dependent result, a fact that is reflected in the ambiguity of the value of the Fourier transform of the dipolar interactions Eq. (2.13b) in the limit k → 0. In the present case, however, such problems do not arise, since the spin-wave mediated interaction is such thatΓ (0) k=0 = 0. This eliminates the potential ambiguity and therefore eliminates the shape dependence. We also do not expect the effective dipolar interactions to induce domains, in contrast to genuine ferromagnets.
Performing an analogous calculation to the one above yields for the antiferromagnet AF(1) Even though the quantitative mapping from the XY Hamiltonian (2.1) to the effective Ising Hamiltonian (3.38) only holds for low densities of impurity bonds, it is useful to study the effective Ising Hamiltonian (3.38) in its own right, i.e., without requiring the impurity bonds to be dilute.
A maximally dense superlattice is defined by For such a lattice, one finds the ferromagnetic (F) and antiferromagnetic [AF(1)] states to be degenerate, since, cf. Eqs. (2.13b) and (2.13c), The identity (4.21) obeyed by the kernel (2.13b) can be used together with the expression (4.4a) and the fact that only q of the form (0, 0, k z ) T enter it, to show that for a maximally dense superlattice all antiferromagnetic states AF(m) are degenerate with the ferromagnet. More generally, it is shown in appendix A that the ferromagnet is degenerate with any Ising configuration in which the spins in every given plane at fixed z coordinate are ferromagnetically aligned, irrespective of the relative orientation of the magnetization of different planes. This degeneracy is, however, lifted at finite dilution, whereby the way in which the dilution is realized is crucial. For example, diluting the impurity density n imp by maintaining a cubic superlattice, but increasing its integer lattice spacing disfavors the ferromagnetic state. This is illustrated in Fig. 8, where we plot the energies per impurity as a function of superlattice spacing . For small , the energy difference is obtained from the representations (4.17) and (4.18). In the dilute limit, n imp = −3 → 0, the reciprocal lattice L only contains small wavevectors, and we may replace 1 − cos k i (where k i = 2π n i / ) in the kernel (2.13b) by (2π n i ) 2 /2 2 with n i ∈ Z for i = x, y, z, i.e., where the exchange anisotropy parameter α was defined in Eq. (3.21), and (4.22b) The sum over n z can be carried out explicitly, 2π α (n 2 x + n 2 y ) sinh 2π α(n 2 x + n 2 y ) . (4.23) Hence, δ(α) is always positive. For the isotropic limit α = 1, one finds δ(1) ≈ 0.1042. Alternatively, one can calculate the antiferromagnetic energy directly in real space using the dipolar form (2.14). This can be used to calculate the energies of other antiferromagnetic states AF(m), which all scale as (4.24) From the results (4.16, 4.22a) it follows that c 1 = δ(1) + 2/3, while one finds the higher c m 's to decrease monotonically with increasing m. From this we conclude that a dilute cubic superlattice orders antiferromagnetically with layer magnetizations that alternate in sign (m = 1).

Dilute tetragonal, face-centered, and body-centered tetragonal superlattices
One readily generalizes the above calculation to tetragonal superlattices L with unit vectors (A , 0, 0) T , (0, A , 0) T , (0, 0, C ) T , where A and C are fixed integers while the integer-valued dilution parameter will be taken to infinity. This case is obtained from that of a cubic lattice by substituting 25) in Eq. (4.22a) and (4.22b). Independently of the ratio C/A of the tetragonal superlattice, the Ising antiferromagnetic state AF(1) is favored over the Ising ferromagnetic state F. However, similarly as in lattice problems of physical electric or magnetic dipoles 26 where the interactions have reversed global sign, a different ground state is found in dilute body-centered or face-centered tetragonal lattices. The difference arises because closest neighbors in these lattices have a stronger tendency to have ferromagnetic interactions than in simple tetragonal lattices. For the face-centered tetragonal lattice, the basis vectors are (A, A, 0), (A, 0, C), and (0, A, C). The corresponding dual basis vectors in reciprocal space are e 1 = π(1/A, 1/A, −1/C), e 2 = π(1/A, −1/A, 1/C), and e 3 = π(−1/A, 1/A, 1/C). Their linear combinations with integer coefficients span the reciprocal lattice L . It is convenient to represent a generic reciprocal lattice vector G ∈ L as G = n 1 e 1 + n 2 e 2 + n 3 (e 2 + e 3 ). With this choice, the asymptotic energy difference between the ferromagnetic and the antiferromagnetic states in the infinite dilution limit n imp → 0 can be written as (4.26b) where for vanishing wavevector we have to set g 0 = 0. Carrying out the sum over n 3 one finds (−1) n 1 −n 2 J ⊥ × π 2α (n 2 1 + n 2 2 ) sinh π 2α (n 2 1 + n 2 2 ) .
For body-centered tetragonal lattices one finds the same expression, but with the replacement 2α → α. The energy difference turns out to be always negative for any value of α, as seen in Fig. 9. Thus, in both these types of superlattices the ferromagnetic state is favored over the layered antiferromagnetic state, whatever the tetragonal aspect ratio.

V. RANDOM IMPURITIES: DILUTE LIMIT
In this section, we study randomly distributed impurities that occupy a fraction n imp of the sites of the cubic host lattice Λ. We assume again that the relevant contenders for the ground state are given by Eqs. (4.13a) and (4.13b). In Eq. (4.13b), we set = 1, since only the lattice constant of the cubic host lattice Λ is relevant. These configurations are expected to come reasonably close to the true ground state and the relevant competing metastable configurations. However, they will differ Independently of the value of α, the antiferromagnetic state has lower energy for simple cubic superlattices while the ferromagnetic state has lower energy for face and body centered superlattices.
in the orientation of a few spins from the simple configurations (4.13a) and (4.13b). The relative fraction of these spins becomes increasingly small as n imp → 0, as discussed below.
If the impurities are distributed randomly according to a Poisson process, the average energy per impurity bond of the trial configurations C = F, AF(m) is given by since any site r of the cubic host lattice Λ is the lower end of an impurity bond with probability n imp , independently of the location of other impurities. From this observation, one might at first conclude that the antiferromagnetic state should dominate again. However, the above consideration does not treat correctly impurities located at short distances from each other. On the one hand, rare pairs of impurities that are located much closer to each other than the average separation n −1/3 imp do not follow the pattern (4.13a) and (4.13b), but simply optimize their mutual interaction energy, irrespective of the global ordering pattern. Since such pairs nevertheless contribute a finite fraction to the total energy estimated above, they must be corrected for, which will turn out to favor the ferromagnetic ordering. This conclusion will become clear below, as a corollary to the discussion of another shortdistance effect, which we will consider first.
Impurity distributions in real materials are usually not simply governed by a Poisson process, but rather, one should expect them to exhibit some short-range correlations. For example, in the case of YBaCuFeO 5 impurity bonds arise due to chemical disorder which occasionally replaces the usual Cu-Fe pairs on bonds along its crystallographic c-axis by impurity configurations consisting in Fe-Fe or Cu-Cu pairs. Fe-Fe pairs differ from Fe-Cu pairs by the sign and magnitude of the resulting magnetic exchange constant. Moreover, both Fe-Fe and Cu-Cu pairs differ from Fe-Cu pairs in their local charge density. The resulting Coulomb repulsion between such impurity configurations thus suppresses the occurence of pairs of impurities at short distances. In a crude manner, we can mimic this effect by a hard constraint on the minimal distance between impurities, excluding distance vectors with |r| ≤ R. With such a constraint the average energy per impurity (5.1) is modified to (5.2) Note that for R = 0 these energies are simply n imp multiplying the energy per impurity ε C L=Λ (R) of a maximally dense system of impurities, recall Eq. (4.14a). As we have shown in the previous section, those energies are all degenerate. Since the sum over r in Eq. (5.2) is dominated by small |r|, even a small R of the order of one lattice constant will have a decisive effect and lifts this degeneracy. In Fig. 10, we plot as a function of R the average energies 2ε F dis (R)/(γ n imp ) and 2ε AF(m=1) dis (R)/(γ n imp ) of the two most relevant competing states. Already, for the smallest effective exclusion radius of R ≥ R c = 1 (in units of the host cubic lattice spacing), we find that the ferromagnetic state (and thus XY spiral order) wins over the antiferromagnetic state (i.e., XY fan order). This numerical result can be un-derstood by recalling that ε F and ε AF(1) are degenerate for R = 0. Upon barring impurities on nearest-neighbor sites on the host cubic lattice, the two states receive a relative energy shift 4n imp Γ r=z = 4n imp × ( A J ), which stabilizes the ferromagnetic state (A ≈ 0.123). Larger exclusion radii tend to reinforce this trend, as shown in Fig. 10. In the limit of large R, the energy per impurity bond of the ferromagnetic state is more favorable than that of the antiferromagnetic one by γ n imp /(3J) in the case of isotropic couplings. This can be understood as follows. For isotropic couplings, the ferromagnetic energy per bond, ε F = −γ n imp /(3J), remains unchanged upon excluding the interactions with a set of sites that is invariant under the cubic symmetry group, as seen in Fig. 10. In contrast, in an antiferromagnetic state, the interactions with the neighbors in thin spherical shells of approximately fixed radius r > R) come with alternating signs. Those tend to cancel the more effectively the larger is R, such that ε AF(1) /(γ n imp ) → 0 as R → ∞.
Even without any repulsive short-range correlations between impurity locations, one expects Ising ferromagnetism to prevail at sufficiently low impurity densities. This is because rare impurities with a neighboring impurity much closer than n −1/3 imp should effectively be taken out of the calculation for the average energy. Indeed, if the close pair is antiferromagnetically coupled, it will anti-align, have no net moment and thus essentially decouples from the global ordering pattern. If instead the pair is ferromagnetically coupled, it forms a bigger spin that can then be incorporated in the consideration like any other typical spin. The net effect of treating such close pairs in this way boils down to considering only original or effective spins with pairwise separations of the order of R eff c n −1/3 imp with some constant c of order 1. The competition for the global ordering pattern then becomes essentially identical to the one of the constrained superlattice above, with R eff now taking the role of the exclusion radius in Eq. (5.2). From these considerations we predict that for sufficiently dilute concentrations n imp (c/R) 3 c the Ising ferromagnetic order prevails.

VI. FINITE-TEMPERATURE TRANSITION TO THE SPIRAL PHASE
The effective Ising model (3.38) undergoes ordering at a critical temperature T Ising ∝ n imp 27 . As long as T Ising lies in the range of low temperatures (3.16), the Ising approximation is well justified. This is certainly the case for n imp 1. Since the reduction to the Ising model neglects some fluctuations, we expect T Ising to be an upper bound to the actual spiral transition temperature T spi . However, the bound should become increasingly tight as the impurity concentration decreases towards n imp → 0.

B. Dipolar approximation
The mean-field theory of the previous section has at least two drawbacks. As usual, the neglect of fluctuations will lead to an overestimate of the critical temperature by a certain factor O(1), which might itself be a function of the ratios between the couplings. This makes it difficult to predict the precise dependence of T Ising on the couplings. A second and more serious drawback of these approximations is the fact that the site-averaged mean field of Eq. (6.4) receives rare, but large contributions from pairs of sites that are nearest neighbors on the underlying lattice Λ. This contribution represents a nonvanishing fraction of the resulting mean field. However, physically it is clear that the Ising spins on very close pairs of sites will lock strongly together and act either as an effective spin with a doubled moment for ferromagnetic pairs, or they essentially decouple from the rest for antiferromagnetic coupling. In either case, these strong short-range couplings have essentially no influence on the long-range ordering, and thus it seems unphysical that such strong couplings should enter in our mean-field estimate of T Ising at all.
Here, we follow a different approach to establish the dependence of T Ising on the couplings. Let ξ be the length scale beyond which we can approximate the interactions J (I) r−r as being anti-dipolar, i.e., given by Eq. (2.14). We assume that we can safely neglect pairs of Ising spins that are within a distance of order ξ of each other [the probability to find another Ising spin a distance ξ from a given one is of order O(ξ d n imp ), a negligible probability as n imp → 0]. If so, we may replace the Ising Hamiltonian The parameter density, upon using the couplings given in Eq. (4.12), see Fig. 1. It is encouraging that our theory overestimates T spi /|Q| only by ≈ 13%, considering the simplifications that go into the modelling of the spin system and the uncertainty in the value of the exchange couplings. As noted above, the ordering temperature is proportional to the concentration of impurity bonds. For the concrete case of YBaCuFeO 5 , our theory predicts that a small fraction of 1% of the oxygen bipyramids realizing the strongly frustrating Fe-Fe magnetic interactions induces a transition to the spiral phase at an ordering temperature of approximately 85 K, see Fig. 1. Note that at those temperatures the constraint of Eq. (3.16) is satisfied by a large margin and the mapping to the effective Ising model is thus well controlled. A subsequent experimental study investigated a family of chemically modified compounds related to YBaCuFeO 5 33 . In those materials the lattice parameters could be altered, which affects the exchange constants, and the concomitant changes to observables such as T spi /|Q| were recorded. Examples of such modifications are the application of uniaxial pressure, or chemical substitution that replaces the atoms between the layers containing the impurity bonds. The latter modifies the interlayer spacing and thus the perpendicular coupling J ⊥ . The experiments of Ref. 33 shows that the ratio T spi /|Q| is only very weakly sensitive to the modification of the interlayer spacing and thus J ⊥ .
These empirical findings can be rationalized by analyzing Eqs. (6.10) and (6.15). In layered materials such as YBaCuFeO 5 , the exchange anisotropy between intra-and inter-layer couplings is large. We model this empirical fact by requiring that α ≡ J /J ⊥ 1 [recall Eq. (3. 21)]. Furthermore, the impurity coupling strength turns out to be large as well, |J imp |/J 1 [recall Eq. (4.12)]. In this limit, J c ≈ 2πJ /(ln α + 2.47) J ⊥ [recall the approximation mentioned in the caption of Fig. 3 ]. The canting angle between the XY spins on either end of an impurity bond comes close to ∆θ ≈ π [recall Eq. (3.15b)]. More precisely, the deviation from π is After dropping this correction, to a first approximation, the ratio between the critical temperature and the spiral twist rate at zero temperature can be approximated by

VII. CONCLUSION AND OUTLOOK
Any three-dimensional lattice hosting XY spins that interact through ferromagnetic nearest-neighbor exchange interactions display a ferromagnetic long-range order below some critical temperature. We have given sufficient conditions under which the replacement of a dilute fraction of the ferromagnetic bonds by antiferromagnetic bonds destabilizes the ferromagnetic order in favor of non-collinear long-range order in the form of a spiral phase. A necessary but not sufficient condition for spiral order is that the antiferromagnetic exchanges along the impurity bonds be sufficiently larger than the ferromagnetic couplings. This induces local canting, which lowers the energy close to the frustrating bond. If this condition is met, a sufficient condition for spiral order is a strong correlation between the impurity bonds such that (i) they all point along a preferred direction and (ii) they are distributed in space such that ferromagnetic interactions dominate between the Ising degrees of freedom associated with the local canting patterns around the impurities. We showed rigorously that (ii) is satisfied for impurities located on Bravais superlattices whose shortest lattice vectors tend to point in directions in which the effective Ising interactions are ferromagnetic, while neighboring impurities along the z-axis, for which the interactions are antiferromagnetic, appear only at larger distance. Small distortions away from a perfectly regular Bravais lattice will not destroy the spiral order. We also argued that completely randomly distributed impurities are prone to stabilize spiral order at low enough impurity density. At higher impurity densities, a short-ranged repulsion among impurity bonds, e.g., due to Coulomb constraints in real materials, has the main effect of reducing the stability of fan states (layered antiferromagnetic orderings of the canting degrees of freedom), and thus also stabilizes spiral order. Hence, once the orientational correlation (i) is ensured, the tendency towards spiral order is rather strong.
On the other hand, if the impurity bonds and their orientations are white-noise correlated in space, the microscopic XY Hamiltonian belongs to the family of three-dimensional XY gauge glasses introduced by Villain. Those host amorphous, glassy order. From this it follows that the zero-temperature phase diagram of two-dimensional XY magnets (as characterized by the strength of the frustrating antiferromagnetic interactions and their spatial correlations) contains at least four stable phases: The ferromagnetic phase, the spiral phase, the fan phase (i.e., ferromagnetic in plane order with ori-entation oscillating from plane to plane), and the gauge glass phase.
From the perspective of the original microscopic XY spins in Hamiltonian (2.1), the phenomenology for small concentrations n imp 1 is the following. Upon lowering the temperature in the XY paramagnetic phase, a continuous phase transition takes place in the threedimensional XY universality class to a ferromagnetic phase at the temperature T XY . This ferromagnetic phase becomes further unstable at the temperature T spi T XY (as estimated by T Ising in Eq. (6.12)), where a XY spiral phase emerges via a continuous phase transition. It is driven by the dilute concentration n imp 1 of impurity bonds that are orientationally correlated. The spiral wavevector Q may serve as an order parameter for this Ising transition. The associated critical exponents are expected to take mean-field values, given the dimensionality and the long-range nature of the dipolar interactions.
What happens as n imp is increased, so that T MF spi ∼ T XY ? In this limit, the effective Ising model (3.38) is not a valid approximation of Hamiltonian (2.1) anymore, so that at this stage we cannot make controlled predictions. However, it seems very likely that at large enough n imp 1, the impurity bonds will dominate the coupling between adjacent a, b-planes, inducing a layered antiferromagnetic state. Upon increasing n imp this state might be reached either continuously, with the spiral wavevector saturating at Q = π, or it appears discontinuously, via a first order transition at some critical value of n imp . In the regime of smaller n imp 1, the critical temperature of the ferromagnetic instability will decrease with increasing n imp , while the spiral instability temperature is expected to continue to increase. They might merge into a single direct transition, if this is not pre-empted by the emergence of a layered antiferromagnetic phase. An approach that is non-perturbative in n imp is needed to address these questions. It remains an open challenge to determine optimal combinations of exchange couplings that would allow to maximize T spi by increasing n imp , and thereby extend the regime of the incommensurate spiral magnetic phase to high temperatures.
In our earlier companion paper Ref. 18, it was argued that YBaCuFeO 5 unites all the essential ingredients of the Hamiltonian discussed in this work, and thus could realize the spiral XY phase described above. The supporting evidence is as follows. On the one hand, Monte Carlo simulations for realistic values of the magnetic exchange couplings in YBaCuFeO 5 yield transition temperatures to the magnetic spiral phase as high as 250 K. On the other hand, it was reported in Ref. 13 that tuning the degree of occupational disorder by changing the annealing procedure of YBaCuFeO 5 affects the transition temperature and the wavevector of the spiral in a way that is qualitatively and quantitatively consistent with Eq. (6.15). Finally, we point out that a mechanism very similar to the one described here might be at work in hole-doped cuprates, where pairs of holes might take the role of the frustrating impurity bonds 36 .

A. Applications to other systems
The main physical mechanism we discussed in this work applies to other systems as well. First, we point out that the restriction to XY spins is not essential. Indeed, we expect that Heisenberg spins with an O(3) symmetry (or any other set of continuous degrees of freedom undergoing spontaneous symmetry breaking) would exhibit essentially the same phenomenology: At low temperatures the unfrustrated system will order ferromagnetically. Frustrating antiferromagnetic impurity bonds induce local canting patterns that are subject to effective pairwise interactions upon integrating out spin waves. A ferromagnetic order of the canting degrees of freedom again imply spiral order for the original Heisenberg spins. If the canting induced by a local impurity bond preserves the coplanarity of the background ferromagnetic order, the problem simply reduces to an effective XY model. This is what we found to happen in the presence of nearest neighbor Heisenberg interactions. However, with more complex interactions, it might occur that the local canting pattern is non-planar. This would imply that the canting does not only have a discrete Ising degree of freedom, but rather a continuous XY −like degree of freedom. Indeed, for an isolated impurity, any rotation of all spins around the direction of the background ferromagnetic magnetization yields an energetically equivalent canting pattern. Upon integrating out spin waves, these effective XY canting degrees of freedom will be coupled through dipole-like interactions, and their ferromagnetic order will again induce a spiral of the original Heisenberg spins.
The phenomenology of magnetic XY spins immediately carries over to superconducting systems, too. There, the role of XY spins is taken by the phase of superconducting islands with a well established amplitude of the superconducting order parameter, and Josephson couplings replace the magnetic exchange couplings. Frustration could be induced by Josephson couplings with a negative sign (based on ferromagnetic materials for example). However, a much simpler way to achieve frustration consists in threading a homogenous magnetic flux through a Josephson junction array. The recent advances in fabrication techniques and nanolithography for such devices should allow to artificially design and control XY systems with a desired spatial pattern of frustrated plaquettes that emulate the presence of the antiferromagnetic impurity bonds in the magnetic analogue. A magnetic spiral phase with ferromagnetic order of the Ising degrees of freedom of the canting patterns then translates into a system of vortices of the same vorticity (sense of circulation), entailing a global supercurrent in the system. This will be explored in future work.