Theoretical Background

Localised rare-earth and actinide ions are well described by $LS$ -coupling with a ground state multiplet determined by Hund's Rules, and characterised by the total angular momentum quantum number $J$ . In a free ion, these $J$ -multiplets containing $2J+1$ levels labeled by $m_J = -J,\ldots,+J$ are degenerate. However, crystal field interactions arising from the breaking of the symmetry of free ions by the regular arrangements of a crystal will break this degeneracy of the free ion $J$ multiplets. This occurs because whereas a free ion will have a Hamiltonian with spherical symmetry, the magnetic ions in a crystal will have a different site symmetry.

Thus the effect of a crystal field can be thought of as a perturbation acting on free magnetic ions due to the mechanisms of crystal bonding. Initially it was attributed to the electric potential due to the charge distribution of surrounding ions on a magnetic ion within the crystal, hence the term crystalline electric field. Expressing this contribution by a multipolar expansion centred on the magnetic ion yields the crystal field potential as a sum over spherical harmonics $Y_{k,q}$

$$V_{CF}(\mathbf{r}) = \sum_{k,q} B_q^k r^k Y_{k,q}(\theta,\phi) %\\$$ (1)

$B_q^k$ are crystal field parameters which determines the strength the interaction. The order $k$ runs over even integers up to $2l+1$ where $l$ is the orbital angular momentum number of the equivalent electrons in question. Thus for rare earth ions, it is 3, and $k=2,4,6$ . The component $q$ runs over $-k,-k+1,\ldots,k$ .

In order to calculate the Hamilton matrix for the crystal field interaction we have to express the CF potential as an operator $\hat{V}_{CF}$ whose matrix elements may be easily calculated. Then the eigenvalue problem $\hat{V}_{CF} \vert\psi\rangle = E \vert\psi\rangle$ may be solved to find the energies $E$ (to an arbitrary constant shift) and wavefunctions $\vert\psi\rangle$ of the magnetic ion in the crystal field. We now come to a problem of notation and normalisation. The first formalism to handle crystal fields was developed by Stevens using operator equivalents [Stevens(1952)]. Stevens first expressed the term $r^k Y_{k,q}$ in Cartesian coordinates, to obtain for example, $r^2 Y_{20}(\theta,\phi) = r^2 (3\cos^2\theta - 1) \equiv (3z^2 - r^2)$ , where we have ignored a constant term in the spherical harmonics. He then used the result that within a manifold of constant $J$ - such as the ground $LS$ state of a rare earth ion - the potential operator $r^k Y_{k,q}$ is equivalent to a similar angular momentum operator formed by taking the Cartesian expression of the potential operator, replacing terms in $x, y, z$ by terms in $J_x, J_y, J_z$ and symmetrising to allow for the non-commutation of $J_x, J_y, J_z$ . The two expressions are then related by an operator equivalent factor, $\alpha, \beta, \gamma$ for $k=2,4,6$ respectively. Thus, for the $q=0$ terms, for example we have:

$$r^2 Y_{20}(\theta,\phi) \equiv r^2 (3\cos^2\theta-1) = (3z^2 - r^2) \equiv \alpha \langle r^2 \rangle (3J_z^2 - J(J+1))$$

$$r^4 Y_{40}(\theta,\phi) \equiv r^4 (35\cos^4\theta-30\cos^2\theta+3) = (35z^4-30r^2z^2+3r^4)$$ (2)

$$\equiv \beta \langle r^4 \rangle (35J_z^4-30J(J+1)J_z^2+25J_z^2-6J(J+1)+3J^2(J+1)^2)$$

$$r^6 Y_{60}(\theta,\phi) \equiv r^6 (231\cos^6\theta-315\cos^4\theta+105\cos^2\theta-5) = (231z^6-316r^2z^4+105r^4z^2-5r^2)$$

$$\equiv \gamma \langle r^6 \rangle (231J_z^6-315J(J+1)J_z^4+735J_z^4+105J^2(J+1)^2J_z^2-525J(J+1)J_z^2$$

$$+294J_z^2-5J^3(J+1)^3+40J^2(J+1)^2-60J(J+1)$$

Stevens then listed the matrix elements for these three operators, and later he and others added matrix elements of all terms with even $q$  [Jones et al.(1959)Jones, Baker, and Pope], which enables magnetic ions in sites of orthorhombic ($C_{2v}$ space group) or higher symmetry to be calculated. The expressions in the angular momentum operators given in equations 2 are termed Steven's operators $O_q^{(k)}$ by later authors, and are listed for all $q$ and $k=0,\ldots,6$ by Smith and Thornley [Smith and Thornley(1966)]. The CF parameters are now specified by $A_{kq}$ , so that the crystal field operator is:

$$\hat{V}_{CF} = \sum_{k,q} A_{kq} \langle r^k \rangle \theta_k O_q^{(k)}$$ (3)

where the operator equivalent factors $\theta_k$ is $\alpha, \beta, \gamma$ respectively for $k=2,4,6$ , after the notation of Judd. These factors were listed for the ground states of the rare earths by Stevens, and a method to calculate the values for other $LS$ states was given by Elliot et al. [Elliot et al.(1957)Elliot, Judd, and Runciman]. Although both $A_{kq}$ and $\langle r^k \rangle$ may be calculated separately from models, we choose to regard both as a single parameter to be varied, and henceforth will refer to $A_{kq}\langle r^k \rangle$ as crystal field parameters under Steven's normalisation. This is because the main use of the SAFiCF program is to fit crystal field parameters to an phenomenological model rather than to compare with theoretical models of the crystal field.

However Steven's methods, whilst effective for hand calculations is not efficient for machine calculations as it would require inputing tables of matrix elements (which may be prone to error, and be hard to debug) and looking up their values. Instead we shall follow the later methods of Wybourne, Judd and others, using the tensor operator techniques due first to Racah [Judd(1998)]. We express the crystal field potential operator in terms of a tensor operator $C_q^{(k)}$ which transforms in the same way as the quantities:

$$C_q^{(k)}(\theta,\phi) = \sqrt{\frac{4\pi}{2k+1}} Y_{k,q}(\theta,\phi)$$ (4)

This defines the normalisation of the crystal field parameters, and is called the Wybourne normalisation by Newman and Ng [Newman and Ng(2000)] ¹. The crystal field operator is now expressed as:

$$\hat{V}_{CF} = \sum_{k,q<0} i B_q^k \left[ C_{\vert q\vert}^{k} -... ...B_0^k O_0^{(k)} + \sum_{k,q>0} B_q^k \left[ C_q^{k} + (-1)^q C_{-q}^{k} \right]$$ (5)

such that $\hat{V}_{CF}$ is hermitian and the parameters $B_q^k$ are real. The matrix elements of $C_q^{(k)}$ may be calculated first by factorising into a part dependent only on $m_J$ and $q$ , and reduced matrix element dependent on $J$ and $k$ by application of the Wigner-Eckart theorem:

$$\langle J,m_J\vert C_q^{(k)} \vert J,m'_J\rangle = (-1)^(J-m_J) \... ..._J & q & m_J \end{array} \right) \langle J\vert\vert C^{(k)}\vert\vert J\rangle$$ (6)

The reduced matrix element $\langle J\vert\vert C^{(k)}\vert\vert J\rangle$ is just a number and is given by:

$$\langle J\vert\vert C^{(k)}\vert\vert J\rangle = (-1)^J (2J+1) \sqrt{2k+1} \left( \begin{array}{ccc} J & k & J 0 & 0 & 0\end{array} \right)$$ (7)

where the large brackets denotes a Wigner $3j$ symbol, which is related to the Clebsch-Gordan coefficients resulting from the coupling of two angular momenta.

However, because the Steven's normalisation is still widely in use, SAFiCF defaults to using it, rather than the Wybourne normalisation. In this case, whilst the $q$ and $m_J$ dependent parts of the matrix elements are the same, we need to replace the reduced matrix element $\langle J\vert\vert C^{(k)}\vert\vert J\rangle$ by another reduced matrix element appropriate to the Steven's operators [Smith and Thornley(1966)]:

$$\langle J\vert\vert O^{(k)}\vert\vert J\rangle = \frac{1}{2^k} \sqrt{\frac{(2J+k+1)!}{(2J-k)!}}$$ (8)

This result may be derived by first taking the $q=0$ components, and constructing the diagonal matrix elements for the operators in equations 2. The Wigner-Eckart theorem then gives:

$$\langle J,m_J\vert O_0^{(k)}\vert J,m_J\rangle = (-1)^{J-m_J} \le... ...J & 0 & m_J \end{array} \right) \langle J\vert\vert O^{(k)}\vert\vert J\rangle$$

Taking $k=2$ , we find that $\langle J,m_J\vert 3J_z^2-J(J+1)\vert J,m_J\rangle = 3m_J^2-J(J+1)$ , and using the algebraic expression for the 3$j$ symbol in the above equation²:

$$\left( \begin{array}{ccc} J & 2 & J -m_J & 0 & m_J \end{array} \right) = (-1)^{J-m_J} \frac{2(3m_J^2 - J(J+1)}{[(2J+3)!/(2J-2)!]^\frac{1}{2}}$$

we arrive the result 8 with $k=2$ . Similar techniques for $k=4$ and $k=6$ yields the $k$ terms in equations 8.

These equations to calculate the crystal field Hamilton matrix within a constant $J$ -manifold is implemented in the function cf_hmltn(J,A2,A4,A6), which by default assumes that the CF parameters, expressed as $2k+1$ component vectors A2,A4,A6, are in Steven's normalisation³.

We now turn to the interaction of the neutron with the electrons of the magnetic ions in a crystal field. Approximating the neutron wavefunction as a plane wave and using Fermi's Golden Rule to obtain the transition probability, then multiplying this by the density of neutron final states and dividing by incident neutron flux (which quantities are proportional to $k'/k$ ), the double differential cross-section is:

$$\frac{d^2 \sigma}{d\Omega dE'} = \frac{k'}{k} \sum_{\sigma \sigma... ...{\kappa})\vert\sigma \mu\rangle \vert^2 \delta(\hbar \omega + E_\mu - E_{\mu'})$$ (9)

where $k$ and $k'$ are the wavevectors of the initial, $\mu$ , and final, $\mu'$ , states of the target ion. The neutron's initial and final polarisation is labeled by $\sigma$ and $\sigma'$ respectively, and $p_\sigma$ is the initial probability distribution of the neutron's spin, whilst $p_\mu$ is the probability distribution of initial target states. $\mathbf{\kappa}=\mathbf{k}-\mathbf{k}'$ is the wavevector transfer, and $\hat{V}(\mathbf{\kappa})$ is the neutron-atom interaction operator.

For magnetic scattering only (ignoring neutron-nuclear interactions), the neutron-electron interaction operator is:

$$\hat{V}_M(\mathbf{\kappa}) = r_0 \hat{\mathbf{\sigma}} \cdot \fra... ...appa^2} [ \mathbf{\kappa} \times \{\hat{\mathbf{Q}} \times \mathbf{\kappa} \} ]$$ (10)

Where $\hat{\mathbf{Q}}$ is the intermediate scattering operator, which is the Fourier transform of the magnetisation density, and has spin and (orbital) momentum dependent parts. In the limit $\kappa \rightarrow 0$ , however, we can use the dipole approximation, and obtain:

$$\hat{\mathbf{Q}}(\mathbf{\kappa}) = \frac{1}{2} \left( \hat{\mathbf{L}} + 2\hat{\mathbf{S}} \right)$$ (11)

SAFiCF does not handle polarised neutrons at present, but we hope to implement it in a future version, so for now, we will omit the sums over $\sigma$ . Furthermore, we take the ideal case where $\kappa = 0$ , and for Russel-Saunders ($LS$ ) coupling replace $\mathbf{L}+2\mathbf{S}$ by $g_J\mathbf{J}$ obtain:

$$\frac{d^2 \sigma}{d\Omega dE'} = \frac{k'}{k} \sum_{\mu \mu'} \fr... ...mathbf{J}\vert\mu\rangle \right\vert^2 \delta(\hbar \omega + E_\mu - E_{\mu'})$$

Where $Z = \sum_\mu \exp(-E_\mu/k_B T)$ is the partition function. We can now average over solid angle to obtain the cross-section at the energies of the transitions between electronic states $\omega = E_\mu-E_{\mu'}$ :

$$\sigma(\omega) = 4\pi \left(\frac{r_0 g_J}{2}\right)^2 \left\{ \f... ...\mu\rangle \vert^2 + \vert\langle \mu'\vert J_z\vert\mu\rangle \vert^2 \right\}$$ (12)

Alternatively, we can calculate the scattering function at $\kappa = 0$ :

$$S(\mathbf{0},\omega) = \frac{1}{N\sigma_b} \left(\frac{r_0 g_J}{2... ...\mu\rangle \vert^2 + \vert\langle \mu'\vert J_z\vert\mu\rangle \vert^2 \right\}$$ (13)

where $N$ is the number of scatterers and $\sigma_b$ is their scattering length. The two equations 12 and 13 are calculated by the function cflvls(Hcf,T) where Hcf is the crystal field Hamilton matrix calculated as described above by cf_hmltn. The matrix elements $\langle \mu'\vert J_\alpha\vert\mu\rangle$ are calculated from the matrix elements of $O_q^{(1)}$ as described above also, by $J_x = O_1^1 - O_1^{-1}$ , $J_y = i (O_1^1 + O_1^{-1})$ and $J_z = O_0^1$ by the function mag_op_J(J).

Finally, we note that this treatment of inelastic neutron scattering by magnetic ions in a crystal field is an approximation which may not agree with measurements. In particular, the magnetic form factor for real measurements with $\kappa \neq 0$ will change the intensity of spectra, and the neutron-electron operator will no longer be proportional to $\mathbf{J}$ . A more sophisticated technique is described by Balcar and Lovesey [Balcar and Lovesey(1989)], again using Racah tensor operator techniques, and will be implemented in future versions which will also include calculations of the crystal field levels of ions in intermediate coupling.