Non-Markovian dynamics of giant emitters beyond the Weisskopf-Wigner approximation

TL;DR

This paper investigates the complex non-Markovian behavior of giant quantum emitters beyond traditional models, providing exact solutions and insights into bound states and superpositions, with potential superconducting circuit implementations.

Contribution

It introduces minimal models of giant emitters coupled to waveguides, solving their dynamics exactly and exploring bound states and superpositions beyond the Weisskopf-Wigner approximation.

Findings

Exact solutions for symmetric giant emitter configurations

Identification of bound states inside and outside the continuum

Feasible implementation in superconducting circuits

Abstract

Giant quantum emitters, whose effective size is comparable to the wavelength of the radiation they couple to, give rise to interference effects and non-Markovian dynamics that lie beyond the scope of the standard Weisskopf-Wigner approximation. Here we study minimal models of giant emitters coupled to cavit$-array waveguides and identify symmetric configurations where the dynamics can be solved exactly. This allows us to capture the emergence of bound states both inside and outside the photonic continuum, and to demonstrate the possibility of engineering their coherent superpositions. Our results provide analytical insight into non-Markovian light-matter interaction and suggest a feasible implementation using superconducting circuit platforms.

Tables

Table 1. Table 1: Experimental parameters of the circuit model taken from Ref.  [ 9 ] . We have set the qubit–line coupling to g 0 ∼ ξ / 5 = 49.8 ​ MHz g_{0}\sim\xi/5=49.8~\mathrm{MHz} for all the calculations presented in the paper.

Parameter	Description	Value
${\omega }_0/2\pi$	Cavity resonance frequency	$5.71\mathrm{GHz}$
$\xi /2\pi$	Hopping strength	$249\mathrm{MHz}$
$g_0/2\pi$	Qubit cavity- array coupling strength	$50\mathrm{MHz}$

Equations

$$H=\delta {\sigma }^{+}{\sigma }^{-}+{\omega }_0\sum _n{a}_n^{†}{a}_n+\xi \sum _n\left(a_{n+1}^{†}{a}_n+a_n^{†}{a}_{n+1}\right)+g_0\left({\sigma }^{+}{a}_0+{\sigma }^{-}{a}_0^{†}\right),$$

$$H=\Delta {\sigma }^{+}{\sigma }^{-}+\sum _k{\omega }_k{a}_k^{†}{a}_k+\frac{g_0}{\sqrt{N}}\sum _k\left({\sigma }^{+}{a}_k+{\sigma }^{-}{a}_k^{†}\right),$$

$$i\dot{\tilde{\alpha }}$$

$$i{\dot{\tilde{\beta }}}_k$$

$$\dot{\tilde{\alpha }}(t)=-\frac{1}{2\pi }{\int }_0^{t}d\tau {\int }_{{\omega }_{\min }}^{{\omega }_{\max }}d\omega J(\omega )e^{-i(\omega -\Delta )\tau }\tilde{\alpha }(t-\tau ),$$

$$J(\omega )=2N|g_k(\omega ){|}^2{\left|\frac{\partial {\omega }_k}{\partial k}\right|}^{-1},$$

$$J(\omega )=\frac{2g_0^2}{\sqrt{4{\xi }^2-{\omega }^2}}.$$

$$\alpha (t)=e^{-(\gamma /2+i\Delta )t},$$

$$\tilde{\alpha }(t)=\frac{e^{i\Delta t}}{2\pi i}{\int }_{\sigma -i\infty }^{\sigma +i\infty }ds\frac{e^{st}}{s+i\Delta +G(s)},$$

$$G(s)=\frac{1}{2\pi }{\int }_{{\omega }_{\min }}^{{\omega }_{\max }}d\omega \frac{J(\omega )}{s+i\omega }.$$

$$\alpha (t)$$

$$y-\Delta +iG(-iy)=0.$$

$$G(s)=\frac{g_0^2}{s\sqrt{1+4{\xi }^2/s^2}}.$$

$$\underset{ϵ\to 0}{\lim }G(\pm ϵ+iy)=\pm \frac{g_0^2}{\sqrt{4{\xi }^2-y^2}},$$

$$\alpha (t)=\frac{4g_0^2}{\pi {\xi }^2}{\int }_{-1}^{1}d\tilde{y}K(\tilde{y})e^{2i\xi \tilde{y}t}+r_{+}{e}^{-i{y}_{+}t}+r_{-}{e}^{-i{y}_{-}t},$$

$$K(\tilde{y})=\frac{\sqrt{1-{\tilde{y}}^2}}{4(2\tilde{y}+\Delta /\xi {)}^2(1-{\tilde{y}}^2)+(g_0/\xi {)}^4},$$

$$r_{\pm }={\left[1+\frac{g_0^2}{y_{\pm }^2(1-4{\xi }^2/y_{\pm }^2{)}^{3/2}}\right]}^{-1}.$$

$${\beta }_n(t)=-i^{|n|+1}{g}_0{\int }_0^{t}d\tau \alpha (t-\tau )e^{-i{\omega }_0\tau }{J}_{|n|}(2\xi \tau ).$$

$$\underset{t\to \infty }{\lim }|{\beta }_n^{\pm }(t){|}^2=\text{expectationvalue}{a_n^{†}{a}_n}_{\pm }=\frac{g_0^2{r}_{\pm }^2}{y_{\pm }^2-4{\xi }^2}{e}^{-2{\kappa }_{\pm }\text{absolutevalue}n},$$

$$H=\delta {\sigma }^{+}{\sigma }^{-}$$

$$H=\Delta {\sigma }^{+}{\sigma }^{-}+\sum _k{\omega }_k{a}_k^{†}{a}_k+\sum _k{\tilde{g}}_k\left({\sigma }^{+}{a}_k+{\sigma }^{-}{a}_k^{†}\right),$$

$${\tilde{g}}_k=\frac{g_0}{N_c\sqrt{N}}\frac{\sin (kd{N}_c/2)}{\sin (kd/2)}.$$

$$J_{\mathrm{eff}}(\omega )=J(\omega )\mathcal{G}(\omega ),\mathcal{G}(\omega )=\frac{1}{N_c^2}\frac{1-\cos (k(\omega )d{N}_c)}{1-\cos (k(\omega )d)}.$$

$$J_{\mathrm{eff}}(\omega )=\frac{g_0^2}{\sqrt{4{\xi }^2-{\omega }^2}}\left[1+\cos (k(\omega )d)\right].$$

$$G(s,d)=\frac{g_0^2}{2}\frac{\mathrm{sign}(\Re \{s\})}{s\sqrt{1+4{\xi }^2/s^2}}\left[1+{\varrho }^d\right],$$

$$\underset{ϵ\to 0}{\lim }G(ϵ+iy,d)-\underset{ϵ\to 0}{\lim }G(-ϵ+iy,d)$$

$$\alpha (t)=\frac{4g_0^2}{\pi {\xi }^2}{\int }_{-1}^{1}d\tilde{y}{K}_{\mathrm{giant}}(\tilde{y},d)e^{i2\xi \tilde{y}t}+\sum _{j=\pm }{r}_j{e}^{-i{y}_{j}t},$$

$$K_{\mathrm{giant}}(\tilde{y},d)=\frac{\sqrt{1-{\tilde{y}}^2}\left[1+\cos (k(-\tilde{y})d)\right]}{8(2\tilde{y}+\Delta /\xi {)}^2(1-{\tilde{y}}^2)+F(\tilde{y},d)},$$

$$\alpha (t)=\frac{4g_0^2}{\pi {\xi }^2}{\int }_{-1}^{1}d\tilde{y}{K}_{\mathrm{giant}}(\tilde{y},d)e^{i2\xi \tilde{y}t}+r_m{e}^{-i{\omega }_{m}t},$$

$$\text{expectationvalue}{a_n^{†}{a}_n}_{\text{BIC}}=\frac{g_0^2{r}_m^2}{2(4{\xi }^2-{\omega }_m^2)}\left[1+\cos (\varphi (|n+d/2|-|n-d/2|))\right],$$

$$\underset{t\to \infty }{\lim }|\alpha (t){|}^2=r_{+}^2+r_m^2+2r_{+}{r}_m\cos (y_{+}-{\omega }_m)t,$$

Full text

I Introduction

The theory of open quantum systems has been applied successfully to a broad variety of problems in quantum optics, solid state physics, quantum information and computation, and even to some dynamical processes in biological systems [1, 2, 3]. This description often relies on the assumption of having a memoryless environment which allows the application of the so called Markovian approximation. A key quantity in this approximation is the system-environment spectral function that can be frequency-dependent in general. In the context of light-matter interaction in free space, this function is almost constant in the relevant range of optical frequencies, a fact that justifies the use of the Markovian approximation. The first theoretical approach to explain the process of spontaneous emission based on this assumption was provided by Weisskopf and Wigner in 1930 [4]. This formalism works perfectly fine and describes the irreversible exponential decay of atoms in most situations in free space. The situation is quite different when we deal with * structured* environments, showing non-trivial frequency behaviors of the spectral function [5]. This is the case of quantum emitters coupling to photonic crystals [6], plasmonic waveguides [7], atomic waveguides [8], transmission line superconducting circuits [9, 10], among others. The dynamics of the emitters turns non-Markovian and a treatment beyond the Weisskopf-Wigner approach is needed. As an example of the breakdown of the Markovian dynamics is the existence of bound states outside the continuum (BOC) when tuning the frequency of an emitter within a band gap [11, 5, 12]. This leads to a partial (non-exponential) emission of the radiation into the medium, exponentially localizing the photonic wavefunction at the position of the emitter. In this case, the spectral function has two van Hove singularities at the band-edges, making the spectral function highly sensitive to frequency variations and therefore not amenable for a Markovian treatment. In order to have an appropriate description of a structured environment it becomes necessary to go beyond the standard Markovian treatment and to take into account the complexity of a realistic spectral function. This is indeed possible for some experimentally relevant one-dimensional structured waveguides as we will see later in this paper.

In recent years artificial atoms known as giant atoms have brought a lot of attention due to their interesting interference properties [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55]. These artificial atoms are considered giant since their effective size is comparable to the wavelength of the electromagnetic radiation they can couple to. As a consequence of this fact they can no longer be treated as point-like matter, and therefore the common dipole approximation becomes invalid. In the most direct and standard experimental implementation of these systems using superconducting circuits they can be realized by coupling a transmon qubit to meandering coplanar waveguides [56]. These transmon qubits are coupled to distant points along the waveguide, so there can be significant phase differences and time delays of the propagating radiation along the waveguide. This situation is highly non-Markovian and gives rise to revivals of the emitter population, frequency-dependent rates, and to the creation of bound states inside the continuum (BIC). BIC states arise from interference and have being experimentally observed in a broad variety of physical systems [57]. A very interesting result is the ability of a giant atom to create persistently oscillating bound states as consequence of the overlap of two bound states inside the continuum [58], a fact that could find potential applications in the encoding and manipulation of quantum information in cavity-QED.

Inspired by the work of Berman and Ford [59], in this paper we study the non-Markovian dynamics of a minimal model for a giant emitter coupled to a cavity-array waveguide beyond the Weisskopf-Wigner approach. Our treatment extends the applicability of the formalism to current interesting and experimentally relevant systems in the context of giant atoms. The paper is organized as follows: In section II we describe the main model employed throughout the paper. It consists of a single emitter coupled to cavity-array waveguide. In section III and IV we describe the general formalism based on the exact inversion of the Laplace transform and its application to the basic model. In section V we show that these methods can be extended to the case of giant emitters in waveguide-QED and present an exactly solvable minimal model. In section VI a full microscopic derivation of the Hamiltonian based on a lumped-element circuit model is given, and finally in section VII we outline the conclusions.

II Model

Let us consider a two-level system with excited state energy $\delta$, coupled to a one-dimensional array of cavities acting as an effective waveguide. The Hamiltonian of this system reads as follows,

[TABLE]

where $\omega_{0}$ is the on-site cavity energy and $\xi$ is the hooping parameter between nearest-neighbor cavities. The two-level system couples to the cavity array trough a single cavity located at the center of the waveguide at $x_{0}=0$. The cavity-array waveguide can be diagonalized by taking advantage of the translational invariance trough the discrete Fourier transform of the bosonic operators, i.e., $a_{k}=N^{-1/2}\sum_{n}e^{ikn}a_{n}$. Applying this transformation to the full Hamiltonian (1) and moving to a rotating frame with respect to the natural frequency of the cavities $\omega_{0}$, we obtain,

[TABLE]

where $\Delta=\delta-\omega_{0}$, and $\omega_{k}=2\xi\cos k$ is the dispersion relation of the cavity-array waveguide. An initially excited two-level system will emit radiation into the waveguide modes as a consequence of the spontaneous emission. In the single-excitation manifold we can describe this process by computing the probability amplitudes associated to finding the atom in the excited sate $\ket{e}$ and the cavity array in the collective vacuum $\ket{\mathbf{0}}=\ket{0...0}$, which we call $\alpha$, and the one corresponding to the atom in the ground state $\ket{g}$ and a photon occupying the $k-$th mode $a^{\dagger}_{k}\ket{\mathbb{0}}$, which we name $\beta_{k}$. Using the Scrhödinger equation we get the following coupled differential equations for the probability amplitudes,

[TABLE]

where we have made the following change of variables: $\tilde{\alpha}=e^{i\Delta t}\alpha$ and $\tilde{\beta}_{k}=e^{i\omega_{k}t}\beta_{k}$. We can decouple eqs.(II) and write down the following integro-differential equation for the atomic probability amplitude as

[TABLE]

where we have introduced the spectral function $J(\omega)\equiv 2\pi\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})$, that characterizes the effects of the bath of harmonic oscillators in Fourier space on the dynamics of the two-level system [1, 3]. Here $g_{k}=g_{0}/\sqrt{N}$, which is independent of $k$. We can explicitly compute the spectral function in the continuum limit by taking $\sum_{k}\to(N/2\pi)\int^{\pi}_{-\pi}dk$ in the definition of the spectral function. This leads to the following expression [3]

[TABLE]

provided that the dispersion relation $\omega_{k}$ has a single branch (see inset figure in fig.1(b)). The spectral function for this case reads,

[TABLE]

In fig.1(b) we show the behavior of the spectral function as a function of the frequency. This spectral function has singularities at the band-edge frequencies $\pm 2\xi$, and it is almost constant for most frequencies near the center of the band dispersion. As we can see from the expression in eq.(4), near the center of the band the function $J(\omega)$ varies slowly with frequency and we have that $J(\omega)\approx J(0)$. We can then write the usual purely exponential decay expression for the atomic probability amplitude as

[TABLE]

being $\gamma=J(0)=g_{0}^{2}/\xi$ the decay rate in the Markovian limit. This approximation is obviously not valid for frequencies close to the band edges, and we must use an exact treatment in order to describe the actual dynamics.

III Exact dynamics for a general spectral function

In this section we briefly review the method employed by Berman and Ford [59] to characterize the dynamics of decay of a two-level system into a continuum of modes. This method allows us to find explicit solutions for the probability amplitudes that go beyond the Wigner-Weisskopf approximation, providing a systematic way of dealing with arbitrary spectral functions. This is specially well suited for the study of non-Markovian or memory effects naturally present in the study of structured baths described by nontrivial spectral functions. Eq.(4) can be formally solved using the Laplace transform. By doing so and taking the formal inverse Laplace transform we can write the formal solution as

[TABLE]

where we have defined the function

[TABLE]

It can be shown that the function $G(s)$ has a discontinuous jump along the imaginary axis, from $-\omega_{\rm{max}}$ to $-\omega_{\rm{min}}$, i.e, there exists a branch cut along this interval. Choosing the integration contour depicted in fig.2(a) and applying the residue theorem we arrive at the following exact solution

[TABLE]

where r_{j}=\left[1+G^{\prime}(s)\big{|}_{s=-iy_{j}}\right]^{-1}, with the poles $y_{j}$ being solutions to the equation

[TABLE]

The existence of such poles depends entirely on the function $J(\omega)$. In the case of the cavity array treated here there are two poles that lie outside the branch cut, bellow and above the branch points. The solution in eq.(III) shows two contributions. The first term represents the contribution from scattering states inside the continuum (fig.1 (c)). In the long time limit this term decays to zero and does not contribute to the steady state. The second term however, obtained from the residues, gives rise to a finite contribution in the steady state. This is always the case if the spectral function $J(\omega)$ is bounded, either from above or bellow [59]. In the case of the the spectral function given in (6) there are two bound states, one bellow and one above of the continuum of states, corresponding to two distinct solutions of the equation of eq.(11). The energy of these states is plotted in fig.1(c) as a function of the coupling where we can see that they move farther away from the continuum as we increase the coupling to the waveguide. In the following section we illustrate the use of this general method and solve the dynamics for the cavity array bath and its spectral function exactly.

IV Exact solution for a cavity-array waveguide

In order to apply the formal solution in (III) we employ the explicit form of the spectral function given in eq.(6), together with the expression for $G(s)$ in this case, which reads

[TABLE]

As mentioned above this function has a branch cut along the imaginary axis, within the interval $[-2\xi,2\xi]$. Using the fact that

[TABLE]

we can obtain the following expression for the excited state atomic probability amplitude,

[TABLE]

where the integral kernel is given by,

[TABLE]

and the residues

[TABLE]

Notice that we have made the change of variable $\tilde{y}=y/2\xi$. The poles $y_{\pm}$ are real solutions to eq.(11) and represent the energy of two bound states outside the continuum band. This solution has been equivalently found in a different way in Ref.[11] by spanning the state vector into bound and scattering eigenstates. One nice thing about the present method is that there is no need to separate the contributions coming from the continuum and the bound modes, since they appear naturally as consequence of the integration around the branch cut and the residues of each pole inside the contour, respectively. In fig.2(b) we show the decay probability of the two-level system as a function of time varying the qubit frequency from $-5\xi/2$ to $5\xi/2$. We can see that if the emitter is on resonance (or close to) with the center of the band the dynamics is indeed very well described by an exponential decay, corroborating the Markovian approximation. Oscillations in the decay probability start to appear as we get detuned from resonance. Close to the band-edges we start to gain some contribution from the bound state (the state corresponding to $y_{+}$ in this case), and the decay becomes partial. This means that part of the excitation is trapped around the position of the emitter, a phenomenon that is well known in this type of systems and has been well studied in the literature of waveguide-QED [12, 11]. The non-zero asymptotic value of the decay probability is given exactly by the square of the residue $r^{2}_{+}$ (dotted lines in the main plot of fig.2(c)). Fig.2(c) also shows the exponentially localized wave functions in position space of the radiation field and their asymmetry as we move towards the upper band limit. In this case the value of the positive pole (energy of the upper bound state) increases as well as the localization length of its wave function.

The probability amplitude of the radiation field in position space can be obtained from the Fourier transform of the one in momentum space, i.e., $\beta_{n}=(1/\sqrt{N})\sum_{k}\beta_{k}e^{-ikn}$. Using the expression for $\beta_{k}$ in eq.(II) we can write,

[TABLE]

For values of $\Delta$ closer to the upper/lower band limit it is possible to compute the steady state of the radiation field as,

[TABLE]

where $\kappa_{\pm}=\lambda_{\pm}^{-1}=\cosh^{-1}\left({|y_{\pm}|/2\xi}\right)$, being $\lambda_{\pm}$ the corresponding localization length for each bound state, respectively. For values of $\Delta$ not that close to the upper/lower limit band and off-resonance from the center of the band, there can be a small contribution from the other bound state. In this case there is a small interference term $\sim r_{+}r_{-}e^{-(\kappa_{+}+\kappa_{-})\absolutevalue{n}}\cos[(y_{+}-y_{-})t]$. In fig.3 we show the dynamics of propagation of the radiation field for two relevant cases. The first situation is when the frequency of the emitter lies inside the band at $\Delta=3\xi/2$. This corresponds to a field that propagates out of the origin at a speed given by the group velocity evaluated at the corresponding value of the momentum $k$. The second case illustrates the dynamical formation of the upper bound state that localizes exponentially around the emitter as we can observe from fig.3(b).

V A single giant emitter in a cavity-array waveguide

Let us now consider a minimal model for a giant emitter. This consists of a two-level system coupled to $N_{c}$ cavities at positions $x_{j}$, as we illustrate in fig.4 for $N_{c}=2$. The Hamiltonian for this model is a straightforward generalization of the Hamiltonian (1), and it is given by,

[TABLE]

Transforming to the Fourier space and assuming equally symmetric spaced couplings with respect to the position of the emitter we get,

[TABLE]

where the couplings are now momentum-dependent, i.e.,

[TABLE]

On can compute the effective spectral function associated to the problem, which reads,

[TABLE]

where $J(\omega)$ is the spectral function corresponding to the emitter coupled to a single cavity in the waveguide (eq.(6)), and $\mathcal{G}(\omega)$ is a function of the frequency and the distance between coupling cavities modifying the spectral function as a consequence of considering a giant emitter [28]. In the following we focus on the particular case of $N_{c}=2$ as a minimal model. This case shows non-trivial dynamics and gives rise to the appearance of bound states inside the continuum (BIC).

V.1 Exact solution for $N_{c}=2$

In order to explore concrete examples and to obtain explicit expressions for the dynamics of the model with a giant emitter we will focus on the case of $N_{c}=2$. This is perhaps the simplest model of a giant emitter having only two equidistantly spaced coupling points to the waveguide. The effective spectral function can be written as,

[TABLE]

Interestingly, this spectral function has $\lceil d/2\rceil$ zeros at specific frequencies given by $\omega_{m}=2\xi\cos(\pi(2m+1)/d),$ with $m\in\mathbb{Z}$. This is shown in fig.5(b) where we show the effective spectral function as a function of frequency for different separation distances $d$. Notice that the function $G(s,d)$ is distance-dependent in this case, and can be obtained in a closed form as,

[TABLE]

where $\varrho(s)=is/2\xi-i\operatorname{sign}(\real\{s\})\sqrt{1+s^{2}/4\xi^{2}}$. The discontinuous jump across the branch cut is given by the difference,

[TABLE]

We can then write the solution as

[TABLE]

being the kernel function

[TABLE]

with $F(\tilde{y},d)=(g_{0}/\xi)^{4}(1+\cos(k(-\tilde{y})d))-4(g_{0}/\xi)^{2}(2\tilde{y}+\Delta/\xi)\sqrt{1-\tilde{y}^{2}}\sin(k(-y)d).$ Let us first discuss the contribution coming from scattering states to the dynamics. In the limit of $d\to 0$ the above kernel simplifies to the kernel for a small emitter given in eq.(15). It is interesting to analyze the behavior of the kernel function as it encodes information about the dynamics. For instance, the sharpest peak in this kernel will dominate the long-time dynamics. This is the case for the small emitter, where the narrow peaks close to the band-edges give rise to a power law decay $\sim t^{-1}$ of the occupation probability at intermediate times, that eventually transitions to a faster algebraic decay $\sim t^{-3}$ at very long times [11]. Let us first describe the situation of $\Delta$ near to the center of the band. In fig.5(a) we plot the kernel function given in eq.(27) along the frequency band. For $d\to 0$ (small emitter) we see that there is a main peak centered at the middle of the band and two singularities at the band-edges located at $\tilde{y}=\pm 1$. Near to these singularities we find two local maxima very close to $\pm 1$. As the emitter gets giant the kernel shows an interesting behavior. For example, for $d=2(2n+1)$ there is no central peak ($F(0,d)=0$) and the kernel behaves smoothly around the center of the band. The kernel function has $\lceil d/2\rceil-1$ zeros along the band located at same frequencies that the ones for the effective spectral function, except for the one at $\tilde{y}=0$, since the kernel function has a finite value at the center of the band. For large distances $d$, more zeros in the kernel function start to appear, making the kernel a highly oscillatory function. Contrary to the case of a small emitter, the initial decay is no longer well approximated by an exponential behavior and the contribution from scattering states gets enhanced as the distance is increased. This behavior is expected due to the fact that there is no Lorentzian-like central peak. For distances such that $d=4n$ the central peak for the small emitter case survives but becomes rounded with faster decaying tails for large $n$. In this case we obtain an approximately initial exponential decay with a rate of $\sim\gamma/2$. Interestingly, after a time $t=d/2\xi$ this decay gets enhanced or accelerated as consequence of the constructive interference of the two initial emitted radiation wavepackets from each distant cavity. The behavior of the decay probability at long times remains similar to the small emitter situation showing a power law behavior $\sim t^{-3}$.

V.2 Bound states inside the continuum

So far we have discussed the dynamics in terms of scattering and bound states that can exist outside the continuum of the allowed modes. It turns out that bound states can also be found inside the continuum and we refer to them as BIC states [60]. In order to search for the existence of such states in our minimal model we need to look at the zeros of the effective spectral function. From our earlier discussion we already encountered that $J_{\rm{eff}}(\omega_{m})=0$ for $\omega_{m}=2\xi\cos(\pi(2m+1)/d).$ If the frequency of emitter is such that $\Delta=\omega_{m}$ we have that the equation for the poles is $y-\omega_{m}+iG(-iy)=0$, which can be satisfied for $y=\omega_{m}$ since $G(-i\omega_{m})=0.$ Therefore, if the frequency of the emitter coincides with a zero of the effective spectral function there will be a non-zero contribution from a pole of the emitter probability amplitude $\alpha(s)$ in the Laplace domain [60]. We can then describe the dynamics in this situation (provided that there is no contribution from the bound states outside the continuum) as following,

[TABLE]

where r_{m}=\left[1+G^{\prime}(s,d)\big{|}_{s=-i\omega_{m}}\right]^{-1}, is the residue corresponding to the pole at $s=-i\omega_{m}$. An exact expression for this residue for arbitrary $d$ is shown in the appendix A. The situation of a giant emitter coupled simultaneously to two cavities is probably a minimal model for the existence of a bound state inside the continuum. A BIC state arises for $d=2$ when the frequency of the emitter is tuned exactly at the center of the band, as shown in figures 6(a) and 6(c). We see that there is fast initial decay given by the contribution from the scattering states followed by the contribution from the residue of the corresponding BIC, being this the dominant term resulting in a nonzero steady state probability at long times. As we increase $d$ more BIC states are allow to exists inside the continuum as we can notice from fig.6(b) for $d=12$. The corresponding field profile for a BIC state can be obtained exactly in the long time limit and it is given by,

[TABLE]

where $\phi=\arcsin(\Delta/2\xi)+\pi/2$. Figure 7(a) shows the formation of a BIC state and its profile at long times. Notice that the field is confined over the effective length of the giant emitter, i.e. between the two distant cavities coupling to the emitter, and is zero outside this region due to the destructive interference. This confinement of radiation in which the giant emitter behaves an an effective cavity has also been found in the case of two distant emitters as a consequence of the formation of BIC states [61].

V.3 Coherent oscillations of bound states

Interestingly, the giant emitter model also allows the coherent superposition of two bound states of different type, inside and outside the resonant band. This persistent coherent oscillation has also been found for BIC or dark states in the continuum of Ohmic waveguides for a giant atom having at least $N_{c}=3$ simultaneous equidistant interactions [58]. Here, however, we show that it is possible to have a coherent superposition of two distinct bound states lying outside and inside the continuum. Such sates take place whenever the frequency of the emitter satisfies the condition $J_{\rm{eff}}(\omega_{m})=0$, while simultaneously lying sufficiently close to the band edge in order to allow the contribution from a pole outside the band. Under such conditions we have that the long time limit of the occupation probability is given by

[TABLE]

where we have chosen the emitter frequency to match the frequency of the BIC state at $\Delta=\omega_{m}$ very close the upper band limit. This dynamics at long times is illustrated in fig.6(c) for $d=12$, where a non-decaying oscillating bound state with large amplitude is shown. We can see that eq.(30) agrees perfectly with the oscillation at long times for $\Delta=\omega_{0}=2\xi\cos(\pi/12).$ Another interesting situation is shown to happen for frequencies of the emitter non exactly matching the BIC condition and also lying sufficiently close of the band edge. In this case it is also possible to generate coherent Rabi-like oscillations lasting for very long times. These states will eventually decay since they are not truly bound states but quasi-bound states. We show one example of such state in fig.6(c) and its corresponding field in fig.7(c) for $d=30$ and $\Delta=2\xi$. In this case the integral over scattering states decays very slowly, having a quasi-steady state value for a long period of time before decaying algebraically at very long times. Strictly speaking, these states are coherent superpositions of slow decaying scattering modes and bound states outside the continuum. These oscillations will die at long times but it can be useful to address and have access to such modes in some particular experimental situations. The field corresponding to the oscillating bound state formed by the superposition of both, the BIC and BOC state is given by,

[TABLE]

where the field profiles for the two bound states are given by $\Phi_{+}(n)=e^{-(\kappa_{+}-i\pi/2)|n-d/2|}+e^{-(\kappa_{+}-i\pi/2)|n+d/2|}$ and $\Phi(n)=e^{i\phi|n-d/2|}+e^{i\phi|n+d/2|}$. This is shown in fig.7(b) for the evolution of the radiation field creating a breathing radiation mode characterizing an oscillating bound state as a combination of two bound states, one inside and one outside the continuum. The quasi-oscillating bound state of the field is shown in fig.7(c). The possibility to create or engineer such infinite oscillating states can be potentially exploited for the design and tailoring of non-stationary cavities supporting particular modes, either for trapping or manipulation, by means of the natural drive field of an oscillating mode, the quantum state of additional qubits interacting with the waveguide.

V.4 Master equation and entanglement entropy

We can generalize above results in order to allow the giant emitter to be initialized in any arbitrary initial state having coherences. The giant emitter density matrix at any time can be written as,

[TABLE]

from which we can obtain the following master equation,

[TABLE]

where $\tilde{\Delta}(t)=-2\Im{\dot{\alpha}(t)/\alpha(t)}$, and $\Gamma(t)=-2\Re{\dot{\alpha}(t)/\alpha(t)}$. The above master equation is time-local but in general non-Markovian, since $\Gamma(t)$ can take negative values [62]. An interesting quantity to analyze is the entanglement entropy resulting from computing the von Neumann entropy of the reduced density matrix (32), i.e., $S(\varrho)=-\Tr{\varrho\log\varrho}$. In the spontaneous emission process, the standard theory of Weisskopf and Wigner predicts that an initially excited atom undergoes an irreversible decay to the ground state. This means that an initial pure product state is mapped onto another pure final state. The entanglement entropy at the steady state is then zero and no atom-field entanglement remains at long times. As we have discussed in the previous sections, the presence of bound states modifies the steady state of the emitter-waveguide system and therefore some amount of entanglement generated in the process survives in the long time limit. In particular, if the steady state is an oscillating bound state like the one in eq.(30), the entanglement entropy is a time-dependent oscillating function with a frequency given by $y_{+}-\omega_{m}$.

VI Superconducting circuits model

Circuit-QED is the current suitable platform for the implementation and simulation of giant atoms since the first demonstration using superconducting circuits [56]. Here we provide the microscopic circuit model derivation of the full Hamiltonian model in (V).

VI.1 Transmission line model for a cavity-array waveguide

In order to simulate a cavity-array waveguide we consider a chain of capacitively coupled LC resonators. We describe such a system trough the classical Lagrangian,

[TABLE]

We assume that the chain is grounded at the boundary nodes with $\Phi_{0}=\Phi_{N+1}=0.$ The conjugated charges associated to each node are obtained via,

[TABLE]

The classical Lagrangian can also be written in matrix form defining the flux and charge vectors: $\mathbb{\Phi}=(\Phi_{1},\Phi_{2},...,\Phi_{N})^{T}$, $\mathbb{Q}=(Q_{1},Q_{2},...,Q_{N})^{T},$ as follows,

[TABLE]

being $\mathbb{C}$ the capacitance matrix with elements $C_{nm}=C_{\Sigma}\delta_{n,n}-C(\delta_{n,m+1}+\delta_{n+1,m}),$ and $\mathbb{L}$ a diagonal inductance matrix, respectively. Using this form the charge vector can be written as $\mathbb{Q}=\mathbb{C}\mathbb{\dot{\Phi}}$ and the Hamiltonian is obtained by performing the Legendre transform,

[TABLE]

The inverse capacitance matrix can be well approximated to first order assuming that $C\ll C_{0}$, resulting in the following approximate Hamiltonian,

[TABLE]

Classical canonical variables are now promoted to quantum hermitian flux and charge operators, satisfying the canonical commutation relation $[\Phi_{n},Q_{m}]=i\delta_{n,m}$. In terms of bosonic creation and annihilation operators they read,

[TABLE]

for which $[a_{n},a^{\dagger}_{m}]=\delta_{nm}.$ The second quantized Hamiltonian in the rotating wave approximation is then given by

[TABLE]

where $\omega_{0}=(L_{0}C_{\Sigma})^{-1/2}$, and $\xi=\omega_{0}C/2C_{\Sigma}$.

VI.2 Coupling to a giant transmon qubit

A giant transmon circuit is now coupled to the above transmission line at nodes positioned at $x_{j}$. This can be done by renormalizing the charge of the transmon Hamiltonian. The full Hamiltonian reads,

[TABLE]

where $C_{g}$ is the coupling capacitance and $C_{\Sigma^{\prime}}=C_{\Sigma}+C_{g}$ is the renormalized capacitance of the transmission line due to the coupling to the transmon. The second-quantized Hamiltonian of the giant transmon coupled to the transmission line up to first order in $C_{g}/C_{\Sigma^{\prime}}$ is given by

[TABLE]

where $\Delta\approx\sqrt{8E_{C}E_{J}}$ and $g_{0}=C_{g}\sqrt{\Delta\omega_{0}/C_{\Sigma_{q}}C_{\Sigma^{\prime}}}/2.$ Limiting our transmon to two lowest energy levels (transmon qubit) we get the Hamiltonian in eq.(V) up to a normalization factor of $N^{-1}_{c}$ in the coupling parameter. It is worth mentioning that in the presence of the transmon qubit, the frequencies $\omega_{0}$ and $\xi$ in $H_{\rm{photons}}$ should be obtained using the renormalized capacitance of the transmission line, $C_{\Sigma^{\prime}}$.

VII Conclusions

We have presented exact analytical solutions for the spontaneous emission dynamics of quantum emitters coupled to structured photonic environments, focusing in particular on giant emitters interacting with cavity–array waveguides. By employing the exact inversion of the Laplace transform, we have gone beyond the standard Weisskopf–Wigner approximation and demonstrated how the interplay of interference and non-Markovian effects leads to qualitatively new behavior. Our results reveal the existence of distinct bound states, both inside and outside the continuum, and show that giant emitters enable their coherent superposition, giving rise to persistent oscillatory dynamics. This provides a minimal and exactly solvable model for understanding interference-induced localization and long-lived quantum coherences in structured reservoirs. We have outlined a feasible implementation using superconducting circuit architectures, where the required cavity–array and giant atom configurations are within reach of current experiments. This suggests that the predicted phenomena, such as bound states in the continuum and their coherent interplay with conventional bound states, could be observed and exploited in near-term devices in circuit-QED. These results may be potentially extended to giant multi-emitter systems, non-Hermitian interactions, and time-dependent driving setups, offering a route toward engineered non-Markovian environments for quantum information processing and quantum simulation with giant emitters.

Acknowledgements.

This work was supported by the Secretaría de Ciencia, Humanidades, Tecnología e Innovación (SECIHTI) Mexico, under grant No. CBF2023-2024-2888, and by DGAPA-PAPIIT-UNAM under grant No. IA104625. The author gratefully acknowledges valuable discussions with Sergi Terradas, Luis M. Moreno, and David Zueco on related projects.

Appendix A Residue for the BIC state

The BIC state corresponds to a pole inside the brach cut in the interval $[-2\xi,2\xi]$. In this case the integration path in fig.2(a) can be deformed in order to account for this contribution. This results in the following residue for eq.(28),

[TABLE]

where,

[TABLE]

being $\eta(u)\equiv\sqrt{1-u^{2}}$, and $u^{-1}=\omega_{m}/2\xi$, $|u|<1$.