Non-perturbative treatment of giant atoms using chain transformations

Superconducting circuits coupled to acoustic waveguides have extended the range of phenomena that can be experimentally studied using tools from quantum optics. In particular giant artificial atoms permit the investigation of systems in which the electric dipole approximation breaks down and pronounced non-Markovian effects become important. While previous studies of giant atoms focused on the realm of the rotating-wave approximation, we go beyond this and perform a numerically exact analysis of giant atoms strongly coupled to their environment, in regimes where counterrotating terms cannot be neglected. To achieve this, we use a Lanczos transformation to cast the field Hamiltonian into the form of a one-dimensional chain and employ matrix-product state simulations. This approach yields access to a wide range of system-bath observables and to previously unexplored parameter regimes.


I. INTRODUCTION
Quantum optical theory provides a solid framework for the study of light-matter interaction.Yet paradigmatic models are based on several approximations, such as the rotating-wave, electric dipole and Born-Markov approximations [1].While the underlying assumptions are typically well justified, recent experimental advances have paved the way for investigations of yet unexplored parameter and physical regimes.Superconducting circuits offer a versatile platform for such studies in which artificial atoms may be efficiently and strongly coupled to electromagnetic and sound waves [2].In particular giant atoms challenge standard approximations and can only be accurately described when taking the finite spatial extent of the artificial atom into account [3], plus a finite propagation speed if coupled to sound waves [4,5] and counter-rotating terms beyond the rotating-wave approximation (RWA) at strong couplings.Recent work has already capitalized on this and demonstrated several intruiging effects that occur in giant atomic setups, including decoherence-free interactions [6], non-exponential atomic decay [7], oscillating bound states [8] and chiral atom-waveguide couplings [9].Still theoretical treatment has so far been restricted to couplings in the realm of the RWA.
At elevated light-matter couplings several physical phenomena can only be accurately captured by taking multiple field modes into account [10][11][12].In this regime unphysical properties of single-mode models become more apparent such as, e.g., causality violations in the form of superluminal signaling [13][14][15][16].In contrast to the single-mode quantum Rabi model (QRM) [17], the corresponding multimode problem is not known to be integrable and requires novel techniques for theoretical treatment.The regime where the coupling strength becomes comparable to the bare resonance frequencies in the system is referred to as the ultra-strong coupling (USC) regime [18].Previous works have established matrix-product state (MPS) simulations as a means to explore quantum optics phenomena of small atoms in the USC regime [19], and they have proven useful for the study of non-Markovian light-matter interactions [20,21].While the USC regime is becoming more and more experimentally accessible, its theoretical study still requires improved analytical and numerical methods, making it a timely research topic.Moreover, at even stronger couplings and within the deep and extremely strong coupling regimes, other non-perturbative methods become available again [22,23].
Here we investigate the low-energy physics and the dynamics of giant atoms beyond the RWA, in the USC regime and with multimode interactions, using a numerically exact, nonperturbative approach.We model the giant atoms as two-level systems.The coupling points we model by a profile function with a finite width thus suppressing the coupling to high frequency modes and motivating a natural UV cutoff.Apart from this UV cutoff our approach requires no further approximations of the model Hamiltonian.While our approach is general, we mainly focus on superconducting qubits coupled to acoustic field modes and the resulting non-Markovian effects which are due to a finite speed of sound.In particular, we investigate the dynamics of a single giant atom coupled to an acoustic waveguide with intrinsic time delay, thus extending the analysis of previously predicted oscillating bound states [8] beyond the single-excitation subspace.Our theoretical treatment of the system-reservoir interaction relies on a so-called chain (also star-to-chain or Lanczos) transformation.This unitary transformation casts the field into the shape of a linear harmonic chain, which is particularly suited for numerical simulation.Rooting back to the numerical renormalization group [24], these methods are widely used in the study of open quantum systems (e.g., see [25][26][27]), but have also proven useful in quantum optics as seen, for example, in [16,[28][29][30].We follow a similar numerical approach as [16], which allows us to go beyond the single-excitation subspace and numerically study system and bath observables using MPS [31,32].
This work is structured as follows.In Sec.II, we introduce the setup and theoretical model of our study.We show how the underlying Hamiltonian can be cast into a form amenable to an efficient numerical analysis even at strong coupling and beyond the RWA, using chain-mapping techniques.In contrast < l a t e x i t s h a 1 _ b a s e 6 4 = " c s h L p q c r y 9 5 r V s c M 5 X C n B 9 A 9 m m k = " > A A A C 6 X i c h V F N S x x B E H 1 O 4 v f X x h y F M L g I H m S Z l U U 9 S j T B i 8 S A q 4 I j S 8 / Y z j Y 7 X / T 0 C r p 4 y x / w F n L N z a v + G f 0 t O e R N O w o q Y g 8 9 V f 3 q 1 e u q r i C P V W E 8 7 2 7 I + f B x e G R 0 b H x i c m p 6 Z r b 2 a W 6 / y P o 6 l O 0 w i z N 9 G I h C x i q V b a N M L A 9 z L U U S x P I g 6 G 2 W 8 Y M z q Q u V p X v m P J f H i Y h S d a p C Y Q h 1 a l / 8 H 4 m M h O s v + 1 1 h / E J F i e j 4 i T B d n Q w u L j u 1 u t f w 7 H J f O 8 3 K q a N a u 1 n t H j 5 O k C F E H w k k U h j 6 M Q Q K f k d o w k N O 7 B g D Y p q e s n G J S 0 w w t 0 + W J E M Q 7 f E f 8 X R U o S n P p W Z h s 0 P e E n N r Z r p Y 5 P 5 u F Q O y y 1 s l / Y L 2 H / e F x a I 3 b x h Y 5 b L C c 9 q A i u N W c Y e 4 Q Z e M 9 z K T i v l Y y / u Z Z V c G p 1 i 3 3 S j W l 1 u k 7 D N 8 0 t l i R B P r 2 Y i L b 5 Y Z U S O w 5 z O + Q E r b Z g X l K z 8 q u L b j E 1 p h r b Q q a a U o q K d p y 9 d n P R x z 8 + V Q X z v 7 K 4 3 m a q P 1 s 1 X f + F o N f A z z W M A S p 7 q G D W x j l 3 W E + I V r 3 O D W 6 T l X z m / n z w P V G a p y P u P Z c v 7 + B 2 e c m 5 Y = < / l a t e x i t > ⌦ ˆ z < l a t e x i t s h a 1 _ b a s e 6 4 = " Y f i R A i q l F o / t v O o U d 8 w 9 U R P 4 W B k = " > A A A C / n i c h V H L a h R B F D 1 p H 3 n 4 m u g y m 8 J B C C h N d w z R 5 e A L N 0 I E J w m k 4 1 B d U + l p p l 9 U 1 w Q m Q 8 D f 8 A f c i V t 3 2 e o n 6 L d k k d N l R 0 h C S D X V 9 9 a 5 5 5 6 6 t 2 5 c Z W l t g + D P n H f j 5 q 3 b 8 w u L S 3 f u 3 r v / o L P 8 c K s u J 0 b p v i q z 0 u z E s t Z Z W u i + T W 2 m d y q j Z R 5 n e j s e v 2 7 i 2 w f a 1 G l Z f L L T S u / l M i n S / V R J S 2 j Q e R 7 F 2 s p B I K J n 0 U h a o Q b B 5 2 g o k 0 Q b 0 Q K h e C q i X N q R y W c j X / l H g 0 4 3  x w H r i P A F P 3 G D X 9 6 1 9 9 X 7 5 n 3 / H + q N V T l L u L e 8 H / 8 A h w u h g A = = < / l a t e x i t > e j j m 0 s + h / n Z 2 Z o q l + e L s F u c 9 h 4 / V h T G M Y 5 p T X c A S 1 r H J O j y c 4 h 4 P e L S q 1 p V 1 Y 9 1 + U K 2 W P G c E 3 5 Z 1 9 w 6 7 k p n Y < / l a t e x i t > |gi < l a t e x i t s h a 1 _ b a s e 6 4 = " p w / O e w 8 f q w h j G M c 2 p L m A J 6 9 h k H R 5 O c Y 8 H P F p V 6 8 q 6 s W 4 / q F Z L n j O C b 8 u 6 e w e 2 t J n W < / l a t e x i t > |ei < l a t e x i t s h a 1 _ b a s e 6 4 = " x A t V b t e j j m 0 s + h / n Z 2 Z o q l + e L s F u c 9 h 4 / V h T G M Y 5 p T X c A S 1 r H J O j y c 4 h 4 P e L S q 1 p V 1 Y 9 1 + U K 2 W P G c E 3 5 Z 1 9 w 6 7 k p n Y < / l a t e x i t > |gi < l a t e x i t s h a 1 _ b a s e 6 4 = " p w / O e w 8 f q w h j G M c 2 p L m A J 6 9 h k H R 5 O c Y 8 H P F p V 6 8 q 6 s W 4 / q F Z L n j O C b 8 u 6 e w e 2 t J n W < / l a t e x i t > |ei < l a t e x i t s h a 1 _ b a s e 6 4 = " x A t V b t e j j m 0 s + h / n Z 2 Z o q l + e L s F u c 9 h 4 / V h T G M Y 5 p T X c A S 1 r H J O j y c 4 h 4 P e L S q 1 p V 1 Y 9 1 + U K 2 W P G c E 3 5 Z 1 9 w 6 7 k p n Y < / l a t e x i t > |gi < l a t e x i t s h a 1 _ b a s e 6 4 = " p w / O e w 8 f q w h j G M c 2 p L m A J 6 9 h k H R 5 O c Y 8 H P F p V 6 8 q 6 s W 4 / q F Z L n j O C b 8 u 6 e w e 2 t J n W < / l a t e x i t > |ei … < l a t e x i t s h a 1 _ b a s e 6 4 = " v c 0 g z 6 < l a t e x i t s h a 1 _ b a s e 6 4 = " I z w N Z t 8 3 p c t E 0 / O e w 8 f q w h j G M c 2 p L m A J 6 9 h k H R 5 O c Y 8 H P F p V 6 8 q 6 s W 4 / q F Z L n j O C b 8 u 6 e w e 2 t J n W < / l a t e x i t > |ei < l a t e x i t s h a 1 _ b a s e 6 4 = " v 5 x q l X U s I 9 D 1 u H j B n 9 w j 7 9 W Z N 1 a v 6 y 7 I d U q j X K W 8 W p Z v 5 8 B 9 G S g d w = = < / l a t e x i t > x q l X U s I 9 D 1 u H j B n 9 w j 7 9 W Z N 1 a v 6 y 7 I d U q j X K W 8 W p Z v 5 8 B 9 G S g d w = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x g a G < l a t e x i t s h a 1 _ b a s e 6 4 = " c s h L p q c r y 9 5 r V s c M 5 FIG. 1. Setup and schematic illustration of chain transformation.(a) Giant atom with M = 3 equidistantly spaced coupling points at distance τ , coupled to a periodic waveguide of length L. Emitter-waveguide couplings are locally described by Gaussian smearing functions f s in Eq. ( 6).(b) Chain transformation maps system-reservoir model with a one-to-all coupling to a linear chain as described by chain parameters α1,...,N and β1,...,N−1 in Eq. ( 15), with Kj = σx ⊗ (fj âj + f * j â † j ).(c) Illustration of a two-atom setup with braided coupling points, which can be chain transformed into a chain with next-to-nearest-neighbor couplings, as indicated in (d).
to earlier works, our description does not rely on the assumption of a point-like emitter-bath coupling, but we promote coupling points to smeared coupling functions with a finite spatial support.We provide estimates for the required values of all characteristic system parameters of an experimental implementation using superconducting circuits at the end of Sec.II.In Sec.III, we present an analysis of the low-energy physics of the system.In particular, we discuss elementary excitations of the ground state as a function of increased emitter-reservoir coupling strength, in analogy with the well-understood quantum Rabi model.In Sec.IV we present a study of the temporal dynamics of a single giant atom coupled to an acoustic waveguide, with an intrinsic time delay, at three coupling points.We showcase and discuss the implications of the breakdown of the RWA at strong coupling.A stability analysis of the findings with respect to experimentally relevant parameters is provided.Finally, we conclude our work in Sec.V, discuss possible future research directions and highlight the wide-ranged applicability of our approach, e.g., to systems with multiple giant atoms and multilevel emitters.The latter is particularly important to realistic implementations in which, depending on the chosen gauge, the two-level approximation is no longer applicable for sufficiently strong couplings [33][34][35].Note that we use natural units ( , c = 1) throughout this work.

II. SETUP AND THEORETICAL FRAMEWORK
In this section, we present our theoretical framework and introduce the chain transformation that we employ for the study of stationary (see Sec. III) and dynamical (see Sec. IV) properties of two-level emitters coupled to a waveguide.
Setup.-A schematic illustration of the setup and the chain transformation is provided in Fig. 1.We treat a single quan-tum emitter as a two-level system coupled to the waveguide modes at M coupling points, cf.Fig. 1(a).For simplicity, we focus on equidistantly spaced coupling points, with a nonzero, significant propagation time τ between neighboring coupling points.Such a system may be realized with a superconducting qubit piezoelectrically coupled to an acoustic waveguide at several locations [7].The interaction between emitter and waveguide modes is usually described by a one-toall coupling, i.e., the emitter couples to all non-interacting field modes.Once brought into the form of a linear chain, cf.Fig. 1(b), well-developed techniques based on MPS can be utilized for efficient numerical studies of various system and bath observables.Note that for setups with multiple emitters, where n emitters couple to one waveguide as schematically depicted in Fig. 1(c), the chain transformation, as reviewed in App.C, casts the field into a linear chain, with each mode coupling to its n nearest neighbors as indicated in Fig. 1(d).
As mentioned, here we model the atom as a two-level system, i.e., we use the two-level approximation (TLA).For couplings above the weak coupling regime, the validity is known to be highly gauge dependent [33][34][35] and only specific gauges still allow for the TLA to be applied beyond weak coupling.Here we chose the TLA interaction Hamiltonian akin to the dipole gauge which, for the quantum Rabi model was found to perform reasonably well in the USC regime [33].
Model.-The total Hamiltonian can be decomposed as the sum of the atomic, the field, and the interaction Hamiltonian, Ĥtot = ĤA + Ĥf + Ĥint . ( Assuming a two-level emitter with frequency Ω, and a massless field in a periodic cavity of length L described by modes with wavenumbers k j = 2πj/L, the non-interacting terms in (1) can be written as with the ground and excited states of the emitter, |g and |e , and the annihilation (creation) operator a ( †) j of field mode j.

The interaction Hamiltonian reads
where λ is a dimensionless coupling constant and π denotes the field momentum.The smearing function f (x) models the emitter-waveguide coupling and has the dimensions of a density.For a giant-atom setup as shown in Fig. 1(a), where M > 1, we consider a sum of single-point couplings of the form with coupling points centered around the positions x 1 , ..., x M .The shape of f s (x) may not directly correspond to the physical shape of a given coupling point, but should be chosen to correctly capture the frequency dependence of the coupling strength (see Eq. ( 9)).In this work, each coupling point is described by a Gaussian profile function with 2d being the effective diameter of each coupling point and dx f s (x) = 1.Other choices for f s (x) can equally be considered, and some examples are discussed in App.B. Note that the choice f s (x) = δ(x) results in a UV divergent coupling which, however, does not occur in physical models [36].Field modes.-Thefield momentum operator π(x), which is equal to the time derivative ∂ t φ(x) of the field amplitude, expressed in terms of field eigenmodes, reads Hence, we can rewrite the interaction Hamiltonian as The coefficients f j for a giant atom with equidistant coupling points follow straightforwardly from the coefficients 2L dx e ikj x f s (x) for a single coupling point.For example, for a giant atom with three coupling points at x l = −τ, 0, τ (l = 1, 2, 3), we find For the Gaussian profile (6), replacing the integral The behavior of |f j | is shown in Fig. 2 for emitters coupling to the waveguide with this smearing function through M = 1 and M = 3 points, respectively.The decay of |f j | for sufficiently large j allows us to introduce a UV cutoff and only consider a finite number of 2N field modes, i.e., we restrict the index to −N ≤ j ≤ N and also discard the zero mode, to which the atom does not couple.Note that this UV cutoff is the only simplification of the original physical model that the present approach requires.In particular, it does not rely on the rotating wave approximation (RWA) or the Wigner-Weisskopf approximation.
Chain modes.-Thedynamics of the Hamiltonian, after the UV cutoff, is straightforward to treat numerically if the coupling is weak, such that the RWA can be applied, and if one restricts attention to the single-excitation subspace of the approximate Hamiltonian.However, in order to treat many excitations within the RWA, and to study USC beyond the domain of the RWA, here we employ a chain transformation of the field modes.Such a chain transformation yields a new basis of field mode operators ĉ0 , ..., ĉ2N−1 , which we refer to as chain modes.These are related to the eigenmodes of Ĥf by a non-mixing Bogoliubov transformation, The front chain mode is chosen as ĉ0 = 1 such that the interaction Hamiltonian takes the form Using Lanczos algorithms (see App. C) the chain modes are appropriately chosen such that they cast the field Hamiltonian into the form of a harmonic chain with nearest-neighbor hopping interactions only, with real coefficients α i , β i ∈ R. Fig. 3 shows a plot of these coefficients for the setup that we will use in our numerical examples (see Tab. I).
Note that if the atom has an even profile function f (x) = f (−x) such as (10), it does not couple to the odd sector of the field modes.Then, by introducing the basis change â(±) j = (â j ± â−j ) / √ 2, we can restrict attention to the N field modes of the even sector and, accordingly, only construct N chain modes as linear combinations of even field modes.
Numerical simulations.-Thelow-energy sector of (1) can be efficiently described using MPS [31,37], once the interaction and field Hamiltonians have been cast into their respective forms (14) and (15).In the following sections, we use both density-matrix renormalization group (DMRG) and time-evolution algorithms to study the stationary and dynamical properties of giant atoms coupled to a waveguide, as a function of the coupling strength λ and the emitter frequency Ω.
Unless stated otherwise, the default configuration that we consider is that of a giant atom coupling to the chain modes at M = 3 coupling points, located at x l = −L/20, 0, L/20, each modelled by the Gaussian smearing (6), with all parameters as specified in Tab.I.The default value for the atom frequency Ω is chosen to be resonant with the 80 th field mode which features the largest coupling coefficient |f j |, as can be seen in Fig. 2.
At strong coupling λ, care must be taken to ensure that the truncation errors associated with increasing chain-mode occupation numbers ĉ † i ĉi are still negligible in the numerical calculations.We find that this is possible even deep in the USC regime, as discussed in Sec.III, using 25 bosons per site.In our numerical calculations using the iTensor software package [38], we also ensure convergence with respect to MPS bond dimension (≈ 200) and chain length N ≤ 1000 at a singular value decomposition (SVD) cutoff of 10 −12 .Experimental considerations.-Experimentally,the present system and its considered initial state can be realized and prepared, e.g., using superconducting qubits coupled to an acoustic cavity [7,[39][40][41].Based on prototypical parameters used in our calculations, cf.Table I, one may choose a qubit frequency of Ω/(2π) = 2.4 GHz.At a typical sound velocity of c = 3 km/s this yields a distance of ≈ 5 µm between the coupling points of the atom.In an implementation, instead of a periodic waveguide with length L, one can consider an open-ended waveguide of length L/2 ≈ 50 µm.These ballpark values are realistic and consistent with recent experimental implementations.As described in Secs.III and IV, we identify the onset of USC at around λ ≈ 0.15, which amounts to an acoustic qubit-waveguide coupling of ≈ 8.5 MHz per coupling point.Considering that the qubit couples through three coupling points, the total coupling between qubit and waveguide is comparable to the estimated free spectral range of ≈ 30 MHz, placing the setup in the strong-multimode regime.This value is comparable with previously reported acoustic coupling strengths, and there are various prospects for state-of-the-art experimental settings to be operated even more deeply in the strong coupling regime, e.g., by appropriate choice of material [40].

III. LOW-ENERGY SPECTRUM AND EIGENSTATES
With the approach presented above, it is possible to investigate giant atoms beyond the realm of the RWA and the singleexcitation subspace.Since the field is not traced out for an effective-system description, the approach also yields full access to field observables such as photon numbers or field energy density, and allows for, e.g., the investigation of virtual photon clouds.As a first step, here we calculate and characterize the ground state and first excited state of our giant atom setup, as a function of the coupling strength.Hereby, we explore the entire USC regime and access the onset of the deep-strong coupling (DSC) limit.
For single-mode models, such as the QRM, the USC and DSC regimes are well understood and characterized [42], and both have been achieved on several experimental platforms.In the following, we will see that the lowest energy eigenstates of our multimode model generally follow the intuition based on the single-mode QRM.Yet the onset of signatures related to the USC and DSC regimes is shifted to smaller coupling strengths, which underlines that the effective emitterfield coupling is enhanced as the emitter simultaneously couples to the field via a multitude of modes.
The absolute values for the ground state in Fig. 4(a) behave monotonically and, thus, make it difficult to distinguish different regimes.However, the energy differences plotted in Fig. 4(b) provide a richer picture: The energy gap of the Hamiltonian ∆H tot starts at ∆H tot (λ = 0) = |k 1 | = 2π/L for λ = 0, then decreases over an intermediate range of 0.5 λ 1.5, and, finally, closes at λ 1.5.The energy differences for the separate terms of the Hamiltonian behave accordingly at low and large λ, but they exhibit prominent peaks in the intermediate region, where the energy gap ∆H tot is closing most rapidly.The behavior of the spectrum in the intermediate range of λ resembles the spectrum of the single-mode QRM [42] in the USC, whereas for λ 1.5, the spectrum resembles the DSC of the QRM.In the QRM, the USC sets on when the ratio of coupling strength to emitter gap is of the order of ∼ 0.1, and DSC sets on at a ratio of ∼ 1.In our approach, analogously, the range of USC can be estimated by considering the ratio of the energy scale of the interaction Hamiltonian λ √ µ 0 to the atom's gap Ω, where we have √ µ 0 /Ω ≈ 0.69.From this comparison, one expects the USC to lie within 0.15 λ 1.5, which agrees well with our numerical findings.In contrast, the coupling strength between emitter and the resonant field mode is only |f 80 |/Ω ≈ 0.076, which wrongfully would suggest the USC regime to lie at much higher λ.This underlines that the emitter couples efficiently to many field modes, cf.Fig. 2, and that a single-mode description would fail.
Structure of eigenstates.-Wecan further investigate the structure of the obtained lowest eigenstates in the different coupling regimes, and compare them to our expectations based on the QRM, by characterizing them in terms of observables such as photon numbers and the atomic population.
In the perturbative regime, as λ → 0, clearly the ground state of the system approaches the free ground state, i.e., the product of the atom ground state and the vacuum |ψ GS → |g, 0 , and the first excited state is obtained by placing a single photon into the first free field mode |ψ ES → â † 1 |ψ GS .In fact, as shown in Fig. 10 of App.D, we find a large overlap | ψ ES | â † 1 |ψ GS | 2 for our numerically obtained eigenstates for sufficiently small coupling.
To characterize the eigenstates in USC and beyond, we consider the atomic population p e = (1 + σz )/2, the total excitation number of the field nfield = j nj , with nj = â † j âj for the jth field mode, and the occupation number of the lowest free field mode n1 .Figs. 5(a)-(c) show these expectation values in the ground state, and their increase as a function of coupling strength λ.The lower panel, Fig. 5(d)-( For the lowest values of λ, we recognize the results of |ψ GS and |ψ ES lying close to |g, 0 and â † 1 |g, 0 , respectively.For large couplings, where we saw above that the two pairs approximately form a degenerate pair, we see that the states also agree in the occupation observables.This pair is characterized by a large number of field excitations that are spread out over many field modes, and by the atom approaching half occupation, p e → 1/2.In fact, this is what one would expect in the DSC where the interaction Hamiltonian Ĥint dominates over the other parts of the total Hamiltonian Ĥtot and, thus, eigenstates of Ĥtot lie close to eigenstates of Ĥint .Eigenstates of Ĥint , however, would be given by the product of the eigenstates of σx = |e g| + |g e|, given by |±X , and eigenstates of the field operator ĉ0 + ĉ † 0 , which can be approximated well by coherent states with a (positive or negative) eigenvalue with respect to ĉ0 , thus allowing for the construction of a degenerate pair.This behavior is reminiscent of the eigenstates of the single-mode QRM.Also there, for small coupling, the ground state of the system contains no photonic excitations, and in the DSC regime, the number of photonic excitations grows linearly while the atomic population saturates at half occupation.

IV. OSCILLATING BOUND STATES
This section studies dynamical properties of giant atoms with a focus on oscillating bound states.These states were recently predicted to arise for giant atoms in the RWA regime under certain, fine-tuned conditions.In our approach we can simulate the dynamics of giant atoms far into the USC regime, and up to times set by the waveguide crossing time, before finite-size effects occur.Here we observe the formation of oscillating bound states and show that they are robust against variations in the coupling parameters.
Oscillating bound states are a fascinating phenomenon of giant atoms: When an initially excited giant atom decays into a waveguide then, under certain resonance conditions, a significant part of the energy may end up oscillating back and forth between the atom and field.The first derivation of this phenomenon in [8], using RWA and δ-coupling points for the atom, identified specific combinations of parameter values for the number of coupling points, the coupling strength and the atom's frequency, at which oscillating bound states appear.In particular, in view of future experimental studies, this raises the questions of whether oscillating bound states can also be expected for finite-width coupling points, whether the appear-  ance of oscillating bound states is robust against deviations in the coupling and frequency parameters, and whether oscillating bound states also appear in the strong coupling regimes.
In the following, we are able to answer these questions in the affirmative.Fig. 6 demonstrates the appearance of an oscillating bound state for the giant atom setup as introduced in Tab.I, whose parameters were chosen to correspond closely to an oscillating bound state configuration of [8].The giant atom is initially in the excited state |e when it is coupled to the waveguide in the vacuum at time t = 0, i.e., it starts out with a population of p e = 1.After the initial decay process, which takes of the order of approximately 5τ to 10τ , the system can realize an oscillating bound state, for certain parameters.These states are characterized by a steady oscillation in the atomic population.Because our setup uses a periodic waveguide, we can only meaningfully describe the atom's dynamics up to t 18τ .After this time, radiation emitted at t = 0 has traversed the waveguide and reaches back to the atom from the other direction [43].
The plots of Fig. 6 suggest that the appearance of oscillating bound states, to a certain extent, is robust against variations both in the atom frequency and the coupling strength.In view of the fact that our approach accounts for a non-zero, realistic width of the coupling points, this observation appears encouraging with respect to experimental implementations.Fig. 6(a), where the atom frequency Ω is varied while the coupling strength is fixed at λ = 0.4, shows regions with oscillating bound states appearing roughly periodically.Fig. 6(b), where the atom frequency is fixed at Ω = 160π/L while the coupling strength is varied, shows oscillating bound states only in the range of 0.3 λ 0.4.Within the RWA [8], one expects oscillatory bound states to appear periodically both in Ω and in λ.However, based on the analysis in the previous section, we would count all data in Fig. 6(a), and all data in Fig. 6(b) with non-trivial dynamics, towards the USC regime, and thus beyond the regime where the RWA is valid.Fig. 7 demonstrates that the parameter regime we consider requires numerically exact calculations, by comparing our MPS results to calculations obtained within the RWA.It also illustrates that the error introduced by the RWA can change unexpectedly, probably due to the multimode couplings of our approach.At the default values of Ωτ /(2π) = 4.0 and λ = 0.4, the agreement between the RWA and MPS results may still appear acceptable.Estimating the dimensionless emitter-waveguide coupling by λ √ µ 0 /Ω, one may thus assume that the agreement between the MPS and RWA calculations will improve when λ is decreased or Ω is increased.As far as the coupling strength is concerned, this is what we observe.However, in the atom frequency, the agreement of the RWA and MPS results is highly non-monotonous.Whereas there is an overall trend for the agreement to improve as Ω is increased, significant oscillations in the quality of the approximation can be observed.In Fig. 7, this is illustrated by the inset for the population dynamics at Ωτ /(2π) = 4.5.
Here, in contrast to the reasonably good agreement between both curves in the presence of the oscillating bound state at Ωτ /(2π) = 4, the RWA results in a significantly different prediction for the atom population.The results of this section show that our approach allows one to explore the dynamics of the system far into the USC regime.In fact, our numerical results indicate that for the setup considered here, intermediate-time evolutions are feasible up to a coupling strength of λ ≈ 1.5, at which point (i) the MPS simulations become too costly and (ii) the system enters the DSC regime, as outlined in Sec.III.Thus, since other, perturbative approaches are more suitable at DSC, we expect our method to be most useful in the intermediate USC regime.Further improvements of our numerical approach can be made by careful choice of parameters, e.g., cavity length L. Increasing L in order to extend the maximal simulation times, or to decrease the free spectral range of the cavity, would result in a larger number of modes in the field below the UV cutoff that need to be taken into account.However, since the resulting chain length scales linearly in the cavity length, the increase in computational costs (cf.[43]) may well be feasible.This could prove useful for further investigations of systems, e.g., motivated by concrete experimental setups.

V. CONCLUSIONS & OUTLOOK
In summary, we have investigated the low-energy sector and time dynamics of a giant atom coupled to a waveguide, beyond the rotating-wave approximation.We have outlined in detail how a prototypical model describing a giant atom coupled to all non-interacting field modes below a physically well-motivated UV cutoff can be conveniently cast into a form which is amenable to efficient numerical treatment using matrix-product states.This approach has enabled us to compute the low-energy spectrum of the system at highly elevated coupling strengths, i.e., going beyond the single-excitation subspace and identifying the onset of different strong-coupling regimes.Based on previous findings in the context of the thoroughly explored standard quantum Rabi model, we have identified these regimes as ultra-strong and deep-strong coupling limits.In contrast to earlier work, we have described the coupling between emitter and waveguide not as a point-like coupling, but using a profile function with a finite spread, suppressing the coupling to high-frequency field modes and allowing for a UV cutoff.Since the presented approach is numerically exact and provides full access to a variety of system and bath observables, we were able to analyze how the contributions to the ground and first-excited state energies are distributed among the system, field and interaction Hamiltonians.Using our numerical toolbox, we have calculated the low-lying eigenstates with up to ≈ 10 excitations in the entire system, including bath and emitter.Based on the relatively low computational costs of these simulations, our study paves the way for further numerical investigations of waveguide quantum electrodynamics with multiple giant atoms in all coupling regimes.Furthermore, we have studied the time evolution of the composite system in an acoustodynamical setting, which may be realized by coupling a superconducting qubit at several locations to an acoustic waveguide.It has previously been suggested that such setups, when operated in the non-Markovian limit, can host bound states characterized by a persistent exchange of energy between the artificial atom and its environment.Here we have explicitly taken into account the significant time delay caused by a finite propagation speed of the acoustic modes, to investigate the pronounced non-Markovian features that arise as a consequence.In contrast to earlier works, we have not been restricted to the single-excitation subspace, and demonstrated the emergence and robustness of oscillating bound states over a wide parameter range.In particular, the breakdown of the rotating-wave approximation can be carefully monitored by applying our ansatz to the models with and without counterrotating terms, respectively.Beyond the scope of this work, which focuses on single giant atoms, the chain transformation approach opens up the opportunity to also non-perturbatively study systems composed of two or several giant atoms coupled to a common environment, and within the ultra-strong coupling regime.In fact, already within the rotating-wave approximation, it can offer advantages since simulations of chain transformed systems, based on matrix-product states, can treat many numbers of excitations in the system without any adjustments, whereas the Hilbert space dimension of direct diagonalization scales unfavorably.
The investigation of systems with several atoms and many excitations is motivated by intriguing phenomena that already arise within the single-excitation subspace and the rotatingwave approximation.App.A and Fig. 8 present two examples of this: The former derives the formation of an oscillating bound state between two giant atoms.The latter shows the emission of radiation from two giant atoms initially prepared as Bell states.Depending on the relative phase of the Bell state, either all energy is radiated away into the waveguide or part of it remains bound in the field between the two atoms.Future research directions include the investigation of superradiance, chiral quantum acoustics with and without an intrinsic time delay, and explicitly time-dependent models (see also [44]), to implement gates between giant atoms.
At sufficiently strong couplings, it also becomes important to go beyond the two-level approximation and consider higher-lying excited states of the emitter.Numerical simulations based on the matrix-product state ansatz can treat fewlevel emitters, and thus the techniques employed in this work can readily be adapted for future studies of non-Markovian dynamics beyond the two-level approximation at strong couplings.An oscillating bound state in the two-atom setup can be found in cases where the dark-state conditions in Table II are fulfilled for various n.In Fig. 9, we show the resulting dynamics in different parameter regimes, but all in the non-Markovian regime where γτ > 1.While Fig. 9(a) displays a fast decay of the initial excitation, Figs.9(b) and Fig. 9(c) show the emergence of dark states.In the long-time limit, these do not decay despite their dissipative environment.The setup corresponding to Fig. 9(b) hosts a symmetric dark state for n = 2 (compare Table II).In Fig. 9(c), one symmetric (n = 10) and two anti-symmetric (n = 9, 11) dark states are present.More explicitly, the long-time limit of the initially excited atom is given by p  The version of the exact simple Lanczos algorithm which is most stable in numerical implementations is: Lanczos algorithm: v 0 = 0; β 0 = 0; v 1 = v/ v ; for j = 1 to n : w = Av j − β j−1 v j−1 ; α j = w † v j ; r j = w − α j v j ; β j = r j ; if β j = 0 : end; v j+1 = r j /β j ; end.
In practical implementations, the break condition can be replaced by β j < for a sufficiently small bound.
In finite precision arithmetic, rounding errors occur in (C2) which can be represented by an error vector, Av j = β j−1 v j−1 + α j v j + β j v j+1 + f j . (C3) Thus, defining ξ k,j = v † k v j as a symbol for the inner products of the iteratively obtained vectors, these no longer fulfill the ideal Kronecker relation ξ k,j = δ k,j .A key point is now that for the Lanczos algorithm to remain stable, it is not necessary to reorthogonalize all vectors, but it is sufficient to keep the v j semi-orthogonal, i.e., max 1≤k≤j−1 |ξ k,j | ≤ √ , for the roundoff unit .Hence, reorthogonalization is only required when this bound is violated at any iteration step of the algorithm.
The growth of the ξ k,j elements is determined by the recurrence relations [48] β j ξ k,j+1 = β k ξ j,k+1 + α k ξ j,k − α j ξ k,j (C4) together with ξ j,j = 1 and ξ k,k−1 = v † k v k−1 .These, however, cannot be exactly calculated in numerical implementations since the error vectors f k are not known.Instead, the idea of partial reorthogonalization is to give an estimate for the terms θ k,j ≡ v † j f k − v † k f j and ξ j,j+1 by simulating them with random numbers, ξ j,j+1 = n β 1 β j Ψ, Ψ ∈ N (0, 0.6), (C5) where N (0, χ) is a zero mean normal distribution with variance χ.These estimates are then used in the original version of the algorithm to determine which vectors, if any, should be reorthogonalized at any given step of the algorithm [48].After a reorthogonalization has occurred, the relevant ξ k,j elements are reset to a normal distribution, ξ k,j+1 = Ξ, Ξ ∈ N (0, 1.5). (C7) For our purpose, we found the following simplified version to be sufficient, applying full orthogonalization to all vectors (we also used wider normal distributions, as in [51]): Lanczos algorithm with partial orthogonalisation: v 0 = 0; β 0 = 0; v 1 = v/ v ; for j = 1 to n : w = Av j ; α j = v † j w; r j = w − α j v j − β j−1 v j−1 ; β j = r j ; Compute ξ k,j+1 for k = 1, . . ., j − 1 using Eq.(C4); Set ξ j,j+1 using Eq.(C5); Set ξ j+1,j+1 = 1; if max 1≤k≤j (|ξ k,j+1 |) ≥ √ : Orthogonalise r j against v 1 , . . ., v j ; Perform orthogonalisation in the next iteration; Reset ξ k,j+1 using Eq.(C7); Recalculate β j = r j ; if β j = 0 : end; v j+1 = r j /β j ; end. (C8) For setups where b different emitters couple to the field, block Lanczos algorithms can be used to transform the field Hamiltonian.The block Lanczos procedure takes a Hermitian matrix A ∈ C n×n and an orthonormal set of complex vectors Q 1 = (v 1 , . . ., v b ) as inputs.The algorithm then iteratively computes a unitary basis Q = (Q 1 , . . ., Q p ) and a block tridiagonal matrix T such that Lanczos algorithm, we get the following procedure: Block Lanczos algorithm: p = n/b; Q 0 , B 0 = 0; for j = 1 to p : Y = AQ j ; Also the block Lanczos algorithm needs to be stabilized in numerical implementations (e.g., see [52]).
8 A O 3 x G U n b J 0 u 2 r V Z d v 4 i w h A l F C b I o V H A 0 s 8 g U f P b R Y g A F b E 9 z I g Z e q m L a x x h i b k T s j Q Z k u i Y / 4 S n 3 R Y t e G 4 0 a 5 e t e E v G b Z g p 8 I T 7 nV O M y W 5 u 1 f R r 2 h P u Q 4 c l V 9 4 w c 8 p N h V P a m I q L T v E D c Y s R G d d l 5 i 3 z r J b r M 5 u u L P b x 0 n W T s r 7 K I U 2 f 6 r / O G 0 Y M s b G L C L x 1 z I Q a s T s f 8 A U K 2 j 4 r a F 7 5 T E G 4 j o e 0 0 l n t V I p W U V L P 0 D a v z 3 o 4 5 v D i U C 8 7 W 2 t + u O G v f 1 z v 9 l 6 1 A 1 / A C h 5 j l V N 9 g R 7 e Y 5 N 1 K H z F M X 7 h t / f F + + Z 9 9 3 7 8 o 3 p z b c 4 j n F v e z 1 P u q 6 H 4 < / l a t e x i t > 0 ĉ † 0 ĉ1 + h.c.< l a t e x i t s h a 1 _ b a s e 6 4 = " / B T KI J w B 2 1 4 M U n Z m 7 d c 0 I o k 9 j c U = " > A A A C + 3 i c h V F N T 9 t A E H 2 Y U j 5 b A h y 5 r B o h c U C R U y H g i I C i X i p R q Q E k D N H a 2 T i b O L a 1 3 i A B 5 W / 0 D / R W 9 d p b r / Q / t L + l h z 5 v D R I g x F r r e f t m 5 u 3 M T p g n u r C + / 3 v M G 3 8 x 8 X J y a n p m d u 7 V 6 / n a w u J h k Y 1 M p F p R l m T m O J S F S n S q W l b b R B 3 n R s l h m K i j c L B b + o / O l S l 0 l n 6 y F 7 k 6 H c o 4 1 V 0 d S U u q X f O D R H W t + D x o 9 0 V g d N w j F s G a C H r S C t n u n w U d G c f K 3 B H t W t 1 v + G 6 J x 6 B Z g T q q d Z D V / i B A B x k i j D C E Q g p L n E C i4 H e C J n z k 5 E 5 x R c 4 Q a e d X u M Y M c 0 e M U o y Q Z A f 8 x z y d V G z K c 6 l Z u O y I t y T c h p k C K 9 z 7 T j F k d H m r I i 5 o / 3 J f O i 5 + 8 o Y r p 1 x W e E E b U n H a K X 4 g b 9 F j x H O Z w y r y t p b n M 8 u u L L r Y c t 1 o 1 p c 7 p u w z u t P Z o 8 e Q G z i P w D s X G V M j d O d z v k B K 2 2 I F 5 S v f K g j X c Y d W O q u c S l o p S u o Z 2 v L 1 W Q / H 3 H w 4 1 M f g 8 G 2 j u d F Y / 7 h e 3 9 6 p B j 6 F Z b z B K q e 6 i W 2 8

5 x a r t a 2 9 o 7 O
r u 6 e 3 r 3 9 g s D A 0 v J P G D e X J s h e H s d p z R S r D I J J l H e h Q 7 i V K i r o b y l 2 3 t p z F d 0

5 x a r t a 2 9 o 7 O
r u 6 e 3 r 3 9 g s D A 0 v J P G D e X J s h e H s d p z R S r D I J J l H e h Q 7 i V K i r o b y l 2 3 t p z F d 0

5 x a r t a 2 9 o 7 O
r u 6 e 3 r 3 9 g s D A 0 v J P G D e X J s h e H s d p z R S r D I J J l H e h Q 7 i V K i r o b y l 2 3 t p z F d 0

5 x a r t a 2 9 o 7 O
r u 6 e 3 r 3 9 g s D A 0 v J P G D e X J s h e H s d p z R S r D I J J l H e h Q 7 i V K i r o b y l 2 3 t p z F d 0

5 x a r t a 2 9 o 7 O
r u 6 e 3 r 3 9 g s D A 0 v J P G D e X J s h e H s d p z R S r D I J J l H e h Q 7 i V K i r o b y l 2 3 t p z F d 0

5 2 m 7 s 9 F
e e 7 v W 6 n b r g a / g H u 7 j E a e 6 i S 5 e o c c 6 Y n z G D / z E r + B T 8 C X 4 G n w 7 o w Z L d c 5 d / L e C 7 3 8 B x U + i c w = = < / l a t e x i t > Periodic waveguide length L < l a t e x i t s h a 1 _ b a s e 6 4 = "

5 x a r t a 2 9 o 7 O
r u 6 e 3 r 3 9 g s D A 0 v J P G D e X J s h e H s d p z R S r D I J J l H e h Q 7 i V K i r o b y l 2 3 t p z F d 0

5 x a r t a 2 9 o 7 O
r u 6 e 3 r 3 9 g s D A 0 v J P G D e X J s h e H s d p z R S r D 3 T 2 9 f / 8 D g 0 P D I 6 N h 4 a W L y I E s 6 y p c 1 P 4 k S d e S J T E Z h L G s 6 1 J E 8 S p U U b S + S h 1 5 r P Y 8 f n k u V h U m 8 r y 9 S e d I W Q R y e h b 7 Q h

d 4 0 6 P
K T E b j A l p u k p G 5 e Y Y Y W 5 I 7 I k G Y L o k P 8 + T 9 c F G v O c a 2 Y 2 O + A t I b d m p o t t 7 l O r 6 J O d 3 y r p Z 7 S / u D 9 b r P

2 < l a t e x i t s h a 1 _ b a s e 6 4 =
c b N C e X o w = = < / l a t e x i t > Ĥint," 3 x b c S p n a 3 + V 6 c / E F q 6 U b 4 I z e QO I = " > A A A C 4 X i c h V F N a 9 t A E H 1 W v 5 y 0 a Z 0 G e s l F x B R 6 C E Y K p s 0 x t E 3 I p e B A b Q f s Y F b y x l 6 s L 1 Z r Q + r 6 B / R W e s 0 t 1 / Q P t b + l h z x t l E J j i l e s Z v b N m 7 c z O 0 E W q d x 4 3 q + K 8 + D h o 8 d P q m v r T 5 9 t P H 9 R 2 3 z Z y d O p D m U 7 T K N U n w Y i l 5 F K Z N s o E 8 n T T E s R B 5 H s B p M P R b w 7 k z p X a f L Z X G T y L B a j R J 2 r U B h C g 9 q r / l g Y 9 3 j Q j 4 U Z 6 3 i u E r P r + o t B r e 4 1 P L v c Z c c v n T r K 1 U p r v 9H H E C l C T B F D I o G h H 0 E g 5 9 e D D w 8 Z s T P M i W l 6 y s Y l F l h n 7 p Q s S Y Y g O u F / x F O v R B O e C 8 3 c Z o e 8 J e L W z H T x m v v I K g Z k F 7 d K + j n t H + 4 v F h v 9 9 4 a 5 V S 4 q 4 0 m l N u o H 8 A P s q m 6 7 6 5 Z + C n x L F z 2 e u p U g Q h l r f O + c e + 6 Z e + e G R a J K 4 3 m / l 5 w H y w 8 f P V 5 Z X V t / 8 v T Z 8 8 b G i 7 0 y n + p I 9 q I 8 y

FIG. 3 .
FIG.3.Coefficients appearing in the chain form of the field Hamiltonian(15) for a giant atom with parameters as in Tab.I.
= 160π/L TABLE I. Giant atom geometry and parameters used as default in figures and numerical results, unless stated otherwise.The (periodic) waveguide's length L sets the overall length scale.

FIG. 4 .
FIG. 4. Lowest eigenenergies of Ĥ as function of coupling λ.(a) Ground-state energy Ĥ (red, triangle) and contributions from atomic Hamiltonian ĤA (orange, triangle), field Hamiltonian Ĥf (blue, circle) and interaction Hamiltonian Ĥint (green, triangle).(b) Gap ∆H between ground and first-excited state energies (red, triangle), and decomposition into contributions as in (a).The labels are the same in (a) and (b).

FIG. 5 .
FIG. 5. Atomic population and occupation of field modes in the ground and first excited state, as a function of coupling strength λ.The upper row shows (a) the expectation values of the total photon number nfield , as well as (b) the number operator of the first field mode n1 = â † 1 â1, and of (c) the atomic occupation.Analogously, (d)-(f) show the increase of the expectation value in the first excited state, i.e., ∆n1 = ψES|n1|ψES − ψGS|n1|ψGS , etc.

FIG. 6 .
FIG. 6. Dynamic of the atomic population for variation (a) of the atom frequency Ω and (b) of the coupling strength λ, whereas the other parameter is kept at its default value of (a) λ = 0.4 and (b) Ωτ /(2π) = 4.The atom starts excited with pe = 1 at t = 0.For certain parameter regions oscillating bound states form.Note that the coupling point separation τ = L/20 is used as the length scale.(Step sizes used for plotting: time ∆t = 5 × 10 −4 , frequency ∆Ω = 0.04, coupling strength ∆λ = 0.0025.)

FIG. 7 .
FIG. 7. Visualization of the Ω−λ-parameter space explored in Fig. 6, together with a comparison of the exact MPS results to RWA calculations.The dotted black lines represent the parameters plotted in Figs.6(a) and (b), with their intersection point corresponding to the default parameters of Tab.I.The insets compare results from Fig. 6 (solid lines) with results obtained using RWA (dashed, red lines).