The effects of thermal and correlated noise on magnons in a quantum ferromagnet

The investigation of spin waves and magnons has lead to the emerging field of magnonics [1–3] which aims to enable magnons as information carriers for both classical [4, 5] and quantum [6] information technology. Since magnons do not carry charge their interaction and dissipation is minimal compared to electronic circuits and thus magnonics could enable information processing with hugely reduced power consumption [3, 4]. Furthermore magnonics could enable transport [7–10], processing [11, 12] and storage of quantum bits in the same platform of spin systems, since single-spin systems, such as the nitrogen- or silicon-vacancy centre in diamond, are strong candidates for future quantum bits. Magnons also enable new approaches to entanglement and quantum processing since (unlike photons) they can be created with a simple spin-gate operation on a spin. A rigorous investigation of their response to quantum noise including spatially correlated noise is therefore an important step towards further progress and miniaturisation of magnonic systems and devices.

Using quantum master equations to calculate macroscopic properties of solid state systems can yield details about the material properties and dynamics. This is a distinctly different application from the calculation of expectation values in small quantum systems of only a few states. For certain systems such as the quantum ferromagnet, analytical solutions can be obtained from master equation approaches. Numerical methods for large systems are possible via mapping of master equations to a quantum jump formalism [13]. The effects of uncorrelated noise in a quantum antiferromagnet have been investigated analytically [14] and shown the applicability and usefulness of this technique. In addition, master equations techniques have been applied in other condensed matter scenarios such as topologically ordered systems [15–17].

The effects of spatially correlated noise are increasingly relevant for quantum technology devices, which aim for an increasing density of controlled quantum systems and thus the simplified model of uncorrelated noise are more likely to break down. Spatially correlated noise has shown to produce interesting and fundamentally different dynamics to uncorrelated noise, even in the Markovian regime. In ion traps for example the occurrence of decoherene-free subspaces has inspired quantum computation solutions [18, 19]. In the field of quantum metrology a re-instatement of the superior Heisenberg precision scaling has been proven possible in the presence of spatial noise correlations [20]. In spin chains and light-harvesting complexes spatial noise correlations have shown to enable robust transport through protected states [21–23]. These results pose the question how correlated noise affects the dynamics of magnons and macroscopic quantities in a quantum ferromagnet

This paper is structured as follows: in section 2 we introduce the spin-wave Hamiltonian, in section 3 we discuss its interaction with a thermal environment, in section 4 we derive and solve a master equation for magnons, in section 5 we discuss the macroscopic magnetisation, in section 6 we point out relevant differences between ferro- and antiferromagnetic behaviour and reach our conclusions in section 7.

We begin by introducing the key concepts and the Heisenberg model Hamiltonian for spins in real space with nearest-neighbour interaction. We then show how this Hamiltonian maps via the Holstein–Primakoff transform with a spin-wave approximation to a bosonic system. Subsequent transformation into k-space via Fourier lattice transform of the operators diagonalises the Hamiltonian and yields the dispersion relation of the system. This is the Hamiltonian describing 'magnons', the elementary collective magnetic excitations.

We start with a Heisenberg model for the spins in the quantum ferromagnet with only nearest-neighbour interaction and a uniform magnetic field B in the negative z-direction:

$$\begin{eqnarray}&&{H}_{S}=J\displaystyle \sum _{\langle {\bf{r}},{\bf{r}}^{\prime} \rangle }{{\bf{S}}}_{{\bf{r}}}\cdot {{\bf{S}}}_{{\bf{r}}^{\prime} }-\displaystyle \sum _{{\bf{r}}}\gamma {{BS}}_{{\bf{r}}}^{z},\end{eqnarray} \tag{ 1 }$$

where J < 0 for a ferromagnet. Antiferromagnetic behaviour (J > 0) has been investigated similarly [14]. Here $\gamma ={g}_{e}{\mu }_{B}/{\hslash }$ is the gyromagnetic ratio and ${{\bf{S}}}_{{\bf{r}}}=({S}_{x},{S}_{y},{S}_{z})$ is the vector spin-operator of at position r = (x, y, z).

This Hamiltonian can be mapped to a system of interacting bosons via the Holstein–Primakoff transformation [24, 25] for spins S ≥ 1/2:

$$\begin{eqnarray}\begin{array}{rcl}{S}_{{\bf{r}}}^{z} & = & S-{n}_{{\bf{r}}}\\ {S}_{{\bf{r}}}^{+} & = & \sqrt{2S}\phi ({n}_{{\bf{r}}}){a}_{{\bf{r}}}\\ {S}_{{\bf{r}}}^{-} & = & \sqrt{2S}{a}_{{\bf{r}}}^{\dagger }\phi ({n}_{{\bf{r}}}),\end{array}\end{eqnarray} \tag{ 2 }$$

where ${n}_{{\bf{r}}}={a}_{{\bf{r}}}^{\dagger }{a}_{{\bf{r}}}$ and:

$$\begin{eqnarray}&&\phi =\sqrt{1-\displaystyle \frac{{n}_{{\bf{r}}}}{2S}}.\end{eqnarray} \tag{ 3 }$$

This function is then approximated by a series expansion to first order in the normal-ordered number operator [25]:

$$\begin{eqnarray}&&\phi \approx 1-(1-\sqrt{1-1/2S}){n}_{{\bf{r}}}\qquad \mathrm{for}\ \,\langle {a}_{{\bf{r}}}^{\dagger }{a}_{{\bf{r}}}\rangle \ll 2S.\end{eqnarray} \tag{ 4 }$$

We note that this linear spin-wave approximation limits the regime to low excitation numbers of magnons relative to their maximum set by the spin $2S$. This is fulfilled both in the ordered ferromagnetic equilibrium state as well as in the equilibration processes considered here, as long as the bath temperature is 'cold' relative to the Debye temperature of the material, i.e. such that the bath fluctuations cannot excite large numbers of magnons into the system. Furthermore, we will be working in a 3D system, leading to spontaneous symmetry breaking at finite temperature further justifying the spin-wave treatment to be qualitatively correct.

With the linear spin-wave approximation the Hamiltonian (1) becomes (for a D-dimensional system):

$$\begin{eqnarray}&&{H}_{{\rm{LSW}}}={E}_{0}+(\gamma B-2{JDS})\displaystyle \sum _{{\bf{r}}}{a}_{{\bf{r}}}^{\dagger }{a}_{{\bf{r}}}+{JS}\displaystyle \sum _{\langle {\bf{r}},{\bf{r}}^{\prime} \rangle }({a}_{{\bf{r}}}^{\dagger }{a}_{{\bf{r}}^{\prime} }+{a}_{{\bf{r}}^{\prime} }^{\dagger }{a}_{{\bf{r}}}),\end{eqnarray} \tag{ 5 }$$

where the constant energy shift ${E}_{0}={{JNDS}}^{2}+\gamma {BSN}$ with the total number of spins N originates from a summation over a constant.

The Hamiltonian is then diagonalized for D = 3 by replacing ${a}_{{\bf{r}}}$ and ${a}_{{\bf{r}}}^{\dagger }$ with their three-dimensional Fourier lattice transform [26–28]:

$$\begin{eqnarray}\begin{array}{rcl}{a}_{{\bf{r}}} & = & {a}_{x,y,z}=\displaystyle \frac{1}{\sqrt{N}}\displaystyle \sum _{{\bf{k}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{a}_{{\bf{k}}}\\ {a}_{{\bf{r}}}^{\dagger } & = & \displaystyle \frac{1}{\sqrt{N}}\displaystyle \sum _{{\bf{k}}}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{a}_{{\bf{k}}}^{\dagger }.\end{array}\end{eqnarray} \tag{ 6 }$$

Taking into consideration a simple cubic lattice as the discrete spin lattice, the detailed expression for the wave vector elements is ${k}_{x}=\tfrac{2\pi }{{N}_{x}{d}_{x}}{n}_{x}$ and ${r}_{x}={d}_{x}{\tilde{n}}_{x}$, where N_x is the number of sites in x-direction, d_x is the lattice constant in x-direction and ${n}_{x},{\tilde{n}}_{x}\in {\mathbb{N}}$ are the integer summation indices. Analogous expressions apply to y- and z-directions. For the Hamiltonian, equation (5), we choose the three-dimensional nearest-neighbour parameterisation ${\sum }_{\langle {\bf{r}},{\bf{r}}^{\prime} \rangle }\to {\sum }_{{\bf{r}}}{\sum }_{\mu =1}^{3}$ and ${\bf{r}}^{\prime} \to {\bf{r}}+{d}_{\mu }{\hat{{\bf{r}}}}_{\mu }$, where ${\hat{{\bf{r}}}}_{\mu }$ is the unit vector in x-, y- or z-direction and d_μ the respective lattice constant. One then diagonalises the Hamiltonian in k-space using ${\sum }_{{\bf{r}}}\exp [{\rm{i}}({\bf{k}}-{\bf{k}}^{\prime} )\cdot {\bf{r}}]={\delta }_{{k}_{x},{k}_{x}^{{\prime} }}\,{\delta }_{{k}_{y},{k}_{y}^{{\prime} }}\,{\delta }_{{k}_{z},{k}_{z}^{{\prime} }}$ and finds:

$$\begin{eqnarray}&&{H}_{{\rm{LSW}}}={E}_{0}+\displaystyle \sum _{{\bf{k}}}\omega ({\bf{k}}){a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mathrm{with}\ \ \quad \omega ({\bf{k}})=\left[\gamma B+2{JS}\displaystyle \sum _{\mu =1}^{3}(\cos \,{k}_{\mu }{d}_{\mu }-1)\right].\end{eqnarray} \tag{ 8 }$$

This is the diagonal linear spin-wave Hamiltonian. The excitations of these uncoupled harmonic oscillators are called magnons as they represent the elementary magnetic excitation. The ${\bf{k}}$-dependent energy of the magnons ω(k) defines their dispersion relation.

In this section we extend our model of the system of a ferromagnet to be an open quantum system, which allows us to model the interaction between the spins and the lattice phonons, which we will assume to be in a thermal equilibrium. Instead of assuming an ad hoc master equation we will show a more rigorous treatment and start by assuming an interaction Hamiltonian between the ferromagnetic system and an environment of bosonic phonon-modes, figure 1. This is the basis for the derivation of a master equation. Any general interaction Hamiltonian [29] can always be written in the mathematical form ${H}_{{\rm{int}}}={\sum }_{j}{s}_{j}{B}_{j}$, where multiple system operators s_j couple to environmental bath operators B_j. Defining these operators s_j and B_j in this form of the coupling is the starting point for a treatment with the Bloch–Redfield formalism, which allows us to write down the master equation directly from this form [30–33].

Zoom In Zoom Out Reset image size

We will start with an interaction Hamiltonian with spin-operators in real space and include an arbitrary spatial correlation function $f(| {\bf{r}}-{\bf{r}}^{\prime} | )$ of the system-environment couplings, which defines how the coupling between a spin and an environmental mode behaves as a function of distance $| {\bf{r}}-{\bf{r}}^{\prime} |$ between them. We then show how this form simplifies by transformation into k-space such that magnons only couple to environmental modes with the same wave vector k. The spatial correlation function is then contained in the k-dependent coupling constant. This simplification is essential in order to solve the resulting Bloch–Redfield equations analytically later.

We note that introducing this spatial decay function $f(| {\bf{r}}-{\bf{r}}^{\prime} | )$ in the system-environment couplings is subtly different to the definition of a spatial decay function in the environmental noise correlations as done elsewhere [30]. In the latter case the decay function is bound by multipartite correlation rules while here the decay function of couplings is unrestricted in its functional form and will always result in a Lindblad form and a physical time evolution.

We assume that our system of spins couples to a large number of bosonic environmental modes, such as phonons. This coupling of the spins to the lattice phonons is the common spin-boson model [34–38] where the spin operator ${\bf{S}}=({S}^{x},{S}^{y},{S}^{z})$ couples linearly to the lattice displacement operator ${\bf{R}}=(X,Y,Z)$. Details of an underlying microscopic model have been discussed [39] as a local exchange interaction between the electron spin and magnetic moment and the local couplings between the electronic charge and lattice displacements. The interaction of the form ${S}^{x}X={S}^{x}({A}^{\dagger }+A)$, with the local position operator X of the environmental mode, is one of the common spin-boson models [34, 35]. We chose this as the relevant example since an interaction ${S}^{y}Y$ yields completely analogous dynamics and an interaction ${S}^{z}Z$ would lead to terms which are negligible in the spin-wave regime in analogy to the antiferromagnetic case [14], see appendix A for details. We adopt capital letters for the environmental creation operators ${A}_{{\bf{r}},j}^{\dagger }$ of mode j at spatial position ${\bf{r}}$ to distinguish them from the system modes ${a}_{{\bf{r}}}^{\dagger }$. At each position there is a collection of phonon modes with different energies ${\omega }_{j}$ and respective coupling strength $g({\omega }_{j})$:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{r}},{\bf{r}}^{\prime} }f(| {\bf{r}}-{\bf{r}}^{\prime} | ){S}_{r}^{x}({A}_{{\bf{r}}^{\prime} ,j}^{\dagger }+{A}_{{\bf{r}}^{\prime} ,j}).\end{eqnarray} \tag{ 9 }$$

After the Holstein–Primakoff transformation with linear spin-wave approximation this interaction reads:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{r}},{\bf{r}}^{\prime} }f(| {\bf{r}}-{\bf{r}}^{\prime} | )({a}_{{\bf{r}}}^{\dagger }+{a}_{{\bf{r}}})({A}_{{\bf{r}}^{\prime} ,j}^{\dagger }+{A}_{{\bf{r}}^{\prime} ,j}).\end{eqnarray} \tag{ 10 }$$

Using the secular approximation this can be simplified further. The secular approximation allows us to neglect coupling or noise terms (i.e. off-diagonal superoperator elements) which are small relative to the difference of two on-diagonal superoperator elements. For a weak system-environment coupling g(ω_j) relative to the magnon energy we can neglect those terms that create both a magnon and an environmental phonon, i.e. terms proportional to ${a}^{\dagger }{A}^{\dagger }$ and ${aA}$, since these terms create small off-diagonal superoperator elements proportional to g(ω_j) corresponding to diagonal elements separated by the magnon energy. We only take into account terms where a magnon is created and an environmental phonon annihilated or vice versa, i.e. ${a}^{\dagger }A$ and ${{aA}}^{\dagger }$ terms. This is more intuitively understood in the interaction picture of the system and environment where the ${a}^{\dagger }{A}^{\dagger }$ and ${aA}$ terms become fast oscillating terms with a frequency corresponding to the sum of magnon and phonon energy. The fast oscillation averages out their effect. The ${a}^{\dagger }A$ and ${{aA}}^{\dagger }$ terms are slow or non-oscillating terms since two counter-rotating factors cancel each other out. In this way the secular approximation can also be regarded as a rotating-wave approximation. Neglecting the ${a}^{\dagger }{A}^{\dagger }$ and ${aA}$ terms is a fairly common procedure in master equation approaches [40, 41], mathematically due to the Riemann–Lebesgue lemma [29]. The interaction Hamiltonian becomes:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{r}},{\bf{r}}^{\prime} }f(| {\bf{r}}-{\bf{r}}^{\prime} | ){a}_{{\bf{r}}}^{\dagger }{A}_{{\bf{r}}^{\prime} ,j}+\ {\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 11 }$$

We note that for this approximation to be generally correct the condition $g({\omega }_{j})\ll \gamma B$ is required since there are otherwise system magnons with an energy that is not large enough to justify the separate scales required for the secular approximation. In this case the exact form of the coupling becomes more relevant and ${a}^{\dagger }{A}^{\dagger }$ terms as well as ${aA}$ terms need to be considered. However, this condition can be fulfilled even for magnetic fields small compared to the spin–spin coupling $J\gg \gamma B\gg g({\omega }_{j})$.

Again we make use of the Fourier lattice transform, equation (6), and rewriting ${{\rm{e}}}^{{\rm{i}}({\bf{k}}^{\prime} \cdot {\bf{r}}^{\prime} -{\bf{k}}\cdot {\bf{r}})}={{\rm{e}}}^{{\rm{i}}{\bf{k}}^{\prime} \cdot ({\bf{r}}^{\prime} -{\bf{r}})}{{\rm{e}}}^{{\rm{i}}({\bf{k}}^{\prime} -{\bf{k}})\cdot {\bf{r}}}$ we find:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{k}},{\bf{k}}^{\prime} }\displaystyle \sum _{{\bf{r}},{\bf{r}}^{\prime} }{{\rm{e}}}^{{\rm{i}}{\bf{k}}^{\prime} \cdot ({\bf{r}}^{\prime} -{\bf{r}})}f(| {\bf{r}}-{\bf{r}}^{\prime} | )\displaystyle \frac{1}{N}{{\rm{e}}}^{{\rm{i}}({\bf{k}}^{\prime} -{\bf{k}})\cdot {\bf{r}}}{a}_{{\bf{k}}}^{\dagger }{A}_{{\bf{k}}^{\prime} ,j}+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 12 }$$

We then substitute ${\bf{u}}={\bf{r}}^{\prime} -{\bf{r}}$ and write:

$$\begin{eqnarray}&&\displaystyle \sum _{{\bf{r}}^{\prime} =1}^{N}{{\rm{e}}}^{{\rm{i}}{\bf{k}}^{\prime} \cdot ({\bf{r}}^{\prime} -{\bf{r}})}f(| {\bf{r}}-{\bf{r}}^{\prime} | )=\displaystyle \sum _{{\bf{u}}=-N/2}^{N/2-1}{{\rm{e}}}^{{\rm{i}}{\bf{k}}^{\prime} \cdot {\bf{u}}}f(| {\bf{u}}| )\end{eqnarray} \tag{ 13 }$$

or, to be more precise, we substitute for each component ${r}_{\mu }^{{\prime} }-{r}_{\mu }={d}_{\mu }{\tilde{n}}_{\mu }^{{\prime} }-{d}_{\mu }{\tilde{n}}_{\mu }={d}_{\mu }{\hat{n}}_{\mu }={u}_{\mu }$, where d_μ is the lattice constant of dimension μ, the variables ${\tilde{n}}_{\mu },{\tilde{n}}_{\mu }^{{\prime} },{\hat{n}}_{\mu }\in {\mathbb{Z}}$ are the integer summation indices and μ takes the values of the three dimensions:

$$\begin{eqnarray}&&\displaystyle \sum _{{\tilde{n}}_{\mu }^{{\prime} }=1}^{{N}_{\mu }}{{\rm{e}}}^{{{\rm{i}}{k}}_{\mu }^{\prime} {d}_{\mu }({n}_{\mu }^{{\prime} }-{n}_{\mu })}f(| {\bf{r}}-{\bf{r}}^{\prime} | )=\displaystyle \sum _{{\hat{n}}_{\mu }=-{N}_{\mu }/2}^{{N}_{\mu }/2-1}{{\rm{e}}}^{{{\rm{i}}{k}}_{\mu }^{{\prime} }{d}_{\mu }{\hat{n}}_{\mu }}f(| {\bf{u}}| ).\end{eqnarray} \tag{ 14 }$$

In doing so we have assumed that the expression is independent of ${\bf{r}}$, and that the summation over ${\bf{r}}^{\prime}$ always runs over all possible values of the difference ${\bf{r}}^{\prime} -{\bf{r}}$. Physically this assumption means that all edge effects are neglected, which occur, when any r_μ is close to 1 or N_μ. We then identify the remaining summation over r as a Kronecker delta and after performing the summation over ${\bf{k}}^{\prime}$, the interaction Hamiltonian then becomes (appendix B):

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{k}}}\mathop{\underbrace{\displaystyle \sum _{{\bf{u}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{u}}}f(| {\bf{u}}| )}}\limits_{F({\bf{k}})}{a}_{{\bf{k}}}^{\dagger }{A}_{{\bf{k}},j}+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 15 }$$

This is the interaction Hamiltonian in k-space in which we find that magnons of the wave vector ${\bf{k}}$ only couple to environmental modes of the same wave vector. The respective coupling strength is given by the three-dimensional Fourier transform of the coupling correlation function $f(| {\bf{r}}-{\bf{r}}^{\prime} | )$. This Fourier transform $F({\bf{k}})$ is explicitly:

$$\begin{eqnarray}&&F({\bf{k}})=\displaystyle \sum _{{\hat{n}}_{x}={N}_{x}/2}^{{N}_{x}/2-1}\displaystyle \sum _{{\hat{n}}_{y}={N}_{y}/2}^{{N}_{y}/2-1}\displaystyle \sum _{{\hat{n}}_{z}={N}_{z}/2}^{{N}_{z}/2-1}f(| u| ){{\rm{e}}}^{{\rm{i}}({k}_{x}{d}_{x}{\hat{n}}_{x}+{k}_{y}{d}_{y}{\hat{n}}_{y}+{k}_{z}{d}_{z}{\hat{n}}_{z})}.\end{eqnarray} \tag{ 16 }$$

For large N and smooth $f(| u| )$ one can take the continuous limit and integrate. Since we neglected edge effects and assumed the difference ${\bf{u}}={\bf{r}}-{\bf{r}}^{\prime}$ to go over all possible values, this is the thermodynamic limit and integration should go from $-\infty$ to $\infty$, where the normalisation of $f(| u| )$ does not need to be adapted if it decays on a small scale relative to the crystal lengths ${d}_{\mu }{N}_{\mu }$.

$$\begin{eqnarray}&&F({\bf{k}})=\iiint {\rm{d}}{\hat{n}}_{x}\,{\rm{d}}{\hat{n}}_{y}\,{\rm{d}}{\hat{n}}_{z}f(| u| ){{\rm{e}}}^{{\rm{i}}({k}_{x}{{\rm{d}}}_{x}{\hat{n}}_{x}+{k}_{y}{{\rm{d}}}_{y}{\hat{n}}_{y}+{k}_{z}{{\rm{d}}}_{z}{\hat{n}}_{z})}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1}{{{d}}_{x}{{d}}_{y}{{d}}_{z}}\iiint {{\rm{d}}{u}}_{x}\,{{\rm{d}}{u}}_{y}\,{{\rm{d}}{u}}_{z}f(| u| ){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{u}}},\end{eqnarray} \tag{ 18 }$$

where ${V}_{d}={d}_{x}{d}_{y}{d}_{z}$ is the volume of the unit cell and ${\bf{u}}$ carries the unit length while ${\hat{n}}_{\mu }$ is dimensionless.

Therefore, the Fourier transform of $f(| u| )$ depends only on the magnitude of the wave vector, $F({\bf{k}})=F(k)$, allowing further simplifications, appendix C. We note that this isotropy follows from the assumed isotropy of the noise correlation function $f(| u| )$ for the whole lattice. It is not connected to any isotropy of the individual site's system-bath coupling operators.

With this the Hamiltonian is:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{k}}}F(k){a}_{{\bf{k}}}^{\dagger }{A}_{{\bf{k}},j}+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 19 }$$

This interaction Hamiltonian in k-space is much simpler as it only couples magnons of wave vector ${\bf{k}}$ to phonons of the same wave vector and thus there is only one summation over ${\bf{k}}$, despite two summations occurring in the spatially correlated real-space form, equation (11). We therefore use this as the starting point for our master equation.

For the standard Bloch–Redfield approach we can now identify the full Hamiltonian $H={H}_{{\rm{sys}}}+{H}_{{\rm{int}}}+{H}_{{\rm{env}}}$, where the system Hamiltonian H_sys is given by the linear spin-wave Hamiltonian from equation (7), the interaction Hamiltonian H_int is given by equation (19) and the environmental Hamiltonian is defined as the bosonic environmental modes:

$$\begin{eqnarray}&&H={H}_{{\rm{sys}}}+{H}_{{\rm{int}}}+{H}_{{\rm{env}}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{H}_{{\rm{sys}}}={H}_{{\rm{LSW}}}={E}_{0}+\displaystyle \sum _{{\bf{k}}}\omega ({\bf{k}}){a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\mathrm{with}\ \ \quad \omega ({\bf{k}})=\left[\gamma B+2{JS}\displaystyle \sum _{\mu =1}^{3}(\cos \,{k}_{\mu }{d}_{\mu }-1)\right]\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{k}}}F(k){a}_{{\bf{k}}}^{\dagger }{A}_{{\bf{k}},j}+{\rm{h}}.{\rm{c}}.,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{H}_{{\rm{env}}}=\displaystyle \sum _{{\bf{k}},j}{\omega }_{{\bf{k}},j}\,{A}_{{\bf{k}},j}^{\dagger }{A}_{{\bf{k}},j}.\end{eqnarray} \tag{ 24 }$$

The interaction Hamiltonian has the form ${H}_{{\rm{int}}}={\sum }_{j}{s}_{j}{B}_{j}$, mentioned above, i.e. multiple products of system operators s_j with bath operators B_j. We label the system operators ${s}_{1{\bf{k}}}={a}_{{\bf{k}}}^{\dagger }$ and Hermitean conjugate ${s}_{2{\bf{k}}}={a}_{{\bf{k}}}$ with the corresponding bath operators ${B}_{1{\bf{k}}}=\sqrt{2S}{\sum }_{j}g({\omega }_{j})F(k){A}_{{\bf{k}},j}$ and ${B}_{2{\bf{k}}}=\sqrt{2S}{\sum }_{j}{g}^{* }({\omega }_{j}){F}^{* }(k){A}_{{\bf{k}},j}^{\dagger }$, respectively, from which the master equation can be directly derived via Bloch–Redfield approach [13, 30, 31]. This is a convenient version of the equivalent master equation techniques [29, 34, 40, 41].

The Bloch–Redfield master equation in k-space is then directly based on calculating the environmental spectral functions using the bath operators ${\tilde{B}}_{1{\bf{k}}},{\tilde{B}}_{2{\bf{k}}}$ in the interaction picture of system and bath. There are four environmental spectral functions due to the four combinations of the two types of bath operators. However, only two out of the four turn out to be non-zero⁶ . These are the coefficients in the master equation responsible for absorption processes and emission processes. The spectral function which corresponds to phonon absorption is given by

$$\begin{eqnarray}&&{C}_{{\prime} }^{{\rm{abs}}}[-\omega ({\bf{k}})]={\int }_{-\infty }^{\infty }{\rm{d}}\tau {{\rm{e}}}^{-{\rm{i}}\omega ({\bf{k}})\tau } \langle {\widetilde{B}}_{1{\bf{k}}}(\tau ){\widetilde{B}}_{{\prime} }(0) \rangle \end{eqnarray}$$

$$\begin{eqnarray}&&=\,2S{\int }_{-\infty }^{\infty }{\rm{d}}\tau {{\rm{e}}}^{-{\rm{i}}\omega ({\bf{k}})\tau }\displaystyle \sum _{j,j^{\prime} }{g}^{* }({\omega }_{j})g({\omega }_{j^{\prime} })| F(k){| }^{2}\langle {\tilde{A}}_{{\bf{k}},j}^{\dagger }(\tau ){\tilde{A}}_{{\bf{k}}^{\prime} ,j^{\prime} }(0)\rangle \end{eqnarray}$$

$$\begin{eqnarray}&&=\,2S{\int }_{-\infty }^{\infty }{\rm{d}}\tau {{\rm{e}}}^{-{\rm{i}}\omega ({\bf{k}})\tau }\displaystyle \sum _{j,j^{\prime} }{g}^{* }({\omega }_{j})g({\omega }_{j^{\prime} })| F(k){| }^{2}{{\rm{e}}}^{{\rm{i}}{\omega }_{j}\tau }\bar{n}({\omega }_{j}){\delta }_{j,j^{\prime} }{\delta }_{{\bf{kk}}^{\prime} }\end{eqnarray}$$

$$\begin{eqnarray}&&=\,2S| F(k){| }^{2}\displaystyle \sum _{j}| g({\omega }_{j}){| }^{2}\delta [\omega ({\bf{k}})-{\omega }_{j}]\,\bar{n}({\omega }_{j})\end{eqnarray}$$

$$\begin{eqnarray}&&=\,2S| F(k){| }^{2}{ \mathcal J }[\omega ({\bf{k}})]\,\bar{n}[\omega ({\bf{k}})],\end{eqnarray} \tag{ 25 }$$

where $\bar{n}(\omega ):= {[\exp ({\hslash }\omega /{k}_{B}T)-1]}^{-1}$ is the Bose–Einstein distribution and ${ \mathcal J }(\omega )={\sum }_{j}{g}^{2}({\omega }_{j})\delta (\omega -{\omega }_{j})$ is the spectral density. In the above calculation the exponential ${{\rm{e}}}^{{\rm{i}}{\omega }_{j}\tau }$ comes from the interaction picture ${\tilde{A}}_{{\bf{k}},j}^{\dagger }(\tau )={{\rm{e}}}^{{{\rm{i}}{H}}_{{\rm{env}}}\tau }\,{A}_{{\bf{k}},j}^{\dagger }\,{{\rm{e}}}^{-{{\rm{i}}{H}}_{{\rm{env}}}\tau }={{\rm{e}}}^{{\rm{i}}{\omega }_{{\bf{k}},j}\tau }{A}_{{\bf{k}},j}^{\dagger }$ due to the environmental Hamiltonian ${H}_{{\rm{env}}}={\sum }_{{\bf{k}},j}{\omega }_{{\bf{k}},j}\,{A}_{{\bf{k}},j}^{\dagger }{A}_{{\bf{k}},j}$ with an environmental dispersion relation ${\omega }_{{\bf{k}},j}$ which as a function of ${\bf{k}}$ is given by the Fourier transform of the spatial coupling function between different environmental bosons. Consistent with phonons in the harmonic approximation, we assume uncoupled environmental bosons and hence ${\omega }_{{\bf{k}},j}={\omega }_{j}$ and $\bar{n}({\omega }_{{\bf{k}},j})=\bar{n}({\omega }_{j})$. In the last step of equation (25) we simply change the argument of the Bose–Einstein distribution; this is allowed because of the delta-function. Analogously we calculate the other spectral function, which describes emission processes:

$$\begin{eqnarray}&&{C}_{2k,1k^{\prime} }^{\mathrm{em}}[\omega ({\bf{k}})]=2S{\int }_{-\infty }^{\infty }{\rm{d}}\tau {{\rm{e}}}^{{\rm{i}}\omega ({\bf{k}})\tau }\langle {\tilde{B}}_{2{\bf{k}}}(\tau ){\tilde{B}}_{1{\bf{k}}^{\prime} }\rangle \end{eqnarray}$$

$$\begin{eqnarray}&&=\,2S{\int }_{-\infty }^{\infty }{\rm{d}}\tau {{\rm{e}}}^{{\rm{i}}\omega ({\bf{k}})\tau }\displaystyle \sum _{j,j^{\prime} }g({\omega }_{j}){g}^{* }({\omega }_{j^{\prime} })| F(k){| }^{2}\,{{\rm{e}}}^{-{\rm{i}}{\omega }_{j}\tau }\langle {A}_{{\bf{k}},j}{A}_{{\bf{k}}^{\prime} ,j^{\prime} }^{\dagger }\rangle \end{eqnarray}$$

$$\begin{eqnarray}&&=\,2S\,| F(k){| }^{2}{ \mathcal J }[\omega ({\bf{k}})](\bar{n}({\omega }_{j})+1).\end{eqnarray} \tag{ 26 }$$

With this the master equation for the density matrix ρ becomes (note that ${\hslash }=1$ throughout):

$$\begin{eqnarray}\begin{array}{rcl}\dot{\rho } & = & {\rm{i}}[\rho ,{H}_{s}]+2S\displaystyle \sum _{{\bf{k}}}| F(k){| }^{2}{ \mathcal J }[\omega ({\bf{k}})]{\bar{n}}_{{\bf{k}}}\left({a}_{{\bf{k}}}^{\dagger }\rho {a}_{{\bf{k}}}-\displaystyle \frac{1}{2}\{{a}_{{\bf{k}}}{a}_{{\bf{k}}}^{\dagger },\rho \}\right)\\ & & +2S\displaystyle \sum _{{\bf{k}}}| F(k){| }^{2}{ \mathcal J }[\omega ({\bf{k}})]({\bar{n}}_{{\bf{k}}}+1)\left({a}_{{\bf{k}}}\rho {a}_{{\bf{k}}}^{\dagger }-\displaystyle \frac{1}{2}\{{a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}},\rho \}\right),\end{array}\end{eqnarray} \tag{ 27 }$$

where the Bose–Einstein distribution ${\bar{n}}_{{\bf{k}}}$ and the spectral density ${ \mathcal J }(\omega )$ are as defined below equation (25) and $\omega ({\bf{k}})=\gamma B\,+\,2{JS}[-D+{\sum }_{\mu =1}^{3}\cos ({k}_{\mu }{d}_{\mu })]$ is the system's dispersion relation, derived in equation (8).

For solid state environments, the spectral density of the bath is usually parameterised in the continuous limit [34, 35] as:

$$\begin{eqnarray}&&{ \mathcal J }(\omega )=\alpha {\omega }^{s}{\omega }_{c}^{s-1}{{\rm{e}}}^{-\omega /{\omega }_{c}},\end{eqnarray} \tag{ 28 }$$

where α accounts for the strength of the coupling and ${\omega }_{c}$ is the cut-off frequency of the bath (typically it would be the Debye frequency of the material). Typically three cases are distinguished, $s\gt 1$ (super-Ohmic), s = 1 (Ohmic) and $s\lt 1$ (sub-Ohmic).

This master equation is of Lindblad-type and we can identify the coefficient as a single magnon decay rate towards thermal equilibration for each k value:

$$\begin{eqnarray}&&{\gamma }_{{\bf{k}}}=2S\,| F(k){| }^{2}{ \mathcal J }[\omega ({\bf{k}})].\end{eqnarray} \tag{ 29 }$$

We solve equation (27) for the expectation value $\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(t)\rangle$ by using the adjoint master equation in the Heisenberg picture⁷ . The solution is given by⁸ :

$$\begin{eqnarray}&&\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(t)\rangle =\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle \,{{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t}+\bar{n}[\omega ({\bf{k}})](1-{{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t}).\end{eqnarray} \tag{ 30 }$$

This result has several very important implications: we have found that magnons decay to their thermal equilibrium with a single exponential rate ${\gamma }_{{\bf{k}}}$ and we have found an analytical expression for this decay rate for any magnon with a given wave vector ${\bf{k}}$ and any given spatial noise correlation function $f(| {\bf{u}}| )$. This magnon decay rate shows how fast different plane spin waves thermalise in the Brillouin zone, and equally applies to single magnon wave packets big enough to behave like plane waves. This magnon decay rate determines their decoherence. This, together with the magnon speed, determines the coherence range of the magnons and is an important indicator for their usefulness as quantum and classical information carriers. Figure 2 shows the magnon decay rate for spatially uncorrelated noise. We will investigate the magnon decay rate for different temperatures, and noise correlation lengths below.

Zoom In Zoom Out Reset image size

Furthermore the thermal equilibrium distribution of magnons as a function of temperature is given by the long-time limit $\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(t=\infty )\rangle =\bar{n}[\omega ({\bf{k}})]$ which gives the average number of magnons for any wave vector ${\bf{k}}$. It is obtained by inserting the system's dispersion relation $\omega ({\bf{k}})$ (which is itself defined by the Fourier transform of the spatial homogeneous spin–spin coupling) into the Bose–Einstein distribution $\bar{n}(\omega )$.

Beyond the investigation of single magnons, our master equation approach also allows us to calculate the time evolution of the average magnetisation of an entire crystal. We will sum over all sites and obtain the magnetisation of the entire ferromagnet rather than the microscopic expectation value of a single site. The ability to obtain the dynamics of a macroscopic solid-state quantity from calculations of a microscopic quantum master equation shows the versatility of this master-equation approach which can be used beyond the usual applications of nanoscale quantities. The magnetisation is furthermore an experimentally more accessible quantity than the time evolution of single magnons.

The macroscopic magnetisation $\langle {m}_{z}\rangle$ is converted to bosonic operators via the Holstein–Primakoff transform:

$$\begin{eqnarray}&&\langle {m}_{z}\rangle =\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{r}}}\langle {S}_{{\bf{r}}}\rangle =S-\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{r}}}\langle {a}_{{\bf{r}}}^{\dagger }{a}_{{\bf{r}}}\rangle \end{eqnarray}$$

$$\begin{eqnarray}&&=\,S-\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{k}}}\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}\rangle .\end{eqnarray} \tag{ 31 }$$

Accordingly the time-evolution of a non-equilibrium state of the magnetisation is given by:

$$\begin{eqnarray}&&\langle {m}_{z}(t)\rangle =S-\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{k}}}[{{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t}\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle +\bar{n}[\omega ({\bf{k}})](1-{{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t})].\end{eqnarray} \tag{ 32 }$$

This describes the relaxation of the system for example from one (system) temperature to another (bath) temperature or from one magnetic field to another. Assuming that we start from an equilibrium state with some magnetic field ${B}_{0}\ne B$ and then evolve to the equilibrium state with magnetic field B, we have $\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle =\bar{n}[{\omega }_{0}({\bf{k}})]=\bar{n}[\gamma {B}_{0}+2{JS}{\sum }_{\mu =1}^{3}(\cos \,{k}_{\mu }{d}_{\mu }-1)]$

$$\begin{eqnarray}&&\langle {m}_{z}(t)\rangle =S-\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{k}}}\bar{n}[\omega ({\bf{k}})]-\displaystyle \sum _{{\bf{k}}}(\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle -\bar{n}[\omega ({\bf{k}})]){{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t}.\end{eqnarray} \tag{ 33 }$$

The discrete summations are explicitly given by a summation over the Brillouin zone, where each k_μ (where $\mu =x,y,z$) is discretised as ${k}_{\mu ,n}=\tfrac{2\pi }{{N}_{\mu }{d}_{\mu }}n$ with the whole numbers $n\in \left\{-\tfrac{{N}_{\mu }}{2},\,\mathrm{...},\tfrac{{N}_{\mu }}{2}-1\right\}$. For large enough N we can change the summation to an integral over the Brillouin zone volume:

$$\begin{eqnarray}&&\langle {m}_{z}(t)\rangle =S-\displaystyle \frac{1}{N}\displaystyle \sum _{{k}_{x}=-\pi /d}^{\pi /d}\displaystyle \sum _{{k}_{y}=-\pi /d}^{\pi /d}\displaystyle \sum _{{k}_{z}=-\pi /d}^{\pi /d}\{\bar{n}[\omega ({\bf{k}})]+(\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle -\bar{n}[\omega ({\bf{k}})]){{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t}\},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\langle {m}_{z}(t)\rangle \approx S-\displaystyle \frac{{{\rm{d}}}^{3}}{{(2\pi )}^{3}}{\int }_{-\pi /{d}}^{\pi /{d}}{{\rm{d}}{k}}_{x}{\int }_{-\pi /{d}}^{\pi /{d}}{{\rm{d}}{k}}_{y}{\int }_{-\pi /{d}}^{\pi /{d}}{{\rm{d}}{k}}_{z}\,\{\bar{n}[\omega ({\bf{k}})]+(\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle -\bar{n}[\omega ({\bf{k}})]){{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t}\}.\end{eqnarray} \tag{ 35 }$$

As an example we choose as the initial state the equilibrium distribution which would occur at high external magnetic field $\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle =\bar{n}[\omega ({\bf{k}})]$ and regard the decay to a new equilibrium position for low magnetic field. For a saturated, i.e. fully magnetised initial state $\langle {a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}(0)\rangle ={\mathrm{lim}}_{B\to \infty }\bar{n}[\omega ({\bf{k}})]=0$ because the magnetic field B aligns all spins to the ground state of zero magnons present. In other words the magnetic field (by appearing in the dispersion relation, equation (8)) effectively increases the energy offset of the system energies so that there are no excited magnon states populated in the thermal equilibrium of infinite magnetic field. This simplifies the initial magnetisation $\langle {m}_{z}(t)\rangle =S$ and the time-dependent magnetisation:

$$\begin{eqnarray}&&\langle {m}_{z}(t)\rangle \approx S-\displaystyle \frac{{{\rm{d}}}^{3}}{{(2\pi )}^{3}}{\int }_{-\pi /{d}}^{\pi /{d}}{{\rm{d}}}^{3}k\,\bar{n}[\omega ({\bf{k}})](1-{{\rm{e}}}^{-{\gamma }_{{\bf{k}}}t}).\end{eqnarray} \tag{ 36 }$$

The continuous limit for the calculation of the magnon decay rate equation (29) and the magnetisation equation (35) facilitates finding an analytical expression. The change from discrete to continuous is justified by smoothness, and smoothness in k-space is guaranteed by large enough number of spins N in real space as this means 'high resolution' in the Brillouin zone in k-space.

However we need to keep in mind, that the Fourier-transform of the spatial correlation function, equation (16) is a separate discrete calculation. This can only be changed to a continuous Fourier transform equation (18) for a correlation length which is not too short (otherwise the correlation function is not smooth enough to change the sum to an integral) and not too long (otherwise the correlations do not decay over the length of the entire crystal and a finite summation cannot be changed to an infinite Fourier transform). Uncorrelated noise and nearest-neighbour correlations are therefore investigated by discrete Fourier transformation, while functional forms are only investigated as continuous Fourier transforms for correlation lengths not shorter than the spin spacing.

Next we will regard the dynamics in the Ohmic case at different temperatures and for several different cases of spatial correlation functions f(u). In each case we give the corresponding Fourier transform F(k), equation (18), the magnon decay rate ${\gamma }_{{\bf{k}}}$, equation (29), and a numerical plot of the decay of the magnetisation to its equilibrium value, equation (35). Note that most variables are dimensionless in the following and the plots hence make relative statements about the effects of spatial correlations. For those spatial correlation functions with a correlation length ξ this length is given in units of the spacing d between the system's (equally spaced $d={d}_{x}={d}_{y}={d}_{z}$) sites, i.e. $\xi =10$ means a correlation length of 10 spins.

In the Ohmic case the spectral density, equation (28), becomes:

$$\begin{eqnarray}&&{ \mathcal J }[\omega ({\bf{k}})]=\alpha \omega ({\bf{k}}){{\rm{e}}}^{-\omega ({\bf{k}})/{\omega }_{c}}\qquad \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&\mathrm{with}\ \,\,\omega ({\bf{k}})=2{JS}\left[-3+\displaystyle \sum _{\mu =1}^{3}\cos ({k}_{\mu }{d}_{\mu })\right].\end{eqnarray} \tag{ 38 }$$

For a sufficiently high cut-off frequency (${\omega }_{c}=100\,| J|$ in the following calculations), the spectral density as a function of ${\bf{k}}$ is essentially proportional to the system's dispersion relation. Note that a large cut-off frequency ${\omega }_{c}$ does not contradict a weak system-environment coupling $g({\omega }_{j})$ mentioned in section 3.

5.1. Thermal noise in spatially uncorrelated environments

For spatially uncorrelated thermal noise we have:

$$\begin{eqnarray}&&f(u)={\delta }_{{\bf{u}},0}^{(3)}={\gamma }_{0}\,{\delta }_{{u}_{x},0}\,{\delta }_{{u}_{y},0}\,{\delta }_{{u}_{z},0},\,\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&F(k)=1,\,\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\gamma }_{{\bf{k}}}=2S\,| F(k){| }^{2}{ \mathcal J }[\omega ({\bf{k}})]\,\end{eqnarray}$$

$$\begin{eqnarray}&&=\,4{S}^{2}\alpha J\left[-3+\displaystyle \sum _{\mu =1}^{3}\cos ({k}_{\mu }{d}_{\mu })\right]{{\rm{e}}}^{-\omega ({\bf{k}})/{\omega }_{c}}.\,\end{eqnarray} \tag{ 41 }$$

Figure 3 shows the decay of magnetisation with time from a fully magnetised ferromagnet to its thermal equilibrium in the absence of external magnetic field. The speed of the dynamics is governed by the noise intensity, which is set by the energy exchange between system and environment α (i.e. the square of their coupling) as well as J, since in the Ohmic case regarded here, the spectral density scales linearly with frequency $\omega ({\bf{k}})$, which in turn is linear in $| J|$ in the absence of external magnetic field. Thus the time is given in units of $\alpha | J|$.

Zoom In Zoom Out Reset image size

Increasing temperature decreases the equilibrium value of the magnetisation as increased thermal noise disturbs the magnetic order increasingly. For low temperatures the equilibrium value reaches close to saturation magnetisation. The relevant scale of temperature ${k}_{B}T$ is set by the system energies which in the absence of external magnetic field are defined by ${\hslash }| J|$. Apart from changing the equilibrium magnetisation, temperature also changes the decay slightly, which can be seen in figure 3 by rescaling the decay curves at different temperatures to start at the same initial and final values. Increasing temperature leads to an overall slightly faster decay of the magnetisation as it increases the thermally populated phonons, which in turn increases the interactions between system and environment.

5.2. Spatial nearest-neighbour correlations in the noise

An uncorrelated noise environment is only found in the limit where the noise correlation length is far below the lattice constant. This is the case when the phononic noise signal varies randomly on those short length scales. We first investigate a correction to this limit, where the spatially decaying interaction between a phononic site and the spins is long enough to reach not just one spin but also the nearest neighbour spin. We therefore compare spatially uncorrelated noise with an environment where the noise of nearest-neighbours of the cubic lattice has a small correlation. This will develop a first understanding of how the introduction of any spatial correlations in the noise changes the decay characteristics of the magnetisation. We introduce the parameter η, which quantifies the amount of nearest neighbour noise correlation, where $\eta =0$ is the limit of uncorrelated noise:

$$\begin{eqnarray}&&f(u)={\delta }_{{\bf{u}},0}^{(3)}+\eta \,{\delta }_{u,1}^{(3)},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&F({\bf{k}})=1+\eta \displaystyle \sum _{\mu =1}^{3}\cos ({k}_{\mu }{d}_{\mu }).\end{eqnarray} \tag{ 43 }$$

Figure 4(A) shows the change to the magnon decay rate by introducing nearest-neighbour correlations in the noise: The decay rate increases for low k-values, while it reduces for high k-values, i.e. for fast magnons. This can be useful in the context of magnonics, where typically fast magnons are used for information transport and slow decay is desirable. For the change of the magnetisation from one equilibrium state to another the increased decay rate will be the dominant effect since the Bose–Einstein distribution is heavily weighted towards low k-values. Indeed figure 4(B) shows that the magnetisation decays faster to its equilibrium value with increasing nearest-neighbour correlations in the noise than in the uncorrelated case. The magnon decay rates in the Brillouin zone change considerably with the introduction of noise correlations. At lower values the decay rates increase, which is the strongest effect as the thermal distribution is weighted towards ${\bf{k}}=0$. At the edges of the Brillouin zone however spatial noise correlations strongly reduce the decay rates of fast magnons. This leads to a maximum of the decay rate for intermediate ${\bf{k}}$-values.

Zoom In Zoom Out Reset image size

5.3. Gaussian spatial decay of noise correlations

Rather than the simple but unrealistic nearest-neighbour correlations, we now consider Gaussian spatial decay of noise correlations with a correlation length ξ. This also allows us to investigate much stronger spatial correlations, since correlations between all pairs of spins are considered and since a correlation length of $\xi =1.3$ spin spacings already implies that nearest neighbours have a correlation above 0.5. The spatial correlation function is:

$$\begin{eqnarray}&&f(u)=\exp (-| {\bf{u}}{| }^{2}/{\xi }^{2}),\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&F(k)={\pi }^{3/2}{\xi }^{3}\exp \left(-\displaystyle \frac{{k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}}{4}{\xi }^{2}\right).\end{eqnarray} \tag{ 45 }$$

We find that F(k) and hence the magnon decay rate becomes peaked around zero with a peak width inversely proportional to the correlation length ξ.

Figure 5 shows the magnon decay rate for spatial noise correlations with different correlation lengths ξ. Since a strong increase of the decay rate is observed for small k-values we have plotted the y-axis with logarithmic scaling. We can see that the decay rate very close to ${\bf{k}}=0$ increases very strongly and the fact that this same decay rate remains close to zero in the case of only nearest-neighbour correlations (figure 4(A)) seems to be resulting from the nearest-neighbour restriction of noise correlations.

Zoom In Zoom Out Reset image size

We have seen that the introduction of spatial correlations and the increase of the correlation length in the phononic noise of a quantum ferromagnet leads to an increased magnon decay for small wave vectors and a decreased decay for large wave vectors in the Brillouin zone. It is interesting to point out that within the approximations of our model the introduction of spatial noise correlations, can be entirely reflected in k-space only by a change in the decay rate γ_k while different k-modes do not experience correlated noise. This shows that fundamentally noise correlations are dependent on the operator basis and that a Fourier lattice transform can mathematically change a system from experiencing correlated noise to uncorrelated noise. This is mathematically parallel to the fact that the system's spatial couplings (i.e. off-diagonal elements in the Hamiltonian) are only reflected in the dispersion relation of the diagonal k-space Hamiltonian, which does mathematically not show any 'coupling' between different k-modes of magnons.

The effect of both increased and decreased decay rates due to spatial noise correlations in a Heisenberg-model ferromagnet can be considered an extension to the effect of super- and sub-radiance in an atomic gas, first discussed by Dicke [42]. Super- and sub-radiance are phenomena with increasingly rich relevance and applications: sub-radiance due to correlated noise has been shown to enable undisturbed classical information transport in spin chains [21]. Recently super-radiance has been measured in ensembles of nitrogen-vacancy centres [43].

The occurrence of super-radiance around small wave vectors and sub-radiance for large wave vectors can be understood by regarding the single-magnon states which are given by (see equation (7.244) in [25]) $| {\bf{k}}\rangle ={\sum }_{r}{{\rm{e}}}^{{\rm{i}}{k}{r}}| r\rangle$, where $| r\rangle$ represents the state with only the one spin at position r flipped against the external magnetic field and all others aligned. These states are the eigenvectors of the Hamiltonian and correspond to the decay rate γ_k. We can therefore regard these states at small and large k values as corresponding to the occurrence of super- and sub-radiance, respectively, in the ferromagnet. These single magnon states in the limit of k ≈ 0 and $k\approx \pm \pi /d$ become very similar to the prototypical super- and sub-radiant states $| \uparrow \downarrow \rangle +| \downarrow \uparrow \rangle$ and $| \uparrow \downarrow \rangle -| \downarrow \uparrow \rangle$, respectively. The sub-radiance originates from the relative phase of ${{\rm{e}}}^{{\rm{i}}{k}{r}}={{\rm{e}}}^{{\rm{i}}\pi }=-1$ which leads to a cancellation in the environmental coupling. In this context it is also interesting to mention that the introduction of classical noise into a quantum evolution has been associated with sustaining a broader momentum distribution in the time evolution [44].

Figure 6 shows the magnetisation as a function of time for different temperatures and correlation lengths. Since the equilibrium magnetisation $\langle {m}_{z}(t=\infty )\rangle$ is dependent on the temperature as shown in the inset of figure 3(A), we plot the relative magnetisation ${ \mathcal M }$ to focus on the spatial correlation effects and compare these at different temperatures.

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{\langle {m}_{z}(t)\rangle -\langle {m}_{z}(\infty )\rangle }{\langle {m}_{z}(0)\rangle -\langle {m}_{z}(\infty )\rangle }.\end{eqnarray} \tag{ 46 }$$

At low temperatures we find a similar result to nearest-neighbour correlations: increasing correlation length leads to a faster decay of the correlations, i.e. super-radiance is dominating. At higher temperatures however, we start to notice a second effect: there is also a very slowly decaying part, which becomes increasingly important with longer correlation length. Correlated noise pushes the magnetisation decay rate into both extremes: a fast decaying part and a slowly decaying part. In figure 5 we saw that super-radiance occurs around ${\bf{k}}=0$ and sub-radiance occurs for large k-values. At low temperatures, where the distribution of magnons is peaked at zero, it therefore makes sense that this super-radiant effect is dominant. At higher temperatures, where the distribution of magnons is wider the sub-radiant effect becomes more prominant. For longer correlation lengths, where the magnon decay rate is peaked with a smaller width eventually sub-radiance dominates at high temperatures.

Zoom In Zoom Out Reset image size

The antiferromagnet has been discussed in detail in [14]. It can be calculated in a similar way, however an additional Bogoliubov transformation is required after the Holstein–Primakoff transform. The resulting dispersion relations and magnon decay rates are different from the ferromagnetic case. We will show below a short comparison in uncorrelated noise environments and want to point out that this comparison is strongly dependent on the type of spectral density. For Ohmic spectral density the magnon decay rate in the Brillouin zone behaves qualitatively the same. Figure 2 showed the ferromagnetic case and figure 7 compares this to the (in the Ohmic case) similar antiferromagnetic behaviour.

Zoom In Zoom Out Reset image size

6.1. Ohmic versus super-Ohmic spectral densities

Figure 7 compares the decay rates of the super-Ohmic and Ohmic case and allows a comparison to the dispersion relation. For ease of display we have plotted all functions as a function of a diagonal line through the Brillouin zone. While the magnon decay rate as a function of the wave vector ${\bf{k}}$ behaves quite similar in Ohmic spectral density, this is no longer the case for super-Ohmic spectral densities, where the magnon decay rate in the antiferromagnetic case goes to zero for magnons with maximal ${\bf{k}}$. Figure 8 shows this behaviour with two k_z-slices of the Brillouin zone. This strongly reduced decay rate of fast magnons in the antiferromagnet is an interesting result. In the context of quantum information transport a fast magnon with a slow decay rate means that information might be transported robustly. Note that this is true independent of temperature and the number of spins.

Zoom In Zoom Out Reset image size

The cause of the different behaviour lies in the antiferromagnetic dispersion relation, which goes to zero for magnons at maximal ${\bf{k}}=(\pm \pi /{d}_{x},\pm \pi /{d}_{y},\pm \pi /{d}_{z})$. The antiferromagnetic dispersion relation is given by [14]:

$$\begin{eqnarray}&&{\omega }_{\mathrm{AF}}({\bf{k}})=2{J}_{\mathrm{AF}}S\sqrt{{D}^{2}-{\left(\displaystyle \sum _{\mu =1}^{3}\cos {k}_{\mu }{d}_{\mu }\right)}^{2}}.\end{eqnarray} \tag{ 47 }$$

This behaviour translates to the decay rate due to the spectral density in the super-Ohmic case. In the Ohmic case the behaviour of the dispersion relation does not translate to the decay rates because the linear scaling of the dispersion relation is counteracted by an additional factor in the antiferromagnetic magnon decay rate [14] which scales as ${(D+{\sum }_{\mu }{k}_{\mu }{d}_{\mu })}^{-1}$ and leads to finite values of the decay rate for maximal ${\bf{k}}$ in the Ohmic case.

In conclusion, we found an analytical solution for the time-evolution of magnons in a quantum ferromagnet with spatially correlated phononic noise. The magnons decay with a single exponential decay rate to their thermal equilibrium. This remains true for both spatially uncorrelated and correlated phononic noise environments and we have found an analytical expression for this decay rate as a function of wave vector ${\bf{k}}$ and the noise correlation function $f(| {\bf{r}}-{\bf{r}}^{\prime} | )$. The introduction of noise correlations causes faster decay (super-radiance) of slow magnons and slow decay (sub-radiance) of fast magnons.

We also calculated the total macroscopic magnetisation of the entire quantum ferromagnet from the microscopic master equation. The magnetisation decays from a fully magnetised state to its thermal equilibrium faster with higher temperature and shows less remaining magnetisation in the equilibrium state with higher temperature. Introducing nearest-neighbour correlations in the noise leads to an overall slightly faster decay of the magnetisation. However, considering the more general model of spatial correlations with a correlation length ξ and a Gaussian profile we find that spatial correlations not only introduce faster decay of the magnetisation but for stronger correlations and higher temperatures there also appears a slowly-decaying part to the magnetisation. Thus strong spatial noise correlations split the decay of the magnetisation into both extremes: instead of one typical time scale there appears a fast and slow decaying part of the magnetisation.

In the comparison between a ferromagnet and an antiferromagnet we find qualitatively similar behaviour in the magnon decay rate in Ohmic environments. However in super-Ohmic noise environments, i.e. when the noise spectral density scale superlinear with the frequency, the antiferromagnet (in contrast to the ferromagnet) shows a magnon decay rate that goes to zero at the largest $| {\bf{k}}|$ values in the Brillouin zone, which may be an interesting regime for magnonics and quantum information transport.

We would like to thank Andrew Greentree for helpful discussions and support. This work was supported in part by the Australian Government through the Australian Research Council (ARC) under projects DP140100375, FT110100225 and CE170100026. It was also supported by computational resources provided by the Australian Government through the National Computational Infrastructure National Facility. JJ acknowledges support by the German Bundesministerium für Bildung und Forschung BMBF. AR and MAMD acknowledge the Spanish MINECO grant FIS2015-67411, the CAM research consortium QUITEMAD+ S2013/ICE-2801, and US Army Research Office through grant W911NF-14-1-0103 for partial financial support.

Here we show in detail why we simplified the form of the system-environment interaction in the main article as ${S}^{x}X$. We take the full interaction Hamiltonian and show that the resulting master equation is equivalent.

The full interaction Hamiltonian is given by ${H}_{{\rm{int}}}\propto {\boldsymbol{S}}\cdot {\boldsymbol{R}}$, where ${\boldsymbol{S}}=({S}^{x},{S}^{y},{S}^{z})$ is the spin vector and ${\boldsymbol{R}}=(X,Y,Z)$ the position operators of the bosonic environment. Here we assume a local environmental model,

$$\begin{eqnarray}\begin{array}{rcl}{H}_{{\rm{int}}} & = & \displaystyle \sum _{j}\displaystyle \sum _{{\boldsymbol{r}}}g({\omega }_{j})[{S}_{{\boldsymbol{r}}}^{x}({A}_{{\boldsymbol{r}},j}^{x}+{A}_{{\boldsymbol{r}},j}^{x\dagger })\\ & & +{S}_{{\boldsymbol{r}}}^{y}({A}_{{\boldsymbol{r}},j}^{y}+{A}_{{\boldsymbol{r}},j}^{y\dagger })+{S}_{{\boldsymbol{r}}}^{z}({A}_{{\boldsymbol{r}},j}^{z}+{A}_{{\boldsymbol{r}},j}^{z\dagger })],\end{array}\end{eqnarray} \tag{ A.1 }$$

here A and ${A}^{\dagger }$ stand for annihilation and creation operators of the environmental boson modes, and we have assumed that the coupling function $g({\omega }_{j})$ is isotropic and the same for every member of the lattice. On the other hand, the Hamiltonian of the environment is

$$\begin{eqnarray}&&{H}_{{\rm{env}}}=\displaystyle \sum _{j}\displaystyle \sum _{{\boldsymbol{r}}}{\omega }_{j}({A}_{{\boldsymbol{r}},j}^{x\dagger }{A}_{{\boldsymbol{r}},j}^{x}+{A}_{{\boldsymbol{r}},j}^{y\dagger }{A}_{{\boldsymbol{r}},j}^{y}+{A}_{{\boldsymbol{r}},j}^{z\dagger }{A}_{{\boldsymbol{r}},j}^{z}),\end{eqnarray} \tag{ A.2 }$$

which is written as

$$\begin{eqnarray}&&{H}_{{\rm{env}}}=\displaystyle \sum _{j}\displaystyle \sum _{{\boldsymbol{k}}}{\omega }_{j}({A}_{{\boldsymbol{k}},j}^{x\dagger }{A}_{{\boldsymbol{k}},j}^{x}+{A}_{{\boldsymbol{k}},j}^{y\dagger }{A}_{{\boldsymbol{k}},j}^{y}+{A}_{{\boldsymbol{k}},j}^{z\dagger }{A}_{{\boldsymbol{k}},j}^{z}),\end{eqnarray} \tag{ A.3 }$$

after taking Fourier transform.

After the Holstein–Primakoff transformation with linear spin-wave approximation, the interaction term reads

$$\begin{eqnarray}\begin{array}{rcl}{H}_{{\rm{int}}} & = & \displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\boldsymbol{r}}\in A}[\sqrt{2S}({a}_{{\boldsymbol{r}}}+{a}_{{\boldsymbol{r}}}^{\dagger })({A}_{{\boldsymbol{r}},j}^{x}+{A}_{{\boldsymbol{r}},j}^{x\dagger })\\ & & -{\rm{i}}\sqrt{2S}({a}_{{\boldsymbol{r}}}-{a}_{{\boldsymbol{r}}}^{\dagger })({A}_{{\boldsymbol{r}},j}^{y}+{A}_{{\boldsymbol{r}},j}^{y\dagger })\\ & & +(S-{a}_{{\boldsymbol{r}}}^{\dagger }{a}_{{\boldsymbol{r}}})({A}_{{\boldsymbol{r}},j}^{z}+{A}_{{\boldsymbol{r}},j}^{z\dagger })].\end{array}\end{eqnarray} \tag{ A.4 }$$

We now can consider several simplifications of this Hamiltonian based on the following facts,

• The term ${a}_{{\boldsymbol{r}}}^{\dagger }{a}_{{\boldsymbol{r}}}$ is negligible in comparison to the others in the regime where the spin-wave theory is valid: $\langle {a}_{{\boldsymbol{r}}}^{\dagger }{a}_{{\boldsymbol{r}}}\rangle \ll 2S$.
• We ignore the terms $S({A}_{{\boldsymbol{r}},j}+{A}_{{\boldsymbol{r}},j}^{\dagger })$ because they are fast oscillators in a rotating-wave approximation which we may neglect in the weak coupling limit.

Therefore, we arrive at

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}\displaystyle \sum _{{\boldsymbol{r}}}g({\omega }_{j})[{a}_{{\boldsymbol{r}}}({A}_{{\boldsymbol{r}},j}^{x}+{A}_{{\boldsymbol{r}},j}^{x\dagger }-{\rm{i}}{A}_{{\boldsymbol{r}},j}^{y}-{\rm{i}}{A}_{{\boldsymbol{r}},j}^{y\dagger })+{a}_{{\boldsymbol{r}}}^{\dagger }({A}_{{\boldsymbol{r}},j}^{x}+{A}_{{\boldsymbol{r}},j}^{x\dagger }+{\rm{i}}{A}_{{\boldsymbol{r}},j}^{y}+{\rm{i}}{A}_{{\boldsymbol{r}},j}^{y\dagger })],\end{eqnarray} \tag{ A.5 }$$

and by taking the Fourier transform this interaction becomes

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}\displaystyle \sum _{{\boldsymbol{k}}}g({\omega }_{j})[{a}_{{\boldsymbol{k}}}({A}_{-{\boldsymbol{k}},j}^{x}+{A}_{{\boldsymbol{k}},j}^{x\dagger }-{\rm{i}}{A}_{{\boldsymbol{-k}},j}^{y}-{\rm{i}}{A}_{{\boldsymbol{k}},j}^{y\dagger })+{a}_{{\boldsymbol{k}}}^{\dagger }({A}_{{\boldsymbol{k}},j}^{x}+{A}_{{\boldsymbol{-k}},j}^{x\dagger }+{\rm{i}}{A}_{{\boldsymbol{k}},j}^{y}+{\rm{i}}{A}_{{\boldsymbol{-k}},j}^{x\dagger })],\end{eqnarray} \tag{ A.6 }$$

Because of the Riemann–Lebesgue lemma [45], for small coupling $g({\omega }_{j})$ we can safely neglect the counter-rotating terms in (A.6) and arrive at:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}\displaystyle \sum _{{\boldsymbol{k}}}g({\omega }_{j})[{a}_{{\boldsymbol{k}}}({A}_{{\boldsymbol{k}},j}^{x\dagger }-{\rm{i}}{A}_{{\boldsymbol{k}},j}^{y\dagger })+{a}_{{\boldsymbol{k}}}^{\dagger }({A}_{{\boldsymbol{k}},j}^{x}+{\rm{i}}{A}_{{\boldsymbol{k}},j}^{y})].\end{eqnarray} \tag{ A.7 }$$

Now the problem becomes equivalent to a collection of uncoupled harmonic oscillators given by their operators ${a}_{{\boldsymbol{k}}}$, which are coupled to two sets of independent environments (x and y) characterised by ${\boldsymbol{k}}$.

We assume that the environments are in the Gibbs state at some temperature T,

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{E} & = & {Z}^{-1}{{\rm{e}}}^{-\beta {H}_{E}}\\ & = & {Z}^{-1}{{\rm{e}}}^{-\beta \displaystyle \sum _{j}\displaystyle \sum _{{\boldsymbol{k}}}{\omega }_{j}({A}_{{\boldsymbol{k}},j}^{x\dagger }{A}_{{\boldsymbol{k}},j}^{x}+{A}_{{\boldsymbol{k}},j}^{y\dagger }{A}_{{\boldsymbol{k}},j}^{y}+{A}_{{\boldsymbol{k}},j}^{z\dagger }{A}_{{\boldsymbol{k}},j}^{z})},\end{array}\end{eqnarray} \tag{ A.8 }$$

where $Z=\mathrm{Tr}({{\rm{e}}}^{-\beta {H}_{E}})$ is the partition function with $\beta =1/{k}_{{\rm{B}}}T$. From now on we shall use natural units ${\hslash }={k}_{{\rm{B}}}=1$.

The standard tools to obtain a master equation for a weak interaction with the environment can be found in [45–48]. The bath correlation functions split in two different types, depending on which kind of process are associated to them, emission or absorbtion. We have

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\mathrm{abs}}(\omega ({\boldsymbol{q}}),{\boldsymbol{k}},{\boldsymbol{q}}) & = & 2S{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{u}\displaystyle \sum _{i,j}{{\rm{e}}}^{-{\rm{i}}[\omega ({\boldsymbol{q}})-{\omega }_{i}]u}g({\omega }_{i})g({\omega }_{j}){\rm{Tr}}[{\rho }_{E}({A}_{{\boldsymbol{k}},i}^{x\dagger }-{\rm{i}}{A}_{{\boldsymbol{k}},i}^{y\dagger })({A}_{{\boldsymbol{q}},j}^{x}+{\rm{i}}{A}_{{\boldsymbol{q}},j}^{y})]\\ & = & 2S{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{u}\displaystyle \sum _{j}{{\rm{e}}}^{-{\rm{i}}[\omega ({\boldsymbol{q}})-{\omega }_{j}]u}{g}^{2}({\omega }_{j}){\rm{Tr}}[{\rho }_{E}({A}_{{\boldsymbol{k}},j}^{x\dagger }{A}_{{\boldsymbol{q}},j}^{x}+{A}_{{\boldsymbol{k}},j}^{y\dagger }{A}_{{\boldsymbol{q}},j}^{y})]\end{array}\end{eqnarray}$$

$$\begin{eqnarray}&&=\,4S{\int }_{-\infty }^{\infty }{\rm{d}}{u}\displaystyle \sum _{j}{{\rm{e}}}^{-{\rm{i}}[\omega ({\boldsymbol{k}})-{\omega }_{j}]u}{g}^{2}({\omega }_{j}){\bar{n}}_{{\boldsymbol{k}}}({\omega }_{j}){\delta }_{{\boldsymbol{k}},{\boldsymbol{q}}},\,\end{eqnarray} \tag{ A.9 }$$

where

$$\begin{eqnarray}&&{\bar{n}}_{{\boldsymbol{k}}}({\omega }_{j}):= {[\exp ({\omega }_{j}/T)-1]}^{-1}\end{eqnarray} \tag{ A.10 }$$

is the number of quanta in the bath labelled by '${\boldsymbol{k}}$' at frequency ${\omega }_{j}$. In the continuous limit, we introduce the spectral density function, ${ \mathcal J }(\omega ):= {\sum }_{j}{g}^{2}(\omega )\delta (\omega -{\omega }_{j})$ (formally), and change the sum by an integral

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\mathrm{abs}}(\omega ({\boldsymbol{k}}),{\boldsymbol{k}},{\boldsymbol{q}}) & = & S\displaystyle \int {\rm{d}}\omega {\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{u}{{\rm{e}}}^{-{\rm{i}}[\omega ({\boldsymbol{k}})-\omega ]u}{ \mathcal J }(\omega ){\bar{n}}_{{\boldsymbol{k}}}(\omega ){\delta }_{{\boldsymbol{k}},{\boldsymbol{q}}}\\ & = & 2\pi S\displaystyle \int {\rm{d}}\omega \delta [\omega ({\boldsymbol{k}})-\omega ]{ \mathcal J }(\omega ){\bar{n}}_{{\boldsymbol{k}}}(\omega ){\delta }_{{\boldsymbol{k}},{\boldsymbol{q}}}\\ & = & 2\pi S{ \mathcal J }[\omega ({\boldsymbol{k}})]{\bar{n}}_{{\boldsymbol{k}}}{\delta }_{{\boldsymbol{k}},{\boldsymbol{q}}},\end{array}\end{eqnarray} \tag{ A.11 }$$

where ${\bar{n}}_{{\boldsymbol{k}}}:= {[\exp (\omega ({\boldsymbol{k}})/T)-1]}^{-1}.$

Similarly we obtain that

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\mathrm{em}}(\omega ,{\boldsymbol{k}},{\boldsymbol{q}}) & = & \displaystyle \frac{S}{2}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{u}\displaystyle \sum _{i,j}{{\rm{e}}}^{{\rm{i}}[\omega ({\boldsymbol{q}})-{\omega }_{i}]u}g({\omega }_{i})g({\omega }_{j})\mathrm{Tr}[{\rho }_{E}({A}_{{\boldsymbol{k}},i}^{x}+{\rm{i}}{A}_{{\boldsymbol{k}},i}^{y})({A}_{{\boldsymbol{q}},j}^{x\dagger }-{\rm{i}}{A}_{{\boldsymbol{q}},j}^{y\dagger })]\\ & = & 2\pi S{ \mathcal J }[\omega ({\boldsymbol{k}})]({\bar{n}}_{{\boldsymbol{k}}}+1){\delta }_{{\boldsymbol{k}},{\boldsymbol{q}}}.\end{array}\end{eqnarray} \tag{ A.12 }$$

We can see that the inclusion of the ${S}^{y}Y$ term in the spin-boson coupling of system and environment has indeed simply led to an additional noise term analogous to the one from the ${S}^{x}X$ term. The ${S}^{z}Z$ term has led to negligible terms in the linear spin-wave regime.

After changing the summation index from ${\bf{r}}^{\prime}$ to ${\bf{u}}$ the interaction Hamiltonian reads:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\sqrt{2S}\displaystyle \sum _{j}g({\omega }_{j})\displaystyle \sum _{{\bf{k}},{\bf{k}}^{\prime} }\displaystyle \sum _{{\bf{u}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}^{\prime} \cdot {\bf{u}}}f(| {\bf{u}}| ){a}_{{\bf{k}}}^{\dagger }{A}_{{\bf{k}}^{\prime} ,j}\displaystyle \sum _{{\bf{r}}}\displaystyle \frac{1}{N}{{\rm{e}}}^{{\rm{i}}({\bf{k}}^{\prime} -{\bf{k}})\cdot {\bf{r}}}+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ B.1 }$$

We then identify the last summation over ${\bf{r}}$ as a Kronecker-delta in each dimension by inserting ${k}_{\mu }=\tfrac{2\pi {n}_{\mu }}{{N}_{\mu }{d}_{\mu }}$ and ${r}_{\mu }={d}_{\mu }{\tilde{n}}_{\mu }$ for each component:

$$\begin{eqnarray}&&\displaystyle \sum _{{\tilde{n}}_{\mu }=1}^{{N}_{\mu }}\displaystyle \frac{1}{{N}_{\mu }}\exp \left[{\rm{i}}\displaystyle \frac{2\pi }{{N}_{\mu }}({n}_{\mu }^{{\prime} }-{n}_{\mu }){\tilde{n}}_{\mu }\right]={\delta }_{{n}_{\mu }^{{\prime} },{n}_{\mu }}\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{\bf{r}}}\displaystyle \frac{1}{N}{{\rm{e}}}^{{\rm{i}}({\bf{k}}^{\prime} -{\bf{k}})\cdot {\bf{r}}}={\delta }_{{k}_{x}^{{\prime} },{k}_{x}}{\delta }_{{k}_{y}^{{\prime} },{k}_{y}}{\delta }_{{k}_{z}^{{\prime} },{k}_{z}}.\end{eqnarray} \tag{ B.3 }$$

The Hamiltonian then becomes equation (15) from the main article after performing the summation over ${\bf{k}}^{\prime}$.

Since the function $f(| u| )$ only depends on the magnitude $| {\bf{u}}|$, we can simplify this three-dimensional Fourier transform to a one-dimensional integral using polar coordinates. The result is also sometimes called a 3D Hankel transform or 3D Fourier-Bessel transform. For this integration we choose without loss of generality the polar direction to coincide with the ${\bf{k}}$ direction, such that ${\bf{k}}\cdot {\bf{u}}={ku}\cos (\theta )$ and then simplify by substitution $v=\cos \theta$:

$$\begin{eqnarray}&&\begin{array}{}F({\bf{k}})=\displaystyle \frac{1}{{V}_{d}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{0}^{\infty }f(u){{\rm{e}}}^{{\rm{i}}{k}{u}\cos \theta }{u}^{2}\sin \theta \,{\rm{d}}{u}\,{\rm{d}}\theta \,{\rm{d}}\phi \end{array}\end{eqnarray}$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{4\pi }{{V}_{d}k}{\int }_{0}^{\infty }f(u)u\sin ({ku})\,{\rm{d}}{u}.\,\end{eqnarray} \tag{ C }$$