Theory of magnetic ordering in the heavy rare earths: ab-initio electronic origin of pair- and four- spin interactions

We describe an ab-initio disordered local moment theory for long period magnetic phases and investigate the temperature and magnetic field dependence of the magnetic states in the heavy rare earth elements (HRE), namely paramagnetic, conical and helical anti-ferromagnetic(HAFM), fan and ferromagnetic (FM) states. We obtain a generic HRE magnetic phase diagram which is consequent on the response of the common HRE valence electronic structure to f-electron magnetic moment ordering. The theory directly links the first order HAFM-FM transition to the loss of Fermi surface nesting as well as providing a template for analysing the other phases and exposing where f-electron correlation effects are particularly intricate. Gadolinium, for a range of hexagonal, close-packed lattice constants, c and a, is the prototype and applications to other HREs are made straightforwardly by scaling the pair and quartic local moment interactions with de Gennes factors and choosing appropriate lanthanide contracted c and $a$ values.

for Gd's half-filled shell through to Lu's complete set of fourteen which causes the lanthanide contraction of the lattice [2]. The magnetism is complicated. When cooled through T c , Gd's paramagnetic (PM) phase undergoes a second order transition to a ferromagnetic (FM) state whereas, at T N , Tb, Dy and Ho form incommensurate, helical antiferromagnetic (HAFM) phases where the magnetization spirals around the crystal c-axis. When the temperature is lowered further both Tb and Dy undergo a first order transition at T t to a FM phase with basal plane orientation and Ho forms a conical HAFM ground state.
Further exotic phases emerge when the metals are subjected to magnetic fields and they have been extensively studied in experiments [3][4][5][6][7][8][9][10]. Below its T c Gd preserves its FM order as the strength of the magnetic field applied along its easy axis is increased. This is in sharp contrast to Tb [3,[11][12][13], Dy [4][5][6][7]14] and Ho [8,10] above T t , which first distort their HAFM order (dis-HAFM) before undergoing a first order transition into a fan magnetic structure followed by a second order transition to a FM state with further increase in the magnetic field. Dy and Ho also exhibit signs of additional spin-flip and vortex transitions associated with subtle changes in measured magnetization curves. In this letter we argue that much of this diversity stems directly from the valence electronic structure that all the HRE elements share.
Extensive experimental and theoretical investigations satisfactorily explain the onset of magnetic order from the PM state [15][16][17][18][19][20][21] via a detailed version of the famous Ruderman-Kittel-Kasuya-Yoshida (RKKY) pairwise interaction. The existence of nesting vectors q nest separating parallel Fermi surface (FS) sheets of the valence electrons provokes a singularity in the conduction-electron susceptibility. This feature results in a q nest -modulated magnetic phase [1,22], identified as a HAFM structure, incommensurate with the underlying lattice.
The lanthanide contraction changes the FS topology [23] and acts as the decisive factor for the emergence of the nesting vectors. This has resulted in the construction of a universal crystallomagnetic phase diagram which links the magnetic ordering that emerges from the PM state to the specific c and a lattice parameters of a heavy lanthanide system [16,18].
The prominence of RKKY interactions in the discussion of lanthanide magnetism promotes a deeper inspection of the common HRE valence electronic structure. As magnetic order among the local f-electron moments of the HREs develops with decreasing temperature and/or strengthening applied magnetic field the valence electron glue spin-polarizes and qualitatively changes. An indication of this was found by Khmelevskyi  lated effective exchange interactions from the FM state to be different from those in the PM state [24]. This effect has a potentially profound feedback on the interactions between the magnetic moments and has wider relevance for other magnetic systems where the physics is also typically couched in RKKY terms, including giant magnetoresistive nanostructures [25], rare earth clusters [26], magnetic semiconductors [27] and spin glasses [28]. In this letter we show how the response and feedback from the heavy lanthanide valence electrons to the ordering of local magnetic moments create multi-site interactions and determine the main features of the magnetic phase diagrams. We establish a reference, summarised in Fig. 1, against which these magnetic properties can be analysed to discriminate specific, subtle f-electron features.
A simple classical spin model with pairwise exchange interactions, magnetic anisotropy contribution, and a Zeeman external magnetic field term describes magnetic field-driven phase transitions in some anti-ferromagnetic insulators [29,30]. If the local moments of a HAFM state are pinned by anisotropy and crystal field effects to spiral around a particular direction, the effect of a magnetic field causes a first order transition to a fan or conical phase where the spins now oscillate about the field direction. In higher fields, the fan or cone angles smoothly decrease to zero in the FM state. Such a model applied to the HREs addresses only part of phenomenon, however, since it fails to reproduce the first order dis-HAFM to FM and second order fan to FM transitions at low temperatures and fields. It misses a tricritical point in consequence. In the seminal work by Jensen and Mackintosh [31] where the formation of field-induced fan and helifan phases was investigated theoretically for the first time, the key aspects of the HRE magnetic phase diagrams were only reproduced if ad hoc temperature dependent pair-wise exchange interactions were incorporated from a fit to spin wave measurements conducted at a series of temperatures. We find instead that much of the magnetic phase complexity is directly traced back to the behavior of the valence electrons. Evenson and Liu [15] maintained that the first order HAFM-FM transition is driven by a magnetoelastic effect. We show rather that while there is a magnetostructural coupling it is not necessary for the transition.
We have extended the ab-initio density functional theory (DFT) based, disordered local moment (DLM) approach [32] to address this issue and construct a generic H-T magnetic phase diagram of the HREs. Gd is a convenient prototype system owing to its seven localized f-electrons per atom in an S-state which form a large moment and the small crystal field and spin-orbit coupling effects that are prevalent. Choosing Gd enables us to abstract the common HRE valence electron effects on the magnetic properties and selecting the c and a lattice constants for other elements makes the analysis appropriate to Tb, Dy, Ho etc. Fig.1 shows the results for Gd using the lattice parameters appropriate to Dy [18].
The DLM-DFT describes the effects of thermally induced 'local moment' fluctuations on the underlying valence electronic structure of a magnet. For many materials such as the HRE's these magnetic excitations can be modeled by allowing the orientations of local, in the case of the HRE's f-electron, moments to vary very slowly on the time-scale of the valence electronic motions. By taking appropriate ensemble averages over their orientational configurations DLM-DFT determines the system's magnetic properties and describes magnetic phase diagrams ab-initio [33][34][35], temperature dependent magnetic anisotropy [36][37][38]  and field and temperature induced metamagnetic transitions [39].
The theory has the advantage that valence electronic structure can be monitored as a function of local moment disorder. This is highly pertinent owing to the recent development of advanced time-dependent spectroscopy techniques. Time-resolved resonant X-ray and ultrafast magneto-optical Kerr studies confirm the central tenet that the dynamics of the HRE core-like f-electrons are on a much longer time-scale than the excitations of the valence electrons [40]. Time-and angle-resolved photoemission (ARPES) studies have demonstrated the differing dynamics of spin-polarized valence states in correlated materials [41][42][43].  [19] angle-resolved photoemission measurements pointing to the magnetic exchange splitting of the FS as the principal mechanism for the fading of the nesting vectors [22] resulting in the stability of the FM phase at the ground state in Tb and Dy.
To follow the repercussions of this insight we specify a generalised Grand Potential Ω({ê n,i }) from DFT in which the local moments are constrained to point along directions . This quantity is averaged over many such configurations with a probability distri- . n and i count over layers stacked along the c-axis (i.e. the z-axis) and sites within a layer respectively. The local aver- (1) J nn 's are interpreted as pair-wise local moment interactions and the quartic coefficients K nn ,n n arise from the mutual feedback between local moment magnetic order and the spin polarized valence electrons illustrated in Fig.2.
For a phase diagram such as Fig.1 we construct the Gibbs free energy G of the system where T and S n (m n ) = −k B P n (ê) ln P n (ê)dê are the temperature and the magnetic entropy of the n-th layer, respectively (k B being Boltzmann's constant). The second term couples the external magnetic field H to the local magnetic moments each with magnitude µ. The last term, F u (m n ) = F 0 (ê n,i ·ẑ) 2 [45], describes a uniaxial anisotropy with strength F 0 [46] and fixes the easy axis. For selected T and H values, the Gibbs free energy where β = 1/k B T . The h n 's are therefore the Weiss fields for this mean field theory [32,34,35,47].
We carried out DLM-DFT calculations [44] for Gd within this framework for the c and a lattice parameters appropriate to Gd, Tb, Dy and Ho and in each case calculated charge and magnetization densities self-consistently for the PM state ({m n = 0}). The self interaction correction was used to capture the strong correlations of the f-electrons [18,48,49]. A local moment of µ ≈7.3 µ B established on each Gd site. We then divided the hexagonal lattice into 10 layer stacks and specified identical sets of {m n } (n = 1, · · · , 10) values for each stack to define the magnetic order parameter for each layer. For each c and a pair, using the effective one electron PM potentials, we calculated the h n = −∇ mnΩ values for each {m n } set. By thoroughly sampling the extensive the m n space we tested and established a method [44] to extract the J nn and K nn ,n n constants of Eq.1. We also checked that higher order terms were vanishingly small [44]. Fig.1 summarises the results. At a value m F M =0.503 the system undergoes a first order transition from a HAFM to FM state in zero field at T t =262K which correlates with the FS topological changes depicted in Fig.2. When the 4 site K nn ,n n 's are neglected in Eq.1 the calculated phase diagram is very different [44] and a FM phase does not appear at low temperatures and fields.
In the presence of long-ranged magnetic order, quantified by m FM , effective pair interactions mediated by the valence electrons affected by the long-range magnetic order can be specified as J ef f. nn = J nn + n n K nn ,n n m 2 FM [44] and they incorporate the influence of the 4-site terms from Eq.1. Fig.3 shows their lattice Fourier transform relevant to Figs.1 and 2 revealing the effect of the valence electron spin polarization. As shown in the inset, for m FM = 0, the interactions have a long-ranged oscillatory nature so that J ef f. (q) peaks at q nest ≈ 0.2 2π cĉ , (full red line). This is a direct consequence of the FS nesting shown in Fig.2(a) and which drives the HAFM magnetic order. We also show J ef f.  By comparing Gibbs free energies of the FM, HAFM, conical, fan, and helifan structures obtained, we constructed the T -H phase diagram. Fig.1 shows the results using c and a values appropriate to Dy when H was applied along the easy direction and continuous/dashed lines correspond to second/first order phase transitions. We imposed a single site uniaxial anisotropy of typical magnitude F 0 =+6.3meV/site [46] which precluded the conical phase when the magnetic field was applied in the easy ab-plane [44].   in Fig. 1). Theoretical estimate of the tricritical point ('A') is also given. Gd has a PM-to-FM second order transition at T C =274K (T C =293K in experiment [1]).
cooling and we find a tricritical point which is marked 'A' in Fig. 1.
We can adapt this Gd-prototype model to a specific heavy lanthanide element or alloy by using suitable lanthanide-contracted lattice constant values [18] and accounting for the specific f-electron configuration. The Gd ion has orbital angular momentum L = 0 and negligible spin orbit coupling effects. LS-coupling, however, is important for the HREs. A simple measure to account for the different total angular momentum values, J, is to scale the J nn and K nn ,n n interactions with the famous de Gennes factor [1] (g J − 1)2J(J + 1), where g J is the Lande g-factor. We applied this treatment to Gd, Tb, Dy and Ho and each metal except Gd had a phase diagram of the form of Fig. 1 consistent with experiment.
In Table I we compare results with those available from experiment for T N and T t , and the values of H for the highest T for the dis-HAFM phase. We also give a theoretical estimate of the tricritical point. The comparison overall shows that the theory correctly captures trends and transition temperature and field magnitudes. When the c and a values are further decreased our model predicts that the FM phase does not appear at low T and fields in accord with some experiments [52]. Discrepancies between the model and experiment can further highlight where f-electron correlation [53] effects are leading to more complicated physics. For example the complex spin slip phases in Ho reported at low temperatures [8,10] are not found in our model. The same applies to the vortex and helifan phases inferred from some experimental studies [7,31].
In summary we claim a dominant role for the valence electrons in the temperaturefield magnetic phase diagrams of the HRE metals. Our ab-initio theory incorporates lattice structural effects coming from the lanthanide contraction on this glue and makes the link be- describes the magnetic order within the n-th ferromagnetic layer. We express the free energy of the system as F =Ω − T k B n S n (m n ), where k B is the Boltzmann constant and the entropy per site in the n-th ferromagnetic layer is calculated as S n = 1 + ln(4π) + ln (sinh(A n )/A n ) − A n coth(A n ). The Gibbs free energy G is then obtained by adding the external magnetic field Legendre transformation to the free energy F and a term to describe a uniaxial anisotropy, which leads to Eq. (2) in the manuscript. An analytical expression for the uniaxial anisotropy can be derived by performing the integral F u (m n ) = F 0 (ê n,i ·ẑ) 2 = F 0 P n (ê n,i )(ê n,i ·ẑ) 2 dê n,i which becomes where F 0 is the strength of the anisotropy term. Prior to minimization of G we need to evaluate the dependence of the internal energyΩ({m n }) on the order parameters {m n } (or {A n }). Our DLM-DFT theory can accomplish this as it explicitly calculates ∇ mnΩ [32].
This requires an extensive exploration of the directional derivative in order to obtain an accurate description ofΩ. In the following section we explain in detail this methodology.

Minimization of the Gibbs Free Energy
Our methodology to minimize G consists in solving Eq. The magnetic structure obtained after this iterative procedure minimizes G, but it does not necessarily correspond to its lowest minimum. Therefore, the construction of the magnetic phase diagram requires a comparison of the Gibbs free energy of the various magnetic structures of interest. The magnetic structure with the lowest value of G is considered as the most stable phase for each (β,H) point. In particular, in the manuscript we have considered ferromagnetic (FM), and helical antiferromagnetic (HAFM), fan, helifan, and conical (CON) structures with a periodicity prescribed by 10 layers. The appropriate initial {A 0 n } arrangement needs to be chosen carefully for the exploration of each magnetic phase.
It is important to make the following point. Our DLM-DFT theory permits a direct calculation of ∇ mnΩ at each step of the iterative procedure. However, for each magnetic phase the calculation of G requires an evaluation of the internal energyΩ also. We achieve this by finding an analytical expression forΩ, i.e., Eq. (1) in the manuscript, that fits satisfactorily our numerical DLM-DFT data of ∇ mnΩ . The quadratic and quartic J nn and K nn ,n n coefficients are extracted as the result of this fitting. In the following section we give detailed information about the form ofΩ and ∇ mnΩ .

Analytical expression forΩ
The magnetic structures formed by the hcp HRE elements are composed of layers with FM alignment within each layer stacked perpendicular to the c-axis. The HAFM structure observed in Tb, Dy, and Ho is modulated along the c-axis and incommensurate with the lattice. The corresponding wave vector associated with the HAFM structure observed in our results (see Fig. 3 in the manuscript) is close to q nest 0.2 2π cĉ , which corresponds to turn angles of about 36 • in good agreement with experiment [19]. To study the period of the HAFM phase we have divided the hcp lattice system into 10 layer stacks. The magnetic order of each of these ferromagnetic layers is specified by the order parameters {m n }. Consequently, our internal energyΩ only depends on the 10 different vectors {m n , n = 1, .., 10} that are repeated along the c-axis. Taking into account the symmetry of the lattice system, we have found the following expression to fit satisfactorily the DLM-DFT data: Hence, we have structured the coefficients J nn and K nn ,n n in the manuscript by the 6 constants {J n , n = 1, .., 6} and the 7 constants {K n , K , n = 1, .., 6} respectively. To corroborate the accuracy of the J nn constants we have additionally applied a linear response theory at the PM state and recalculated them in the reciprocal space [18, 32? ]. We verified that the PM state is unstable to the formation of a HAFM structure prescribed by a 10 layer periodicity. The long-range constants (n − n > 6) shown in the inset of Fig. 3 in the manuscript have been calculated within this approach.
We have obtained DLM-DFT data for the Gd prototype within this framework. As explained in the manuscript, we have carried out calculations for Gd with the appropriate lattice spacings to mimic Tb, Dy, and Ho [18]. = − q,nn J nn + q ,n n K nn ,n n e −iq ·(R n −R n ) (m(q ) · m(−q )) e −iq·(Rn−R n ) (m(q) · m(−q)) , which directly defines some effective pair interactions in the presence of long-range magnetic order m(q ), i.e., shows mode-mode coupling, When q = 0, i.e. m(0) = m FM , J ef f. nn = J nn + n n K nn ,n n m 2 FM as mentioned in the manuscript. In the language of Eq. (6) these effective quadratic coefficients are expressed as The modified lattice Fourier transforms shown in Fig. 3 in the manuscript have been calculated from the effective pair interactions shown in the expression above.   obtained from the fitting of the DLM-DFT data of Gd as a magnetic prototype of Tb, Dy, and Ho.

Magnetic phase diagrams of Tb, Dy, and Ho
We show in Table I the results of the fitting of the DLM-DFT calculation of Gd with the lattice attributes of Tb, Dy, and Ho. Fig.4 shows the corresponding magnetic phase diagrams constructed from our results. The de Gennes factor [1] has been used to take into account the specific f-electron configurations. We remark that for the quartic coefficients we have used the square of the de Gennes factor. As stated in the manuscript the common features shown in Fig.4 are determined by the mutual feedback between the magnetic ordering and the valence electrons. An inspection of Fig.4 confirms that the lanthanide contraction favors the stability of the HAFM structure. In fact, for the lattice spacing of Ho, where the contraction is greater, the FM-HAFM transition disappears.
In Fig. 5 we show the magnetic phase diagram to accompany Fig. 1 in the manuscript (Dy lattice spacing) when the magnetic field is applied along the hard direction. The figure shows that the cone structure is stabilized if the magnetic field is applied parallel to the  hard axis. The transition from the cone to the FM/PM state is of second (first) order at high (low) temperatures. Fig. 1 in the manuscript is radically altered and qualitatively at odds with experimental results [7] if we neglect the quartic terms, i.e. set K nn ,n n =0, and so omit the feedback between the valence electronic structure and lanthanide local f-electron magnetic moment order. This is shown in Fig. 6.