Multiband and array effects in matter-wave-based waveguide QED

Experimenting with the light-matter interface at the quantum level Hood et al. (2016); Roy et al. (2017); Carusotto et al. (2020) has led to the discovery of new exciting phenomena, including implementations of photonic quantum matter Chang et al. (2018); Carusotto and Ciuti (2013), effective long-range interactions between neutral atoms Douglas et al. (2015), chiral quantum optics Mitsch et al. (2014), nondestructive photon counting Malz and Cirac (2020) and quantum logic gates Zheng et al. (2013); Paulisch et al. (2016) among others. The plethora of these developments contrasts with the conceptual simplicity of coupling one or many quantum emitters to a photonic waveguide (a.k.a. waveguide QED); photons are the ideal carriers of quantum information as they can travel large distances without colliding with one another, while atoms can access and store this information by absorbing the photons.

Implementations using ultracold atoms in optical lattices, originally proposed by de Vega et al. (2008), challenge this paradigm as they reproduce the same physics while switching the roles of matter and light which now become the radiation and emitters, respectively. This comes with advantages, such as accessibility to different parameter regimes and a high tunability of the system Navarrete-Benlloch et al. (2011); Stewart et al. (2017a) with the potential to enable state-of-the-art applications including the simulation of giant atoms González-Tudela et al. (2019), perfect subradiance González-Tudela and Cirac (2018) and topological effects Bello et al. (2019). The experimental realization of these systems has been achieved recently Krinner et al. (2018); Stewart et al. (2020) resulting in the observation of phenomena such as tunable Lamb shifts, bound-state beats, non-Markovian decay dynamics and time-of-flight pictures of the bound states.

In the analysis of waveguide-QED systems Rzazewski et al. (1982); Lewenstein et al. (1988); Kofman et al. (1994); Lambropoulos et al. (2000); Berman and Ford (2010), a key idea is that the emission properties can be drastically altered by manipulating the mode distribution of the radiation Bykov (1975); Yablonovitch (1987); John (1987). However, due to the difficulty of a multiband analysis Bykov (1975), the emission of matter wave radiation has only been analyzed for the cases of a single emitter coupled to modes with one parabolic band edge de Vega et al. (2008); Stewart et al. (2017a); Krinner et al. (2018) or one sinusoidal band Calajó et al. (2016); Stewart et al. (2020). In this paper, we go beyond these restrictions and present how to evaluate the exact lattice functions required to solve the dynamics. Using techniques of complex analysis, we find that this can be done in an efficient way by using infinite products. These allow us to solve the system with broader generality (i) by considering the action of multiple emitters in the system; (ii) by accounting for the full band dispersion relation for the radiated modes; (iii) going beyond the Markovian regime, as we solve for arbitrarily high couplings; and (iv) also accounting for the finite size of the emitters. For simplicity, and in agreement with the conditions of the experiments Krinner et al. (2018); Stewart et al. (2020), we restrict ourselves to 1D non-interacting systems within the single excitation subspace.

This paper is organized as follows. In Sec. II, we give an introduction to the ultracold atom platform and its connection with the Hamiltonian for spontaneous emission of matter waves. In Sec. III, we derive the formal solution of the equations of motion for multiple emitters. These equations can be evaluated using the infinite-product representations presented in Sec. IV. We focus on the case of a single emitter in Sec. V, where we describe the system spectrum in detail and compare the analysis results with the experiments. In Sec. VI we study the form of the bound states, consisting of matter waves dynamically anchored to the emitter. Finally, in Sec. VII we review the case of an infinite array of emitters and the emergence of polaritons that result from the hybridization between the emitter array and the potential in which the matter waves propagate.

We consider a two-level atom (states “red” $\left|a\right\rangle$ and “blue” $\left|b\right\rangle$) of mass $m$ in a 1D state-dependent optical lattice of recoil momentum $k$ and potential $\mathbf{V}_{a(b)}(z)=V_{a(b)}\sin^{2}(kz)$, as experimentally realized in Krinner et al. (2018); Stewart et al. (2020). When the atom is in $\left|a\right\rangle$, the lattice potential is deep ($V_{a}\gg E_{r}=(\hbar k)^{2}/2m$) and the atom is locally confined to the harmonic-oscillator ground state $\phi_{j}(z)$ in site $j$, with negligible tunneling between sites. An oscillating external field of frequency $\omega_{\mu}$ couples the atom to $\left|b\right\rangle$, which experiences a much shallower potential $\lvert V_{b}\rvert\ll V_{a}$, with tunneling $J_{b}/\hbar$.

The motional states $\left|1^{b}_{q}\right\rangle$ of the $\left|b\right\rangle$ atoms are the Bloch waves $\psi_{q}(z)$ with lattice momentum $q\in(-\infty,+\infty)$ which are solutions of the Mathieu equation Arscott et al. (2014). We choose to normalize these to a single Brillouin zone ($BZ$) but work in the extended zone scheme for a more direct connection with the free-particle case. We take their total energy $\hbar\omega_{b,q}$ as the sum of the band energy $\varepsilon_{q}$ (see Fig. 1(b)) of $\mathbf{V}_{b}(z)$ and the internal energy $\hbar\omega_{b}$. We note that for diverging energy, the dispersion relation $\varepsilon_{q}$ approximates that of a free particle subject to the constant average potential $V_{b}/2$; this will become important in Section IV.

If the atom is initially localized in a single well, during short evolution times $t\ll\hbar/J_{b}$ the system is a simulator for spontaneous emission of a photon by an isolated quantum emitter in a photonic crystal de Vega et al. (2008); Stewart et al. (2017a). However, for longer times, the $\left|b\right\rangle$ atom propagates sufficiently far to be reabsorbed by neighboring lattice sites of $\mathbf{V}_{a}$, such that dynamical array effects become noticeable, as was experimentally observed in Krinner et al. (2018); Stewart et al. (2020). These array effects are detrimental to the Markovianity of the system Cascio et al. (2019) and they are especially acute in the ultracold platform due to the strong retardation between emitters.

For such a system of $N$ lattice sites the uncoupled Hamiltonian is given by

where $\hat{a}_{j}^{\dagger}=|1^{a}_{j}\rangle\left\langle 0\right|$ creates an $\left|a\right\rangle$ atom in the $j$-th site from the vacuum $\left|{0}\right\rangle$, $\hat{b}_{q}^{\dagger}=|{1^{b}_{q}}\rangle\left\langle 0\right|$ creates a $\left|b\right\rangle$ atom with lattice momentum $q$, and the explicit sum in modes is $\sum_{q}\equiv\int_{-\infty}^{+\infty}\frac{\mathrm{d}\;\!\!q}{2k}$ for the chosen normalization and zone scheme. The coupling part of the Hamiltonian (in the interaction picture) takes the form

where $\Omega$ is the coupling between the two internal states, $z_{j}=j\pi/k$ is the position of the $j$-th site, and $\gamma_{q}=\left\langle\psi_{q}\left|\right.\phi_{0}\right\rangle$ (see Fig. 1(b)) is the Franck–Condon overlap between the tightly confined wavefunction and the Bloch wave at $q$, and $\Delta=\omega_{\mu}+\omega_{a}-\omega_{b}$ is the detuning of the coupling field Note (1). The Hamiltonian thus takes the form of a system of Weisskopf-Wigner Hamiltonians Weisskopf and Wigner (1930), in which $\hbar\Delta$ plays the role of the excitation energy of a quantum emitter.

Generally, the state of an atom in the state-dependent lattice can be expressed as a linear combination of Wannier states (harmonic-oscillator ground states) and Bloch waves

The time evolution in this picture, $i\hbar{\partial_{t}}|\psi(t)\rangle={\hat{H}}_{ab}|\psi(t)\rangle$, is given by

Following a standard approach in the field Kofman et al. (1994); Lambropoulos et al. (2000); de Vega et al. (2008); Navarrete-Benlloch et al. (2011); Quang et al. (1997); Stewart et al. (2017a) and extending it to a system of sites, the system of differential equations can be solved via Laplace transform $\tilde{A}(s)\equiv\mathcal{L}\{A(t)\}$. Considering processes of spontaneous emission ($B_{q}(t=0)=0$), Eqns. (4) transform into

Substituting the second equation into the first yields

with the bath correlation function

To proceed, it is convenient to express the Bloch waves $\psi_{q}$ entering the calculation of the Franck–Condon factor $\gamma_{q}$ in the basis of plane waves

and to replace the sum over lattice momenta $q$ with an integration over the associated energy $E$. After integration in the complex plane (see Appendix A), this yields

where lattice momentum $q(E)$ is the positive branch of the inverse function of the analytic extension of $\varepsilon_{q}$, $a_{ho}$ is the harmonic oscillator length (i.e. $\phi_{0}(z)=\exp[-z^{2}/(2a_{ho}^{2})]/\sqrt[4]{\pi a_{ho}^{2}}$), ${\rho(E)=2\mathrm{d}\;\!\!q(E)/\mathrm{d}\;\!\!E}$ is the density of Bloch states (DoS), with the prefactor of 2 accounting for the existence of left and right-movers, and $\operatorname{erfc}(x)=2/\sqrt{\pi}\int_{x}^{\infty}\exp(-y^{2}){\mathrm{d}\;\!\!y}$ is the complementary error function.

We note that, because of the Gaussian envelope $\exp\left[{-{(q(E)+k(m+r))^{2}a_{ho}^{2}}}\right]$ the infinite sum in Eqn. (9) converges rapidly, and for a tightly confining emitter ($\lvert a_{ho}q(E)\rvert,\>a_{ho}k\ll 1$), the expression above simplifies greatly. It follows from (8) that $\sum_{m}\upsilon_{m}(E)=1$, and the $\operatorname{erfc}(x)$ in the formula approximates a step function. Hence,

The linear system of equations (6) can be solved for finite $N$ with basic linear algebra (the infinite case $N=\infty$ is treated in detail in Sec. VII), with the solution in matrix form written as $\tilde{\mathbf{A}}(s)=[s\mathbf{I}+\tilde{\mathbf{G}}(i\hbar s+\hbar\Delta)]^{-1}\mathbf{A}(0)$. Using the Wick-rotated variable $E=i\hbar s+\hbar\Delta$ in the inverse Laplace transform then yields

The integral can be further simplified using the residue theorem, for which is important to know the singularities of the integrand. They include poles, i.e. the solutions of the equation

as well as square-root branch points at the band edges $E_{A0}$, $E_{A1}$, $E_{B1}$, $E_{B2}$, etc. introduced in the next section.

Using the Laplace transform to solve the system dynamics (see Eqn. (11)) makes it necessary to evaluate different lattice-related functions in the complex energy plane. This is challenging since some of these functions, such as the lattice momentum $q(E)$, are multivalued and some others, such as the Bloch waves $\psi_{q(E)}(0)$, are physically defined only up to a phase factor. In this section we develop an efficient way for doing so via an infinite-product representation.

We start by analyzing the band structure of the matter-wave vacuum using a standard textbook formula Ashcroft and Mermin (1976),

relating the lattice momentum $q(E)$ to the transmission coefficient $t(E)$ for a plane wave of energy $E$ going through the isolated potential barrier described by $V(z)=\mathbf{V}_{b}(z)$ for $z\in[0,\pi/k]$ and 0 otherwise. Importantly, we will see that $t(E)$ does not need to be evaluated explicitly, but instead just some special points of $T(E)$ will be necessary for the analysis, as they fully determine the physical properties of the system.

These special points, shown in Fig. 2(a), are the energies where $T(E)$ takes zero, unity, or extreme values. In particular, the energies with $\lvert T\rvert=1$ correspond to the band edges of the dispersion relation, which host Bloch waves carrying an integer multiple of the recoil momentum and definite parity. Correspondingly, we label the band edges with even Bloch waves (Mathieu cosines) as $\{E_{An}\}_{n=0}^{\infty}$ and those with odd waves (Mathieu sines) as $\{E_{Bn}\}_{n=1}^{\infty}$ Frenkel and Portugal (2001). Finally, we label the zeroes $\{E_{Cn}\}_{n=1}^{\infty}$ and the extrema $\{E_{Dn}\}_{n=1}^{\infty}$.

The ordering of these energies for different potential depths, shown in Fig. 2(b), is easily understood for free-space motion $V_{b}=0$ and the limit $V_{b}\to\pm\infty$, as the lattice spectrum becomes that of a quantum harmonic oscillator. We note that flipping $V_{b}\mapsto-V_{b}$ (as experimentally done in Stewart et al. (2020), see Fig. 4(c)) leaves the band structure unchanged, but swaps the parity of the edges belonging to every other energy gap. Physically, this transformation is equivalent to displacing the emitter by half a lattice period.

With these definitions, we can efficiently perform an analytical extension of several functions into the complex energy plane by using infinite products. For instance, consider $T(E)$. Since $t(E)\to 1$ as $\lvert E\rvert\to\infty$, the asymptotic expression for this function readily follows from Eqn. (13), and already resembles $T(E)$ quite well (see Fig. 2(a)). A perfect match can be achieved if one “corrects” the zeroes of the approximation by first dividing through them and then multiplying with the actual zeroes of T(E), which gives

This identity is in fact guaranteed by Liouville’s theorem of complex analysis Fine and Rosenberger (2012) (assuming that $T(E)$ has no complex singularities), since the quotient between the two sides of the equation is, by construction, an entire function that tends to 1.

Equation (14) in return allows for the analytic extension both of the lattice momentum, via $q(E)=\frac{k}{\pi}\arccos T(E)$, and of the energy bands $\varepsilon_{n}(q)=T^{-1}_{n}\left(\cos\frac{\pi q}{k}\right)$. In particular, we note that $q(E_{Dn})$ are the branch points where the $n$ and $n+1$ bands cross.

An analogous reasoning can be applied for analytical extensions of other lattice functions such as $\rho(E)$, $\psi_{q(E)}(0)$ and $\psi^{\prime}_{q(E)}(0)$. Since at extreme energies ($\lvert E\rvert\to\infty$) an atom travelling through the lattice potential $\mathbf{V}_{b}(z)$ behaves like a free particle subject to the average constant potential $V_{b}/2$, it is straightforward to find their asymptotic expressions, which can again be corrected further via infinite products. In the case of the density of states (DoS), this results in

For the Bloch waves $\psi_{q(E)}(0)$ taken at the emitter position, the asymptotic value for high energy is $\sqrt{k/\pi}$, and the value has to vanish at the odd band edges $E_{Bn}$ due to parity. Furthermore, Eqns. (7,10) (or equivalently Eqn. (31)) imply that $\rho(E)\psi^{2}_{q(E)}(0)=\frac{2k}{i\pi}\sum_{q}{\lvert\psi_{q}(0)\rvert^{2}}/(\varepsilon_{q}-E)$ and, given that the RHS of this expression is analytic for all $E$ outside the bands, the singularities of $\psi^{2}_{q(E)}(0)$ in this region have to be simple poles matching the zeroes $E_{Dn}$ of the DoS. This leads to the expression

and similarly

Finally, we note that Eqns. (16,17) can be used as initial conditions in the Mathieu equation to obtain the value of the Bloch wave at any other point through $\psi_{q(E)}(z)=\psi_{q(E)}(0)C(E,z)+i\psi^{\prime}_{q(E)}(0)S(E,z)$ where $C$ ($S$) is an entire function in both of their arguments, (anti-)symmetric in $z$ corresponding to the unnormalized Mathieu cosine (sine) function.

The decay of an isolated emitter for the case of coupling to a single, sinusoidal band or a parabolic band edge has been analyzed in earlier theoretical works Stewart et al. (2017a); Lombardo et al. (2014); Calajó et al. (2016). Using our formalism, we now generalize this treatment to include infinitely many (non-)sinusoidal bands, arising for arbitrary lattice depth $V_{b}$. For simplicity, we consider the case $V_{b}>0$ and disregard band gaps beyond an arbitrary cutoff index $\Lambda$ (in practice $\Lambda$ can always be chosen to be small, given that the band gaps quickly get narrower).

As seen in the previous section, the dynamics $A(t)$ of a tightly confining emitter is governed by the bath correlation function $\tilde{G}_{j\text{-}J}(E)\equiv\tilde{G}_{0}(E)$ obtained from Eqns. (10), (15) and (16). To evaluate $\tilde{G}_{0}(E)$, we first introduce the (truncated) products $\Pi_{A}(E)=\prod_{n=0}^{\Lambda}(E-E_{An})$ and $\Pi_{B}(E)=\prod_{n=1}^{\Lambda}(E-E_{Bn})$, as well as their ratio $\Pi_{B/A}(E)=\Pi_{B}(E)/\Pi_{A}(E)$. Using these products, the bath correlation function can then be approximated as $\hbar\tilde{G}_{0}(E)\approx\kappa\sqrt{\Pi_{B/A}(E)}$, with the coupling constant $\kappa=(\Omega/2)^{2}\hbar a_{ho}\sqrt{2\pi m}$. To avoid ambiguity, we consider that all square roots $\sqrt{\cdots}$ in this section give back complex numbers with argument in $(-\pi/2,+\pi/2]$.

To find the poles of the inverse Laplace transform  (11), we multiply Eqn. (12) with its algebraic conjugate; the poles then correspond to zeroes of the polynomial $(E-\hbar\Delta)^{2}\Pi_{A}(E)+\kappa^{2}\Pi_{B}(E)$. By considering the degree of this polynomial and its changes in sign between band edges, it is easy to see that there are no bound states in the continuum (BIC) Hsu et al. (2016), whereas each band gap (including $E<E_{A0}$) contains at least one pole. This differs from the multi-emitter case $1<N<\infty$, where subradiance and retardation effects can lead to a BIC Sinha et al. (2020); Calajó et al. (2019). For sufficiently weak couplings ($\kappa\ll E_{r}^{3/2}$), these in-gap poles can be approximated as $E_{n}\approx E_{An}-\kappa^{2}\Pi_{B}(E_{An})\>/\>[(E_{An}-\hbar\Delta)^{2}\prod_{m\neq n}(E_{An}-E_{Am})]$. There are two additional poles, approximately $\hbar\Delta\mp i\kappa\sqrt{\Pi_{B/A}(\hbar\Delta)}$, which lie in the same gap if they are real; otherwise one of them ($E_{M}$) has negative imaginary part and can lead to Markovian (exponential) decay of the emitter, and the other ($E_{M}^{*}$) is its complex conjugate.

In determining the spectral decay properties, see Fig. 3, we see that not all of these poles contribute towards the residue theorem: due to the square root singularities at the band edges, we can visualize the integrand domain as a Riemann surface consisting of an ‘upper sheet’ where the integration paths in the complex plane are located and a ‘lower sheet’ on the other side of the branch cuts. Only the poles on the upper sheet will contribute towards the residue theorem and have a physical interpretation.

For positive lattice depths $V_{b}>0$ and Markovian couplings $\Omega\ll\min_{q\in\mathbb{R}}|\Delta-\varepsilon_{q}/\hbar\rvert$, the pole $E_{n}$ is in the upper sheet only when $\operatorname{sign}\{\hbar\Delta-E_{An}\}=(-1)^{n}$. Of the two extra poles, only one is in the upper sheet; in particular $E_{M}^{*}$ cannot be there because its positive imaginary part would lead to exponential growth of the population (see Eqn. (18)).

For larger couplings, there is still a change of sheets for one of the poles living in the $n^{th}>0$ gap as $\hbar\Delta$ crosses the value $E_{Bn}$. On the other hand, by increasing the coupling it is possible to make two lower poles co-located in a gap combine into a double pole and then split into a Markovian pole $E_{M}$ and its conjugate $E_{M}^{*}$. Whereas $E_{M}^{*}$ remains always unphysical, $E_{M}$ can make it to the upper sheet, as depicted in Fig. 3(c).

This figure reveals behavioral differences between the decay next to a band edge hosting even Bloch waves and one hosting odd ones. Whereas the former behaves as expected from previous studies Stewart et al. (2017a), the latter displays an increase in the Markovian component of the decay at non-Markovian couplings. One would naively expect that reabsorption and emission scale equally with the vacuum coupling; however this is not the case here as BS formation is suppressed for these parameters (see Fig. 3(b)). Despite the phenomenological similarities, this ultra-Markovian emission is not to be confused with superradiance, as a single emitter is enough to create this effect. For illustration purposes, let us consider the following example. Under the same conditions as Fig. 3(b,c), a quantum emitter with detuning $\hbar\Delta=E_{B1}+\tfrac{1}{5}E_{r}$ and coupling $\hbar\Omega=\tfrac{5}{2}E_{r}$ emits half of its population at only $t=0.52\hbar/E_{r}$ and it emits more than $90\%$ ultimately. In contrast, while the initial decay is very similar at $\hbar\Delta=E_{A0}+\tfrac{1}{5}E_{r}$ and $\hbar\Omega=\tfrac{5}{2}E_{r}$, most of the radiation is reabsorbed at $t=2.5\hbar/E_{r}$ and the red population subsequently oscillates with large amplitude as the system cycles between emission and reabsorption.

The integration path of Eqn. (11) can be adapted to the singularities described in Fig. 3(a) by circling around the physical poles and branch cuts of the integrand Note (2), leading to the time evolution

where

denotes the residue of the poles and the branch contribution

is well defined unless $E_{M}$ stands on the branch cut. The branch contribution $I\left(E,t\right)$ tends to zero as the time $t$ increases but in a non-exponential fashion, making the decay non-Markovian. The only contributions persistent in time, i.e. the bound states (BS), are the real upper poles residing in the band gaps.

The general features of the resulting time evolution are shown in Fig. 4(a). It is mostly Markovian for detunings deep inside the bands and non-Markovian around the band edges; there is no decay deep inside the gaps. We note that that emission is suppressed at the band edges $E_{Bn}$ of odd parity with respect with the edges $E_{An}$, whose Bloch waves have the same parity as the emitter.

After the emission, the emitter population can oscillate by the beating of various bound states. This effect, which has also been measured experimentally Stewart et al. (2020), is most noticeable in the center of the first energy band. We compare our calculation with the experimental data in Fig. 4(b).

Although the one-emitter model ($N=1$) matches the observed decay dynamics at short times, it underestimates the amount of subsequent reabsorption seen in the experiment, in which the optical lattice provided an array of emitters. As already discussed in Stewart et al. (2020), the subsequent oscillations are dominated by reabsorption as the emitted radiation spreads across the emitter array. Using the formalism developed in this paper, the presence of neighboring ground-state emitters (i.e. empty lattice sites) surrounding an excited emitter can readily be taken into account in Eqn. (11) via an analogous approach, and already working with a lattice of three sites ($N=3$) shows a marked improvement for the second oscillation. While it generally gets harder to analyze and numerically solve larger arrays, a special case is $N=\infty$, analyzed in Sec. VII. As seen in the figure 4(b), the overall agreement with the experiment qualitatively improves further to longer time scales, but residual deviations persist. They are likely due to differences in the initial state (the experiment worked with a sparsely populated array with more than one excitation) and collisions between atoms.

An analogous formula to (18) can be written for the time evolution $B_{q}(t)$ of the emitted modes,

The main difference is that the integration displays an additional pole at energy $\varepsilon_{q}$. Importantly, this real pole does not correspond to a bound state (discussed further in the next section), but to a mode that has completely abandoned the emitter and keeps travelling free indefinitely. This causes similarity between the emission spectrum and the dispersion relation of the medium (see Fig. 4(c)). We note that at smaller couplings, the applicability of the $N=1$ theory extends to longer times, which was the case for the parameters of Fig. 4(c). More generally, a multi-emitter analysis of the emitted modes is accessible through Eqn. (4).

As seen in the previous section, the long-time dynamics of the quantum emitter is dictated by the presence of bound states. BS have been broadly studied in the literature Shi et al. (2016) and in 1D are limited to be short ranged Sánchez-Burillo et al. (2020), which is the reason why they are often depicted as decaying exponentially in space Stewart et al. (2017a); Shi et al. (2018, 2016); de Vega et al. (2008); Calajó et al. (2016) while the detailed features are often overlooked. However, the reader might expect that the presence of a lattice potential for the radiated modes induces a corrugation of this evanescent wave Bykov (1975). Even more strikingly, recent experimental work Stewart et al. (2020) found that time-of-flight distributions of these BS can possess two sharp peaks at opposite momenta, suggesting common features with stationary waves despite the absence of boundaries. In this section we clarify these apparent inconsistencies by computing the exactly spatial distribution of the BS and presenting a simple physical picture that encompasses all these effects.

The BS wavefunction $\left|\psi\right\rangle(t)$ is defined in the interaction picture (see Eqn. (3)) by

where $A(0)=\left[1+\sum_{q}\lvert\gamma_{q}\rvert^{2}\left(\frac{\hbar\Omega/2}{E_{BS}-\varepsilon_{q}}\right)^{2}\right]^{-1/2}$ is the normalization constant. The red part of this expression follows from isolating the contribution of a single pole in Eqn. (18), whereas the blue part and the equation (12) for the BS energy $E_{BS}$ follow from integrating (4).

The spatial distribution of this bound state is better understood in the Schrödinger picture $\left|\psi(t)\right\rangle_{S}=e^{-i\hat{H}_{0}t/\hbar}\left|\psi(t)\right\rangle$, where position is a time-independent operator. The blue part of $\left|\psi(t)\right\rangle_{S}$, has a position-space wave function $B_{S}(z,t)=\sum_{q}B_{S,q}(t)\psi_{q}(z)$ which can be integrated in a very similar way to the derivation of (9), leading to

This expression might seem cumbersome, but we note that it is real at $t=0$ and symmetric; that the prefactor is regular also when $E_{BS}=E_{Dn}$, due to the cancellation of the zero of $\rho(E_{BS})$ with the pole of $\psi_{q(E_{BS})}^{2}(0)$; that the time evolution is that of an eigenstate of the system (even though the Hamiltonian is not time-independent due to the external coupling, which causes the red part of the BS to have a frequency lower than the blue by $\omega_{\mu}$); and that the dominant modes of the sums are $\upsilon_{0}(E_{BS})$ and $\upsilon_{-n}(E_{BS})=\upsilon_{0}^{*}(E_{BS})$ (with $n$ numbering the gap of $E_{BS}$) when the blue lattice is shallow as they correspond with the momenta that are closest to the free-particle dispersion relation (see Fig. 5). The BS is localized at the center of the site and decays to both sides with an inverse exponential decay length $\operatorname{Im}q(E_{BS})$ and a wavenumber $\operatorname{Re}q(E_{BS})=nk$, and it is consistent with those presented in Refs. Stewart et al. (2017a); Bykov (1975). In summary, expression (23) shows that the BS is a linear combination of evanescent waves whose momenta match the analytic extension to the band gaps of the blue lattice dispersion relation, as presented in Fig. 5.

Looking back at the spectrum of the system described in the previous section, these bound states may be created adiabatically by choosing a detuning $\Delta$ resonant with the $n^{th}$ band gap and increasing the vacuum coupling strength $\Omega$ adiabatically, as empirically shown in Refs. Krinner et al. (2018); Stewart et al. (2020). The adiabatic condition consists of the change in these parameters being slower than the time scale defined by the energy difference between the instantaneous BS energy and the nearest energy edge. This imposes two restrictions for adiabatic excursions in the $\Delta$-$\Omega$ plane (see Fig. 3(b)): $\Omega$ must not be brought back to 0 while $\Delta$ is in some region other than the $n^{th}$ band gap since the BS energy would cross $E_{An}$, and $\Delta$ must not cross the corresponding odd band edge $E_{Bn}/\hbar$ since the BS energy would also cross $E_{Bn}$.

The procedure of adiabatic creation allows for testing experimentally the BS probability distributions derived in this section. The lattice-momentum distribution $\lvert B_{q}(t)\rvert^{2}=\lvert B_{S,q}(t)\rvert^{2}$ of Eqn. (22) matches with the band-map measurements of Refs. Krinner et al. (2018) (for the case $V_{b}=0$; see also Stewart et al. (2017a)) and Ref. Stewart et al. (2020) (case $V_{b}=2.5E_{r}>0$). Furthermore, signatures of the real-space distribution $\lvert B_{S}(z,t)\rvert^{2}$ of Eqn. (23) might be accessed by highly resolving in situ imaging techniques Bakr et al. (2009).

In general, as the number of emitters in the array increases, the spectrum of the system becomes more elaborate, but the case of infinitely many emitters is an exception to this. As $N\to\infty$, the Hamiltonian regains the lattice periodicity and states with different lattice momentum $q$ decouple via Bloch’s theorem. In our formalism this is reflected in the fact that vectors of the form $\mathbf{A}(0)=(\ldots,e^{iq\pi(j-1)/k},e^{iq\pi j/k},e^{iq\pi(j+1)/k}\ldots)^{\operatorname{T}}$ are eigenstates of $\tilde{\mathbf{G}}(E)$ while being independent of $E$. We will refer to these states as matter-wave polaritons J. Kwon et al. (2021), in analogy to of their quantum optical counterparts Carusotto and Ciuti (2013); Shi et al. (2018); Basov et al. (2021). While photons are also involved in this new type of quasiparticle by constituting the optical lattice, the usual role of a photon in a polariton is taken over by matter waves. In the regime $\lvert V_{b}\rvert\ll E_{r}\ll\lvert V_{a}\rvert$, the dispersion relation of the radiated modes is nearly quadratic and matter-wave polaritons are similar to exciton polaritons Carusotto and Ciuti (2013); Schneider et al. (2017), whereas for $E_{r}\ll\lvert V_{b}\rvert\ll\lvert V_{a}\rvert$ the now strongly confining optical lattice for blue atoms forms the analogue of a coupled cavity array in circuit QED Hartmann et al. (2006); Greentree et al. (2006); Noh and Angelakis (2016); Hartmann (2016); Fitzpatrick et al. (2017); Ma et al. (2019).

Taking the Fourier transform $\mathcal{F}\hskip-2.0ptA_{\hskip 1.0ptq}=\sum_{j}A_{j}e^{-i\pi jq/k}$ as a change to the decoupled basis, Eqn. (11) simplifies to

where the eigenvalue

follows immediately by applying the definition of eigenvalue to (7) and using the identity $\sum_{j=-\infty}^{+\infty}e^{i(q^{\prime}-q)\pi j/k}=2k\sum_{n=-\infty}^{+\infty}\delta(q^{\prime}-(q+2kn))$. The integrand has thus become a meromorphic real function, free of the branch cuts at the band edges and free of complex poles. This indicates that radiation cannot escape the emitter array, which is to be expected since the array is infinite.

Since all of the poles are real, we follow a suggestion in  Bykov (1975) and visualize Equation (25) as shown in Fig. 6(a), in order to locate all of the solutions $\mathcal{E}_{n}(q)$ that physically correspond to the polariton energy bands. Some simple properties that follow from this graph are that polariton bands neither cross each other nor the original energy bands ($\mathcal{E}_{1}(q)<\hbar\Delta$, $\varepsilon_{n-1}(q)<\mathcal{E}_{n}(q)<\varepsilon_{n}(q)$ for $n\geq 2$), although they might cross the detuning at the points $\tilde{g}_{q}^{-1}(0)$ where the couplings to different bands cancel mutually. Excited polariton bands tend to these points in the limit of very large coupling; otherwise they soon approximate the energy bands ($\mathcal{E}_{n\gg 1}(q)\approx\varepsilon_{n-1}(q)$).

The resulting band structure (purple lines in Fig. 6(b)) is exotic and cannot be obtained by a simple periodic potential. An indicator of this is the positive effective mass of both the ground and the first excited polariton band near $q=0$.

The motional properties associated with these bands follow from applying the residue theorem to (24), which leads to

where the residues are given by

The residues for bands far away from the detuning are negligible, rendering the structure of higher bands/gaps irrelevant for the dynamics. Alternatively, these results may also be obtained by decoupling the Hamiltonian Shi et al. (2018) or as an ansatz in the Schrödinger equation, while taking into consideration that these residues satisfy the normalization condition ($r_{n}\in[0,1]$ and $\sum_{n=1}^{\infty}r_{n}(q)=1$), energy conservation ($\sum_{n=1}^{\infty}r_{n}(q)\mathcal{E}_{n}(q)=\hbar\Delta$) and $\sum_{n=1}^{\infty}r_{n}(q)/(\varepsilon_{m}(q)-\mathcal{E}_{n}(q))=0$. With this solution of the system, we can explain the dynamical behaviour of the experimental data Krinner et al. (2018); Stewart et al. (2020) at longer evolution times $t\geq\hbar/J_{b}$ (see Fig. 4(b), case $N=\infty$).

Moreover, as different quasi-momenta are decoupled unless they differ by an even multiple of the recoil momentum $k$, one can define a periodic momentum-dependent detuning $\Delta\equiv\Delta(q)$ to account more accurately for the dispersion relation $\varepsilon_{n}^{\scalebox{0.5}{$(a)$}}(q)$ for the $\left|a\right\rangle$ states in the $n^{\text{th}}$ band, while also modifying the Franck–Condon overlap into $\gamma_{q}=\left\langle\psi^{\scalebox{0.5}{$(b)$}}_{q}\left|\right.\psi^{\scalebox{0.5}{$(a)$}}_{n,q}\right\rangle$. This allows more customization of the resulting polariton bands $\mathcal{E}_{n^{\prime}}(q)$, whose hopping rates

can be freely tuned. An extreme example of this is shown in Fig. 6(c,d) where the first band of the $\left|b\right\rangle$ states is coupled with the second band of the $\left|a\right\rangle$ in a way that the polaritons hop two lattice sites at a time, without going through the intermediate site at all ($\bar{J}_{1}=0$, $\bar{J}_{2}\neq 0$). This opens the possibility for the analog simulation of $J_{1}$-$J_{2}$ quantum spin models Chandra and Doucot (1988); Dagotto and Moreo (1989) with ultracold bosons in 1D, by inducing effective spin interactions Kuklov and Svistunov (2003); Altman et al. (2003); Liao et al. (2021). We note that frustration Balents (2010), a core feature of the $J_{1}$-$J_{2}$ model, in our system has kinetic origin and is generated by coupling bands of opposite effective mass Sträter and Eckardt (2015).

We have shown the importance of having multiple emitters coupled to multiple band modes in ultracold-atom realizations of waveguide QED. Already on the single emitter level, the intricate band structure enriches the spectrum, partially due to the presence of energy edges hosting odd Bloch waves (which present an unexpected enhancement of the Markovian part of the decay in the non-Markovian regime), and partially due to the presence of multiple bound states in different band gaps. Our analysis was made possible by what to our knowledge are novel, simple analytical expressions for multiple well-known functions describing the dynamical properties of a particle in a lattice potential. The fast convergence of these expressions and their applicability on the whole complex energy plane make them relevant in contexts beyond the scope of this work, although there is no trivial generalization to higher dimensional lattices. And finally we have studied polariton formation in the matter-wave context, showing that this system not only can act as an analog simulator of photonic phenomena, but also as a wider platform for studying low-dimensional frustrated systems.