Magnetostriction Reveals Orthorhombic Distortion in Tetrahedral Gd-compounds

We report detailed thermal expansion and magnetostriction experiments on GdCoIn$_5$ and GdRhIn$_5$ single crystal samples that show a sudden change in the dilation at a field B$^\ast$ for temperatures below the N\'eel transition temperature TN. We present a first-principles model including crystal-field effects, dipolar and exchange interactions, and the dependence of the exchange couplings with lattice distortions in order to fully account for the magnetostriction and magnetic susceptibility data. The mean-field solution of the model shows that a transition between metastable states occurs at the field B$^\ast$. It also indicates that two degenerate phases coexist in the sample at temperatures below TN. This allows to explain the lack of observation, in high resolution x-ray experiments, of an orthorhombic distortion at the N\'eeel transition even though the magnetic structure breaks the tetragonal symmetry and the magnetoelastic coupling is significant. These conclusions could be extended to other tetragonal Gd-based compounds that present the same phenomenology.

We report detailed thermal expansion and magnetostriction experiments on GdCoIn5 and GdRhIn5 single crystal samples that show a sudden change in the dilation at a field B for temperatures below the Néel transition temperature TN . We present a first-principles model including crystal-field effects, dipolar and exchange interactions, and the dependence of the exchange couplings with lattice distortions in order to fully account for the magnetostriction and magnetic susceptibility data. The mean-field solution of the model shows that a transition between metastable states occurs at the field B . It also indicates that two degenerate phases coexist in the sample at temperatures below TN . This allows to explain the lack of observation, in high resolution x-ray experiments, of an orthorhombic distortion at the Néel transition even though the magnetic structure breaks the tetragonal symmetry and the magnetoelastic coupling is significant. These conclusions could be extended to other tetragonal Gd-based compounds that present the same phenomenology.

I. INTRODUCTION
Rare-earth magnetic compounds are among the strongest permanent magnets and present the highest magnetostrictive responses ever recorded. These remarkable properties stem from the large magnetic moments of the rare-earth ions with partially filled f-shells and the magnetic anisotropy associated with crystal-field effects and spin-orbit couplings. The magnetic structure of these compounds is mainly determined by the Ruderman-Kittel-Kasuya-Yosida (RKKY) exchange interactions between the magnetic moments of the rare earth ions and by the crystal field, which dominate over the dipolar interaction.
In many materials the magnetoelastic couplings lead to a significant spontaneous lattice distortion concomitant with the magnetic order at zero applied magnetic field. In some systems, the associated changes in the lattice constants can be as large as a few percent [1]. When the symmetry of the magnetic order is lower than the lattice symmetry, a reduction of the latter is expected at the magnetic transition. This can occur, e.g., when the magnetoelastic couplings generate changes in the lattice parameters which do not preserve the lattice symmetry.
In the rare earth series, Gadolinium stands as different. In solid state compounds it is generally found as a trivalent ion Gd 3+ for which Hund's rules indicate the maximum spin allowed S = 7/2, and zero angular momentum L = 0. As a consequence of the latter, the magnetic coupling to the lattice via crystal-field effects is expected to be weak. The magnetic anisotropy observed in Gd compounds is therefore usually attributed to the dipolar interaction [2]. The dependence of the exchange couplings on the relative distance between ions can, however, give rise to large magnetoelastic couplings.
Although a coupling between lattice and magnetic orders is observed on most Gd compounds there are few re-ports of lattice symmetry breaking at the magnetic transition (see Ref. [1] for a review). To the best of our knowledge, such lattice symmetry breaking in Gd compounds has only been confirmed for the ferromagnetic GdZn compound where the lattice distorts from cubic to orthorhombic at the magnetic transition [3].
Many Gd compounds with tetragonal lattice show antiferromagnetic (AFM) order [2,4]. Competing magnetic couplings can result in an AFM order other than the trivial G-AFM, where every pair of first neighbors is antiparallel. For instance, a second-neighbor AFM coupling can lead to a C-AFM order where chains of parallel-aligned moments order antiparallel (antiferromagnetically) between them (see Fig. 1) [5,6]. If these chains are aligned along the basal plane, it is expected that below the AFM ordering temperature T N , the lattice lowers its symmetry to orthorhombic [7]. Although symmetryconserving distortions at the AFM transition have been easily detected [4], high-precision x-rays experiments (upto |a − b|/a ∼ 2 × 10 −4 ) do not show any difference between a and b lattice parameters [7,8]. Such an intriguing absence of lattice symmetry breaking at the Néel transition in Gd-based AFM systems has been referred to as the magnetoelastic paradox [7].
In this work we present very sensitive dilation experiments across the Néel transition in GdCoIn 5 (T N ≈30 K) [9] and GdRhIn 5 (T N ≈40 K) [10]. Both systems show, for temperatures below T N , an abrupt change of the longitudinal linear forced magnetostriction in an external field of ∼ 1 Tesla. We perform detailed calculations which indicate that the dipolar interaction, the crystal electric field, and the strain dependence of the magnetic exchange couplings are all essential to account for the observed lattice distortions and magnetic structure. We also show that the dilation data is compatible with the existence of a tetragonal to orthorhombic distortion of the lattice at the Néel transition. The predicted or-thorhombic distortions are below the x-ray resolution and result from the competition between the different magnetoelastic couplings.
The rest of this paper is organized as follows. Section II presents the main experimental and theoretical results for the magnetostriction and thermal expansion data for the GdCoIn 5 and GdRhIn 5 compounds. Section III presents the model, while Sec. IV presents the numerical simulations and their interpretation. Finally in Sec. V we present the conclusions.

II. MAIN RESULTS
We present below the main experimental and theoretical results on the magnetostriction and thermal expansion data for the GdCoIn 5 and GdRhIn 5 compounds. Both materials crystallize in the tetragonal HoCoGa 5 structure (see Fig. 1). The magnetic structure inferred by x-ray experiments in GdRhIn 5 is of C-type [8] (see Fig. 1). High quality single crystals of GdCoIn 5 (GCI) and GdRhIn 5 (GRI) were grown by the self-flux technique and characterized as described elsewhere [9]. The specific heat and the magnetic susceptibility data are compatible with S = 7/2 spins at the Gd 3+ ions coupled by exchange interactions mediated by the conduction electrons. The main difference between GdCoIn 5 and GdRhIn 5 is the larger exchange coupling along the c-axis in the latter which leads to a higher T N [11].
Platelet-shaped crystals of typical size 1 × 1 × 0.4 mm 3 were selected for the dilation experiments which were performed with a high resolution (∆L ≤ 1Å) capacitive dilatometer [12]. All dilation experiments under magnetic field were carried out in the longitudinal configuration, i.e. with the magnetic field B parallel to the sample dimension L being measured. Figure 2 summarizes the most important experimental observations of this work. It displays the forced magnetostriction (lattice dimension change driven by an external magnetic field) along theâ-axis for both GCI and GRI. At temperatures below T N , ∆L a /L a shows a sudden increase at an in-plane field B ∼ 1 T. Above this field, ∆L a /L a becomes field independent for GCI and shows a weak increase with increasing magnetic field for GRI. No hysteresis effects are observed. At temperatures above T N the magnetostriction becomes negligible and no sudden change is observed. This type of change in ∆L for T < T N is usually seen on ferromagnetic materials and attributed to the change from the zero magnetization to the saturated magnetic state (see chapter 1 on Ref. [13]). As we will see in the following Section, in these antiferromagnetic systems an applied magnetic field induces a spin-flop transition (at the field B ) to a metastable state where the spins point mainly along theĉ-axis [14,15]. To account for the value B and the absence of hysteresis both the dipole interaction and the effect of the crystal field in second order perturbation theory needs to be considered in the calculations. The theoretical results (dashed lines in Fig. 2), were obtained using fitting interaction parameters constrained to the range of estimated values (see Sec. III). Along theĉ-axis, the field dependence of the forced magnetostriction (MS) is quadratic at all temperatures as it is shown in Fig. 2 for GCI. GRI mimics this behaviour (not shown here). The spontaneous and the forced magnetostrictions are a measure of the strength of the magnetoelastic couplings which can be extracted from the thermal expansion data. The main panel of Fig. 3 shows theâ-axis thermal expansion ∆L a /L a of GCI at zero field and at B = 5 T > B * . A small kink signposts the magnetic transition at T N . In the paramagnetic state (T > T N ), ∆L a /L a ∝ T 2 as seen in the inset of Fig. 3. This quadratic non-magnetic thermal-expansion background is also included in the main panel and it extrapolates at low temperatures to ∆L a (non-magnetic, T → 0)/L a = −2.5 × 10 −5 . We take the quadratic fit as the non-magnetic thermal expansion, and the difference between the zero field thermal expansion curve and that fit corresponds to the spontaneous MS. Accordingly, the difference between the finite field and the zero field data is the forced MS. The total MS is the sum of the forced and the spontaneous MSs and can be obtained substracting to the finite field curve the non magnetic fit. According to this analysis, theâ-axis spontaneous MS is positive in GCI, resulting in a zero temperature expansion ∆L a (B = 0, T → 0)/L a − ∆L a (non-magnetic, T → 0)/L a = 2.5 × 10 −5 , while the forced MS is also positive giving ∆L a (B =5 T, T → 0)/L a − ∆L a (B = 0, T → 0)/L a = 2 × 10 −5 . The zero-temperature spontaneous and forced MS's of both GCI and GRI are summarized in Table I. An equivalent analysis can be performed from the thermal expansion data along theĉ-axis. In this case, however, the spontaneous MS has the opposite sign The theoretical description of the MS data needs also to account for the pronounced anisotropy observed in the magnetic susceptibilities below T N as it was reported pre- viously [10,16]. Interestingly, the observed difference between theâ-axis andĉ-axis susceptibilities is rapidly suppressed as the magnetic field is raised above ∼ 1 T [16,17].

III. MAGNETOELASTIC HAMILTONIAN
In this Section we present the model used to describe the experimental data. We present the different magnetic interaction terms in the Hamiltonian and analyze their coupling to the lattice degrees of freedom.

A. Exchange interactions
The magnetic exchange interactions between localized magnetic moments at the Gd 3+ ions were determined by magnetic susceptibility and specific heat experiments combined with first principles calculations [9,11]. The exchange interactions up to the fifth nearest neighbor are presented in Fig. 4. These couplings do not determine unambiguously the magnetic ground state as the energy only depends on the relative orientation of the magnetic moments. ). Figure 4 presents four different magnetic moment arrangements having the same exchange energy.
The exchange couplings J i , shown on Fig. 4, are modified when the lattice is distorted where δa, δb, and δc are uniform lattice distortions along theâ,b, andĉ axis, respectively. The magnetic exchange interaction between magnetic moments S i at the Gd 3+ ions can be written as where the J i,j couplings are equal to J 0 for first nearest neighbors, J 1 for second nearest neighbors, etc. Although this interaction and its dependence with the lattice distortions is the largest, it is not enough to explain the observed ground state configuration, nor the magnetoelastic data or the magnetic anisotropy. For example, the Hamiltonian given by Eq. (3) leads to the same energy for the C-AFM aa and the C-AFM ac configurations. As we show below, the dipolar interactions break this degeneracy.

B. Dipolar interactions
The dipolar interactions have an explicit dependence with the distance between spins: where K is Kelvin scale unit, a B is the Bohr radius, r i is the position of the i-th spin and Note that here and in what follows we take k B = 1 and use K for the energy units. H D introduces a magnetic anisotropy. Since the distance between nearest-neighbor spins is larger along theĉ-axis than along theâ orb axes, the dipolar configurations with the lowest energy have the spins in the a-b plane. Still those states remain highly degenerate as a continuum of configurations with in-plane second nearest neighbors antiparallel have the same exchange and dipolar energies [see e.g. in Figs. 4a), 4d), and 6]. The experimentally observed order are however the ones shown in Figs. 4a and 4d). The first state, Fig. 4a), becomes the lowest lying state under a distortion a → a + δa and b → b + δb such that the lattice parameters a and b become different (orthorhombic distortion). The uniaxial magnetic anisotropy along theĉ-axis introduced by the dipolar interaction is however larger than what is inferred from the value of B (see Appendix A). An additional source of magnetic anisotropy, which is due to crystal-field effects, needs to be considered to explain the value of B .

C. Crystal-field effects
The Gd 3+ ion is, according to Hund's rules, in a 4f 7 state with L = 0. The spin-orbit i σ i . l i coupling, however, mixes this L = 0, S = 7/2 state with a higher energy multiplet with L = 1, S = 5/2 and J = 7/2 (see Appendix B for details), which is affected by the tetragonal crystal field. As a consequence crystal field (CF) effects, although small, are present and need to be considered.
The total crystal-field effect can be written as where S i is the component of the i-th spin along theˆaxis. The first term in Eq. (5) is the intrinsic CF [18] and the last term is induced by distortions between the a and b lattice parameters (see Ref. [19]). The contribution due to c deformations is negligible with respect to the intrinsic contribution B 2 S 2 c . The latter combined with the dipolar contribution to the anisotropy determine B .

D. Elastic energy
The elastic energy for a uniform distortion can be approximated as where the elastic constants C ab el ∼ C c el ∼ C el = 70000K/Å 2 are estimated from the elastic properties of materials of the same family of compounds [20,21] and Density Functional Theory (DFT) results (see Appendix C).
To evaluate this energy we approximate the large S = 7/2 spins on the Gd 3+ ions with classical magnetic moments.
We consider a lattice of L × L × L sites. We find that L 12 is enough to obtain L-independent results. We minimize the energy considering uniform deformations δa, δb, δc, and spin rotations restricted to magnetic orders that preserve the 8-site magnetic unit cell of Fig.  4.
The coupling parameters used in the simulations of GdCoIn 5 and GdRhIn 5 are presented in Table II. The exchange coupling parameters were obtained from DFT calculations [11]. The remaining parameters were obtained fitting the magnetostriction and the magnetic susceptibility. Reference values for the rate of change of the exchange coupling with the lattice parameter changes were obtained from DFT calculations [? ]. Note that for simplicity only the variation of the largest exchange couplings (J 0 and J 1 ) with the lattice distortions was considered.
A. Ground state for zero magnetic field (B = 0) The total energy as function of δa and δb is shown on Fig. 5. Two degenerate minima are obtained with magnetic orders C-AFM aa (see Fig. 4) and C-AFM bb which is related to C-AFM aa by a rotation of the lattice by 90 • around theĉ-axis. In the C-AFM bb configuration the spins and the spin chains are along theb-axis. The distortions associated with these minima satisfy δa aa = δb aa and δa aa = δb bb , δb aa = δa bb .
The fact that δa aa = δb aa indicates that the ground state crystal symmetry is reduced from tetragonal to orthorhombic. The magnitude of these distortions is however below the precision of high resolution XRD.
From here on we assume that the real system is composed by a mixture of these two states and consider the average of both distortions δa = (δa aa + δa bb )/2 (8) δb = (δb aa + δb bb )/2 (9) to compare with the experimental results. Assuming a homogeneous distortion of the lattice we have ∆L a /L a = δa/a. In the zero field case (B = 0) we haveδa/a ≡ δb/b, which corresponds to the spontaneous MS, since the model does not consider non-magnetic distortions. These distortions are presented in Table III and show a good agreement with the experimental results (see Table  I).

B. Magnetostriction
To analyze the effect of an external magnetic field on the striction we include the Zeeman coupling Under a magnetic field along theâ-axis, the energy of the C-AFM aa spin configuration remains unchanged while the energy of the C-AFM bb configuration [see Fig. 4  d)] is reduced. A large energy barrier separates however these two configurations. To change from C-AFM aa to C-AFM bb at zero field, the distortions need to be interchanged δa ↔ δb, the spins rotated 90 • interchanging nearest neighbour correlations from antiferromagnetic along theâ-axis to ferromagnetic, and vice-versa along theb-axis. Although the C-AFM ac order has a larger energy than C-AFM bb , it is much closer in configuration space to C-AFM aa . The spin-spin correlations do not change and the distortions are only slightly modified as they are determined to a large extent by these correlations. The exchange coupling only depends on the relative orientation of the spins. Our numerical results show that for a field B > B the configuration CAFM ac has a lower energy than CAFM aa and a transition between these two metaestable phases occurs. The critical field B is determined by the dipolar energy and the intrinsic crystal field parametrized by B 2 . For these compounds the CF reduces the critical field as it lowers the energy of the CAFM ac configuration. At B = B an energy barrier separates the C-AFM bb and the C-AFM ac configurations. In our scenario, the regions of the sample which at zero field were in the C-AFM aa configuration, with distortions δa aa , δb aa have, at B = B , a sudden change to the C-AFM ac spin configuration with different distortions δa ac and δb ac (see Table III). The regions of the sample in the C-AFM bb configuration, however, remain in it and the field only tilts the spins slightly in its direction. This leads to a forced MS at a field B > B  Table I for GdCoIn 5 and GdRhIn 5 at B = 5 T.
For low fields, B < B , only the deformations associated with the C-AFM bb order change. This is reflected on a small change of the total deformation. Around the critical field there is a sudden change of the deformations following the change of the magnetic order from C-AFM aa to CAFM ac . Figure 6 presents the energy of the system for a path in the space of spin configurations that passes through the C-AFM aa , C-AFM ac , and C-AFM bb phases. To go from one configuration to the other, each spin is rotated by the same angle and the energy is minimized with respect to distortions δa, δb and δc. In the absence of an external field, there are two degenerate minima for the C-AFM aa and C-AFM bb configurations. While the energy E aa of the C-AFM aa configuration does not depend on the magnetic field, both the C-AFM ac and C-AFM bb lower their energy with increasing magnetic field. At a field B ∼ 0.85 T the C-AFM aa configuration becomes unstable and there is a transition to the C-AFM ac phase. This transition has a very small hysteresis loop of width ∼ 0.01 T which is not observed in the experiments. Even for fields B ∼ 16 T a barrier separates the C-AFM ac and C-AFM bb phases.
We present in Appendix A a minimal model that captures the main physical ingredients and allows to obtain the model parameters from the experimental results.
The magnetostriction with the magnetic field parallel to theĉ-axis is simpler to understand because, in this case, there is no spin-flop transition. As the magnetic field is increased, the local spins tilt in theĉ direction. This leads to a change in the spin-spin correlations of antiparallel spins where M (B) M S is the uniform magnetization along theĉ-axis and M S is its saturation value. Since M (B) ∝ B, this results, to leading order, in a c lattice parameter change δc ∝ B 2 .

V. SUMMARY AND CONCLUSIONS
We analyzed the magnetoelastic properties of GdCoIn 5 and GdRhIn 5 . We measured the thermal expansion and the longitudinal magnetostriction on single crystals of both compounds using a high resolution capacitive dilatometer and constructed from first principles a model to account for the observed data.
These compounds present a number of intriguing properties. The observed magnetic order below T N is a C-type antiferromagnet which has a lower point symmetry than the lattice, but the expected tetragonal to orthorhombic distortion is not observed in high resolution XRD experiments. The in-plane dilation presents a sudden change at a field B ∼ 1 T and temperatures below the Néel transition which is not observed along theĉ-axis (for fields along the same axis). Contrary to expectations for a spin-flop transition [22,23], no hysteresis effects are observed at B . Although crystal-field effects are expected to be negligible in the Hund's rule ground-state multiplet of the Gd 3+ ion, a magnetic anisotropy is clearly observed in these compounds.
To understand the observed magnetic structure and reproduce the magnetostriction and magnetic susceptibility data we find it necessary to consider the spin-spin exchange interactions and their dependence on the lattice distortions, the dipolar interactions, and the crystal-field effects due to the mixing of the terms 8 S and 6 P as a consequence of the spin-orbit coupling.
The exchange couplings and their dependence on lattice distortions were estimated from first principles DFT based calculations while the crystal-field model and parameters were obtained from second order perturbation theory. The final model parameters were obtained for each compound through a fitting procedure of the magnetostriction and magnetic susceptibility data. As a consistency check of the model parameters, we estimated the change in the Néel temperature with an applied hydrostatic pressure which shows and excellent agreement with the expected value from Ehrenfest's thermodynamic equations (see Appendix D).
Our model can fully account for the observed experimental data, including the observed spin-flop transition and the absence of evidence of tetragonal symmetry breaking in these compounds. The main assumption is that the magnetically ordered state is a spatially inhomogeneous mixture of the two possible degenerate ground states. This assumption is needed to explain the absence of asymmetry, in the magnetic susceptibility and magnetostriction, between theâ-axis and theb-axis measurements.
Interestingly, there are other examples in the literature of tetragonal Gd-based compounds that show a similar behaviour, i.e. antiferromagnetic order with magnetic anisotropy below T N and a sudden change of the forced magnetostriction under a moderate magnetic field: GdNi 2 B 2 C (Refs. [7,24]), GdAg 2 (Ref. [25]) and GdRu 2 Si 2 (Ref. [26]) among others. The model discussed here could apply also to these cases. In this appendix we solve, at the mean field level, a simplified model that captures the main magnetoelastic properties of the system in the ordered phase. We focus the analysis on the modifications of the lattice parameters in the a-b plane which present sudden changes when a strong enough longitudinal magnetic field is applied. In the absence of an external magnetic field and assuming a C-AFM correlated state the model reads: Here, the first two terms stem from the crystal field and the dipolar interaction. The third and fourth terms are due to the dependence of the exchange interactions on the lattice parameters. The fifth term is the elastic energy and in the last term, λ is a Lagrange multiplier included to enforce the spin normalization | S| 2 = 1. The parameters of the model are given by: is the difference in dipolar energy between the C-AFM ac and C-AFM aa phases. γ 2 = A is given by the variation of the crystal field with lattice distortions [see Eq. (5)].
db )S 2 and β 2 = dJ1 da S 2 are given by the variation of the nearest neighbor and in-plane diagonal magnetic couplings with the lattice distortions.
There are four sets of distortions and spin projections that satisfy ∂F/∂ = 0 for all ∈ {δa, δb, S a , S b , S c , λ}: The ground state is doubly degenerate (C-AFM aa and C-AFM bb ) while C-AFM ac and C-AFM bc are degenerate higher energy states since γ 1 > 0 and γ 2 β 1 > 0. C-AFM ac (C-AFM bc ) is however unstable with respect to a tilt of the spins, increasing S a (S b ) . This fact preserves the spin-spin correlations (see Fig. 6 in the main text): A magnetic field along theâ-axis produces a magnetization Mâ in the states C-AFM bb , C-AFM ac and C-AFM bc by tilting the spins along the same axis. This leads to a reduction of the energy of these states by ∼ B 2 a /T N compared to the C-AFM aa state which remains unchanged. At fields larger than the energy of the C-AFM ac state becomes lower than the energy of the C-AFM aa state. Additionally the C-AFM aa state becomes unstable with respect to an increase in S c : (A7) This explains the transition between the metastable states C-AFM aa and C-AFM ac at B = B . Here E = 3K(1 − 2ν) is the Young modulus assuming an isotropic material where ν is Poisson's ratio. This results in anisotropic elastic constants along theâ −b andĉ axes, C ab el = Ec and C c el = Ea 2 /c, respectively. Using ν ∼ 0.22 (see below), we obtain C a el ∼ 71000K/Å 2 and C c el ∼ 27000K/Å 2 . The DFT calculations were performed using a supercell of 2 × 2 × 2 unit cells. On Fig. 7 we show the change on the total energy as one of the unit cell lattice parameters (a or c) is changed. The quadratic behaviour allows us to obtain the young modulus in both situations resulting C a el = 52000K/Å 2 ; C c el = 8000K/Å 2 .
Appendix D: Predicted effect of pressure on TN As a consistency check of the model parameters we calculate here the change of T N with external pressure. We find a good agreement between the model results and the inferred from Ehrenfest's thermodinamic equations (approximately ∼ 0.8K/GPa and ∼ 1.2K/GPa respectively, see below).
An external hidrostatic pressure P produces changes δb/b = δa/a = δc/c = −P/E on the lattice parameters. This modifies the Néel temperature for the C-AFM order through the exchange coupling constants J i dependence on the lattice parameters. For our model parameters (Table II), we find that the changes in T N are dominated by the dependence of J 1 on the lattice parameters For an isotropic Young modulus of 131GPa (see appendix C), we obtain δT N P ∼ 0.8K/GPa.
Ehrenfest's equations allow to relate the specific heat and the lattices changes on a second order transition For GdCoIn 5 , using V m = 155.85Å 3 , the lattice parameters of Table III, and the changes of specific heat and dilation reported in [9], we obtain ∆α V ∼ 6 × 10 −6 /K, ∆c P ∼ 15 J molK ,