# High fidelity state reconstruction of a qubit via dynamics of a dissipative
resonator

G. Di Bello1 L. M. Cangemi2 V. Cataudella3,4 G. De Filippis3,4 A. Nocera5 C.
A. Perroni3 1Dip. di Fisica E. Pancini - Università di Napoli Federico II -
I-80126 Napoli, Italy  2Department of Chemistry Bar-Ilan University, Ramat-Gan
52900, Israel  3SPIN-CNR and Dip. di Fisica E. Pancini - Università di Napoli
Federico II - I-80126 Napoli, Italy  4INFN, Sezione di Napoli - Complesso
Universitario di Monte S. Angelo - I-80126 Napoli, Italy  5Department of
Physics and Astronomy and Stewart Blusson Quantum Matter Institute, University
of British Columbia, Vancouver, B.C., Canada, V6T 1Z1

###### Abstract

Using a dissipative quantum Rabi model, we study the dynamics of a slow qubit
coupled to a fast quantum resonator embedded in a viscous fluid from weak to
strong and ultra-strong coupling regimes. Solving the quantum Heisenberg
equations of motion, perturbative in the internal coupling between qubit and
resonator, we derive functional relationships directly linking the qubit
coordinates on the Bloch sphere to resonator observables. We then perform
accurate Matrix Product State simulations using the time dependent variational
principle, and compare our results both with the analytical solutions obtained
via the Heisenberg equations of motion, and with numerical solutions of a
Lindblad master equation, perturbative in the external coupling between
resonator and environment. The aim is to extend the standard qubit readout,
limited to a closed system for a qubit in a state along one quantization axis,
to the case of a spin directed along an arbitrary direction as a part of an
open system. We demonstrate that a qubit state pointing to generic coordinates
on the Bloch sphere can be fully reconstructed, up to the strong coupling
regime with a damped resonator. Indeed, by monitoring the qubit fidelity with
respect to free limit, we point out that, up to the strong coupling regime,
the qubit state accurately fulfills the functional relationships of the
Heisenberg equations of motion. Interestingly, a weak to intermediate
resonator coupling to the bath is able to anticipate the state reconstruction
of the qubit with respect to the same procedure for the closed system.
Moreover, we find a parameter regime in which, even if the standard qubit
readout fails, the resonator is still able to monitor the qubit dynamics.
Finally, in the ultra-strong coupling regime, non-Markovian effects become
robust and the dynamics of qubit and resonator are inextricably entangled
making the qubit state reconstruction difficult. Our work is not only
interesting from a fundamental point of view, but also for quantum control
issues, aimed to provide information about system dynamics through
measurements.

##  I Introduction

The interaction between matter and electromagnetic fields is one of the most
fundamental processes occurring in nature. This interaction, specified to
quantized electromagnetic fields, is the central focus of investigation in
cavity quantum electrodynamics (QED) and in circuit QED, with artificial atoms
and on-chip resonators. This research field is really important for both
atomic physics and quantum optics Gu  _et al._ (2017a). The quantum Rabi model
Rabi (1936) describes the most used setup, in which a two-level system (qubit)
is coupled to a single-mode resonator (cavity) field. This model gives a
mathematical framework for various quantum phenomena, inter alia,
nonclassical-state generation, quantum state transfer and topological physics
using photons. In particular, it is more general than the Jaynes-Cummings
model Jaynes and Cummings (1963), since the Rabi model explicitly includes the
counter-rotating terms in the qubit-resonator interaction Zueco  _et al._
(2009). Hence, it is crucial for describing strong and ultra-strong coupling
regimes in which the qubit-cavity coupling becomes about an order of magnitude
smaller than the characteristic frequency of both qubit and cavity. These
models help to interpret the experiments that have reached the strong coupling
regime for various QED setups Wallraff  _et al._ (2004) and the ultra-strong
one in superconducting QED. In the latter case the coupling between artificial
atoms and resonators is obtained, for instance, by inductively coupling a flux
qubit and an LC oscillator via Josephson junctions Yoshihara  _et al._ (2017).
The ultra-strong coupling regimes of light-matter interaction will keep on
expanding at the frontier of quantum optics and quantum physics. Thus, topics
related to these regimes will remain a prominent field in the foreseeable
future Forn-Díaz  _et al._ (2019). The other issue to address for a realistic
description of such systems is assessing quantitatively the role of
dissipation and decoherence Gu  _et al._ (2017a); Perroni  _et al._ (2016);
Nocera  _et al._ (2016); Perroni and Benenti (2021).
Coupled qubit-resonator systems are currently employed in circuit QED setups
since the qubit readout through the Stark dynamical shift of the resonator
frequency provides information on the spin component along the quantization
axis Gu  _et al._ (2017a). This procedure is strictly valid only in the weak
coupling limit between qubit and resonator and when interactions with the
environment are neglected. Actually, in the weak coupling regime, the closed
Rabi model reduces to the closed Jaynes-Cummings model Jaynes and Cummings
(1963), for which it is possible to show that the splitting of the energy
eigenvalues depends on qubit state, allowing its readout through the resonator
state Zagoskin (2011). It is still unclear if efficient qubit readout is
possible also when the inevitable coupling to the environment is considered in
the qubit-resonator model.
Recently, the readout procedure has been generalized in some experiments
Hacohen-Gourgy  _et al._ (2016); Campagne-Ibarcq  _et al._ (2016); Ficheux
_et al._ (2018) showing that it is possible to follow simultaneously non-
commuting spin components of a qubit through resonator observables. In
particular, in Ref. Ficheux  _et al._ , 2018, the procedure for a full state
reconstruction of a superconducting qubit has been improved. These recent
experiments have inspired the theoretical work presented in this paper. We
will discuss in the concluding section IV how results exposed in this paper
can be compared to those obtained in the experiments Hacohen-Gourgy  _et al._
(2016); Campagne-Ibarcq  _et al._ (2016); Ficheux  _et al._ (2018).

In this paper, we make an accurate theoretical analysis of static and dynamics
properties of the open quantum Rabi system, from the weak up to the ultra-
strong coupling regimes. The cavity frequency is greater than the qubit one,
so that the system is in the anti-adiabatic regime. We use both numerically
exact and approximate methods for solving the dynamics with the aim of
establishing the validity limits of different methods Breuer  _et al._ (2016);
Weimer  _et al._ (2021). Firstly, we solve the Heisenberg equations of motions
(HEM) in the weak coupling limit between the qubit and the cavity, in the
presence of Ohmic bath. These allow us to derive analytical functional
relationships directly linking cavity and qubit coordinates. Hence, to test
the validity of the derived relationships, we perform extensive Matrix Product
State (MPS) simulations. We further corroborate our MPS numerical results by
solving Lindblad master equations (LME) for the coupled qubit-cavity system.
Only MPS simulations provide correct results in every parameter region, being
variationally exact in both the qubit-cavity coupling and coupling to the
environment. Indeed, LME is obtained in the weak coupling limit between the
cavity and the environment and so it is valid under the Markovian
approximation Zagoskin (2011). Therefore, in the following, LME results help
us assess the non-Markovianity of the solutions obtained with MPS numerics.
Up to strong coupling regime, we retrieve the typical qubit readout by making
a Fourier transform of the average excitation number of the oscillator for a
qubit in the up state. Then, we extend the readout to a generic superposition
of up and down qubit states including the effects of the bath. Furthermore, we
analyse the resonator dynamics in order to reconstruct the qubit state during
its evolution. In fact, in addition to the number of bosons, we compute the
mean values of the oscillator coordinates. We find that a full qubit state
reconstruction can be achieved from weak to strong coupling regime if the
resonator is not strongly damped. We quantify the quality of the quantum state
reconstruction by monitoring the qubit fidelity with respect to free limit. We
show that, up to the strong coupling regime, we can infer accurately the qubit
state from the functional relationships of the Heisenberg equations of motion.
Interestingly, in the strong coupling regime, when the qubit-cavity coupling
is not too strong, the standard readout fails to provide the z spin component,
but a full state reconstruction is still possible by following the resonator
coordinates and excitation number over time. Quite surprisingly, results
reveal that dissipation and decoherence can facilitate the prediction of qubit
dynamics in the parameter regime from weak to strong coupling.
Finally, in the ultra-strong coupling regime, the qubit dynamics cannot be
separated from the oscillator evolution due to an enhancement of their
entanglement. As a result, the qubit state reconstruction from the resonator
dynamics becomes prohibitive. The Wigner quasi-probability distribution for
the resonator shows that it is in a damped momentum squeezed state, whose
squeezing parameter depends on both the internal and the external couplings.
Moreover, from the comparison with MPS simulations, we observe how Lindblad
equation fails to provide correct results due to the non-Markovian effects of
the interaction.
The results presented in this work are not limited to the anti-adiabatic
regime. Actually, we will show how a full quantum state reconstruction can be
extended to a parameter regime where the qubit energy scale is of the order of
that of the oscillator in the weak to strong coupling regime.
The present work can be useful for better understanding the dissipative
quantum Rabi model from a fundamental point of view. Furthermore, it can help
to address quantum control issues, which are aimed to gain information about
system dynamics through measurements enabling more powerful performance in
computing, sensing, and metrology.

The paper is organised as follows: in Sec. II we will describe the dissipative
quantum Rabi model and the methods used to solve the dynamics of this open
quantum system; in Sec. III, we will discuss our results for three regimes of
parameters: weak, strong, and ultra-strong coupling, which are summarized in
the Fig. 12 in the concluding section IV.

##  II  Dissipative quantum Rabi model and solution methods

In this section, we describe the quantum Rabi model, where a quantized single-
mode field is coupled to a two-level system (TLS) through a quantum dipolar
light-matter interaction:

$\mathcal{H}_{closed}=\,\frac{\Delta}{2}\sigma_{z}+\omega_{0}a^{\dagger}a+g\sigma_{x}(a+a^{\dagger}),$
(1)

where $\Delta$ is the tunneling matrix element, $a\,(a^{\dagger})$ is the
annihilation (creation) operator for the oscillator with frequency
$\omega_{0}$, $g$ stands for the strength of the qubit-cavity coupling and
$\sigma_{z}$ and $\sigma_{x}$ are the canonical Pauli matrices. The effects of
decoherence and dissipation induced by the electromagnetic environment are
taken into account via a linear coupling à la Caldeira-Leggett to a collection
of $N$ independent bosonic modes at zero temperature:

$\mathcal{H}_{bath}=\,\sum_{j=1}^{N}\left[\omega_{j}a^{\dagger}_{j}a_{j}-(a+a^{\dagger})\left|c_{j}\right|(a_{j}+a^{\dagger}_{j})+\frac{x_{0}^{2}}{2}M_{j}\omega_{j}^{2}\right].$
(2)

The bosonic modes have frequencies $\omega^{2}_{j}=k_{j}/M_{j}$, coordinates
and momenta given by $x_{j}$ and $p_{j}$, respectively; furthermore $x_{0}$
denotes the position operator of the oscillator with mass $m$ and frequency
$\omega_{0}$: $x_{0}=\sqrt{\frac{\hbar}{2m\omega_{0}}}(a+a^{\dagger})$. Units
are such that $\hbar=k_{B}=1$. The coupling constants to the bath are
$\left|c_{j}\right|=\sqrt{\frac{k_{j}\omega_{j}}{4m\omega_{0}}}$. Moreover, we
have neglected the energy shift $\sum_{j=1}^{N}\omega_{j}/2$, which does not
affect the dynamics. The last term ensures that the total energy is bounded
from below, irrespective of the strength of the coupling to the bath. Hence,
we can rewrite the Hamiltonian of the system plus the environment, by defining
a renormalized oscillator frequency
$\tilde{\omega}_{0}=\sqrt{\omega_{0}^{2}+\sum_{j=1}^{N}M_{j}\omega_{j}^{2}/m}$
and renormalized coupling strengths
$\tilde{g}=g\sqrt{\frac{\omega_{0}}{\tilde{\omega}_{0}}}$, between qubit and
resonator, and
$\left|\tilde{g}_{j}\right|=\sqrt{\frac{k_{j}\omega_{j}}{4m\tilde{\omega}_{0}}}$,
between oscillator and each bath bosonic mode. The total Hamiltonian then
reads:

$\displaystyle\mathcal{H}=\,$
$\displaystyle\frac{\Delta}{2}\sigma_{z}+\tilde{\omega}_{0}b^{\dagger}b+\tilde{g}\sigma_{x}(b+b^{\dagger})$
(3) $\displaystyle+\,$
$\displaystyle\sum_{j=1}^{N}\omega_{j}a^{\dagger}_{j}a_{j}-(b+b^{\dagger})\sum_{j=1}^{N}\left[\left|\tilde{g}_{j}\right|(a_{j}+a^{\dagger}_{j})\right],$

where $b\,(b^{\dagger})$ is the annihilation (creation) operator for the
renormalized oscillator with frequency $\tilde{\omega}_{0}$ and coordinates
$x=\sqrt{\frac{1}{2m\tilde{\omega}_{0}}}(b+b^{\dagger})$ and
$p=i\sqrt{\frac{m\tilde{\omega}_{0}}{2}}(b^{\dagger}-b)$. The bath is then
represented by an Ohmic spectral density:
$J(\omega)=\sum_{j=1}^{N}\left|\tilde{g}_{j}\right|^{2}\delta(\omega-\omega_{j})=\frac{\alpha}{2}\omega\Theta(\omega_{c}-\omega)$,
$\omega_{c}$ is the cutoff frequency, and $\theta(x)$ Heaviside function. Here
the dimensionless parameter $\alpha$ measures the strength of the resonator-
bath coupling (see Fig. 1).
This model can also be mapped Zueco and García-Ripoll (2019) in such a way
that, by including the resonator as a further oscillator of the bath, the
qubit is coupled to the N+1 bath modes as follows:

$\mathcal{H}_{mapped}=\frac{\Delta}{2}\sigma_{z}+\sum_{l=1}^{N+1}\tilde{\omega}_{l}b^{\dagger}_{l}b_{l}+\sigma_{x}\sum_{l=1}^{N+1}\left[\left|\beta_{l}\right|(b_{l}+b^{\dagger}_{l})\right],$
(4)

where now the N+1 bosonic modes are defined by the operators $b_{l}$ and
$b^{\dagger}_{l}$ and have frequencies $\tilde{\omega}_{l}$. The couplings
between the qubit and each bosonic mode $\left|\beta_{l}\right|$ are defined
such that we can describe the bath in terms of an effective spectral density:

$J_{eff}(\omega)=\sum_{l=1}^{N+1}\left|\beta_{l}\right|^{2}\delta(\omega-\tilde{\omega}_{l})=\frac{2g^{2}\omega_{0}^{2}\omega\alpha}{(\omega_{0}^{2}-\omega^{2})^{2}+(\omega_{0}\alpha\pi\omega)^{2}},$
(5)

which is Ohmic at low frequencies:
$J_{eff}(\omega)\approx\frac{2g^{2}\alpha}{\omega_{0}^{2}}\omega$. Therefore,
the qubit is coupled to an oscillator bath through an effective constant
proportional to $g^{2}\alpha/\omega_{0}^{2}$. Indeed, in the anti-adiabatic
regime ($\omega_{0}$ large), when the qubit-resonator coupling $g$ is small,
the qubit is only weakly coupled to its surroundings.

Figure 1: Dissipative quantum Rabi model described by the Hamiltonian in Eq.
(1) for a qubit of frequency $\Delta$ coupled through $g$ to a cavity of
frequency $\omega_{0}$, in turn coupled to a bath via a coupling strength
$\alpha$.

In the following subsections we are going to show the approaches useful to
describe the time-evolution of the systems of interest (qubit and oscillator)
when coupled with the environment through the cavity. First, we will find the
solution of the HEM for position and momentum of the cavity, perturbative in
the coupling strength $g$. Then, in the second subsection, we present the
dynamic simulations through MPS method. In Appendix, we expose the solution of
global LME for the composite system qubit plus oscillator, that we will use to
compare with HEM and MPS results. We underline the regimes in which the three
methods provide a correct description. In particular, we are interested in the
cases that allow to restore and extend the typical qubit readout. Moreover,
thanks to the MPS simulations, we will be able to verify the functional
relationships provided by HEM, useful to unveil the regimes where a full qubit
state reconstruction is possible, which are summarized in the phase diagram in
Fig. 12. Actually, MPS allows us to go beyond the other two perturbative
methods, correctly analysing strong and ultra-strong coupling regimes, when
the effect of the bath is not negligible.

###  II.1 Perturbative solution of HEM

This approach is possible since our bath is made up of bosonic modes linearly
coupled to the system made of qubit plus resonator. Hence, we can exactly
eliminate the dynamics of the bath from HEM for the variables pertaining to
the reduced system (qubit and oscillator). Starting from the Hamiltonian
$\mathcal{H}$ in Eq. (3) and evaluating the HEM
$\dot{A}(t)=i\left[\mathcal{H},A(t)\right]$ for a generic observable $A$
belonging to the qubit-resonator-bath system, we find:

$\begin{cases}\dot{x}_{j}(t)=\frac{p_{j}(t)}{M_{j}}\\\
\dot{p}_{j}(t)=-M_{j}\omega_{j}^{2}x_{j}(t)+\lambda_{j}x(t)\,;\end{cases}\\\ $
(6) $\begin{cases}\dot{x}(t)=\frac{p(t)}{m}\\\
\dot{p}(t)=-m\tilde{\omega}_{0}^{2}x(t)-\tilde{g}\sigma_{x}(t)\sqrt{2m\tilde{\omega}_{0}}+\sum_{j=1}^{N}\lambda_{j}x_{j}(t).\end{cases}$
(7)

The first system of equations refers to the bath positions and momenta (Eq.
(6)), while the second one is for the oscillator observables (Eq. (7)). The
coefficients that appear in both the systems are
$\lambda_{j}=-\tilde{g}_{j}\sqrt{2m\tilde{\omega}_{0}}\sqrt{2M_{j}\omega_{j}}$.
These differential equations are coupled, such that we can express for example
$x(t)$ from the first system as a function of the $x_{j}(t)$ and substitute it
in the second one. Therefore, we obtain an equation for $x(t)$ (and similarly
for $p(t)$) as follows:

$\ddot{x}(t)+\tilde{\omega}_{0}^{2}x(t)+\tilde{g}\sigma_{x}(t)\sqrt{\frac{2\tilde{\omega}_{0}}{m}}+\frac{d}{dt}\int_{0}^{t}ds\gamma(t-s)x(s)=\frac{B(t)}{m}.$
(8)

It may be viewed as the quantum analogue of a classical stochastic
differential equationBreuer  _et al._ (2002), involving a damping kernel

$\gamma(t-s)=4\tilde{\omega}_{0}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega}\cos\left[\omega\left(t-s\right)\right]$
(9)

and a force operator

$B(t)=\sum_{j=1}^{N}\frac{\lambda_{j}}{\sqrt{2M_{j}\omega_{j}}}\left(a_{j}e^{-i\omega_{j}t}+a^{\dagger}_{j}e^{i\omega_{j}t}\right),$
(10)

that on average becomes a stochastic force, whose statistics depends on the
initial bath states distribution. In our case, when the Ohmic spectral density
has an infinite cutoff, $\omega_{c}\rightarrow\infty$, we get
$\gamma(t-s)\rightarrow 4\gamma\delta(t-s)$, where
$\gamma=\pi\frac{\alpha}{2}\tilde{\omega}_{0}$ is the cavity decaying rate. In
all the approaches used in this work, we consider the limit where $\omega_{c}$
is the largest frequency scale. This approximation allows us to simplify the
derivative of the integral with the damping kernel to $2\gamma\dot{x}(t)$.
Thus, we perform the average over the initial distribution $\rho(0)$ and find
the HEM for the oscillator coordinates:

$\displaystyle\left\langle\ddot{x}(t)\right\rangle+\tilde{\omega}_{0}^{2}\left\langle
x(t)\right\rangle+2\gamma\left\langle\dot{x}(t)\right\rangle+\tilde{g}\left\langle\sigma_{x}(t)\right\rangle\sqrt{\frac{2\tilde{\omega}_{0}}{m}}=$
$\displaystyle\,0$ (11)
$\displaystyle\left\langle\ddot{p}(t)\right\rangle+\tilde{\omega}_{0}^{2}\left\langle
p(t)\right\rangle+2\gamma\left\langle\dot{p}(t)\right\rangle-\tilde{g}\left\langle\sigma_{y}(t)\right\rangle\sqrt{2\tilde{\omega}_{0}m}=$
$\displaystyle\,0\,,$ (12)

where the third terms are friction forces and the fourth ones are “external”
forces due to the interaction with the qubit. Clearly the solution of these
coupled differential equations ($\left\langle
x(t)\right\rangle=\frac{\left\langle p(t)\right\rangle}{m}$) depends on the
time evolution of the mean values of the qubit Pauli matrices. A way to find
an analytical solution in the anti-adiabatic regime
($\Delta\ll\tilde{\omega}_{0}$) is to consider the ratio between the qubit-
oscillator coupling and the resonator frequency,
$\tilde{g}/\tilde{\omega}_{0}$ as a weak perturbation. By neglecting terms
greater than or equal to the first order in $\tilde{g}/\tilde{\omega}_{0}$ in
the resonator equations, the HEM for the Pauli matrices can be solved at zero
order in $\tilde{g}/\tilde{\omega}_{0}$, describing the well-known Rabi
oscillations for a free evolving qubit. Apart from a decaying homogeneous
solution on a transient time of the order of $~{}\gamma^{-1}$, the solutions
for the mean dimensionless oscillator position
$\tilde{x}=x\sqrt{2m\tilde{\omega}_{0}}$ and momentum
$\tilde{p}=p\sqrt{2/(m\tilde{\omega}_{0})}$ are the following:

$\displaystyle\left\langle\tilde{x}(t)\right\rangle\rightarrow$
$\displaystyle\,-2\tilde{\omega}_{0}k\left\langle\sigma_{x}(t)\right\rangle$
(13) $\displaystyle\left\langle\tilde{p}(t)\right\rangle\rightarrow$
$\displaystyle\,2\Delta k\left\langle\sigma_{y}(t)\right\rangle,$ (14)

where the coefficient $k$ is

$k=\frac{\tilde{g}(\tilde{\omega}_{0}^{2}-\Delta^{2})}{\Delta^{4}+\tilde{\omega}_{0}^{4}+4\gamma^{2}\Delta^{2}-2\Delta^{2}\tilde{\omega}_{0}^{2}}.$
(15)

At the zero-order in the coupling $\tilde{g}$, the mean Pauli matrices are
those of a free spin

$\displaystyle\left\langle\sigma_{x}(t)\right\rangle=$
$\displaystyle\,\left\langle\sigma_{x}(0)\right\rangle\cos(\Delta
t)-\left\langle\sigma_{y}(0)\right\rangle\sin(\Delta t)$ (16)
$\displaystyle\left\langle\sigma_{y}(t)\right\rangle=$
$\displaystyle\,\left\langle\sigma_{y}(0)\right\rangle\cos(\Delta
t)+\left\langle\sigma_{x}(0)\right\rangle\sin(\Delta t)$ (17)
$\displaystyle\left\langle\sigma_{z}(t)\right\rangle=$
$\displaystyle\,\left\langle\sigma_{z}(0)\right\rangle.$ (18)

In Sec. III, we will show that, in the weak coupling regime, this treatment of
the mean values is sufficient to get accurate results.
In the next section, thanks to the MPS simulations, we will show that Eqs.
(13-14) are valid up to the strong coupling regime. Moreover, we will show
that the main effects of the non-zero $\tilde{g}$ are to introduce a
renormalized frequency ($\tilde{\Delta}_{x}$ and $\tilde{\Delta}_{y}$ for
$\left\langle\sigma_{x}(t)\right\rangle$ and
$\left\langle\sigma_{y}(t)\right\rangle$ respectively) and a decay rate
($\kappa_{x}$ and $\kappa_{y}$ for $\left\langle\sigma_{x}(t)\right\rangle$
and $\left\langle\sigma_{y}(t)\right\rangle$ respectively) in Eqs. (16-17):

$\displaystyle\left\langle\sigma_{x}(t)\right\rangle=$
$\displaystyle\,\left(\left\langle\sigma_{x}(0)\right\rangle\cos(\tilde{\Delta}_{x}t)-\left\langle\sigma_{y}(0)\right\rangle\sin(\tilde{\Delta}_{x}t)\right)e^{-\kappa_{x}t}$
(19) $\displaystyle\left\langle\sigma_{y}(t)\right\rangle=$
$\displaystyle\,\left(\left\langle\sigma_{y}(0)\right\rangle\cos(\tilde{\Delta}_{y}t)+\left\langle\sigma_{x}(0)\right\rangle\sin(\tilde{\Delta}_{y}t)\right)e^{-\kappa_{y}t},$
(20)

while for the z component, a good fit is obtained by using a linear regression
with slope $a_{z}$ ($<0$) and intercept $b_{z}$:

$\left\langle\sigma_{z}(t)\right\rangle=\,a_{z}t+b_{z}.$ (21)

We will show that the relationships above fit very well the results obtained
with MPS approach in Sec. III.
There remains the problem of reconstructing the z-component of the qubit
trajectory on the Bloch sphere. An approximate time evolution can be written
starting from the average excitation number for the oscillator in terms of its
coordinates and their standard deviations $\delta\tilde{x}$ and
$\delta\tilde{p}$. By using Eqs. (13-14), one obtains:

$\displaystyle\left\langle\tilde{n}(t)\right\rangle=$
$\displaystyle\,\tilde{\omega}_{0}^{2}k^{2}(1-\left\langle\sigma_{z}(t)\right\rangle^{2})+k^{2}(\tilde{\omega}_{0}^{2}-\Delta^{2})\left\langle\sigma_{y}(t)\right\rangle^{2}$
$\displaystyle+$
$\displaystyle\frac{1}{4}\left((\delta\tilde{x}(t))^{2}+(\delta\tilde{p}(t))^{2}\right)-\frac{1}{2}.$
(22)

Since $\left\langle\sigma_{y}(t)\right\rangle$ can be derived from previous
Eqs. (13-14), we get access to $\left\langle\sigma_{z}(t)\right\rangle$
through the resonator observables over time. In particular, as found in the
next section, the third term in Eq. (22) is only a small correction, not
strongly dependent on time. Moreover, the mean number of bosons depends on the
square of spin components. Therefore, this quantity provides the same results
if the qubit is initially either in the up or down state.
In the case of a qubit state initially up along z direction, Eqs. (16-17), but
also Eqs. (19-20), vanish, therefore, Eq. (22) provides a direct link between
the average excitation number and the z component of the spin. This is
equivalent to the standard qubit readout. On the other hand, if the qubit is
initially in a linear combination of up and down states along z direction,
Eqs. (13-14) are different from zero, hence the mean number of bosons in Eq.
(22) displays a more complex time behaviour.
One of the aims of this paper is to extend the readout to the case of a qubit
initially in a linear combination of up and down states in the presence of a
bath. Ultimately, the approximate but analytical solutions in Eqs. (13-14-22)
are the functional relationships linking qubit and oscillator observables
which will guide us in the remaining part of the paper. These equations
explicitly show how it is possible with proper tuning to follow the qubit
state through the oscillator time evolution. We will follow essentially this
procedure in the next section: we determine oscillator observables from
independent MPS simulations, then we will use the functional HEM Eqs.
(13-14-22) to fit the MPS data for the oscillator observables by using MPS
results for the spin components. In particular, the parameters of Eqs. (19-20)
are obtained from MPS simulations. In the next section, we will call this
procedure as $HEM+fit$. We will find that, up to the strong coupling regime,
the agreement is excellent, implying that, by inverting Eqs. (13-14-22), the
oscillator observables provide access to all three spin components, therefore
to qubit state reconstruction.

####  II.1.1 Some remarks on the perturbative analytical solution

Let us have a closer look to the solution of the HEM derived above. We stress
two points which will be relevant in the next sections:

  * 1)

dissipation and decoherence induced by the bath on the oscillator are not
detrimental to the full qubit state reconstruction;

  * 2)

in the anti-adiabatic regime, with increasing oscillator-qubit coupling, the
oscillator is in a well approximated squeezed state, depending on the qubit
state, with a time dependent squeezing parameter $r(\alpha,\tilde{g})$.

The first consideration is due to the fact that the homogeneous solution of
the HME is proportional to the factor $e^{-\gamma t}$, hence, for
$\gamma\rightarrow\infty$ it vanishes. We expect that as the coupling to the
bath increases, the qubit state reconstruction can start earlier, on a time
scale of the order of $\tau\approx 3/\gamma$.
For the second point, it is useful to consider a standard canonical
transformation for the cavity at time $t$ such that the new annihilation
operator is $c=b+\frac{\tilde{g}}{\tilde{\omega}_{0}}\sigma_{x}$ and therefore
the new oscillator Hamiltonian is
$\mathcal{H}_{o}=\tilde{\omega}_{0}c^{\dagger}c-\frac{\tilde{g}^{2}}{\tilde{\omega}_{0}}$.
We obtain just the Hamiltonian of a coherent state centred at
$-\frac{\tilde{g}}{\tilde{\omega}_{0}}\sigma_{x}$. Within the anti-adiabatic
regime, the qubit is much slower than the cavity, therefore it can be followed
by the oscillator in a coherent state moving according to the
$\left\langle\sigma_{x}\right\rangle$ of the qubit. Indeed, the HEM for qubit
and oscillator are coupled so that $\sigma_{x}$ can be written as a term to
the zero-order in $\tilde{g}$ plus a term which depends on
$\tilde{g}^{2}\tilde{x}$. As a consequence, a term to the third order in
$\tilde{g}$ proportional to $\tilde{x}^{2}$ appears in the oscillator
Hamiltonian, indicative of a squeezed state, with squeezing parameter
$r(\alpha,\tilde{g})$. In particular, being the coupling to the qubit in the
oscillator position, we expect a momentum squeezed state.

###  II.2 Matrix Product States simulations

It is well known that three main issues complicate the study of the dynamics
of open quantum systems:

  * •

no exact analytical solution generally exists;

  * •

perturbation theory fails at strong coupling between the system and the bath;

  * •

size of the Hilbert space grows exponentially in the system size, making exact
diagonalization prohibitive.

Figure 2: MPS chain of sites representing the Hamiltonian in Eq. (3). The
first site is occupied by the qubit, the second one by the oscillator with
Hilbert space of dimension $N_{o}$ and the sites from the third to the
$(N+2)$_th form the Ohmic bath of bosonic modes, each with Hilbert space of
dimension $N_{ph}$. Qubit and oscillator are coupled via $\tilde{g}$, while
the oscillator and the N bath modes are coupled by long range interactions of
strength $\tilde{g}_{j}$.

For these reasons, accurate and efficient numerical studies are needed.
We solve the Hamiltonian in Eq. (3) using the time-dependent MPS simulations.
In particular, we adopt the star geometry depicted in Fig. 2 to describe the
long-range interactions between the qubit, the cavity mode, and the bath
modes.
Because of the long-range character of the interactions, we adopt two
different methods for the solution of the time-dependent Schrödinger equation.
The first one was developed in Ref. Zaletel  _et al._ , 2015 and consists in a
first order approximation of the unitary time-evolution operator in terms of a
Matrix Product Operator (MPO). This method, to which we refer to $W^{I}$ in
the following, has an error per site which diverges with the system size $L$,
while giving a time-step error of $O(dt^{2})$.
The second method we use is the time-dependent variational principle (TDVP)
Haegeman  _et al._ (2011, 2016); Paeckel  _et al._ (2019), where the time-
dependent Schrödinger equation is projected to the tangent space of the MPS
manifold of fixed bond dimension at the current time. In this work we employ
the two-site TDVP (2TDVP in Ref. Paeckel  _et al._ , 2019) using the second
order integrator by sweeping left-right-left with half time step $dt/2$. The
main advantage of this method is that it has a smaller time-step error
$O(dt^{3})$, and its accuracy is controlled only by the MPS bond dimension and
the threshold to terminate the Krylov series. In this work, we stop the Krylov
vectors recurrence when the total contribution of two consecutive vectors to
the matrix exponential is less than $10^{-12}$.
In the star-geometry considered in this work, depicted in Fig. 2, we have
placed the qubit on the first site, on the second one the resonator with
Hilbert space of dimension $N_{o}$ and on the remaining sites the collection
of $N$ bosonic modes of the bath, each with Hilbert space of dimension
$N_{ph}$. The coupling between qubit and oscillator is $\tilde{g}$, while the
oscillator and the bath experiment long-range interactions with couplings
$\tilde{g}_{j}$.
Using the ITensor library Fishman  _et al._ (2020), we start the time-
evolution from a product state with the cavity and the bath modes in their
vacuum state (zero temperature), while the qubit is placed in a generic point
on the surface of the Bloch sphere. In order to reduce the simulation time and
simultaneously reach a longer dynamics, we tested the $W^{I}$ and TDVP methods
against each other on the exactly solvable closed model. We observed that the
TDVP method reproduces the analytical exact solution by using a time step two
or three orders of magnitude larger than that needed by $W_{I}$. Moreover, as
expected, $W^{I}$ has shown a much smaller accuracy than the TDVP method.
Thus, we decided to use the TDVP method for our simulations in the presence of
the interaction with the bath modes. An important observation about our
numerical simulations is in order: we did not use more sophisticated
approaches like local basis optimization Zhang  _et al._ (1998); Bursill
(1999); Friedman (2000); Wong and Chen (2008); Brockt  _et al._ (2015); Jansen
_et al._ (2022). We instead have converged our simulations in the number of
Fock states in the cavity and bath modes, finding the best compromise between
the smallest bond dimension and longest simulation times. Our truncation error
has kept below $10^{-12}$ requiring a maximum bond dimension of
$D_{\rm{max}}=20$. At the same time, this optimal maximum bond dimension has
allowed us to reach a time size for the simulations of
$t_{\rm{final}}=35/\Delta$, about 5 periods of the qubit. All the optimal
parameters used in the MPS simulations are specified in the captions of the
figures of the section III. A similar convergence analysis has already been
done in previous works for the applyMPO method and for one De Filippis  _et
al._ (2020) or more De Filippis  _et al._ (2021) qubits interacting with a
bath.
We finally note that, in the star geometry, one could also adopt the TEBD
method with swap gates. It was recently shown, however, in Ref. Liu  _et al._
, 2022, that it usually requires larger bond dimensions compared to 2TDVP,
despite giving smaller accumulated errors for long time evolutions.

In the next section, we compare the results of our MPS simulations against the
perturbative solutions shown in the previous sections. Indeed, we have used
our MPS simulations to explore the parameter regions where LME or HEM are
expected to fail.

Figure 3:  Cavity position in $(a)$, cavity momentum in $(b)$ and cavity
number in $(c)$ as a function of time for two different methods: MPS
simulations (blue solid line) and the HEM approach (green dashed line)
(particular solutions in Eqs. (13-14)). Parameters:
$\omega_{0}=10\Delta,g=0.1\Delta,\alpha=0.05,\omega_{c}=50\Delta$; for MPS
simulations, $N_{o}=10,N=500,N_{ph}=3,D_{\rm{max}}=20$.

##  III Comparison of different methods

In this section we showcase our results by applying the previous methods to
the dynamics of our model, in the anti-adiabatic regime
($\Delta\ll\tilde{\omega}_{0}$), from the weak to the ultra-strong coupling
regime. Specifically, up to the strong coupling regime, for weak to
intermediate couplings of the resonator to the bath, we find that the
resonator allows the typical qubit readout (up or down initial qubit state).
Furthermore, through the time evolution of resonator observables, starting
from a generic qubit initial state, we reconstruct the qubit state by
evaluating the qubit fidelity in time with respect to its free evolution (see
Fig. 12 for the sum up of the results). For parameters from the strong to the
ultra-strong coupling regime, the interactions (between them and with the
bath) are so important that the qubit behaviour differs from the free one as a
result of the resonator back-action upon the qubit. However, our analysis
allows us to go beyond the Lindblad approximation, observing non-Markovian
features of the dynamics, even through the Wigner quasi-probability
distribution for the resonator state.
This section is organised as follows. In Subsec. III.1 we will show how the
state reconstruction of the qubit works in the weak coupling regime for the
three techniques (HEM, LME, MPS); we evaluate the qubit fidelity, and clarify
in what sense dissipation and decoherence are useful for the qubit state
reconstruction. In Subsec. III.2 we will identify the range of bath-resonator
couplings for which a typical readout is feasible in the strong coupling
regime, for both the up state and the generic one; moreover, we will continue
the analysis for the qubit state reconstruction. In Subsec. III.3, finally, we
will discuss the features of the dynamics in the ultra-strong coupling regime,
where entanglement and non-Markovianity are non-negligible and neither qubit
readout nor state reconstruction are achievable.

###  III.1 Weak coupling regime

We set the parameters such that the qubit frequency $\Delta$ is taken as
unity, hence, for $\omega_{0}=10\Delta$, $g=0.1\Delta$ and $\alpha=0.05$,
$g/\Delta\leq\Gamma/\Delta\,\bigcap\,\Gamma/\Delta\leq 0.1\omega_{0}/\Delta$,
with $\Gamma=\alpha\omega_{0}\pi/2$ the bare resonator decay rate, involving
that the quantum Rabi system is in the weak coupling regime Gu  _et al._
(2017b) (Fig. 12). We remark that, in the limit of small $\alpha$, the
agreement between HME, LME and MPS results is excellent. In this subsection,
we set the coupling $\alpha=0.05$, an intermediate value of the resonator-bath
coupling, since it allows to understand the general method used in the next
subsections.

We present our implementation of qubit state reconstruction, by following the
oscillator state. Starting from the state $\ket{0}\otimes\ket{\nearrow}$, with
$\ket{\nearrow}=(\cos(\theta/2)\ket{\uparrow}+\sin(\theta/2)\,e^{-i\phi}\ket{\downarrow})$,
where $\theta=0.3\pi$ and $\phi=1.2\pi$, and an empty bath, we use MPS and the
complete solutions of HEM (only particular solutions in Eqs. (13-14)). As
shown in Fig. 3, we find a good agreement between the two approaches. The
cavity starts following the qubit after $\tau\approx 3/\gamma\approx
3.12/\Delta$. The agreement is good for the oscillator coordinates, plotted in
the panels (a) and (b) of Fig. 3. There is a small difference in the amplitude
of the oscillations which is not quite well reproduced by the HEM approach,
because it is proportional to the qubit-resonator coupling constant. For the
number, instead, panel (c) of Fig. 3 shows an excellent agreement. We notice
how the number becomes constant, only shifted from $0$ to $0.0257$, after a
rapid oscillation.

Both LME and MPS show the possibility to follow the qubit and are consistent
with each other. They also confirm that the HEM solution can be improved by
using the fit for the mean values of spin observables given in Eqs. (19-20).
In fact, in panel $(a)$ of Fig. 4, we notice this agreement for the time
evolution of the oscillator position. For LME and HEM there is a difference
with MPS in the amplitude of the oscillations due to a more accurate treatment
of the coupling.

We assess the quality of the quantum state reconstruction by evaluating the
qubit fidelity $F_{q}(t)$ with respect to its free evolution, when decoupled
from the cavity:

$\rho^{free}(t)=\frac{1}{2}\left(\mathbf{I}+\left\langle\vec{\sigma}^{free}(t)\right\rangle\cdot\vec{\sigma}^{free}\right),$
(23)

where the free Bloch vector describes the Rabi oscillations and it is a pure
state. Thus, the qubit fidelity reads:

$\displaystyle F_{q}(t)=$
$\displaystyle\left(\Tr\sqrt{\sqrt{\rho^{free}(t)}\rho(t)\sqrt{\rho^{free}(t)}}\right)^{2}$
(24) $\displaystyle=$
$\displaystyle\bra{\psi_{\rho^{free}(t)}}\rho(t)\ket{\psi_{\rho^{free}(t)}}$
(25) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+\left\langle\vec{\sigma}^{free}(t)\right\rangle\cdot\left\langle\vec{\sigma}(t)\right\rangle\right).$
(26)

A fidelity above the 90 percent is assumed reliable. Indeed, in this coupling
regime, the state reconstruction is robust, as shown in panel $(b)$ of Fig. 4.

Figure 4:  Cavity position in $(a)$ and qubit fidelity in $(b)$ as a function
of time (in units of $1/\Delta$). In panel $(a)$ the cavity position is
computed with three methods: MPS simulations (blue solid line), LME solution
(red dashed line) and HEM approach (green dotted line), improved by the fit
given in Eqs. (19-20). Parameters:
$\omega_{0}=10\Delta,g=0.1\Delta,\alpha=0.05,\omega_{c}=50\Delta$; for MPS
simulations, $N_{o}=10,N=500,N_{ph}=3,D_{\rm{max}}=20$.

This panel also shows a comparison between MPS and Lindblad methods by
following a linear fit for both: $F_{q}(t)\approx a_{n}-\kappa_{n}t$, where
$n=M$ for fitting of MPS curve, $n=L$ for fitting of LME plot. The Lindblad
average fidelity decreases with a slope smaller in absolute value than that
computed with MPS. Anyway, both MPS and Lindblad solutions give a fidelity
close to $100\%$ up to $t\approx 5T$, where $T$ is the qubit period (see Fig.
12 and Table 1 in the concluding section IV).
With increasing $\alpha$, as expected, the HEM homogeneous solution
effectively drops earlier since the transient $\tau\approx 3/\gamma$ gets
reduced up to very strong couplings, when the interactions break the qubit-
resonator coherence. Furthermore, starting from the vacuum state, the
oscillator evolves in a momentum squeezed state, not so different from the
vacuum, due to the low value of $\tilde{g}$, as underlined in II.1.1. In
particular, the figure Fig. 5 shows that the position uncertainty
$\delta\tilde{x}(t)>1$, the momentum one is instead $\delta\tilde{p}(t)<1$,
indication of a momentum squeezed state, so that the product
$\delta\tilde{x}(t)\delta\tilde{p}(t)\geq 1$, as ruled by the uncertainty
principle.

Figure 5:  Cavity position uncertainty (blue solid line), momentum one (red
dashed line) and uncertainty product (green dotted line) as a function of time
(in units of $1/\Delta$) from MPS simulations. Parameters:
$\omega_{0}=10\Delta,g=0.1\Delta,\alpha=0.05,\omega_{c}=50\Delta$; for MPS
simulations, $N_{o}=10,N=500,N_{ph}=3,D_{\rm{max}}=20$.

###  III.2 Strong coupling regime - non-perturbative state reconstruction

We choose again the qubit frequency $\Delta$ as unity and
$\omega_{0}=10\Delta$, by tuning $g$ and $\alpha\in[0.01,0.1]$ such that
$g/\Delta\geq\Gamma/\Delta\,\bigcap\,g/\Delta\leq
0.1\omega_{0}/\Delta\,\bigcap\,\Gamma/\Delta\leq 0.1\omega_{0}/\Delta$ or only
$g/\Delta\leq 0.1\omega_{0}/\Delta$ if $\Gamma/\Delta\geq
0.1\omega_{0}/\Delta$ (see Fig. 12). Hence, the closed quantum Rabi system is
in the strong coupling regime. At the same time, the ratio between the
coupling and the sum of the qubit and oscillator frequencies of the closed
system is $g/(\Delta+\omega_{0})\approx 10^{-2}$, so that it is possible to
rewrite the Hamiltonian of the Rabi model (Eq. (1)) as the Bloch-Siegert one
Bloch and Siegert (1940)

$\displaystyle\mathcal{H}_{BS}$
$\displaystyle=\frac{\Delta}{2}\sigma_{z}+\omega_{0}a^{\dagger}a+\omega_{BS}\left[\sigma_{z}\left(n+\frac{1}{2}\right)-\frac{1}{2}\right]$
(27) $\displaystyle+G(n)a^{\dagger}\sigma_{-}+a\sigma_{+}G(n),$ (28)

with the Bloch-Siegert frequency $\omega_{BS}=g^{2}/(\Delta+\omega_{0})$ and
$G(n)=g[1-n\omega_{BS}/(\Delta+\omega_{0})]$. This Hamiltonian can be block
diagonalised at fixed total number of excitations $N$ obtaining the
eigenvalues Klimov and Chumakov (2009)

$\displaystyle E^{BS}_{N,\pm}$
$\displaystyle=(\Delta+\omega_{0})N\pm\Gamma_{N}-g\epsilon,$ (29)
$\displaystyle\Gamma_{N}$
$\displaystyle=\sqrt{\frac{(\Delta-\omega_{0}+2g\epsilon
N)^{2}}{4}+g^{2}N(1-\epsilon^{2}N)^{2}},$ (30)

where $\epsilon=g/(\Delta+\omega_{0})$. As for the Jaynes-Cummings model, the
splitting of the energy eigenvalues depends on the qubit state ($\pm$), which
can be read through the resonator state.
Another way to determine the excitation energies of the closed system is to
directly apply the stationary perturbation theory for low $g$ on the
Hamiltonian in Eq. (1). The unperturbed Hamiltonian reads:

$\mathcal{H}_{q}+\mathcal{H}_{o}\ket{n,\pm}^{(0)}=\pm\frac{\Delta}{2}\sigma_{z}+\omega_{0}n\ket{n,\pm}^{(0)}=E_{n,\pm}^{(0)}\ket{n,\pm}^{(0)},$
(31)

so that the eigenvalues up to the second order in $g$ are the following:

$E_{n,\pm}\approx\pm\frac{\Delta}{2}+\omega_{0}n\pm
g^{2}\left(\frac{n+1}{\Delta\mp\omega_{0}}+\frac{n}{\Delta\pm\omega_{0}}\right).$
(32)

Hence the excitation energies can be easily evaluated as:

$\displaystyle E_{u}$ $\displaystyle=E_{1,-}-E_{0,+}$ (33) $\displaystyle
E_{d}$ $\displaystyle=E_{1,+}-E_{0,-},$ (34)

where $E_{u}$ is the energy corresponding to the up state of the qubit, lesser
than the cavity energy, and $E_{d}$ that of the down one, greater than the
cavity energy.
The excitation energies are obviously in agreement with the Bloch-Siegert
Hamiltonian (Eq. (29)) approach. In fact, similarly to the Jaynes-Cummings
model Zagoskin (2011), the oscillator frequency $\omega_{0}$ is shifted
depending on the qubit state as follows:

$\omega_{0}\rightarrow\omega_{0}\pm\left[\Delta+g^{2}\left(\frac{1}{\Delta+\omega_{0}}+\frac{1}{\Delta-\omega_{0}}\right)\right],$
(35)

so that the ($-$) sign refers to the excitation energy $E_{u}^{BS}$, while the
($+$) to $E_{d}^{BS}$.

Figure 6:  FFT of the mean cavity number $\tilde{n}$ computed at $g=0.6$ as a
function of frequency (in units of $\Delta$) through MPS simulations for three
values of $\alpha$ and for the two initial states: ($\ket{\uparrow}$ (left
panel), $\ket{\nearrow}$ (right panel)). Parameters:
$\omega_{0}=10\Delta,g=0.6\Delta,\alpha\in[0.01,0.03],\omega_{c}=50\Delta$;
for MPS simulations, $N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$.

We first verify if the standard qubit readout through the dissipative cavity
works also in this coupling regime. In particular, we perform the same
analysis for a qubit in an up state ($\ket{0}\otimes\ket{\uparrow}$) or in a
generic state ($\ket{0}\otimes\ket{\nearrow}$), and an empty cavity in contact
with the bath at zero temperature. In order to locate the excitation energies
of the open system and read the qubit state from the resonator, we analyse the
peaks of the Fast Fourier Transform (FFT) of the mean oscillator number.
As expected, we observe in Fig. 6 only a peak for the up state, with frequency
smaller that the cavity one, and two peaks for the generic qubit state due to
the presence of up (left peak) and down (right peak) contributions. With
increasing the resonator-bath coupling $\alpha$, we find that the peaks become
less resolved making the qubit readout more difficult. In both panels of Fig.
6, there is a critical value $\alpha_{c}$ of $\alpha$, about $0.03$, above
which the qubit readout is not feasible. We see how the FFT below the critical
value for the generic initial state has memory of the down contribution.

We then pursue the full qubit state reconstruction choosing the initial state
as $\ket{0}\otimes\ket{\nearrow}$ analysing in detail the system behaviour for
$g=0.6$ (the same as in Fig. 6) and $g=1.0$, with $\alpha\in\\{0.01,0.1\\}$.
For $\alpha=0.01$, results from LME and MPS are still in agreement and the
full state reconstruction is possible (see Fig. 7).

Figure 7:  Cavity position (in panels $(a)$ and $(b)$) and qubit fidelity (in
panels $(c)$ and $(d)$) as a function of time (in units of $1/\Delta$) for
$g=0.6$ and $g=1.0$ at $\alpha=0.01$ using different methods: MPS simulations
(blue solid line), the LME solution (red dashed line) and the HEM approach
(green dotted line), improved by the fitting with spin observables obtained by
MPS simulations as in Eqs. (19-20). Parameters:
$\Delta=1,\omega_{0}=10\Delta,g=\\{0.6\Delta,1.0\Delta\\},\alpha=0.01,\omega_{c}=50;$
for MPS simulations $N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$.

At $g=0.6$, the three approaches give the correct description, because the
regime is still perturbative in $g$ and the effect of the bath negligible. At
$g=1.0$, there is a difference in the amplitude of the oscillations of the
cavity position obtained from HEM. Indeed, this solution is strictly valid
only for small couplings $g$, unable to significantly influence the qubit
dynamics. For the fidelity ($(c)$ and $(d)$ in Fig. 7) we compare MPS and
Lindblad methods. For $g=0.6$, they overlap almost completely with each other,
ever above the $95\%$ for $5$ qubit periods. This robust fidelity together
with the two peaks of the FFT in Fig. 6 shows how it is possible to still use
the HEM approach for oscillator coordinates and energy to reconstruct the
qubit state in time. Instead, for $g=1.0$ the system ends up in the limiting
case when the qubit state reconstruction is feasible for both LME and MPS
($F_{q}(5T\Delta)\approx 90\%$). Therefore, it is possible through the
resonator to read the qubit state which is still close to the free behaviour.

Figure 8:  Cavity number (in panels $(a)$ and $(b)$), cavity position (in
panels $(c)$ and $(d)$) and qubit fidelity (in panels $(e)$ and $(f)$) as a
function of time (in units of $1/\Delta$) at $g=0.6$ and $g=1.0$,
respectively, for $\alpha=0.1$ by using different methods: MPS simulations
(blue solid line), the LME solution (red dashed line) and the HEM approach
(green dotted line), improved by the fitting of the spin observables as in
Eqs. (19-20). Parameters:
$\Delta=1,\omega_{0}=10\Delta,g=\\{0.6\Delta,1.0\Delta\\},\alpha=0.1,\omega_{c}=50;$
for MPS simulations $N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$.

For $\alpha=0.1$, we find interesting behaviours with increasing $\tilde{g}$.
We start again from $g=0.6$, which has been considered both in Fig. 6 and 7.
We recall that, for $\alpha=0.1$, the standard readout is not effective to get
the spin component along the z-axis. Therefore, we investigate if a state
reconstruction can still be achieved. As shown in Fig. 8, the time evolution
of the resonator number agrees with that predicted from the HEM improved by
the fitting of the spin z-component (Eq. (21)), while the LME remains
perturbative and does not significantly change the state from the initial
vacuum. For the position, instead, the agreement is good at $g=0.6$ between
LME and HEM, while there is a slightly different amplitude in the MPS
position. In fact, LME fails at high coupling to the bath while HEM at high
renormalization of the qubit dynamics. Thus, these perturbative approaches are
unable to properly transfer the effect of the environment to the qubit through
the interaction with the cavity. To understand this point, one can consider
the mapped model in Eq. (4) where the effective spectral density is Ohmic at
low frequencies with coupling proportional to $g^{2}\alpha/\omega_{0}^{2}$.
For $g=1.0$, we observe also a discrepancy between LME and HEM, due to the
increased $g$. Moreover, they are phase-shifted with respect to the MPS due to
the fact that a high value of the coupling in the mapped model Zueco and
García-Ripoll (2019) means a qubit strongly coupled to the effective bath.
The panels $(e)$ and $(f)$ of Fig. 8 show qubit fidelity for $\alpha=0.1$,
comparing MPS and LME methods. For $g=0.6$, LME gives an incorrect behaviour,
while the MPS simulations provide the limiting case $F_{q}(5T\Delta)\approx
90\%$ in which the state reconstruction can be performed. Indeed, one can
still determine x, y and z spin components from the behaviour of oscillator
position, momentum and number as predicted by HEM solutions. In panel Fig.
8$(f)$, for $g=1.0$, the difference between LME and MPS becomes large,
indicative of the fact that in this regime LME fails, while for MPS
$F_{q}(5T\Delta)\approx 60\%$. In fact, LME would provide still a perturbative
description, with an acceptable fidelity; MPS, however, shows a fidelity
almost linearly decreasing so that it results $F_{q}(5T\Delta)\approx 60\%$.
Summarizing, we have found that for low $\alpha$ in the strong coupling regime
both the typical readout for an up state and for a generic state can be
achieved making the qubit state reconstruction possible. For increasing
$\alpha$, anyway, the readout is unachievable, while for $g\leq 0.6$ a state
reconstruction is still attainable (see Fig. 12 and Table 1 in the concluding
section IV).
In analogy with weak coupling regime discussed in the previous subsection, we
find that the effect of the decoherence helps the qubit state reconstruction
by reducing the initial transient and the oscillator is in a moving momentum
squeezed state as pointed out in II.1.1.

###  III.3 Ultra-strong coupling regime - interactions prevent state
reconstruction

Lastly, with the qubit frequency $\Delta$ taken as unity and
$\omega_{0}=10\Delta$, by focusing on $g/\Delta=5$ and tuning
$\alpha\in\\{0.01,0.1\\}$ such that $g/\Delta\geq 0.1\omega_{0}/\Delta$ (see
Fig. 12), the quantum Rabi system is in the ultra-strong coupling regime.

As in the previous regime, if one chooses $\ket{0}\otimes\ket{\uparrow}$ as
initial state, for low values of $\alpha$ it is possible to follow the peaks
of the FFT of the mean resonator number. In Fig. 9 we see, as expected, that,
for $\alpha=0.01$, a peak is prominent, while, for $\alpha\geq 0.1$, the FFT
is almost flat. However, even for $\alpha=0.01$, the frequency corresponding
to the peak is no more in perfect agreement with the Bloch-Siegert shift given
in Eq. (35) because the system is in the non-perturbative region of ultra-
strong coupling regime Rossatto  _et al._ (2017)
($g\approx\Delta+\omega_{0}$). Moreover, also the peak corresponding to the
down contribution is emerging. Therefore, the qubit readout is not feasible
because of the marked back-action of the resonator. Below we will show that,
in this coupling regime, qubit and resonator are indeed strongly entangled.

Figure 9:  FFT of the mean cavity number $\tilde{n}$ computed through MPS
simulations as a function of the frequency (in units of $\Delta$) for two
values of $\alpha$. Parameters:
$\Delta=1,\omega_{0}=10\Delta,g=5\Delta,\alpha=\\{0.01,0.1\\},\omega_{c}=50;$
for MPS simulations $N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$. Figure 10:
Cavity position (in panels $(a)$ and $(b)$) and qubit fidelity (in panels
$(c)$ and $(d)$) as a function of time (in units of $1/\Delta$) for
$\alpha=0.01$ and $\alpha=0.1$, respectively, using three methods: MPS
simulations (blue solid line), LME solutions (red dashed line) and HEM
approach (green dotted line), improved by the fitting of spin observables as
in Eqs. (19-20). Parameters:
$\Delta=1,\omega_{0}=10\Delta,g=5\Delta,\alpha=\\{0.01,0.1\\},\omega_{c}=50;$
for MPS simulations $N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$.

For the initial state $\ket{0}\otimes\ket{\nearrow}$ and the bath empty,
regardless of the value of $\alpha$, qubit state reconstruction is no more
feasible (see Fig. 10). In particular, the agreement between different
computational methods is good for $\alpha=0.01$ at short times, while there is
a dephasing in LME solution due to the fact that a high value of the coupling
in the mapped model Zueco and García-Ripoll (2019); De Filippis  _et al._
(2022) means a qubit strongly coupled to the effective bath. For larger
$\alpha$ also HEM fails because the fit is not enough to describe the qubit
dynamics. In Fig. 10$(c)$ and Fig. 10$(d)$ we compare the fidelity for the MPS
and LME approaches. For $\alpha=0.01$, the results are very similar,
oscillating around $0.5$ due to the growing entanglement, while for
$\alpha=0.1$ they have an exponential decaying behaviour, with an
overmodulation for MPS. Anyway, during the first $5$ qubit periods, both MPS
and LME solutions give a minimum fidelity close to $20\%$. We emphasize that a
state reconstruction is still possible, at least for low $\alpha$ values, but
what we are reconstructing is the qubit state as a part of the coupled system.
As a consequence, it is very entangled, so its evolution is far from the free
one and the fidelity is low (see Fig. 12 and Table 1 in the concluding section
IV).

This set of parameters allows us to observe the moving of the squeezed state
(being the amplitude of the displacement proportional to $\tilde{g}$) (Fig.
11). The position of the resonator oscillates around zero as shown in Fig. 10,
while its momentum remains close to the zero. The cavity starts in the vacuum
and evolves in a momentum squeezed state whose center oscillates with damping
around zero until stopping at the origin of the quantum phase space, as
underlined in II.1.1.

Figure 11:  Reduced density matrix of the cavity computed through MPS
simulations at $8$ times $t\Delta$ during its evolution. From this, we obtain
the Wigner quasi-probability distribution by using the proper function of the
package QuTip Johansson  _et al._ (2012). Parameters:
$\Delta=1,\omega_{0}=10\Delta,g=5\Delta,\alpha=0.1,\omega_{c}=50;$ for MPS
simulations $N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$.

##  IV Discussion and conclusions

In the present work, we have provided an exhaustive description of a method
for reconstructing the qubit state through the dynamics of its coupled
resonator. In particular, we have tested the functional relationships derived
from the perturbative (in $g$) solutions of HEM, by simulating the quantum
Rabi model with MPS numerical method. We have compared it also with a
numerical solution of the LME, which is perturbative in the oscillator-bath
coupling and unable to describe the non-Markovian effects of the bath. We have
found the critical values of the coupling to the bath for the usual case of
the qubit readout (up state) and for a generic qubit state since we have
extended the procedure to an open system. Moreover, we have implemented full
state reconstruction of the qubit, by computing the time evolution of the
cavity observables. We have found the parameter regimes where the procedure
does not work, taking the qubit fidelity with respect to its free evolution as
a measure of the goodness of the state reconstruction. Table 1 shows the qubit
fidelity computed at time $t_{\rm{final}}=35/\Delta$ for weak and strong
coupling regimes, while for the ultra-strong regime the minimum value is
reached at time $t=7/\Delta$ for the MPS method and $t=19.5/\Delta$ for the
LME one.

Table 1: Minimum qubit fidelity $F_{q}(t)$ (%) for the three regimes analysed in the paper  Method | Weak | Strong | Ultra-strong
---|---|---|---
| $g=0.1,\alpha=0.05$ | $g=0.6,\alpha=0.1$ | $g=5.0,\alpha=0.1$
MPS | 99.8 | 61.2 | 20.4
Lindblad | 99.9 | 88.0 | 18.7
Figure 12: Phase diagram of the quantum Rabi model for varying $\alpha$ and
$g$. On the background there are the three regions of parameters analysed in
the paper: weak coupling (light grey), strong coupling (grey) and ultra-strong
coupling (black). The two methods under study are shown using different
patterns: qubit readout with blue “+” and qubit state reconstruction with red
“//”. The grey vertical line at $\alpha=2/(10\pi)$ underlines that above this
value the regime is again the strong coupling one because
$\Gamma/\Delta>0.1\omega_{0}/\Delta$. Parameters:
$\Delta=1,\omega_{0}=10\Delta,g\in[0\Delta,2\Delta],\alpha\in[0,2/(10\pi)],\omega_{c}=50,N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$.

In Fig. 12 we show the validity of the two methods analysed in detail in the
paper, that are the standard qubit readout and the full qubit state
reconstruction via oscillator dynamics in the three regimes: weak, strong and
ultra-strong coupling. We stop the plot at $\alpha=2/(10\pi)$ that corresponds
to $\Gamma/\Delta=0.1\omega_{0}/\Delta$. The latter indicates that the
estimation of the perturbative decay rate $\Gamma$ for the definition of the
parameter regions no longer holds. Hence, above this value of $\alpha$ the
regime is again the strong coupling one. We emphasize that for internal
coupling $g\leq 0.8\Delta$ the state reconstruction is possible for coupling
strengths $\alpha$ bigger than those necessary for the standard readout.
Furthermore, for increasing internal couplings, even if the LME solution would
give a reliable qubit fidelity, regardless of the coupling to the bath, the
MPS simulations have shown that it is already much lesser than the acceptance
threshold of $90\%$. Nevertheless, the state reconstruction is still possible,
but the interactions have modified significantly the qubit dynamics with
respect to the free one. In the ultra-strong regime the fidelity oscillates
around the value of $50\%$, that can be easily interpreted by looking at the
qubit entanglement in time.

Figure 13: Qubit Von Neumann entropy ($S_{q}=-\Tr\\{\rho_{q}\log\rho_{q}\\}$)
and purity ($P_{q}=1-\Tr\\{\rho_{q}^{2}\\}$) computed through MPS simulations
in panels $(a)$ and $(b)$ respectively, for $g=0.1,\,\alpha=0.05$, i.e. weak
coupling (blue solid line), $g=0.6,\,\alpha=0.1$ i.e. strong coupling (red
dashed line) and $g=5.0,\,\alpha=0.1$ i.e. ultra-strong coupling (green dash-
dotted line). Parameters:
$\Delta=1,\omega_{0}=10\Delta,g=\\{0.1\Delta,0.6\Delta,5\Delta\\},\alpha=\\{0.05,0.1\\},\omega_{c}=50;$
for MPS simulations $N_{o}=20,N=500,N_{ph}=3,D_{\rm{max}}=20$.

In Fig. 13, the qubit Von Neumann entropy
($S_{q}=-\Tr\\{\rho_{q}\log\rho_{q}\\}$) and the purity
($P_{q}=1-\Tr\\{\rho_{q}^{2}\\}$) are displayed for the three different
regimes of parameters examined in the paper. We clearly observe how they
increase towards the maximum values (dotted black lines $S_{q}^{max}=\log 2$
and $P_{q}^{max}=1/2$, respectively) in the ultra-strong coupling regime.
Actually, for the coupled system qubit-resonator, entanglement becomes bigger
and bigger. This is the main reason why the qubit fidelity with respect to its
free evoution becomes smaller and smaller.

Finally, we point out how the method utilised in the work to reconstruct the
qubit state is not restricted to the anti-adiabatic regime
($\Delta\ll\tilde{\omega}_{0}$), but the description is valid for a wider
range of parameters. In Fig. 14, we depict the forbidden region of the set of
three parameters ($\tilde{\omega}_{0}/\Delta,\alpha,\tilde{g}/\Delta$), while
$\Delta=1$, by limiting the $\tilde{g}/\Delta$ values to $2$ and the $\alpha$
ones to $1$ (very high coupling strengths) and $\tilde{\omega}_{0}/\Delta$ to
$10$, of the order of the value chosen in the paper. This region is obtained
by performing the union among the inequalities that define the conditions for
neglecting the spurious terms and obtaining the analytical particular
solutions for $\tilde{x}(t)$ and $\tilde{p}(t)$ (Eqs. 13-14). For the mean
position, at fixed initial conditions $\langle\sigma_{x}(0)\rangle$ and
$\langle\sigma_{y}(0)\rangle$ the inequalities to be both verified read:

$\displaystyle\frac{2\pi\alpha\Delta\tilde{\omega}_{0}^{2}\tilde{g}\langle\sigma_{x}(0)\rangle}{\left(\tilde{\omega}_{0}^{2}-\Delta^{2}\right)^{2}+\pi^{2}\alpha^{2}\tilde{\omega}_{0}^{2}\Delta^{2}}$
$\displaystyle\ll 1$ (36)
$\displaystyle\frac{2\pi\alpha\Delta\tilde{\omega}_{0}^{2}\tilde{g}\langle\sigma_{y}(0)\rangle}{\left(\tilde{\omega}_{0}^{2}-\Delta^{2}\right)^{2}+\pi^{2}\alpha^{2}\tilde{\omega}_{0}^{2}\Delta^{2}}$
$\displaystyle\ll 1,$ (37)

and similarly those for the momentum are:

$\displaystyle\frac{2\pi\alpha\Delta^{2}\tilde{\omega}_{0}\tilde{g}\langle\sigma_{x}(0)\rangle}{\left(\tilde{\omega}_{0}^{2}-\Delta^{2}\right)^{2}+\pi^{2}\alpha^{2}\tilde{\omega}_{0}^{2}\Delta^{2}}$
$\displaystyle\ll 1$ (38)
$\displaystyle\frac{2\pi\alpha\Delta^{2}\tilde{\omega}_{0}\tilde{g}\langle\sigma_{y}(0)\rangle}{\left(\tilde{\omega}_{0}^{2}-\Delta^{2}\right)^{2}+\pi^{2}\alpha^{2}\tilde{\omega}_{0}^{2}\Delta^{2}}$
$\displaystyle\ll 1.$ (39)

Figure 14: Our method is not valid for the values of
($\tilde{\omega}_{0}/\Delta,\alpha,\tilde{g}/\Delta$) in the red 3D region, by
taking $\Delta$ as unity. The region is the result of the union among the
inequalities (Eqs. 36-39) that define the conditions for obtaining (Eqs.
13-14) for the oscillator dynamical variables.

At this stage, as recalled in the introduction I, we can thoroughly discuss
the results of the experiments in Ref. Ficheux  _et al._ , 2018. Indeed, in
these experiments, a full qubit state reconstruction is realised by reading
the quadratures of a coupled resonator via homodyne and heterodyne
measurements. The unraveling of quantum trajectories has been used in order to
follow the qubit dynamics. The method is also known as the method of quantum
jumps, exactly the same as the Lindblad operators Breuer  _et al._ (2002). The
resulting master equation is then an appropriately weighted stochastic average
over all of the different times at which the jumps could occur, and all of the
different types of jumps that can occur Daley (2014). Hence, this approach is
equivalent to the Lindblad approach for the qubit in the limit of a
sufficiently large number of experimental realizations. It is remarkable that
a full quantum state reconstruction can be obtained at any time by raw
averaging measurement outcomes of many realizations of a single experiment
despite the incompatibility of the three spin components of the qubit.
In our work the idea of how performing the qubit state reconstruction is
similar to that in Ref. Ficheux  _et al._ , 2018, that is to read the
resonator dynamics. In fact, the functional relationships linking the qubit
dynamics to the resonator evolution, that we have derived from HEM, bear some
resemblance to those used in Ref. Ficheux  _et al._ , 2018 to interpret
measurements. The main difference is that in Ref. Ficheux  _et al._ , 2018 the
bath has not been explicitly considered, while, in this paper, we have
exploited its effects on the qubit state reconstruction. Moreover, the LME
used in this work takes qubit and resonator on the same footing. Finally, we
have introduced dissipation and decoherence through the presence of an Ohmic
bath, without describing a measurement apparatus. A generalization of our
study to a more realistic model, including the measurement process
description, deserves future investigations Breuer  _et al._ (2002); Walls
_et al._ (1985).

The study carried out in the present work can be extended to the systems where
there are more qubits and more resonators. However, while these
generalizations are straightforward in HME and LME approaches, they would be
clearly more challenging within a MPS framework. Such architectures are
typical of the recent noisy intermediate-scale quantum processors.
Furthermore, the resonator could be investigated in more detail, since the
interest about the bosonic behaviour, specially in circuit QED, is growing. In
fact, quantum information can be encoded into subspaces of a bosonic
superconducting cavity mode with long coherence time Ma  _et al._ (2021).

## References

  * Gu  _et al._ (2017a) X. Gu, A. F. Kockum, A. Miranowicz, Y.-x. Liu, and F. Nori, Phys. Rep. 718, 1 (2017a).
  * Rabi (1936) I. Rabi, Phys. Rev. 49, 324 (1936).
  * Jaynes and Cummings (1963) E. T. Jaynes and F. W. Cummings, Proc. IEEE 51, 89 (1963).
  * Zueco  _et al._ (2009) D. Zueco, G. M. Reuther, S. Kohler, and P. Hänggi, Phys. Rev. A 80, 033846 (2009).
  * Wallraff  _et al._ (2004) A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, Nat. 431, 162 (2004).
  * 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._ (2019) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, RMP 91, 025005 (2019).
  * Perroni  _et al._ (2016) C. Perroni, D. Ninno, and V. Cataudella, J. Condens. Matter Phys. 28, 373001 (2016).
  * Nocera  _et al._ (2016) A. Nocera, C. A. Perroni, V. M. Ramaglia, and V. Cataudella, Beilstein J. Nanotechnol. 7, 439 (2016).
  * Perroni and Benenti (2021) C. A. Perroni and G. Benenti, Advances in Thermoelectricity: Foundational Issues, Materials and Nanotechnology 207, 115 (2021).
  * Zagoskin (2011) A. M. Zagoskin,  _Quantum Engineering: Theory and Design of Quantum Coherent Structures_ (Cambridge University Press, 2011).
  * Hacohen-Gourgy  _et al._ (2016) S. Hacohen-Gourgy, L. S. Martin, E. Flurin, V. V. Ramasesh, K. B. Whaley, and I. Siddiqi, Nat. 538, 491 (2016).
  * Campagne-Ibarcq  _et al._ (2016) P. Campagne-Ibarcq, P. Six, L. Bretheau, A. Sarlette, M. Mirrahimi, P. Rouchon, and B. Huard, Phys. Rev. X 6, 011002 (2016).
  * Ficheux  _et al._ (2018) Q. Ficheux, S. Jezouin, Z. Leghtas, and B. Huard, Nat. Commun. 9, 1 (2018).
  * Breuer  _et al._ (2016) H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, RMP 88, 021002 (2016).
  * Weimer  _et al._ (2021) H. Weimer, A. Kshetrimayum, and R. Orús, RMP 93, 015008 (2021).
  * Zueco and García-Ripoll (2019) D. Zueco and J. García-Ripoll, Phys. Rev. A 99, 013807 (2019).
  * Breuer  _et al._ (2002) H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Zaletel  _et al._ (2015) M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B 91, 165112 (2015).
  * Haegeman  _et al._ (2011) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
  * Paeckel  _et al._ (2019) S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Ann. Phys. 411, 167998 (2019).
  * Fishman  _et al._ (2020) M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor software library for tensor network calculations (2020), arXiv:2007.14822 .
  * Zhang  _et al._ (1998) C. Zhang, E. Jeckelmann, and S. R. White, Phys. Rev. Lett. 80, 2661 (1998).
  * Bursill (1999) R. J. Bursill, Phys. Rev. B 60, 1643 (1999).
  * Friedman (2000) B. Friedman, Phys. Rev. B 61, 6701 (2000).
  * Wong and Chen (2008) H. Wong and Z.-D. Chen, Phys. Rev. B 77, 174305 (2008).
  * Brockt  _et al._ (2015) C. Brockt, F. Dorfner, L. Vidmar, F. Heidrich-Meisner, and E. Jeckelmann, Phys. Rev. B 92, 241106 (2015).
  * Jansen  _et al._ (2022) D. Jansen, J. Bonča, and F. Heidrich-Meisner, arXiv preprint arXiv:2206.00985 (2022).
  * De Filippis  _et al._ (2020) G. De Filippis, A. de Candia, L. Cangemi, M. Sassetti, R. Fazio, and V. Cataudella, Phys. Rev. B 101, 180408 (2020).
  * De Filippis  _et al._ (2021) G. De Filippis, A. de Candia, A. Mishchenko, L. Cangemi, A. Nocera, P. Mishchenko, M. Sassetti, R. Fazio, N. Nagaosa, and V. Cataudella, Phys. Rev. B 104, L060410 (2021).
  * Liu  _et al._ (2022) K. T. Liu, D. N. Beratan, and P. Zhang, Phys. Rev. A 105, 032406 (2022).
  * Gu  _et al._ (2017b) X. Gu, A. F. Kockum, A. Miranowicz, Y.-x. Liu, and F. Nori, Phys. Rep. 718, 1 (2017b).
  * Bloch and Siegert (1940) F. Bloch and A. Siegert, Phys. Rev. 57, 522 (1940).
  * Klimov and Chumakov (2009) A. B. Klimov and S. M. Chumakov,  _A group-theoretical approach to quantum optics: models of atom-field interactions_ (John Wiley & Sons, 2009).
  * Rossatto  _et al._ (2017) D. Z. Rossatto, C. J. Villas-Bôas, M. Sanz, and E. Solano, Phys. Rev. A 96, 013849 (2017).
  * De Filippis  _et al._ (2022) G. De Filippis, A. de Candia, G. Di Bello, C. Perroni, L. Cangemi, A. Nocera, M. Sassetti, R. Fazio, and V. Cataudella, arXiv preprint arXiv:2205.11555 (2022).
  * Johansson  _et al._ (2012) J. R. Johansson, P. D. Nation, and F. Nori, Comput. Phys. Commun. 183, 1760 (2012).
  * Daley (2014) A. J. Daley, Adv. Phys. 63, 77 (2014).
  * Walls  _et al._ (1985) D. Walls, M. Collet, and G. Milburn, Phys. Rev. D 32, 3208 (1985).
  * Ma  _et al._ (2021) W.-L. Ma, S. Puri, R. J. Schoelkopf, M. H. Devoret, S. M. Girvin, and L. Jiang, Sci. Bull. 66, 1789 (2021).
  * Lindblad (1976) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * Manzano and Hurtado (2018) D. Manzano and P. Hurtado, Adv. Phys. 67, 1 (2018).
  * González  _et al._ (2017) J. O. González, L. A. Correa, G. Nocerino, J. P. Palao, D. Alonso, and G. Adesso, OSID 24, 1740010 (2017).
  * Farina  _et al._ (2020) D. Farina, G. De Filippis, V. Cataudella, M. Polini, and V. Giovannetti, Phys. Rev. A 102, 052208 (2020).
  * Pang (1999) T. Pang, An introduction to computational physics (1999).

###### Acknowledgements.

We acknowledge interesting discussions with R. Fazio. A.N. acknowledges
support from the Max Planck-UBC-UTokyo Center for Quantum Materials and the
Canada First Research Excellence Fund (CFREF) Quantum Materials and Future
Technologies Program of the Stewart Blusson Quantum Matter Institute (SBQMI).

##  Appendix A Lindblad master equation

The most straightforward approach to the dynamics of an open quantum system is
based on the solution of a quantum master equation which generalises the
Liouville - Von Neumann equation for the reduced density operator of the
system of interest. In fact, when the system is isolated, described by a time-
independent Hamiltonian, or closed, described by a time-dependent Hamiltonian
(e.g. a driving field), its dynamics is fully characterised by the Liouville -
Von Neumann equation:

$\frac{d\rho_{I}(t)}{dt}=-\frac{i}{\hbar}[\mathcal{H}_{I}(t),\rho_{I}(t)]=\mathcal{L}\rho_{I}(t),$
(40)

written in the interaction picture. The subscript $I$ stays for “interaction”
and the whole system Hamiltonian is
$\mathcal{H}(t)=\mathcal{H}_{0}+\mathcal{H}_{I}(t)$, with $\mathcal{H}_{0}$
the unperturbed energy of the separated quantum systems and
$\mathcal{H}_{I}(t)$ the interaction between them. Moreover, $\mathcal{L}$ is
the Liouville superoperator that acting on the density matrix gives its time
evolution.
The famous Lindblad master equation Lindblad (1976) (LME) is a Markovian
quantum master equation that has the advantage of being easy to solve
numerically through suitable tools, and analytically in some special
situations. An open quantum system is said to be Markovian when its behaviour
at a given time is independent of its behaviour in the past. The main drawback
of Markovian framework is the lack of the back-action upon the system of
interest due to the environment, that is indeed neglected due to the
hypothesis of the so-called Born-Markov approximation, valid for systems
weakly coupled to the environment Manzano and Hurtado (2018). In fact,
according to this approximation, the environment relaxes before the system
changes its state, $\rho(t)=\rho_{S}(t)\otimes\rho_{B}(0)$ so that it cannot
act on the system (Born approximation). Furthermore, Markov approximation,
i.e. $\tau_{R}\gg\tau_{B}$, where $\tau_{R}$ is the relaxation time of the
system $S$ and $\tau_{B}$ is the correlations time of the bath, means that the
dynamics over a time of the order of $\tau_{B}$ cannot be resolved. The Born-
Markov approximation does not warrant the complete positivity of the dynamical
map. Therefore, we need to perform one further approximation, the secular or
rotating wave approximation (RWA). If $\tau_{S}$ is the time scale of the
system evolution and $\tau_{S}\ll\tau_{R}$, the non-secular terms may be
neglected, since they oscillate very rapidly during the time over which
$\rho_{S}$ varies appreciably. If we separate the Hermitian and non-Hermitian
parts of the dynamics of the system, and we return to the Schrödinger picture
diagonalising the matrix formed by the rate coefficients, we obtain the LME
for $\rho_{S}(t)$:

$\begin{split}&\frac{d\rho_{S}}{dt}\equiv\mathcal{L}\rho_{S}=-\frac{i}{\hbar}[\mathcal{H}_{S}+\mathcal{H}_{LS},\rho_{S}]\\\
&+\sum_{\omega,k}\gamma_{k}(\omega)\left(L_{k}(\omega)\rho_{S}L_{k}^{\dagger}(\omega)-\frac{1}{2}\left\\{L_{k}^{\dagger}(\omega)L_{k}(\omega),\rho_{S}\right\\}\right),\end{split}$
(41)

where the first term is the unitary evolution with $\mathcal{H}_{LS}$ the Lamb
and Stark shift Hamiltonian, whose role is to renormalize the system energy
levels due to the interaction with the environment. The sum is the Dissipator
with rates $\gamma_{k}$ and $L_{k}$ the Lindblad or jump operators.

When dealing with multipartite systems as in our Hamiltonian given in Eq. (3),
the key point is how to eigendecompose the system observables occurring in the
interaction Hamiltonian. One way (Local) consists in neglecting the “internal”
interaction among the parties González  _et al._ (2017). In this case the
eigenstates are factorised and one can only add the dissipators pertaining the
different operators in order to describe the dissipation. The other one
(Global) instead is based on the use of the correct eigenstates of the
Hamiltonian of the system surrounded by the bath. Moreover, it is known that
the global approach provides an accurate description of the system evolution
at long timescales Farina  _et al._ (2020). Thus, in the paper we follow this
approach, because we aim to properly describe the interaction between qubit
and oscillator up to times of the order of some qubit periods. Hence, we
diagonalise the Hamiltonian of qubit plus oscillator system (Eq. (3) without
the bath) at each time step $dt$, we eigendecompose the system observables and
we evaluate the Lindbladian $\mathcal{L}$. Eq. (41) for the reduced density
matrix of the bipartite system $\rho_{S}(t)$, vector in the Liouville space,
is solved numerically through the very well-known fourth-order Runge-Kutta
algorithm, derived from two different Taylor expansions of the dynamical
variable $\rho_{S}$ and its derivative Pang (1999):

$\rho_{S}(t+dt)=\rho_{S}(t)+\frac{1}{6}(c_{1}+2c_{2}+2c_{3}+c_{4}),$ (42a)
$\displaystyle c_{1}=\,\,$ $\displaystyle
dt\,\mathcal{L}\left(\rho_{S}(t)\right),$ $\displaystyle c_{2}=\,\,$
$\displaystyle dt\,\mathcal{L}\left(\rho_{S}(t)+\frac{c_{1}}{2}\right),$
$\displaystyle c_{3}=\,\,$ $\displaystyle
dt\,\mathcal{L}\left(\rho_{S}(t)+\frac{c_{2}}{2}\right),$ $\displaystyle
c_{4}=\,\,$ $\displaystyle dt\,\mathcal{L}\left(\rho_{S}(t)+c_{3}\right).$
(42b)

We study the solution of LME in order to compare a Markovian description with
other analytical and numerical methods, such as the MPS approach introduced in
Sec. II.2.