Multiferroic magnetic spirals induced by random magnetic exchanges

Multiferroism can originate from the breaking of inversion symmetry caused by magnetic-spiral order. The usual mechanism for stabilizing a magnetic spiral is competition between magnetic exchange interactions differing by their range and sign, such as nearest-neighbor and next-nearest- neighbor interactions. Since the latter are usually weak the onset temperatures for multiferroism via this mechanism are typically low. By considering a realistic model for YBaCuFeO$_5$ we propose an alternative mechanism for magnetic-spiral order, and hence for multiferroism, that occurs at much higher temperatures. We show using Monte-Carlo simulations and electronic structure calculations based on density functional theory that the Heisenberg model on a geometrically non-frustrated lattice with only nearest-neighbor interactions can have a spiral phase up to high temperature when frustrating bonds are introduced randomly along a single crystallographic direction as caused, e.g., by a particular type of chemical disorder. This long-range correlated pattern of frustration avoids ferroelectrically inactive spin glass order. Finally, we provide an intuitive explanation for this mechanism and discuss its generalization to other materials.

Insulators with magnetic spiral order are of particular interest because of their associated multiferroism [1][2][3][4][5] in which the breaking of inversion symmetry by the magnetic spiral drives long-range ferroelectric order. The magnetic order can then be manipulated by an applied voltage with minimal current dissipation due to the insulating nature offering potential for low loss memory devices.
For many insulators, such as the orthorhombic manganites RMnO 3 (R=Dy,Tb) [6][7][8][9][10], spiral order results from a competition between nearest-neighbor and further-neighbor magnetic exchanges of comparable strength. As a consequence the onset temperature T spi is set by the rather low energy scale of a typical further-neighbor exchange, strongly limiting the multiferroic temperature range. Alternative routes to stabilizing magnetic spirals at higher temperatures are, therefore, of fundamental and technological interest.
The phenomenology of the spiral magnet YBaCuFeO 5 , which has one of the highest critical temperatures among the magnetically driven multiferroics [11,12], suggests that a particular type of chemical disorder might provide such a route. As the temperature is lowered below T N ∼ 440 K in YBaCuFeO 5 , the paramagnetic state undergoes a transition to a commensurate magnetic order with wave vector q N = ( 1 2 , 1 2 , 1 2 ). Then, below T spi < T N , a multiferroic magnetic spiral phase sets in, with a propagation wave vector along the c crystallographic axis, q spi = ( 1 2 , 1 2 , 1 2 − Q). The value of Q increases smoothly from Q(T spi ) = 0 as temperature is decreased. Impor-tantly, the reported values of T spi range from 180 K to 310 K [11][12][13][14][15][16] depending on the preparation conditions, and it was recently shown [16] that T spi and Q increase systematically with Fe 3+ /Cu 2+ occupational disorder. These observations suggest that chemical disorder plays an essential role in stabilizing the magnetic spiral motivating our search for a microscopic mechanism by which disorder facilitates, or even drives, magnetic spiral order.
In this paper, we introduce a classical Heisenberg spin model for YBaCuFeO 5 in which spiral order is indeed induced by chemical disorder. We describe the local moments of Cu 2+ and Fe 3+ as classical Heisenberg spins, S r , localized at the positions r of a lattice. We assume only nearest-neighbor exchange interactions, and use the magnitudes calculated from the local spin density approximation including an effective Hubbard U correction (LSDA+U) for YBaCuFeO 5 [17]. Without chemical disorder, the magnet is unfrustrated and establishes commensurate antiferromagnetic order at T N ≈ 300 K. Frustration is introduced through dilute impurity bonds with enhanced exchange couplings of opposite sign. The structure of YBaCuFeO 5 is such that the chemical disorder introduces only collinear impurity bonds parallel to the c axis. We show that the induced frustration results in a local canting of the antiferromagnetic order parameter around the impurity bond, which spontaneously breaks the local inversion symmetry around each impurity so that S r imp +δr = S r imp −δr . The Goldstone modes of the antiferromagnet mediate a long-range coupling be-tween the cantings of different impurities and establish long-range order among the cantings that induces a continuous twist of the antiferromagnetic order parameter in the direction parallel to the impurity bonds. Related mechanisms are likely responsible for the numerical results obtained in the two-dimensional systems studied by Ivanov et al. [18] and by Capati et al. [19]. There, too, impurity bonds were orientationally correlated. We note that randomly oriented impurity bonds would have a spin glass solution [20][21][22][23], whose magnetic order does not couple to a net electric polarization and thus does not lead to multiferroism. This mechanism results in a T spi of the order of a typical exchange coupling. Our Monte Carlo simulations for YBaCuFeO 5 yield T spi as high as 250 K depending on the concentration of the impurity bonds and their strength, in a manner that is consistent with the experimentally observed dependence of T spi and q spi on the amount of Fe 3+ /Cu 2+ occupational disorder [16].
Microscopic origin of spin-spiral state in YBaCuFeO 5 . YBaCuFeO 5 forms a vacancy ordered perovskite structure for which planes of Y ions separate bilayers of BaCuFeO 5 . The latter consist of corner-sharing FeO 5 /CuO 5 bipyramids, as depicted in Fig. 1(a).
A recent ab initio study [12] revealed that the lowest energy arrangements of Fe and Cu ions contain only Fe 3+ -Cu 2+ bipyramids as shown in the left panel of Fig. 1(a). The exchange interactions between the magnetic ions making up these low-energy structures were also calculated using LSDA+U. The exchange interactions within the Fe-Cu bipyramids were found to be uniform in sign (ferromagnetic), and so unfrustrated, and thus cannot explain the emergence of a magnetic spiral. Bipyramids of Cu 2+ -Cu 2+ and Fe 3+ -Fe 3+ , shown in the right panel of Fig. 1(a), are energetically more costly, but they nevertheless occur as local defects that form during preparation and thermal treatment of the sample. The Fe 3+ -Fe 3+ bipyramids introduce strongly frustrating antiferromagnetic couplings along the c axis [24]. In the following, we show that a small concentration of Fe 3+ -Fe 3+ bipyramids is at the microscopic origin of the spiral order.
We model the magnetic ordering of YBaCuFeO 5 using the classical Hamiltonian Here, the classical spin S r is a three-component unit vector located at the site r of a cubic lattice isomorphic to the tetragonal lattice of YBaCuFeO 5 . To distinguish between the lattice vectors of the real material and those of the spin model we label the latter as x, y, and z. Their length is equal to a = b ≈ c/2 where a , b , and c are the lattice vectors of YBaCuFeO 5 as shown in Fig. 1 (with c z). To facilitate comparison with experiments, we express the wave vectors in units of the reciprocal vectors of the crystallographic unit cell of YBaFeCuO 5 throughout. We retain only nearest-neighbor exchange interactions. The ab planes host spins on a square lattice, coupled antiferromagnetically with the strong exchange J = −28.9 meV [25]. These planes are stacked along the c axis in bilayers (Fig. 1). Within a bilayer, the nearest neighbors along the c axis, that is the spins in a bipyramid, are coupled with the weak ferromagnetic exchange interaction J ⊥ = 2.7 meV > 0. Adjacent bilayers are coupled with the relatively weak antiferromagnetic coupling J ⊥ = −5.5 meV. Finally, the defect bipyramids are modelled by a small concentration, n imp , of randomly located "impurity bonds" lying within the bilayers. They substitute J ⊥ > 0 with the strong antiferromagnetic coupling J imp = −95.8 meV calculated for Fe 3+ -Fe 3+ bipyramids and locally frustrate the ferromagnetic intra-bilayer couplings inducing a local canting for |J imp | sufficiently larger than J . At the same time, we assume impurities to be sufficiently dilute, (n imp J imp J ⊥ ), so that intralayer antiferromagnetic alignment does not become favored within bilayers. The impurity bonds are long-range correlated in the sense that they are always oriented parallel to the c axis. We do not include the Cu 2+ -Cu 2+ bipyramids (which stoichiometry implies to be as abundant as the Fe 3+ -Fe 3+ bipyramids) since their inter-layer exchange is substantially smaller than all other couplings.
The second term in Eq. (1) describes a small easy plane anisotropy with ∆ = 0.5 meV which favors magnetic moments in the plane perpendicular to the direction defined by the unit vectorn. The spiral order parameter is defined by where N is the number of spins in the lattice. Neutron scattering shows thatn =ĉ, which ensures a cycloidal component of the spiral and thus a coupling to the electric polarization [26]. For the calculation of the thermodynamics, the orientation ofn plays no role. In our Monte Carlo simulations, we use a cubic superlattice of linear dimension L×L×2L (where L ≤ 28) with periodic boundary conditions in plane and open boundary conditions in the z direction. Impurity bonds are randomly located, but subject to the constraint that any two impurity bonds in adjacent bilayers are further apart than a cutoff in-plane distance ξ min (tuned to be 2.5 or 4 in units of the lattice spacing). This cutoff embodies the main effect of the strong Coulomb repulsion between Fe 3+ -Fe 3+ bipyramids. We consider small impurity bond concentrations n imp ≤ 0.04 ∼ J ⊥ /|J imp | per unit cell. Monte Carlo results include a disorder averaging over 16 configurations of the impurity bonds. A more complete description of the model, the method and detailed Monte The magnetic ions are depicted as pink spheres. They correspond to either Cu 2+ or Fe 3+ . They form stacked bilayers (regions enclosed by the rectangles) obtained by stacking a pair of x-y square lattices (corresponding to YBaCuFeO5 bilayers) along the z axis. The antiferromagnetic exchange coupling J within a layer is assumed to be independent of whether it couples Fe 3+ with Fe 3+ , Cu 2+ with Cu 2+ or Fe 3+ with Cu 2+ . The inter-layer coupling J ⊥ within a bilayer is ferromagnetic for Cu 2+ -Fe 3+ bipyramids (see text). Impurity bonds have a strong antiferromagnetic exchange, J imp , similar in magnitude to the ab initio value for Fe 3+ -Fe 3+ bipyramids. The impurity bonds are randomly distributed, up to the constraint that two impurity bonds in adjacent bilayers cannot be closer than a minimal distance ξ min . Carlo data including system-size dependence are given in the supplemental material [27].
For the clean, unfrustrated case (n imp = 0), we find a transition at T N 300 K = O(J ) from the hightemperature paramagnetic phase to a low-temperature collinear antiferromagnetic phase with the ordering wave vector ( 1 2 , 1 2 , 1 2 ). This is consistent with the antiferromagnetic ordering observed experimentally at high temperature in Ref. [12] and shown in Fig. 1(b). The calculated specific heat C at constant volume is shown in the first panel of Fig. 2(a). It shows a typical λ-peak at T N , characteristic of a continuous transition.
With a finite concentration of impurity bonds, n imp > 0, the peak in C broadens, while its position remains almost constant as long as n imp ≤ 0.04. However, the magnitude of the collinear order parameter, m AF , shown in the second panel of Fig. 2(a), is strongly suppressed below T = T spi 150 K for n imp = 0.02. This fact suggests the onset of spiral order. Simultaneously, the electric polarization (estimated as P 2 ) becomes nonzero. The associated susceptibility χ P exhibits a peak, which seems to diverge as the system size increases, as shown in the supplemental material [27]. In Fig. 2(b), we show the spin-structure factor S(q c ) as a function of the wave vector q c along the c axis (averaged over q a and q b ) and temperature, for four values of n imp . At n imp = 0.02, the propagation wave vector q c of the magnetic order decreases smoothly from π below T spi , suggesting a con-tinuous transition from the antiferromagnetic phase to a spiral-ordered phase, consistent with the experimental observations reported in Ref. [12]. A small residual m AF below T spi remains. This might be due to either finite-size effects or to a coexistence of the spiral and antiferromagnetic order [27]. Figure 2(c) shows T spi , estimated from the peak of χ P (see Supplemental Material), as a function of impurity concentration. At large impurity concentration, n imp 0.04, a direct transition from the paramagnet to an incommensurate spiral state with q spi < 1/2 is also compatible with the finite resolution of our data.
We carried out Monte Carlo simulations for two values of the minimal in-plane distance between impurity bonds ξ min = 2.5 and ξ min = 4 for each value of n imp . As n imp increases, T spi increases and almost reaches T N (as estimated from the peak in specific heat) at n imp = 0.04. Note that at n imp = 0.04, S(q c ) has a maximum at q c < π even at high temperature. This behavior does not rule out a direct transition from the paramagnetic to the spiral phase for larger (but not too large) values of n imp . Note that our values of T spi constitute lower bounds due to finite size effects in the Monte Carlo simulations. Indeed, the peak values for χ P are still moving to higher temperatures for the largest linear system sizes L that were used in the Monte Carlo simulations. Hence, the finite size effects on T spi are expected to be sizable for the values of L used in our Monte Carlo simulations. Finally, we discuss the short range repulsion among the impurity bonds, which we model by imposing ξ min . We find that when we neglect to include this effect, by setting ξ min = 0, the low temperature state does not support any spiral order over the range 0.02 ≤ n imp ≤ 0.04 that we studied. Instead, a "fan state" in which the intralayer ferromagnetic order oscillates is stabilized, as shown in panel (c) of Fig. 4; this state will be discussed in more detail elswehere [28]. Comparing the phase diagrams obtained for ξ min = 4 and ξ min = 2.5 in Fig. 2, we see that T spi decreases with decreasing ξ min . Indeed, the spiral order is favored if the presence of impurity bonds is suppressed in directions forming small angles with the c axis, where the coupling has antiferromagnetic sign.
In summary, our simulations reveal that in the range 0.02 ≤ n imp ≤ 0.04 the low temperature state is a spiral with a temperature-dependent wave vector provided the impurity bonds obey the condition that ξ min 2.5. Both the spiral ordering temperature T spi and the component 1/2 − q spi ·ĉ of the ordering wave vector depend on the concentration n imp of the impurity bonds and are proportional to the latter in the limit of small n imp .
Mechanism for spiral stabilization -Finally, we analyze the limit of low concentration of impurity bonds n imp and low temperature T , to provide insight into the mechanism for stabilization of the spiral. As discussed in detail in Ref. [28], impurity bonds induce cantings which are coupled through the spin stiffnes of the hosting spin lattice. For systems with continuous symmetry such effective coupling is dipolar and analogous to the Coulomb coupling that occurs for vortices in the two dimensional XY model. Moreover, the net value of dipoles couples linearly to the winding of the spins of the hosting lattice along the z axis. Therefore, distributions of impurity bonds having a ground state with a non-vanishing net dipole stabilize a spin spiral. At T ∆, the spins become coplanar XY spins with S ·n = 0. For n imp = 0, the antiferromagnetic ground state has all spins parallel in the xy-plane, e.g., parallel to the x axis. A single impurity bond with exchange of magnitude 1) renders the ground state two-fold degenerate, see Fig. 3. Far away from the impurity bond, the staggered magnetization parallel to the x axis is restored. Close to the impurity bond, the coplanar spins are canted away from the x axis with the two ground states differing in the local sense of rotation of spins, as given by P loc (r) = Sr ∧ Sr +ẑ ·n. At finite but small n imp > 0 the impurity bonds are well separated. Thus, the low energy configurations can be labelled by local Ising variables σr = sgn P loc (r) describing the sign of the local cantings, while the remaining degrees of freedom are well captured by spin waves. Integrating them out yields an effective Hamiltonian H eff for the Ising degrees of freedom for a system with N lattice sites, where we simplified the model slightly by setting |J ⊥ | = J ⊥ . We label the impurity bonds by the siter at their lower end and by L the ensemble of such sites for a given distribution of impurity bonds. Here, Q is the saddle-point value of the spiral wave vector and is proportional to the net value of the Ising spins, α = (|J imp | + |J |) |P imp loc | where P imp loc encodes the canting angle at an isolated impurity bond [28]. As shown in Ref. [28], at large distances the kernel Γ takes the (anti-) dipolar form Note that Γr ,r is ferromagnetic if ∆r 2 x + ∆r 2 y > √ 2|J | ∆r 2 z /|J ⊥ | and antiferromagnetic otherwise, as illustrated in Fig. 3(b). In particular, the ratio of the in-plane nearest-neighbor interaction to the inter-plane one scales like (J /J ⊥ ) 3/2 1. We thus expect ferromagnetic order in-plane. The inter-layer order might be either ferromagnetic or antiferromagnetic. Now, the Coulomb repulsion between Fe 3+ /Fe 3+ bipyramids amounts to a suppression of antiferromagnetically coupled pairs of nearest-neighbor impurities. A mean-field calculation shows that this constraint stabilizes ferromagnetic inter-layer order. In the numerical simulation, we retain the essential part of the short-distance repulsion by forbidding impurity bonds on adjacent bilayers to be closer than ξ min in plane. Increasing ξ min indeed enhances the tendency toward a ferromagnetic Ising ground state with a finite Q = n imp |α/J ⊥ | (as opposed to competing inter-layer antiferromagnetic order, which has vanishing Q).
The latter immediately translates into a spiral order of the spins S r , with the characteristic wave vector q spi = ( 1 2 , 1 2 , 1 2 − Q). The linear dependence of Q on n imp is in qualitative agreement with the n imp -dependence of the peak of the low T structure factor, as calculated with Monte Carlo simulations and presented in Fig. 2 (b).
The ground state of the effective Ising Hamiltonian (3) can be solved analytically when the set of impurity bonds forms a superlattice. Following Ref. [28], we assume that the impurity bonds form a Bravais lattice with the basis A, B, and C [29], such that n imp = a·(b∧c) A·(B∧C) 1. In Figs. 4(a) and 4(b) we compare the ground state of Eq. (3a) with the ground state obtained by Monte Carlo simulations of Eq. (1). For both calculations we take |J ⊥ | = |J ⊥ | = 4.1 meV to be the average of the actual |J ⊥ | and |J ⊥ | given in Fig. 1. In agreement with the magnetic structure obtained by the refined analysis of the elastic neutron scattering data, the obtained magnetic spiral has the feature that the rotation of magnetic moments mainly happens between neighboring ablayers coupled by impurity bonds (e.g., bilayers of the YBaCuFeO 5 structure). Finally, Fig. 4(d) shows the concentration dependence of the polarization P (see Eq. 2) calculated both numerically from Eq. (1) and using the effective Ising Hamiltonian (3) for various superlattices for which the ground state is a spiral. The values of polarization predicted by the effective model match well the values of polarization obtained numerically for the ground state of a finite system. Furthermore, the value of P increases proportionally to n imp for small values of P as predicted by the effective model.

Conclusions.
We have presented a mechanism in which impurity bonds, when sufficiently strong, induce a local frustration in an otherwise non-frustrated lattice of classical spins. The associated cantings become long-range  correlated at low temperature, resulting in a spiral magnetic order in which the sense of rotation of the spiral is spontaneously chosen. The spiral formation is in contrast to the spin-glass phase expected for white-noise correlated random magnetic exchange interactions [20][21][22][23]30], and is caused by orientational long-range correlations, which align the impurity bonds along one crystallographic direction of the lattice. The critical temperature below which the spiral order develops is controlled by nearest-neighbor exchange couplings, which can be sizable, so the mechanism is relevant for engineering high temperature multiferroics.
When the model is applied to YBaCuFeO 5 with realistic couplings, it predicts a multiferroic phase that is consistent with several distinct features observed experimentally. (i) It captures the transition to a commensurate antiferromagnet phase with wave vector q N ≡ ( 1 2 , 1 2 , 1 2 ) at the Néel temperature T N , followed by a transition to an incommensurate spiral phase with wave vector q spi ≡ ( 1 2 , 1 2 , 1 2 − Q) at the lower critical temperature T spi . (ii) It gives rise to a magnetic spiral phase, in which the rotation of the magnetic moments neighboring along c occurs mainly when they belong to the same bilayer. (iii) The observed temperature dependence of Q in YBaCuFeO 5 is well reproduced. (iv) Finally, the dependence of both T spi and q spi on the concentration n imp is in qualitative agreement with the dependence on the annealing conditions of the YBaCuFeO 5 samples: The faster the quench ( and thus the higher the expected defect concentration), the larger the measured T spi and q spi .
We close by mentioning other compounds with magnetic spiral order whose origin is not understood so far, but to which the mechanism presented in this paper might apply. For solid solutions of Cr 2 O 3 and Fe 2 O 3 [31] and for the doped hexaferrite Ba (1−x) Sr x Zn 2 Fe 12 O 22 [32], the wave vector of the spiral order is known to change smoothly from a commensurate value at high temperature to an incommensurate value at low temperature. In all cases, the maximal value of the concentration dependent T spi is high (T max spi = 148 K for solid solutions of Cr 2 O 3 and Fe 2 O 3 [31] and T max spi higher than 294 K for Ba (1−x) Sr x Zn 2 Fe 12 O 22 [33]). Moreover, the wave vector and the transition temperature to the spiral state are dependent on the dopant concentration. Since in both compounds cation substitution can introduce impurity bonds, they are likely candidates to realize "spiral order by disorder". Supplemental Materials: Multiferroic magnetic spirals induced by random magnetic exchanges The algorithm used for the Monte Carlo simulations is described in more details. First, we describe the model and the way in which thermodynamic averages are computed. Second, we compare the results obtained for the largest size of the system considered when the minimal distance ξ min is set to 4 and 2.5 in units of in-plane lattice vectors. Finally, the dependence of the thermodynamic averages on the system size is discussed.

MONTE CARLO (MC) SIMULATION
We model the magnetic properties of YBaCuFeO 5 (YBCF) as follows. We use the crystallographic unit cell labelled by the vector r ∈ R 3 and comprising two magnetic ions aligned along the c direction and belonging to the same bilayer. Therefore, in what follows, the position of each unit cell is defined by where m, n, l ∈ Z and a , b , and 2c is a basis for the crystallographic lattice vectors. We distinguish the two inequivalent magnetic sublattices with the index µ = 1, 2. As described in the main text, the classical spins S µ,r are unit vectors in R 3 and a single-ion easy-plane anisotropy with the characteristic energy scale ∆ > 0 is included. All together, the classical Hamiltonian [Eq. (1) in the main text] is There are three characteristic energy scales. The antiferromagnetic in-plane exchange coupling J ≤ 0, the intra-unit cell (intra-layer along the c direction) ferromagnetic exchange coupling J ⊥ ≥ 0, and the inter-unit cell (inter-layer) antiferromagnetic exchange coupling J ⊥ ≤ 0 (also along the c direction). They enter Eq. (S1b) according to the rules, J 1,2 r,r = J ⊥ δ r,r + J ⊥ δ r,r +2c , (S1c) J 2,1 r,r = J ⊥ δ r,r + J ⊥ δ r,r −2c , (S1d) The subscript L for the Hamiltonian (S1b) indicates its dependence on the choice made for the set L of impurity bonds. This dependence arises from the term We employ the exchange Monte Carlo [34] and overrelaxation [35] method using 200 temperature (T ) points in the range of 10 ≤ k B T ≤ 30 (meV −1 ) for L = 28. Replicas are exchanged between two neighbouring temperatures. The distribution of the temperatures is optimized so that the exchange rate is independent of T . Thermodynamic observables are measured using the reweighting method [36]. Periodic boundary conditions in the ab plane and open boundary conditions along the c direction are used in order to deal with a magnetic spiral order with an incommensurate wave vector.
Monte Carlo simulations are performed for several sets L chosen randomly. As described in the main text, the randomly chosen impurity bonds satisfy the constraint that there be a minimal lateral distance ξ min between impurity bonds separated by 2c along the c axis of the pseudocubic lattice. We denote with · · · MC the Monte Carlo average and with · · · L the average over randomly chosen sets L. Denoting with Eq. (S1a) the position of the unit cell, we define the collinear antiferromagnetic order parameter m AF as the one that would occur in the absence of impurity bonds The specific heat is defined by The order parameter for the spiral state is defined to be where generalizes Eq. (2) in the main text. The structure factor at We start by discussing the ξ min dependence of the results. Figure S1(a) shows MC data obtained for L = 28 at n imp = 0, 0.02, 0.03, 0.04 for ξ min = 4. Figure S1(b) shows the results for the case ξ min = 2.5. For both values of ξ min , the low-temperature values of the collinear order parameter m AF (upper panel) are suppressed as n imp is increased. Moreover, for n imp = 0, m AF shows the re-entrant behavior described in the main text as the temperature is decreased. Interestingly, for the case of ξ min = 2.5 the collinear antiferromagnetic order parameter m AF seems to not vanish completely at low temperatures for the lowest impurity concentration. This might either arise from finite size effects (too few windings of the spiral within the finite lattice) or indicate a possible coexistence of the magnetic spiral state with a collinear spin structure. The middle panels of Figs. S1(a) and S1(b) show the temperature dependence of the specific heat C. The peak position in C coincides well with the onset of the collinear order parameter m AF and is approximately unaffected by n imp and ξ min . The polarization P (bottom panels) becomes nonvanishing below T spi . For T < T spi , P increases continuously with decreasing T to some saturation value at T = 0. Hence, T spi is interpreted as the onset of the spiral phase through a continuous phase transition provided n imp > 0. The temperature below which the magnetic spiral state appears is better identifiable by the positions of the peaks of the dielectric susceptibility χ P , defined as and shown in the insets of the lower panels in Fig. S1. Figure S2 shows the dependence on q c in Eq. (S5) of the spin-structure factor (S6) computed for eight different configurations of impurities. The size of the system and the minimal distance are set to L = 28 and ξ min = 4,  respectively. One clearly sees a peak in the dependence of S(q c ) on q c at a value of q c that deviates from 0.5 as n imp is increased. This deviation approximately increases linearly with the concentration of impurity bonds n imp , in agreement with the predictions of the low-temperature effective model.
To study the approach to the thermodynamic limit, we investigate the system-size dependence of the MC results for n imp = 0.02. As shown in Fig. S3, the dependence of m AF on the system size is not monotonic around 120 K: m AF has a small but non-vanishing residual value even at L = 28. Therefore, as the system size is not large enough to obtain convergent results for m AF , we cannot exclude the possibility of the coexistence of the spin collinear order and the spin spiral order at low temperatures. On the other hand, P converges with increasing L rapidly for the lowest temperature considered to a nonvanishing and finite value. For example, at T 23 K, P does not change appreciably when increasing L from L = 12 to L = 28. This observation establishes the existence of the spin-spiral state at low T . However, the peak position in χ P drifts toward high temperatures as the system size is increased. This indicates that our calculations might actually underestimate the value of T spi . Figure S4 shows the system-size dependence of T N (peak position in C) and T spi (peak position in χ P ). This clearly shows that T spi is still monotonically increasing for L increasing from 12 to 28, while T N is approximately constant for L increasing from 12 to 28.
Although our MC results for given large L × L × L strongly suggest the existence of a finite T spi , extrapolating T spi in the thermodynamic limit L → ∞ with n imp held fixed is difficult. Furthermore, the transition may be sensitive to the geometry of the system since the effective interactions between impurity bonds are long ranged. Thus, we systematically study the geometry dependence of the spin-spiral transition. Figure S5 compares T spi obtained for different geometries for n imp =0.02 and ξ min =4. We adopt tetragonal shapes by which the dimensions L × L × L are substituded by the dimensions L a × L b × L c and we vary the ratio L a /L c = L b /L c from 1 to 8. As the ratio L a /L c becomes larger, T spi tends to increase for the same numbers of spins N . This trend is consistent with the dipole effective interaction, which is ferromagnetic in transverse directions a and b . All the series of T spi for different values of the L a /L c ratio extrapolate to T spi 200 -250 K in the thermodynamic limit.
To confirm this result, we present the MC simulation for the Binder cumulant g P of the the polarization P , where Π is defined by Eq. (S4b). The Binder cumulant is first computed for each configuration and then averaged over different configurations of impurity bonds. We show the temperature dependence of g P for various values of the ratio L a /L c in Fig. S6. One clearly sees an intersection of the data at T spi 225 K, which is lower than T c . The