Dynamical Susceptibility and Excitations - Formalism for mcdisp

This section is describes the formalism used in the calculation of the magnetic excitations. Because the procedure is not standard, we list the most important formulas.

We assume a quantum mechanical system that can be described by the Hamiltonian

  $$\hat \mathcal H=\sum_{n=1}^N \hat \mathcal H(n) -\frac{1}{2} \su... ..._{n'} - \mathbf R_n) \hat \mathcal I_{\alpha}^n \hat \mathcal I_{\beta}^{n'}.$$ (237)

The first term $\hat \mathcal H(n)$ denotes the Hamiltonian of a subsystem $n$ (e.g. an ion, or cluster of ions). The second term describes a bilinear interaction between different subsystems through the operators $\hat \mathcal I_{\alpha}^n$, with $\alpha = 1,2,...,m$. The operators $\hat \mathcal H(n)$ and $\hat \mathcal I_{\alpha}^n$ act in the subspace $n$ of the Hilbert space, i.e. $[\hat \mathcal I_{\alpha}^n,\hat \mathcal I_{\alpha}^{n'}]=0$, $[\hat \mathcal H(n),\hat \mathcal I_{\alpha}^{n'}]=0$ and $[\hat \mathcal H(n),\hat \mathcal H(n')]=0$ for $n \neq n'$³⁶. For example, in the case of a Heisenberg exchange between magnetic ions we would identify the set of operators with $\alpha=1,2,3$ with the three components of the spin: $\hat \mathcal I_1 \leftrightarrow \hat S_x, \hat \mathcal I_2 %%@ \leftrightarrow \hat S_y, \hat \mathcal I_3 \leftrightarrow \hat S_z$. The beauty of the analysis which follows is that it can be applied to almost any Hamiltonian of the form (253). The analysis of complex magnetic systems can thus be attempted by starting from a simple form such as the Heisenberg model and by introducing, step-by-step, more complexity into the model. For example, anisotropy and interactions with extended range can be introduced by modifying ${\mathcal J}_{\alpha\beta}(\mathbf R_n - \mathbf R_{n'})$, higher order operators can be introduced by extending the index range for $\alpha$, and a complex single-ion term $\hat \mathcal H(n)$ may be added. Another example for a Hamiltonian (253) is the problem of lattice dynamics, which can be treated in the framework of this formalism by identifying the operators $\hat \mathcal I^n_{\alpha}$ with the atomic displacements $u^{n}$. Here the index $\alpha$ is not necessary and $n$ refers to both, the atomic position index and the spatial coordinate of the displacement, $n=(1,x),(1,y),(1,z),(2,x),(2,y),(2,z), ...$. Note that this can be done, because the three spatial components of the displacement operators commute with each other (in contrast to the components of the spin) and each displacement component acts in its own subspace of the Hilbert space. The kinetic energy will be part of the single ion term $\hat \mathcal H(n)$. Allowing more complexity to the system, both the spin and lattice degrees of freedom can be introduced and spin-phonon interactions can be handled by the theory.

The main limitation of the approach is that it neglects fluctuations associated with phase transitions and quantum disorder. We are primarily concerned, therefore, with excitations associated with a well-ordered ground state.

The translational symmetry of the system is represented by a Bravais lattice (which, in general, will be a superlattice of a crystal lattice). The position of subsystem $n$ can be specified by a lattice vector $\ensuremath{\boldsymbol\ell}$ and a basis vector $\mathbf b_s$. The latter is the position of $n$ relative to $\ensuremath{\boldsymbol\ell}$. The index $s$ ($s=1,2,...,N_b$) labels the subsystems within the unit cell.

The calculation of the excited states of the system starts from a mean-field model for the ground-state order. We define a mean field acting on each subsystem by

  $$H_{\alpha}^s=\sum_{\ensuremath{\boldsymbol\ell}'s'\beta} {\mathc... ...ath{\boldsymbol\ell}-\mathbf b_s) \langle \hat \mathcal I^{s'}_{\beta}\rangle,$$ (238)

where $\langle \hat \mathcal I^{s'}_{\beta}\rangle$ represents the thermal expectation value at a temperature $T$ in the mean field acting on subsystem $s'$. Note that the mean field is periodic in the lattice, so does not depend on $\ensuremath{\boldsymbol\ell}$. The mean-field Hamiltonian for subsystem $s$ is then given by

  $$\hat \mathcal H^{\rm MF}(s) =\hat \mathcal H(s) - \sum_{\alpha=1}^m H_{\alpha}^s \hat \mathcal I^s_{\alpha}$$ (239)

The mean-field ground state is obtained from the self-consistent solution of (254) and (255). This iterative procedure is illustrated in fig. 32. The mean field Hamiltonian (255) for the subsystem $s$ is used to calculate the thermal expectation values $\langle \hat \mathcal I^{s'}_{\beta}\rangle$ for the initial mean field acting on all subsystems $s'=1,2, ...,m$. Equation (254) is then used to calculate a new set of mean fields. These are again used in (255), and the procedure repeated until convergence is reached to within some specified precision. The free energy of the mean field ground state is evaluated and compared to that of other solutions obtained at the same temperature (computed from other initial states and superlattices). The solution with the lowest free energy corresponds to the stable ground state.

We now turn to the excited states. From linear response theory it can be shown [1, page 143] that the excited states are poles of the dynamical susceptibility, which is defined by

  $$\chi_{BA}^{}(\omega)=\lim_{\varepsilon\to0^+}\Bigg[\sum_{aa'}^{E_... ..._a^{}-n_{a'}^{}) +\frac{i\varepsilon}{\omega+i\varepsilon}\chi_{BA}'(el)\Bigg]$$ (240)

where

  $$\chi_{BA}'(el)=\frac{1}{kT} \sum_{aa'}^{E_a= E_{a'}}\langle a\ver... ... a'\rangle \langle a'\vert\hat{A}-\langle \hat{A}\rangle\vert a\rangle\,n_a^{}$$ (241)

and

  $$n_a^{}=\frac{\exp(-E_a^{}/kT)}{\sum_{a'}\exp(-E_{a'}/kT)};\qquad \langle \hat{A}\rangle=\sum_a\langle a\vert\hat{A}\vert a\rangle\,n_a^{}$$ (242)

Here the energy levels and eigenstates of the Hamiltonian (253) are denoted by $E_{{a}}$ and $\vert{a} \rangle$, respectively. $n_{{a}}$ is the corresponding Boltzmann occupation probability. $\hat A$ and $\hat B$ are quantum mechanical operators describing the perturbation to the Hamiltonian and the response of the system according to the general concept of linear response theory [1]. The expression (256) is based on a system with well defined energy levels implying that the poles of $\chi_{BA}^{}(z=\omega+i\varepsilon)$ are all lying on the real axis, or that the absorptive part of the response function

  $$\chi_{BA}''(\omega)\equiv\lim_{\varepsilon\to0^+}\frac{1}{2i} \left[\chi_{BA}^{}(z)-\chi_{AB}^{}(-z^\ast)\right]$$ (243)

becomes a sum of $\delta$-functions, which are only non-zero when $\hbar\omega$ is equal to the excitation energies $E_{a'}-E_a$. If the susceptibility contains an elastic contribution, then the function $\chi_{BA}''(\omega)/\omega=\pi\chi_{BA}'(el)\delta(\omega)$ in the zero frequency limit. In any realistic, interacting system the energy levels are no longer discrete states and fluctuations will cause spontaneous transitions between the different levels. Nonzero probabilities for such transitions may be accounted for in a phenomenological way by replacing $E_{a'}-E_a$ in equation (256) by $E_{a'}-E_a-i\Upsilon_{a'a}$, where $\Upsilon_{a'a}\ge0$. In this approximation the response becomes a sum of Lorentzians

  $$\chi_{BA}''(\omega)\simeq \sum_{aa'} \frac{\langle a\vert\hat{B}\... ...+\frac{\hbar\omega\,\Upsilon_0^{}}{(\hbar\omega)^2+\Upsilon_0^2}\chi_{BA}'(el)$$ (244)

The same result is obtained if $\varepsilon$ is kept as a nonzero positive quantity in (256) instead of taking the limit $\varepsilon\to0^+$, i.e. if assuming $\varepsilon=\Upsilon_{a'a}/\hbar$ in the different terms in the sum and $\varepsilon=\Upsilon_0/\hbar$ in the elastic term.

Because of the periodicity of our system we define generalized susceptibilities $\chi_{\alpha\beta}^{ss'}({\mathbf Q},\omega)\equiv\chi_{BA}$ by choosing the Fourier transform operators

	$$\hat A$$	$=$	$$\frac{1}{\sqrt{N}}\exp({\rm i}{\mathbf Q} \cdot \mathbf b_{s'})\s... ...{\boldsymbol\ell}')\hat \mathcal I^{(\ensuremath{\boldsymbol\ell}' s')}_{\beta}$$	(245)
	$$\hat B$$	$=$	$$\frac{1}{\sqrt{N}}\exp({\rm i}{\mathbf Q} \cdot \mathbf b_{s})\su... ...th{\boldsymbol\ell})\hat \mathcal I^{(\ensuremath{\boldsymbol\ell}s)}_{\alpha},$$	(246)

where $N$ is the number of unit cells. It will also be convenient to introduce the Fourier transform of the two-body interaction

  $${\mathcal J}_{\alpha\beta}^{ss'}({\mathbf Q})=\sum_{\ensuremath{\... ...i}\mathbf Q \cdot (\ensuremath{\boldsymbol\ell}'+\mathbf b_{s'}-\mathbf b_s)\}$$ (247)

We have arbitrarily chosen ${\ensuremath{\boldsymbol\ell}} = 0$ since ${\mathcal J}({\mathbf Q})$ is the same for all $\ensuremath{\boldsymbol\ell}$ due to the translational symmetry.

Figure 32: Illustration of the iterative flow used for solving the mean-field Hamiltonian eq. (255)

The calculation of the dynamical susceptibility³⁷ $\overline{\chi}({\mathbf Q},\omega)$ from the Hamiltonian (5) is carried out within the mean field – random phase approximation (MF–RPA) [1,74]. This approximation neglects correlations in the differences $\hat \mathcal I^n(t) - \langle \hat \mathcal I^n (t) \rangle$ of different subsystems $n$. In this approach the dynamical susceptibility $\overline{\chi}({\mathbf Q},\omega)$ for a primitive lattice ($s=s'=0$) can be calculated from the solution to

  $$1=\left [ \overline{\chi}^0(\omega)^{-1}-\overline{\mathcal J}({\mathbf Q}) \right]\overline{\chi}({\mathbf Q},\omega),$$ (248)

where $\overline{\chi}^0(\omega)$ is the usual single ion magnetic susceptibility tensor. This equation can be written for the more general case of several subsystems ( $s=1,2,3,4,...,N_B$) as

  $$\delta_{ss'}= \sum_{s''=1}^{N_B}\left[ \delta_{ss''}[\overline{\c... ...l J}}^{ss''}({\mathbf Q}) \right] \overline{\chi}^{s''s'}({\mathbf Q},\omega),$$ (249)

or, in index notation, to

  $$\delta_{ss'}\delta_{\alpha\beta}= \sum_{s''=1}^{N_B}\sum_{\delta=... ...a}^{ss''}({\mathbf Q}) \right] \chi_{\delta\beta}^{s''s'}({\mathbf Q},\omega),$$ (250)

where $\chi^{s}_{\alpha\beta}(\omega)$ is the subsystem susceptibility (the same for all $\ensuremath{\boldsymbol\ell}$), given by

  $$\chi^{s}_{\alpha\beta}(\omega)= \sum_{jj'} \frac{\langle j\vert\h... ...ta}\rangle\vert j\rangle}{\epsilon_{j'}-\epsilon_j-\hbar \omega} (p_j-p_{j'}).$$ (251)

where for the sake of simplicity we omit the $s$ index on all quantities on the right-hand side. Here $\epsilon_j$ and $\epsilon_{j'}$ are energy levels of the subsystem $s$ as calculated self-consistently within the mean-field theory using the Hamiltonian (255), $\vert j\rangle$ and $\vert j'\rangle$ denote the corresponding eigenstates and $p_j$ the corresponding population numbers:

  $$p_{j} = \frac{\exp(-\epsilon_{j}/{k}T)}{\sum_{j'}\exp(-\epsilon_{j'}/{k}T)}.$$ (252)

The writing of (267) has been simplified in two ways. The obvious one is that $\omega$ should read $\omega+i \varepsilon$ where $\varepsilon\to0^+$. Secondly, the elastic contribution is included in (267) by assuming the use of the following convention: $\epsilon_{j'}$ is being replaced by $\epsilon_{j'}+d$ in all terms where $\epsilon_{j'}=\epsilon_j$. The shift in energy introduced is $\vert d\vert\ll k^{}T$ and hence $p_j-p_{j'}=p_j\,d/k^{}T$ to leading order. Notice that the matrix elements of the thermal expectation values in (267) are only nonzero in the special cases of $j=j'$. Using the two conventions equation (267) becomes equivalent to (256) in the limit of $d\to0$ (after taking the limit $\varepsilon\to0^+$). Since the expectation values are only needed in (267) when considering the elastic contribution, we may use this fact to signal that the second convention has to be applied whenever the expectation values are subtracted from the operators.

In order to evaluate equations (264)-(267) without producing a numerical divergence it is necessary to add to $\omega$ a small imaginary constant $\omega \rightarrow \omega+i\epsilon$ and insert this into equation (267). Technically this is done in the jjjpar::chi0() function in source file bin/src/jjjpar_basmodfunc.cpp. If the option -r $\epsilon$ is used, the program McDisp calculates the above expression for every energy and stores the result in ./results/mcdisp.dsigma. Also program singleion can be used to evaluate the magnetic subsystem susceptibility according to (267), see section 19.

If the option -r is not used, the program mcdisp uses only the extremely fast DMD (Dynamical Matrix Diagonalisation) algorithm[39,75] to calculate excitation energies and intensities and store the result in mcdisp.qom, mcdisp.qei, etc.. The flowing chart of such a calculation is shown in fig. 33 and the formalism is outlined hereafter: