# Unconventional mechanism of virtual-state population through dissipation

Alejandro Vivas-Viaña  Departamento de Física Teórica de la Materia Condensada
and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid,
28049 Madrid, Spain  Alejandro González-Tudela  Institute of Fundamental
Physics IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain  Carlos Sánchez
Muñoz<EMAIL_ADDRESS>Departamento de Física Teórica de la
Materia Condensada and Condensed Matter Physics Center (IFIMAC), Universidad
Autónoma de Madrid, 28049 Madrid, Spain

###### Abstract

Virtual states are a central concept in quantum mechanics. By definition, the
probability of finding a quantum system in a virtual state should be
vanishingly small at all times. In contrast to this notion, we report a
phenomenon occurring in open quantum systems by which virtual states can
acquire a sizable population in the long time limit, even if they are not
directly coupled to any dissipative channel. This means that the situation
where the virtual state remains unpopulated can be metastable. We describe
this effect by introducing a two-step adiabiatic elimination method, that we
termed hierarchical adiabatic elimination, which allows one to obtain
analytical expressions of the timescale of metastability in general open
quantum systems. We show how these results can be relevant for practical
questions such as the generation of stable and metastable entangled states in
dissipative systems of interacting qubits.

Introduction.— The concept of virtual in quantum mechanics is of paramount
importance, e.g. in the context of virtual transitions between coherently
unconnected states [1, 2, 3] or in the description of scattering processes in
QFT where interactions are mediated by virtual particles [4]. In situations
where strongly off-resonant “virtual” states mediate interactions between
quasi-resonant “real” states, an adiabatic elimination over the fast degrees
of freedom—the virtual ones—allows one to reduce the dimensionality of the
problem and obtain an effective description of the slow degrees of freedom,
i.e. the real states. This technique of adiabatic elimination, which can be
formulated in several alternatives ways—e.g., the Schrieffer-Wolff
transformation [2]—is ubiquitous in the description and design of quantum
phenomena, e.g. quantum optical applications in atomic physics [5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15] or exotic dynamics in the ultrastrong coupling regime
of cavity QED [16, 17, 18]. A significant effort has been made to establish
the mathematical foundations of this technique [19, 20, 21] and its extension
to dissipative contexts for its application in open quantum systems [22, 23,
24, 25].

The fundamental underlying assumption for the adiabatic elimination of a
virtual state is that the coupling between the real subspace $\mathcal{H}_{R}$
and the virtual subspace $\mathcal{H}_{V}$ is perturbative; i.e., the coherent
coupling rate is much smaller than the energy difference between subspaces; as
a result, one can obtain an effective Hamiltonian acting only in
$\mathcal{H}_{R}$. Consequently, when this approximation applies, any initial
state in $\mathcal{H}_{R}$ will remain within that subspace, and
$\mathcal{H}_{V}$ will not be populated. In this work, we show that the
situation can be radically different when there is also dissipative dynamics,
even if dissipative process only take place within $\mathcal{H}_{R}$. We
unveil an unconventional mechanism by which, in the long time limit, virtual
states acquire a sizable occupation probability, comparable to that of the
real states. These findings can have great importance in the understanding and
engineering of interactions between quantum systems in driven-dissipative
contexts [26, 27, 28, 29].

In order to study this phenomenon, we start discussing what is arguably the
simplest scenario that can be described in terms of virtual states [see
Fig.1(a)]: two quasi-resonant “real” states, effectively interacting through
the mediation of a third, strongly off-resonant “virtual” state. Crucially, we
enable a spontaneous decay between the real states, which can be provided, for
instance, by the coupling to a surrounding environment in a Markovian regime.
Contrary to the familiar intuition, the situation in which the virtual state
remains “virtual” is, in this case, only metastable [30, 31], and, in the long
time limit, the system eventually relaxes to a stationary state where the
virtual state has a sizable population.

This process of _de-virtualization_ occurs through an unconventional mechanism
of population enabled by dissipation. In this work, we introduce a technique
of _hierarchical adiabatic elimination_ to obtain analytical approximations of
the time-dependent elements of the system density matrix and expressions for
the characteristic metastability timescales. We show how our novel technique
can be used to described metastable dynamics in different systems involving
two interacting qubits, where the phenomenon reported has strong implications
for the generation of stable and metastable entanglement via dissipation.

Model.— The first model we study consists of a three-level system
configuration, sketched in Fig. 1(a). The Hilbert space spans a basis
$\left\\{|1\rangle,|2\rangle,|V\rangle\right\\}$, where the states $|1\rangle$
and $|2\rangle$ represent two real states, and $|V\rangle$ will play the role
of a virtual state, being strongly detuned from $|1\rangle$ and $|2\rangle$.
We define lowering operators as $\hat{\sigma}_{i,j}\equiv|i\rangle\langle j|$
$(i,j={1,2,V})$. The real states are coupled to $|V\rangle$ with a coupling
rate $\Omega$, and there is an irreversible decay process within the real
subspace, with state $|2\rangle$ decaying towards $|1\rangle$ with a decay
rate $\Gamma$. This specific $\Lambda$ model could be motivated, for instance,
by the description of a quantum-optical system consisting of two interacting
qubits coherently excited at the two-photon resonance with a Rabi frequency
$\Omega$ [32, 33, 34, 35], where the ground state $|gg\rangle$ corresponds to
the real state $|1\rangle$, the doubly-excited state $|ee\rangle$ is the
excited real state $|2\rangle$, and the symmetric single-excitation state
$|S\rangle=\frac{1}{\sqrt{2}}(|eg\rangle+|ge\rangle)$—detuned from the two-
photon transition energy due to the interaction between qubits—corresponds to
the virtual state $|V\rangle$. The two-photon decay channel can be enabled,
for instance, by a cavity in resonance with the two-photon transition [36,
37]. A change to the rotating frame of the drive would directly yield the
configuration shown in Fig. 1(a).

Figure 1: a) Scheme of the system: two quasi-resonant “real” quantum states,
interacting via a third, strongly off-resonant “virtual” state. There is
spontaneous decay between the real states. (b) The system in (a) can describe
two interacting two-level systems under coherent driving at the two-photon
resonance, in the rotating frame of the drive. Decay between real states can
be engineered with a cavity.

The resulting time-independent Hamiltonian is
$\hat{H}=\hat{H}_{0}+\hat{H}_{d}$, where $\hat{H}_{0}$ is the bare Hamiltonian
($\hbar=1$),

$\hat{H}_{0}=\Delta_{1}|1\rangle\langle 1|+\Delta_{2}|2\rangle\langle
2|+\Delta_{V}|V\rangle\langle V|,$ (1)

and $\hat{H}_{d}$ is the Hamiltonian of the driving/coupling term

$\hat{H}_{d}=\Omega\left(\hat{\sigma}_{1,V}+\hat{\sigma}_{2,V}+\text{H.c}\right),$
(2)

where $\Delta_{i}$ ($i=1,2,V$) stand for the free energy parameters, where we
will assume that $\Delta_{V}\gg\Delta_{2},\Delta_{1},\Omega$ and
$\Delta_{2}\approx\Delta_{1}\approx 0$. We assume that the evolution of the
system is governed by a quantum master equation [38],

$\frac{d\hat{\rho}}{dt}=-i[\hat{H},\hat{\rho}]+\frac{\Gamma}{2}\mathcal{L}_{\hat{\sigma}_{21}}[\hat{\rho}]+\frac{\Gamma_{V}}{2}\mathcal{L}_{\hat{\sigma}_{1V}}[\hat{\rho}],$
(3)

where the Lindblad term $\mathcal{L}_{\hat{O}}\equiv
2{\hat{O}}\hat{\rho}{\hat{O}}^{\dagger}-\left\\{{\hat{O}}^{\dagger}{\hat{O}},\hat{\rho}\right\\}$
describes processes of spontaneous decay. Unless stated otherwise, we will
consider $\Gamma_{V}=0$, i.e., we assume there is only one process of
spontaneous decay, from $|2\rangle$ to $|1\rangle$ (the case $\Gamma_{V}\neq
0$ will be considered later only for comparison). The dynamics of the system
can be studied straightforwardly by numerically solving Eq. (3). Figure 2(a)
shows the occupation probability of the excited state,
$\rho_{2,2}\equiv\langle 2|\hat{\rho}|2\rangle$ and the virtual state
$\rho_{V,V}\equiv\langle V|\hat{\rho}|V\rangle$ versus time. One can clearly
appreciate the existence of two distinct relaxation timescales. Within the
first relaxation timescale ($t\sim 1/\Gamma$), the system behaves according to
the standard intuition regarding virtual states: $|V\rangle$ remains
unpopulated, mediating the interaction between $|1\rangle$ and $|2\rangle$,
which gives rise to coherent Rabi oscillations between these two states with a
two-photon Rabi frequency $\Omega_{2\mathrm{p}}=\Omega^{2}/\Delta_{V}$, damped
by spontaneous emission of rate $\Gamma$ into a stationary state. This
situation can be described simply in terms of a coherently driven two-level
system spanned by $|1\rangle$ and $|2\rangle$. This stationary regime is,
however, metastable, and in a much longer timescale, which for this particular
choice or parameters is $t\sim 10^{4}/\Gamma$, $\rho_{V,V}$ develops a
population comparable to $\rho_{2,2}$. Clear evidences of this metastable
behaviour in open quantum systems can be found in the spectrum of eigenvalues
of the Liouvillian superoperator $\mathcal{L}$ [30, 31]. All these eigenvalues
$\\{\lambda_{k},k=1,2,\ldots\\}$—ordered here by its real values, so that
$\real(\lambda_{k})\geq\real(\lambda_{k+1})$—have a negative real part, and
the eigenvalue with the largest real part is necessarily equal to zero,
$\lambda_{1}=0$, its corresponding eigenstate being the steady state of the
system. The second largest real value of the Liouvillian spectrum,
$\real(\lambda_{2})$, is the Liouvillian gap [39], and it gives the relaxation
time necessary to reach the steady state, $\tau_{2}=1/|\real(\lambda_{2})|$.
Metastability results when $\lambda_{2}$ is well separated from the rest of
eigenvalues by a second gap, so that $\real(\lambda_{3})\ll\real(\lambda_{2})$
[30] (here, we assume for simplicity that, as in the case of our model, there
is only one metastable state, rather than a manifold). Then, the system
relaxes to a metastable state in a timescale
$\tau_{3}=1/|\real(\lambda_{3})|$, which will eventually evolve into the
actual steady state in a time $\tau_{2}\gg\tau_{3}$. The system we consider
here exhibits precisely this clustering of eigenvalues characteristic of
metastability, as can be seen in Fig. 2(b) where we confirm that $\tau_{3}\sim
1/\Gamma$, and $\tau_{2}\sim 10^{4}/\Gamma$. The steady state value of the
virtual state occupation probability can be computed analytically, yielding

$\rho_{V,V}^{\mathrm{ss}}=\frac{\Omega^{2}\left(\Gamma^{2}+4\Omega^{2}\right)}{2\Omega^{2}\left(\Gamma^{2}+6\Omega^{2}\right)+\Gamma^{2}\Delta_{V}^{2}}.$
(4)

In the limit $\Omega^{2}\gg\Gamma\Delta_{V}$, this expression indeed yields a
sizable population $\rho_{V,V}^{\mathrm{ss}}\approx 1/3$, clearly establishing
that the virtual state will get populated in the long time limit.

Figure 2: (a) Dynamics of the excited real state and the virtual state. Solid
lines are numerical calculations, dashed lines are analytical predictions from
the hierarchical adiabatic elimination. Dot-dashed represents a case where
$\Gamma_{V}\neq 0$, showing the stabilization of the metastable state. (b)
Liouvillian spectrum, featuring the defining characteristic of metastability
in open quantum systems: a metastability gap between $\lambda_{2}$ and
$\lambda_{k>2}$. (c) Same as in (b), with $\Gamma_{V}=\Gamma$. In this system,
metastability is no longer present. (d) Liouvillian eigenvalues versus
$\Gamma$, for $\Gamma_{V}=0$ (solid, black) and
$\Gamma_{V}=10^{-5}\Omega_{2\mathrm{p}}$. Our analytical prediction for the
value of the Liouvillian gap $\lambda_{2}$ for $\Gamma_{V}=$ is shown in
dashed-red. If $\Gamma_{V}\neq 0$, it is seen that metastability disappears
for when $\Gamma_{c}<\Gamma_{V}$. Parameters: $\Omega/\Delta_{V}=0.01$; in
(a—c), $\Gamma/\Delta_{V}=10^{-5}$. $\Gamma_{V}$ is zero unless indicated
otherwise, in which case $\Gamma_{V}=\Gamma$.

Hierarchical Adiabatic elimination (HAE).—In order to have an estimate of the
survival time of the metastable state, it would be desirable to obtain an
analytical expression of $\lambda_{2}$. A direct analytical solution for the
time evolution of the density matrix through the diagonalization of
$\mathcal{L}$ is not readily available. Nevertheless, it is clear from our
previous discussion the existence of a hierarchy of timescales, which suggests
that a series of adiabatic elimination techniques could be applied. _(i)_ The
shortest timescale is clearly governed by Hamiltonian dynamics, evidenced by a
fast oscillatory evolution of the density matrix elements. This oscillatory
dynamics stabilizes into a steady state in a timescale of the order
$\tau_{3}\sim 1/\Gamma$. In these timescales, the population of $|V\rangle$
plays the role of a fast variable: $|V\rangle$ mediates effective interactions
within the real subspace, i.e., it plays the role of a virtual state that can
be eliminated within a purely Hamiltonian evolution. This is the first
adiabatic elimination that we will perform. _(ii)_ The longest timescale is
characterized by a very slow evolution of $|V\rangle$. In this long timescale,
the relaxation of the real variables in a time $1/\Gamma$ occurs almost
instantaneously, meaning that one can treat the real variables as the _fast_
variables in a dissipative sense, i.e., they relax quickly into a time-
dependent quasi-steady state that follows the slow evolution of $\rho_{V,V}$.
This can be described in terms of a second adiabatic elimination. Note that a
direct application of standard adiabatic elimination techniques in dissipative
context, e.g. the projection-operator method [40, 25], would directly
eliminate the virtual subspace and, therefore, it would fail to capture the
mechanism that populates this state in the long time limit. We present instead
a two-step method, that we label hierarchical adiabatic elimination.

First adiabatic elimination.—Our starting point is the set of differential
equations describing the evolution of the elements of the total density
matrix, obtained from Eq. (3) as

$\displaystyle\dot{\rho}_{V,V}$
$\displaystyle=2\Omega\imaginary\left[\rho_{1,V}-\rho_{V,2}\right],$ (5a)
$\displaystyle\dot{\rho}_{2,2}$
$\displaystyle=-\Gamma\rho_{2,2}+2\Omega\imaginary\left[\rho_{V,2}\right],$
(5b) $\displaystyle\dot{\rho}_{1,V}$
$\displaystyle=i\Delta_{V}\rho_{1,V}+i\Omega\left[1+\rho_{1,2}-2\rho_{V,V}-\rho_{2,2}\right],$
(5c) $\displaystyle\dot{\rho}_{1,2}$
$\displaystyle=-(\Gamma/2)\rho_{1,2}+i\Omega\left[\rho_{1,V}-\rho_{V,2}\right],$
(5d) $\displaystyle\dot{\rho}_{V,2}$
$\displaystyle=-\left(i\Delta_{V}+\Gamma/2\right)\rho_{V,2}-i\Omega\left[\rho_{1,2}+\rho_{2,2}-\rho_{V,V}\right].$
(5e)

Based on the assumption that $|\Delta_{V}-\Delta_{i}|\gg|\Omega|$
($i\in\\{1,2\\}$), we perform an adiabatic elimination consisting in setting
$\dot{\rho}_{1,V}=\dot{\rho}_{V,2}=0$. In the limit $\Delta_{V}\gg\Gamma$, the
resulting effective equations that govern the dynamics of the real subspace
become

$\displaystyle\dot{\rho}_{2,2}$
$\displaystyle\approx-\Gamma\rho_{2,2}-2\Omega_{2\mathrm{p}}\,\text{Im}[\rho_{1,2}]+\frac{\Gamma\,\Omega_{2\mathrm{p}}}{\Delta_{V}}\,\rho_{V,V},$
(6a) $\displaystyle\dot{\rho}_{1,2}$
$\displaystyle\approx-\Gamma/2\,\rho_{1,2}+i\Omega_{2\mathrm{p}}(2\rho_{2,2}+\rho_{V,V}-1),$
(6b)

where we defined a two-photon Rabi frequency,
$\Omega_{2\mathrm{p}}\equiv\Omega^{2}/\Delta_{V}$. These equations can be
solved considering $\rho_{V,V}$ as a time-independent parameter with a fixed
value (i.e. $\dot{\rho}_{V,V}=0$). A natural choice would be to set
$\rho_{V,V}=0$. In that case, Eqs. 6a and 6b simply describe the dynamics of
two resonant levels coupled via a second-order process with a Rabi frequency
of $\Omega_{2\mathrm{p}}$ with a standard decay; e.g., the regime of coherent
two-photon driving of the transition $|gg\rangle\leftrightarrow|ee\rangle$ in
the case of a two-atom system depicted in Fig.1(b). Such a two-level system
dynamics describes accurately the short-timescale oscillatory dynamics of Fig.
2(a), where the initial state was set to be $|1\rangle$. From now on, we focus
on the strong coupling limit $\Gamma\lesssim\Omega_{2\mathrm{p}}$; otherwise,
the system is overdamped and will basically remain in the ground state
$|1\rangle$.

_Second adiabatic elimination—_ While the usual approach when eliminating a
virtual state is to indeed assume $\rho_{V,V}=0$ for all times, we have
already seen that this approach is eventually bound to fail, since
$\rho_{V,V}$ develops a sizable population within a characteristic timescale
$\tau_{2}\gg 1/\Gamma$ which, crucially, is orders of magnitude longer than
the relaxation time of Eqs. 6a and 6b. This suggest we can make a second
adiabatic elimination based on this separation of timescales. From the first
adiabatic elimination conditions ($\dot{\rho}_{V,i}=0$) and Eq. (S1a), we can
obtain a differential equation for $\rho_{V,V}$ which is a function of itself
and the real-subspace elements, i.e.
$\dot{\rho}_{V,V}(t)=f[\rho_{V,V}(t);\rho_{1,2}(t);\rho_{2,2}(t)]$ (see
Supplemental Material for a full expression). The second adiabatic elimination
consist then in substituting $\rho_{1,2}(t)$ and $\rho_{2,2}(t)$ in that
equation by their steady state solutions of Eqs. 6a and 6b obtained for a
given $\rho_{V,V}=\rho_{V,V}(t)$, yielding a dynamical equation that only
depends on $\rho_{V,V}(t)$, i.e.
$\dot{\rho}_{V,V}(t)=f[\rho_{V,V}(t);\rho_{1,2}^{\mathrm{ss}}(\rho_{V,V}(t));\rho_{2,2}^{\mathrm{ss}}(\rho_{V,V}(t))]$.
Here, $\rho_{2,2}$ and $\rho_{1,2}$ act as fast variables that relax into a
time-dependent stationary state that follows the slow evolution of
$\rho_{V,V}$. Solving this differential equation, one obtains
${\rho_{V,V}(t)\approx\rho_{V,V}^{\mathrm{ss}}\left(1-e^{-\Gamma_{c}t}\right)}$,
where we have defined the relaxation rate

$\Gamma_{c}\approx\frac{3\Gamma\Omega^{2}}{2\Delta_{V}^{2}},$ (7)

obtained under the assumption $\Omega_{2\mathrm{p}}\gg\Gamma$ (a full, more
cumbersome expression that does not require that assumption is provided in the
Supplementary Material). Equation (7) is the desired expression that gives us
the survival time of the metastable regime, and thus, it must correspond to
the Liouvillian gap, $\Gamma_{c}=|\real(\lambda_{2})|$. We have checked that
this is indeed the case in Fig. 2(d), which depicts the spectrum of
eigenvalues of $\mathcal{L}$ as a function of $\Gamma$, showing a perfect
match between our analytical expression of $\Gamma_{c}$ and $\lambda_{2}$. We
also note the perfect matching of the analytical solutions of the dynamics in
Fig. 2(a) (full expressions in Supplemental Material).

Mechanism of population.—We will now provide some insights into the
dissipative mechanism that results in the population of the virtual state. In
order to do so, we perform an analysis from the perspective of quantum
trajectories using the method of quantum jumps [41, 42, 43]. Inspection of
individual trajectories—see e.g. the example of Fig. 3(a)— reveals that the
virtual state gets populated through non-Hermitian evolution between quantum
jumps—the effect of a jump is, in fact, to strongly decrease the population of
the virtual state—. This can be understood if one considers the information
about the system leaked to the environment during a time interval with no
jumps [44]. When no jump occurs, the system is logically more likely to be in
a state that cannot emit, i.e., either $|1\rangle$ or $|V\rangle$. However,
since $|1\rangle$ is resonantly coupled to $|2\rangle$ to second-order in
perturbation theory, a system in $|1\rangle$ will eventually evolve into
$|2\rangle$ and lead to a jump. In other words, as a period without a jump
becomes longer, $|V\rangle$ becomes the most likely state, and the system is
updated accordingly, increasing its population. This purely-dissipative
mechanism will slowly accumulate over time, explaining the population buildup
of the virtual state in our system. This intuition is further confirmed by
computing a conditional density matrix for the particular trajectory in which
no jumps occur at all, see Fig. 3(b). In this particular scenario (whose
probability naturally decreases over time), the population of the virtual
state saturates to its maximum possible value.

Figure 3: (a) Population of the virtual state through non-Hermitian evolution
between quantum jumps. (b) Conditional evolution when no jumps take place. Red
line: probability of no-jump. Blue dashed line: population of the virtual
state conditioned to no jumps, showing that, in this case, it reaches its
maximum possible value. Blue straight line: population of the virtual state
for the general evolution. Parameters: $\Omega/\Delta_{V}=0.1,\
\Gamma/\Delta_{V}=10^{-3},\ \Delta_{2}=\Delta_{1}=0.$

This mechanism would be completely disrupted if there were additional
dissipative channels involving the virtual state. We can consider this
situation by setting $\Gamma_{V}\neq 0$ in Eq. (3), thus including a channel
of spontaneous emission from the virtual state to the ground state. The ratio
between decay rates, $\Gamma_{V}/\Gamma_{c}$, will determine whether virtual
state population occurs or not. When $\Gamma_{V}\gtrsim\Gamma_{c}$,
dissipation from $|V\rangle$ outcompetes the mechanism of population, and one
recover the simple dynamics in terms of a driven two-level system, as can be
seen in Fig. 2(a). Consistently with this, there is no longer a Liouvillian
eigenvalue corresponding to a metastable state—as shown in Fig. 2(c-d)—, as
$\Gamma_{V}$ becomes comparable to $\Gamma_{c}$, $\lambda_{2}$ is pulled
towards values $\sim\Gamma$ and the metastability gap disappears.

Generality of the HAE and its implications for entanglement generation.—The
HAE method introduced in this work can be a valuable tool to obtain analytical
insights about metastable dynamics in open quantum sytems. This can have
strong implications for quantum technological applications such as the
generation of entanglement in dissipative quantum systems. In order to
illustrate this, we now apply the HAE to describe entanglement generation in
two different systems displaying metastability.

First, we consider entanglement generation in the system of two interacting
qubits already described in Fig. 1(b), which, as noted before, maps into the
the three level model discussed so far in this text, assuming that the
occupation of the antisymmetric state
$|A\rangle\propto(|eg\rangle-|eg\rangle)$ is completely decoupled from the
dynamics and remains equal to zero. This map between models allows us to use
the density matrix elements estimated with the HAE—c.f. Fig. 2(a)—to compute
the concurrence and quantify the degreee of entanglement between the two
qubits [45, 46, 47, 48]. The results are shown in Fig. 4(a), evidencing the
formation and stabilization of entanglement at short timescales $t\sim
1/\Gamma$, due to the coherence built between the states $|ee\rangle$ and
$|gg\rangle$ via the two-photon drive. Notably, this entanglement is long-
lived, but metastable, and its survival time is given by $1/\Gamma_{c}$, i.e.,
the relaxation rate obtained in Eq. (7) via the HAE. The unconventional
population of the virtual state $|S\rangle$ thus destroys entanglement in the
long time limit.

Next, we consider another two-qubit system, sketched in the inset of Fig.
4(b). In this case, qubits experience collective decay with rate $\Gamma$,
inducing transitions $|ee\rangle\rightarrow|S\rangle$ and
$|S\rangle\rightarrow|gg\rangle$. Furthermore, each qubit is driven with the
same Rabi frequency $\Omega$, but each of them is detuned from the drive
frequency by an absolute value $\delta$ with opposite sign. This system was
introduced in Refs. [49, 50] in the context of chiral waveguides; the
Hamiltonian is $\hat{H}=\sqrt{2}\Omega(|S\rangle\langle gg|+|ee\rangle\langle
S|)+(\delta-i\Delta\gamma/2)|A\rangle\langle S|+\mathrm{H.c.}$, with
$\Gamma=2(\gamma_{R}+\gamma_{L})$ and
$\Delta\gamma\equiv\gamma_{R}-\gamma_{L}$, where $\gamma_{R}$ and $\gamma_{L}$
describe decay into right and left propagating modes respectively. This
configuration was shown to stabilize in the long time limit to the fully-
entangled dark state $|A\rangle$ provided $\Omega\gg\Delta\gamma,\delta$, as
we show explicitly in Fig. 4(b). The application of HAE to describe this
system follows exactly the same reasoning detailed above, with $|A\rangle$
playing the role of the “virtual” state that gets populated over time. Using
this technique, we are able to establish the timescale of formation of the
entangled state, which is given by $\tau\approx
24\Omega^{2}/[\Gamma(4\delta^{2}+\Delta\gamma^{2})]$. The analytical results
obtained from the HAE method match perfectly the exact calculations, as shown
in Fig. 4(b). In the limit $(\delta,\Delta\gamma)\rightarrow 0$ we find
$\tau\rightarrow\infty$, meaning that the metastable state becomes the steady
state, as reported, e.g., in Ref [51].

Figure 4:  Entanglement formation in systems of two coupled qubits. (a)
Formation of metastable entangled state in the system sketched in Fig. 1(b),
surviving for a time $\tau_{2}\sim 1/\Gamma_{c}$. Parameters:
$\Gamma/\Delta_{V}=10^{-5},\ \Delta_{2}/\Delta_{V}=0.$ (b) Stabilization of
the entangled antisymmetric state in the two-qubit system sketched in the
inset. The evolution is perfectly described by the HAE. Parameters:
$\Omega/\Gamma=1$, $\delta/\Gamma=0.01$, $\Delta\gamma/\Gamma=0.01$.

_Conclusion—._ We have shown that, in open quantum systems, off-resonant
virtual states can get populated in the long time limit even if they are not
connected to any dissipative channel, meaning that the regime where the
virtual state is not populated is metastable. We introduce a method of
hierarchical adiabatic elimination that approximates the dynamics and provides
analytical expressions of the lifetime of the metastable state. Our method can
be applied in a variety of metastable open quantum systems to obtain valuable
insights in questions such as the dissipative stabilization of entangled
states.

###### Acknowledgements.

The authors are thankful to C. Navarrete-Benlloch for insightful discussions.
C.S.M. acknowledges that the project that gave rise to these results received
the support of a fellowship from la Caixa Foundation (ID 100010434), from the
European Union’s Horizon 2020 Research and Innovation Programme under the
Marie Sklodowska-Curie Grant Agreement No. 847648, with fellowship code
LCF/BQ/PI20/11760026. A.G.T. acknowledges support from CSIC Research Platform
on Quantum Technologies PTI-001 and from Spanish project
PGC2018-094792-B-100(MCIU/AEI/FEDER, EU). All authors acknowledge financial
support from Proyecto Sinérgico CAM 2020 Y2020/TCS-6545 (NanoQuCo-CM).

## References

  * Cohen‐Tannoudji  _et al._ [1997] C. Cohen‐Tannoudji, J. Dupont‐Roc, and G. Grynberg,  _Photons and Atoms_ (Wiley, 1997).
  * Cohen‐Tannoudji  _et al._ [1998] C. Cohen‐Tannoudji, J. Dupont‐Roc, and G. Grynberg,  _Atom—Photon Interactions_ (Wiley, 1998).
  * Sakurai and Napolitano [2017] J. J. Sakurai and J. Napolitano,  _Modern Quantum Mechanics_ (Cambridge University Press, 2017).
  * Ryder [1996] L. H. Ryder,  _Quantum Field Theory_ (Cambridge University Press, 1996).
  * Gaubatz  _et al._ [1990] U. Gaubatz, P. Rudecki, S. Schiemann, and K. Bergmann,  _Population transfer between molecular vibrational levels by stimulated Raman scattering with partially overlapping laser fields. A new concept and experimental results,_ The Journal of Chemical Physics 92, 5363 (1990).
  * Bergmann  _et al._ [1998] K. Bergmann, H. Theuer, and B. W. Shore,  _Coherent population transfer among quantum states of atoms and molecules,_ Reviews of Modern Physics 70, 1003 (1998).
  * Lütkenhaus  _et al._ [1998] N. Lütkenhaus, J. I. Cirac, and P. Zoller,  _Mimicking a squeezed-bath interaction: Quantum-reservoir engineering with atoms,_ Physical Review A 57, 548 (1998).
  * Warszawski and Wiseman [2000] P. Warszawski and H. M. Wiseman,  _Adiabatic elimination in compound quantum systems with feedback,_ Physical Review A 63, 013803 (2000).
  * Brion  _et al._ [2007] E. Brion, L. H. Pedersen, and K. Mølmer,  _Adiabatic elimination in a lambda system,_ Journal of Physics A: Mathematical and Theoretical 40, 1033 (2007).
  * Dimer  _et al._ [2007] F. Dimer, B. Estienne, A. S. Parkins, and H. J. Carmichael,  _Proposed realization of the Dicke-model quantum phase transition in an optical cavity QED system,_ Physical Review A 75, 013804 (2007).
  * Burgarth  _et al._ [2019] D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa,  _Generalized Adiabatic Theorem and Strong-Coupling Limits,_ Quantum 3, 152 (2019).
  * Gamel and James [2010] O. Gamel and D. F. V. James,  _Time-averaged quantum dynamics and the validity of the effective Hamiltonian model,_ Physical Review A 82, 052106 (2010).
  * Damanet  _et al._ [2019] F. Damanet, A. J. Daley, and J. Keeling,  _Atom-only descriptions of the driven-dissipative Dicke model,_ Physical Review A 99, 033845 (2019).
  * Kaufman  _et al._ [2020] B. Kaufman, T. Rozgonyi, P. Marquetand, and T. Weinacht,  _Adiabatic elimination in strong-field light-matter coupling,_ Physical Review A 102, 063117 (2020).
  * Burgarth  _et al._ [2021] D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa,  _Eternal adiabaticity in quantum evolution,_ Physical Review A 103, 032214 (2021).
  * Garziano  _et al._ [2015] L. Garziano, R. Stassi, V. Macrì, A. F. Kockum, S. Savasta, and F. Nori,  _Multiphoton quantum Rabi oscillations in ultrastrong cavity QED,_ Physical Review A 92, 063830 (2015).
  * Garziano  _et al._ [2016] L. Garziano, V. Macrì, R. Stassi, O. Di Stefano, F. Nori, and S. Savasta,  _One Photon Can Simultaneously Excite Two or More Atoms,_ Physical Review Letters 117, 043601 (2016).
  * Stassi  _et al._ [2017] R. Stassi, V. Macrì, A. F. Kockum, O. Di Stefano, A. Miranowicz, S. Savasta, and F. Nori,  _Quantum nonlinear optics without photons,_ Physical Review A 96, 023818 (2017).
  * Comparat [2009] D. Comparat,  _General conditions for quantum adiabatic evolution,_ Physical Review A 80, 012106 (2009).
  * Mirrahimi and Rouchon [2009] M. Mirrahimi and P. Rouchon,  _Singular Perturbations and Lindblad-Kossakowski Differential Equations,_ IEEE Transactions on Automatic Control 54, 1325 (2009).
  * Paulisch  _et al._ [2014] V. Paulisch, H. Rui, H. K. Ng, and B.-G. Englert,  _Beyond adiabatic elimination: A hierarchy of approximations for multi-photon processes,_ The European Physical Journal Plus 129, 12 (2014).
  * Santos and Sarandy [2021] A. C. Santos and M. S. Sarandy,  _Generalized transitionless quantum driving for open quantum systems,_ Physical Review A 104, 062421 (2021).
  * Reiter and Sørensen [2012] F. Reiter and A. S. Sørensen,  _Effective operator formalism for open quantum systems,_ Physical Review A 85, 032111 (2012).
  * Azouit  _et al._ [2016] R. Azouit, A. Sarlette, and P. Rouchon,  _Adiabatic elimination for open quantum systems with effective Lindblad master equations,_ in  _2016 IEEE 55th Conference on Decision and Control (CDC)_ (IEEE, 2016) pp. 4559–4565.
  * Finkelstein-Shapiro  _et al._ [2020] D. Finkelstein-Shapiro, D. Viennot, I. Saideh, T. Hansen, T. Pullerits, and A. Keller,  _Adiabatic elimination and subspace evolution of open quantum systems,_ Physical Review A 101, 042102 (2020).
  * Plenio  _et al._ [1999] M. B. Plenio, S. F. Huelga, A. Beige, and P. L. Knight,  _Cavity-loss-induced generation of entangled atoms,_ Physical Review A 59, 2468 (1999).
  * Diehl  _et al._ [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller,  _Quantum states and phases in driven open quantum systems with cold atoms,_ Nature Physics 4, 878 (2008).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac,  _Quantum computation and quantum-state engineering driven by dissipation,_ Nature Physics 5, 633 (2009).
  * Chang  _et al._ [2018] D. E. Chang, J. S. Douglas, A. González-Tudela, C.-L. Hung, and H. J. Kimble,  _Colloquium : Quantum matter built from nanoscopic lattices of atoms and photons,_ Reviews of Modern Physics 90, 031002 (2018).
  * Macieszczak  _et al._ [2016] K. Macieszczak, M. GuţǍ, I. Lesanovsky, and J. P. Garrahan,  _Towards a Theory of Metastability in Open Quantum Dynamics,_ Physical Review Letters 116, 240404 (2016).
  * Macieszczak  _et al._ [2021] K. Macieszczak, D. C. Rose, I. Lesanovsky, and J. P. Garrahan,  _Theory of classical metastability in open quantum systems,_ Physical Review Research 3, 033047 (2021).
  * Varada and Agarwal [1992] G. V. Varada and G. S. Agarwal,  _Two-photon resonance induced by the dipole-dipole interaction,_ Physical Review A 45, 6721 (1992).
  * Hettich  _et al._ [2002] C. Hettich, C. Schmitt, J. Zitzmann, S. Kühn, I. Gerhardt, and V. Sandoghdar,  _Nanometer resolution and coherent optical dipole coupling of two individual molecules._ Science 298, 385 (2002).
  * Haakh and Martín-Cano [2015] H. R. Haakh and D. Martín-Cano,  _Squeezed Light from Entangled Nonidentical Emitters via Nanophotonic Environments,_ ACS Photonics 2, 1686 (2015).
  * Vivas-Viaña and Sánchez Muñoz [2021] A. Vivas-Viaña and C. Sánchez Muñoz,  _Two-photon resonance fluorescence of two interacting nonidentical quantum emitters,_ Physical Review Research 3, 33136 (2021).
  * del Valle  _et al._ [2010] E. del Valle, S. Zippilli, F. P. Laussy, A. Gonzalez-Tudela, G. Morigi, and C. Tejedor,  _Two-photon lasing by a single quantum dot in a high-Q microcavity,_ Physical Review B 81, 035302 (2010).
  * Ota  _et al._ [2011] Y. Ota, S. Iwamoto, N. Kumagai, and Y. Arakawa,  _Spontaneous Two-Photon Emission from a Single Quantum Dot,_ Physical Review Letters 107, 233602 (2011).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Kessler  _et al._ [2012] E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D. Lukin, and J. I. Cirac,  _Dissipative phase transition in a central spin system,_ Physical Review A 86, 012116 (2012).
  * Haken [2004] H. Haken,  _Synergetics_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004).
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight,  _The quantum-jump approach to dissipative dynamics in quantum optics,_ Reviews of Modern Physics 70, 101 (1998).
  * Brun [2002] T. A. Brun,  _A simple model of quantum trajectories,_ American Journal of Physics 70, 719 (2002).
  * Gerry and Knight [2004] C. Gerry and P. Knight,  _Introductory Quantum Optics_ (Cambridge University Press, 2004).
  * Haroche and Raimond [2006] S. Haroche and J.-M. Raimond,  _Exploring the Quantum: Atoms, Cavities, and Photons_ (Oxford University Press, 2006).
  * Wootters [1998] W. K. Wootters,  _Entanglement of Formation of an Arbitrary State of Two Qubits,_ Physical Review Letters 80, 2245 (1998).
  * Wootters [2001] W. K. Wootters,  _Entanglement of formation and concurrence,_ Quantum Information and Computation 1, 27 (2001).
  * Plenio and Virmani [2007] M. B. Plenio and S. Virmani,  _An introduction to entanglement measures,_ Quantum Information and Computation 7, 1 (2007).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki,  _Quantum entanglement,_ Reviews of Modern Physics 81, 865 (2009).
  * Ramos  _et al._ [2014] T. Ramos, H. Pichler, A. J. Daley, and P. Zoller,  _Quantum spin dimers from chiral dissipation in cold-atom chains,_ Physical Review Letters 113, 237203 (2014).
  * Pichler  _et al._ [2015] H. Pichler, T. Ramos, A. J. Daley, and P. Zoller,  _Quantum optics of chiral spin networks,_ Physical Review A 91, 042116 (2015).
  * Gonzalez-Tudela  _et al._ [2011] A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal,  _Entanglement of Two Qubits Mediated by One-Dimensional Plasmonic Waveguides,_ Physical Review Letters 106, 020501 (2011).

Supplementary Material

##  I Further details on the Hierarchical Adiabatic Elimination

This section elaborates all the steps involved in the hierarchical adiabatic
elimination technique presented in the main text. Let us remind that, as an
example, we consider a Hilbert space of three states,
$\mathcal{H}=\left\\{\mathinner{|{1}\rangle},\mathinner{|{2}\rangle},\mathinner{|{V}\rangle}\right\\}$.
Here, $|1\rangle$ and $|2\rangle$ represent the “real” states, and $|V\rangle$
is an off-resonant state that plays the role of a virtual state that mediates
the interactions between $|1\rangle$ and $|2\rangle$. The evolution of the
elements of the density matrix are given by the master equation of Eq. (3) in
the main text,

$\displaystyle\dot{\rho}_{V,V}$
$\displaystyle=2\Omega\imaginary\left[\rho_{1,V}-\rho_{V,2}\right],$ (S1a)
$\displaystyle\dot{\rho}_{2,2}$
$\displaystyle=-\Gamma\rho_{2,2}+2\Omega\imaginary\left[\rho_{V,2}\right],$
(S1b) $\displaystyle\dot{\rho}_{1,V}$
$\displaystyle=i\Delta_{V}\rho_{1,V}+i\Omega\left[1+\rho_{1,2}-2\rho_{V,V}-\rho_{2,2}\right],$
(S1c) $\displaystyle\dot{\rho}_{1,2}$
$\displaystyle=-(\Gamma/2)\rho_{1,2}+i\Omega\left[\rho_{1,V}-\rho_{V,2}\right],$
(S1d) $\displaystyle\dot{\rho}_{V,2}$
$\displaystyle=-\left(i\Delta_{V}+\Gamma/2\right)\rho_{V,2}-i\Omega\left[\rho_{1,2}+\rho_{2,2}-\rho_{V,V}\right].$
(S1e)

###  I.1 First adiabatic elimination

The first stage of the dynamics can be completely described within the real
Hilbert subspace,
$\mathcal{H}_{R}=\left\\{\mathinner{|{1}\rangle},\mathinner{|{2}\rangle}\right\\}$
after an adiabatic elimination of the virtual state. More specifically, this
adiabatic elimination consists in setting
$\dot{\rho}_{1,V}=\dot{\rho}_{V,2}=0$. This is done under the assumption that
the energy difference between the real and virtual subspaces is much larger
than the coupling rate, i.e. $\Delta_{V}\gg\Omega$. Thus, by substituting the
virtual coherence terms by their steady state values, the system gets
described by the following differential equations

$\displaystyle\dot{\rho}_{1,2}$
$\displaystyle\approx\left(-\frac{\Gamma}{2}+\frac{\Gamma\Omega^{2}}{(i\Gamma-2\Delta_{V})\Delta_{V}}\right)\rho_{1,2}(t)-\frac{i\Omega^{2}}{\Delta_{V}}+\left(-\frac{2\Omega^{2}}{\Gamma+2i\Delta_{V}}+\frac{i\Omega^{2}}{\Delta_{V}}\right)\rho_{2,2}(t)+\left(\frac{2\Omega^{2}}{\Gamma+2i\Delta_{V}}+\frac{2i\Omega^{2}}{\Delta_{V}}\right)\rho_{V,V}(t),$
(S2a) $\displaystyle\dot{\rho}_{2,2}$
$\displaystyle\approx\left(-\Gamma-\frac{4\Gamma\Omega^{2}}{\Gamma^{2}+4\Delta_{V}^{2}}\right)\rho_{2,2}(t)-\frac{2\Omega^{2}}{\Gamma+2i\Delta_{V}}\rho_{1,2}(t)-\frac{2\Omega^{2}}{\Gamma-2i\Delta_{V}}\rho_{2,1}(t)+\frac{4\Gamma\Omega^{2}}{\Gamma^{2}+4\Delta_{V}^{2}}\rho_{V,V}(t).$
(S2b)

In the limit $\Delta_{V}\gg\Gamma$, we can simplify them and obtain a more
familiar set of equations,

$\displaystyle\dot{\rho}_{1,2}$
$\displaystyle\approx-\frac{\Gamma}{2}\rho_{1,2}(t)-i\Omega_{2\mathrm{p}}(1-2\rho_{2,2}(t))+i\Omega_{2\mathrm{p}}\rho_{V,V}(t),$
(S3a) $\displaystyle\dot{\rho}_{2,2}$
$\displaystyle\approx-\Gamma\rho_{2,2}(t)-2\Omega_{2\mathrm{p}}\imaginary[\rho_{1,2}(t)]+\frac{\Gamma}{\Delta_{V}}\Omega_{2\mathrm{p}}\rho_{V,V}(t).$
(S3b)

These formulas show a well-known structure, since they correspond to two
resonant levels coupled via a second-order process with a Rabi frequency,
$\Omega_{2\mathrm{p}}\equiv\Omega^{2}/\Delta_{V}$ along with a standard decay,
$\Gamma$. The only addition is an extra term related to the virtual
population, $\rho_{V,V}$. However, this term evolves in a much slower
timescale than $\rho_{1,2}$ and $\rho_{2,2}$, so in these equations it can be
treated as a time-independent parameter with a fixed value. At the beginning
of the evolution, we may set $\rho_{V,V}=0$. The time-dependent analytical
solutions of these equations are well known and can be found in any quantum
optics textbook,

$\displaystyle\rho_{2,2}(t)$
$\displaystyle\approx\frac{4\Omega_{2\mathrm{p}}^{2}}{\Gamma^{2}+8\Omega_{2\mathrm{p}}^{2}}\left[1-e^{-3\Gamma
t/4}\left(\cosh(\kappa t)+\frac{3\Gamma}{4\kappa}\sinh(\kappa
t)\right)\right],$ (S4a) $\displaystyle\rho_{1,2}(t)$
$\displaystyle\approx\frac{-2i\Omega_{2\mathrm{p}}\Gamma}{\Gamma^{2}+8\Omega_{2\mathrm{p}}^{2}}\left[1-e^{-3\Gamma
t/4}\left(\cosh(\kappa
t)+\left(\frac{\kappa}{\Gamma}+\frac{3\Gamma}{16\kappa}\right)\sinh(\kappa
t)\right)\right],$ (S4b) $\displaystyle\rho_{i,V}(t)$
$\displaystyle=\rho_{V,V}(t)\approx 0\quad(i=1,2),$ (S4c)

where
$\kappa\equiv\frac{1}{2}\sqrt{\frac{\Gamma^{2}}{4}-16\Omega_{2\mathrm{p}}^{2}}$.
One can clearly see that, for these equations, the relaxation time towards a
stationary state occurs in a timescale $\sim 1/\Gamma$.

###  I.2 Second adiabatic elimination

In much longer timescales than $1/\Gamma$, $\rho_{2,2}$ and $\rho_{1,2}$ can
be considered as “fast” variables, since they relax to a steady state in a
very short time. This allows us to perform a second adiabatic elimination:
from Eqs. (S3), it is clear that, if we assume that $\rho_{V,V}$ will be
virtually unchanged in a timescale $\sim 1/\Gamma$, we could take it as a
time-independent parameter and obtain a stationary solution for $\rho_{2,2}$
and $\rho_{1,2}$ that is dependent on $\rho_{V,V}$. This quasi-steady state
will adiabatically follow any slow change of $\rho_{V,V}$. The expression of
this $\rho_{V,V}$-dependent steady state can be obtained by solving a linear
system equations of the form $M.\vec{\rho}+\vec{b}=0$ for the vector
$\vec{\rho}=\left\\{\rho_{2,2}^{SS},\rho_{2,1}^{SS},\rho_{1,2}^{SS}\right\\}$,
where $M$ and $\vec{b}$ are given by

$M=\begin{bmatrix}-\Gamma\left(1+\frac{4\Omega^{2}}{\Gamma^{2}+4\Delta_{V}^{2}}\right)&-\frac{2\Omega^{2}}{\Gamma-2i\Delta_{V}}&-\frac{2\Omega^{2}}{\Gamma+2i\Delta_{V}}\\\
-\frac{2\Omega^{2}}{\Gamma-2i\Delta_{V}}-i\Omega_{2\mathrm{p}}&-\Gamma\left(\frac{1}{2}-\frac{i\Omega_{2\mathrm{p}}}{\Gamma-2i\Delta_{V}}\right)&0\\\
-\frac{2\Omega^{2}}{\Gamma+2i\Delta_{V}}+i\Omega_{2\mathrm{p}}&0&-\Gamma\left(\frac{1}{2}-\frac{\Gamma\Omega_{2\mathrm{p}}}{i\Gamma-2\Delta_{V}}\right)\end{bmatrix}$
(S5)

and

$\vec{b}=\begin{pmatrix}\frac{4\Gamma\Omega^{2}}{\Gamma^{2}+4\Delta_{V}^{2}}\\\
\frac{2\Omega^{2}}{\Gamma-2i\Delta_{V}}-2i\Omega_{2\mathrm{p}}\\\
\frac{2\Omega^{2}}{\Gamma+2i\Delta_{V}}+2i\Omega_{2\mathrm{p}}\end{pmatrix}\rho_{V,V}(t)+\begin{pmatrix}0\\\
i\Omega_{2\mathrm{p}}\\\ -i\Omega_{2\mathrm{p}}\end{pmatrix}.$ (S6)

By solving this linear system we obtain a set of equations for the quasi-
stationary values of $\rho_{2,2}$ and $\rho_{1,2}$ that depend on the
population of the virtual state at any given time,

$\rho_{2,2}^{SS}[\rho_{V,V}(t)]=\frac{16\Omega^{4}(\Delta_{V}^{2}+\Omega^{2})}{\Gamma^{4}\Delta_{V}^{2}+32\Omega^{2}(\Delta_{V}^{2}+\Omega^{2})+4\Gamma^{2}(\Delta_{V}^{4}+3\Delta_{V}^{2}\Omega^{2}+\Omega^{4})}+\\\
\frac{4\left[\Gamma^{2}\Delta_{V}^{2}\Omega^{2}-4\Omega^{4}(\Delta_{V}^{2}+\Omega^{2})\right]}{\Gamma^{4}\Delta_{V}^{2}+32\Omega^{2}(\Delta_{V}^{2}+\Omega^{2})+4\Gamma^{2}(\Delta_{V}^{4}+3\Delta_{V}^{2}\Omega^{2}+\Omega^{4})}\rho_{V,V}(t),$
(S7)

$\rho_{1,2}^{SS}[\rho_{V,V}(t)]=\frac{2\Omega^{2}\left(-i\Gamma\Delta_{V}\left(\Gamma^{2}+4\Delta_{V}^{2}\right)-2\Gamma\Omega^{2}(\Gamma+4i\Delta_{V})-8\Omega^{4}\right)}{4\Gamma^{2}\left(3\Delta_{V}^{2}\Omega^{2}+\Delta_{V}^{4}+\Omega^{4}\right)+\Gamma^{4}\Delta_{V}^{2}+32\Omega^{4}\left(\Delta_{V}^{2}+\Omega^{2}\right)}+\\\
\frac{4\Omega^{2}\left(2\Gamma\Omega^{2}(\Gamma+4i\Delta_{V})+\Gamma\Delta_{V}(\Gamma+i\Delta_{V})(2\Delta+i\Gamma)+12\Omega^{4}\right)}{4\Gamma^{2}\left(3\Delta_{V}^{2}\Omega^{2}+\Delta_{V}^{4}+\Omega^{4}\right)+\Gamma^{4}\Delta_{V}^{2}+32\Omega^{4}\left(\Delta_{V}^{2}+\Omega^{2}\right)}\rho_{V,V}(t).$
(S8)

After the first adiabatic elimination, the differential equation that governs
the dynamics of $\rho_{V,V}$ became:

$\dot{\rho}_{V,V}(t)\approx-\frac{4\Gamma\Omega^{2}}{\Gamma^{2}+4\Delta_{V}^{2}}\rho_{V,V}+\frac{\Omega_{2\mathrm{p}}\Gamma(i\Gamma+2\Delta_{V})}{\Gamma^{2}+4\Delta_{V}^{2}}\rho_{1,2}+\frac{\Omega_{2\mathrm{p}}\Gamma(-i\Gamma+2\Delta_{V})}{\Gamma^{2}+4\Delta_{V}^{2}}\rho_{2,1}+\frac{4\Gamma}{\Gamma^{2}+4\Delta^{2}}\rho_{2,2}.$
(S9)

Substituting the pseudo-stationary values of $\rho_{1,2}$ and $\rho_{2,2}$
into this equation, we obtain a differential equation for $\rho_{V,V}(t)$
which is function of itself, i.e.,
$\dot{\rho}_{V,V}(t)=f[\rho_{V,V}(t);\rho_{1,2}^{SS}[\rho_{V,V}(t)];\rho_{2,2}^{SS}[\rho_{V,V}(t)]]$.
Namely, the differential equation for the virtual state population becomes:

$\dot{\rho}_{V,V}(t)=\frac{4\Gamma\Omega^{4}(\Gamma^{2}+4\Omega^{2})}{\Gamma^{4}\Delta_{V}^{2}+32\Omega^{2}(\Delta_{V}^{2}+\Omega^{2})+4\Gamma^{2}(\Delta_{V}^{4}+3\Delta_{V}^{2}\Omega^{2}+\Omega^{4})}-\frac{4\left[12\Gamma\Omega^{6}+\Gamma^{3}\Omega^{2}(\Delta_{V}^{2}+2\Omega^{2})\right]}{\Gamma^{4}\Delta_{V}^{2}+32\Omega^{2}(\Delta_{V}^{2}+\Omega^{2})+4\Gamma^{2}(\Delta_{V}^{4}+3\Delta_{V}^{2}\Omega^{2}+\Omega^{4})}\rho_{V,V}(t),$
(S10)

obtaining

$\rho_{V,V}(t)=\rho_{V,V}^{SS}\left(1-e^{-\Gamma_{c}t}\right),$ (S11)

where $\rho_{V,V}^{SS}$ stands for the virtual steady state population (its
expression will be given in the next section) and $\Gamma_{c}$ stands for the
relaxation rate, which corresponds to the Liouvillian gap:

$\Gamma_{c}=\frac{4\left[12\Gamma\Omega^{6}+\Gamma^{3}\Omega^{2}(\Delta_{V}^{2}+2\Omega^{2})\right]}{\Gamma^{4}\Delta_{V}^{2}+32\Omega^{2}(\Delta_{V}^{2}+\Omega^{2})+4\Gamma^{2}(\Delta_{V}^{4}+3\Delta_{V}^{2}\Omega^{2}+\Omega^{4})}.$
(S12)

We can reduce this formula under the assumption $\Delta_{V}\gg\Omega,\Gamma$
and the effective strong coupling regime $\Omega_{2\mathrm{p}}\gg\Gamma$. In
this situation, the relaxation rate reduces to

$\Gamma_{c}\approx\frac{3\Gamma\Omega^{2}}{2\Delta_{V}^{2}}.$ (S13)

###  I.3 Summary of analytic expressions for the time-dependent density matrix
elements

Once we know the analytic expression for the time-dependent virtual state
population, the remaining formulas are easily computed:

$\displaystyle\rho_{V,V}(t)$
$\displaystyle\approx\rho_{V,V}^{SS}\left[1-e^{-\Gamma_{c}t}\right],$ (S14a)
$\displaystyle\rho_{2,2}(t)$
$\displaystyle\approx\rho_{2,2}^{SS}\left[1+\frac{4\Omega^{4}\left(\Gamma^{2}+4\Delta_{V}^{2}\right)-\Gamma^{4}\Delta_{V}^{2}+16\Omega^{6}}{4\Gamma^{2}\left(3\Delta_{V}^{2}\Omega^{2}+\Delta_{V}^{4}+\Omega^{4}\right)+\Gamma^{4}\Delta_{V}^{2}+32\Omega^{4}\left(\Delta_{V}^{2}+\Omega^{2}\right)}e^{-\Gamma_{c}t}\right],$
(S14b) $\displaystyle\rho_{1,2}(t)$
$\displaystyle\approx\rho_{1,2}^{SS}\left[1-\frac{2i\Omega^{2}\left(\Gamma^{2}+4\Omega^{2}\right)\left(2\Gamma\Omega^{2}(\Gamma+4i\Delta_{V})+\Gamma\Delta_{V}(\Gamma+i\Delta_{V})(2\Delta_{V}+i\Gamma)+12\Omega^{4}\right)}{\Gamma\Delta_{V}\left(4\Gamma^{2}\left(3\Delta_{V}^{2}\Omega^{2}+\Delta_{V}^{4}+\Omega^{4}\right)+\Gamma^{4}\Delta_{V}^{2}+32\Omega^{4}\left(\Delta_{V}^{2}+\Omega^{2}\right)\right)}e^{-\Gamma_{c}t}\right],$
(S14c) $\displaystyle\rho_{1,V}(t)$
$\displaystyle\approx\rho_{1,V}^{SS}\left[1+\frac{2\Omega^{2}\left(\Gamma^{2}+4\Omega^{2}\right)\left(-2i\Gamma\Omega^{2}\left(\Gamma^{2}+6i\Gamma\Delta_{V}+2\Delta_{V}^{2}\right)+\Gamma^{2}\Delta_{V}\left(\Gamma^{2}+4\Delta_{V}^{2}\right)+8\Omega^{4}(3\Delta_{V}-2i\Gamma)\right)}{\Gamma\left(\Gamma\Delta_{V}-2i\Omega^{2}\right)\left(4\Gamma^{2}\left(3\Delta_{V}^{2}\Omega^{2}+\Delta_{V}^{4}+\Omega^{4}\right)+\Gamma^{4}\Delta_{V}^{2}+32\Omega^{4}\left(\Delta_{V}^{2}+\Omega^{2}\right)\right)}e^{-\Gamma_{c}t}\right],$
(S14d) $\displaystyle\rho_{2,V}(t)$
$\displaystyle\approx\rho_{2,V}^{SS}\left[1-\frac{\left(\Gamma^{2}+4\Omega^{2}\right)\left(\Gamma^{2}\Delta_{V}^{2}(\Gamma+2i\Delta_{V})-4\Omega^{4}(\Gamma-6i\Delta_{V})+4i\Gamma\Delta_{V}\Omega^{2}(\Gamma+i\Delta_{V})\right)}{4\Gamma^{3}\left(3\Delta_{V}^{2}\Omega^{2}+\Delta_{V}^{4}+\Omega^{4}\right)+\Gamma^{5}\Delta_{V}^{2}+32\Gamma\Omega^{4}\left(\Delta_{V}^{2}+\Omega^{2}\right)}e^{-\Gamma_{c}t}\right],$
(S14e)

where $\rho_{i,j}^{SS}$ $(i,j=1,2,V)$ are the steady state density matrix
elements,

$\displaystyle\rho_{V,V}^{SS}$
$\displaystyle=\frac{\Omega^{2}\left(\Gamma^{2}+4\Omega^{2}\right)}{\Gamma^{2}\left(\Delta_{V}^{2}+2\Omega^{2}\right)+12\Omega^{4}},\qquad\rho_{2,2}^{SS}=\frac{4\Omega^{4}}{\Gamma^{2}\left(\Delta_{V}^{2}+2\Omega^{2}\right)+12\Omega^{4}},$
(S15a) $\displaystyle\rho_{1,2}^{SS}$
$\displaystyle=-\frac{2i\Gamma\Delta_{V}\Omega^{2}}{\Gamma^{2}\left(\Delta_{V}^{2}+2\Omega^{2}\right)+12\Omega^{4}},\qquad\rho_{1,V}^{SS}=\frac{\Gamma\Omega\left(-\Gamma\Delta_{V}+2i\Omega^{2}\right)}{\Gamma^{2}\left(\Delta_{V}^{2}+2\Omega^{2}\right)+12\Omega^{4}},\qquad\rho_{2,V}^{SS}=-\frac{2i\Gamma\Omega^{3}}{\Gamma^{2}\left(\Delta_{V}^{2}+2\Omega^{2}\right)+12\Omega^{4}}.$
(S15b)

The analytic expression within the metastability regime are just the one for a
single two-level system,

$\rho_{2,2}^{M}=\frac{4\Omega_{2\mathrm{p}}^{2}}{\Gamma^{2}+8\Omega_{2\mathrm{p}}^{2}},\quad\rho_{1,2}^{M}=\frac{2i\Gamma\Omega_{2\mathrm{p}}}{\Gamma^{2}+8\Omega_{2\mathrm{p}}^{2}},\quad\rho_{i,V}^{M}\approx
0\quad\forall i=1,2,V.$ (S16)

##  II Calculation of the concurrence

We now consider that the Hilbert space
$\mathcal{H}=\left\\{\mathinner{|{1}\rangle},\mathinner{|{2}\rangle},\mathinner{|{V}\rangle}\right\\}$
is mapping the Hilbert space of two interacting quantum emitters, in such a
way that $\left\\{|1\rangle,|2\rangle,|V\rangle\right\\}$ now stand for the
triplet
$\left\\{|gg\rangle,|ee\rangle,\mathinner{|{S}\rangle}=1/\sqrt{2}(|ge\rangle+|eg\rangle\right\\}$,
where the antisymmetric state is neglected due to be completely disconnected
from the dynamics. In this scenario, we make the assumption that $\rho$ has
the following structure

$\rho\approx\begin{pmatrix}\rho_{1,1}&0&0&\rho_{1,2}\\\ 0&\rho_{V,V}&0&0\\\
0&0&0&0\\\ \rho_{2,1}&0&0&\rho_{2,2}\end{pmatrix},$ (S17)

so that the square roots of the eigenvalues of the density matrix
$\rho\sigma_{y}\otimes\sigma_{y}\rho^{*}\sigma_{y}\otimes\sigma_{y}$ are given
by

$\displaystyle\lambda_{1}$
$\displaystyle=\sqrt{\left|\rho_{1,2}\right|{}^{2}+2\left|\rho_{1,2}\right|\sqrt{-\rho_{2,2}\left(\rho_{2,2}+\rho_{V,V}-1\right)}-\rho_{2,2}\left(\rho_{2,2}+\rho_{V,V}-1\right)},$
(S18) $\displaystyle\lambda_{2}$
$\displaystyle=\sqrt{\left|\rho_{1,2}\right|{}^{2}-2\left|\rho_{1,2}\right|\sqrt{-\rho_{2,2}\left(\rho_{2,2}+\rho_{V,V}-1\right)}-\rho_{2,2}\left(\rho_{2,2}+\rho_{V,V}-1\right)}$
(S19) $\displaystyle\lambda_{3}$ $\displaystyle=\rho_{V,V}.$ (S20)

This expressions allows us to compute analytically the expression for the
concurrence. This can be done taking into account the analytic formulas shown
before, Eq. (S16). The metastable value of non-zero concurrence that is
achieved is thus given by:

$\mathcal{C}^{\text{m}}(\rho)=\frac{2\sqrt{2}|\Omega_{2\mathrm{p}}|}{\Gamma^{2}+8|\Omega_{2\mathrm{p}}|^{2}}\left(\sqrt{2|\Omega_{2\mathrm{p}}|^{2}+\Gamma^{2}+\Gamma\sqrt{\Gamma^{2}+4|\Omega_{2\mathrm{p}}|^{2}}}-\sqrt{2|\Omega_{2\mathrm{p}}|^{2}+\Gamma^{2}-\Gamma\sqrt{\Gamma^{2}+4|\Omega_{2\mathrm{p}}|^{2}}}\right).$
(S21)

This expression allows us to establish the optimum value of
$|\Omega_{2\mathrm{p}}|$ that maximizes the metastable value of the
concurrence,

$|\Omega_{2\mathrm{p}}|^{opt}=\frac{\Gamma}{2\sqrt{2}}.$ (S22)