Spontaneous emission of Schrödinger cats in a waveguide at ultrastrong coupling

Photons describe the granular structure of the electromagnetic field radiated by single coherent sources, such as atoms and quantum dots. These quanta of light constitute well-defined monochromatic excitations because the spontaneous emission rate Γ is much smaller than the transition frequency Δ of an emitter. More precisely, for atomic decay in three-dimensional space [1], the width of the transition lines, $\alpha ={\rm{\Gamma }}/{\rm{\Delta }}$, also defines a dimensionless coupling constant $\alpha ={({ \mathcal P }/e\lambda )}^{2}{\alpha }_{\mathrm{QED}}$ involving two small factors: the ratio of the atomic dipole ${ \mathcal P }$ to the photon wavelength λ, and the vacuum fine structure constant ${\alpha }_{\mathrm{QED}}\simeq 1/137$. Typically ${ \mathcal P }/e\lambda \lesssim {\alpha }_{\mathrm{QED}}$, so that the light–matter coupling to a 3D continuum, $\alpha \lesssim {({\alpha }_{\mathrm{QED}})}^{3}\simeq {10}^{-6}$, is vanishingly small.

With the advent of microwave quantum optics in superconducting one-dimensional waveguides [2–9], the coupling constant can reach the still small value $\alpha \simeq {\alpha }_{\mathrm{QED}}\simeq {10}^{-2}$. Thus, a coupling of order one can only be attained by tweaking the fine structure constant itself. This requires the use of superconducting waveguides with high impedance, since ${\alpha }_{\mathrm{QED}}={Z}_{\mathrm{vac}.}/2{R}_{{K}}$ can be interpreted as the ratio of the vacuum impedance ${Z}_{\mathrm{vac}.}=1/{\epsilon }_{0}c\simeq 376\ {\rm{\Omega }}$ to the quantum of resistance ${R}_{{K}}={e}^{2}/h\simeq 25812\ {\rm{\Omega }}$. Long chains of Josephson junctions [10–13] constitute low-loss metamaterials for the propagation of microwave photons with characteristic impedance Z up to the order of R_K, allowing to obtain values of the coupling constant $\alpha =Z/{R}_{K}$ of order one. This physical regime reveals anomalous scattering properties of photons [14–18] impinging on a non-trivial dressed vacuum [19]. This Letter aims at answering a seemingly simple question: is the radiation from a single emitter still described by a discrete photon in the ultrastrong coupling regime?

Clearly, the radiation released for $\alpha \simeq 1$ is strongly non-monochromatic, because ultrafast energy relaxation processes gives a large broadening of the atomic transition lines. In addition to modified spectral properties, drastic effects also occur in the quantum statistics of the radiation. Indeed, in contrast to spontaneous decay at weak coupling, where the photon state emerges as an n = 1 Fock excitation [1], the radiated field at ultrastrong coupling contains states with photon number larger than one. This surprising result can be understood as follows: photons with energy much smaller or much larger than the atomic transition frequency can be excited because one is not limited to resonant transitions in the non-perturbative regime of quantum electrodynamics (QED). But for low energy photons, the large available energy and the strong coupling constant make favorable the creation of multiple photon excitations. Infact, we find that the full quantum state of the electromagnetic field possesses a simple structure given by an optical Schrödinger cat, namely a superposition of strongly displaced coherent states. We stress that cat states obtained by spontaneous decay in a large impedance environment are very different from the cats obtained in cavities [20], for two reasons. First, as already mentioned above, the radiation is spectrally broad at ultrastrong coupling, and thus the resulting cat states are rather localized in the time domain rather than in frequency. Secondly, we find that these cat states display an intrinsic loss of coherence at longer times than the sudden time scale T₁ associated to energy relaxation, but shorter than the long decoherence time T₂ for complete memory loss of the initial atomic state. The separation of time scales ${T}_{2}\gg {T}_{1}$ at ultrastrong coupling, see figure 1, illustrates again the stark difference with the weak-coupling regime of quantum optics, where ${T}_{2}=2{T}_{1}$, and will be one of the main focus of the rest of the paper.

Zoom In Zoom Out Reset image size

Having presented the physics at play, we must emphasize that the computational aspects of QED at ultrastrong coupling are far from trivial, due to the breakdown of perturbation theory [19, 21, 22]. In contrast to cavity-QED at ultrastrong coupling [23–25], the Hilbert space is unmanageable for high impedance waveguides because: (i) a large number of electromagnetic modes are involved, up to ${N}_{\mathrm{modes}}=1000$ for long chains of Josephson junctions; (ii) the average total number of photons ${\bar{n}}_{\mathrm{tot}}$ is larger than one. For the case $\alpha =1.0$ considered in the following, one obtains up to ${\bar{n}}_{\mathrm{tot}}=6$ photons, so that the full quantum mechanical problem requires to tackle more than ${({N}_{\mathrm{modes}})}^{{\bar{n}}_{\mathrm{tot}}}\simeq {10}^{18}$ quantum states, well beyond the reach of brute force diagonalization. Based on the physical idea that coherent states are the most stable quantum states at ultrastrong coupling [26–28], we propose here a new methodology for harnessing cat states which relies on an efficient computational algorithm. This technique is not only conceptually simpler than state of the art methods [29–33], but also very powerful.

We will specify our study to the standard model of waveguide-QED, described by the Hamiltonian:

$$\begin{eqnarray}&&H=\displaystyle \frac{{\rm{\Delta }}}{2}{\sigma }_{x}+\displaystyle \sum _{k}{\omega }_{k}{a}_{k}^{\dagger }{a}_{k}-{\sigma }_{z}\displaystyle \sum _{k}\displaystyle \frac{{g}_{k}}{2}({a}_{k}^{\dagger }+{a}_{k}),\end{eqnarray} \tag{ 1 }$$

defining the light–matter coupling g_k of a given mode with frequency ${\omega }_{k}$ to the two-level atom, described by Pauli matrices ${\sigma }_{j}$ and transition frequency Δ. At variance with quantum optics conventions, we have intentionally written the atomic splitting as a ${\sigma }_{x}$ term, in order to emphasize the natural selection of coherent states caused by the ${\sigma }_{z}$ light–matter coupling. Hence, the bare atomic ground state is $| g\rangle =[| \uparrow \rangle -| \downarrow \rangle ]/\sqrt{2}$ in our notation. The dimensionless coupling strength α can be encapsulated from the spectral density $J(\omega )=\pi {\sum }_{k}{g}_{k}^{2}\delta (\omega -{\omega }_{k})=2\pi \alpha \omega {{\rm{e}}}^{-\omega /{\omega }_{{\rm{p}}}}$, with ${\omega }_{{\rm{p}}}$ the plasma frequency. We assume here a linear dispersion relation ${\omega }_{k}=k$, with unit speed of light in the medium, which is justified for atomic transitions well below ${\omega }_{{\rm{p}}}$. In addition, we do not consider here the important effect of coupling the high impedance finite-size Josephson waveguide to 50 Ω lines, which will restrict the spectrum to a set of discrete resonances (slightly broadened by the contacts). However, for a long-enough Josephson array, the limits in resolution will be well-below the characteristic scales of the system. Finally, we consider purely unitary dynamics, assuming that extrinsic losses, apart from the coupling to external contacts, are negligible.

As discussed above, the state vector is represented at all times by superposed coherent states [26–28]:

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =\displaystyle \sum _{m=1}^{{N}_{\mathrm{cs}}}[{p}_{m}(t)| {f}_{m}(t)\rangle | \uparrow \rangle +{q}_{m}(t)| {h}_{m}(t)\rangle | \downarrow \rangle ].\end{eqnarray} \tag{ 2 }$$

Here the set of amplitudes $\{{p}_{m},{q}_{m}\}$ are complex and time dependent, and a set of discrete multimode coherent states are introduced: $| {f}_{m}(t)\rangle ={{\rm{e}}}^{{\sum }_{k=1}^{{N}_{\mathrm{modes}}}[{f}_{m}(k,t){a}_{k}^{\dagger }-{f}_{m}^{* }(k,t){a}_{k}]}| 0\rangle$. This decomposition allows in principle to target an arbitrary state of the full Hilbert space (with an exponential cost). One must note a strong difference here with the standard Glauber–Sudarshan decomposition, which relies on a continuous expansion of the state vector onto coherent states. Because of the large number of modes involved in waveguide QED at ultrastrong coupling, the continuous representation is not suitable for numerical purposes. Our discrete expansion (2) can however be understood as a discretized version of the continuous integral representation of an arbitrary wavefunction onto coherent states. Most importantly, we argue in this paper that physical states obtained from standard protocols (such as spontaneous emission) can be efficiently simulated with a discrete set of coherent states, showing only a polynomial cost in the number of modes N_modes and coherent states N_cs. This computational gain was previously demonstrated for the full ground state of (1), and is here extended to the dynamics, allowing to simulate for the first time large systems, up to thousands of electromagnetic modes.

Technically, the exact Schrödinger dynamics controlled by the Hamiltonian (1) can be obtained from the real Lagrangian density:

$$\begin{eqnarray}&&{ \mathcal L }=\langle {\rm{\Psi }}(t)| \displaystyle \frac{{\rm{i}}}{2}\overrightarrow{{\partial }_{t}}-\displaystyle \frac{{\rm{i}}}{2}\overleftarrow{{\partial }_{t}}-{ \mathcal H }| {\rm{\Psi }}(t)\rangle ,\end{eqnarray} \tag{ 3 }$$

by applying the time-dependent variational principle [34], $\delta \int {\rm{d}}t{ \mathcal L }=0$, upon arbitrary variations of the state vector (2). This results in Euler–Lagrange equations $\tfrac{{\rm{d}}}{{\rm{d}}{t}}\tfrac{\partial { \mathcal L }}{\partial \dot{v}}=\tfrac{\partial { \mathcal L }}{\partial v}$ for the set of variables $v=\{{p}_{m},{q}_{m},{f}_{m}(k),{h}_{m}(k)\}$, which can be solved by numerical integration [35–38] using a specially devised algorithm (see appendix for details).

We provide here clear proof of the good convergence of our algorithm. First, the left panel of figure 2 shows how the time-dependent error $\mathrm{Err}({t})$ vanishes with the number of coherent states, for the deep ultracoupling regime $\alpha =1.0$. Here, the error is defined by the squared norm $\mathrm{Err}(t)\equiv \langle {\rm{\Phi }}(t)| {\rm{\Phi }}(t)\rangle$ of the auxiliary state $| {\rm{\Phi }}(t)\rangle \equiv ({\rm{i}}{\partial }_{t}-H)| {\rm{\Psi }}(t)\rangle$, which is zero for the exact Schrödinger dynamics. Clearly the error is already small at all times for a single coherent state ${N}_{\mathrm{cs}}=1$, and goes quickly to zero with increasing number of terms in the decomposition equation (2). The precise scaling of the algorithm with N_cs is demonstrated in the right panel of figure 2. We find that the error (here computed at a finite and fixed time T = 80) decreases typically with an inverse square power ${[{N}_{\mathrm{cs}}]}^{-2}$, independently of the coupling strength α. This shows that our methodology is based on a physically well-motivated decomposition of the state vector, and is not tied to a particular regime of the spin-boson model. The coherent state expansion of the time-dependent state vector (2) thus provides numerically accurate results in all regimes of coupling for a small computational effort.

Zoom In Zoom Out Reset image size

We consider first the important issue of the time scales governing the physics in the ultrastrong coupling regime. For this purpose, we prepare the initial state as $| {\rm{\Psi }}(t=0)\rangle =| 0\rangle \otimes | \uparrow \rangle =| 0\rangle \otimes [| g\rangle +| e\rangle ]/\sqrt{2}$ in a superposition of the two bare atomic levels, and with the environment in its vacuum. Since the two-level splitting is described by a ${\sigma }_{x}$ coupling in equation (1), one expects precession and decay to zero of the transverse spin component $\langle {\sigma }_{z}(t)\rangle$ on a time scale T₂, and relaxation of the longitudinal term $\langle {\sigma }_{x}(t)\rangle$ towards its finite equilibrium value on a time scale T₁. This standard behavior is well obeyed in the weak-coupling regime $\alpha \ll 1$, with ${T}_{2}=2{T}_{1}$, as seen from the dot-dashed green curve in the first two panels of figure 1. For the intermediate value $\alpha =0.4$ (full blue curve), precession of $\langle {\sigma }_{z}(t)\rangle$ is nearly overdamped, as is well established [21].

Remarkably, the ultrastrong coupling regime $\alpha =1.0$ (dashed red curve) shows a striking decoupling between the decoherence time T₂ and the energy relaxation time T₁, with ${T}_{2}\gg {T}_{1}$. The underlying physics can be anticipated: energy relaxation, related to emission of radiation, occurs on the short time scale T₁, because the excited atomic level is strongly damped at large α. This is clearly seen from the rapid saturation of $\langle {\sigma }_{x}(t)\rangle$ in the middle panel of figure 1 at $\alpha =1.0$. However, rapid decoherence is prohibited because the atom is dressed by its electromagnetic environment [19] on a large spatial scale ${L}_{K}={({\omega }_{{\rm{p}}}/{\rm{\Delta }})}^{\alpha /({\alpha }_{c}-\alpha )}$ (here ${\alpha }_{c}=1+{\rm{\Delta }}/{\omega }_{{\rm{p}}}$, corresponding to the threshold of full localization). Decoherence thus takes a considerable time ${T}_{2}\simeq {L}_{K}$ for the complete relaxation of all quantum correlations between the dressed atom and the is larger than one,radiated field. This effect is seen by the very slow decay of $\langle {\sigma }_{z}(t)\rangle$ (dashed red curve in the left panel of figure 1). Indeed, for $\alpha =1.0$ and ${\rm{\Delta }}/{\omega }_{{\rm{p}}}=0.2$, ${T}_{2}\simeq 3000$, which was not reached at the final time of our simulations. As we will see in the following, this separation of time scales is key to understanding the physics at play, due to its strong impact on the structure of the emitted light.

We turn to the emission protocol, taking now the bare excited level of the atom $| {\rm{\Psi }}(t=0)\rangle =| 0\rangle \otimes [| \uparrow \rangle +| \downarrow \rangle ]/\sqrt{2}=| 0\rangle \otimes | e\rangle$ as the initial state. According to Wigner–Weisskopf theory [1] valid at weak coupling $\alpha \ll 1$, one expects spontaneous emission of a single photon and decay of the atom towards its bare ground state. Although our theory (2) is based on coherent states and not Fock excitations, we can show that it does recover the standard quantum optics results at $\alpha \to 0$, while providing a simple physical picture in the ultrastrong coupling regime. For this purpose, we plot in figure 3 the real part of the coherent state amplitudes ${f}_{m}(x,T)$ and ${h}_{m}(x,T)$ (for m = 1, 2) for ${N}_{\mathrm{cs}}=2$, as a function of the spatial separation x along the waveguide (the atom is located at x = 0). Here we choose a time $T\gg {T}_{2}\simeq {L}_{K}$ long enough that the system has fully relaxed to its dressed ground state. The emitted wavepacket is clearly seen at a distance $x\simeq T$, due to ballistic propagation of the wavefront, and the displacements satisfy the following relations: ${f}_{1}(x,t)\simeq -{f}_{2}(x,t)\simeq -{h}_{1}(x,t)\simeq {h}_{2}(x,t)\equiv {f}_{\mathrm{wp}}(x-t)$ (in addition, ${p}_{1}=-{p}_{2}={q}_{1}=-{q}_{2}\equiv p$ (not shown)). Within the spatial region $x\lt {L}_{K}$ where the atom is dressed by its static surrounding cloud, a different set of relations is observed: ${f}_{1}(x,t)\simeq {f}_{2}(x,t)\simeq -{h}_{1}(x,t)\simeq -{h}_{2}(x,t)\equiv {f}_{\mathrm{cl}}(x)$. Factoring each coherent state in real space $| {f}_{1}\rangle =| {f}_{\mathrm{cl}}\rangle \otimes | {f}_{\mathrm{wp}}\rangle$ simplifies the wavefunction (2):

$$\begin{eqnarray}| {\rm{\Psi }}(T)\rangle & \simeq & [{p}_{1}| {f}_{1}\rangle +{p}_{2}| {f}_{2}\rangle ]| \uparrow \rangle +[{q}_{1}| {h}_{1}\rangle +{q}_{2}| {h}_{2}\rangle ]| \downarrow \rangle \\ & \simeq & [p| {f}_{\mathrm{cl}}\rangle | {f}_{\mathrm{wp}}\rangle -p| {f}_{\mathrm{cl}}\rangle | -{f}_{\mathrm{wp}}\rangle ]| \uparrow \rangle +[p| -{f}_{\mathrm{cl}}\rangle | -{f}_{\mathrm{wp}}\rangle -p| -{f}_{\mathrm{cl}}\rangle | {f}_{\mathrm{wp}}\rangle ]| \downarrow \rangle \\ & \simeq & [| {f}_{\mathrm{cl}}\rangle | \uparrow \rangle -| -{f}_{\mathrm{cl}}\rangle | \downarrow \rangle ]\otimes p[| {f}_{\mathrm{wp}}\rangle -| -{f}_{\mathrm{wp}}\rangle ].\end{eqnarray} \tag{ 4 }$$

The last line in (4) is straightforwardly interpreted as the absence of correlations between the dressed qubit $| {{\rm{\Psi }}}_{\mathrm{cl}}\rangle \equiv [| {f}_{\mathrm{cl}}\rangle | \uparrow \rangle -| -{f}_{\mathrm{cl}}\rangle | \downarrow \rangle ]/\sqrt{2}$ (which comprises both the atom and its entangled neighboring cloud in the waveguide) and the emitted wavepacket $| {{\rm{\Psi }}}_{\mathrm{wp}}\rangle \equiv \sqrt{2}p[| {f}_{\mathrm{wp}}\rangle -| -{f}_{\mathrm{wp}}\rangle ]$. This is physically expected as the system relaxes at long times to a unique dressed ground state, while emitting a stream of electromagnetic radiation carrying the excess energy but no quantum correlations with the atom. We stress that the above expressions for $| {{\rm{\Psi }}}_{\mathrm{cl}}\rangle$ and $| {{\rm{\Psi }}}_{\mathrm{wp}}\rangle$ are only approximate, as small quantum corrections arise at increasing α [26, 27], which are accounted for by the terms $n\gt 2$ in the expansion (2).

Zoom In Zoom Out Reset image size

The approximate wavefunction (4) nicely recovers the result of Wigner–Weisskopf theory in the quantum optics regime $\alpha \to 0$. Indeed, in this case the displacements in the dressing cloud are vanishingly small, $| {f}_{\mathrm{cl}}(x)| \ll 1$, and thus $| {{\rm{\Psi }}}_{\mathrm{cl}}\rangle =| 0\rangle [| \uparrow \rangle -| \downarrow \rangle ]/\sqrt{2}\equiv | 0\rangle | g\rangle$, so that the atom has correctly relaxed to its bare ground state. The quantum state describing the emitted light also simplifies, since the displacements in the wavepacket are small, $| {f}_{\mathrm{wp}}| \ll 1$. A first-order Taylor expansion of the coherent state gives $| {{\rm{\Psi }}}_{{\rm{wp}}}\rangle =2\sqrt{2}p{\sum }_{k}{f}_{{\rm{wp}}}(k){{\rm{e}}}^{{\rm{i}}{k}{t}}{a}_{k}^{\dagger }| 0\rangle$, which is the expected one-photon Fock state. The monochromatic nature of the emitted light can be seen from the spatiotemporal oscillations of the fields in the upper panel of figure 3, with an envelope controlled by the underdamped dynamics of the qubit (first panel of figure 1 for $\alpha =0.1$). Moving towards the ultrastrong coupling regime for increasing α values (middle and lower panel in figure 3), two major changes occur. First, the displacements in the entanglement cloud $x\lt {L}_{K}$ penetrate deeper and deeper within the waveguide, due to the increase of the screening length L_K with α, as shown by the shaded area in figure 3. Second, the emitted wavepacket becomes very localized temporally, with associated displacements that clearly grow in magnitude. This indicates that the radiation is spectrally broad, and that the number of emitted photons grows with increasing α, as was anticipated.

Let us finally analyze in more detail the quantum properties of the emitted radiation, starting with the approximate expression $| {{\rm{\Psi }}}_{{\rm{wp}}}\rangle =\sqrt{2}p[| {f}_{{\rm{wp}}}\rangle -| -{f}_{{\rm{wp}}}\rangle ]$ for the emitted wavepacket. Although this state is a Fock state when the average number of emitted photons is close to one, as shown above, it clearly turns into an odd parity Schrödinger cat when the displacements grow large, as seen for $\alpha =0.8$ in the lower panel in figure 3. In order to check this idea more precisely, we compute the Wigner distribution of the emitted wavepacket. In contrast to photons released into a cavity [20, 39], we stress again that the radiation is not purely monochromatic, due to the significant damping effect on the atom caused by the strong coupling to the waveguide. Accordingly, following standard practice [40], one uses an optimized temporal filter w(t) of the output signal. The filter function w(t) is defined as to match precisely the shape of any of the wavefronts in figure 3 (which show the same amplitude), in the spatial domain where the emitted signal is decoupled from the short distance cloud (here in the approximate range $100\lt x\lt 250$). For the values of $\alpha \leqslant 0.8$ in figure 3, a measurement time T = 200 is sufficient to ensure this decoupling, but for the computation of figure 4 with $\alpha =1.0$, a longer time T = 900 was required, due to the increase in T₂.

Zoom In Zoom Out Reset image size

Owing to the linear dispersion in the waveguide, the averaging can be performed spatially by defining an effective creation operator ${b}^{\dagger }={\sum }_{x}w(x){a}^{\dagger }(x){\rm{\Theta }}(x-{L}_{K})$, in which a Θ-function is used to take out the static bound component of the screening cloud. The global scale of the filter function w(t) is set by imposing standard commutation relation for the effective mode, namely $[b,{b}^{\dagger }]=1$. The Wigner distribution [20] is then computed using the standard expression $W(\lambda )=\int ({{\rm{d}}}^{2}\beta /{\pi }^{2}){C}_{{\rm{s}}}(\beta ){{\rm{e}}}^{\lambda {\beta }^{* }-{\lambda }^{* }\beta }$, with the symmetrized correlation function ${C}_{{\rm{s}}}(\beta )=\langle {\rm{\Psi }}| {{\rm{e}}}^{\beta {b}^{\dagger }-{\beta }^{* }b}| {\rm{\Psi }}\rangle$. For $\alpha =0.1$, the phase space distribution in the left panel of figure 4 shows the characteristic circular form of the n = 1 Fock state with negative quasi-probability $-2/\pi$ at the origin.

At ultrastrong coupling for $\alpha =1.0$, the emitted radiation undergoes radical changes, as shown in the right panel of figure 4. The Wigner distribution now presents two positive lobes (signature of the two classical configurations of the cat), but also a negative region near the origin, fingerprint of the characteristic quantum interference, or 'whiskers', of a Schrödinger cat. Surprisingly, the maximum negative amplitude does not reach the expected value $-2/\pi$ of a perfectly coherent cat, although our system does not present any extrinsic source of decoherence for the optical modes (such as leaks into a 3D continuum). We argue that the physical source of decoherence is the dressed qubit itself, a very unusual feature. This phenomenon can be understood from the qubit dynamics shown in figure 1, in relation to the separation of time scales ${T}_{2}\gg {T}_{1}$. Indeed, the cat state is emitted on a short scale ${T}_{1}\simeq 1/{\omega }_{{\rm{p}}}$ in the ultrastrong coupling regime, due to the sudden release of energy. But the atom maintains its coherence on a longer time scale ${T}_{2}\simeq {L}_{K}$ due to its long-range spatial entanglement with the waveguide, since $\langle {\sigma }_{z}(t)\rangle$ does not decay. From the no-cloning theorem [1], the quantum information stored in the dressed atom state cannot be transfered to the wavepacket for times t ≪ T₂, and thus the coherence of the emitted cat state is only partial at intermediate timescales. This explains why the negative lobe of the cat state in the right panel of figure 4 does not quite reach the maximal value $-2/\pi$.

We have to point out that a definite phase is seen in the cat state of figure 4, which evolves as the state propagates between the emitter and the measurement setup. Its absolute value is tied to the original form of the Hamiltonian (1), where the qubit is seen to couple to the first quadrature of the field. Moreover, it must be stressed that a true measurement setup will take place outside the waveguide, which must be adapted to a low impedance environment. It may be that phase information is ultimely lost in the final output field, as pointed out previously in the context of anisotropic dielectrics [41], hence perhaps modifying the cat structure. Interestingly, this problem bears some similarity to transport in interacting quantum wires (Luttinger liquids [42]), which can be described by squeezed plasmonic modes. The interaction fingerprints visible in the conductance of an infinite wire are indeed suppressed once the wire is smoothly connected to non-interacting leads [43, 44]. The question of impedance matching in ultra-strongly coupled waveguides has not been addressed to our knowledge in the recent circuit-QED literature [14–18].

We conclude our analysis of the emitted cat by displaying its actual photon content. In the left panel of figure 5, we show the spectrally resolved photon number density, which is obtained as previously by cutting out the static part of the field tied to the atom. In the quantum optics regime $\alpha =0.01$, the lineshape is a narrow Lorentzian centered at the bare transition frequency Δ. For increasing α, the peak is shifted to lower frequencies due to the dressing of the atom by the bosonic bath, but more remarkably, its lineshape becomes spectrally very broad. This illustrates again the fact that the cat state generated by spontaneous emission has a rather localized character in the time domain. We stress that the number of emitted photons has to vanish for zero frequency modes, because the coupling constant g_k vanishes at $k\to 0$. This is also true for the curve $\alpha =0.9$, although the downturn of the curve can barely be seen on the graph, due to the strong renormalization of the emission frequency to tiny frequencies. The size of the cat can be measured by integrating the photon density curves, giving the photon numbers shown in the right panel of figure 5, both for the wavepacket (solid line) and the total wavefunction (dashed line) comprising wavepacket and entanglement cloud. At ultrastrong coupling for $\alpha =1.0$, the wavepacket contains ${\bar{n}}_{\mathrm{wp}}\simeq 3$ photons, in agreement with the cat size in figure 4, with an average total photon number of ${\bar{n}}_{\mathrm{tot}}={\bar{n}}_{\mathrm{wp}}+{\bar{n}}_{\mathrm{cl}}=6$, which highlights the complexity of the quantum problem at play, involving more than ${1000}^{{\bar{n}}_{\mathrm{tot}}}\simeq {10}^{18}$ quantum states.

Zoom In Zoom Out Reset image size

To summarize, we have found that Schrödinger cat states are spontaneously radiated by a single emitter in an infinite high impedance medium at ultrastrong coupling. These cats show unusual properties in comparison to standard quantum optics protocols [20, 45]: they are spectrally very broad and partially quantum coherent at intermediate stages of the dynamics, due to a strong separation in time scales between the slow decoherence and the fast energy relaxation of the emitter. After complete information loss of an arbitrary initial state of the qubit $[u| g\rangle +v| e\rangle ]\otimes | 0\rangle$, the qubit relaxes to a unique many-body ground state $| \mathrm{GS}\rangle$, associated to the polarization in the waveguide. Quantum information is then preserved in a generic cat state $| \mathrm{GS}\rangle \otimes [(u+v)| +f\rangle +(u-v)| -f\rangle ]$ of the radiated field at long times. However, impedance matching of the waveguide to a low-impedance external measurement apparatus may drastically affect the quantum correlations of spontaneous emission (as well as the structure of inelastically scattered light [15, 17]), a generic problem that will be addressed in a future study.

We thank H Baranger, A Chin, T Grall, T Meunier, A Nazir, J Puertas-Martinez, M Schiró and I Snyman for useful discussions, and especially N Roch for stimulating several aspects of this work. Financial support from the Nanoscience Foundation is also recognized.

The starting point of the method is to represent the state vector at an arbitrary time by an expansion onto coherent states:

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =\displaystyle \sum _{n=1}^{{N}_{\mathrm{cs}}}[{p}_{n}(t)| {f}_{n}(t)\rangle | \uparrow \rangle +{q}_{n}(t)| {h}_{n}(t)\rangle | \downarrow \rangle ],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&| {f}_{m}(t)\rangle ={{\rm{e}}}^{{\displaystyle \sum }_{k=1}^{{N}_{\mathrm{modes}}}[{f}_{m}(k,t){a}_{k}^{\dagger }-{f}_{m}^{* }(k,t){a}_{k}]}| 0\rangle ,\end{eqnarray} \tag{ 6 }$$

where the set of variables $v=\{{p}_{n},{q}_{n},{f}_{n}(k),{h}_{n}(k)\}$ are complex and time dependent, with $n=1\ldots {N}_{\mathrm{cs}}$. Following the time-dependent variational principle [34], we define the real Lagrangian density ${ \mathcal L }=\langle {\rm{\Psi }}(t)| \tfrac{{\rm{i}}}{2}\overrightarrow{{\partial }_{t}}-\tfrac{{\rm{i}}}{2}\overleftarrow{{\partial }_{t}}-{ \mathcal H }| {\rm{\Psi }}(t)\rangle$, from which arbitrary variations of the state vector give the simple Euler–Lagrange equations $\tfrac{{\rm{d}}}{{\rm{d}}{t}}\tfrac{\partial { \mathcal L }}{\partial \dot{v}}=\tfrac{\partial { \mathcal L }}{\partial v}$ for our set of variables. For the spin-boson model [21], the equations explicitly read:

$$\begin{eqnarray}&&-{\rm{i}}\displaystyle \frac{\partial E}{\partial {p}_{j}^{* }}=\displaystyle \frac{1}{2}\displaystyle \sum _{m}(2{\dot{p}}_{m}-{p}_{m}{\kappa }_{{mj}})\langle {f}_{j}| {f}_{m}\rangle ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}-{\rm{i}}\displaystyle \frac{\partial E}{\partial {f}_{j}^{k* }} & = & \displaystyle \sum _{m}\left[{p}_{m}{p}_{j}^{* }{\dot{f}}_{m}^{k}\langle {f}_{j}| {f}_{m}\rangle -\displaystyle \frac{1}{4}(2{\dot{p}}_{m}-{p}_{m}{\kappa }_{{mj}}){p}_{j}^{* }({f}_{j}^{k}-2{f}_{m}^{k})\langle {f}_{j}| {f}_{m}\rangle \right.\\ & & \left.+\displaystyle \frac{1}{4}(2{\dot{p}}_{m}^{* }-{p}_{m}^{* }{\kappa }_{{mj}}^{* }){p}_{j}{f}_{j}^{k}\langle {f}_{m}| {f}_{j}\rangle \right],\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\kappa }_{{mj}}=\displaystyle \sum _{k^{\prime} \gt 0}[{\dot{f}}_{m}^{k^{\prime} }{f}_{m}^{k^{\prime} * }\,+\,{\dot{f}}_{m}^{k^{\prime} * }{f}_{m}^{k^{\prime} }-2{f}_{j}^{k^{\prime} * }{\dot{f}}_{m}^{k^{\prime} }].\end{eqnarray} \tag{ 9 }$$

Identical equations (up to a minus sign) are obtained for the variables q_n and h^k_n. We have denoted here $E=\langle {\rm{\Psi }}| H| {\rm{\Psi }}\rangle$ the average energy, whose explicit expression is:

$$\begin{eqnarray}E & = & \displaystyle \frac{{\rm{\Delta }}}{2}\displaystyle \sum _{n,m}({p}_{n}^{* }{q}_{m}\langle {f}_{n}| {h}_{m}\rangle +{p}_{m}{q}_{n}^{* }\langle {h}_{n}| {f}_{m}\rangle )+\displaystyle \sum _{n,m}({p}_{n}^{* }{p}_{m}\langle {f}_{n}| {f}_{m}\rangle {W}_{{nm}}^{f}+{q}_{n}^{* }{q}_{m}\langle {h}_{n}| {h}_{m}\rangle {W}_{{nm}}^{h})\\ & & -\displaystyle \frac{1}{2}\displaystyle \sum _{n,m}({p}_{n}^{* }{p}_{m}\langle {f}_{n}| {f}_{m}\rangle {L}_{{nm}}^{f}-{q}_{n}^{* }{q}_{m}\langle {h}_{n}| {h}_{m}\rangle {L}_{{nm}}^{h}),\end{eqnarray} \tag{ 10 }$$

where we have defined ${W}_{{nm}}^{f}={\sum }_{k\gt 0}{\omega }_{k}{f}_{n}^{k* }{f}_{m}^{k}$, ${W}_{{nm}}^{h}={\sum }_{k\gt 0}{\omega }_{k}{h}_{n}^{k* }{h}_{m}^{k}$, ${L}_{{nm}}^{f}={\sum }_{k\gt 0}{g}_{k}({f}_{n}^{k* }\,+\,{f}_{m}^{k})$, ${L}_{{nm}}^{h}={\sum }_{k\gt 0}{g}_{k}({h}_{n}^{k* }+{h}_{m}^{k})$.

In contrast to the dynamics with a single coherent state [36, 38], one encounters here a computational difficulty [35], because time-derivatives ${\dot{f}}_{m}^{k^{\prime} }$ of all possible coherent state amplitudes enter the dynamical equation ruling a given field f_j^k in equation (8) through the parameter ${\kappa }_{{mj}}$ in equation (9). Indeed, for stability reasons it is crucial to formulate the dynamical equations in an explicit form ${\dot{f}}_{j}^{k}=F[v]$, where F is only a functional of the variables $v=\{{p}_{n},{q}_{n},{f}_{n}^{k},{h}_{n}^{k}\}$ without reference to their time derivatives. Numerical inversion in order to bring the system into explicit form is however prohibitive (unless the number of modes is small, for instance in the case of the Wilson discretization [37], which is not adapted to study the bath dynamics), as it would cost ${({N}_{\mathrm{modes}}\times {N}_{\mathrm{cs}})}^{3}$ operations. It turns out that a very convenient trick allows to make the inversion in ${({N}_{\mathrm{cs}})}^{6}$ operations, which is favorable provided ${N}_{\mathrm{cs}}\ll {N}_{\mathrm{modes}}$, as is the case for very long Josephson arrays or broadband environments. First, we invert equations (7) and (8) by expliciting the ${\kappa }_{{nj}}$ dependence:

$$\begin{eqnarray}&&{\dot{p}}_{i}=\displaystyle \sum _{{jn}}{A}_{{ijn}}{\kappa }_{{nj}}+{G}_{i},\ \ \ \ {\dot{f}}_{i}^{s}=\displaystyle \sum _{{jn}}{B}_{{ijns}}{\kappa }_{{nj}}+{H}_{{is}},\end{eqnarray} \tag{ 11 }$$

defining the compact notation ${A}_{{ijn}}=\tfrac{1}{2}{M}_{{ij}}^{-1}{p}_{n}\langle {f}_{j}| {f}_{n}\rangle$, ${G}_{i}={\sum }_{j}({M}_{{ij}}^{-1}{P}_{j})$, ${B}_{{ijns}}=\tfrac{1}{2}{({N}^{-1})}_{{ij}}{p}_{n}{f}_{n}^{s}\langle {f}_{j}| {f}_{n}\rangle -{\sum }_{{lm}}{({N}^{-1})}_{{il}}{A}_{{mjn}}{f}_{m}^{s}\langle {f}_{l}| {f}_{m}\rangle$, and ${H}_{{is}}={\sum }_{j}{({N}^{-1})}_{{ij}}{F}_{j}^{s}-{\sum }_{{jm}}{({N}^{-1})}_{{ij}}{G}_{m}{f}_{m}^{s}\langle {f}_{j}| {f}_{m}\rangle$, with ${P}_{j}=-{\rm{i}}\partial E/\partial {p}_{j}^{* }$, ${F}_{j}^{k}=-{\rm{i}}\partial E/\partial {f}_{k}^{j* }\,+\,(1/2)[{P}_{j}{p}_{j}^{* }\,+\,{P}_{j}^{* }{p}_{j}]{f}_{k}^{j}$, and the overlap matrices ${M}_{{jm}}=\langle {f}_{j}| {f}_{m}\rangle$ and ${N}_{{jm}}={p}_{m}\langle {f}_{j}| {f}_{m}\rangle$. The evolution equation (11) is still not in explicit form, because the parameters ${\kappa }_{{nj}}$ explicitly depend on time derivatives in equation (9). However, we can now use equation (11) to replace all the ${\dot{f}}_{i}^{k}$ terms in equation (9), which gives a closed equation for the κ matrix:

$$\begin{eqnarray}{\kappa }_{{im}} & = & \displaystyle \sum _{{jns}}[({f}_{i}^{s* }-2{f}_{m}^{s* }){B}_{{ijns}}{\kappa }_{{nj}}+{f}_{i}^{s}{B}_{{ijns}}^{* }{\kappa }_{{nj}}^{* }]\\ & & +\displaystyle \sum _{s}[{f}_{i}^{s* }{H}_{{is}}+{f}_{i}^{s}{H}_{{is}}^{* }-2{f}_{m}^{s* }{H}_{{is}}].\end{eqnarray} \tag{ 12 }$$

Inverting this linear systems with ${({N}_{\mathrm{cs}})}^{2}$ parameters provides the final ${({N}_{\mathrm{cs}})}^{6}$ scaling of our algorithm. Of course the gain is important only provided that N_cs stays small during the unitary time evolution. We show now that the dynamics indeed converges rapidly for a surprisingly small number of coherent states, demonstrating that the decomposition (5) is physically well motivated.