Decay rates and energies of free magnons and bound states in dissipative XXZ chains

TL;DR

This paper investigates how nonlocal dissipation affects the energies and decay rates of free magnons and bound states in dissipative XXZ chains, providing analytical solutions and proposing an experimental realization with Rydberg atoms.

Contribution

It analytically derives bound state solutions in dissipative XXZ chains with nonlocal dissipation and explores their potential for experimental realization.

Findings

Bound state energies can lie outside the two-magnon continuum.

Decay rates of bound states are confined within the two-magnon decay continuum.

Nonlocal dissipation enables more flexible engineering of bound state dispersions.

Abstract

Chains of coupled two-level atoms behave as 1D quantum spin systems, exhibiting free magnons and magnon bound states. While these excitations are well studied for closed systems, little consideration has been given to how they are altered by the presence of an environment. This will be especially important in systems that exhibit nonlocal dissipation, e.g. systems in which the magnons decay due to optical emission. In this work, we consider free magnon excitations and two-magnon bound states in an XXZ chain with nonlocal dissipation. We prove that whilst the energy of the bound state can lie outside the two-magnon continuum of energies, the decay rate of the bound state has to always lie within the two-magnon continuum of decay rates. We then derive analytically the bound state solutions for a system with nearest-neighbour and next-nearest-neighbour XY interaction and nonlocal dissipation, finding that the inclusion of nonlocal dissipation allows more freedom in engineering the energy and decay rate dispersions for the bound states. Finally, we numerically study a model of an experimental set-up that should allow the realisation of dissipative bound states by using Rydberg-dressed atoms coupled to a photonic crystal waveguide (PCW). We demonstrate that this model can exhibit many key features of our simpler models.

Equations

$$\begin{array}{r}\dot{\hat{\rho }}(t)=\frac{i}{\hslash }\left[\hat{\rho }(t),\hat{H}\right]+\sum _{i,l}^N\frac{{\Gamma }_{il}}{2}\left(2{\hat{\sigma }}_i^{-}\hat{\rho }(t){\hat{\sigma }}_l^{+}-\left\{{\hat{\sigma }}_l^{+}{\hat{\sigma }}_i^{-},\hat{\rho }(t)\right\}\right),\end{array}$$

$$\hat{H}=\hslash \Delta \sum _i^N{\sigma }_i^z+\sum _{i\ne l}^N\hslash V_{il}{\hat{\sigma }}_i^{+}{\hat{\sigma }}_l^{-}+\frac{\hslash J_z}{2}\sum _i^N{\hat{\sigma }}_i^z{\hat{\sigma }}_{i+1}^z.$$

$$\begin{array}{rl}& \frac{dG(ij,t)}{dt}-{\delta }_{ij}\delta (t)=-i\Delta G(ij,t)-\frac{\Gamma }{2}G(ij,t)+\\ & 4i{J}_{z}G(ij,t)-i\sum _{p\ne j}^N\left(V_{pj}-i\frac{{\Gamma }_{pj}}{2}\right)G(ip,t),\end{array}$$

$$G(k,\omega )=\underset{ϵ\to 0}{\lim }i(\omega -E(k)+iϵ{)}^{-1},$$

$$E(k)=-4J_z+\Delta -\frac{i\Gamma }{2}+\sum _{l=1}^N\left(2V_{l0}-i{\Gamma }_{l0}\right)\cos (kl)$$

$$\begin{array}{rl}& \frac{dG(ij,lm;t)}{dt}-(1-{\delta }_{ij})\delta (t)({\delta }_{il}{\delta }_{jm}+{\delta }_{im}{\delta }_{jl})=\\ & \left(-2i\Delta +8i{J}_z-\Gamma -4i{J}_z{\delta }_{m,l+1}\right)G(ij,lm;t)\\ & -i\sum _{p\ne l}^N{J}_{pl}G(ij,pm;t)-i\sum _{p\ne m}^N{J}_{pm}G(ij,pl;t)\\ & +2i{\delta }_{lm}\sum _{p\ne m}^N{J}_{pm}G(ij,pm;t),\end{array}$$

$$\begin{array}{rl}& G(r,r^{\prime },Q,\Omega )=\Gamma (r,r^{\prime },Q,\Omega )h(r)\\ & -\sum _{r^{\prime \prime }}^{N}K(r,r^{\prime \prime },Q,\Omega )G(r^{\prime \prime },r^{\prime },Q,\Omega ).\end{array}$$

$$\begin{array}{rl}& K(r,r^{\prime };Q,\Omega )=\frac{2i}{N}\sum _{q\in BZ}\frac{\cos (q{r}^{\prime })}{\Omega -S(q,Q)}\text{crossproduct}\\ & \left[4i{J}_z\cos (q)-2i\left(V(r)-i\frac{\Gamma (r)}{2}\right)\cos (Qr/2)\right],\\ & \Gamma (r,r^{\prime };Q,\Omega )=-\frac{2i}{N}\sum _{q\in BZ}\frac{\cos (q{r}^{\prime })\cos (qr)}{\Omega -S(q,Q)}.\end{array}$$

$$\begin{array}{rl}& S(q,Q)=E(Q/2+q)+E(Q/2-q)\\ & =-8J_z+2\Delta -i\Gamma +\sum _{j=1}^N\left(4V_{j0}-2i{\Gamma }_{j0}\right)\cos (jQa/2)\cos (jqa),\end{array}$$

$$\det \left[\delta (r,r^{\prime \prime })+K(r,r^{\prime \prime },Q,\Omega )\right]=0.$$

$$\begin{array}{rl}& \left(1-\frac{1}{N}\sum _{q\in BZ}\frac{8J_z{\cos }^2(qa)}{\Omega -S(q,Q)}\right)\left(1+\frac{1}{N}\sum _{q^{\prime }\in BZ}\frac{S(q^{\prime },Q)-t}{\Omega -S(q^{\prime },Q)}\right)\\ & +\frac{8J_z}{N^2}\sum _{q,q^{\prime }\in BZ}\frac{\cos (qa)\cos (q^{\prime }a)(S(q^{\prime },Q)-t)}{\left[\Omega -S(q,Q)\right]\left[\Omega -S(q^{\prime },Q)\right]}=0,\end{array}$$

$$\begin{array}{r}(\Omega +t)\left[\frac{I_0(\Omega ,Q)}{8J_z}-I_0(\Omega ,Q)I_2(\Omega ,Q)+I_1(\Omega ,Q{)}^2\right]=0,\end{array}$$

$$I_m(\Omega ,Q)={\int }_{-\pi }^{\pi }\frac{{\cos }^m(q)}{\Omega -S(q,Q)}\frac{dq}{2\pi }.$$

$$\begin{array}{r}\dot{\hat{\rho }}(t)=i\left[\hat{\rho }(t),\hat{H}\right]+\sum _k\left(2{\hat{J}}_k^{-}\hat{\rho }(t){\hat{J}}_k^{+}-\left\{{\hat{J}}_k^{+}{\hat{J}}_k^{-},\hat{\rho }(t)\right\}\right).\end{array}$$

$$\begin{array}{rl}& |Q⟩=\sum _{ij}^N{\alpha }_Q{f}_Q(|r_i-r_j|)e^{iQ(r_i+r_j)/2}{\hat{\sigma }}_i^{+}{\hat{\sigma }}_j^{+}|0⟩,\end{array}$$

$$\begin{array}{rl}& \frac{d{\rho }_Q(t)}{dt}\approx -\sum _{k}8{\gamma }_k|{\alpha }_{Q}F(Q/2-k){|}^2{\rho }_Q(t),\end{array}$$

$$F(Q/2-k)=\sum _{r\ne 0}{f}_Q(r)e^{ir(Q/2-k)}$$

$$\sum _{k}2|{\alpha }_{Q}F(Q/2-k){|}^2=1,$$

$${\tilde{\gamma }}_Q=\sum _{k}2{\gamma }_k|{\alpha }_{Q}F(Q/2-k){|}^2<{\gamma }_{\text{min}},$$

$$\sum _{k}2({\gamma }_k-{\gamma }_{\text{min}})|{\alpha }_{Q}F(Q/2-k){|}^2<0.$$

$$\sum _{k}2({\gamma }_k-{\gamma }_{\text{max}})|{\alpha }_{Q}F(Q/2-k){|}^2>0,$$

$$\begin{array}{rl}& \Re [E(k)]=-4J_z+\Delta +2V_{12}\cos (ka)\\ & |\text{imaginary}[E(k)]|=\frac{\Gamma }{2}+{\Gamma }_{12}\cos (ka)\\ & \Re [S(q,Q)]=-8J_z+2\Delta +4V_{12}\cos (Qa/2)\cos (qa)\\ & |\text{imaginary}[S(q,Q)]|=\Gamma +2{\Gamma }_{12}\cos (Qa/2)\cos (qa).\end{array}$$

$$\begin{array}{r}\Omega (Q)=-4J_z+2\Delta -i\Gamma +\frac{(2V_{12}-i{\Gamma }_{12}{)}^2}{4J_z}{\cos }^2(Qa/2),\end{array}$$

$$\begin{array}{rl}& \Re [\Omega (Q)]=-4J_z+2\Delta +\frac{4V_{12}^2-{\Gamma }_{12}^2}{4J_z}{\cos }^2(Qa/2)\\ & |\text{imaginary}[\Omega (Q)]|=\Gamma +\frac{V_{12}{\Gamma }_{12}}{J_z}{\cos }^2(Qa/2).\end{array}$$

$$\begin{array}{rl}& \Re [\Omega (Q)]=-4J_z+2\Delta \\ & |\text{imaginary}[\Omega (Q)]|=\Gamma \mp 2J_z{\cos }^2(Qa/2),\end{array}$$

$$\begin{array}{rl}\Re [E(k)]=& -4J_z+\Delta +2V_{12}\cos (ka)+2V_{13}\cos (2ka)\\ |\text{imaginary}[E(k)]|=& \frac{\Gamma }{2}+{\Gamma }_{12}\cos (ka)+{\Gamma }_{13}\cos (2ka)\\ \Re [S(q,Q)]=& -8J_z+2\Delta +4V_{12}\cos (Qa/2)\cos (qa)\\ & +4V_{13}\cos (Qa)\cos (2qa)\\ |\text{imaginary}[S(q,Q)]|=& \Gamma +2{\Gamma }_{12}\cos (Qa/2)\cos (qa)\\ & +2{\Gamma }_{13}\cos (Qa)\cos (2qa).\end{array}$$

$$\begin{array}{rl}& \Omega (Q)=-8J_z+2\Delta -i\Gamma +4J_{13}\cos (Qa)+\frac{J_{12}^2{\cos }^2(Qa/2)}{J_z}\\ & +\frac{J_{12}^2{\cos }^2(Qa/2)J_{13}\cos (Qa)}{2J_z^2}+\frac{8J_z^2}{2J_z+J_{13}\cos (Qa)},\end{array}$$

$$\begin{array}{rl}& \Re [\Omega (Q)]=-8J_z+2\Delta +\frac{4V_{12}^2-{\Gamma }_{12}^2}{4J_z}{\cos }^2(Qa/2)\\ & +\frac{V_{13}(4V_{12}^2-{\Gamma }_{12}^2)-2{\Gamma }_{13}{\Gamma }_{12}{V}_{12}}{8J_z^2}\cos (Qa){\cos }^2(Qa/2)\\ & +4V_{13}\cos (Qa)+\frac{16J_z^2(4J_z+2V_{13}\cos (Qa))}{(4J_z+2V_{13}\cos (Qa){)}^2+({\Gamma }_{13}\cos (Qa){)}^2}\\ & |\text{imaginary}[\Omega (Q)]|=\Gamma +\frac{V_{12}{\Gamma }_{12}}{J_z}{\cos }^2(Qa/2)+2{\Gamma }_{13}\cos (Qa)\\ & +\frac{{\Gamma }_{13}(4V_{12}^2-{\Gamma }_{12}^2)+8V_{13}{\Gamma }_{12}{V}_{12}}{16J_z^2}\cos (Qa){\cos }^2(Qa/2)\\ & -\frac{16J_z^2{\Gamma }_{13}\cos (Qa)}{(4J_z+2V_{13}\cos (Qa){)}^2+({\Gamma }_{13}\cos (Qa){)}^2}.\end{array}$$

$$\begin{array}{rl}\Re [\Omega (Q)]& \approx -4J_z+2\Delta +2V_{13}\cos (Qa)\\ |\text{imaginary}[\Omega (Q)]|& \approx \Gamma +{\Gamma }_{13}\cos (Qa),\end{array}$$

$$A_{ij}=\frac{J_{xy}{e}^{iK|r_{ij}|}}{2\sqrt{1-(\delta /(2J)+i{\gamma }_c/(4J){)}^2}}.$$

Full text

I Introduction

One very interesting direction of recent research on ultracold atomic or molecular gases involves the study of the collective quantum dynamics of internal excitations of the atoms (or molecules) positioned in ordered arrays. Such systems behave as strongly coupled two-level quantum systems (i.e. spin-1/2 systems), and can explore fundamental issues in the quantum dynamics of many-body systems subject to strong interparticle interactions Yan et al. (2013); Labuhn et al. (2016); Fukuhara et al. (2013)

A famous example of a strong-interaction phenomenon in quantum spins systems is provided by magnon bound states, first proposed by Bethe Bethe (1931) more than 80 years ago. In this work, it was shown that magnon bound states could form in 1D spin-1/2 Heisenberg chain with nearest-neighbour interactions, lowering their energy compared to free magnons in the system. Subsequent work then extended this result to higher dimensions, anisotropic spin chains and arbitrary spin including solitons Wortis (1963); Haldane (1982); Southern et al. (1994); Schneider (1981) and spin chains with long-range interactions Torrance and Tinkham (1969); Majumdar (1969); Ono et al. (1971); Letscher and Petrosyan (2018). Furthermore, magnon bound states have been studied in systems with frustration Kecke et al. (2007), topological structure Qin et al. (2017, 2018) and in Floquet systems Agarwala and Sen (2017); Kudo et al. (2009). They have also recently been observed experimentally Fukuhara et al. (2013) and shown to have an important role in magnetisation switching Barker et al. (2013), transport Krimphoff et al. (2017); Ganahl et al. (2012) and to have interesting effects on entanglement entropy Mölter et al. (2014).

One key aspect in all of these studies is that the system is closed and so the question of bound state decay rates is not considered. However, if the system is coupled to an external environment, then the excitations will eventually decay and so it is natural to ask how long lived these excitations can be. For a system with local dissipation, the decay rate of both free excitations and bound states will be given by $m$ times the local decay rate Longo and Evers (2014a) where $m$ is the number of excitations. However, for systems involving radiative decay, the dissipation typically becomes nonlocal, where a range of decay rates to the environment exist, which are either superradiant (greater than the local decay rate) or subradiant (smaller than the local decay rate). In these scenarios, the relative decay rates of the free excitations and bound states becomes unclear. For example, is it possible for the decay rate of the bound states to be smaller than that of the free magnons?

In this work, we address the question of bound state decay rates in systems with nonlocal dissipation. We look at three models with a nearest-neighbour Ising interaction, which is crucial for the bound states to form, and different forms of XY interaction and nonlocal dissipation. The first two models are a nearest-neighbour and next-nearest-neighbour XY interaction for which we can obtain analytical results. The final model is an experimentally achievable setting in which to observe our results with Rydberg dressed atoms coupled to a photonic crystal waveguide.

The layout of the paper is as follows. In Sec II, we derive the general equations needed to obtain the energy and decay rate of the free excitations and bound states. In Sec III, we show that in general the decay rate of the bound state lies within the two-magnon decay rate continuum. Then in Sec IV, we obtain the energies and decay rates for the three models. In Sec V we discuss our results and experimental implementation before drawing conclusions in Sec VI.

II Model

We consider a macroscopic number, $N$, of two-level systems fixed in position on a 1D optical lattice with spacing, $a$, and periodic boundary conditions. The atoms interact with an electromagnetic field which acts as an environment for the system. We assume the Markovian and Born approximations, which are valid provided the coupling between the system and environment is weak. These allow us to describe the system using a master equation approach. We will later discuss the validity of this approximation in relation to our results. The resultant master equation is given by

[TABLE]

where the square brackets represent a commutator and curly brackets represent the anti-commutator. The spin operators are defined as $\hat{\sigma}^{z}_{i}=\ket{e_{i}}\bra{e_{i}}-\ket{g_{i}}\bra{g_{i}}$, $\hat{\sigma}^{-}_{i}=\ket{g_{i}}\bra{e_{i}}$ and $\hat{\sigma}^{+}_{i}=\ket{e_{i}}\bra{g_{i}}$, where $\ket{e_{i}}$ and $\ket{g_{i}}$ are the excited and ground states of the atom respectively. We require that the eigenvalues of the matrix $\Gamma_{il}$ are all greater than or equal to zero, in order for Eq. (1) to describe decay of the excited state, driven by the operators $\hat{\sigma}^{-}_{i}$. Then the steady state density matrix is given by $\rho_{ss}=\ket{0}\bra{0}$ where $\ket{0}=\prod_{i}^{N}\ket{g_{i}}$. The Hamiltonian is given by

[TABLE]

For the rest of the paper, we will work in units of $\hbar=1$. The Hamiltonian in Eq. (2) conserves the number of excitations in the system whilst the dissipator allows the excitations to decay. We can therefore talk about the dynamics of few-magnon excitations. To compute the energies and decay rates of one- and two-magnon excitations in our system, we employ a Green’s function method.

We first start with the single magnon Green’s function, defined as $G(ij,t)=\Tr[\hat{\sigma}^{-}_{i}(t)\hat{\sigma}^{+}_{j}(0)\hat{\rho}(0)]\Theta(t)=\bra{0}\hat{\sigma}^{-}_{i}(t)\hat{\sigma}^{+}_{j}(0)\ket{0}\Theta(t)$, choosing the initial condition, $\rho(0)$, to be the pure state $\ket{0}\bra{0}$. The single magnon Green’s function obeys the following equation

[TABLE]

where $\Gamma=\Gamma_{ii}$. Fourier transforming Eq. (3) gives the spectrum of the single magnon states from the poles of

[TABLE]

where

[TABLE]

is the single magnon dispersion, with the real part corresponding to the energy and the magnitude of the imaginary part corresponding to the decay rate.

For two magnons, we consider the Green’s function $G(ij,lm;t)=\Tr[\hat{\sigma}^{-}_{i}(t)\hat{\sigma}^{-}_{j}(t)\hat{\sigma}^{+}_{l}\hat{\sigma}^{+}_{m}\rho(0)]\Theta(t)$, which obeys the following equation

[TABLE]

where $J_{pl}=V_{pl}-i\Gamma_{pl}/2$. This equation can be rewritten as a matrix equation and partially Fourier transformed with $G(r,r^{\prime},Q,\Omega)=\sum_{R}e^{-iRQ}\int_{-\infty}^{\infty}G(ij,lm,t)e^{i\Omega t}dt$ to give (see Appendix A)

[TABLE]

where $r=r_{i}-r_{j}$, $r^{\prime}=r_{l}-r_{m}$ and

[TABLE]

The momenta $q$ and $Q$ in Eq. (7) and Eq. (8) are the difference and sum of momenta, defined by $q=(k_{1}-k_{2})/2$ and $Q=k_{1}+k_{2}$, where $k_{1}$ and $k_{2}$ are the momenta of the individual magnons. The momenta $q$ are summed over the Brillouin zone denoted by BZ. The function in the denominator of Eq. (8), $S(q,Q)$, is the dispersion of two free magnons, given by

[TABLE]

which determines the poles of $\Gamma(r,r^{\prime},Q,\Omega)$, whilst the two-magnon bound states are given by solutions to the determinant equation

[TABLE]

Because of the nearest-neighbour Ising coupling, this determinant equation can be simplified to (see Appendix B)

[TABLE]

where $t=8J_{z}-2\Delta+i\Gamma$. In the limit $N\rightarrow\infty$, we can rewrite Eq. (11) as

[TABLE]

where

[TABLE]

In Section IV, we shall find the energies and decay rates of the bound states by solving Eq. (12) (or Eq. (11) where appropriate) for three specific forms of the XY interaction and nonlocal dissipation: a nearest-neighbour model, next-nearest-neighbour model and a photonic crystal waveguide model. Note that $\Omega=t=-8J_{z}+2\Delta-i\Gamma$ is always a solution to Eq. (12). However, this solution always lies within the two-magnon energy continuum. In general, we will dismiss any solutions that lie inside the two-magnon energy continuum where the bound state is no longer well defined because it can scatter into the continuum states and become a resonance. While it is possible to have bound states that exist in the scattering continuum Hsu et al. (2016), these usually occur when the system has certain symmetries that protect the state, which we are not aware of existing in our models.

III General Decay Rates of Bound States

We first show that in general, for any model with nonlocal dissipation of the form given in the master equation, Eq. (1), the decay rate of the bound state always lies within the two-magnon decay rate continuum, i.e. the bound state cannot decay more quickly or slowly than its constituent parts. To show this, we consider Eq. (1) rewritten in diagonal form

[TABLE]

Here, $\hat{J}_{k}$ is a decay operator for mode $k$, given by $\hat{J}_{k}^{-}=\sqrt{\gamma_{k}}\sum_{i}^{N}c_{i}^{k}\hat{\sigma}^{-}_{i}$, where $c_{i}^{k}$ is the $i^{th}$ component of the $k^{th}$ eigenvector of $\Gamma_{il}/2$ and $\gamma_{k}$ is the corresponding eigenvalue. For a periodic or large enough system, the eigenvector components are given by $c_{i}^{k}=e^{ikr_{i}}/\sqrt{N}$. To determine the decay rate of the bound state, we focus on the initial dynamics of the pure state density matrix, $\hat{\rho}(0)=\ket{Q}\bra{Q}$ where $\ket{Q}$ is the wavefunction of a bound state with momentum $Q$, given by

[TABLE]

where $f_{Q}(r)$ is some localised function that determines the spatial decay of the bound state, with $r=|r_{i}-r_{j}|$, and $\alpha_{Q}$ is a normalisation constant given by $\alpha_{Q}=1/\left(2N\sum_{r\neq 0}|f_{Q}(r)|^{2}\right)$. The equation of motion for a pure bound state density matrix at short initial times is given by

[TABLE]

where $\rho_{Q}(t)=\bra{Q}\hat{\rho}(t)\ket{Q}$ and

[TABLE]

is the Fourier transform of the localised function. At later times, there can be the population of coherences between the bound state and scattering states, which we have neglected. We can see that the bound state density matrix has a decay rate of $4\tilde{\gamma}_{Q}$, where $\tilde{\gamma}_{Q}\equiv\sum_{k}2\gamma_{k}|\alpha_{Q}F(Q/2-k)|^{2}$, which is weighted sum of all single magnon decay rates. Note that $\tilde{\gamma}_{Q}$, is equivalent to the decay rate we will obtain from our Green’s function method.

For local dissipation where $\gamma_{k}=\gamma\equiv\Gamma/2$, the sum over $k$ in $\tilde{\gamma}_{Q}$ can be completed to give

[TABLE]

and so the decay rate of the bound state wavefunction (which is half the decay rate of the pure density matrix) is $2\gamma$ as expected. For nonlocal dissipation, in order to have a bound state decay rate that exists below the two-magnon decay rate continuum, we would need

[TABLE]

where $\gamma_{\text{min}}$ is the smallest decay rate for a single magnon. However, using Eq. (18), we can rewrite this condition as

[TABLE]

Both $|\alpha_{Q}F(Q/2-k)|^{2}$ and $\gamma_{k}-\gamma_{\text{min}}$ are always positive, which means this condition can never be fulfilled. The lowest decay rate that could possibly be achieved for the bound state is the lowest decay rate that can be achieved for two free magnons, although this may not always obey the bound state equation. The same argument applies for showing that the bound state cannot have a decay rate above the two-magnon decay rate continuum, such that

[TABLE]

where $\gamma_{\text{max}}$ is the largest decay rate in the system. Again $|\alpha_{Q}F(Q/2-k)|^{2}>0$, but $\gamma_{k}-\gamma_{\text{max}}<0$, so this condition can not be satisfied and the bound state decay rate must always lie within the two-magnon decay rate continuum.

IV Results

IV.1 Nearest-Neighbour Model

Having now shown in general that the decay rate of the bound state always lies within the two-magnon decay rate continuum, we now look at three specific models for dissipative bound states. The first model we consider is where all interactions are nearest-neighbour (NN). The energies and decay rates of the one and two free magnon states are given by

[TABLE]

Solving Eq. (12) gives the following bound state solution (see Appendix C)

[TABLE]

which can be written in terms of the energy and decay rate as

[TABLE]

These expressions first appeared in Ref. Longo and Evers (2014b), although we analyse them in more detail here. For the expressions in Eqs. (24), there are limits to the parameters we can choose for the solutions to satisfy the bound state equation, Eq. (12). However, provided we choose $V_{12}$ and $\Gamma_{12}$ such that the energy term in Eq. (24) lies below the two-magnon energy continuum, then we find the bound state equation is always satisfied. We also have to impose $\Gamma/2\geq|\Gamma_{12}|$ in order for the dissipator to always give decay.

Comparing the bound state solution Eq. (24) to the free magnon dispersions in Eq. (22), we see the energy and decay rate of the bound state depend on a mixture of the interaction and dissipation. The presence of nonlocal dissipation creates a negative shift in energy compared to the XY interaction, which means that the bound state energy is shifted further from the two-magnon energy continuum than in a closed system. This is important as the effects of nonlocal dissipation will not only cause the bound state to decay, but will alter its dynamics travelling through the lattice meaning that even if the bound state has a very small decay rate, it is not sufficient to ignore environmental effects. Furthermore, due to nonlocal dissipation, there is more freedom to engineer the bound state energy and decay than in a closed system. For example, the bound state energy band can be made entirely flat by choosing $V_{12}=\Gamma_{12}/2$. Also, by choosing $V_{12}=0$ such that there is no XY interaction, the bound state experiences only local dissipation, with a decay rate of $\Gamma$, whereas the one and two free magnons still experience nonlocal dissipation. Finally, in the limit where $V_{12},\Gamma_{12}\ll J_{z}$, the effects of the XY interaction and nonlocal dissipation become negligible, with the energy of the bound state tending to $-4J_{z}$ and the decay rate tending to $\Gamma$ which would be expected for an Ising model with local dissipation.

The relative signs of the XY interaction, nonlocal dissipation and Ising interaction allow the bound state decay rate to be tuned such that it is either entirely subradiant or superradiant, with the most super- or subradiant decay at $Qa=0$ and a decay rate of $\Gamma$ at the band edge, $Qa=\pm\pi$. To find how subradiant or superradiant it is possible to make the bound state, we extremise the decay rate of the bound state with respect to the parameters $V_{12}$ and $\Gamma_{12}$, while still obeying the constraint that the bound state energy must lie below the two-magnon energy continuum. We also maintain a fixed decay rate $\Gamma$ (otherwise there is always a trivial minimal decay rate with $\Gamma=\Gamma_{12}=0$). We find the extremal decay rates and corresponding energies are given by

[TABLE]

where the negative sign gives the maximal (minimal) decay rate and the positive sign gives the minimal (maximal) decay rate for $J_{z}<0$ ($J_{z}>0$). The largest values for $\Gamma_{12}$ and $V_{12}$ occur when the bound state makes contact with the energy continuum at $Qa=0$. In figure 1, we show the minimal decay rate solution for $J_{z}<0$ and $\Gamma=2|\Gamma_{12}|$.

The bound state decay rate lies in the two-magnon decay rate continuum as expected and is smaller than half the free magnon decay rates at $Qa=\pi$ and $2/3$ of the continuum at $Qa=0$, with the lowest energy bands of the two-magnon continuum having the smallest decay rates. For the maximal decay rate solution, the results are similar to Figure 1, but the decay rates reverse, with the lowest energy bands having the highest decay rates and the bound state solution having a larger decay rate than most of the two-magnon decay rate continuum.

IV.2 Next-Nearest-Neighbour Model

The NN model studied in the previous section demonstrated many features of dissipative bound states, but also missed some qualitative features of bound states with longer range hopping. We therefore consider a next-nearest-neighbour (NNN) model, finding that the inclusion of additional site interactions produces important differences in the properties of the bound state compared to a NN model. The one and two free magnon energies and decay rates are given by

[TABLE]

The bound state solution is given by (see Appendix D)

[TABLE]

where $J_{12}=V_{12}-i\Gamma_{12}/2$ and $J_{13}=V_{13}-i\Gamma_{13}/2$. Writing in terms of the energy and decay rate gives

[TABLE]

As for the NN model, there is a constraint on the values of the dissipative couplings to ensure the magnons always decay, which is $\Gamma/2\geq|\Gamma_{12}+\Gamma_{13}|$. Likewise, we have to choose parameters that satisfy the bound state condition Eq. (12), finding again that provided the energy of the bound state lies below the continuum, then Eq. (12) is satisfied. Our NNN bound state solution is the same as that found in Ref. Letscher and Petrosyan (2018) but with a complex XY interaction. This is also true of our NN result in Eq. (24), which can be obtained by taking the bound state result in Ref. Wortis (1963) with a complex XY interaction.

The inclusion of an additional site in the XY interaction and nonlocal dissipation results in a more complex bound state solution than in the NN model. Looking at the terms in Eq. (28) in more detail, we see that the NN solution in Eq. (24) can be recovered by letting $V_{13},\Gamma_{13}=0$, and that now we have additional terms due to two-site hopping processes and a term that mixes the NN and NNN parameters. Because of the new magnon hopping terms, the decay rate of the bound state is no longer fixed to be $\Gamma$ at $Qa=\pm\pi$ as was the case for NN interactions, and the smallest and largest decay rates do not have to occur at $Qa=0$ anymore. Therefore the inclusion of NNN interactions allows more freedom in choosing at what momenta $Q$ the bound state can have its highest or smallest decay rate. However, we can now no longer engineer an entirely flat energy band due to the presence of both $\cos(Qa/2)$ and $\cos(Qa)$ terms (unless trivially the NNN couplings are set to zero). Looking at the limit of $V_{13},\Gamma_{13},V_{12},\Gamma_{12}\ll J_{z}$, we find Eq. (28) simplifies to

[TABLE]

We find that there is now always a contribution to the decay rate from the NNN interactions, that means even tightly confined bound states still experience the effects of nonlocal dissipation, which was not the case for the NN model. We can also see that the smallest decay rate will occur at $Qa=0$ ($Qa=\pm\pi$) and largest decay rate at $Qa=\pm\pi$ ($Qa=0$) for $\Gamma_{13}<0$ ($\Gamma_{13}>0$).

We now extremise the NNN bound state decay rate for a fixed $\Gamma$ with respect to the parameters $V_{12}$, $V_{13}$, $\Gamma_{12}$ and $\Gamma_{13}$ to find the smallest and largest decay rates the bound state can have while its energy remains separate from the two-magnon energy continuum. Due to the complexity of Eqs. (28), we solve this numerically, finding that the solution with minimal (maximal) decay rate occurs when $V_{12}=\pm 1.135J_{z}$, $V_{13}=-0.293J_{z}$, $\Gamma_{12}=\pm 1.926J_{z}$ and $\Gamma_{13}=0.578J_{z}$, and the maximal (minimal) solution occurs when $V_{12}=\mp 1.135J_{z}$, $V_{13}=-0.293J_{z}$, $\Gamma_{12}=\pm 1.926J_{z}$ and $\Gamma_{13}=-0.578J_{z}$ for $J_{z}<0$ ($J_{z}>0$), where in both cases, we are free to choose the positive or negative sign. The largest values of all parameters occur when the bound state energy makes contact with the two-magnon energy continuum at $Qa=0$ as was the case for the NN interactions. In figure 2, we show the minimal solution with $J_{z}<0$ and $\Gamma=2|(\Gamma_{12}+\Gamma_{13})|$.

Again, we find the decay rate of the bound state lies within the two-magnon decay rate continuum, with the bound state having a smaller decay rate than $30\%$ of the continuum at $Qa=\pi$ and up to $70\%$ of the continuum at $Qa=0$. We should note there is a second minimal (maximal) decay rate solution with parameters $V_{12}=\Gamma_{12}=0$, $\Gamma_{13}=+0.402J_{z}$ and $V_{13}=-0.827J_{z}$ and maximal (minimal) solution for $V_{12}=\Gamma_{12}=0$, $\Gamma_{13}=-0.402J_{z}$ and $V_{13}=-0.827J_{z}$ for $J_{z}<0$ ($J_{z}>0$). However, we have not shown this solution as it is more unphysical due to the absence of the NN terms.

IV.3 Photonic Crystal Waveguide Model

We now study one final model which should be an experimentally realisable set-up to study dissipative bound states. We consider Rydberg dressed two-level atoms that are coupled to a photonic crystal waveguide (PCW). Systems of two-level atoms where one state is a Rydberg state or Rydberg dressed are already well studied as realisable quantum simulators Schauss (2018); Weimer et al. (2010); Whitlock et al. (2017); Nguyen et al. (2018); Glaetzle et al. (2015). Likewise, PCWs are also gaining attention as a method for quantum simulation and quantum information processing due to the high tunability of the interactions between coupled quantum emitters Hartmann (2016); Hood et al. (2016); Goban et al. (2014); Douglas et al. (2015); González-Tudela et al. (2015). For atoms coupled to a PCW, photons emitted from the atoms can propagate to other atoms along the chain, which mediates an effective XY interaction and nonlocal dissipation. For a single mode in a dissipative PCW, the XY interaction and nonlocal dissipation are given by Calajó et al. (2016) $V_{il}=\Im{A_{ij}}$ and $\Gamma_{il}=2\Re{A_{ij}}$, where $A_{ij}$ is of the form

[TABLE]

The parameter $J_{xy}$ is the coupling of the atoms to the PCW, $J$ is an energy scale determining the PCW bandwidth, and $K_{wq}a=\pi-\arccos{(\delta/(2J)+i\gamma_{c}/(4J))}=k_{wg}a+i\kappa_{wg}a$ is the PCW wavevector. The PCW wavevector depends on the detuning, $\delta=(\omega_{eg}-\omega_{wg})$, of the atomic transition frequency, $\omega_{eg}$, from the photon mode frequency, $\omega_{wg}$, and also the loss rate of photons from the PCW, $\gamma_{c}$. If $|\delta/J|<2$, then the photon lies within the bandwidth and can propagate along the PCW with a group velocity given by the denominator of Eq. (30), $2\sqrt{1-(\delta/(2J)+i\gamma_{c}/(4J))^{2}}$. However, if $|\delta/J|>2$, then the photon cannot propagate and instead exponentially decays along the PCW.

In order for bound states to form, we also need an Ising interaction. This can be engineered by dressing Glaetzle et al. (2015) either the excited state, $\ket{e}$ or ground state, $\ket{g}$, of an atom with a Rydberg state $\ket{r}$, giving a new state $\ket{\tilde{e}}=\ket{e}+\beta\ket{r}$ where $\beta=\Omega_{d}/2\Delta_{d}$, set by the drive $\Omega_{d}$ and detuning $\Delta_{d}$ that couple $\ket{e}$ to $\ket{r}$. The atoms then interact with an Ising interaction of the form

[TABLE]

where $U_{0}=\hbar\Omega_{d}^{4}/8\Delta_{d}^{3}$ and $R_{c}$ is some cut off length to the interaction. For small $R_{c}$, this is a good approximation to a NN Ising interaction. The sign and magnitude of $U_{0}$ can be fixed by the laser detuning and it is also possible to add additional XY interactions between the atoms which gives more freedom in tuning $V_{ij}$ separately from $\Gamma_{ij}$.

For the PCW system, the one and two free magnon energies and decay rates are given by

[TABLE]

where $\Delta=V_{11}/2+\delta/2+\delta_{\text{add}}$, with $\delta_{\text{add}}$ being an additional detuning to those from the waveguide, and

[TABLE]

For the rest of this section, we will choose the additional detuning, $\delta_{\text{add}}$ such that $\Delta=0$ and so we can ignore the contributions to energy from the onsite term, $V_{11}$ and detuning from the waveguide mode $\delta$. We will also work with $J_{z}<0$.

In figure 3, we plot the energy and decay rate of the single magnon dispersion for $\gamma_{c}/J=2$, $\delta/J=0$ and $J_{xy}/|J_{z}|=3$. If $|\delta/J|<2$ and $\gamma_{c}/J$ is small, then about the points $k=\pm k_{wg}$, the decay rate is well modelled by two Lorentzians with a width of $4\sinh(\kappa_{wg}a/2)$ and maximum value of $\Gamma/[4\tanh(\kappa_{wg}a/2)]$. Similarly, the energy of the magnon is well described by the derivative of a Lorentzian with width $4\sinh(\kappa_{wg}a/2)$ and maximal (minimal) values given by $\pm\Gamma/[8\sinh(\kappa_{wg}a/2)]$. As $\gamma_{c}/J$ decreases (and so $\kappa_{wg}\rightarrow 0$), the energies of the magnons and decay rates about $k=\pm k_{wg}$ diverge within the photonic bandwidth ($|\delta/J|<2$). However, outside the bandwidth ($|\delta/J|>2$), the energy of the magnon is bounded and its decay rate drops to zero as $\gamma_{c}\rightarrow 0$, leaving the system effectively closed. The single magnon dispersions can be thought of as the hybridisation of a photon propagating through the waveguide with a dispersion $\omega_{k}=\omega_{wg}-J\cos(k)$ and momentum $k$, and a single atom with energy $\omega_{eg}$.

We now look at the bound state solutions in the PCW and discuss their properties. The bound state condition, Eq. (12), is too complex to be solved analytically, so we instead tackle the problem numerically for finite sized systems by solving Eq. (11). In Figure 4, we plot some typical solutions of Eq. (11) for a system size of $N=99$, with $\gamma_{c}/J=2$, $J_{xy}/|J_{z}|=3$ and for $\delta/J=(-3,-1.5,0,1.5,3)$.

We see that bound state decay rate lies within the two-magnon decay rate continuum as expected, and is smaller than the decay rate of the lowest energy bands of the continuum for $\delta/J<-2$, but larger than the decay rate of the lowest energy bands of the continuum for $\delta/J>2$. For intermediate detunings, whether the bound state decay rate is smaller or larger than the lowest energy bands depends on the momentum of the bound state. As for the NNN model, we find the minimal and maximal decay rate of the bound state is no longer constrained to occur at $Qa=0$ and that the decay rate at $Qa=\pi$ is not given by $\Gamma$ as a consequence of the long-range interactions. If $\kappa_{wg}$ is large enough, then the bound state solutions are well modelled by the NNN analytics due to the exponential decay of the PCW interaction. This can be seen by the close agreement between the NNN and PCW bound state solutions when $\delta/J=\pm 3$, which gives the largest $\kappa_{wg}$. For intermediate detunings, the agreement is not as good, but can be made increasingly better for larger $\gamma_{c}/J$.

In figure 5, we plot the momentum for which the bound state has the smallest decay rate as a function of $\delta/J$ and $\gamma_{c}/J$. We find that there is a transition between the bound state having the smallest decay rate at $Qa=0$ when $|\delta/J|<1.4$ to $Qa=\pi$ when $|\delta/J|>1.4$. This transition can be explained by looking at the weak XY limit of the NNN bound state solutions given by Eq. (29). In the weak limit, we find that the momentum where the decay rate of the bound state is smallest transitions from $Qa=0$ to $Qa=\pi$ when $\Gamma_{13}$ changes sign. We show when $\Gamma_{13}=0$ in figure 5 by the red dashed lines, and find it agrees well with the transition in the PCW, with $\Gamma_{13}<0$ when $|\delta/J|>1.4$. The transition moves to larger values of $|\delta/J|$ as $\gamma_{c}/J$ increases, and also becomes sharper as the NNN solution becomes a better approximation to the PCW results.

Finally, we discuss how the bound state formation depends on $\delta/J$ and $\gamma_{c}/J$. Figure 6 shows where the bound state rejoins the two magnon energy continuum as a function of $\delta/J$ and $\gamma_{c}/J$.

We find there is a region inside the bandwidth that extends along the $\gamma_{c}/J$ axis where the bound state joins the continuum and that, as $J_{xy}/|J_{z}|$ increases, this region also increases in size. The reason the bound state starts to rejoin the continuum for small $\gamma_{c}/J$ inside the bandwith is due to the diverging strength of the single magnon energy around $k=\pm k_{wg}$. For increasingly large systems, more momentum modes around these points are allowed and so the energy range of the two magnon continuum grows until the bound state is absorbed. However, outside the bandwidth and in the small $\gamma_{c}/J$ limit, the bound state energy can remain separate from the two-magnon energy continuum for any value of $J_{xy}/|J_{z}|$ provided $\delta/J$ is large enough. This is because the two-magnon energy continuum is now bounded as $\gamma_{c}/J\rightarrow 0$ and so bound states can remain separate from the continuum. As mentioned in our discussion of the single magnon dispersion, the imaginary part of the PCW interaction, Eq. (30), becomes negligible in this limit, and so the system becomes closed, with the decay rate of the bound state dropping to zero. When $\gamma_{c}/J$ becomes large, or when $|\delta/J|\gg 2$, the XY interaction becomes increasingly shorter ranged due to the exponential decay, until eventually it is negligible compared to the Ising interaction. In this limit, the bound state is well separated from the two-magnon energy continuum with the bound state energy tending to $-4J_{z}$ and the decay rate tending to $\Gamma$.

Our analysis of a PCW has shown how many features of dissipative bound states can be obtained for a single photonic mode and how, for large $\kappa_{wg}a$, the PCW is well described by the NNN analytics. For a single mode, it is not possible to obtain the NN results, no matter how large $\kappa_{wg}a$ is. To see why this is the case, we look at the NNN bound state solution in Eq. (27). We can see that for an exponentially decaying function, $J_{13}\sim J_{12}^{2}/J_{z}$, which means that there is always a NNN contribution to the bound state solution that is of the order of the NN parts, so the NNN contribution cannot be ignored. However, it could be possible to engineer more exotic XY interactions by combining many modes or coupling to more than one waveguide. This could also be done in parallel with different Rydberg dressing schemes or allowing other interactions, such as dipole interactions, to occur between atoms.

V Discussion

We have shown that two-magnon bound states can generally form in dissipative spin chains with XY and Ising interactions. We find the inclusion of nonlocal dissipation not only gives the bound state a momentum dependent decay rate, but also alters the bound state energy compared to a closed system or system with local dissipation. Nonlocal dissipation also allows for a greater degree of freedom in engineering the energy and decay rate of the bound state. We have shown that the decay rate of the bound state cannot be smaller or larger than its constituent free magnons. Nevertheless, it is still possible to achieve bound states that have a decay rate much lower than a large proportion of the two-magnon decay rate continuum.

We now discuss the experimental set-up of the PCW in more detail. To engineer the bound states, we need to choose an appropriate scheme for Rydberg dressing for the atoms. Rydberg dressing has already been achieved experimentally Zeiher et al. (2016) with $\prescript{87}{}{\text{Rb}}$ atoms, taking the Rubidium hyperfine states $\ket{g}=\ket{1,-1}$ and $\ket{e}=\ket{2,-2}$ and dressing with a suitable Rydberg state of $\ket{r}=\prescript{31}{}{P_{1/2}}$. Therefore, it should be possible to engineer suitable Ising-like interactions with NN or even beyond NN range. The PCW itself can be realised with a SiO alligator waveguide Goban et al. (2014); Douglas et al. (2015) with high tunability over the allowed modes and loss processes. Previous experiments with cold atoms in waveguides have used Caesium, but it should be possible to engineer a waveguide suitable for Rubidium Perrella et al. (2018). When studying the bound states, one has to be careful not to violate the Markovian approximation. For the Markovian approximation to be valid, it is required that the time for a photon to travel down the PCW is negligible compared to the decay rate of the atoms Calajó et al. (2016). This gives the condition

[TABLE]

which is satisfied provided the coupling of the atoms to the waveguide, $J_{xy}$, is weak and also that the detuning is away from the band edge at $\delta=\pm 2J$ when $\gamma_{c}/J$ is small. The expression Eq. (34) also shows that the system needs to be finite to not violate the Markovian approximation. However, we have checked and found that there are bound state solutions with similar properties to those in the main text for finite size systems with open boundary conditions. Therefore, it should be possible to observe many of our bound states results for large enough finite sized systems with open boundary conditions or periodic boundary conditions.

Finally, measurement of the bound state decay rate and energy should be possible by observing the emission when the bound state decays. Following the steps outlined in Ref. Longo and Evers (2014b), the emission properties of the bound state are given by the correlator $g(t,\boldsymbol{r})=\langle\boldsymbol{\hat{E}}^{(-)}(t,\boldsymbol{r})\boldsymbol{\hat{E}}^{(+)}(t,\boldsymbol{r})\rangle$ which can be calculated from the electric field, $\boldsymbol{\hat{E}}^{(-)}(t,\boldsymbol{r})$. For decay of a pure bound state, $\hat{\rho}(0)=\ket{Q}\bra{Q}$, the correlator $g(t,\boldsymbol{r})$ is given by

[TABLE]

where $t_{r}\equiv t-r/c$, $\Delta_{0}^{k}=\text{Re}[E(k)]$, $\Delta_{k}^{Q}=\text{Re}[\Omega(Q)]-\text{Re}[E(k)]$, $\eta=\omega_{eg}^{2}/(4\pi\epsilon_{0}c)$ and $W(\boldsymbol{r})=\boldsymbol{d}/r-\boldsymbol{r}(\boldsymbol{d}.\boldsymbol{r})/r^{3}$ is the far-field dipole emission profile. There are two contributions to the emission of the bound state; one from the decay of the bound state to a single magnon with momentum $k$, and one from the decay of a single magnon to the ground state. The delta functions determine the emission angle $\beta$ for each of these decay processes in terms of the momentum and energy of the bound state and single magnons, where $\beta$ is defined from the perpendicular axis from the spin chain. The total emission is then a sum over all these processes. The quantity $|\alpha_{Q}F(Q/2-k)|^{2}$ that determined the decay rate of the bound state also plays a crucial role in the angular dependence of the emission, which was noted in Longo and Evers (2014b). By examining the spatial and temporal emission of the bound state, it should be possible to determine its energy and decay rate for a given momentum $Q$.

In future work, it would be interesting to extend our results to $m$ magnon-bound states and to see how the decay rates of different magnon sectors compare to one another. Given our proof that the two-magnon bound state decay rate must lie within the continuum of decay rates, it seems likely that this would also be true for $m$ magnon states, and possibly also true for magnon states with larger spin and in systems of higher dimension. It would also be interesting to study different forms of dissipators and find systems where the bound state can have a decay rate that lies outside the two-magnon continuum.

VI Conclusions

We have studied the energies and decay rates of one and two free magnons and two-magnon bound states in an XXZ model with nonlocal dissipation. We have proved that in general the decay rate of the bound state must lie within the decay rate continuum of two free magnons. We have then examined three examples of dissipative bound states in more detail, first looking at two forms of the XY interaction analytically; a nearest-neighbour model and next-nearest-neighbour model. We have found that the inclusion of nonlocal dissipation leads to momentum dependent decay rates and changes in the energy of the bound state compared to a closed system or a system with local dissipation. The nonlocal dissipation also allows a higher degree of tunability in the energies and decay rates of the bound states. Finally, in our third example, we have numerically studied an experimentally realisable model to observe dissipative bound states using Rydberg dressed atoms coupled to a photonic crystal waveguide, which demonstrates many key features of our simpler models and can also be used to obtain our next-nearest-neighbour results within certain parameter regimes.

VII Acknowledgements

This work was supported by EPSRC Grant Nos. EP/K030094/1 and EP/P009565/1 and by the Simons Foundation. Statement of compliance with EPSRC policy framework on research data: All data accompanying this publication are directly available within the publication.

Appendix A Deriving the Bound State Determinant Equation

Below, we outline the steps to obtain the bound state equation in Eq. (11). For an open quantum system, provided the Liouvillian operator is time independent, any Heisenberg operator will obey the adjoint master equation, given by Breuer and Petruccione (2007)

[TABLE]

Therefore, the Green’s function $\Tr(\hat{A}(t)\hat{B(0)}\hat{\rho}(0))=\bra{0}\hat{A}(t)\hat{B}\ket{0}$, with the initial condition $\hat{\rho}(0)=\ket{0}\bra{0}$, will obey

[TABLE]

For the two-magnon Green’s function, $G(ij,lm;t)=\bra{0}\hat{\sigma}^{-}_{i}(t)\hat{\sigma}^{-}_{j}(t)\hat{\sigma}^{+}_{l}\hat{\sigma}^{+}_{m}\ket{0}\Theta(t)$, this gives

[TABLE]

where $J_{pl}=V_{pl}-i\Gamma_{pl}/2$. In order to solve Eq. (38), it will be useful to view it as a matrix equation Majlis (2000) given by $(\boldsymbol{\mathcal{L}}+\boldsymbol{\delta\mathcal{L}})\boldsymbol{G}=\boldsymbol{\mu}\boldsymbol{h}$, where the matrices are defined as

[TABLE]

To solve Eq. (38), we now follow the same steps taken by Wortis Wortis (1963) by introducing the function $\Gamma(ij,lm;t)=G(il;t)G(jm;t)+G(im;t)G(jl;t)$, where $G(jl;t)$ is the single magnon Green’s function. We find that $\Gamma(ij,lm;t)$ obeys Eq. (38) without the last two terms and no $1-\delta_{ij}$ term. Viewed in terms of matrices, this means $\boldsymbol{\mathcal{L}}\boldsymbol{\Gamma}=\boldsymbol{\mu}$ and so we can write $\boldsymbol{\mathcal{L}}=\boldsymbol{\mu}\boldsymbol{\Gamma}^{-1}$. This allows Eq. (38) to be rewritten as

[TABLE]

where in the last line we have defined

[TABLE]

In order to obtain the bound state solutions, we now need to partially Fourier transform Eq. (40). The Fourier transform of $\Gamma(ij,ab;t)$ is given by

[TABLE]

By using the definition of $\Gamma(ij,ab;t)$ and the Fourier transform of the single magnon Green’s function, this can be written as

[TABLE]

where $r_{ia}=r_{i}-r_{a}$. We now rewrite the momentum sums using the sum and difference of momenta, $Q=k_{1}+k_{2}$ and $q=(k_{1}-k_{2})/2$, and also the sum and difference of coordinates $R=(r_{i}+r_{j})/2$, $r=r_{i}-r_{j}$ and $R^{\prime}=(r_{a}+r_{b})/2$, $r^{\prime}=r_{a}-r_{b}$. Once we evaluate the frequency integrals, we then obtain

[TABLE]

where $S(q,Q)$ is the two free magnon dispersion, defined in Eq. (9) in the main text. Similarly, we can Fourier transform and rewrite $K(lm,pq;t)$ as

[TABLE]

where $R^{\prime\prime}=(r_{p}+r_{v})/2$ and $r^{\prime\prime}=r_{p}-r_{v}$. Transforming Eq. (40) by inserting the results of Eq. (44) and Eq. (45) gives

[TABLE]

This equation is obeyed provided we set the integrand to zero such that

[TABLE]

The bound state solutions are found when the determinant of the matrix $\delta(r^{\prime},r^{\prime\prime})+K(r^{\prime},r^{\prime\prime},Q,\Omega)$ is singular, which means $G(r,r^{\prime\prime},Q,\Omega)$ cannot be written as the sum of two free magnon solutions. The bound state solutions are therefore solutions to

[TABLE]

If the Ising interaction is nearest-neighbour such that $U_{il}=J_{z}\delta_{l,i+1}$, we can simplify the determinant in Eq. (48) to obtain Eq. (11) in the main text.

Appendix B Simplifying the Determinant Condition

We first define the Ising and XY matrices,

[TABLE]

where

[TABLE]

This allows us to rewrite the determinant condition, Eq. (48), as

[TABLE]

The determinant can be simplified by subtracting the last column from all the other columns, $\text{C}_{1}-\text{C}_{N}$, $\text{C}_{2}-\text{C}_{N}$,…$\text{C}_{N-1}-\text{C}_{N}$, giving

[TABLE]

where we partially Laplace expand the determinant. For the first determinant, we can swap the first and last column, $C_{1}\leftrightarrow C_{N}$ and then swap the first and last row, $R_{1}\leftrightarrow R_{N}$. In the second determinant, we can carry out the row-swap operation, $R_{N}\leftrightarrow R_{N-1}$, followed by $R_{N-1}\leftrightarrow R_{N-2}$, $R_{N-2}\leftrightarrow R_{N-3}$ etc. until the last row becomes the first row. This then gives

[TABLE]

which are the determinants of arrowhead matrices, where an arrowhead matrix is a matrix of the form

[TABLE]

Using the Sherman-Morrison-Woodbury formula, we can evaluate the determinant of the arrowhead matrix by rewriting Eq. (54) as

[TABLE]

where

[TABLE]

Using this gives a determinant of

[TABLE]

Substituting the values of $a$, $b_{i}$, $c_{i}$ and $d_{i}$ for the two arrowhead matrices in Eq. (53), we obtain the determinant equation

[TABLE]

Once we plug in the definitions of $\boldsymbol{ZZ}$ and $\boldsymbol{XY}$ into Eq. (58), we obtain Eq. (11) in the main text.

Appendix C Nearest-Neighbour Bound State Solution

Here we derive the analytic expression for the bound state energy and decay rate given by Eq. (24) when the XY interaction and nonlocal dissipation is nearest-neighbour. We can evaluate the integrals as defined in Eq. (13) using contour integration. Substituting $z=\exp(iq)$, the integral transforms into

[TABLE]

where we have defined $\alpha=(2V_{12}+i\Gamma_{12})\cos(Qa/2)$ and $t=8J_{z}-2\Delta+i\Gamma$. The integral has a pole of order $m$ at $z=0$ and simple poles at $z_{\pm}=(\Omega+t)/2\alpha\pm\sqrt{((\Omega+t)/2\alpha)^{2}-1}$. The two poles only coincide at $|z|=1$, so the case of double poles can be ignored for the derivation. Evaluating the integrals gives

[TABLE]

where the $\pm 1$ sign depends on whether $z_{+}$ or $z_{-}$ lie in the contour. Substituting these solutions into the bound state equation, Eq. (12), we obtain the equation

[TABLE]

which gives the solution $\Omega+t=4+\alpha^{2}/(4J_{z})$.

Appendix D Next-Nearest-Neighbour Bound State Solution

To derive the analytic expression for the next-nearest-neighbour bound state solution given by Eq. (27), we use the substitution $z=e^{iq}$ to transform the integral in Eq. (13) into the following contour integral

[TABLE]

where $\beta=(2V_{13}-i\Gamma_{13})\cos(Qa)$, $t=8J_{z}-2\Delta+i\Gamma$ and $\alpha=(2V_{12}-i\Gamma_{12})\cos(Qa/2)$. The quartic in the denominator is palindromic, which means the solutions obey a quadratic in $(z+1/z)$. Therefore, if $z$ is a solution to the quartic, then so too is $1/z$, and this immediately indicates that only two of the four roots can exist inside the contour. We also find that the residue of the roots $1/z$ and $z$ only differ by a sign. The integrals in Eq. (62) can therefore be evaluated to give

[TABLE]

where

[TABLE]

The sign of $F_{1/2}$ depends on whether the root $z_{1/2}$ or its inverse lies inside the contour. Substituting the integral solutions into the bound state equation, Eq. (12), gives

[TABLE]

We can now solve Eq. (65) to obtain the solution given in Eq. (27) in the main text. There is also the possibility of a double root when $\Omega+t=2\beta+\alpha^{2}/(4\beta^{2})$. When this is the case, the denominator the integrals in Eq. (12) can be simplified to $(4\beta\cos(q)-\alpha)^{2}/(4\beta)$. We can then evaluate the NNN integrals without using contour integration, but find these solutions do not obey the bound state solution.