Resolution of superluminal signalling in non-perturbative cavity quantum electrodynamics

Recent technological developments have made it increasingly easy to access the non-perturbative regimes of cavity quantum electrodynamics known as ultrastrong or deep strong coupling, where the light–matter coupling becomes comparable to the bare modal frequencies. In this work, we address the adequacy of the broadly used single-mode cavity approximation to describe such regimes. We demonstrate that, in the non-perturbative light–matter coupling regimes, the single-mode models become unphysical, allowing for superluminal signalling. Moreover, considering the specific example of the quantum Rabi model, we show that the multi-mode description of the electromagnetic field, necessary to account for light propagation at finite speed, yields physical observables that differ radically from their single-mode counterparts already for moderate values of the coupling. Our multi-mode analysis also reveals phenomena of fundamental interest on the dynamics of the intracavity electric field, where a free photonic wavefront and a bound state of virtual photons are shown to coexist. Quantum Rabi model is a standard tool for describing cavity quantum electrodynamics, but the potential shortcomings of its single-mode version are usually neglected. Here, the authors show that, in the ultrastrong coupling regime, a multimode Rabi model is mandatory in order to avoid unphysical results.

From the theoretical side, the investigation of these nonperturbative regimes proceeded through the analysis of archetypical Hamiltonians, adapted to model different physical implementations and parameter regimes. The quantum Rabi model, describing a single two-level system (TLS) coupled to a single mode of the electromagnetic field, stands out as the simplest and the most iconic of them. Presently well understood for arbitrary values of the coupling 20 , it has been successfully employed to model the first observation of strong coupling 21 and, with some tweaks, of deep strong coupling 16 . Its mathematical properties 22 and the possible implementations with synthetic models 23,24 have also become object of interest.
To what extent any particular physical implementation is faithfully described by the quantum Rabi model depends largely upon how well it satisfies two assumptions: the emitter behaves effectively as a TLS, and only a single mode of the electromagnetic field significantly couples with it. The validity of the latter assumption is far from universal, and it has often been recognized that when the coupling is large enough to significantly hybridize the emitter with higher-lying photonic modes, those should be included in the Hamiltonian description 17,[25][26][27][28][29][30][31][32] .
The first major result of this paper will be to show, exploiting a simple gedanken experiment, that, at least in the case of cavities with an harmonic multi-mode structure, there is actually an intrinsic problem in the description of a emitter-cavity system in terms of the single-mode quantum Rabi model, which becomes unphysical in the deep strong coupling regime since it allows for superluminal signalling. In order to better understand the practical relevance of such a problem, we will then perform a rigorous analysis of the multi-mode version of the quantum Rabi model, exploiting both numerical and analytical approaches. Such analysis will reveal that the failure to consider higher-lying photonic modes has a profound impact already in the ultrastrong coupling regime, that is, for values of the coupling nowadays routinely achieved in experiments. So far, such observations have mainly consisted of transmission experiments probing the low-energy spectrum of the system 33 . It is worth noticing that, in the kind of systems we are focussing on, one can obtain a low-energy spectrum of the single-mode description that does not differ greatly from the full, multi-mode case if one uses distinct fitting parameters. However, in contrast to these previous works, our analysis reveals that the different nature of the eigenstates and their degeneracy have critical consequences on the system dynamics.

Results
The problem of superluminal signalling. We will focus most of our discussion on the simple physical system sketched in Fig. 1a: a perfect, one-dimensional cavity of length L coupled to a single TLS of frequency ω x placed at its centre. When only the coupling to the lowest mode of frequency ω c = πc/L is considered, such a system is perfectly described by the standard Rabi Hamiltonian (we take hereafter ħ = 1): In order to see how this Hamiltonian allows for superluminal signalling when g ' ω x ; ω c , let us consider the situation sketched in Fig. 1b, with an observer placed close to one of the mirrors and the system initialized in a factorized state, with the TLS either in its ground g j i or excited e j i energy level and the cavity field in its vacuum state. Such a configuration can be prepared performing only local operations on the TLS, i.e. by non-adiabatically switching on its coupling to the cavity 34,35 .
After a timescale τ R ≈ 2πg −1 , the Hamiltonian in Eq. (1) will lead to an evolution of the cavity field, conditional on the initial state of the TLS. The cavity mode is delocalized along the cavity and the observer can thus, measuring the local field, acquire an information on the initial state of the TLS, placed at a distance L 2 . Unless τ R ) L 2c , the observer can thus measure the state of the TLS, placed at a distance L 2 , in a time smaller than L 2c . The above inequality can be expressed in terms of coupling and bare frequencies as ω c ) g, showing that the parameter regime in which superluminal signalling becomes possible coincides with the non-perturbative coupling regimes of cavity QED.
Multi-mode quantum Rabi model. In order to better understand the impact of the single-mode approximation, we will study the same model of Eq. (1) but now considering the full, real-space Fig. 1 The problem of superluminal signalling in the single-mode Rabi model. a Schematic view of a qubit embedded in a perfect 1D cavity, together with the depiction of the three lowest cavity modes. When the qubit is only coupled to the fundamental mode, the system is described by the Rabi Hamiltonian. b Violation of relativistic causality by the single-mode Rabi model in regimes where g ≈ ω c . An observer placed close to the cavity edge can retrieve information about the initial state of the TLS before light is able to reach its position. c A multi-mode description is able to capture the spatio-temporal structure of the light field necessary to comply with causality electric field inside the cavity: where we have taken into account a single relevant polarization along the z axis. Here A is the transverse area of the cavity, and without any loss of generality, we have taken periodic boundary conditions to simplify the numerical analysis. By defining the symmetric modes: the dipolar coupling interaction H int = −d · E, where the dipole operator is d = μσ x u z , yields the multi-mode Rabi Hamiltonian: ðn þ 1Þω c a y n a n À i ffiffiffiffiffiffiffiffiffiffiffi n þ 1 p gσ x a n À a y with g ffiffiffiffiffiffiffi ffi 2ω c p μ= ffiffiffiffiffiffiffiffiffiffiffiffi 2ε 0 LA p and N the total number of modes included in the description. Equation (4) is well defined in the electric dipolar approximation and the low-energy part of its spectrum converges in the limit of an ideal multi-mode cavity N → ∞, when the TLS frequency ω x includes the N-dependent renormalization due to the dipole self-interaction in the Power-Zienau-Woolley gauge 36,37 .
In the standard Coulomb gauge in which ω x is microscopically independent from N, convergence would require instead to consider the diamagnetic A 2 term in the Hamiltonian 27, 28 . Recent works have proved that this remains true also in the case of superconducting circuits 31,38,39 , assuring that our results are applicable also to this important class of systems. Given that we consider ω x to be an experimentally measured value, we will not explicitly mark its dependency upon N.
In general, the total number of modes N involved will depend on the specific physical implementation of the quantum Rabi model, e.g. due to the finite size of the emitter, with several tens of them being a typical figure 31 . Even for these finite values of N, computing the dynamics of Eq. (4) for large g/ω c is a computationally formidable task, because even in the groundstate each photonic mode contains a finite population of virtual photons 1 . As explained in the Methods section, we thus adopt the approach of refs. 40,41 , recasting the Hamiltonian into the form of a chain with nearest neighbour interactions, which can then be efficiently solved by using matrix product states (MPS) [42][43][44] .
System dynamics. In Fig. 2a, we plot the time evolution of the TLS population versus g/ω c , with the TLS initially in its excited state and zero photons in the cavity, ψð0Þ j i¼ e j i 0 j i, obtained, respectively, solving Eq. (1) (single-mode) and Eq. (4) (multimode). This initial configuration is a superposition of excited states of the coupled light-matter system, which could be initialized by applying a π pulse in a decoupled system and then by non-adiabatically switching on the coupling 34,35 . As an alternative approach to obtain an initial excited configuration, one could also apply a suitable pulse to the coupled system in its ground state 45 . In any case, the effects that we report here appear as long as the system is initially in some superposition of excited states. Figure 2b shows a plot along the dashed lines in Fig. 2a, corresponding to g/ω c = 0.6. It is clear that the single-mode approximation drastically fails as the system enters the nonperturbative region, with completely different physics taking place already for values of the coupling well below the boundary of the deep strong coupling regime. While for the considered values of the coupling the Rabi oscillations are distorted in the single-mode case, for the multi-mode Hamiltonian the TLS relaxes immediately and remains most of the time in a superposition of g j i and e j i yielding a population of 1/2, experiencing a sequence of sharply peaked revivals that bring it back to the excited state at times multiple of the cavity roundtrip time, 2π/ω c . Even for lower values of the coupling-before these revival peaks are fully formed-one can observe a perturbation of the Rabi oscillations taking place at those specific times. In Fig. 2c, Single mode (SM) 10 we plot the amplitude of the electric field inside the cavity, x ∈ (−L/2, L/2) as a function of time, for the case g/ω c = 0.6. The electric field features the coexistence between two distinct components; (i): a localized cloud bound at the position of the TLS, and (ii): a free wavefront propagating at the speed of light. The free wavefront is backscattered at the edges of the cavity and returns at the position of the emitter at times 2πn/ω c , when all the light is perfectly absorbed by the TLS-see inset of Fig. 2cyielding the revival peaks in its population.
In order to gain further insight into the dynamical features of the multi-mode quantum Rabi model in the non-perturbative regime, we perform now an analysis similar to the one applied in ref. 2 to the single-mode case. To do so, we split the Hamiltonian into two parts, H = H I + H II , with H II ¼ ω x 2 σ z , and start by studying the action of H I alone. While in the single-mode case neglecting H II is a good approximation only in the limit ω x ≈ 0 of the deep strong coupling regime 2 , we will show that it is enough to describe the features that we have reported for the multi-mode model even at the resonant condition ω x ≈ ω c and in the ultrastrong coupling regime. Let us consider that H I is acting on a wavefunction whose matter component is one of the eigenstates of σ x , ± j i. In that case, H I takes the form of a collection of driven harmonic oscillators: ðn þ 1Þω c a y n a n Ç i ffiffiffiffiffiffiffiffiffiffiffi n þ 1 p g a n À a y The evolution under this Hamiltonian can be readily solved by means of a unitary transformation where D n β ð Þ ¼ exp βa y n À β Ã a n Â Ã is a displacement operator acting on mode n with a sign that depends on the state of the TLS, and β 0 = ig/ω c . This transformation gives a Hamiltonian without the driving term, n¼0 ðn þ 1Þω c a y n a n À g 2 =ω c Â Ã . We can write the evolution of an initial state with no photons ψð0Þ j i ± = Q n 0 j i n ± j i under the effect of H I as: where Çξ N ðtÞ j i≡ Q NÀ1 n Çβ n ðtÞ . Here β n ðtÞ represents a coherent state in the nth cavity mode, with β n (t) given by: The corresponding trajectories in phase space for each cavity mode are depicted in Fig. 3. The single-mode case was already introduced in ref. 2 ; it features circular trajectories corresponding to oscillations around the centre of an harmonic oscillator displaced by β 0 . The period of these oscillations is given by 2π/ω c , and it is associated with the revivals in the probability of the initial state, corresponding to those times when the state in phase space crosses the (0, 0) point. In the multi-mode case, this picture is extended, with each mode of frequency ω c (n + 1) following a circular trajectory, whose radius and period depend on n as 1= ffiffiffiffiffiffiffiffiffiffiffi n þ 1 p and 1/(n + 1), respectively. With high-energy modes oscillating faster than low-energy ones, the total period of the dynamics is fixed, as in the single-mode case, by the period of the fundamental mode, ω c . The revival probability of the initial state is given by: If the TLS is initially in an excited state e j i = þ j i À À j i ð Þ = ffiffi ffi 2 p , as in the case we solved numerically, the resulting wavefunction consists of a superposition: with a revival probability given as well by Eq. ( . The exponent is given by a sum that diverges logarithmically with N for all t except for t = 2πn/ω c : This means that the overlap decays quickly to some stationary value O N that goes to zero with increasing N as O N % 1= 2e γ ðN þ 1Þ ½ 4g 2 =ω 2 c (with γ the Euler-Mascheroni constant) and then experiences sharp revivals at multiples of the cavity roundtrip time. In contrast to the single-mode case, where the width of the revival peaks is given by g/ω c , these decays and revivals occur on a short timescale τ ≈ 2π/(Nω c ), which justifies the approximation of neglecting H II as long as (i): the decay is fast enough, Nω c ) ω x ; and (ii): the stationary value of the overlap after the decay is small enough, ω x O N ( g. This sets two conditions on N and g for the multi-mode physics to become relevant and the effect of light propagation that we report to manifest, breaking down the single-mode Rabi physics. We have observed that, for ω x = ω c , values of N ∈ [10, 100] and g=ω c ⪆ 0.25 are sufficient to fulfil these conditions, meaning that these effects will be relevant already in the ultrastrong coupling regime for systems involving only several tens of cavity modes. A more detailed analysis of the implications of a finite N is provided in Supplementary Notes 1 and 2. Interestingly, these results show that the multi-mode Rabi model can work as a dynamical description of wavefunction collapse based only on the Schrödinger equation. This is related to previous efforts [46][47][48] , which, in the spirit of the many-worlds theory, describe the wavefunction reduction as a unitary evolution that includes the measurement device as part of the quantum system 49,50 .

Multi mode
As we showed before numerically, the revivals can also manifest in the population of the TLS, which within our approximation is trivially related to the overlap O N (t) as: This expression reproduces perfectly the extremely sharp revival profiles that we report in Fig. 2b that were numerically computed for N = 50. Furthermore, it is easy to show how the collection of circular trajectories of the multi-mode case gives rise to the spatial profile of the electric field that we obtained numerically. The amplitude of the electric field is given by: which, when plotted, shows a perfect agreement to the profile in Fig. 2c. This is explicitly shown in Fig. 4a, which depicts a comparison between numerical calculations and Eq. (12) at a given time. Equation (12) can be decomposed into a timedependent term, corresponding to (i) the part of the field that is emitted from the TLS and propagates freely towards the ends of the mirror, and (ii) a time-independent term, corresponding to the part of the field that remains bound to the TLS at the centre of the cavity. These terms have their origin in the time-dependent and -independent parts of the coherent amplitude β n (t) of each of the cavity modes, see Eq. (7), and the ratio between them will depend on the initial state (being 1/2 in our particular case).
Propagative and bound photons. The plot of the electric field in Fig. 2c seems to clearly attribute the regular peaks in Fig. 2a, b with period 2π/ω c to a rather trivial propagative effect of photons bouncing back and forth, and as such it had already been described in ref. 51 within the rotating wave approximation, which a priori excludes the presence of any non-perturbative effect. Still our analysis shows that those peaks have the same origin as those reported in ref. 2 for the single-mode quantum Rabi model in the deep strong coupling regime, in which, of course, the concept of propagation is non-relevant. Here we have shown that these two seemingly unrelated phenomena are effectively the same and that, in the multi-mode case, it is intimately related to light propagation and thus relativistic causality. This provides an intuitive physical understanding of why this phenomenon manifests at much lower coupling rates than actually predicted by the single-mode model: it is linked to a propagation that cannot be neglected when the coupling frequency becomes comparable to the cavity roundtrip, since it would allow for superluminal signalling.
In order to understand the second component of the dynamics (the localized cloud bound at the position of the TLS), let us recall that we expressed the Hamiltonian as a collection of displaced harmonic oscillators. Therefore, the absolute value of the timeindependent part of β n (t) describes a coherent state at the equilibrium position of the nth displaced oscillator, β 0 = ffiffiffi n p , i.e. its vacuum state. We can then understand the time-independent part of the wavefunction as a set of displaced oscillators in vacuum, which corresponds to the ground state of the system. We have verified this by numerically computing the ground state using imaginary-time evolution, see Fig. 4a. The results obtained confirm that the ground state of a TLS non-perturbatively coupled to a cavity is indeed constituted by a localized cloud of photons around the TLS, which is in a superposition with a population corresponding to that observed in the revivals n σ h i ¼ 1=2. Those virtual photons have been demonstrated to exist also in lossy systems 52 , although once the coupling with the environment is properly considered 35, 53-56 their non-radiative nature becomes apparent. Our results provide a more transparent way to understand them as a localized, bound state of photons; in future works, the methods that we present here might be applied to study their properties in lossy systems. Bound states have already been documented in the context of ultrastrong coupling of a quantum emitter to open lines 57,58 , and there is much literature discussing their existence in boson impurity models in the single photon 59, 60 and, more relevant to our discussion, multiphoton case 61 . They are associated with eigenstates of the to the ground state of the system for g/ω c = 0.6. Solid, blue (dasheddotted, yellow): numerical (analytical) calculation of the electric field for an initial state e j i 0 j i after evolving for a time t = π/2ω c , confirming that the dynamics of the system is given by the independent evolution of two freely propagating wavepackets plus a localized cloud of photons corresponding to the ground state of the light-matter system. b Low-energy spectrum of the single-mode (red) and multi-mode (blue) Rabi Hamiltonian as a function of the coupling rate. For each value of g, the eigenvalues are expressed with respect to the ground state. Here, ω x = ω c system whose energy lie outside the energy spectrum of the bath, which in this case would be constituted by the infinite set of cavity modes.
For a given set of parameters, the spectrum of eigenvalues obtained from the multi-mode Rabi Hamiltonian strongly differs from the result given by the single-mode one, see Fig. 4b. In the large-coupling limit, both models feature a series of equispaced energy levels similar to the bare ones, a result shown above in the derivation of H ′ I and well known for the single-mode case 20,27,62,63 . However, the results predicted by both models differ substantially in a range of couplings approximately delimited by 0.1 ≲ g/⍵ c ≲ 2 for the low-energy eigenstates. The results shown in Fig. 4b evidence that transition energies should be fitted with a multi-mode Rabi Hamiltonian in order to obtain a proper description of the system; the use of a single-mode Rabi Hamiltonian might lead to a qualitatively similar prediction for the low-energy transitions but yielding an incorrect estimation of the system parameters. Owing to this possibility, an unambiguous evidence of the breakdown of the single-mode Rabi model physics enforced by causality should come from the analysis of the dynamics of observables, such as the TLS population, that, as we have shown, carry unequivocal signatures of the propagation of light inside the cavity.

Discussion
We have performed a thorough theoretical analysis of a single emitter coupled to a photonic resonator. Our first result has been that, at least for resonators with harmonic spectra, like standard λ/2 cavities, the single-mode quantum Rabi model is incompatible with relativistic causality. By means of quasi-exact numerical calculations using MPS, we have then studied the multi-mode version of the quantum Rabi model confirming that, beyond certain values of the coupling rate, the single-mode model fails to describe the physics of a TLS coupled to the electric field inside a cavity. The failure of the model occurs in the regime of ultrastrong coupling, well before reaching the limit of deep strong coupling, and where the single-mode Rabi model is often invoked. This failure does not only manifest in the spectrum of eigenvalues, which differs from the one given by the single-mode model, but most importantly in the dynamics, which features freely propagating photonic wavepackets inside the cavity that coexist with a bound state of virtual photons corresponding to the ground state of the system.
Our theoretical analysis is most timely. Advances in superconducting circuits in fact not only recently led to the first observation of the deep strong coupling regime in a single-mode setup 16 , but multi-mode effects in the ultrastrong coupling have also been recently reported 17 . Although this work primarily deals with the failure of the single-mode approximation, we verified that our results are not qualitatively affected by the breakdown of the TLS approximation. In Supplementary Note 3, we in fact extend our investigations beyond the quantum Rabi model, considering as matter degree of freedom a bosonic field with a small Kerr nonlinearity. We found that, in this situation, although higher modes are also involved in the dynamics, our conclusions remain valid.
These results bring a deeper understanding of a system of central importance in quantum mechanics and therefore are very relevant for the design of new technologies aiming to exploit the physics of light-matter coupling in the ultrastrong coupling regime.

Methods
Computation of system dynamics with MPS. We make use of the approach presented in refs. 40,41 and define a new set of operators by means of an unitary transformation b i ¼ P N n¼0 U i;n a n to recast the Hamiltonian in Eq. (4) into another with nearest neighbour interactions: with U i;n ≡ ffiffiffiffiffiffiffiffiffiffiffi n þ 1 p Q i ðn; 1; 0; NÞρ À1 i (Q i being the Hahn polynomials); t i ≡ −A i ρ i +1 /ρ i ; ω i ≡ 1 + A i + C i and A i ¼ ði þ 2Þ 2 ðN À iÞ 2ði þ 1Þð2i þ 3Þ ; ð15Þ where we used the Pochhammer symbol (z) i = z(z + 1)…(z + i − 1). Writing the Hamiltonian in this form allows us to compute its dynamics very efficiently using the MPS method.
Data availability. The data that support the findings of this study are available from the corresponding author upon request.