Microscopic theory of magnetic detwinning in iron-based superconductors with large-spin rare earths

Detwinning of magnetic (nematic) domains in Fe-based superconductors has so far only been obtained through mechanical straining, which considerably perturbs the ground state of these materials. The recently discovered non-mechanical detwinning in EuFe$_{2}$As$_{2}$ by ultra-low magnetic fields offers an entirely different, non-perturbing way to achieve the same goal. However, this way seemed risky due to the lack of a microscopic understanding of the magnetically driven detwinning. Specifically, the following issues remained unexplained: (i) ultra-low value of the first detwinning field of $\sim$0.1$\,$T, two orders of magnitude below that of BaFe$_{2}$As$_{2}$, (ii) reversal of the preferential domain orientation at $\sim$1$\,$T, and restoration of the low-field orientation above 10$-$15$\,$T. In this paper, we present, using published as well as newly measured data, a full theory that quantitatively explains all the observations. The key ingredient of this theory is a biquadratic coupling between Fe and Eu spins, analogous to the Fe-Fe biquadratic coupling that drives the nematic transition in this family of materials.

Introduction One of the most admirable experimental feats in studies of Fe-based superconductors (FeBS) is the mechanical detwinning of the low-temperature phases of the parent compounds in the 122 families [1,2]. This allowed impressive insight into the physics of spin-driven nematicity, a phenomenon that arguably rivals the superconductivity itself in these materials. In this connection, one of the most intriguing and unexpected findings was that this nematic physics is ensured by a sizable biquadratic magnetic interaction, something unheard of in localized magnetic moment systems, and never investigated in itinerant magnetic metals. This phenomenon was first discovered computationally [3] and later shown to provide the only physically meaningful description of spin dynamics in FeBS [4]. There is growing evidence that it is not limited to FeBS, but occurs also in other itinerant systems [5].
Mechanical straining is not the only way to detwin FeBS. It was shown that a static magnetic field of ∼ 15 T leads to partial detwinning [6] and pulsed fields of ∼ 30 T to nearly complete detwinning [7]. Later we will analyze these facts in more detail, but at the moment we emphasize that these are relatively large fields, even though the in-plane magnetic anisotropy energy of the FeAs planes was experimentally shown to be of the order of 0.5 meV and, therefore, sizable compared to e.g. elemental Fe, where it is only a few µeV. Against this background, it came as a complete surprise when it was discovered [8,9] that substituting Ba with Eu lowers the field needed for full detwinning by two orders of magnitude. In principle, this magnetic detwinning allows for a virtually nonintrusive (the energy scale associated with this field is less than 20 mK) investigation of the physics of the nematic state. However, this seemingly exciting opportunity was met with limited enthusiasm for the simple reason that no plausible microscopic explanation could be found for the detwinning itself, given the minuscule amplitude of the required field. Even more striking was the discovery that by increasing the magnetic field gradually one can switch the sign of detwinning twice: initially, twin domains orient in such a way that the Fe-Fe ferromagnetic bonds along the crystallographic b are parallel to the applied field (we call this the b-twin, see Fig. 1a). This process is essentially complete at H 0 ∼ 0.5 T. Then, at H 1 ∼ 1 T domains spontaneously rotate in-plane by 90 o , and at H 2 ∼ 10 T start to turn back to the detwinned state that was initially generated at H 0.5 T, see also Fig. 1(b-h). With such a complex phase diagram, and no theoretical understanding of the underlying phenomena, it is indeed worrisome to embark on systematic studies of nematicity with the risk that unknown magnetic physics may affect the findings. The goal of this paper is to remedy this situation and present a full and quantitative theory explaining all the above observations. It appears that magnetically induced detwinning is intimately related to the nature of nematicity itself, namely, it is also driven by a sizable biquadratic interaction, which, in turn, is the consequence of the Janusian itinerant-localized nature of the FeBS.
Formalism We start with the (simpler) case of BaFe 2 As 2 . The minimal approximation is the nearestneighbor (n.n) Heisenberg model with single site anisotropy: where i, j label Fe sites, f i is the unit vector directed along the Fe magnetic moment at site i,M is its absolute value, and the summation is over all inequivalent n.n bonds [10].
H is the external field in the corresponding energy units.
Here and below we use tildes over symbols for Fe-related parameters.
It is known experimentally that the Fe moments lie within the ab plane and are oriented along the longer a-axis (i.e.,D > 0). Neutron scattering provides esti- mates for the parameters asJ ∼ 30 meV,D ≈ 0.25 meV [11]. Minimizing Eq. (1) with respect to f, for H||b (no linear susceptibility appears for H||a) we observe that H generates a canting of the Fe spins away from the a axis by an angleα = sin −1 (M H/4J), which corresponds to an energy gain of E = (M H) 2 /4J per formula unit (f.u). There are two ways how the system may take advantage of this energy gain, even if the field is along a. First, all Fe spins may rearrange (abruptly) from being aligned along a to being predominantly aligned along b. This process is called spin-flop and occurs whenJ >D. UsingM ≈ 1 µ B and the above parameters, we esti-mateH flop = 2 2DJ/M ≈ 11 meV ≈ 135 T. Thus, the spin-flop never occurs at typically lab-accessible fields. The other way is to switch an entire "a-twin" (i.e., a domain with antiferromagnetic bonds along the crystallographic a in field direction) to a "b-twin"(in other words, to rotate the crystal structure by 90 o around c, keeping the moments aligned along the magnetic easy a axis. This process is associated with an unknown energy barrier ∆. Given that in the field of ∼ 15 T the Stanford group has observed partial [6], and in ∼ 30 T nearly full detwinning [7], we deduce ∆ ∼ E = (M H) 2 /4J ≈ 0.025 meV/f.u.
One can verify these deductions against mechanical detwinning. The latter is an indirect process wherein it is difficult to access the microscopic strain and stress.
Reported values for the latter differ from σ = 6-20 MPa [12,13]. Assuming σ = 10 MPa and taking the elastic modulus to be 10 GPa [14,15], we derive a strain of ε ≈ 0.1%. Finally, using the calculated dependence of the total energy on the microscopic strain [16], we find that this strain corresponds to ∆ ≈ 0.01 meV/f.u, which is in agreement with our magnetic estimate.
Yet another estimate can be obtained by considering the stress σ on the unit cell during mechanical detwinning and calculate the energy associated with the displacement δ from a σ to b σ. Using the reported lattice parameters of EuFe 2 As 2 [17], we arrive at δ = 1.6×10 −12 m. Assuming σ = 6 MPa, this leads also to an energy of ∆ ≈ 0.01 meV/f.u.
Given that the detwinning energy, according to this theory, depends quadratically on the magnetic moment, and inversely on the exchange constant, one may naively assume that the same mechanism will be operative in EuFe 2 As 2 , given that the ordering temperature in the Eu sublattice is much smaller, T N = 19 K, and the moment much larger, M = 7 µ B , than for the Fe sublattice. Such a phenomenology was adapted in Ref. [9] in order to parametrize the observed effect. However, it is easy to see that it is microscopically untenable. Indeed, while it is possible and reasonable to write down interactions inside the Eu sublattice in Heisenberg form, where i, j label Eu sites, i + the opposing Eu site in the next layer, e i the unit vectors directed along the Eu magnetic moment M at site i, and the ferromagnetic J < 0 and antiferromagnetic J ⊥ > 0 constants determine the in-plane and interplanar ordering, respectively, it is not possible to describe the interaction between the Fe and Eu subsystems in the same manner, for the simple reason that the Heisenberg exchange field induced by the Fe planes on the Eu sites is zero by symmetry. In fact, any bilinear coupling between the Fe stripes and in-plane Eu spins is zero by symmetry, including Heisenberg, Dzyaloshinky-Moria and dipole interaction (in the last case the field induced by the Fe planes on the Eu sites is non-zero, but is directed strictly along z, see supplement). Note that we did not include any single-site anisotropy in Eq. (2), because Eu adopts a valence state of +2 in this system. Due to the closed f −shell, with 7 electrons in the spin-majority channel, Eu 2+ has zero angular momentum and negligible magnetocrystalline anisotropy. This is confirmed by first principles calculations, presented in the supplement to this article. Without an interaction between Fe and Eu, there is no physical mechanism by which the Eu spin dynamics may affect the detwinning.
Another intriguing problem, possibly related to this one, is the fact that even the basic magnetic properties of EuFe 2 As 2 cannot be explained within a simple Heisenberg model. Indeed, the magnetic susceptibility of EuFe 2 As 2 above T N is dominated by the Eu spins and well described by a nearly isotropic Curie-Weiss law, in accordance with the previous paragraph. Eq. (2) suggests that k B T CW = 3 7 (4J − 2J ⊥ ) (the quantum prefactor 3M/(M + 2) would have been 1 if M were 1 µ B , but for M = 7 µ B it becomes 3/7) The effective moment has been determined to be M = 7 µ B with T CW ≈ −20 K [18]. Thus the Néel temperature appears to be equal to the mean-field transition temperature of the individual Eu planes. In other words, each Eu plane orders magnetically at the mean-field temperature, not at all suppressed by fluctuations, and immediately at the transition the antiferromagnetic stacking of the individual planes along c is acquired.
At this point, it is instructive to look at first principles calculations and what they tell us about J ⊥ and J . The former appears to be very small and decreases with the value of the Hubbard U used on the Eu f orbitals. This is not surprising, because it is set by the competition between superexchange, proportional to t 2 ⊥ /U, and Schrieffer-Wolff driven [19] double exchange, proportional to (t 2 Eu-Fe /U ) 2 N (0), with N (0) the density of states, and t ⊥ and t Eu-Fe the effective Eu-Eu and Eu-Fe hopping across the planes. For U − J = 0 we get J ⊥ = 0.14 meV (1.6 K) and J = −2 meV (23 K), while for U − J = 5 eV we find J ⊥ = 0.26 meV (3 K) and J =-0.8 meV (9 K). The LDA+U results, which we believe are closer to reality, correspond to T CW ∼ −13 K, and T N ≈ 4.27J /[3.12 + log(J /J ⊥ )] = 9 K for the Heisenberg model [20]. The ratio T CW /T N is close to 1.4 rather than the experimental 1.13, and one can see that in order to reduce it to 1.13 one needs to increase the J ⊥ /J ratio to ∼ 0.5, a rather unrealistic 3-dimensionality.
The situation could be remedied if one were to assume a finite single-site anisotropy for Eu of D ≈ 1 K, because in such a case one has to replace J ⊥ with 21,22], which leads to a sufficient increase of the Néel temperature and consequently to T CW /T N = 1.13. But, as we have argued above, Eu 2+ has no single-site anisotropy. One of the results of our paper, however, is that a biquadratic coupling K between the Eu and Fe subsystems is operative in EuFe 2 As 2 , which plays the key role in its magnetic detwinning. This interaction acts as an effective anisotropy for the Eu subsystem, which tries to imprint the orientation of the Fe spins on the Eu spins. The corresponding number (8K) that we extract from the experiment amounts to D eff ∼ 0.6 K. Substituting this into the formulas above we get T CW /T N ≈ 1.18, which is within the error margin of the experimental ratio. Of course, given the model character of these calculation, a discrepancy of ≈ 25% is not too alarming. However, it is noteworthy that after adding the experimentally determined biquadratic term the discrepancy virtually disappears completely.
Note that, unlike correlated insulators, FeBS have a considerable biquadratic interaction within the Fe subsystem, which plays a crucial role in their nematic behavior [23]. Whether this is a universal property of correlated magnetic metals, which simply had not been given proper attention before, or is unique for FeBS is unknown. With this in mind, we have combined the Hamiltonians from Eqs. (1) and (2) in the following way: where the Greek subscripts label the Fe sites, Latin the Eu sites, and the last summation runs over all n.n Fe-Eu pairs. We have then estimated the biquadratic Eu-Fe coupling term K from LDA+U calculations (see SM), and obtained an estimate of K ∼ 0.4 meV. Given the uncertainty in the calculations, this should be taken as evidence that K is not negligible; later we will determine the actual amplitude of K directly from the experiment.
In the Supplementary Information we present detailed derivations of the solutions of Eq. (3) for both possible orientations of the external field with respect to the crystallographic axes, and for all relevant field regimes. The discussion below omits less relevant parts of the full the-ory, concentrating on rationalizing the actual experimental observations. Domain energetics We start our discussion in the low field regime, that is 0 < H H 1 ∼1 T, before considering higher fields. In this regime the Fe single-site anisotropy dominates and Fe spins are always oriented along the crystallographic a axis, and so are, initially, Eu spins. We will distinguish two cases: first, when H is applied perpendicular to the initial orientation of Eu spins (b-twin) or, second, parallel to those (a-twin). In this regime the former is always lower in energy, since the latter has formally zero spin susceptibility. We will characterize the orientation of the Eu and Fe spins by their respective angles ϕ andφ with respect to the external field. Figure 1(a) is drawn for a b-twin domain, i.e., a ⊥ H. Equation (3) can then be rewritten in terms of these angles (per one formula unit) as where α's are the Fe angles measured from the magnetic easy a-axis.
For the b-twin domains and low fields, i.e.φ 2 =φ 1 = π/2, ϕ 1 = ϕ 2 = ϕ,α 2 =α 1 = 0 this yields where we have expressed everything in terms of p = cos ϕ The equilibrium tilting angle and energy are given by The biquadratic term in the definition of E 0 is always trying to minimize the angle between the Fe spins and Eu spins. Thus, in simplified terms (see supplement for details), if π/4 < ϕ < π/2, Fe spins prefer to orient perpendicular to the field (φ = π/2), and the b-twin domain is always lower in energy. When ϕ becomes smaller than π/4, it becomes more favorable to orient Fe spins parallel to the field (φ = 0), and this is the first critical field H 1 at which the first domain reorientation from b-twins to a-twins takes place (note that in this field regime onlỹ ϕ = π/2, or, in case of the a-twin, 0 is allowed, any intermediate value ofφ is severely punished by the Fe-Fe exchange and single-site anisotropy).
After this reorientation has occurred, the total energy is expressed as and consequently p min Thus the first reorientation field H 1 is defined by E min At the field H sat a = 4(J ⊥ −4K)/M the angle ϕ becomes 0, that is to say, Eu spins are perfectly aligned with the field. Further increase of the field does not change the total energy (aside from the Zeeman term −M H), because from that point on the differential spin susceptibility of the Eu subsystem becomes zero. However, while theoretically important, this field does not manifest itself as a change of regime in domain dynamics.
H sat a is too small to incur any Fe spin dynamics, but, in principle, with further increase of H one needs to include the spin susceptibility of the Fe subsystem. The latter is zero as long as Fe spins are parallel to Eu spins, satisfying the biquadratic coupling. Yet, in a sufficiently strong field, i.e. at H H 2 a potential energy gain from allowing Fe spins to screen the field outweighs the loss of the biquadratic interaction. Mathematically, the latter, in this regime H > H sat a , is reduced to an effective single-site anisotropy for the Fe subsystem, subtracting from the actual anisotropy, and the transition in question is described by the same formulas as a typical spin-flop transition. We can find the corresponding field value in the same way as one derives the spin-flop field in textbooks: we need the energy gain from the Fe spin-flop to overcome the energy loss due to noncollinearity (the loss occurs both due to the Fe-Fe exchangeJ and because the Fe site-anisotropy is much larger than the biquadratic coupling).
In this case, the total energy (since now ϕ = 0,α = where we have now usedp = cosφ and E 0 . Minimizing with respect top yields: with the energy gain compared to Eq. (8) with p = 1 being This expression changes sign at . This is the second critical field at which the reorientation back to the b-twin domains is initiated.
The domain dynamics, with the initial detwinning to b-twins at H 0 , first reorientation to a-twins at H 1 , and second reorientation back to b-twins at H 2 is illustrated in Fig. 1(b-h) and also in a supplemental movie (url). There, we depict how Eu and Fe spins (and the structural axes follow the latter) rotate in an external field. One can see that, despite the physical simplicity, the actual dynamics are rather complicated.
Determination of the Coupling Constants While the H 0,1,2 defined above directly manifest themselves in the experiment, the actual expression for the energy difference between the two types of domains (which is needed to describe domain dynamics, as opposed to the thermodynamic equilibrium) is a complicated piecewise function of the field. The full derivation can be found in the Supplement, where explicit expressions for all critical fields are obtained. In the relevant field range for the experiments discussed here this function is given by In the second case dE > 0 changes to dE < 0 at The coupling constants J ⊥ and K can be determined from experiment. Magnetization, magnetostriction, neutron and magneto-transport measurements, for instance, all can be used to estimate the domain population ratio. In the following we will determine the coupling constants using new magnetization data [9], and then use them to calculate this ratio as a function of field in order to compare it to the experiment.
When measured along the tetragonal [110] T direction, a roughly linear magnetization M (H) was reported in Ref. [9], interrupted only by two pronounced jumps around 0.1 T and 0.6 T. Saturation sets in above 0.9 T with M ≈ 5 µ B . Due to the uncharacteristically small saturated moment, which does not agree well with previously published data [18] and the theoretically expected value of M = 7 µ B , we have re-measured the magnetization of a EuFe 2 As 2 sample from Zapf et al. [9]. The results are shown in Fig. 2 for decreasing magnetic field. The overall behavior is similar, but we found a saturation magnetization of 6.82 µ B , in good agreement with the theoretical expectations and Ref. [18].
The step-like increase of the magnetization around 0.65 µ B T was associated with a spin-flip transition of the Eu 2+ moments in Ref. [9]. However, magnetostriction, magneto-transport and unpublished neutron diffraction data indicate that the feature is associated with the reorientation of domains, rather than an intrinsic spin-flip of the a-twin domain and must, therefore, be identified with H 1 . Furthermore, the Eu saturation field of the btwin domain can be extrapolated to H sat b ≈ 1.14 T from the slope of the low-field region between 0.2 T and 0.5 T (Fig. 2); note that this field needs to be extrapolated, and cannot be measured directly because at H 1 T virtually all domains are a-twins (Fig. 1g). The constants can be extracted from these two fields using the following set of equations: which yields J ⊥ = 0.093(11) meV and K = 0.0049 (18) meV for this particular sample. The determination of the coupling constants via a different set of equations is discussed in the supplement. The respective results are summarized there in Tab. (SII). We have cross-checked the results by determining the parameters from various samples grown with different methods and investigated with various measurements techniques like magnetostriction, magneto-transport and neutron diffraction data. The results appear to be consistent between measurements, but we found a noticeable sample dependence, which also seems to be related to the synthesis methods of the single crystals, see also Tab. (SII). The averaged values are J ⊥ = 0.121(24) meV, K = 0.0072 (22) meV and J ⊥ /K ≈ 18(4). However, since each Eu atom is surrounded by 8 Fe atoms, 8K = 0.057 meV is the more representative quantity to gauge the biquadratic coupling strength in the system.
Energy barrier and domain dynamics The domain dynamics are driven by the energy difference dE between the domains. In Fig. 3 we show dE on a semilog plot as a function of the reduced variable M H/J ⊥ for a representative ratio of K/J ⊥ = 0.05. Positive values correspond to the b-twin domain being the ground state, while for negative ones the a-twin is the ground state. The phase diagram in thermodynamic equilibrium over a large parameter space in reduced coordinates {K/J ⊥ , M H/J ⊥ } is shown in Fig. 4. The remainder of the formulas used in constructing this phase diagram can be found in the supplement.
Until now we have assumed that the reorientation of twin domains has no energy cost. In reality, however, there is an unknown energy barrier ∆ associated with the reorientation, i.e. the energy difference dE = E a-twin − E b-twin between the two twin variants needs to exceed a certain threshold before reorientation occurs. In the following, we will assume ∆ for various domain walls to be log-normal distributed, i.e. log ∆ is normal [24]. This is a typical distribution e.g. of grain sizes in polycrystalline matter [25]. The cumulative distribution function F X (x) of a positive log-normal distributed variable , with the location and the scaling parameters µ and σ. Due to the fact that dE changes sign as a function of applied field, the log-normal distribution in our case leads to the following fit function to the (a-twin) domain population: with n 0 the fraction of a-twins at H = 0 and dn = n sat − n 0 the difference between the saturated a-twins n sat and n 0 . Together with Eq. (12) this function fully describes the domains population. Prior to performing the fit, the domain population needs to be extracted from the available measurements. The domain population n(H) can be determined in a variety of ways. Arguably the most exact data can be extracted from field-dependent neutron diffraction measurements on a free standing sample, which will be published soon [26]. In the following we will use magnetostriction (MS) data by Zapf et al. [9], which agree with the preliminary data of Ref. [26]. The assuming full detwinning at the observed minimum around H 0 , with b||H, i.e. n(H 0 ) = 0. This assumption is justified, as preliminary neutron data [26] indicate a domain distribution with n(H 0 ) = 0.06. Furthermore, the significant pressure of the dilatometer in field direction, 1.35 MPa [9] aids the n = 0 alignment at small magnetic fields (but hinders it at larger fields, when n → 1). The extracted domain population is shown in Fig. 5. The solid red line represents the fit to the data (blue symbols). The energy barrier for this particular sample and measurement technique was determined to ∆ = 0.01 meV. Among the other investigated samples ∆ ranges roughly between 10 −3 meV and 10 −2 meV, in agreement with the estimates presented in the introduction.

Consistency check
Using the averaged constants from the previous paragraph we find H 1 = 0.85 T and H 2 = 35 T, in very good agreement with experiment. Utilizing the knowledge about the energy barrier, we can also calculate the initial detwinning field from the condition dE = ∆, through For the determined ∆-range between 0.001 meV and 0.01 meV this yields 0.09 T to 0.28 T which is also in excellent agreement with the experimental evidence.
Going even further, the presented theory allows us to calculate our new magnetization data, by weighting the Eu magnetization of the twin sublattices with the domain population n(H) (dashed line in Fig. 2), that we calculated using ∆ = 0.004 meV, n i = 0.22 and dn = 0.71.  Fig. 3(b). The a-twin domain population is reduced, due to the pressure of the dilatometer, which favors the b-twin domain in this orientation.
The total magnetization is given by where M a and M b are the magnetization of the sublattices given by M cos ϕ, i.e.
{H sat a , ∞}, The result is depicted by the solid red line in Fig. 2. Summary We present a microscopic, physically meaningful and quantitative description of the observed magnetic detwinning effect in FeBS, with all its complexity. In particular, the following mysteries have been resolved: (i) strong detwinning in minuscule fields despite absence of spin-orbit coupling effects in Eu 2+ ions; (ii) coupling of Eu spin orientation to the Fe sublattice, despite any bilinear interaction canceling out by symmetry, and (iii) double reversal of the preferential domain orientation with the increase of the external field. We show that all these issues find a natural explanation if a biquadratic Eu-Fe coupling is included in the model. We also show that such a term does actually appear in first-principles calculation, with an amplitude even stronger than needed to explain the experimental data. Furthermore, we were able to describe quantitatively not only the thermodynamic phase diagram, but even the dynamics of detwinning (as deduced from our new magnetization data), assuming that the detwinning energy barriers are distributed according the log-normal law (quite typical in crystal morphology).
Although we focus on the Eu-based 122-system, our findings should be of great interest for other large-spin rare earths as well, such as the Gd-based 1111-compound, which, like EuFe 2 As 2 , also features large S=7/2 spinonly moments, or the recently discovered Eu-based 1144systems. Furthermore, the theory also captures the physics of the high-field (H 15 T) detwinning, which occurs even in non Eu-based iron pnictides. The microscopic understanding of the phenomenon of magnetic detwinning in ultra-low magnetic fields, as deduced above, opens new avenues for the experimental investigation of spin-driven nematicity in Fe-based superconductors. Acknowledgments We are grateful to Shibabrata Nandi and Yinguo Xiao for giving us access to their neutron diffraction data, to Makariy Tanatar for sharing with us his data on mechanical detwinning and to Ian Fisher for discussing with us the concept of this work. We also like to convey our gratitude to Shuai Jiang, Christian Stingl and Jeevan H.S. for permitting us access to their samples and magnetostriction data and Sina Zapf and Martin Dressel for collaboration in the early stage of this project (Ref. [9]), and to James Glasbrenner for verifying our biquadratic calculations using an all-electron method. J.M. thanks Patrick Seiler for support and discussions. J.M. and P.G. are supported by DFG through SPP 1458. I.I.M. is supported by ONR through the NRL basic research program.

A. Computations
In order to estimate computationally the key parameters of the theory (keeping in mind that the final values would be extracted from the experiment), we have performed calculations using the plane wave code VASP [1] with the setup as described in Ref. [2], and using the Generalized Gradient Approximation (GGA) to the Density Functional Theory (DFT). Our first goal was to verify the conjecture based on the general considerations presented in the main text, that the single-site anisotropy (SSA) for Eu can be neglected. This is not a trivial task: rotating both Eu and Fe spins together incurs a much larger Fe SSA, while rotating only Eu entangles Eu SSA with the Fe-Fe biquadratic coupling (and may be unreliable for different reasons, as discussed below). We can assess the SSA anisotropy of Eu indirectly though, bypassing the problem of Eu-Fe coupling entirely. We can replace Fe 2+ with the isovalent, but nonmagnetic Zn 2+ . The EuAs trilayers then remain intact, and if there is any coupling of Eu magnetic energy with the orthorhombic distortion of the Eu sublattice, it should manifest itself in EuZn 2 As 2 as well. Everything else is interaction between Eu spins and Fe spins, which can be described as a combination of a Heisenberg (bilinear) and a biquadratic coupling. We have performed this calculation using both the VASP code and the all-electron Linear Augmented Plane Waves WIEN2k code [3], and obtain the SSA of 0± 0.007 meV, with the error window an order of magnitude smaller than required to explain the experiment. We can get an independent estimate of the Eu SSA by going back to our compound and fixing the Fe spins along the c axis, while rotating Eu spins in the xy plane. This way, the biquadratic coupling between Eu and Fe remains always zero. In these calculations we found the total energy to be independent of the Eu orientation within a similar error bar, thus confirming that the Eu closed shell does not allow for any SSA.
We also note that possible anisotropic exchange interactions of the form W (M x ij M x i+1,j − M y ij M y i,j+1 ) contribute to the calculated energy differences in both these procedures at the same footing as SSA (since the Eu spins remain collinear), so the null results of the calculations also excludes these type of effects. Finally, the Dzyaloshinskii-Moria interaction in this system is also ruled out, because the midpoint between two Eusites is an inversion center. These two interaction exhaust relevent anisotropic exchange term allowed by the square-planar symmetry.
We have also tried to calculate the biquadratic coupling K from first principles. To this end, we have performed two sets of calculations: one with Fe spins along the easy axis a, and Eu spins either parallel or perpendicular to the former, and then taking the energy difference, or with Fe spins along b, and Eu spins either parallel or perpendicular to the former, again taking the energy difference. Not unexpectedly, the result depends on the Hubbard U, but is consistent between the two sets for the same U . In particular, for U − J = 7 eV, we found K ≈ 0.4 meV. We have also performed calculations gradually rotating the Eu moments with respect to the Fe ones, using the constrained-magnetic-moment formalism, incorporated in VASP. We have obtained the total energy roughly proportional to cos 2 φ, where φ is the angle between the Fe and Eu spins, with the coefficient consistent with the above estimate. Compared to the values extracted from the experiment, these num-bers are too large, partially because the GGA underestimates the wave function localization, and the biquadratic term is small and very sensitive to localization. However, the discrepancy is still too large. It should be noted that while collinear LDA+U+SO calculations have a solid history by now, and have been tested for many materials, including those containing rare earths, the same cannot be said about noncollinear calculations. Proper account of the nonspherical part of the double-counting term in LDA+U is very challenging [4] . Because of that, VASP developers recommend using a spherically averaged potential, replacing U with U − J and setting J to zero. This is what we did, and, in fact, we are not aware of any publicly available DFT package that would allow for noncollinear magnetic calculation and properly included nonspherical LDA+U. With this in mind, we have tested the stability of the above-mentioned coefficient in front of cos 2 φ with respect to U . We found that it depends unphysically strongly on U . This means that while calculations can be taken as an indication that the experimental, or even larger biquadratic coupling is consistent with the band structure calculations, its amplitude cannot be reliably estimated computationally.
Estimating J ⊥ does not require noncollinear calculations and is more reliable. we compared the total energies of ferromagnetic Eu planes, stacked ferromagnetically or antiferromagnetically along c. Finally, for estimating J || we compared ferromagnetic and checkerboard-antiferromagnetic in-plane arrangements of Eu ions, always allowing stripe antiferromagnetism for Eu. All these patterns are illustrated in Fig. S1.

The role of bilinear coupling between Fe and Eu
The fact that the Heisenberg coupling induces zero exchange field at the Eu sites (which project onto the centers of Fe 4 plaquettes) from the stripe-ordered Fe planes is quite obvious. One can show that not only Heisenberg exchange, but any symmetry-allowed bilinear spin-spin coupling generates at the Eu sites, at best, an exchange field perpendicular to the Fe planes. In fact, since the symmetry positions of the As atoms with respect to the Fe planes are exactly the same as the Eu atoms, the following derivation exactly follows the well known theory of the As hyperfine field in BaFe 2 As 2 in the NMR community (see, for instance Ref. [5]) .
Let us assume an arbitrary bilinear coupling between Fe and Eu, of the form, for a single bond with e and f the Eu and Fe spins, respectively and summation over repeated Cartesian indices. The structure of T depends on the direction of the bond. The Heisenberg exchange corresponds to T = c · δ αβ , dipole interaction to c · (3d α d β − d 2 δ αβ ), the DM interaction to c · ε αβγ d γ (equivalent to d · e × f ), where c is a constant and d the Eu-Fe bond vector. In general there are symmetry constraints on allowed elements of the T matrix, but our derivation below will not be using them explicitly.
We consider a plaquette of four Fe atoms at {±1, ±1, 0} (site labels starting from {1,1,0} counting counterclockwise), with their spins along x as {1, 0, 0}, {−1, 0, 0}, {−1, 0, 0}, {1, 0, 0}, and Eu at {0, 0, h}. The exchange field generated by the Fe 1 at the Eu site is In the following we show how the exchange fields of the other sites can be expressed in terms of the Fe 1 site with spin direction either along {1,0,0} or along {0,1,0}. The Fe 2 site e.g. is obtained by rotating a {0,1,0} spin at the Fe 1 site (corresponding to h α = T (1) αy ) around the z axis by π/2. This rotation is given by the matrix After the rotation we get a spin along {−1, 0, 0} = f (2) at site 2, generating the field βy . By symmetry, this must be the same as the field generated by f (2) through the matrix T (2) , namely T (2) αx . Thus, T βy . Using similar arguments, we find that Adding up the four fields, leads to: zy )}, which does not have any component in the xy plane. If the interaction in question is dipole, this is similar to the well-known derivation of the NMR hyperfine field on As. In case of the DM interaction, the last term also cancels and the entire field is zero.

B. Experiment
Single crystals of EuFe 2 As 2 were grown using three different techniques, namely, Fe-As self-flux [6], Sn-flux and the Bridgman method [7]. The respective sample preparation methods were similar to those in the corresponding references. Samples were characterized using energy-dispersive x-ray analysis and x-ray diffraction to confirm composition and structure. Afterwards, they were oriented using a Laue camera and subsequently cut along the tetragonal [110] T -direction with an electric spark erosion technique. Typical dimensions after cutting were approximately 2×1×0.1 mm 3 .
Magnetization measurements were performed in a standard Magnetic Property Measurements System (MPMS), while DC-transport data was recorded using a standard Physical Property Measurement System (PPMS) in the Montgomery geometry.

II. MODEL DETAILS
In the main text, we have determined the equilibrium phase diagram. Here, we shall derive the energies of each type of domain for all fields, so as to be able to deduce the energy cost/gain of any reorientation. Eq. 12 in the main text is based on this derivation. We will also present, for completeness, the energetics at ultra-high fields, not accessible in current experiments.
We start from Eq. 4 of the main text, and will first calculate the energy of the b-twin (H b), and then of the atwin, in all fields. We will consider separately two different regimes: first, small fields compared to H 2 defined in the main text, and, second, large fields comparable with, or larger than H 2 . In the former regime we can safely neglect the Fe spins' canting (φ 2 =φ 1 = π/2,α 2 =α 1 = 0, while in the latter we can use the fact that the Eu magnetization is already saturated and ϕ 1 = ϕ 2 = 0.
A. b-Domain: H b.

Small Fields
At small fields we have,φ 2 =φ 1 = π/2,α 2 =α 1 = 0 and ϕ = ϕ 1 = ϕ 2 , which, when inserted into Eq. (4) of the main text leads to with p = cos ϕ and E 0 = −2J −2D−J ⊥ −8K. Minimizing Eq. (1) with respect to p, leads for M H < 4(J ⊥ + 4K) to the equilibrium tilting angle and energy given in Eq. (6-7) of the main text. For sufficiently small magnetic fields this solution is always below E 0 and Eu 2+ moments can screen the field continuously. The tilting angle ϕ changes gradually according to At fields larger than H sat b = (4J ⊥ + 16K)/M (derived from p = 1) the moments are fully aligned with the external magnetic field and the energy changes to

Large Fields
At considerably higher fields the Fe moments will start to deflect (α = π/2−φ) from the ground state orientation, while the Eu 2+ moments are already saturated. In this case Eq. (4) of the main text yields where we have now introducedp = cosφ and E 0 . The equilibrium energy is now determined by minimizing with respect top, which then reads i.e. Eq. (3) with an additional term stemming from the canting of the Fe moments away from the easy axis.
B. a-Domain: H a

Small Fields
At small fields we have the following conditions:α 1 = ϕ 1 = 0,α 2 =φ 2 = π and ϕ = ϕ 1 = ϕ 2 , leading to again with p = cos ϕ and E 0 , which after minimizing w.r.t. p results for J ⊥ > 4K and M H < 4J ⊥ − 16K in Eqs. (9-10) of the main text. We note the sign change in the denominator and the additional 8K-term. Due to this term this state will only be populated once H is large enough. At smaller fields the Eu 2+ moments will remain in their ground state configuration. The field at which the moments spin-flop H flop a = 8/M K(J ⊥ − 4K) is determined from the condition E min a = E 0 . In the spin-flop phase the Eu 2+ moments gradually rotate towards saturation according to However, if the biquadratic coupling K is strong compared to the interplanar coupling J ⊥ between Eu 2+ moments, they directly flip from the ground state into the fully saturated state, skipping Eq.

Large field
If the field is large enough the Fe moments can become subject to a spin-flop (but not a spin-flip, becausẽ J D ). As estimated in the main text this should happen around 130 T. We can easily incorporate this effect for the sake of completeness, by setting ϕ = 0 and α =φ, which leads toẼ

C. Eu spin angles at H1
In the main text we implied, for simplicity, that ϕ(H 1 ) = π/4 when the system changes from b-twins to a-twins at H 1 . However, this holds only in the limit of K → 0. At finite K the Eu spin angle changes discontinuously from a value ϕ > π/4 to a value ϕ < π/4, as can be seen when inserting the expression for H 1 , Eq. (11) in the main text, into Eqs. (2) and (4). This corresponds to a discontinuous decrease of ϕ from roughly 55 o to 25 o for our averaged findings of J ⊥ and K. The respective jump in magnetization M (H 1 ) can easily be calculated using Eq. (15) of the main text and is in agreement with the experiment.

D. Detwinning Dynamics
The full energy map for both types of domains in all fields allows us to calculate the energy cost/gain of switching from a-twin domains to b-twin domains, dE, as a function of the reduced field, M H/J ⊥ , and reduced Table SI. The energy difference dE in different regimes, assuming 12K/J ⊥ < 1. In case of 12K/J ⊥ > 1 and 8K/J ⊥ < 1, H1 changes to H sat 1 , which lies in the interval {H sat a , H sat b } (Fig. S3) Since the domain dynamics are driven by the energy difference dE between the domains we show in Fig. S4 dE = E a − E b on a semi-log color map as a function of K/J ⊥ and M H/J ⊥ . The dotted line represents K/J ⊥ = 0.05, for which dE was depicted in the main text (Fig. 3). Positive and negative values correspond to the b-twin and a-twin domains being the ground state. All fields and regimes are summarized in Table SI for the case 12K/J ⊥ < 1.