A new theoretical approach to the strain dependence of magnetic crystal-field anisotropy

We report on the derivation of analytical equations for ab-initio calculations of the strain dependence of crystal-electric-field (CEF) parameters for arbitrary deformations. The calculation is based on the fundamental assumption that the charge distribution deforms in the same way as the crystal. Based on this deformed-charge model, simple formulas for the practical usage are given for various site symmetries of cubic lattices under uniform strain. These formulas can be used to predict the change of the magnetic crystal-field anisotropy under strain, which is important for the design of magnetic materials and devices. As an example for the power of the method, we present a calculation of the magnetic contribution to the thermal expansion in some rare-earth-based materials.


Introduction
The strain dependence of magnetic anisotropy is of great importance in the design of novel materials and needs to be tuned to achieve technological applications of magnetic materials. For example, advanced ultrathin spintronic nanodevices with customized electronic and magnetic properties will need a local-strain engineering [1]. Historically, this problem has been recognized in investigations of magnetoelastic parameters in rareearth materials [2][3][4]. Some recent research with focus on transition-metal oxides has been using cluster methods [5] and spin-polarized bandstructure calculations [6]. However, a full understanding and quantitative prediction of the strain dependence of magnetic interactions is missing.
Research on this matter started with the theoretical analysis of magnetoelastic constants [2][3][4]. Since then, the huge problem of calculating reliable crystal-electricfield (CEF) and exchange parameters in magnetic materials by ab-initio methods have been in the focus of research for many years.
In many studies of magnetoelastic effects, the comparison of theory and experiment involved the fitting of one or more magnetoelastic parameters to the data [7,8].
In this work, we make an attempt to calculate the strain dependence of the magnetic crystal-field anisotropy for an arbitrary deformation. In many cases, these single-ion effects play an essential role in the description of the anisotropic elastic behava e-mail: thomas.stoeter@tu-dresden.de (corresponding author) ior. Other anisotropies, such as those of the magnetic exchange (e.g., dipolar or Dzyaloshinskii-Moriya interactions, higher-order coupling effects), are not considered here. However, these can be additionally included using the analysis of the distortion dependence of the exchange parameters. Thus, it is in principle possible to investigate changes in physical properties caused by a symmetry breaking, such as elimination of latticeinduced magnetic frustrations. Macroscopic magnetic anisotropies, such as shape or surface anisotropy, are outside the scope of our microscopic considerations. For example, the considered magnetic deformations can be induced in magnetic devices by (i) external tensile or (ii) compressive strain or by (iii) chemical pressure. Experimentally, the effect on the anisotropy can be measured by high-sensitivity magnetization and other bulk techniques or spectroscopic methods [9]. Moreover, we show how the strain dependence of magnetic anisotropy may be used to predict the magnetic contribution to the thermal expansion.
This paper is made up as follows: In Sect. 2, we give a thorough introduction to the formulas important for the magnetoelastic coupling. We use different representations of the strain tensor to derive explicit formulas of the thermal expansion for various strain modes. In Sect. 3, we give the main result of this paper: the strain-derivative of a quantity closely related to the CEF parameters. From this quantity, the change of the CEF parameters under strain may easily be calculated. We evaluate the main formula for all possible strain modes in Sect. 4. Finally, we show the calculated CEF volume striction of several compounds in Sect. 5. We 125 Page 2 of 9 Eur. Phys. J. B (2021) 94 :125 close with an outlook and a summary of the results in Sects. 6 and 7.

A recapitulation of the theory
The strain is a 3×3 matrix, a second-rank tensor, that represents the deformation of a crystal. In the following, we consider only homogeneous strain, i.e., the strain is constant in all of the samples. We use two different ways to represent the strain tensor. On the one hand, we decompose in irreducible representations of the space group of the crystal and we denote the components of the strain tensor as Γ j i . Here, Γ denotes the irreducible representation of the space group of the crystal, the index j differentiates between multiple representations of the same type Γ , and the index i denotes the different components of the (multi-dimensional) irreducible representation. In this decomposition, the elastic energy takes on a simple form. On the other hand, we decompose in terms of irreducible tensors, i.e., irreducible representations of the rotation group, that transform into each other via a simple law under rotations. We denote the components with JM , with J and M denoting the irreducible representation and its 2J +1 components, respectively; these are analogous to the angular momentum formalism in quantum mechanics. Since the strain tensor is symmetric, both these representations have six components.
The simplest possible expression of the harmonic elastic energy of a unit cell is given by [2,3,10]: where mn are the Cartesian components of the strain tensor , c mnkl are the Cartesian components of the elasticity tensor characteristic of each compound and obey certain symmetries depending on the lattice; V is the volume of a unit cell of the crystal. Depending on the lattice symmetry, we can rewrite the elastic energy in terms of strains adapted to the lattice symmetry Γ,j i [2]: The c Γ jj are the elastic constants of the symmetryadapted strain modes Γ,j i of the ith component of the jth one-, two-or three-dimensional irreducible representation Γ ; these strain modes are given for each lattice symmetry in the classic paper [2]. In particular, the elastic energies of the cubic and hexagonal crystal with the conventional superscripts are given by [3]: Besides the elastic energy, we use the CEF energy to investigate magnetoelastic coupling in this work. Below the ordering temperature of a magnetic compound (if existing), exchange effects are dominating the magnetoelastic behavior of a compound. However, also in that region the CEF-striction mechanism is contributing; however, only above the ordering temperature and at zero field, the CEF striction can be isolated from exchange striction. The CEF HamiltonianĤ CEF of one unit cell projected on the J-manifold of the ground state of the free rare-earth ion up to first order in the symmetrized strains Γ,j i is given by: where B l m (f, f ( )) andÔ l m (J f ) are Stevens parameters and operators of the f th ion in the unit cell under the strain f in the local coordinate system at that position, respectively. The operatorsÔ l m (J f ) are represented in the angular-momentum operators of the local coordinate system at the f th ion. The first term on the right is the CEF Hamiltonian of the undeformed crystal and the second term is the magnetoelastic coupling Hamiltonian. The coupling constants are, therefore, ∂ Γ,j i B l m (f, f ( )), the strain dependence of the CEF parameters, the quantity we are interested in.
The components of f are related to the components of the global strain via rotation matrices. The latter is the strain in the global reference frame parallel to the symmetry axes of the crystal and the former is represented in the local coordinate system used for the CEF Hamiltonian. To get simple transformation laws between the components of the strain tensor, we avoid using the Cartesian standard basis and instead use irreducible tensors as basis. The components in this basis are denoted with f JM and they relate to the components JM in the global coordinate system via: where D J M ,M (α f , β f , γ f ) denotes the Wigner-D-matrix of the rotation of Euler angles α f , β f and γ f between the local and global coordinate system [11]. The strain tensor, being a symmetric tensor, is made up of the six components with J = 0, M = 0 and J = 2, M = −2, −1, . . . , 2. The irreducible tensors in Cartesian coordinates are given in the appendix. Hence, the derivative of the CEF parameter B l m at ion f can be written as: where α f , β f and γ f are the Euler angles that rotate the global reference frame to the local reference frame at ion f . For convenience in the following formulas, we use for all angles, which leads to simple evaluation of the uniform strain mode.
For high symmetry lattices such as cubic, hexagonal, trigonal and tetragonal lattices, the symmetrized strains of the lattice coincide with the irreducible tensors, i.e., the last factor in Eq. (12) is one or zero. For the other lattices, simple linear relations between the symmetrized strains of the lattice and the irreducible tensor components can be found by decomposing the first in terms of the latter, as shown in the Appendix A.
Before calculating the strain dependence of the CEF parameters, we show how it can be used to calculate the equilibrium strain. The elastic energy and the CEF energy change with deformation and compete with each other. In equilibrium, the free energy is minimized, i.e., the derivative of the free energy with respect to Γ j i vanishes [2,3,10]. The strain-dependent part of the free energy is simply the sum of the CEF energy and the elastic energy: where Ô l m denotes the thermal average of the expectation value of the operatorÔ l m . The equilibrium condition yields a system of linear equations in Γ,j i [2]. For example, for the cubic lattice (omitting the index j = 1): and for the hexagonal lattice with Δ = c α 11 c α 22 − c α 12 c α 12 : Having given the explicit formulas of the equilibrium strain, we now come to the evaluation of the strain derivative of the CEF parameters, which will be given in the next section hosting the main result of this paper.

Strain derivative of the CEF parameters
The CEF parameters B l m (f ) are linear functions of ω l m (f ), which are expressed by: where ρ f is the charge distribution in the crystal centered on the f th ion in the unit cell in the local frame of reference, Y l m denotes the spherical harmonics, Ω R the solid angle and R the position vector. The B l m are defined by the tesseral (i.e., real spherical) harmonics and related to ω l m via: where β l m = −|e|p l m r l 4f Θ l . In order to go transform back and forth between B l m and ω l m , we also need the inverse relation: The main assumption of our model is that the charge distribution deforms in the same way as the lattice. This assumption seems natural and is usually justified a posteriori. In the vicinity of structural phase transitions, this assumption may be violated, for example near the variant conversion transition in magnetic shape memory systems such as RCu 2 [12]. For the general formulas, we derive here, we parametrize the specific strain modes by the linear expansion = Δl/l such that the tensor = M , where M is a constant 3 × 3 matrix. To take full advantage of the formalism of spherical harmonics in [11], we will represent the components of the tensor in the spherical basis [11] and not Cartesian; in this basis, we use the indices σ and σ for and the matrix M. Suppose the lattice deformation in the local reference frame is given by the matrix O = E + f M with E the 3 × 3 unit matrix and M a symmetric matrix. Then, we assume that the deformed charge distribution ρ f is related to the undeformed charge distribution ρ f with the deformation matrix O −1 via: where the factor det(O −1 ) is needed to ensure the invariance of the total charge: where to ω l m but of a higher multipolar order of maximal order 8.
In general, not all necessary quantities can be obtained from the unperturbed system. Especially, the last term makes evaluation for small arbitrary deformations difficult. For the special case of the volume changing mode f 00 , however, evaluation is simple. Yet, it is only in systems with cubic lattice symmetry, that the volume changing mode does not couple to other modes when considering the equilibrium strain. Therefore, even though we give explicit formulas for all modes f JM , the calculated example afterwards is focused on the volume changing mode f 00 in cubic systems.

Explicit formulas
In this section, we evaluate the previous formula for the different JM with their specific M in spherical coordinates using the Kronecker delta function δ nm . We start with the simplest deformation mode: where the formula follows from straightforward evaluation or the identity αβ C cγ aαbβ C c γ aαbβ = δ cc δ γγ [11], since the second term in Eq.
From the linearization of ω l m , we can obtain the linearization of B l m via Eq. (24), for example, for the mode f 00 in cubic crystals: Depending on the CEF symmetry at the relevant rareearth site, only few Stevens parameters are non-zero. For the mode α , we now also give the equilibrium strain, which is particularly simple to express, because 3 00 independent of the local coordinate system at position f . Hence, the equilibrium strain in an arbitrary CEF symmetry and for a cubic Bravais lattice reads: where N is the number of rare-earth ions in the unit cell and B l m (1, 0) and Ô l m (1) are evaluated for the first ion as representation of all others in the unit cell. Thus, we see that for the calculation of volume expansion in cubic crystals, the knowledge of the CEF parameters is sufficient. However, the Stevens parameters have very complicated influence on the expansion because they are not only explicitly used in the formulas but also implicitly contained in the CEF levels and wave functions used for thermal average Ô l m . The formulas of the Hamiltonians and expansions for the several important point symmetries are given in Table 1.
We want to mention the effects of anisotropic deformations that are not important for the case of zero field and the paramagnetic state in cubic crystals; following [2] for symmetry reasons the anisotropic modes vanish. However, if behavior in field, axial stress (as in recently available strain cells) or the thermal expansion of crystals with lower symmetry than cubic are considered also the anisotropic modes 2M with M = 0, ±1, ±2 will become relevant and the according expressions above may be used. For example, the hexagonal, trigonal and tetragonal Bravais lattices have one more symmetry conserving mode γ,1 that decomposes in irreducible tensors of type 2M with M = 0, ±1, ±2.
Coming back to the case of the thermal CEF volume striction of cubic crystals, we use Eq. (45) in the examples in the next section.

Example: CEF parameters and CEF volume striction of pyrochlore titanates
In this section, to illustrate the theoretical part, we show how to apply the derived equations from the last section, to calculate the CEF volume striction given only the CEF parameters, the elastic constants and the lattice parameters. The thermal averages of the expectation values of the Stevens operators Ô l m are calculated with McPhase [13]. We show how to derive the CEF parameters for members of an isostructural series when the CEF parameters of one of them is known and the lattice parameters of the whole series is known. The calculated CEF is then used to calculate the CEF volume striction for all compounds of the isostructural series. We use the isostructural series of rare-earth titanate pyrochlores as an example, crystallizing in the cubic space group F d3m (Nr. 227).
The CEF parameters of rare-earth compounds are calculated via [10]: where −|e|p l m is a factor independent of the rare-earth ion and the lattice, γ l m are the coefficients of the multipole expansion of the charge distribution ρ and a linear combination of the ω l m so that the same behavior under volume change applies that we have derived earlier; r l 4f Θ l is a factor dependent only on the ground state and the kind of the rare-earth ion. Θ l is α J , β J or γ J for l = 2, 4, 6, respectively, and tabulated, cf. [10]. If the substiitution of the rare-earth ion R by another rareearth ion R' in a compound does not change the lattice symmetry but only the lattice parameter, then there are two parts contributing to changes of the Stevens parameters: First, the linear expansion or contraction of the lattice influences the CEF parameters according to Eq. (41), with = √ 3 f 00 . And second, the substitution of the rare-earth ion has a large influence shown in Eq. (46). Putting both together, we get: In the following example, we use this formula to calculate the CEF parameters of compounds in the same isostructural series and using these we calculate the CEF volume striction. The rare-earth ion resides at Wyckoff position 16d with local CEF symmetry D 3d . The cubic crystal structure is stable for the relevant temperature range. Using Eq. (45) in the special case of a point symmetry D 3d (cf. Table 1), we can calculate the contribution of the CEF to the thermal expansion. For all rare-earth titanate pyrochlores, the CEF parameters of the rare-earth ions were calculated from the ones of Ho 2 Ti 2 O 7 [14] using Eq. (49). The values compare very well to the ones given in [14] and were then used to calculate the CEF volume striction shown in Fig. 1.
We calculate the CEF volume striction with McPhase using the lattice parameters from [15] and an estimate of the elastic parameter c α = 568 GPa measured for Ho 2 Ti 2 O 7 [16]. Figure 1 shows that the CEF has a significant contribution to the thermal expansion between ΔL/L = CEF = 4× 10 −4 for Tb 2 Ti 2 O 7 and 7.5 × 10 −4 for Tm 2 Ti 2 O 7 between 0 and 500 K. Compared to phononic thermal expansion in the order of 3 × 10 −3 , this would be a measurable effect. Essentially, the order of magnitude is the same for all titanate pyrochlores and largest for Tm 2 Ti 2 O 7 .

Outlook
Several possible ways to work with the results presented in this paper come to mind.
First, for high-symmetry crystals, an intensive study of Eq. (28) might allow evaluation of arbitrary homogeneous deformations. The authors expect many of these values to vanish for high-symmetry crystals. Second, an investigation of low-symmetry crystals seems promising, since significantly larger CEF striction effects can be expected from our model. The reason for this is that the elastic constants, especially shear elastic constants are considerably softened in non-centrosymmetric crystals [17]. Third, an extension of the theory to the strain dependence of the exchange parameters would also be possible. This would take a thorough investigation of the symmetry of the two quantum states of neighboring atoms and is a natural next step to consider. Important theoretical work has already been done in the review [18]. Finally, our formulas may be used as an equationof-state for numerical simulations: Strictly speaking, a crystal is an object that is invariant under some (discrete) translation group; inhomogeneous deformations would break this translation invariance, as would statistically distributed defects. However, given that the quantum state at a specific atom in the crystal is influenced only by a small region surrounding it, the local symmetry of that environment at the position of the atom needs to be considered. That symmetry group may change within the crystal to allow symmetry breaking by defects or non-homogeneous deformations, as possible in non-centrosymmetric crystals [19]. However, the solution of such a problem in general is not accessible with analytic methods and one needs to use numerical methods.

Summary
In this paper, we have shown new formulas of the effect of a deformation on the CEF parameters associated with Stevens operators. We started from a socalled deformed charge model, where the lattice deformation affects the charge distribution one-to-one and even charge clouds are deformed. This ansatz allowed us to derive an analytical formula of the change of the CEF parameters for small deformation . We applied this formula to different deformation modes and found in the simplest case of a pure volume changing mode a simple formula. The resulting formulas can be used to calculate the change of the CEF parameters in new strain experiments on bulk materials as well as thin films. We also gave a recapitulation of the general theory of CEF striction with explicit formulas for cubic and hexagonal crystals. With these formulas, it is easy to expand the theory given in this paper if the formulas of the strain derivative of the Stevens parameters can be evaluated. Finally, we used the formulas for several cases of current scientific interest, such as the rare earth pyrochlores including the spin-ice compounds Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 .
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Irreducible tensors in Cartesian coordinates
The irreducible tensors in Cartesian coordinates are given for the convenience of the reader and may be found in literature [11]: The N × N square matrices form an N 2 -dimensional vector space. Together with the Frobenius inner product as scalar product, we can project out components of matrices in the same way as components of vectors. The Frobenius inner product of two N × N matrices A and B is given by A · B = ijĀ ij Bij. With this scalar product, the irreducible tensors above form an orthonormal basis for the subspace of symmetric matrices like the strain tensor. We can use this basis to uniquely decompose arbitrary symmetric N × N matrices. To illustrate the procedure of decomposing a matrix into its irreducible components, we perform this on the symmetric matrix: (55) For the sake of completeness: To decompose arbitrary N × N matrices or antisymmetric tensors, which we do not need to do in this paper, the irreducible tensors 11 and 1−1 are needed, as well.

B Derivation of the crystal field parameters
Our goal is to expand ω l m ( ) linearly for small deformations of the charge distribution ρ with deformation matrix O = E + M.

(64)
So considering the product Y l+1 μ Y 1 σ there are three possible L = l, l + 1, l + 2 of which only L = l or L = l + 2 contribute since C L0 l+1010 is zero for odd L. As l is even, we get Inserting this into Eq. (61), we obtain Eq. (63), the main result of this paper. This formula is used in the main text to calculate the change of the CEF parameters during a deformation of the lattice.