Topological properties of a dense atomic lattice gas

TL;DR

This paper explores topological phases in a dense 2D atomic lattice gas, revealing how long-range interactions and dissipation influence edge transport and can break traditional bulk-boundary relations.

Contribution

It demonstrates the existence of topological phases in a dense atomic lattice with long-range interactions, highlighting robustness of edge transport despite dissipation and defects.

Findings

Long-range interactions cause one-sided divergences in energy spectra.

Topological edge transport remains robust despite dissipation.

Standard bulk-boundary relation can be broken by long-range interactions.

Abstract

We investigate the existence of topological phases in a dense two-dimensional atomic lattice gas. The coupling of the atoms to the radiation field gives rise to dissipation and a non-trivial coherent long-range exchange interaction whose form goes beyond a simple power-law. The far-field terms of the potential -- which are particularly relevant for atomic separations comparable to the atomic transition wavelength -- can give rise to energy spectra with one-sided divergences in the Brillouin zone. The long-ranged character of the interactions has another important consequence: it can break of the standard bulk-boundary relation in topological insulators. We show that topological properties such as the transport of an excitation along the edge of the lattice are robust with respect to the presence of lattice defects and dissipation. The latter is of particular relevance as dissipation and coherent interactions are inevitably connected in our setting.

Equations

$$H$$

$$\mathcal{D}(\rho )$$

$${\overline{V}}_{jl}$$

$${\overline{\Gamma }}_{jl}$$

$$A_{jl}$$

$$B_{jl}$$

$$A_{jl}^{\prime }$$

$$B_{jl}^{\prime }$$

$$\overline{V}_{jl}\!=\!\left(\!\!\begin{array}[]{cc}\left(A_{jl}+\Delta\delta_{jl}\right)\!\!\left(1-\varepsilon^{2}\right)&\!\!B_{jl}e^{-2i\phi_{jl}}\!\!\left(1-\frac{\varepsilon^{2}}{2}\right)\\ B_{jl}e^{2i\phi_{jl}}\!\!\left(1-\frac{\varepsilon^{2}}{2}\right)&A_{jl}-\Delta\delta_{jl}\end{array}\!\!\right)\!,$$

$$H=\hslash \sum _{\mathbf{k}}{\psi }_{\mathbf{k}}^{†}\cdot {\overline{V}}_{\mathbf{k}}\cdot {\psi }_{\mathbf{k}},$$

$${\overline{V}}_{\mathbf{k}}=\sum _{\mathbf{r}\ne 0}{e}^{i\mathbf{k}\cdot \mathbf{r}}{\overline{V}}_{\mathbf{r}}$$

$$C=\frac{1}{2\pi i}{\int }_{T^2}d\mathbf{k}{F}_{xy}(\mathbf{k}),$$

$$q=\frac{1}{N_{\text{edge}}}\text{Im}\left[\text{ln}\sum _{l=1}^{N_{\text{edge}}}\frac{{\mathbf{d}}_l^{†}{\mathbf{d}}_{l+1}}{|{\mathbf{d}}_l^{†}{\mathbf{d}}_{l+1}|}\right],$$

$$\overline{W}_{jl}=\left(\begin{array}[]{ccc}A_{jl}&B_{jl}e^{-2i\phi_{jl}}&0\\ B_{jl}e^{2i\phi_{jl}}&A_{jl}&0\\ 0&0&C_{jl}\end{array}\right)$$

$$C_{jl}$$

$$\overline{\gamma}_{jl}=\left(\begin{array}[]{ccc}A^{\prime}_{jl}&B^{\prime}_{jl}e^{-2i\phi_{jl}}&0\\ B^{\prime}_{jl}e^{2i\phi_{jl}}&A^{\prime}_{jl}&0\\ 0&0&C^{\prime}_{jl}\end{array}\right),$$

$$C_{jl}^{\prime }$$

$$\mathbf{r}={\mathbf{\tau }}_{\alpha \beta }+\mathbf{T},$$

$$V(\mathbf{k})=\sum _{\mathbf{r}\ne 0}{\mathrm{e}}^{i\mathbf{kr}}v(\mathbf{r})$$

$$1={}^S\xi (r)+{}^L\xi (r),$$

$$V(\mathbf{k})={}^{S}V(\mathbf{k})+{}^{L}V(\mathbf{k})-{}^{\mathrm{self}}V(\mathbf{k}),$$

$${}^{S}V(\mathbf{k})$$

$${}^{L}V(\mathbf{k})$$

$${}^{\mathrm{self}}V(\mathbf{k})$$

$${}^{L}V(\mathbf{k})$$

$$\sum _{\mathbf{T}}{\delta }^2(\mathfrak{r}-\mathbf{T})=\sum _{\mathbf{G}}{\mathrm{e}}^{i\mathbf{G}\mathbf{r}}$$

$${}_{1,n}v(\mathbf{r})$$

$${}_{2,n}v(\mathbf{r})$$

$${}_{3,n}v(\mathbf{r})$$

$${}_{1,n}\mathfrak{V}$$

$${}_{1,n}\mathfrak{V}={}_{1,n}{\mathfrak{V}}_0^c+{}_{1,n}{\mathfrak{V}}_c^{\infty },$$

$${}_{1,n}{\mathfrak{V}}_a^b={\kappa }^{-1}{\int }_a^b\mathrm{d}r\left(1-{\mathrm{e}}^{-\gamma r}\right)\cos (\kappa r)J_n(\tilde{G}r).$$

$$J_n(\tilde{G}r)\approx \sqrt{\frac{2}{\pi \tilde{G}}}\frac{1}{\sqrt{r}}\cos \left(\tilde{G}r-{\varphi }_n\right),$$

$${}_{1,n}{\mathfrak{V}}_c^{\infty }\approx I_0^{\infty }-I_0^c,$$

$$I_a^b=\frac{1}{\kappa }\sqrt{\frac{2}{\pi \tilde{G}}}{\int }_c^d\mathrm{d}r\frac{1-{\mathrm{e}}^{-\gamma r}}{\sqrt{r}}\cos (\kappa r)\cos (\tilde{G}r-{\varphi }_n).$$

Full text

Introduction. Recently, the pursuit of topological phases in quantum many-body systems has been the focus of intense research. The potential application of these topological states for robust quantum computation Kitaev (2003); Nayak et al. (2008) is one of the driving forces for this increased interest. So-called topological insulators are usually characterized by bulk bands separated by a gap and the presence of gapless edge states whose properties are topologically protected against local perturbations such as external disorder or noise Hasan and Kane (2010); Yao et al. (2013a). Paradigmatic examples of these include the integer and fractional quantum Hall effects, which were initially realized on two-dimensional electron gases subject to strong magnetic fields Klitzing et al. (1980); Laughlin (1981); Tsui et al. (1982); Laughlin (1983). Since their discovery, several lattice models that do not require external magnetic fields have been proposed and some realized experimentally Haldane (1988); Ruostekoski et al. (2002); Kane and Mele (2005); Novoselov et al. (2007); Tsukazaki et al. (2007); König et al. (2007); Raghu et al. (2008); Neupert et al. (2011); Wang et al. (2011); Tang et al. (2011); Sun et al. (2011); Wang et al. (2012); Dauphin et al. (2012); Cooper and Dalibard (2013); Barkeshli et al. (2015); Andrijauskas et al. (2015).

Due to the high degree of experimental control that is achievable nowadays, cold atoms and molecules have been proposed as platforms for the exploration of novel topological phases of quantum matter Goldman et al. (2016). In particular, many-body systems that display long-range interactions Zhang et al. (2014); Vodola et al. (2014); Gong et al. (2016); Behrmann et al. (2016) – such as polar molecules Peter et al. (2012); Yao et al. (2012); Manmana et al. (2013); Yao et al. (2013b); Peter et al. (2015), atoms with large magnetic dipoles Yao et al. (2015), and Rydberg atoms Maghrebi et al. (2015) – have been shown to feature topologically non-trivial flat bands and fractional quantum Hall states.

In this paper, we explore the topological properties of a two-dimensional lattice of atoms where long-range interactions arise intrinsically via coherent light scattering between internal atomic states Bettles et al. (2015). We consider the full interaction – going thus beyond the usually employed near-field approximation – which is of relevance for atomic separations comparable to the atomic transition wavelength. This scenario is currently studied in various contexts, e.g., the exploration of collective light scattering and super- and sub-radiant decay Svidzinsky et al. (2010); Keaveney et al. (2012); Bienaimé et al. (2013); Pellegrino et al. (2014); Kwong et al. (2014); Bettles et al. (2015); Bromley et al. (2016); Jenkins et al. (2016). We study two different lattice geometries (a square and a honeycomb lattice) and find that this simple system can support topologically non-trivial phases. The often neglected far-field terms of the interactions lead to one-sided divergences in the single-particle spectrum. We explore the consequences of the long-range character of the interaction on the relation between the topological properties of the bulk and the boundary of a finite size system. Remarkably, we find that the standard bulk-boundary relation Hatsugai (1993a, b), which is well established for short-range interactions, does not generally hold. Furthermore, we find that neither the presence of lattice defects nor dissipation destroy the topological properties of the system. The latter is of particular importance, as the considered system is inevitably open: both dissipation and coherent exchange interactions originate from the coupling of the atoms to the radiation field.

The system. We consider a two-dimensional optical lattice with $N$ sites lying in the $xy$ plane (lattice spacing $a$), where each site is occupied by a single atom. The $j$-th atom is located at position ${\mathbf{r}}_{j}$ with internal levels $\ket{g}_{j}$ (the ground state) and $\ket{-}_{j},\ket{0}_{j},\ket{+}_{j}$ (the excited states) corresponding to the $J=0$ and $J=1$ total angular momentum manifolds, respectively [see Fig. 1(a)]. This level structure is naturally available in a variety of systems such as alkaline-earth-metal atoms Gorshkov et al. (2009, 2010); Olmos et al. (2013), dysprosium atoms Yao et al. (2015), polar molecules Peter et al. (2015); Gorshkov et al. (2013); Syzranov et al. (2014) or Rydberg systems Maghrebi et al. (2015).

The coupling of the atoms to the quantized multimode radiation field results in an effective long-range exchange interaction and collective dissipation Lehmberg (1970); Agarwal (1970); James (1993). Within the dipole and Born-Markov approximations, the dynamics of the atomic system is described by the master equation $\dot{\rho}=-\frac{i}{\hbar}\left[H,\rho\right]+{\cal D}(\rho)$, with

[TABLE]

Here, $\mathbf{d}_{j}=\left(\left|g\right>_{j}\!\left<+\right|,\left|g\right>_{j}\!\left<-\right|\right)^{T}$ and

[TABLE]

Note, that we have used here the fact that the dynamics of the $\ket{0}$ level and the $\{\ket{-},\ket{+}\}$ manifold are decoupled, see Appendix A for details. The coefficients in Eqs. (5) and (8) read

[TABLE]

where $\Gamma$ is the single atom decay rate and $\kappa_{jl}\equiv 2\pi r_{jl}/\lambda$, with $\lambda$ being the wavelength of the transition from the ground to the excited state manifold [see Fig. 1(a)], $r_{jl}=|\mathbf{r}_{j}-\mathbf{r}_{l}|$, and $\phi_{jl}={\rm arg}(\mathbf{r}_{j}-\mathbf{r}_{l})$ is the polar angle between the $j$-th and $l$-th atom.

The many-body Hamiltonian (1) conserves the number of excitations in the system and describes their exchange among the atoms. The coefficients in this Hamiltonian [Eqs. (9) and (10)] decay as $1/r^{3}$ for short distances ($r\leq\lambda$) and as $1/r$ for $r\gg\lambda$ [see Fig. 1(b)]. We will study the topological properties of the system in an intermediate regime, where the full potential needs to be considered. In particular, in the remainder of the paper we will fix the ratio $a/\lambda=0.1$ Olmos et al. (2013); Pellegrino et al. (2014) (note that the effects we show are not constrained to this specific value) and systematically investigate two specific lattice geometries, the square lattice (SL) and the honeycomb lattice (HL) [see Figs. 1(c) and (d), respectively]. We will proceed by analyzing first the band structure of the Hamiltonian (1) and then studying the effect of the dissipation (2) on the edge transport in the topologically non-trivial phases.

In its current formulation, the Hamiltonian (1) is symmetric under time reversal and hence is topologically trivial. In order to possibly find topologically non-trivial phases we consider two ways of breaking the time reversal symmetry. First, we lift the degeneracy of the states $\left|\pm\right>$, which can be achieved by means of a magnetic field such that the two states’ energies are shifted by $\pm\Delta$ [see Fig. 1(a)]. Secondly, we couple off-resonantly the state $\left|+\right>$ via a microwave field to an auxiliary hyperfine state, $\left|2\right>$ Peter et al. (2015). Up to second order in $\varepsilon\approx\Omega/2\delta\ll 1$, where $\Omega$ and $\delta$ are the Rabi frequency and the detuning of the microwave field, respectively, the $\left|+\right>$ state gets dressed and the potential (5) is modified as

[TABLE]

where $\delta_{jl}$ is the Kronecker delta symbol.

Band structure and divergences. To obtain the band structure we express (1) in the reciprocal space as

[TABLE]

where $\psi_{\mathbf{k}}=\frac{1}{\sqrt{N}}\sum_{l}\mathbf{d}_{l}e^{i\mathbf{k}\cdot\mathbf{r}_{l}}$ and

[TABLE]

with $\mathbf{r}$ being the relative distance between any pair of atoms. Note that the potential (13) (due to the terms that decay as $1/r$) makes the sum in Eq. (15)) not convergent. Potentials with similar features have been shown to lead to interesting physical effects, such as supersonic spreading of the correlations Raghu et al. (2008); Neupert et al. (2011); Wang et al. (2011); Tang et al. (2011); Sun et al. (2011); Wang et al. (2012); Dauphin et al. (2012); Cooper and Dalibard (2013); Storch et al. (2015). Such potentials are also encountered in self-gravitating systems Campa et al. (2014) or in electrons in solids subject to the Coulomb force Kantorovich (2004).

To evaluate the sum (15) we employ the Ewald summation technique Bonsall and Maradudin (1977); Kantorovich (2004), see Appendix B for details. As a result, we find that the potential (15) features one sided divergences in the spectrum occurring for wavevectors lying on a circle of radius $\kappa=2\pi/\lambda$ in reciprocal space. More specifically, the potential diverges as $\lim_{k\to\kappa^{+}}\braket{V_{\mathbf{k}}}=-\infty$, while $\lim_{k\to\kappa^{-}}\braket{V_{\mathbf{k}}}$ remains finite (see Appendix B). This situation is depicted in Fig. 1(e) (SL) and Fig. 1(f) (HL). Note that in this work we only study cases where the divergence falls inside the Brillouin zone, $a/\lambda<1/2$ for SL and $a/\lambda<1/3$ for HL.

Chern numbers and edge states. In the language of differential geometry, the topological properties are studied in terms of differentiable fibre bundles assuming a differentiable Hamiltonian map $H:T^{2}\rightarrow\mathcal{M}$, mapping the Brillouin zone (represented as a two-dimensional torus $T^{2}$) to some target space $\mathcal{M}$ Nakahara (2003). Here the topology is characterized by the Chern number defined as

[TABLE]

where $F_{xy}(\mathbf{k})=\partial_{x}A_{y}(\mathbf{k})-\partial_{y}A_{x}(\mathbf{k})$, $A_{\mu}(\mathbf{k})=\langle n(\mathbf{k})|\partial_{\mu}|n(\mathbf{k})\rangle$ is the Berry connection, $\ket{n(\mathbf{k})}$ is an eigenstate of the Hamiltonian (14), $\partial_{\mu}=\partial/\partial k_{\mu},\;\mu=\{x,y\}$, and ${\mathbf{k}}\in T^{2}$.

In the present case the assumption of differentiability is not satisfied due to the discontinuity of the potential (15). However, formally it is still possible to evaluate the Chern number using the algorithm of Ref. Fukui et al. (2005), where special care is taken to avoid the points in the Brillouin zone where $\overline{V}_{\mathbf{k}}$ diverges. This corresponds to evaluating (16) with an effectively bounded Hamiltonian, which in turn yields integer values of $C$. Note also that the divergence appears as a consequence of considering the potential in the thermodynamic limit, a situation that is ill-defined as the finite propagation time of the radiation modes mediating the exchange of photons in the system is neglected Lehmberg (1970); Milonni and Knight (1974).

The resulting phase diagram is shown in Fig. 2(a) (SL) and Fig. 2(b) (HL). We note that in the SL case it is required that $\varepsilon\neq 0$ in order to access the topologically non-trivial region and that in HL several topologically distinct regions can be accessed by tuning $\Delta$ alone (i.e., $\varepsilon=0$).

A hallmark of the non-trivial topology of the bulk is the appearance of edge states in finite-size systems. Considering an infinite strip keeping the system finite in the $y$ direction we plot the corresponding band structure in Fig. 2(c) (SL) and Fig. 2(d) (HL), where the edge states appear in the bulk gaps. Here, the size of the filled circles corresponds to the weight of the given eigenmode on one (red) or the other (blue) edge of the infinite strip.

We note that in the SL case the edge states are nearly degenerate and ”energetically hidden” Weber (2015); Peter et al. (2015) making it difficult to access them experimentally. In the HL case, relatively well separated bands form. Moreover, it is easy to verify that in the example shown indeed the Hall conductivity or, equivalently, the sum of the Chern numbers of the filled bands below a band gap, determines the net number of edge states crossing that gap Hatsugai (1993a); Liu et al. (2012).

Finite system: quasi-momentum and driven-dissipative dynamics. One important consequence of the long-ranged character of the interactions in the context of topological systems is that the bulk-edge correspondence, well established for systems with short-range potentials does not necessarily hold, as we exemplify in the following.

In Fourier space, the edge modes can be seen as modes spanning between different energy bands. In a finite system, to obtain a similar measure of whether any given mode has a chiral edge conductivity associated with it, we define a quasi-momentum following a similar approach to Hafezi et al. (2013); Weber (2015): When an excitation hops from site to site, it accumulates a phase. If this phase of hopping is constant around the edge of the lattice, the hopping occurs around the edge in a given direction. If the hopping phase is random between different lattice sites, however, there is no preferred direction of hopping. We therefore define the quasi-momentum $q$ as the average phase difference between neighbouring lattice sites around the edge of the lattice:

[TABLE]

where $N_{\text{edge}}$ is the number of atoms in the outermost edge of the lattice. In Fig. 3(a) we show the quasi-momentum spectrum for a HL of an hexagonal shape with a total of 486 atoms 111The effects we show in this work can also be observed in smaller lattices.. We have overlayed the band structure results for an infinite strip with the same number of atoms in the $y$ direction for comparison.

Another defining characteristic of an edge state is that the vast majority of its weight is on the physical edge of the system. Writing an eigenstate of the system as $\ket{\varphi}=\sum_{\mu,j}c^{j}_{\mu}\ket{\mu}_{j}$, where $\mu\in\{-,+\}$ and $j$ index the internal degrees of freedom and the spatial position, respectively, we quantify the weight of the eigenstate on the edge of the system as $w=\sum_{\mu,j\,\in\,{\rm edge}}|c^{j}_{\mu}|^{2}$. The blue and red dots in Fig. 3(a) correspond to states where $w<0.95$ and $w>0.95$, respectively. The spatial configuration of two possible edge states from the middle and lowest bulk gap is shown in Figs. 3(i) and (ii) [which we denote as states (i) and (ii) in the following]. A bulk state is shown for comparison in Fig. 3(iii). Here, the size of the red circles correspond to the weight of the state at a given site.

In order to further characterize the edge states in the finite system, we model the dynamics of an excitation in the lattice under the Hamiltonian (1) and weak driving. Specifically, we drive a single atom on the leftmost corner of the HL [marked by a star in the leftmost panels of Figs. 3(b) and (c)] with weak driving strength and we choose the detuning to be $-2\Gamma$ and $-10\Gamma$ [Figs. 3(b) and (c), respectively], so that we resonantly address the states (i) and (ii), respectively. We observe that initially in both (b) and (c) the excitation dynamics is chiral with clockwise propagation before the steady state is reached.

Next, as one of the most appealing features of the edge states is their robustness against local perturbations, we simulate the excitation dynamics including $5\%$ of uniformly distributed lattice defects, i.e., empty sites due to an imperfect loading of atoms [fourth panels in Figs. 3(b,c)]. We observe a clear difference in the excitation propagation when driving the state (i), where the propagation is clearly interrupted, and the state (ii), where the excitations still populate a significant part of the edge. This strongly indicates that only the second one, state (ii), is an actual edge state. This is in contrast to what one might expect from the equilibrium bulk analysis, where the edge state conductivity, proportional to the sum of the Chern numbers of the filled lower lying bands, is the same for both states (i) and (ii) [see the red line in Fig. 2(b)].

Finally, as the mechanism that is responsible for the existence of coherent interactions is also responsible for the dissipation in the system, we simulate the dynamics above including the dissipative term (2). Remarkably, we find that while in the case (i) the dissipation makes the excitation decay after only a few sites, when exciting the edge state (ii) the excitation becomes sub-radiant, i.e., it decays with a rate much smaller than the single atom decay rate $\Gamma$.

Conclusions and outlook. We have found the existence of topologically non-trivial phases in a dense atomic two-dimensional lattice system coupled to the radiation field. We show that one can excite edge states that allow for the transport of an excitation over long distances along the edge of the lattice that are robust to the presence of defects. Moreover, these edge states are remarkably long-lived due to the collective character of the dissipation (a detailed analysis of this effect will be the aim of further investigations). Finally, we have found that, due to the long-range character of the interactions, the bulk-boundary relations, well established for topological insulators with short-range interactions, are not generally valid in our setting.

Appendix A Coupling to the multimode radiation field

The atomic system is coupled to the quantized multimode radiation field. Within the dipole approximation, the Hamiltonian that describes the dynamics of the atomic ensemble and the radiation field reads $H_{\mathrm{af}}=\sum_{j=1}^{N}\hbar\omega_{\mathrm{a}}\mathbf{b}_{j}^{\dagger}\cdot\mathbf{b}_{j}+\sum_{\mathbf{q}\lambda}\hbar\omega_{\mathbf{q}}a_{\mathbf{q}\lambda}^{\dagger}a_{\mathbf{q}\lambda}+i\hbar\sum_{j=1}^{N}\sum_{\mathbf{q}\lambda}\mathbf{g}_{\mathbf{q}\lambda}\cdot\left(a_{\mathbf{q}\lambda}^{\dagger}\mathbf{s}_{j}e^{-i\mathbf{q}\cdot\mathbf{r}_{j}}-\mathbf{s}_{j}a_{\mathbf{q}\lambda}e^{i\mathbf{q}\cdot\mathbf{r}_{j}}\right)$. The first term represents the atomic energy, with $\hbar\omega_{\mathrm{a}}=2\pi\hbar c/\lambda$ being the energy difference between the ground ($\left|g\right>$) and the three degenerate $J=1$ manifold states ($\left|-\right>$, $\left|0\right>$ and $\left|+\right>$ for $m_{J}=-1,0$ and $+1$, respectively) and the vector transition operator for the $j$-th atom (located at $\mathbf{r}_{j}$) being defined as $\mathbf{b}_{j}=\left(\left|g\right>_{j}\!\left<+\right|,\left|g\right>_{j}\!\left<-\right|,\left|g\right>_{j}\!\left<0\right|\right)$. The second term of the Hamiltonian represents the radiation field, and here $\hbar\omega_{\mathbf{q}}$ is the energy of a photon with momentum $\mathbf{q}$ and polarization $\lambda$ and $a_{\mathbf{q}\lambda}$ is the annihilation operator of such a photon ($[a_{\mathbf{q}\lambda},a_{\mathbf{q}^{\prime}\lambda^{\prime}}^{\dagger}]=\delta_{\mathbf{q}\mathbf{q}^{\prime}}\delta_{\lambda\lambda^{\prime}}$). The last term represents the coupling between the two systems, with $\mathbf{s}_{j}=\mathbf{b}_{j}^{\dagger}+\mathbf{b}_{j}$ and the coupling coefficient being given by $\mathbf{g}_{\mathbf{q}\lambda}=d\sqrt{\frac{\omega_{\mathbf{q}}}{2\epsilon_{0}\hbar V}}\hat{e}_{\mathbf{q}\lambda}$, where $V$ is the quantization volume, $d$ the modulus of the transition dipole moment between the ground and the three degenerate excited states, and $\hat{e}_{\mathbf{q}\lambda}$ the unit polarization vector of the photon ($\mathbf{q}\cdot\hat{e}_{\mathbf{q}\lambda}=0$).

Within the dipole and Born-Markov approximations, and tracing out the environment, one obtains the quantum master equation for the atomic density matrix $\rho$, $\dot{\rho}=-\frac{i}{\hbar}\left[H,\rho\right]+{\cal D}(\rho)$, which describes the dynamics of the atomic degrees of freedom Lehmberg (1970); Agarwal (1970); James (1993). The many-body Hamiltonian that describes the coherent evolution of the system reads $H=\hbar\sum_{j\neq l}\mathbf{b}^{\dagger}_{j}\cdot\overline{W}_{jl}\cdot\mathbf{b}_{l}$, characterized by the coefficient matrix

[TABLE]

with $A_{jl}$ and $B_{jl}$ being defined in the main text [Eqs. (5) and (6), respectively] and

[TABLE]

for $j\neq l$. The dissipator ${\cal D}(\rho)$ reads ${\cal D}(\rho)=\sum_{jl}\mathbf{b}_{j}\cdot\overline{\gamma}_{jl}\cdot\rho\mathbf{b}_{l}^{\dagger}-\frac{1}{2}\left\{\mathbf{b}_{j}^{\dagger}\cdot\overline{\gamma}_{jl}\cdot\mathbf{b}_{l},\rho\right\}$, where

[TABLE]

with $A^{\prime}_{jl}$ and $B^{\prime}_{jl}$ being defined in the main text [Eqs. (7) and (8), respectively] and

[TABLE]

It is clear from the form of Eqs. (18) and (20) that the dynamics of the internal level $m_{J}=0$ is disconnected from the other two. We will consider thus only the dynamics on the subspace formed by the states $m_{J}=\pm 1$.

Appendix B Ewald summation

B.1 General notions

We briefly review the Ewald summation technique Bonsall and Maradudin (1977); Kantorovich (2004) introduced originally to cope with the evaluation of the electrostatic potential energy of a gas of electrons interacting through Coulomb force. Here we wish to evaluate the discrete sum (15). Considering two-dimensional lattices with unit cell defined by two lattice vectors $(\mathbf{a}_{1},\mathbf{a}_{2})$, a vector ${\mathbf{r}}$ connecting two atoms $\alpha$ and $\beta$ can be written as

[TABLE]

where ${\bf T}=(m\mathbf{a}_{1},n\mathbf{a}_{2}),\;m,n\in\mathbb{Z}$ is the lattice translation vector, ${\boldsymbol{\tau}}_{\alpha\beta}={\boldsymbol{\tau}}_{\beta}-{\boldsymbol{\tau}}_{\alpha}$ and ${\boldsymbol{\tau}}_{\alpha}$ denotes the position of the $\alpha$-th atom in the unit cell. We also define the reciprocal lattice vectors in the usual way as ${\mathbf{G}}=(m\mathbf{b}_{1},n\mathbf{b}_{2}),\;m,n\in\mathbb{Z}$, where $\mathbf{b}_{j},\;j=1,2$ are the primitive reciprocal lattice vectors.

Let us consider an algebraically decaying function of the atomic separation $v({\mathbf{r}})$ and say we want to evaluate the sum

[TABLE]

The Ewald’s trick consists of separating the contributions of $v$ to the sum to the short (S) and long (L) ranged part by introducing a regulating function $\xi$ through the identity

[TABLE]

so that

[TABLE]

where

[TABLE]

The function $\xi$ is chosen such that ${}^{S}\xi(r)v({\bf r})$ decays exponentially fast as $r\rightarrow\infty$ and (26c) is finite. Due to the fast decay of ${}^{S}\xi$, the short-range part can be readily obtained by evaluating the sum in (26a) in real space. We are thus left with the evaluation of (26b).

We proceed as follows. Using (22), (26b) becomes

[TABLE]

where we in the third line we have used the identity (normalizing the volume of the unit cell to 1)

[TABLE]

and in the last line we made a substitution of variables ${\bf r}={\boldsymbol{\tau}}_{\alpha\beta}+{\bm{\mathfrak{r}}}$ and defined $\widetilde{\bf G}={\bf k+G}$.

Next, due to the form of (5),(9),(10), we consider the following functions $v$

[TABLE]

where $n=0,2$ and $\kappa=2\pi/\lambda$.

B.2 Analytical approach to the $1/r$ term: the divergence

We start with the evaluation of the integral in (27) with the function ${}_{1,n}v({\bf r})$

[TABLE]

where we have chosen ${}^{L}\xi(r)=1-{\rm e}^{-\gamma r}$, $\gamma>0$ is an arbitrary control parameter, $\eta$ is the polar angle of $\widetilde{\bf G}$ and $J_{n}$ is the $n$-th Bessel function of the first kind. Next, we split the integration as

[TABLE]

where

[TABLE]

Finally, we substitute the asymptotic expression for the Bessel function (for large $\widetilde{G}r$)

[TABLE]

where $\phi_{n}=\frac{\pi}{4}(2n+1)$, to ${}_{1,n}\mathfrak{V}_{c}^{\infty}$ so that

[TABLE]

where

[TABLE]

Importantly, we note, that only the term $I_{0}^{\infty}$ can give rise to the divergence in (25), which we now analyze in detail. Substituting the identity

[TABLE]

and using

[TABLE]

we get

[TABLE]

This amounts to evaluation of the integral

[TABLE]

where $\mathcal{G}=\widetilde{G}\pm\kappa$. The final integration in (39) can be done by methods of contour integration along the contour depicted in Fig. 4, i.e.

[TABLE]

Using the property of the branch cut, the integration along $C_{2}$ becomes

[TABLE]

Omitting the c.c. term, we express the leftmost part of (39) as

[TABLE]

It can be immediately seen, that $\int_{R}\rightarrow 0$ since $g(z)\rightarrow R^{-\frac{3}{2}}$ in the limit of large $|z|$. Next, since $\gamma>0$, assuming $\mathcal{G}\neq 0$ and substituting $z=\epsilon{\rm e}^{i\varphi}$, we get from the estimation lemma

[TABLE]

We are thus left with

[TABLE]

where we have defined the function $f$ and the poles

[TABLE]

Plugging (B.2) back to (39) and evaluating the complex conjugated term we get

[TABLE]

Next, we note that this result can be further simplified in specific cases: for $\mathcal{G}\in\mathbb{R}$ and $n=0,2$ respectively, we get the following options

[TABLE]

Lets first consider the $n=0$ case. Using (48), the relation (47) simplifies to

[TABLE]

and similar expression is obtained for $n=2$. Those results can be combined to a single expression

[TABLE]

where $\theta$ is the Heaviside step function. Remarkably, the function (50) is finite in the limit $\mathcal{G}\rightarrow 0^{-}$, but diverges as $\mathcal{G}\rightarrow 0^{+}$, which constitutes the main analytic result of this section.

Remark: The discontinuity and one-sided divergence could be seen already in (39) by setting $\mathcal{G}=0$. This would result in the diverging integral in the vicinity of the origin, namely the part $\int_{\epsilon}\propto\int{\rm d}t\frac{1}{t\sqrt{t}}$ in (40), which, according to the estimation lemma, scales as $\propto\frac{1}{\sqrt{\epsilon}}$.

B.3 Numerical implementation

In the previous section we have found that (26a),(26c) can be evaluated numerically, while (26b) can be written as

[TABLE]

Here, the expressions ${}_{1,n}\mathfrak{V}_{0}^{c}$ and $I_{0}^{c}$ are finite and can by obtained numerically. In principle similar decomposition can be obtained when considering the functions (29b),(29c), however we found it more convenient to evaluate the corresponding element (23) by direct summation in real space.

The final subtlety is related to the fact, that for ${}^{L}\xi=1-{\rm e}^{-\gamma r}$ the self-interaction term (26c) is not defined for ${}_{1,n=2}v$ as it depends on the direction in which the origin is approached (which is due to the dependence of ${}_{1,n=2}v$ on the polar angle of $\bf r$). This problem can be circumvented by considering ${}^{L}\xi=1-{\rm e}^{-\gamma r^{2}}$ instead in which case the self-interaction term vanishes. Here, it is still possible to evaluate the integral in (35) symbolically, which we then use in (51).

In (51) the tuning parameters $\gamma,c$ are chosen such that the approximations used (eq. (33)) are well satisfied.