The total Crystal field Phonon Interaction Hamiltonian

Also the crystal field may be written and expanded in terms of the strain $\epsilon$ and the $\vec 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 (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 in a crystal $j=N_{\rm nuclei}+1,... N_{\rm nuclei}+N_{\rm magnetic ions}$)

$$H = H_{\rm ph} +\sum_{j,\gamma} B_{\gamma}(j,\vec U_1,...,\vec U_N) O_{\gamma}(\vec J_i)$$ (97)

$$= H_{\rm ph} +\sum_{j,\gamma} B_{\gamma}(j,0,\dots,0) O_{\gamma}(\vec J_j) + H_{\rm cfph}$$

$$H_{\rm cfph} = \sum_{i<j,\gamma}\nabla_{\vec U_i} B_{\gamma}(j) \vec U_i O_{\gamma}(\vec J_j)$$

$$= \sum_{i<j,\gamma}\nabla_{\vec U_i} B_{\gamma}(j) (\bar a \vec R_i + \vec P_i)O_{\gamma}(\vec J_j)$$

$$= \sum_{i<j,\alpha\beta=1,2,3,\gamma=1,...} R_i^{\beta} \frac{\part... ...{\gamma}(j)}{\partial U_i^{\alpha}} \epsilon_{\alpha\beta} O_{\gamma}(\vec J_j)$$

$$+\sum_{i<j,\gamma}\nabla_{\vec U_j} B_{\gamma}(j) \vec P_i O_{\gamma}(\vec J_j)$$

$$= -\sum_{j,\alpha=1,..,6,\gamma=1,...} G_{\rm cfph}^{\alpha\gamma}(j) \epsilon_{\alpha} O_{\gamma}(\vec J_j)$$

$$-\sum_{i<j,\alpha=1,2,3,\gamma=1,...} \Gamma^{\alpha\gamma}(ij) P_i^{\alpha}a_0^{-1} O_{\gamma}(\vec J_j)$$ (98)

The definition of the static magnetoelastic constants writing explicitely 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^{\beta} \frac{\partial B_{\gamma}(j)}{\... ...^{\alpha}} + R_i^{\alpha} \frac{\partial B_{\gamma}(j)}{\partial U_i^{\beta}} )$$ (99)

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

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

Note, that equation (103) makes use of the displacement derivatives of the crystal field parameters and not the strain derivatives found in literature [32,33,34]. 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$.

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

$$H = \sum_{i,\gamma} B_{\gamma}(i,0,\dots,0) O_{\gamma}(\vec J_i) + \sum_{i} H_{\rm E}(i) +$$ (101)

$$+ \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) u_{i}^{\... ...j,\alpha=1-6,\gamma}\Gamma^{\alpha\gamma}(ij) u^{\alpha}_i O_{\gamma}(\vec J_j)$$

$$- \sum_{i,\alpha=1-6,\gamma=1,2,3} G_{mix}^{\alpha\gamma}(i) \epsilon_{\alpha} u_{i}^{\gamma} -$$

$$- \sum_{i,\alpha=1-6,\gamma=1,...} G_{\rm cfph}^{\alpha\gamma}(i) \epsilon_{\alpha} O_{\gamma}(\vec J_i)$$

The first line in (106) contains the single ion Hamiltonian (crystal field, phonon), and the second line the elastic energy (also a ”single ion” term), 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 should be possible to determine (i) $\epsilon$, (ii) $\langle O_l^m(i) \rangle$ and (iii) $\langle \vec u_i \rangle$. This should be able to 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 (102) 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\vec R_{ij}\vert^2} (\vec R_{ij}^T\bar \epsilon \vec R_{ij})R_{ij}^{\alpha}R_{ij}^{\beta}$$ (102)

$$+ \frac{c_T(ij)}{2} (R_{ij}^{\alpha} (\bar \epsilon \vec R_{ij})^{\beta}+R_{ij}^{\beta} (\bar \epsilon R_{ij})^{\alpha})$$

$$+2\sum_{ij} \frac{c_L(ij)-c_T(ij)}{\vert\vec R_{ij}\vert^2} R_{ij}^{\alpha}R_{ij}^{\beta} \vec R_{ij}^T \langle \vec P_{i} \rangle$$

$$+2 c_T(ij) \vec R_{ij}^{\beta} \langle P_{i}^{\alpha} \rangle + \... ...ma}(i)}{\partial U_j^{\alpha}} R_j^{\beta} \langle O_{\gamma}(\vec 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 (107) can be evaluated numerically. For each mean field iteration the strain $\epsilon_{\alpha\beta}$ components can be calculated from equation (107) and inserted into (102) 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_{\gamma=1,..6} c^{\alpha\gamma} \epsilon_{\gamma} = \sum_{i,\gamma=1,2,3} G_{mix}^{\alpha\gamma}(i)\langle u_{i}^{\gamma} \rangle$$ (103)

$$+ \sum_{i,\gamma=1,...} G_{\rm cfph}^{\alpha\gamma}(i) \langle O_{\gamma}(\vec J_i) \rangle$$