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

TL;DR

This paper develops a microscopic theory explaining magnetic detwinning in EuFe₂As₂, revealing how biquadratic coupling between Fe and Eu spins enables ultra-low field detwinning and domain reversal.

Contribution

It introduces a comprehensive microscopic model accounting for non-mechanical magnetic detwinning phenomena in EuFe₂As₂, including the effects of Fe-Eu spin coupling.

Findings

Quantitative explanation of ultra-low detwinning field (~0.1 T).

Understanding of domain orientation reversal at ~1 T.

Restoration of low-field domain orientation above 10-15 T.

Abstract

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.

Equations

$${\mathcal{H}}_{\text{Fe}}=-\tilde{D}\sum _i{f}_{i,a}^2+\tilde{J}\sum _{⟨ij⟩}{\mathbf{f}}_i\cdot {\mathbf{f}}_j-\tilde{M}\sum _i{\mathbf{f}}_i\cdot \mathbf{H}$$

$${\mathcal{H}}_{\text{Eu}}=J_{\parallel }\sum _{⟨ij⟩}{\mathbf{e}}_i\cdot {\mathbf{e}}_j+J_{\perp }\sum _i{\mathbf{e}}_i\cdot {\mathbf{e}}_{i^{+}}-M\sum _i{\mathbf{e}}_i\cdot \mathbf{H}$$

$$\mathcal{H}=$$

$$-M\sum _i{\mathbf{e}}_i\cdot \mathbf{H}+J_{\perp }\sum _i{\mathbf{e}}_i\cdot {\mathbf{e}}_{i^{+}}$$

$$-K\sum _{\alpha ,n}{\left({\mathbf{f}}_{\alpha }\cdot {\mathbf{e}}_{\alpha +n}\right)}^2,$$

$$\mathcal{E}=$$

$$-MH(\cos {\phi }_1+\cos {\phi }_2)/2+J_{\perp }\cos ({\phi }_2+{\phi }_1)$$

$$-2K[{\cos }^2({\phi }_2-{\tilde{\phi }}_2)+{\cos }^2({\phi }_2-{\tilde{\phi }}_1)+{\cos }^2({\phi }_1-{\tilde{\phi }}_2)+{\cos }^2({\phi }_1-{\tilde{\phi }}_1)],$$

$$E_b=-MHp(2J_{\perp }+8K)p^2+E_0,$$

$$p_b^{\text{min}}=\frac{MH}{4J_{\perp }+16K}$$

$$E_b^{\min }=-\frac{(MH{)}^2}{8(J_{\perp }+4K)}+E_0.$$

$$E_a=-MHp+(2J_{\perp }-8K)p^2+E_0+8K$$

$$p_a^{\min }=\frac{MH}{4J_{\perp }-16K}$$

$$E_a^{\min }=-\frac{(MH{)}^2}{8(J_{\perp }-4K)}+E_0+8K.$$

$$H_1=\frac{2}{M}\sqrt{2(J_{\perp }^2-16K^2)}$$

$${\tilde{E}}_b=$$

$$-MH+E_0+2J_{\perp }+8K,$$

$${\tilde{p}}_b^{\text{min}}=\frac{\tilde{M}H}{2(2\tilde{J}+\tilde{D}-4K)},$$

$$dE=+\frac{(\tilde{M}H{)}^2}{2(2\tilde{J}+\tilde{D}-4K)}-8K.$$

$$dE=\left\{\begin{array}[c]{ll}\frac{M^{2}H^{2}}{8(J_{\perp}+4K)},&\{0,H_{a}^{\text{flop}}\}\\ K(8-\frac{M^{2}H^{2}}{J_{\perp}^{2}-16K^{2}}),&\{H_{a}^{\text{flop}},H_{a}^{\text{sat}}\}\\ -MH+2J_{\perp}+\frac{M^{2}H^{2}}{8(J_{\perp}+4K)},&\{H_{a}^{\text{sat}},H_{b}^{\text{sat}}\}\\ -8K+\frac{\tilde{M}^{2}H^{2}}{2(2\tilde{J}+\tilde{D}-4K)},&\{H_{b}^{\text{sat}},H_{2}\}.\end{array}\right.$$

$$J_{\perp }$$

$$K$$

$$n(dE)=\left\{\begin{array}[c]{lr}n_{0}\cdot\frac{1}{2}\operatorname{erfc}\left(\frac{\log(dE/\Delta)}{\sqrt{2}}\right),&dE>0\\ dn\cdot\frac{1}{2}\operatorname{erfc}\left(-\frac{\log(-dE/\Delta)}{\sqrt{2}}\right)+n_{0},&dE<0\end{array},\right.$$

$$n(H)=\frac{\Delta L}{L_0}\left(\frac{b}{a-b}+n_0\right),$$

$$n_0=\frac{a\cdot n(H)-b[n(H)+\Delta L/L_0]}{(a-b)(1+\Delta L/L_0)},$$

$$H_0=\frac{2}{M}\sqrt{2\Delta (J_{\perp }+4K)}.$$

$$M(H)=n(H)\cdot M_a+(1-n(H))\cdot M_b,$$

$$M_a$$

$$M_b$$

Full text

**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 familiesFisher et al. (2011); Tanatar et al. (2009). 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 computationallyYaresko et al. (2009) and later shown to provide the only physically meaningful description of spin dynamics in FeBSWysocki et al. (2011). There is growing evidence that it is not limited to FeBS, but occurs also in other itinerant systemsZhang et al. (2017).

Mechanical straining is not the only way to detwin FeBS. It was shown that a static magnetic field of $\sim 15$ T leads to partial detwinningChu et al. (2010) and pulsed fields of $\sim 30$ T to nearly complete detwinningRuff et al. (2012). 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 $\mu$eV. Against this background, it came as a complete surprise when it was discoveredXiao et al. (2010); Zapf et al. (2014) 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 non-intrusive (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}\sim 0.5$ T. Then, at $H_{1}\sim 1$ T domains spontaneously rotate in-plane by 90${}^{o},$ and at $H_{2}\sim 10$ T start to turn back to the detwinned state that was initially generated at $H\lesssim 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 BaFe2As${}_{2}.$ The minimal approximation is the nearest-neighbor (n.n) Heisenberg model with single site anisotropy:

[TABLE]

where $i,j$ label Fe sites, $\mathbf{f}_{i}$ is the unit vector directed along the Fe magnetic moment at site $i,$ $\tilde{M}$ is its absolute value, and the summation is over all inequivalent n.n bondsnot . $\mathbf{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.$, $\tilde{D}>0)$. Neutron scattering provides estimates for the parameters as $\tilde{J}\sim 30$ meV, $\tilde{D}\approx 0.25$ meV Wang et al. (2013). Minimizing Eq. (1) with respect to **f, **for $\mathbf{H||}b$ (no linear susceptibility appears for $\mathbf{H||}a)$ we observe that $H$ generates a canting of the Fe spins away from the $a$ axis by an angle $\tilde{\alpha}=\sin^{-1}(\tilde{M}H/4\tilde{J})$, which corresponds to an energy gain of $E=(\tilde{M}H)^{2}/4\tilde{J}$ 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 when $\tilde{J}>\tilde{D}.$ Using $\tilde{M}\approx 1\,\mu_{\text{B}}$ and the above parameters, we estimate $\tilde{H}^{\text{flop}}=2\sqrt{2\tilde{D}\tilde{J}}/\tilde{M}\approx 11\,\textnormal{meV}\approx 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 90o around $c$, keeping the moments aligned along the magnetic easy $a$ axis. This process is associated with an unknown energy barrier $\Delta.$ Given that in the field of $\sim 15$ T the Stanford group has observed partialChu et al. (2010), and in $\sim 30$ T nearly full detwinningRuff et al. (2012), we deduce $\Delta\sim E=(\tilde{M}H)^{2}/4\tilde{J}\approx 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 $\sigma=$ 6–20 MPaJiang et al. (2013); Blomberg et al. (2012). Assuming $\sigma=10$ MPa and taking the elastic modulus to be 10 GPaShein and Ivanovskii (2009); Yoshizawa and Simayi (2012), we derive a strain of $\varepsilon\approx$ 0.1%. Finally, using the calculated dependence of the total energy on the microscopic strainJesche et al. (2008), we find that this strain corresponds to $\Delta\approx$ 0.01  meV/f.u, which is in agreement with our magnetic estimate.

Yet another estimate can be obtained by considering the stress $\sigma$ on the unit cell during mechanical detwinning and calculate the energy associated with the displacement $\delta$ from $a\parallel\sigma$ to $b\parallel\sigma$. Using the reported lattice parameters of EuFe2As2Xiao et al. (2009), we arrive at $\delta=$ 1.6$\times 10^{-12}$ m. Assuming $\sigma=6$ MPa, this leads also to an energy of $\Delta\approx 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 EuFe2As${}_{2},$ given that the ordering temperature in the Eu sublattice is much smaller, $T_{\textnormal{N}}=19$ K, and the moment much larger, $M=7\,\mu_{\textnormal{B}},$ than for the Fe sublattice. Such a phenomenology was adapted in Ref. Zapf et al. (2014) 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,

[TABLE]

where $i,j$ label Eu sites, $i^{+}$ the opposing Eu site in the next layer, $\mathbf{e}_{i}$ the unit vectors directed along the Eu magnetic moment $M$ at site $i$, and the ferromagnetic $J_{\parallel}<0$ and antiferromagnetic $J_{\perp}>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, Eu2+ 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 EuFe2As2 cannot be explained within a simple Heisenberg model. Indeed, the magnetic susceptibility of EuFe2As2 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_{\text{B}}T_{\text{CW}}=\frac{3}{7}(4J_{\parallel}-2J_{\perp})$ (the quantum prefactor $3M/(M+2)$ would have been 1 if $M$ were $1\,\mu_{\text{B}}$, but for $M=7\,\mu_{\text{B}}$ it becomes $3/7$) The effective moment has been determined to be $M=7\,\mu_{B}$ with $T_{\text{CW}}\approx-20$ K Jiang et al. (2009). 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_{\perp}$ and $J_{\parallel}.$ 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_{\perp}^{2}/U,$ and Schrieffer-Wolff drivenGlasbrenner et al. (2014) double exchange, proportional to ($t_{\textnormal{Eu-Fe}}^{2}/U)^{2}N(0),$ with $N(0)$ the density of states, and $t_{\perp}$ and $t_{\textnormal{Eu-Fe}}$ the effective Eu-Eu and Eu-Fe hopping across the planes. For $U-J=0$ we get $J_{\perp}=0.14$ meV (1.6 K) and $J_{\parallel}=-2$ meV (23 K), while for $U-J=5$ eV we find $J_{\perp}=0.26$ meV ($3$ K) and $J_{\parallel}=$-0.8 meV (9 K). The LDA+$U$ results, which we believe are closer to reality, correspond to $T_{\textnormal{CW}}\sim-13$ K, and $T_{\textnormal{N}}\approx 4.27J_{\parallel}/[3.12+\log(J_{\parallel}/J_{\perp})]=9$ K for the Heisenberg modelYasuda et al. (2005). The ratio $T_{\textnormal{CW}}/T_{\textnormal{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_{\perp}/J_{\parallel}$ ratio to $\sim 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\approx 1$ K, because in such a case one has to replace $J_{\perp}$ with $J_{\text{eff}}=J_{\perp}+D+\sqrt{D^{2}+2J_{\perp}D}$ Irkhin and Katanin (1998, 1999), which leads to a sufficient increase of the Néel temperature and consequently to $T_{\textnormal{CW}}/T_{\textnormal{N}}$ = 1.13. But, as we have argued above, Eu2+ 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 EuFe2As2, 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_{\text{eff}}\sim 0.6$ K. Substituting this into the formulas above we get $T_{\textnormal{CW}}/T_{\textnormal{N}}\approx 1.18$, which is within the error margin of the experimental ratio. Of course, given the model character of these calculation, a discrepancy of $\approx 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 behaviorGlasbrenner et al. (2015). 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:

[TABLE]

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\sim 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 theory, concentrating on rationalizing the actual experimental observations.

Domain energetics We start our discussion in the *low field *regime, that is $0<H\lesssim H_{1}\sim$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 $\mathbf{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 $\varphi$ and $\tilde{\varphi}$ with respect to the external field. Figure 1(a) is drawn for a b-twin domain, i.e., $a\perp H$. Equation (3) can then be rewritten in terms of these angles (per one formula unit) as

[TABLE]

where $\alpha$’s are the Fe angles measured from the magnetic easy $a$-axis.

For the b-twin domains and low fields, i.e. $\tilde{\varphi}_{2}=\tilde{\varphi}_{1}=\pi/2,$ $\varphi_{1}=\varphi_{2}=\varphi,$ $\tilde{\alpha}_{2}=\tilde{\alpha}_{1}=0$ this yields

[TABLE]

where we have expressed everything in terms of $p=\cos\varphi$ and $E_{0}=-2\tilde{J}-2\tilde{D}-J_{\perp}-8K$. The equilibrium tilting angle and energy are given by

[TABLE]

[TABLE]

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 $\pi/4<\varphi<\pi/2,$ Fe spins prefer to orient perpendicular to the field ($\tilde{\varphi}=\pi/2)$, and the b-twin domain is always lower in energy. When $\varphi$ becomes smaller than $\pi/4,$ it becomes more favorable to orient Fe spins parallel to the field $(\tilde{\varphi}=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 only $\tilde{\varphi}=\pi/2$, or, in case of the a-twin, 0 is allowed, any intermediate value of $\tilde{\varphi}$ is severely punished by the Fe-Fe exchange and single-site anisotropy).

After this reorientation has occurred, the total energy is expressed as

[TABLE]

and consequently

[TABLE]

[TABLE]

Thus the first reorientation field $H_{1}$ is defined by $E_{a}^{\min}=E_{b}^{\min},$ or

[TABLE]

At the field $H_{a}^{\text{sat}}=4(J_{\perp}-4K)/M$ the angle $\varphi$ 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 $-MH)$, 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_{a}^{\text{sat}}$ 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\gtrsim 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_{a}^{\text{sat}}$, 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 exchange $\tilde{J}$ and because the Fe site-anisotropy is much larger than the biquadratic coupling).

In this case, the total energy (since now $\varphi=0,$ $\tilde{\alpha}=\pi/2-\tilde{\varphi})$ is

[TABLE]

where we have now used $\tilde{p}=\cos\tilde{\varphi}$ and $E_{0}$. Minimizing with respect to $\tilde{p}$ yields:

[TABLE]

with the energy gain compared to Eq. (8) with $p=1$ being

[TABLE]

This expression changes sign at $H_{2}=4\sqrt{2\tilde{J}K+\tilde{D}K-4K^{2}}/\tilde{M}\approx 4\sqrt{2\tilde{J}K}$/$\tilde{M}.$ 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

[TABLE]

In the second case $dE>0$ changes to $dE<0$ at $H_{1}$, which lies between $H_{a}^{\text{flop}}$ and $H_{a}^{\text{sat}}$. The fields $H_{1},H_{a}^{\text{sat}}$ and $H_{2}$ have been derived above. The transformations associated with $H_{a}^{\text{flop}}$ and $H_{b}^{\text{sat}}$ occur in “wrong”, or minority domains, which are thermodynamically unstable, but occur kinetically. Their expressions, as derived in the Supplement, are $H_{b}^{\text{sat}}=4(J_{\perp}+4K)/M$ and $H_{a}^{\text{flop}}=8/M\sqrt{K(J_{\perp}-4K)}.$

The coupling constants $J_{\perp}$ 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 dataZapf et al. (2014), 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]_{\text{T}}$ direction, a roughly linear magnetization $M(H)$ was reported in Ref. Zapf et al. (2014), interrupted only by two pronounced jumps around 0.1 T and 0.6 T. Saturation sets in above 0.9 T with $M\approx 5\,\mu_{\text{B}}$. Due to the uncharacteristically small saturated moment, which does not agree well with previously published dataJiang et al. (2009) and the theoretically expected value of $M=7\,\mu_{\text{B}}$, we have re-measured the magnetization of a EuFe2As2 sample from Zapf et al. Zapf et al. (2014). 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\,\mu_{\text{B}}$, in good agreement with the theoretical expectations and Ref. Jiang et al. (2009).

The step-like increase of the magnetization around $0.65\,\mu_{\text{B}}\text{T}$ was associated with a spin-flip transition of the Eu2+ moments in Ref. Zapf et al. (2014). 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 b-twin domain can be extrapolated to $H_{b}^{\text{sat}}\approx 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\gtrsim 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:

[TABLE]

which yields $J_{\perp}=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_{\perp}=0.121(24)$ meV, $K=0.0072(22)$  meV and $J_{\perp}/K\approx 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 semi-log plot as a function of the reduced variable $MH/J_{\perp}$ for a representative ratio of $K/J_{\perp}=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_{\perp}$, $MH/J_{\perp}\}$ 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 $\Delta$ associated with the reorientation, i.e. the energy difference $dE=E_{\text{a-twin}}-E_{\text{b-twin}}$ between the two twin variants needs to exceed a certain threshold before reorientation occurs. In the following, we will assume $\Delta$ for various domain walls to be log-normal distributed, i.e. $\log{\Delta}$ is normalLimpert et al. (2001). This is a typical distribution e.g. of grain sizes in polycrystalline matterKurzydlowski and Ralph (1995). The cumulative distribution function $F_{X}(x)$ of a positive log-normal distributed variable $X$ is given by $F_{X}(x)=\frac{1}{2}\operatorname{erfc}\left(-\frac{\log{x}-\mu}{\sigma\sqrt{2}}\right)$, with the location and the scaling parameters $\mu$ and $\sigma$. 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:

[TABLE]

with $n_{0}$ the fraction of a-twins at $H=0$ and $dn=n_{\text{sat}}-n_{0}$ the difference between the saturated a-twins $n_{\text{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 soonMaiwald et al. . In the following we will use magnetostriction (MS) data by Zapf et al. Zapf et al. (2014), which agree with the preliminary data of Ref. Maiwald et al. .

The MS is defined as $\Delta L(H)/L_{0}=(L(H)-L_{0})/L_{0}$, where $L_{0}$ and $L(H)$ are the initial and field-dependent average unit cell length, respectively. They can be expressed by $L(H)=n(H)a+[1-n(H)]b$ and $L_{0}=n_{0}a+(1-n_{0})b,$ again with the (initial) domain population $n_{0}$ at $H=0$ and orthorhombic lattice parameters $a$ and $b$. Solving for $n$ leads to

[TABLE]

from which follows

[TABLE]

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 dataMaiwald et al. indicate a domain distribution with $n(H_{0})=0.06$. Furthermore, the significant pressure of the dilatometer in field direction, 1.35 MPaZapf et al. (2014) aids the $n=0$ alignment at small magnetic fields (but hinders it at larger fields, when $n\rightarrow 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 $\Delta=0.01$ meV. Among the other investigated samples $\Delta$ 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=\Delta$, through

[TABLE]

For the determined $\Delta$-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 $\Delta=0.004$ meV, $n_{i}=0.22$ and $dn=0.71$. The total magnetization is given by

[TABLE]

where $M_{a}$ and $M_{b}$ are the magnetization of the sublattices given by $M\cos{\varphi}$, i.e.

[TABLE]

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 Eu2+ 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 EuFe2As2, also features large $S$=7/2 spin-only moments, or the recently discovered Eu-based 1144-systems. Furthermore, the theory also captures the physics of the high-field ($H\gtrsim 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. Zapf et al. (2014)), 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.