# Waveguide QED in the Dipole Gauge

Sergi Terradas-Briansó  Instituto de Nanociencia y Materiales de Aragón
(INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain  Luis Martín-
Moreno  Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-
Universidad de Zaragoza, Zaragoza 50009, Spain  David Zueco  Instituto de
Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza,
Zaragoza 50009, Spain

###### Abstract

In recent studies on ultrastrong coupling between matter and light in
cavities, the significance of gauge choice when employing the widely-used two-
level approximation has been highlighted. Expanding upon these investigations,
we extend the analysis to waveguide QED, where we demonstrate that truncations
performed in the dipole gauge also yield accurate results. To illustrate this
point, we consider the case of a dipole coupled to a cavity array. Various
numerical and analytical techniques have been employed to investigate the low-
energy dynamics of the system. Leveraging these theoretical tools, we argue
that single photon scattering is an ideal method for investigating gauge-
related issues. Our findings reveal two novel effects in the scattering
spectra, which cannot be reproduced in a truncated model using the Coulomb
gauge. Firstly, the primary resonance is modified due to a Lamb shift
contribution. Secondly, we observe asymmetric transmission amplitudes
surrounding this resonance, reflecting the asymmetry of the spectral density
in this model. Additionally, we explore other features in the scattering
spectra resulting from ultrastrong couplings, such as the emergence of Fano
resonances and inelastic channels. Finally, we propose an experimental test of
our ideas in the context of circuit QED.

##  I Introduction

Photons usually interact weakly with matter, which has lead to the adoption of
several approximations in quantum optics. The most common ones include the
rotating-wave approximation, truncating the matter description to its lowest
energy levels, neglecting the $A^{2}$ term, and the Markovian approximation
for computing spontaneous emission, among others. However, numerous
experiments have demonstrated that discrete quantum emitters can couple to
light beyond the limitations of perturbative coupling. This breakthrough has
been achieved by coupling these emitters to cavities Niemczyk  _et al._
(2010); Forn-Díaz  _et al._ (2010); Yoshihara  _et al._ (2017) and waveguides
Forn-Díaz  _et al._ (2017); Puertas Martínez  _et al._ (2019); Léger  _et al._
(2019); Vrajitoarea  _et al._ .

It has been demonstrated that, in certain cases, the interaction energy can be
comparable to the energies of light and matter. As a consequence, many
approximations break down, indicating that the coupling strength has entered
the ultrastrong regime (USC). One of the primary noticeable effects of
entering the USC is the significant emergence of processes that go beyond the
interchange of a single photon and matter excitation.

As a result, the widely used rotating-wave approximation (RWA) for the
interaction loses its validity, leading to the renormalization of the bare
emitter parameters and the emergence of a nontrivial ground state. Several
interesting phenomena have been discussed in relation to the latter, including
the possibility of converting virtual photons into real ones through ground
state perturbations Ciuti  _et al._ (2005); Stassi  _et al._ (2013); He  _et
al._ ; Gheeraert  _et al._ (2018); Sánchez-Burillo  _et al._ (2019); Cirio
_et al._ (2017), the localization-delocalization transition Peropadre  _et
al._ (2013); Shi  _et al._ (2018), and the potential for performing nonlinear
optics at the limits of single and zero photons Sanchez-Burillo  _et al._
(2014); Sánchez-Burillo  _et al._ (2015); Stassi  _et al._ (2017); Kockum  _et
al._ (2017a); Chen  _et al._ (2017); Kockum  _et al._ (2017b). To gain
comprehensive insights into light-matter interactions in the USC regime, see
Frisk Kockum  _et al._ (2019); Forn-Díaz  _et al._ (2019).

More recently, it has been discovered that the USC regime introduces
additional complications for conventional approaches used to describe light-
matter systems. The commonly employed two-level approximation (TLA) for the
matter subsystem has been found to lack gauge invariance in certain
descriptions. In fact, truncating a momentum-like coupling operator has been
shown to cause significant inconsistencies between complete and truncated
models De Bernardis  _et al._ (2018); Di Stefano  _et al._ (2019); Ashida  _et
al._ (2022). Instead, employing position-based interactions yields more
reliable results when combined with the TLA. This research area has sparked
numerous studies focusing on different approaches to truncating the matter
level subsystems Savasta  _et al._ (2021), and also on the truncation of
photonic levels or the number of modes Taylor  _et al._ (2022); Sánchez Muñoz
_et al._ (2018). Correctly applying the two-level approximation has
significant implications for a wide range of system properties. It affects not
only the energy levels De Bernardis  _et al._ (2018), but also the ground
state Settineri  _et al._ (2021), emission spectra Mercurio  _et al._ (2022);
Salmon  _et al._ (2022), spectral density Salmon  _et al._ (2022), and various
other observables Settineri  _et al._ (2021); Salmon  _et al._ (2022).

In this work, we build upon previous approaches to ensure gauge invariance in
cavity quantum electrodynamics and apply them to the domain of waveguide QED.
We propose that, for a dipole coupled to a waveguide, employing the dipole
gauge is more suitable for matter truncation. To illustrate this, we focus on
the case of a single dipole coupled to a cavity array and conduct both
numerical and analytical calculations. Notably, we demonstrate that scattering
experiments are ideal for testing gauge-related issues. Specifically, we
reveal that the transmittance minima exhibits a red-shift as the coupling
strength increases, even within the lower range of the ultrastrong coupling
regime. This red-shift effect arises from the contribution of Lamb shifts,
which contrasts with the constant resonance observed when truncation is
performed in the Coulomb gauge. This key result highlights a qualitative
distinction between truncation approaches carried out in two different gauges.
Furthermore, our analysis reveals that including counterrotating terms enables
the occurrence of Fano resonances and inelastic scattering processes.

The rest of the manuscript is organized as follows. In Sec. II, we present a
general gauge-invariant description of a single particle interacting with the
electromagnetic field. We apply this description to the case of waveguide QED
and particularize it to a cavity array waveguide. Section III introduces the
framework used to address the scattering processes in these systems and
discusses how gauge invariance affects their description. The different
theoretical methods employed in the scattering computations are described in
Sec. IV. In section V, we present and analyse the numerical results for the
scattering spectra both within and beyond the rotating-wave approximation.
Sec. VI proposes an implementation of this setup in a circuit QED platform.
Finally, we conclude in section VII. The manuscript also includes four
appendices: Appendix A analyses gauge invariance in the weak coupling limit.
Appendix B provides the analytical computation of the self-energy in the RWA.
Appendix C discusses the convergence of Coulomb gauge and dipole gauge
simulations. Appendix D deals with the continuous limit of the dipole gauge
model.

##  II Gauge invariant formulation of waveguide QED

###  II.1 Preliminaries

In this section, we present a formalism that describes the interaction between
a single emitter and the electromagnetic field. A convenient approach to
introducing gauge-invariant light-matter systems can be found in Refs. Di
Stefano  _et al._ (2019); Savasta  _et al._ (2021); Dmytruk and Schiró (2021).
In the Coulomb gauge, the Hamiltonian is expressed as:

$H_{C}=H_{\rm ph}+U^{\dagger}H_{\rm m}U\;.$ (1)

Here, $H_{\rm ph}$ represents the quantized Hamiltonian for the
electromagnetic field, which will be specified below, and $H_{\rm m}$ refers
to the matter Hamiltonian. The unitary operator $U$, which represents a gauge
transformation itself, is given by

$U=e^{iq{\bf A(x)}\cdot{\bf x}/\hbar}\;,$ (2)

where ${\bf A(x)}$ is the vector potential and $q$ denotes the emitter charge.
The variable ${\bf x}$ corresponds to the position operator of the emitter. It
is worth noting that ${\bf A(x)}$ can explicitly depend on this position.

Given $H_{\rm m}={\bf p}^{2}/2m+V({\bf x})$, Eqs. (1) and (2) result in the
minimal coupling Hamiltonian

$H_{C}=\frac{1}{2m}\Big{(}{\bf p}-q{\bf A}\Big{)}^{2}+V({\bf x})\;.$ (3)

This way of expressing $H_{C}$ also offers a straightforward understanding of
the dipole gauge. In fact, the Hamiltonian in the dipole gauge can be
represented by an gauge transformation, equivalent to the Power-Zienau-Woolley
(PZW) one, as [Cf. Eq. (1)],

$H_{D}=H_{\rm m}+UH_{\rm ph}U^{\dagger}\;.$ (4)

Within this formulation, it is evident that $H_{D}=UH_{C}U^{\dagger}$,
ensuring gauge invariance. Furthermore, the light-matter coupling transforms
the _bare_ matter Hamiltonian in the Coulomb gauge, while in the dipole gauge,
it transforms the _free_ electromagnetic field Hamiltonian. Up to this point,
there should be no ambiguity about working in one gauge or the other. However,
in practical calculations, we often need to make approximations in the
Hamiltonians, especially when focusing on low-energy dynamics. In such cases,
we typically truncate the matter Hamiltonian to some minimum energy states,
using methods like the two-level approximation or the single band limit, among
others. The problem of gauge ambiguities arises when making these truncations,
since applying approximations on the matter subsystem in one gauge may lead to
different results than in another. This discrepancy can be understood by
examining Eqs. (1) and (4): matter and photonic operators are not the same in
$H_{C}$ as in $H_{D}$ [cf. Eq. (2)]. For instance, in single-mode cavity QED,
where ${\bf A(x)}={\bf A_{0}(x_{0})}(a+a^{\dagger})$, it has been argued that
the correct starting point for performing the truncation is $H_{D}$, while
truncation in $H_{C}$ yields incorrect results in the USC regime. This is
because in the dipole gauge, the coupling is ${\bf E(x)\cdot x}$ (${E}({\bf
x})$ is the electric field), while in the Coulomb gauge, it goes as ${\bf
A(x)\cdot p}$. As the emitter wave functions are typically localized in space,
the momentum-${\bf p}$ matrix elements can not be neglected even between well-
separated (in energy) eigenstates De Bernardis  _et al._ (2018). The same is
expected to occur in waveguide QED, as discussed below.

###  II.2 Waveguide QED Hamiltonian

EM quantization is convenient in the Coulomb gauge. In what follows, we will
assume that the _position_ of the emitters is fixed. This is the case in the
majority of situations and facilitates the discussion. To be more precise, we
will assume that the emitter position can be written as ${\bf x}={\bf
x_{0}}+\delta{\bf x}$ (${\bf x_{0}}$ is not longer an operator but a vector
position). For all the relevant energy scales ${\bf A(x)}\cong{\bf A(x_{0})}$,
which is the long-wavelength approximation. In this scenario, ${\bf A}$ acts
only on the Hilbert space of the photons (Román-Roche and Zueco, 2022, Sect.
3).

${\bf A}_{\perp}({\bf
x_{0}})=\frac{1}{\sqrt{L}}\sum_{k}\Big{(}{\boldsymbol{\lambda}}_{k}(x_{0},y_{0})\;a_{k}e^{ikz_{0}}+{\rm
h.c.}\Big{)}\;.$ (5)

We have added the suffix $\perp$ to emphasize that the potential vector has
only transverse components in the Coulomb gauge. In this work, we focus on
waveguide QED, which fixes a propagation direction, denoted as $z$ (cf. the
exponentials and the scalar character for the wavevector $k$). Additionally,
this constrains ${\boldsymbol{\lambda}}_{k}$ to the $xy$-plane, which can be
expressed as:

${\boldsymbol{\lambda}}_{k}\equiv\sqrt{\frac{\hbar}{2\omega_{\bf
k}\epsilon_{0}}}{\bf u}_{\bf k}\;,$ (6)

with ${\bf u}_{k}$ normalized functions, satisfying
$\int_{\mathbb{R}_{2}}dxdy\,|{\bf u}_{k}|^{2}=1$, and providing the space
dependence of ${\bf A}$ around the waveguide. The waveguide is assumed to be
surrounded by vacuum, hence the appearance of $\epsilon_{0}$ above. The
frequencies of the different waveguide modes are denoted as $\omega_{k}$:

$H_{\rm ph}=\sum_{k}\hbar\omega_{k}a_{k}^{\dagger}a_{k}\;.$ (7)

Up to this point, we have assumed a waveguide of length $L$, using a discrete
number of modes $k=m\times\pi/L$, where $m=0,\pm 1,...$. However, the
continuum limit can be obtained by taking $L\to\infty$.

At this point, it becomes evident that continuing to use the Coulomb gauge
comes with certain obstacles. The most evident is the fact that the minimal
coupling, introduces the ${\bf A}^{2}$-term which, after Eq. (5), couples all
the waveguide modes, resulting in terms like
$a_{k}a_{k^{\prime}}+a_{k}^{\dagger}a_{k^{\prime}}+{\rm h.c.}$. As a
consequence, the photonic part should be diagonalized using a Bogolioubov
transformation.

Alternatively, we can switch to the dipole gauge (4). To achieve this, we just
need to know how $a_{k}$ transforms under $U$ in (2),

$Ua_{k}U^{\dagger}=a_{k}-iq\frac{1}{\sqrt{L}}{\boldsymbol{\lambda}}_{k}e^{-ikz_{0}}\cdot{\bf
x}$ (8)

Therefore, the waveguide QED Hamiltonian in the dipole gauge can be written
as:

$H_{D}=H_{\rm
m}+\sum_{k}\hbar\omega_{k}a_{k}^{\dagger}a_{k}-iq\frac{1}{\sqrt{L}}{\bf
x}\sum_{k}\Big{(}{\boldsymbol{\lambda}}_{k}\hbar\omega_{k}e^{ikz_{0}}a_{k}-{\rm
h.c.}\Big{)}+\frac{q^{2}}{L}\sum_{k}\hbar\omega_{k}\Big{|}{\boldsymbol{\lambda}}_{k}e^{ikz_{0}}\cdot{\bf
x}\Big{|}^{2}.$ (9)

The last term in Eq. (9) can be absorbed into $H_{\rm m}$, resulting in a
modified matter Hamiltonian ${H}^{\prime}_{\rm m}=H_{\rm
m}+\frac{q^{2}}{L}\sum_{k}\hbar\omega_{k}\Big{|}{\boldsymbol{\lambda}}_{k}e^{ikz_{0}}\cdot{\bf
x}\Big{|}^{2}$. If ${H}^{\prime}_{\rm m}$ can be described using its two
lowest states, say $\\{|{0}^{\prime}\rangle,|{1}^{\prime}\rangle\\}$, $H_{D}$
can be expressed in the form of a spin-boson model:

${\mathcal{H}}_{D}=\frac{\hbar{\Delta}^{\prime}}{2}\sigma_{z}+\sum_{k}\hbar\omega_{k}a_{k}^{\dagger}a_{k}+\sigma_{x}\sum_{k}\Big{(}\hbar
g_{k}a_{k}+{\rm h.c.}\Big{)}$ (10)

where

$g_{k}=\frac{\omega_{k}}{\sqrt{L}}\langle{0}^{\prime}|{\bf
d}|{1}^{\prime}\rangle\cdot{\boldsymbol{\lambda}}_{k},$ (11)

with $\hbar{\Delta}^{\prime}$ being the transition energy between the two
states and ${\bf d}=q{\bf x}$. Here and throughout the text, we use the
notation $\mathcal{H}_{D}$ to distinguish the truncated Hamiltonian in Eq.
(10) from the full model $H_{D}$ in (9). It is customary to define the
spectral density for spin-boson models as:

$J_{D}(\omega)=2\pi\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})\;.$ (12)

Here, we have used the suffix $D$ to emphasize that the spectral density is
gauge-dependent since, it depends on the light-matter coupling (11).

As mentioned earlier, the ${\bf A}^{2}$-term in the Coulomb gauge prevents us
from explicitly expressing the Hamiltonian in the form of (9) and/or (10).
However, at low and intermediate couplings, the ${\bf A}^{2}$ is usually
neglected. By doing so, one can obtain a spin-boson model, but with a
different transition frequency $\Delta$ (the last term of (9) has been
absorbed into $H_{\rm m}$) and different a spectral density. Since the spin-
boson model is determined by its the spectral density, the Coulomb and dipole
gauges are not equivalent after truncation, similar to the case in cavity QED.
Despite this, it is possible to show that,

$\lim_{|{\boldsymbol{\lambda}}_{k}|\to 0}J_{C}(\Delta)\to J_{D}(\Delta)$ (13)

The proof is provided in Appendix A.

Thus, truncation can be safely done in both gauges in the usual scenario of
waveguide QED, in which the light-matter coupling is _small_. However,
entering the USC regime requires more caution, and as we will discuss in
detail later, the dipole gauge proves to be quite convenient. In addition, we
will explore the physical consequences of working in the USC regime and how
they manifest in standard experiments, such as single photon scattering.

###  II.3 The Cavity array case

Figure 1: Schematic representation of a one-dimensional array of coupled
cavities via the hopping parameter $\xi$. The dipole lies within the centre
cavity with a relative distance between its charges $x$.

Up to this point, we have not specified the particular waveguide or emitter we
are considering. Moving forward, we will analyze a model that allows for both
exact numerical treatments and analytical estimates. To describe the emitter,
we will assume it to be spherically symmetric dipole, which enables us to use
a one-dimensional model with position $x$ and momentum $p$ operators to
represent its Hamiltonian. On the other hand, we will focus on a cavity array
for the waveguide, where the emitter is coupled to a single cavity
(specifically, the $n=0$ cavity). We will also assume that the cavities in the
array are single-mode, with the vector potential aligned with the dipole in
the Coulomb gauge. This setup ensures that the model remains fully one-
dimensional. A schematic representation of this system is shown in Fig. 1. As
we will discuss, these simplifications can be justified within the circuit QED
architecture, where the USC regime has been reached, and our proposed ideas
may be implemented. Then, the resulting _full_ light-matter Hamiltonian in the
Coulomb gauge can be expressed as:

$\displaystyle H_{C}$
$\displaystyle=\frac{\left[p-q{A}_{0}(a_{0}+a^{\dagger}_{0})\right]^{2}}{2m}+V(x)$
$\displaystyle+\hbar\omega_{c}\sum_{n}^{N}a^{\dagger}_{n}a_{n}+\hbar\xi\sum_{n}^{N}(a_{n}a^{\dagger}_{n+1}+a^{\dagger}_{n}a_{n+1})$
(14)

In this model, we consider identical $N$ cavities, each with a resonance
frequency $\omega_{c}$. The cavities are coupled to their left and right
neighbors with a strength $\xi$ in a tight-binding fashion (last term of
equation (II.3)). Regarding the dipole, we describe it using a double-well
potential:

$V(x)=-\mu\frac{x^{2}}{2}+\lambda\frac{x^{4}}{4}.$ (15)

Moving to the momenta space $a_{k}=1/\sqrt{N}\sum_{n}a_{n}e^{ikn}$ and
transforming to the dipole gauge following Eqs. (4), (8), we obtain the
waveguide Hamiltonian with $H_{\rm
m}^{\prime}=p^{2}/2m+V(x)+\hbar\omega_{c}q^{2}{A}^{2}_{0}x^{2}$ [Cf. Eqs. (9)
and (II.3)] and,

$\lambda_{k}=A_{0}$ (16)

Our simulations consider both the non-truncated (which we call the full model)
and the truncated version where only the two lowest states,
$\\{|0^{\prime}\rangle,|1^{\prime}\rangle\\}$ of $H_{\rm m}^{\prime}$ are
retained. In the latter case, the model to consider is the spin-boson one (10)
with,

$\displaystyle g_{k}=\frac{\omega_{c}qA_{0}|\langle
1^{\prime}|{x}|0^{\prime}\rangle|}{\sqrt{N}}\frac{\omega_{k}}{\omega_{c}},$
(17)

then from Eq. (12)

$\displaystyle
J_{D}(\omega)=\frac{2g^{2}}{\sqrt{4\xi^{2}-(\omega-\omega_{c})^{2}}}\frac{\omega^{2}}{\omega^{2}_{c}},$
(18)

to quantify the coupling regime of the system in the truncated dipole gauge,
we define the coupling strength $g=q{A}_{0}\omega_{c}|\langle
1^{\prime}|x|0^{\prime}\rangle|$.

Moreover, we find it useful to express the dipole gauge Hamiltonian in
position space, as we will employ this representation in our scattering
calculations. Using $a_{n}=1/\sqrt{N}\sum_{k}a_{k}e^{-ikn}$, it yields,

$\displaystyle H_{D}=$ $\displaystyle
H_{m}^{\prime}+\hbar\omega_{c}\sum_{n}a^{\dagger}_{n}a_{n}+\hbar\xi\sum_{n}(a^{\dagger}_{n}a_{n+1}+{\rm
h.c.})-i\hbar\omega_{c}q{A}_{0}(a^{\dagger}_{0}-a_{0})x-i\hbar\xi
q{A}_{0}\left[(a^{\dagger}_{1}-a_{1})+(a^{\dagger}_{-1}-a_{-1})\right]x.$ (19)

The gauge transformation intertwines light and matter; hence, in this new
gauge, the dipole also couples to the adjacent cavities $n=-1,1$, which, in
the Coulomb gauge, are the ones coupled to cavity $n=0$ via the hopping term.

###  II.4 Truncation of matter levels in both gauges

In general, Hamiltonians (9), (II.3) or even the truncated version (10) are
non-integrable. Besides, the numerical calculation of the spectrum is rather
challenging. Hence, performing a benchmark of the truncation in different
gauges with the full model is not viable. To illustrate how the truncation
affects gauge invariance, we present the eigenenergies of a minimial system
comprising a dipole coupled to three cavities. These are obtained in a gauge
invariant non-truncated representation and within the two-level approximation
in both the Coulomb and dipole gauges in Fig. 2.

Refs. De Bernardis  _et al._ (2018); Di Stefano  _et al._ (2019) introduce a
convenient dimensionless notation for the matter Hamiltonian, which we also
follow in the diagonalization procedure to obtain the spectra of Fig. 2. By
defining a length scale $x_{0}=(\hbar^{2}/(\lambda m))^{1/6}$, we work with
the dimensionless variable $z=x/x_{0}$, allowing us to rewrite the bare matter
Hamiltonian $H_{m}$ as:

$\displaystyle H_{m}=E_{d}\left[\frac{p^{2}_{z}}{2}-\frac{\beta
z^{2}}{2}+\frac{z^{4}}{4}\right],$ (20)

where $\beta=m\mu x^{4}_{0}/\hbar^{2}$, $p_{z}=-i\hbar\partial/\partial z$.

Fig. 2 demonstrates that, similar to the single cavity QED case De Bernardis
_et al._ (2018); Di Stefano  _et al._ (2019), the truncation fails in the
Coulomb gauge, while it accurately matches the full model in the dipole gauge.
This observation highlights that the truncation in the dipole gauge, with an
interaction described by the position operator, preserves the energy levels of
the Hamiltonian. In contrast, calculations within the Coulomb gauge are
profoundly affected by the truncation, which is consistent with the results
obtained the single cavity scenario De Bernardis  _et al._ (2018); Di Stefano
_et al._ (2019).

Figure 2:  The energy spectra of an array of three cavities with the central
one coupled to the dipole, are shown for the full model (solid black line),
the truncated dipole gauge (dashed red line) and the truncated Coulomb gauge
(dashed-dotted blue line). In this plot, we choose $\beta=3.8$ and
$E_{dip}/(\hbar\omega_{c})=63.812$ such that at zero coupling, there is a
resonance between the cavity energy $\omega_{c}$ and the first bare dipole
transition $\Delta$. The hopping parameter is given by
$\xi/\omega_{c}=-1/\pi$. For the full model diagonalization, we consider 18
dipole levels and 18 photonic levels in each of the three cavities, while the
truncated cases use the same number of photonic excitations and 2 dipole
levels.

##  III Scattering: setting the problem

A caveat in USC is that the distinction between “light” and “matter”
subsystems becomes less clear. This is not only due to the creation of
strongly correlated light-matter states but also because the definition of
these subsystems is gauge-dependent. While physical observables remain gauge-
invariant, the choice of matter or light observables can be ambiguous. To
overcome these subtleties and provide _clear and measurable signatures of
truncation issues in different gauges_ , this paper focuses on studying
scattering phenomena. Scattering serves as an ideal testbed for highlighting
gauge-related problems. The scattering problem can be simplified as follows.
The input state, representing our initial condition, is chosen to be the non-
normalized quantum state,

$|\Psi_{in}\rangle=(a_{\phi}^{\dagger})^{N}|\mathrm{GS}\rangle,\qquad
a_{\phi}^{\dagger}=\sum_{x}\phi^{\rm in}_{x}a_{x}^{\dagger}\,,$ (21)

where $|\mathrm{GS}\rangle$ is the ground state of the system and
$\phi_{x}^{\rm in}$ is a Gaussian wavepacket centred in $x_{\rm in}$ with
spatial width $\theta$,

$\phi_{x}^{\rm in}=\exp\left(-\frac{(x-x_{\rm in})^{2}}{2\theta^{2}}+ik_{\rm
in}x\right)\;.$ (22)

Typically, we consider $x_{\rm in}$ located on the left-hand side of the
scatterer, with the wavepacket moves to the right towards it. The wavepacket
(22) is exponentially localized around $k_{\rm in}$, with a width of
approximately $\theta^{-1}$.

The wave packet evolves in time as

$|\Psi(t)\rangle=U(t,0)|\Psi_{\rm in}\rangle={\rm e}^{-iH_{\rm
tot}t}|\Psi_{\rm in}\rangle\;.$ (23)

A final time $t_{\mathrm{out}}$ is chosen, which must be sufficiently large to
allow the photons to move freely along the waveguide after interacting with
the scatterers. The evolution is then described by the $S$-matrix, defined as:

$|\Psi_{\rm out}\rangle=S|\Psi_{\rm in}\rangle\,.$ (24)

The scattering matrix is characterized by its momentum components:

$S_{p_{1}...p_{N^{\prime}},\,k_{1}...k_{N}}=\langle\mathrm{GS}|a_{p_{1}}...a_{p_{N^{\prime}}}\,S\,a_{k_{1}}^{\dagger}...a_{k_{N}}^{\dagger}|\mathrm{GS}\rangle\;.$
(25)

Some comments are pertinent here. The ground state wavefunction
$|\mathrm{GS}\rangle$ appears in the definition of $S$. In the ultrastrong
coupling regime, as discussed before, the ground state differs from the vacuum
state and contains a non-zero number of excitations Sánchez-Burillo  _et al._
(2019). In this regime, the number of excitations is not conserved, so
$N^{\prime}\neq N$ in general in the Eq. (25). However, we can expect some
simplifications to occur. We specialized our discussion to single photon
wavepackets, searching for computational simplicity, and also because in the
experiments, when using low-power coherent classical input/output field
states, the transmittance amplitudes coincide with the single photon
scattering amplitudes. Moreover, it has been shown Sánchez-Burillo  _et al._
(2018) that for wavepackets far away from the scatterer, even in the
ultrastrong coupling regime, we can approximate
$a_{\phi}^{N}|\mathrm{GS}\rangle\cong a_{\phi}^{N}|{\rm vac}\rangle$, where
$|{\rm vac}\rangle$ represents the trivial vacuum of the waveguide with
$a_{n}|{\rm vac}\rangle=0$ for all $n$. Lastly, it has been numerically tested
that the probability of having more of one photon in the output field is
negligible Sanchez-Burillo  _et al._ (2014); Sánchez-Burillo  _et al._ (2015).
As a consequence, the single photon amplitudes can be related to the number of
photons as follows:

${\frac{\langle{\Psi_{\rm out}}|a_{k}^{\dagger}a_{k}|\Psi_{\rm
out}\rangle}{\langle{\Psi_{\rm in}}|a_{k}^{\dagger}a_{k}|\Psi_{\rm
in}\rangle}}=|S_{kk}|^{2}\equiv t_{k}.$ (26)

This equation defines the _transmission amplitude_ $t_{k}$, which is a key
quantity in this work.

Within this framework, it is easy to understand how scattering is free from
ambiguities in the abovementioned sense. Of course, being an observable, the
transmission amplitude is gauge invariant. What makes scattering “special” is
that both the input and output fields have support only on regions well
separated from the scatterer. Consequently, they can be considered as free
photon wavepackets created over the QED vacuum of the waveguide, which remains
the same in all the gauges: $U|\Psi_{\rm in}\rangle=|\Psi_{\rm in}\rangle$
(the same with $|\Psi_{\rm out}\rangle$). Therefore, one can start with the
same initial conditions in both gauges, allow the system to evolve, compute
the amplitude and compare the results after performing truncations in both the
dipole and Coulomb gauges.

##  IV Theoretical methods.

Let us outline the theoretical methods we employed to compute the scattering
spectra in the USC regime. Our approaches involve a Polaron-like
transformation, numerical simulations based on Matrix Product States, and
matching techniques. The readers already familiar with these methods or not
interested in technical details, should better proceed to the next section.

###  IV.1 Polaron picture: effective single excitation dynamics

The polaron formalism provides an effective description of the low-energy
sector of spin-boson models such as (10). This approach is based on applying a
unitary transform that disentangles the light and matter subsystems.

$U_{\mathrm{P}}=\exp\left(-\sigma_{x}\sum_{k}f^{*}_{k}a_{k}-f_{k}a^{\dagger}_{k}\right),$
(27)

where $f_{k}$ represents a set of variational parameters. These parameters are
found by minimising the ground state energy, which is given by the state with
zero excitations
${|\Psi^{\mathrm{GS}}_{\mathrm{P}}\rangle}={|0\rangle}\otimes{|0_{k}\rangle}$,
see _e.g._ Silbey and Harris (1983); Bera  _et al._ (2014); Díaz-Camacho  _et
al._ (2016); Sánchez-Burillo  _et al._ (2019).

After applying (27), the effective Hamiltonian in the polaron picture is given
by

$\displaystyle\mathcal{H}_{P}$
$\displaystyle=\sum_{k}\omega_{k}a^{\dagger}_{k}a_{k}+\frac{\Delta_{r}}{2}\sigma_{z}+2\Delta_{r}\sum_{k}(f_{k}a_{k}\sigma^{+}+\mathrm{h.c.})$
$\displaystyle-2\Delta_{r}\sigma_{z}\sum_{k,p}f_{k}f^{*}_{p}a^{\dagger}_{k}a_{p},$
(28)

where

$\Delta_{r}=\exp\left(-2\sum_{k}\lvert f_{k}\rvert^{2}\right),$ (29)

and

$f_{k}=\frac{g_{k}}{(\Delta_{r}+\omega_{k})}.$ (30)

Here and throughout the manuscript, we use the convention $\hbar=1$ for
simplicity. The variational parameters $f_{k}$ and $\Delta_{r}$ are obtained
self-consistently from Eqs. (29) and (30). The renormalized frequency
$\Delta_{r}$ is a well-known result from the spin-boson model. It tends to
vanish as the coupling increases, eventually leading to the localization-
delocalization quantum phase transition depending on $J(\omega)$ Legget  _et
al._ (1987). The main advantage of the effective Hamiltonian (IV.1) is that it
conserves the number of excitations, making the single-particle dynamics
relatively straightforward.

###  IV.2 Matrix product states

Tensor networks have proven to be effective tools for numerically simulating
light-matter systems in the ultrastrong coupling regime Peropadre  _et al._
(2013); Sánchez-Burillo  _et al._ (2015); Shi  _et al._ (2018). Specifically,
for our one-dimensional chain with nearest-neighbor interactions described by
(19), we utilize matrix product states in conjunction with a time evolution
block decimation (TEBD) algorithm García-Ripoll (2006) to simulate wavepacket
scattering, as discussed in Section III. By simulating the scattering process
for a sufficiently long duration $t_{out}$, we can obtain the outgoing
eigenstate (24), as explained in the previous section. Furthermore, with the
matrix product state representation of ${|\Psi_{\rm in}\rangle}$ and
${|\Psi_{\rm out}\rangle}$, the computation of the transmission amplitude
transmittance described by Eq. (26) can be efficiently implemented.

###  IV.3 Matching techniques

Figure 3: The coupled waveguide is divided into three regions (I), (II), and
(III), which are employed in computing the transmission via matching methods.
Regions (I) and (III) are chosen to have negligible ground-state photonic
population, even in the ultrastrong coupling regime.

Lastly, we also use a matching technique to obtain the transmittance and
reflectance amplitudes. The fundamental concept in this framework involves
dividing the coupled waveguide into three distinct regions, as depicted in
Figure 3. The central region (II) encompasses the dipole, the central cavity,
and a sufficient number of additional cavities to contain the virtual photonic
excitations resulting from the dressing of the dipole in the ultrastrong
coupling regime. Ground state photons exhibit exponential localization around
the emitter, with $\langle a_{n}^{\dagger}a_{n}\rangle_{\rm gs}\sim
e^{-\kappa_{\rm GS}|n|/2}$, where $\kappa_{\rm GS}$ denotes the localization
length (which has been calculated in Román-Roche  _et al._ (2020) but is
irrelevant for our purposes). Consequently, the cavity regions (I) and (III)
can be considered, to a good approximation, to not have photons in the ground
state. In other words, in regions (I) and (III), we assume that $\langle
a_{n}^{\dagger}a_{n}\rangle_{\rm GS}=0$.

Bringing everything together, in regions (I) and (II), we describe single-
photon transport using the following state

${|\Psi\rangle}=\sum_{n,\alpha}\phi_{n,\alpha}a^{\dagger}_{n}{|0_{ph},\alpha_{sc}\rangle}+\sum_{\alpha}f_{\alpha}{|0_{ph},\alpha_{sc}\rangle}\;.$
(31)

Here, ${|0_{ph},\alpha_{sc}\rangle}$ represents the state with zero
excitations in regions (I) and (III). The label $\alpha$ denotes the
eigenstates in the region (II), which can be computed numerically. For the
parameters considered in this work, we have found that up to 5 or 7 cavities
in region (II) are sufficient to obtain an accurate description of the
transmittance.

The quantity $\phi_{n,\alpha}$ denotes the amplitude of having a photon in the
n-th cavity while the scatterer remains in state $\alpha$:

$\phi_{n,\alpha}(k)=\begin{cases}e^{ikn}+r_{k,\alpha}e^{-ikn},&\mathrm{(I)}\\\
t_{k,\alpha}e^{ikn},&\mathrm{(III)}\end{cases}\;.$ (32)

By employing the ansatz (31), the time-independent Schrödinger equation can be
solved, resulting in the determination of the transmittance and reflectance
amplitudes.

##  V Results

The complete transmittance spectrum is quite intricate, as we will soon
discover. To better understand the gauge issues and the impact of truncation
in different gauges, we will initially focus on calculating the scattering
coefficients within the rotating-wave approximation. This approach simplifies
single photon scattering and allows for fully analytical solutions.
Subsequently, we will proceed to solve the full model and comprehensively
discuss the complete transmittance spectrum.

###  V.1 Transmission under the RWA approximation

In order to isolate the effects of the truncation, we will apply the following
approximations. We apply the two-level truncation and the RWA such that
$\mathcal{H}_{D}=\mathcal{H}_{0}+\mathcal{H}_{I}$, with

$\mathcal{H}_{0}=\Delta\sigma^{+}\sigma^{-}+\omega_{c}\sum_{n}a_{n}^{\dagger}a_{n}+\xi\sum_{n}(a^{\dagger}_{n}a_{n+1}+{\rm
h.c.})$ (33)

and

$\mathcal{H}_{I}=g(i\sigma^{-}a_{0}^{\dagger}+{\rm h.c.})+\frac{\xi
g}{\omega_{c}}\left(\sigma^{-}(ia_{1}^{\dagger}+ia_{-1}^{\dagger})+{\rm
h.c.}\right)$ (34)

$\sigma^{\pm}$ are the two level system ladder operators and
$g=\omega_{c}q{A}_{0}\lvert\langle{0^{\prime}|x}|{1^{\prime}}\rangle\rvert$.
We emphasize that the rotating-wave approximation is valid only in the weak
coupling regime, but here we extend our computations to larger couplings to
compare the RWA results directly with those obtained from the full model,
thereby revealing the role of the counter-rotating terms.
Furthermore, we neglect the terms of second order in the couplings,
specifically the $x^{2}$ term in the dipole Hamiltonian. As mentioned, this
term introduces a dependency on the coupling strength $g$ to the TLS
transition $\Delta^{\prime}$ in Eq. (10). However, as shown in Figure 4,
$\Delta^{\prime}$ increases with $g$ but we note that this change is
relatively small compared to the other effects discussed in this section.
Therefore, for the RWA calculations, we will not consider this effect to gain
a qualitative understanding.

Figure 4: Renormalization of the first dipole transition energy in the dipole
gauge $\Delta^{\prime}$ with increasing coupling, presented in units of the
bare energy. The parameters used are the same as in Fig. 2.

The rotating-wave approximation offers several advantages. First, the ground
state becomes trivial, with $|\mathrm{GS}\rangle=|0\rangle$, implying zero
excitations in both the waveguide and the dipole. Moreover, within the RWA, we
can work within the single excitation manifold since the number of excitations
is conserved, i.e., $[\mathcal{H}_{D},N]=0$, where
$N=\sum_{n}a^{\dagger}_{n}a_{n}+\sigma^{+}\sigma^{-}$. In this case, the
matching method introduced in Sec. IV.3 becomes exact using the general
_ansatz_ for a quantum state in the single excitation manifold.

${|\psi\rangle}=\sum_{n}\phi_{n}(k)a^{\dagger}_{n}{|0\rangle}+\phi_{q}\sigma^{+}{|0\rangle}\;.$
(35)

The eigenvalue problem, $\mathcal{H}_{D}{|\psi\rangle}=E{|\psi\rangle}$, leads
to a dipole excited state amplitude given by:

$\phi_{q}=\frac{g}{E-\Delta}\left[\phi_{0}(k)+\frac{\xi}{\omega_{c}}(\phi_{1}(k)+\phi_{-1}(k))\right].$
(36)

In order to find the scattering eigenstates, we use the _ansatz_ for the
photonic amplitude at sites $n\neq 0$ as follows Zhou  _et al._ (2008),

$\phi_{n}(k)=\begin{cases}e^{ikn}+r_{k}e^{-ikn},&n<0\\\
t_{k}e^{ikn},&n>0\end{cases}\;.$ (37)

Using the continuity of the wave function, $\phi_{0^{+}}(k)=\phi_{0^{-}}(k)$
the solution for the transmittance amplitude is:

$\displaystyle
t_{k}=\frac{\Delta-\omega_{k}-\frac{g^{2}}{\omega^{2}_{c}}(\omega_{c}+\omega_{k})}{\omega_{k}-\Delta-i\frac{g^{2}\omega^{2}_{k}/\omega^{2}_{c}}{\sqrt{(2\xi)^{2}-(\omega_{k}-\omega_{c})^{2}}}+\frac{g^{2}}{\omega^{2}_{c}}(\omega_{k}+\omega_{c})}.$
(38)

Figure 5: Transmittance spectra of a single photon within the rotating-wave
approximation in both the truncated Coulomb gauge (a) and the truncated dipole
gauge (b). The dipole gauge resonant frequency obtained in Eq. (40) is shown
as a dotted cyan line. The parameters defining the system are the same as in
Fig. 2.

An equivalent computation can be performed for the Coulomb gauge, assuming the
same approximations mentioned in II.2 to obtain a spin-boson model, this is,
neglecting the $A^{2}$ term and considering the RWA. By applying these
transformations to Hamiltonian (II.3), we obtain a model that also conserves
the number of excitations. Further details about the resulting model can be
found in Appendix A. With this derived number-conserving description of our
system in the Coulomb gauge, we can use the _ansatz_ (35) to obtain $t_{k}$.
It yields the following expression (compare with Eq. (38))

$t_{k}=\frac{\Delta-\omega_{k}}{\omega_{k}-\Delta-\frac{ig^{2}_{C}}{\sqrt{(2\xi)^{2}-(\omega_{k}-\omega_{c})^{2}}}},$
(39)

where $g_{C}=qA_{0}\langle{0\lvert p}|{1}\rangle/m$ is the coupling strength
in the Coulomb gauge, see App. A for more details.

Figure 5 compares both formulas (38) and (39), with the transmittance
$T=\lvert t_{k}\rvert^{2}$, as a function of the coupling strength $g$ in the
dipole gauge and the incoming photon frequency $\omega_{k}$. In the truncated
Coulomb gauge, it is a known result Zhou  _et al._ (2008); Sanchez-Burillo
_et al._ (2014) that the resonance frequency always occurs at the dipole
transition, and the width of the transmittance minima increases as $g^{2}$ as
plotted in Fig. 5 (a).

The equivalent computation in the truncated dipole gauge is presented in
Figure 5 (b). A notable feature in the dipole gauge is that the resonance
moves to lower frequencies as $g$ increases. The resonant frequency can be
obtained by imposing $t_{k}=0$ in Eq. (38), leading to the formula:

$\omega^{\mathrm{RWA}}_{\mathrm{res}}=\omega_{c}\frac{\Delta\omega_{c}-g^{2}}{\omega^{2}_{c}+g^{2}}\;.$
(40)

This dependence is represented by the dotted cyan line in Figure 5 (b). The
resonant frequency shift is a consequence of the coupling-dependent term in
the numerator of Eq. (38), which arises from the couplings to the adjacent
cavities in the real-space Hamiltonian (19). This frequency shift will be
confirmed with the full model in the next section. Another effect of the
correct truncation is the modification of the width of the resonance, as given
by the imaginary term in (38). In Figure 5 (b), we can observe a transmittance
imbalance at both sides of the red-shifted minima. For a given coupling value,
the transmittance is higher at frequencies below the resonance than at
frequencies above.

To gain a deeper understanding of the resonance shift, we employ the resolvent
operator method Cohen-Tannoudji  _et al._ (1998). This approach allows us to
identify the change in resonance as a Lamb shift. The self-energy of our spin-
boson model within the RWA can be expressed as:

$\Sigma(E)=\sum_{k}\frac{\lvert
g_{k}\rvert^{2}}{E-\omega_{k}}=\frac{g^{2}}{N\omega^{2}_{c}}\sum_{k}\frac{\omega^{2}_{k}}{E-\omega_{k}}\;.$
(41)

In the continuum limit, this summation can be written in terms of known
integrals Economou (2006); Lombardo  _et al._ (2014). Further details on this
computation can be found in Appendix B. After performing some manipulations,
the self-energy of the system can be expressed as follows:

$\Sigma(E)=-\frac{g^{2}}{\omega^{2}_{c}}(\omega_{c}+E)+i\frac{g^{2}}{\sqrt{(2\xi)^{2}-(E-\omega_{c})^{2}}}\frac{E^{2}}{\omega^{2}_{c}}\;.$
(42)

The real part of the self-energy represents the Lamb shift, causing the red-
shift proportional to the square of the normalized coupling strength
$g/\omega_{c}$, as given in Eq. (38). Similarly, the imaginary part of (42)
corresponds to half of the spontaneous emission rate, which can also be
obtained by evaluating the spectral density $J_{D}(\omega)$ at the dipole
transition ${\Delta}$, as shown in Eq. (18).

The shift in the resonance frequency (40) is one of the primary findings of
this study. Although it has been calculated under the rotating-wave
approximation, we will observe a similar effect in the full model.
Importantly, this shift is a measurable quantity that sheds light on the
issues associated with truncating the matter in the Coulomb gauge.

###  V.2 Beyond RWA

Having established that the truncation within the Coulomb gauge cannot account
for the Lamb shift and the resonant frequency dependence on the coupling, we
now turn our attention to the full model without applying the RWA. This allows
us to study the transmission in a wide range of parameters.

To ensure the accuracy of our results and interpretations, we performed
simulations using the matching technique by truncating to different numbers of
levels in the dipole gauge. Additionally, we conducted MPS simulations, as
presented in Appendix C, where it can be observed that truncating the dipole
gauge to two levels yields the correct results. Therefore, in the main text,
we will only describe this specific case. Without the RWA, the number of
excitations is no longer conserved. As a result, processes that involve
different numbers of photons between the input and output fields become
possible. However, after conducting thorough numerical investigations, we have
found that these processes have negligible magnitudes. Therefore, in this
study, we can focus on the single photon transmission without significant loss
of accuracy.

Figure 6:  (a) Elastic transmittance spectra for a single-photon propagating
in a coupled cavity array, including all the interaction terms of the
truncated Coulomb gauge, as a function of the dipole gauge coupling strength
and the incoming photon frequency. The red dashed line indicates the scatterer
transition relevant to the photon transport in this gauge. (b) Equivalent
spectra computed in the dipole gauge. The red and orange dashed lines depict
the transition energies of the scatterer in the dipole gauge that play a role
in the transmission process. Additionally, we plot the resonance predicted in
the rotating-wave approximation in a dotted cyan line, as given in equation
(40). The resonance predicted utilising polaron techniques, as described in
section V, is depicted with a dashed white line.

The complete elastic transmittance spectra obtained in the Coulomb and dipole
gauges are depicted in Figures 6 (a) and (b), respectively. In the Coulomb
gauge spectra shown in Figure 6 (a), the transmittance minima remains
constant, much like in the number-conserving transmittance spectra in Figure 5
(a). One notable feature arising from the counterrotating terms in the Coulomb
gauge is a Fano resonance near the center of the band for all coupling
strengths.

The corresponding transmittance spectra in the dipole gauge are plotted in
Figure 6 (b). Several features can be observed in this plot. First, let’s
discuss the resonance frequency shift, which already occurs within the RWA, as
detailed above. To understand this shift in the full model, we can employ the
polaron picture and the effective Hamiltonian $\mathcal{H}_{P}$ (IV.1). This
transformed Hamiltonian is number-conserving, enabling us to compute the self-
energy similarly to the previous section. In this case, the self-energy is
given by:

$\displaystyle\Sigma_{P}(E)=\sum_{k}\frac{4\Delta^{2}_{r}\lvert
f_{k}\rvert^{2}}{(E-\omega_{k}-2\Delta_{r}\sum_{l,p}f_{l}f_{p})}\,.$ (43)

The equation for the resonance is now given by [Cf. with eq. (40)]:

$\omega_{k}-\Delta_{r}+\Re\left(\Sigma_{P}(\omega_{k})\right)=0\,.$ (44)

Here, $\Delta_{r}$ represents the renormalized transition frequency (29). In
Figure 6 (b), the frequency shift is depicted as a white dashed line, matching
the numerical results, and for comparison, we also plot the RWA result (40)
with a cyan dotted line. The shift is smaller than in the RWA case, which is
due to the competition of two effects. On the one hand, we have the
renormalized dipole frequency $\Delta_{r}$ from (29). On the other hand, there
is $\Re\left(\Sigma_{P}(\omega_{k})\right)$, which tends to shift towards
higher frequencies.

The main result of this paper is confirming the resonance shift in the full
model, which provides a qualitatively measurable feature to test gauge issues
related to the truncation of dipole energies.

Furthermore, Fano resonances also appear in the spectra. The first resonance
occurs at intermediate frequencies and spans most of the range of $g$ values.
Another resonance appears at larger coupling strengths and for the higher
frequencies of the band. These Fano resonances can be explained by considering
the interaction between the flying photon and the dipole, which is allowed to
access subspaces with a larger number of excitations due to the inclusion of
counterrotating terms Sanchez-Burillo  _et al._ (2014); Vrajitoarea  _et al._
. As a result, the flying photon can resonate with eigenstates of the
scatterer (defined as region (II) in Sec. IV.3) having different numbers of
excitations. In Figure 7, we plot the eigenstates of region II as a function
of the coupling strength. Two localized eigenstates with odd parity,
corresponding to three and five excitations, are identified and denoted as
$E_{3}$ and $E_{5}$, respectively. We also plot the energy differences between
these states and the ground state, i.e., $E_{3}-E_{0}$ and $E_{5}-E_{0}$, in
Figure 6 (b) as dashed lines in red and orange, respectively. A similar
analysis can be done in the Coulomb gauge, where a three-excitation eigenstate
can be associated with the Fano resonance observed in that spectrum. The
agreement with the numerical results confirms our argument explaining the
presence of the Fano resonances.

Figure 7: Lowest eigenenergies of Region (II), which includes the dipole and
the 5 central cavities, as introduced in IV.3.

Lastly, the full model exhibits inelastic scattering. The presence of bound
states originating from the ultrastrong coupling allows for scattering
processes that leave the scatterer in an excited state. In Figure 8, we
present the inelastic transmittance spectra obtained from the dipole gauge
computation, which reaches a maximum value of $0.25$, a fundamental bound as
explained in Sanchez-Burillo  _et al._ (2014).

These inelastic processes correspond to Raman scattering, leaving the dressed
dipole in an excited bound state. To delineate the parameter range where
inelastic transmission occurs, we use the energy conservation condition:

$\omega^{in}_{k}+E_{0}=E_{n}+\omega^{out}_{k},$ (45)

where $\omega^{out}_{k}\in[\omega_{c}-2\xi,\omega_{c}+2\xi]$, and $E_{0}$ and
$E_{n}$ are the energies of the ground state and bound states for the dressed
dipole. The solid red line in Figure 8 (c) represents the minimum energy of
the incoming photon for which inelastic scattering is possible, given by Eq.
(46):

$\omega^{min}_{ine}=E_{2}-E_{0}+\omega_{c}-2\xi.$ (46)

Our numerical computations verify this condition, providing the separation
line beyond which inelastic transmission is possible.

Figure 8: Inelastic transmittance $(1-T-R)/2$ in the dipole gauge. The orange
dashed line indicates the same Fano resonance as in 6. The red solid line
indicates the minimum energy of the incoming photon required for inelastic
scattering to occur, as given by Eq. (46).

##  VI Implementation

Finally, we propose a circuit QED architecture that provides an experimental
platform to test the ideas developed in this paper. Our proposed circuit,
illustrated in Figure 9, consists of an array of $N$ LC circuits coupled
inductively in series. The central circuit contains a superconducting qubit
(such as a transmon) capacitively coupled to the LC circuit, as shown in the
same figure.

By considering the Kirchhoff equations of motion and selecting $\phi_{0}$ and
$\phi_{q}$ as variables, we can express the current through the capacitor
$C_{r}$ as $C_{r}(\ddot{\phi}_{0}-\ddot{\phi}_{q})$, which leads to a light-
matter coupling via momentum, analogous to the minimal coupling in Eq. (1).
The Hamiltonian for this circuit is given by Devoret (1997); García Ripoll
(2022),

$\displaystyle{H}_{Ch}=$
$\displaystyle\sum_{n}\left[\frac{Q^{2}_{n}}{2C_{r}}+\frac{\phi^{2}_{n}}{2L_{\Sigma}}+\frac{\phi_{n}\phi_{n-1}}{L_{c}}\right]$
$\displaystyle+\frac{(Q_{0}+Q_{q})^{2}}{2C_{J}}+E_{J}\cos(\phi_{q})\;,$ (47)

where $1/2L_{\Sigma}=1/2L_{r}+1/L_{c}$. The suffix $Ch$ indicates that this
Hamiltonian is obtained in the charge gauge, resulting from the choice of
dynamical variables $\phi_{0}$ and $\phi_{q}$. The presence of the minimal
coupling-like feature in this circuit makes it analogous to the light-matter
Hamiltonian in the Coulomb gauge (II.3).

Figure 9: Cavity array implemented in a superconducting circuit QED platform.
Individual LC oscillators are coupled inductively to their neighbours, while
the transmon is directly connected to the capacitance of the central
oscillator, labelled $0$.

Similar to dipolar systems, a unitary transformation can also be applied in
this case (cf. Eq. (2)) De Bernardis  _et al._ (2018):

$U=e^{i\phi_{q}Q_{0}/\hbar}\;.$ (48)

Physically, this transformation changes the dynamical variables (for the
Kirchhoff equations) from $\\{\phi_{0},\phi_{q}\\}$ to
$\\{\delta\phi,\phi_{q}\\}$, where $\delta\phi=\phi_{0}-\phi_{q}$.

By utilizing the mode fluxes and charges,
$\phi_{k}=\frac{1}{\sqrt{N}}\sum_{n}\phi_{n}e^{ikn}$,
$Q_{k}=\frac{1}{\sqrt{N}}\sum_{n}Q_{n}e^{ikn}$, and
$\alpha_{k}=\left(\frac{1}{2L_{\Sigma}}+\frac{\cos(k)}{L_{c}}\right)$, along
with their quantization Devoret (1997); García Ripoll (2022):

$\displaystyle Q_{k}$
$\displaystyle=-i\sqrt{\frac{\hbar\omega_{k}C_{r}}{2}}(a^{\dagger}_{k}-a_{k}),$
$\displaystyle\phi_{k}$
$\displaystyle=\sqrt{\frac{\hbar\omega_{k}}{4\alpha_{k}}}(a^{\dagger}_{k}+a_{k})\;,$

the transformed Hamiltonian $H_{Fl}=UH_{Ch}U^{\dagger}$ can be written as:

$\displaystyle
H_{Fl}=\frac{Q^{2}_{q}}{2C_{J}}+E_{J}\cos\left(\frac{\phi_{q}-\phi_{ext}}{\Phi_{0}}\right)+\sum_{k}{\hbar\omega_{k}}a^{\dagger}_{k}a_{k}-\hbar\sum_{k}g_{k}(a_{k}+a^{\dagger}_{k})\phi_{q}+\frac{1}{N}\sum_{k}\alpha_{k}\lvert\phi_{q}\rvert^{2}\;.$
(49)

Here, the dispersion relation is $\omega_{k}=\omega_{r}+2\xi_{r}\cos(k)$,
where $\omega_{r}=(L_{\Sigma}C_{r})^{-1/2}$ and
$\xi_{r}=\omega_{r}L_{\Sigma}/L_{c}$. The coupling constants are given by:

$\hbar g_{k}=\sqrt{\frac{\hbar}{2L_{\Sigma}\omega_{c}N}}\omega_{k}.$ (50)

It is worth noting that $g_{k}\sim\omega_{k}$ as shown in Eq. (17). Moreover,
the bare matter Hamiltonian, defined by the first two terms in (49), is also
modified by a term scaling with the square of the flux operator. Therefore,
this circuit Hamiltonian is equivalent to (4).

##  VII Conclusions

In this work, we have extended the study of gauge issues in light-matter
coupled systems to the realm of waveguide QED. Gauge problems present
significant technical challenges when studying waveguide QED, particularly in
the ultrastrong coupling regime. We have employed various theoretical methods,
including numerical techniques (matching and MPS) and analytical approaches
(polaron transformation), to describe the system dynamics and compare
truncation in different gauges. We argued that scattering, a natural quantity
in waveguide QED, is ideal for testing different gauges, and it holds
relevance from an experimental perspective.

Our investigations have confirmed that the transmittance spectrum exhibits
both qualitative and quantitative differences when truncating in different
gauges. Numerical results have provided evidence that the dipole gauge is
well-suited for truncation, allowing for accurate transmittance spectra over a
wide parameter range. On the other hand, the Coulomb gauge is found to be
unsuitable for truncation. Figure 6 presents a clear visual representation of
these significant differences.

The main features of correct transmittance spectra include the coupling-
dependent resonant frequency shift, the emergence of Fano-like resonances, and
the occurrence of non-elastic scattering. These findings constitute important
aspects of our study. Furthermore, we have concluded the article by proposing
an experimental implementation of this physics using circuit QED.

##  VIII Acknowledgments

The authors acknowledge funding from the Spanish Government Grants
PID2020-115221GB-C41/AEI/10.13039/501100011033 and TED2021-131447B-C21 funded
by MCIN/AEI/10.13039/501100011033 and the EU “NextGenerationEU”/PRTR, the
Gobierno de Aragón (Grant E09-17R Q-MAD) and the CSIC Quantum Technologies
Platform PTI-001.

##  Appendix A Gauge Invariance at weak coupling

In this section, we demonstrate that in the weak coupling limit, both the
Coulomb gauge and dipole gauge formulations can be represented as spin-boson
models [see Eq. (10)]. Although their spectral densities are not identical,
under the Wigner-Weisskopf approximation, their spontaneous emission rates are
equal, as indicated in Eq. (13).

By substituting the expression of the vector potential in the Coulomb gauge,
as given in the main text equation (5), into the minimal coupling Hamiltonian
(3), we obtain a similar expression to that in the dipole gauge (9)

$H_{C}=H_{\rm
m}+\sum_{k}\omega_{k}a_{k}^{\dagger}a_{k}-\frac{q}{\sqrt{L}}\frac{{\bf
p}}{m}\sum_{k}\Big{(}{\boldsymbol{\lambda}}_{k}e^{ikz_{0}}a_{k}-{\rm
h.c.}\Big{)}+\frac{q^{2}}{2mL}\Big{(}\sum_{k}{\boldsymbol{\lambda}}_{k}e^{ikz_{0}}a_{k}+{\rm
h.c.}\Big{)}^{2}.$ (51)

To derive effective models in the weak coupling regime we consider the limit
$\boldsymbol{\lambda}_{k}\to 0$ in both (9) and (51) retaining only the first-
order terms in $\boldsymbol{\lambda}_{k}$. Therefore neglecting the shift in
$H_{m}$ for the dipole gauge and the ${\bf A}^{2}$-term in the Coulomb gauge.

After taking this limit, we project the effective Hamiltonians into the
subspace spanned by the first two energy levels of the bare matter Hamiltonian
$H_{m}$, denoted as ${{|0\rangle},{|1\rangle}}$. This procedure results in two
spin-boson models, one for each gauge, which can be written in a form similar
to the full model presented in the main text (10).

In this weak coupling limit, the waveguide modes $\omega_{k}$ and the two-
level system transition frequencies are equal for both models. The difference
between the weak coupling Coulomb and dipole spin-boson models lies in their
mode couplings. In the dipole gauge, the expression for the full model
couplings is given in (11). For the approximated model in the dipole gauge, we
have

$g_{{\bf x},k}=\frac{q\omega_{k}}{\sqrt{L}}\langle 0|{\bf
x}|1\rangle\cdot{\boldsymbol{\lambda}}_{k}.$ (52)

The corresponding couplings in the Coulomb gauge can be written as

$g_{{\bf p},k}=\frac{q}{m\sqrt{L}}\langle 0|{\bf
p}|1\rangle\cdot{\boldsymbol{\lambda}}_{k}$ (53)

We can now compute the spectral densities of these two models using the
definition given by Eq. (12) in the main text, obtaining

$\displaystyle J_{D}(\omega)$
$\displaystyle=\frac{q^{2}\omega^{2}}{L}\frac{\big{|}\langle 0|{\bf
x}|1\rangle\cdot{\boldsymbol{\lambda}_{k}}\big{|}^{2}}{\sqrt{4\xi^{2}-(\omega-\omega_{c})^{2}}},$
(54) $\displaystyle J_{C}(\omega)$
$\displaystyle=\frac{q^{2}}{m^{2}L}\frac{\big{|}\langle 0|{\bf
p}|1\rangle\cdot{\boldsymbol{\lambda}_{k}}\big{|}^{2}}{\sqrt{4\xi^{2}-(\omega-\omega_{c})^{2}}}.$
(55)

Utilizing the general relation

$\langle n|{\bf p}|l\rangle=im\Delta_{nl}\langle n|{\bf x}|l\rangle,$ (56)

where $\Delta_{nl}$ is the transition frequency between states ${|n\rangle}$
and ${|l\rangle}$, we can derive the following expression for the spectral
density in the Coulomb gauge

$J_{C}(\omega)=\frac{q^{2}\Delta^{2}}{L}\frac{\big{|}\langle 0|{\bf
x}|1\rangle\cdot{\boldsymbol{\lambda}_{k}}\big{|}^{2}}{\sqrt{4\xi^{2}-(\omega-\omega_{c})^{2}}}.$
(57)

This proves that in the weak coupling limit, the spontaneous emission rate
$J(\Delta)$ is gauge invariant

$\lim_{{\boldsymbol{\lambda}_{k}}\to
0}J_{C}(\Delta)=\lim_{{\boldsymbol{\lambda}_{k}}\to 0}J_{D}(\Delta).$ (58)

##  Appendix B Analytical computation of the Self-Energy of the coupled dipole

In this appendix, we provide the derivation of the exact self-energy of the
dipole excited level in the RWA (42). We begin the derivation from the
expression (41) and assume that the energy $E$ lies within the band
$E\in[\omega_{c}-2\xi,\omega_{c}+2\xi]$. However, the computation can be
easily extended to energies outside the band.

As mentioned in the main text, expanding the expression of the dispersion
relation in Eq. (41) yields:

$\displaystyle\Sigma(E)$
$\displaystyle=\frac{g^{2}}{N\omega^{2}_{c}}\sum_{k}\frac{\omega^{2}_{k}}{E-\omega_{k}}=$
$\displaystyle=\frac{g^{2}}{N\omega^{2}_{c}}\sum_{k}\Bigg{[}\frac{\omega^{2}_{c}}{E-\omega_{k}}+\frac{2\xi\omega_{c}(e^{ik}+e^{-ik})}{E-\omega_{k}}$
$\displaystyle+\frac{\xi^{2}\left(e^{ik}+e^{-ik}\right)^{2}}{E-\omega_{k}}\Bigg{]}.$
(59)

We can convert the summations over modes into integrals by taking the
continuum limit. ($N\to\infty$).

$\displaystyle\Sigma(E)$
$\displaystyle=\frac{g^{2}}{2\pi\omega^{2}_{c}}\int^{\pi}_{-\pi}\Bigg{[}\frac{(\omega^{2}_{c}+2\xi^{2})dk}{E-\omega_{k}}+\frac{2\xi\omega_{c}(e^{ik}+e^{-ik})dk}{E-\omega_{k}}$
$\displaystyle+\int^{\pi}_{-\pi}\frac{\xi^{2}(e^{2ik}+e^{-2ik})dk}{E-\omega_{k}}\Bigg{]}.$
(60)

The entire computation in (B) relies one valuating integrals of the form:

$\displaystyle I(E,n)=$
$\displaystyle\int^{\pi}_{-\pi}\frac{e^{ink}dk}{E-[\omega_{c}-\xi(e^{ik}+e^{-ik})]}.$
(61)

By performing the variable change $z=e^{ik}$, transforms into a closed
integral around a circle of radius $\lvert z\rvert=1$ in the complex plane

$I(E,n)=\frac{1}{\xi}\oint_{\lvert z\rvert=1}\frac{z^{n}dz}{z^{2}+2az+1}$ (62)

where we defined $a=(E-\omega_{c})/(2\xi)$. Since the integral is around a
closed circle, we have $I(E,n)=I(E,-n)$.

Integral (62) can be solved using the Residue Theorem. If the energy lies
within the band ($\lvert a\rvert<1$), the corresponding poles are
$z_{\pm}=-a\pm i\sqrt{1-a^{2}}$. These poles correspond to the limits
$\lim_{\eta\to 0}E\pm i\eta$.

$\lim_{\eta\to 0}I(E\pm i\eta,n)=\mp 2\pi i\frac{\left(-a\pm
i\sqrt{1-a^{2}}\right)^{\lvert n\rvert}}{\sqrt{4\xi^{2}-(E-\omega_{c})^{2}}}.$
(63)

Writing the self-energy (B) in terms of the integrals (62) gives:

$\displaystyle\Sigma(E)=i\frac{g^{2}}{2\pi\omega^{2}_{c}}\Big{[}(\omega^{2}_{c}+2\xi^{2})I(E,0)+$
$\displaystyle 4\xi\omega_{c}I(E,1)+2\xi^{2}I(E,2)\Big{]}.$ (64)

From where, after introducing the solution found in (63), we can obtain the
self-energy as in Eq. 42 of the main text. The real part of the self-energy
gives the Lamb shift of the resonant frequency, whereas the imaginary part
provides half its spectral width.

$\text{Re}(\Sigma(E))=-\frac{g^{2}}{\omega^{2}_{c}}(\omega_{c}+E).$ (65a)
$\text{Im}(\Sigma(E))=\frac{g^{2}}{\omega^{2}_{c}}\frac{E^{2}}{\sqrt{4\xi^{2}-(E-\omega_{c})^{2}}}=\frac{J(E)}{2}.$
(65b)

##  Appendix C Transmission in the Coulomb gauge

In this appendix, we demonstrate the effect of truncation in both gauges by
explicitly computing different transmission spectra while considering various
numbers of dipole levels $N_{d}$. As mentioned in the main text, applying the
two-level approximation in the dipole gauge can accurately approximate the
full model, whereas doing so in the Coulomb gauge does not yield the expected
results. However, we were able to perform some calculations in the Coulomb
gauge for $N_{d}>2$ at the intermediate range of the USC regime.

Figure 10: Transmittance spectra in the Coulomb gauge for a different number
of dipole levels $N_{d}$. (a) and (b) correspond to the spectra obtained using
the two-level approximation and for four levels, respectively. For a higher
number of dipole levels, such as in (c) $N_{d}=5$ and (d) $N_{d}=8$, we also
plot the RWA resonance as a dotted cyan line [Cf. Eq. (40)] and the USC
resonance line predicted with the polaron transformation in the dipole gauge.

In Fig. 10, we plot the transmittance spectra in the Coulomb gauge, including
all coupling terms, for different numbers of dipole levels. Fig. 10 (a) shows
the spectra obtained using the two-level approximation in the Coulomb gauge,
where the resonance is fixed at $\omega_{k}=\Delta$ and an additional Fano
resonance is observed for all couplings. In Fig. 10 (b), we increase the
number of dipole levels, resulting in a slight red-shift of the main resonance
and an increase in the transmittance at lower frequencies around the main
resonance while maintaining the Fano resonance. As the number of dipole levels
increases further in Fig. 10 (c) and (d), we observe a strong red-shift in the
transmittance minima and an imbalance of the transmittance on both sides of
the main resonance. These effects are consistent with what we found and
explained in the main text for the dipole gauge.

In both Fig. 10 (c) and (d), we also plot the minima of transmittance obtained
in the dipole gauge, as shown in the main text in Fig. 6. The blue dotted line
represents the analytical solution in the RWA, as given in Eq. (40), while the
white dashed line gives the solution in the effective polaron picture. It is
evident that capturing the main features obtained in the dipole gauge with the
two-level approximation in the Coulomb gauge requires a significantly higher
number of dipole levels. While the main resonance seems to be well described,
the Fano resonances do not coincide yet. The Fano resonances directly come
from the eigenenergies of the scatterer, which are well approximated in the
dipole gauge, as introduced in Sec. II.4.

However, obtaining an accurate two-level approximation in the Coulomb gauge
from the truncated dipole model is possible, as observed in the single cavity
case Di Stefano  _et al._ (2019).

Applying a truncated version of the Power-Zineau-Wooley transformation (2),
$\mathcal{U}=\exp(ig/\omega_{c}(a_{0}+a^{\dagger}_{0})\sigma_{x})$, we can
recover a Coulomb gauge description from
$\mathcal{H}_{C}=\mathcal{U}\mathcal{H}_{D}\mathcal{U}^{\dagger}$ [Cf. Eqs.
(33), (34)]

$\displaystyle\mathcal{H}_{C}$
$\displaystyle=\sum_{n}\omega_{c}a^{\dagger}_{n}a_{n}+\xi\sum_{n}(a^{\dagger}_{n+1}a_{n}+a^{\dagger}_{n}a_{n+1})$
(66)
$\displaystyle+\frac{\Delta^{\prime}}{2}\left[\cos\left(\frac{2g}{\omega_{c}}(a_{0}+a^{\dagger}_{0})\right)+\sin\left(\frac{2g}{\omega_{c}}(a_{0}+a^{\dagger}_{0})\right)\right].$

##  Appendix D Continuous limit of the system

Cavity array systems can be experimentally implemented in various platforms,
such as photonic crystals and superconducting systems Liu and Houck (2017);
Vrajitoarea  _et al._ , as discussed in Sec. VI. Moreover, these types of
models can also be seen as discretizations of general waveguide models.

In this appendix, we derive the continuous real-space description of our
system in the dipole gauge. This model differs from other standard continuous
waveguide QED models Shen and Fan (2009) due to the couplings to adjacent
cavities.

To begin, we split the momentum space into left and right propagating momenta
as follows:

$H_{ph}=\sum_{j=\\{L,R\\}}\int{dk_{j}}\omega({k_{j}})a^{\dagger}(k_{j})a(k_{j})$
(67)

where
$[a(k_{i}),a^{\dagger}(k^{\prime}_{j})]=\delta(k_{j}-k^{\prime}_{i})\delta_{ij}$.

Next, we introduce creation and annihilation operators for right or left
propagating photons at position $x$ as the Fourier transform of their momentum
counterparts

$\displaystyle a(k_{R}):=\int a_{R}(r)e^{-ik_{R}r}dr,$ (68) $\displaystyle
a(k_{L}):=\int a_{L}(r)e^{-ik_{L}r}dr.$ (69)

By linearizing the dispersion relation around a probe wavevector $k_{0}$,

$\displaystyle\omega_{L}(k)\simeq\omega_{L}(k_{0})-v_{g}(k_{0})k_{L},$ (70)
$\displaystyle\omega_{R}(k)\simeq\omega_{R}(k_{0})+v_{g}(k_{0})k_{R},$ (71)

with $k_{L}=k-k_{0}$ and $k_{R}=k+k_{0}$, we can find an expression for
Hamiltonian (67) in real space as

$\displaystyle H_{ph}$ $\displaystyle=\int
dra^{\dagger}_{R}(r)\left(\omega_{R}(k_{0})-iv_{g}\frac{\partial}{\partial
r}\right)a_{R}(r)$ $\displaystyle+\int
dxa^{\dagger}_{L}(r)\left(\omega_{L}(k_{0})-iv_{g}\frac{\partial}{\partial
r}\right)a_{L}(r).$ (72)

where we have used the relation $k_{R}e^{ik_{R}r}=-i\partial_{r}e^{ik_{R}r}$.

After applying the same linearization to the interaction term the full
Hamiltonian reads

$\displaystyle H_{D}$ $\displaystyle=\sum_{j}\int
dra^{\dagger}_{j}(r)\left(\omega_{j}(k_{0})-iv_{g}\frac{\partial}{\partial
r}\right)a_{j}(r)+H_{m}$ $\displaystyle+i{q{A}_{0}}x\sum_{j}\int
dk_{j}\left(\omega_{j}(k_{0})-iv_{g}\frac{\partial}{\partial
r}\right)a_{j}(r)$ $\displaystyle-i{q{A}_{0}}x\sum_{j}\int
dk_{j}\left(\omega_{j}(k_{0})-iv_{g}\frac{\partial}{\partial
r}\right)a^{\dagger}_{j}(r).$ (73)

## References

  * Niemczyk  _et al._ (2010) T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, A. Marx, and R. Gross, Nature Physics 6, 772 (2010).
  * Forn-Díaz  _et al._ (2010) P. Forn-Díaz, J. Lisenfeld, D. Marcos, J. J. García-Ripoll, E. Solano, C. J. Harmans, and J. E. Mooij, Phys. Rev. Lett. 105, 3 (2010).
  * Yoshihara  _et al._ (2017) F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Nat. Phys. 13, 44 (2017).
  * Forn-Díaz  _et al._ (2017) P. Forn-Díaz, J. J. García-Ripoll, B. Peropadre, J.-L. Orgiazzi, M. A. Yurtalan, R. Belyansky, C. M. Wilson, and A. Lupascu, Nat. Phys. 13, 39 (2017).
  * Puertas Martínez  _et al._ (2019) J. Puertas Martínez, S. Léger, N. Gheeraert, R. Dassonneville, L. Planat, F. Foroughi, Y. Krupko, O. Buisson, C. Naud, W. Hasch-Guichard, S. Florens, I. Snyman, and N. Roch, npj Quantum Inf. 5, 19 (2019).
  * Léger  _et al._ (2019) S. Léger, J. Puertas-Martínez, K. Bharadwaj, R. Dassonneville, J. Delaforce, F. Foroughi, V. Milchakov, L. Planat, O. Buisson, C. Naud, W. Hasch-Guichard, S. Florens, I. Snyman, and N. Roch, Nat. Commun. 10, 5259 (2019).
  * (7) A. Vrajitoarea, R. Belyansky, R. Lundgren, S. Whitsitt, A. V. Gorshkov, and A. A. Houck, arXiv:2209.14972 .
  * Ciuti  _et al._ (2005) C. Ciuti, G. Bastard, and I. Carusotto, Phys. Rev. B 72, 115303 (2005).
  * Stassi  _et al._ (2013) R. Stassi, A. Ridolfo, O. Di Stefano, M. J. Hartmann, and S. Savasta, Phys. Rev. Lett. 110, 1 (2013).
  * (10) Q.-K. He, Z. An, H.-J. Song, and D. L. Zhou, arXiv:1810.04523v2 .
  * Gheeraert  _et al._ (2018) N. Gheeraert, X. H. Zhang, T. Sépulcre, S. Bera, N. Roch, H. U. Baranger, and S. Florens, Phys. Rev. A 98, 1 (2018).
  * Sánchez-Burillo  _et al._ (2019) E. Sánchez-Burillo, L. Martín-Moreno, J. J. García-Ripoll, and D. Zueco, Phys. Rev. Lett. 123, 013601 (2019).
  * Cirio  _et al._ (2017) M. Cirio, K. Debnath, N. Lambert, and F. Nori, Phys. Rev. Lett. 119, 1 (2017).
  * Peropadre  _et al._ (2013) B. Peropadre, D. Zueco, D. Porras, and J. J. García-Ripoll, Phys. Rev. Lett. 111, 1 (2013).
  * Shi  _et al._ (2018) T. Shi, Y. Chang, and J. J. García-Ripoll, Phys. Rev. Lett. 120, 1 (2018).
  * Sanchez-Burillo  _et al._ (2014) E. Sanchez-Burillo, D. Zueco, J. J. Garcia-Ripoll, and L. Martin-Moreno, Phys. Rev. Lett. 113, 1 (2014).
  * Sánchez-Burillo  _et al._ (2015) E. Sánchez-Burillo, J. García-Ripoll, L. Martín-Moreno, and D. Zueco, Faraday Discuss. 178, 335 (2015).
  * Stassi  _et al._ (2017) R. Stassi, V. Macrì, A. F. Kockum, O. Di Stefano, A. Miranowicz, S. Savasta, and F. Nori, Phys. Rev. A 96, 023818 (2017).
  * Kockum  _et al._ (2017a) A. F. Kockum, V. Macrì, L. Garziano, S. Savasta, and F. Nori, Sci. Rep. 7, 5313 (2017a).
  * Chen  _et al._ (2017) Z. Chen, Y. Wang, T. Li, L. Tian, Y. Qiu, K. Inomata, F. Yoshihara, S. Han, F. Nori, J. S. Tsai, and J. Q. You, Phys. Rev. A 96, 1 (2017).
  * Kockum  _et al._ (2017b) A. F. Kockum, A. Miranowicz, V. Macrì, S. Savasta, and F. Nori, Phys. Rev. A 95, 1 (2017b).
  * Frisk Kockum  _et al._ (2019) A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys. 1, 19 (2019).
  * Forn-Díaz  _et al._ (2019) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 25005 (2019).
  * De Bernardis  _et al._ (2018) D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl, Phys. Rev. A 98, 1 (2018).
  * Di Stefano  _et al._ (2019) O. Di Stefano, A. Settineri, V. Macrì, L. Garziano, R. Stassi, S. Savasta, and F. Nori, Nat. Phys. 15 (2019).
  * Ashida  _et al._ (2022) Y. Ashida, T. Yokota, A. İmamoğlu, and E. Demler, Phys. Rev. Res. 4, 023194 (2022).
  * Savasta  _et al._ (2021) S. Savasta, O. Di Stefano, A. Settineri, D. Zueco, S. Hughes, and F. Nori, Phys. Rev. A 103, 053703 (2021).
  * Taylor  _et al._ (2022) M. A. D. Taylor, A. Mandal, and P. Huo, Opt. Lett. 47, 1446 (2022).
  * Sánchez Muñoz  _et al._ (2018) C. Sánchez Muñoz, F. Nori, and S. De Liberato, Nat. Commun. 9 (2018).
  * Settineri  _et al._ (2021) A. Settineri, O. Di Stefano, D. Zueco, S. Hughes, S. Savasta, and F. Nori, Phys. Rev. Res. 3, 1 (2021).
  * Mercurio  _et al._ (2022) A. Mercurio, V. Macrì, C. Gustin, S. Hughes, S. Savasta, and F. Nori, Phys. Rev. Res. 4, 023048 (2022).
  * Salmon  _et al._ (2022) W. Salmon, C. Gustin, A. Settineri, O. Di Stefano, D. Zueco, S. Savasta, F. Nori, and S. Hughes, Nanophotonics 11, 1573 (2022).
  * Dmytruk and Schiró (2021) O. Dmytruk and M. Schiró, Phys. Rev. B 103, 1 (2021).
  * Román-Roche and Zueco (2022) J. Román-Roche and D. Zueco, SciPost Phys. Lect. Notes , 50 (2022).
  * Sánchez-Burillo  _et al._ (2018) E. Sánchez-Burillo, A. Cadarso, L. Martín-Moreno, J. J. García-Ripoll, and D. Zueco, New Journal of Physics 20, 013017 (2018).
  * Silbey and Harris (1983) R. Silbey and R. A. Harris, J. Chem. Phys. 80, 2615 (1983).
  * Bera  _et al._ (2014) S. Bera, A. Nazir, A. W. Chin, H. U. Baranger, and S. Florens, Phys. Rev. B 90, 1 (2014).
  * Díaz-Camacho  _et al._ (2016) G. Díaz-Camacho, A. Bermudez, and J. J. García-Ripoll, Phys. Rev. A 93, 1 (2016).
  * Sánchez-Burillo  _et al._ (2019) E. Sánchez-Burillo, L. Martín-Moreno, J. J. García-Ripoll, and D. Zueco, Phys. Rev. Lett. 123, 1 (2019).
  * Legget  _et al._ (1987) A. J. Legget, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * García-Ripoll (2006) J. J. García-Ripoll, New Journal of Physics 8, 305 (2006).
  * Román-Roche  _et al._ (2020) J. Román-Roche, E. Sánchez-Burillo, and D. Zueco, Phys. Rev. A 102, 023702 (2020).
  * Zhou  _et al._ (2008) L. Zhou, Z. R. Gong, Y. X. Liu, C. P. Sun, and F. Nori, Phys. Rev. Lett. 101, 1 (2008).
  * Cohen-Tannoudji  _et al._ (1998) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, “Nonperturbative calculation of transition amplitudes,” in  _Atom—Photon Interactions_ (John Wiley & Sons, Ltd, 1998) Chap. 3, pp. 165–255.
  * Economou (2006) E. N. Economou, “Green’s functions for tight-binding hamiltonians,” in  _Green’s Functions in Quantum Physics_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 77–110.
  * Lombardo  _et al._ (2014) F. Lombardo, F. Ciccarello, and G. M. Palma, Phys. Rev. A 89, 053826 (2014).
  * Devoret (1997) M. H. Devoret, in  _Fluctuations Quantiques/Quantum Fluctuations_ , edited by S. Reynaud, E. Giacobino, and J. Zinn-Justin (1997) pp. 351–386.
  * García Ripoll (2022) J. J. García Ripoll,  _Quantum Information and Quantum Optics with Superconducting Circuits_ (Cambridge University Press, 2022).
  * Liu and Houck (2017) Y. Liu and A. A. Houck, Nat. Phys. 13, 48 (2017).
  * Shen and Fan (2009) J. T. Shen and S. Fan, Phys. Rev. A 79, 23818 (2009).