The total Crystal field Phonon Interaction Hamiltonian

To calculate crystal field parameters $B_{\gamma}(i)=B_l^m$ ($\gamma$ is a short hand notation for $lm$) for the ion at site $i$ in the crystal with Stevens convention we use the well known expressions of the point charge model:

  $$B_{\gamma}(i)=B_l^m=-\vert e\vert\theta_l^i\langle r^l \rangle \gamma_{lm}(i) p_{lm}$$ (96)

Here, the $\gamma_{lm}(i)$ are computed from the relative position $R_{ji},\Omega_{ji}$ (in spherical coordinates) of the point charges $q_j$ using

  $$\gamma_{lm}(i) = \sum_j \frac{q_j}{2l+1}\frac{Z_{lm}(\Omega_{ji})}{\epsilon_0R_{ji}^{l+1}}$$ (97)

Equation (104) is based on the adiabatic approximation and the dependence of the crystal field parameters on the displacements $\hat \mathbf U_i$ is not considered. Going beyond the adiabatic approximation this dependence may be considered and the crystal field parameters may be expanded in terms of the strain $\epsilon$ and the $\hat \mathbf P_i$ leading to the crystal field phonon interaction. Minimizing the Energy with respect to the strain tensor $\epsilon$ leads to expressions for the strain $\epsilon$ in terms of expectation values of Stevens Operators and displacement operators.

The Hamiltonian can be written as a sum of phonon, crystal field and Zeeman contributions (writing $\gamma$ instead of $lm$ for the crystal field parameter indices, the $i$ index counts nuclei carrying with them the charge producing the crystal field $i=1,...,N_{\rm nuclei}$, the $j$ index runs over magnetic ions,i.e. the only partially filled f or d electron shells in a crystal $j=N_{\rm nuclei}+1,... N_{\rm nuclei}+N_{\rm magnetic ions}$, this is why we writ $i<j$ in the following sums)

	$$\hat H_{\rm tot}$$	$=$	$$\hat H_{\rm phon} +\sum_{i,\gamma} B_{\gamma}(i,\hat \mathbf U_1,... ... O_{\gamma}(\hat \mathbf J_i) - \sum_i g_{J_i} \mu_B \hat \mathbf J_i \mathbf H$$	(98)
		$=$	$$\hat H_{\rm phon} +\sum_{i,\gamma} B_{\gamma}(i,0,\dots,0) O_{\gamma}(\hat \mathbf J_i) + \hat H_{\rm cfph}$$

	$$\hat H_{\rm cfph}$$	$=$	$$\sum_{i<j,\gamma}\nabla_{\hat \mathbf U_i} B_{\gamma}(j) \hat \mathbf U_i O_{\gamma}(\hat \mathbf J_j)$$
		$=$	$$\sum_{i<j,\gamma}\nabla_{\hat \mathbf U_i} B_{\gamma}(j) (\bar a \mathbf R_i + \hat \mathbf P_i)O_{\gamma}(\hat \mathbf J_j)$$
		$=$	$$\sum_{i<j,\alpha\beta=1,2,3,\gamma=1,...} R_i^{\beta} \frac{\part... ...\partial \hat U_i^{\alpha}} \epsilon_{\alpha\beta} O_{\gamma}(\hat \mathbf J_j)$$
			$$+\sum_{i<j,\gamma}\nabla_{\hat \mathbf U_j} B_{\gamma}(j) \hat \mathbf P_i O_{\gamma}(\hat \mathbf J_j)$$

In the last line of this equation we have made use of the invariance of the total crystal field energy under rotations and therefore substituted $\bar a$ with the strain $\bar \epsilon$. Abbreviating notation we arrive at the final result for the crystal field phonon interaction:

	$$\hat H_{\rm cfph}$$	$=$	$$-\sum_{j,\alpha=1,..,6,\gamma=1,...} G_{\rm cfph}^{\alpha\gamma}(j) \epsilon_{\alpha} O_{\gamma}(\hat \mathbf J_j)$$
			$$-\sum_{i<j,\alpha=1,2,3,\gamma=1,...} \Gamma^{\alpha\gamma}(ij) \hat P_i^{\alpha}a_0^{-1} O_{\gamma}(\hat \mathbf J_j)$$	(99)

Note, that the first part of this equation denotes the coupling of the strain to the crystal field and is commonly known as magnetoelastic interaction [33]. The definition of the static magnetoelastic constants writing explicitly the Voigt notation of the first index $\alpha=(\alpha\beta)$ is

  $$G_{\rm cfph}^{(\alpha\beta)\gamma}(j) = -\frac{1}{2}\sum_{i}(... ... + R_i^{\alpha} \frac{\partial B_{\gamma}(j)}{\partial \hat U_i^{\beta}} )$$ (100)

The dynamic magnetoelastic constants (the crystal field phonon coupling constants) are

  $$\Gamma^{\alpha\gamma}(ij)= -a_0\frac{\partial B_{\gamma}(j)}{\partial \hat U_i^{\alpha}}$$ (101)

Note, that equation (108) makes use of the displacement derivatives of the crystal field parameters and not the strain derivatives found in literature [34,35,36]. In section 12.9 a procedure to calculate these derivatives is described.

Summarizing, and remembering the dimensionless phonon displacement operators $\hat \mathbf u = \hat \mathbf P / a_0$, we can write the total Hamiltonian as

	$$\hat H_{\rm tot}$$	$=$	$$\sum_{i,\gamma} B_{\gamma}(i,0,\dots,0) O_{\gamma}(\hat \mathbf J_i) - \sum_i g_{J_i} \mu_B \hat \mathbf J_i \mathbf H+$$	(102)
		$+$	$$\sum_{i} \hat H_{\rm E}(i) + \frac{1}{2}\sum_{\alpha\gamma=1-6} c^{\alpha\gamma} \epsilon_{\alpha}\epsilon_{\gamma} -$$
		$-$	$$\frac{1}{2}\sum_{i\ne j,\alpha\beta} K_{\alpha\beta}(ij) \hat u_{... ...\gamma}\Gamma^{\alpha\gamma}(ij) \hat u^{\alpha}_i O_{\gamma}(\hat \mathbf J_j)$$
		$-$	$$\sum_{i,\alpha=1-6,\gamma=1,2,3} G_{mix}^{\alpha\gamma}(i) \epsilon_{\alpha} \hat u_{i}^{\gamma} -$$
		$-$	$$\sum_{i,\alpha=1-6,\gamma=1,...} G_{\rm cfph}^{\alpha\gamma}(i) \epsilon_{\alpha} O_{\gamma}(\hat \mathbf J_i)$$

The first line in (110) contains the single ion Hamiltonian (crystal field, Zeeman)

and the second line the phonon Einstein oscillator and elastic energy (also a ”single ion” terms), the third line the interaction terms (phonon, crystal field phonon), the forth line the mixing term and the last line the magnetoelastic term. In a selfconsistent solution it is possible to determine (i) $\epsilon$, (ii) $\langle \hat O_l^m(i) \rangle$ and (iii) $\langle \hat \mathbf u_i \rangle$. This will produce multipolar phase diagrams including the magnetostrictive properties without the need of a detailed investigation of the symmetry adapted Hamiltonian.

Setting zero the derivative of the expectation value of the Hamiltonian (105) with respect to $\bar a$ (i.e. minimizing the energy with respect to strain and rotation) yields the following relations

	$$0$$	$=$	$$\sum_{ij} \frac{c_L(ij)-c_T(ij)}{\vert\mathbf R_{ij}\vert^2} (\mathbf R_{ij}^T\bar \epsilon \mathbf R_{ij})R_{ij}^{\alpha}R_{ij}^{\beta}$$	(103)
			$$+ \frac{c_T(ij)}{2} (R_{ij}^{\alpha} (\bar \epsilon \mathbf R_{ij})^{\beta}+R_{ij}^{\beta} (\bar \epsilon R_{ij})^{\alpha})$$
			$$+2\sum_{ij} \frac{c_L(ij)-c_T(ij)}{\vert\mathbf R_{ij}\vert^2} R_{ij}^{\alpha}R_{ij}^{\beta} \mathbf R_{ij}^T \langle \mathbf P_{i} \rangle$$
			$$+2 c_T(ij) \mathbf R_{ij}^{\beta} \langle \hat P_{i}^{\alpha} \ra... ...ial \hat U_j^{\alpha}} R_j^{\beta} \langle O_{\gamma}(\hat \mathbf J_i) \rangle$$

The index $i$ needs only to go over the atoms in the unit cell, because the crystal structure is periodic. There are 9 equations for $\alpha,\beta=1,2,3$ for nine components $a_{\alpha\beta}$. Thus the coefficients of the strain in equation (111) can be evaluated numerically. For each mean field iteration the strain $\epsilon_{\alpha\beta}$ components can be calculated from equation (111) and inserted into (105) until selfconsistency is achieved.

Considering only the strain $\bar \epsilon$ we can make use of the elastic constants in forming the derivative of the Hamiltonian, we get 6 equations for $\alpha=1,...,6$, from which the 6 strain components can be determined.

	$$\sum_{\beta=1,..6} c^{\alpha\beta} \epsilon_{\beta}$$	$=$	$$\sum_{i,\delta=1,2,3} G_{mix}^{\alpha\delta}(i)\langle \hat u_{i}^{\delta} \rangle$$	(104)
		$+$	$$\sum_{i,\gamma=1,...} G_{\rm cfph}^{\alpha\gamma}(i) \langle O_{\gamma}(\hat \mathbf J_i) \rangle$$