# Density dependence of the excitation gaps in an undoped Si/SiGe double-
quantum-well heterostructure

D. Chen  Department of Physics, University of Florida, Gainesville, FL 32611,
USA  S. Cai  Department of Physics, University of Florida, Gainesville, FL
32611, USA  N.-W. Hsu  Department of Electrical Engineering and Graduate
Institute of Electronic Engineering, National Taiwan University, Taipei 10617,
Taiwan  S.-H. Huang  Department of Electrical Engineering and Graduate
Institute of Electronic Engineering, National Taiwan University, Taipei 10617,
Taiwan  Y. Chuang  Department of Electrical Engineering and Graduate Institute
of Electronic Engineering, National Taiwan University, Taipei 10617, Taiwan
E. Nielsen  Sandia National Laboratories, Albuquerque, New Mexico 87185, USA
J.-Y. Li  Department of Electrical Engineering and Graduate Institute of
Electronic Engineering, National Taiwan University, Taipei 10617, Taiwan
Taiwan Semiconductor Research Institute, Hsinchu 30078, Taiwan  C. W. Liu
Department of Electrical Engineering and Graduate Institute of Electronic
Engineering, National Taiwan University, Taipei 10617, Taiwan  T. M. Lu
Sandia National Laboratories, Albuquerque, New Mexico 87185, USA  D. Laroche
Email of Author to whom correspondence should be addressed<EMAIL_ADDRESS>Department of Physics, University of Florida, Gainesville, FL 32611, USA

###### Abstract

We report low-temperature magneto-transport measurements of an undoped Si/SiGe
asymmetric double quantum well heterostructure. The density in both layers is
tuned independently utilizing a top and a bottom gate, allowing the
investigation of quantum wells at both imbalanced and matched densities.
Integer quantum Hall states at total filling factor $\nu_{\text{T}}=1$ and
$\nu_{\text{T}}=2$ are observed in both density regimes, and the evolution of
their excitation gaps is reported as a function of density. The
$\nu_{\text{T}}=1$ gap evolution departs from the behavior generally observed
for valley splitting in the single layer regime. Furthermore, by comparing the
$\nu_{\text{T}}=2$ gap to the single particle tunneling energy,
$\Delta_{\text{SAS}}$, obtained from Schrödinger-Poisson (SP) simulations,
evidence for the onset of spontaneous inter-layer coherence (SIC) is observed
for a relative filling fraction imbalance smaller than ${\sim}50\%$.

††preprint: AIP/123-QED

Owing to their compatibility with the widespread complementary metal-oxide-
semiconductor (CMOS) technology enabling a scalable and affordable fabrication
process, Si-based heterostructures continue to attract a lot of attention for
both fundamental and applied research.Ng (2002) In particular, the substantial
low-temperature mobility improvements achieved in SiGe-based heterostructures
Lu  _et al._ (2009) have enabled various studies of low-disorder phenomena,
such as the fractional quantum hall effect,Lu  _et al._ (2012) valley
splitting Weitz  _et al._ (1996); Neyens  _et al._ (2018); McJunkin  _et al._
(2021), quantum dotsXue  _et al._ (2019); Mills  _et al._ (2019); Lawrie  _et
al._ (2020), as well as long lifetime and high fidelity electron spin
qubitsSimmons  _et al._ (2011); Russ  _et al._ (2018) that are predominantly
studied alongside GaAs-Petta  _et al._ (2005); Sarma and Pinczuk (2008) and
graphene-Bolotin  _et al._ (2009); Shimazaki  _et al._ (2015)based devices.

The vast majority of the experimental efforts on quantum transport in SiGe-
based heterostructures have focused on the development and characterization of
single layer devices, while coupled Si/SiGe bilayers Malissa  _et al._ (2006);
Laroche  _et al._ (2015) have received little attention. However, bilayer
systems present a promising platform in their own right for studying
interlayer coherence and the onset of indirect exciton condensation.Eisenstein
_et al._ (1992); Tutuc, Shayegan, and Huse (2004); Tiemann  _et al._ (2008);
Burg  _et al._ (2017); Li  _et al._ (2017) In addition, the evolution of
valley splitting in Si/SiGe bilayers is of crucial interest for the
realization of long-lived silicon spin qubits.Goswami  _et al._ (2007) Undoped
Si/SiGe double quantum wells have only been investigated in top-gated
devices,Laroche  _et al._ (2015) severely limiting the density tunability of
the devices’ bottom well. Unlike these prior studies, the device presented in
this letter hosts a bilayer system that is dual-gated, allowing independent
control of the density in both layers simultaneously. Utilizing this enhanced
control, the evolution of the density dependence of the excitation gaps at
total filling fractions $\nu_{\text{T}}=1$ and $\nu_{\text{T}}=2$ is reported
both at mismatched and matched densities.

Closely following previous heterostructure designs,Laroche  _et al._ (2015)
the double-quantum well structure presented in this letter was designed to
contain a $5\ \text{nm}$ top quantum well, a $2\ \text{nm}$
$\text{Si}_{\text{0.72}}\text{Ge}_{\text{0.28}}$ interwell barrier and a $15\
\text{nm}$ bottom quantum well. In contrast to the previous structure, the
interwell barrier’s Ge content was increased to 28% to suppress tunnelling
between the two wells. The heterostructure was grown in an ultra-high-vacuum
chemical vapor-deposition system with Si$\text{H}_{\text{4}}$ and
Ge$\text{H}_{\text{4}}$ as precursors on a p-type Si substrate. As depicted in
Fig. 1(a), the layer composition is, from bottom to top: a graded SiGe virtual
substrate with a maximal Ge composition of 14$\%$, a ${\sim}3\
\upmu\textrm{m}$ thick relaxed
$\text{Si}/\text{Si}_{\text{0.86}}\text{Ge}_{\text{0.14}}$ spacer layer, the
(strained) double quantum well structure described above, an additional
${\sim}100\ \upmu\textrm{m}$
$\text{Si}/\text{Si}_{\text{0.86}}\text{Ge}_{\text{0.14}}$ spacer layer and a
$2\ \text{nm}$ Si cap.

Figure 1: (a) Complete layer schematic of the device. (b) Top Panel: Total
density dependence vs top gate voltage with the bottom gate voltage set to 0.
Bottom Panel: Total density dependence vs bottom gate voltage with a top gate
voltage of 0.58 V. (c) Total 2D density vs mobility dependence at fixed
$V_{\text{bot}}=0$ V. A similar plot over a wider density range is shown in
the supplementary material. Top Left (Right) inset: SP simulation results
describing the band diagram and 3D electron density distribution below (above)
$n_{\text{crossover}}$, the blue dotted line in the main panel. (d) Hall
resistance (blue) and longitudinal resistivity (black) as a function of
magnetic field for $V_{\text{top}}$ = 0.63 V and $V_{\text{bot}}$ = $-22.51$
V. The top and bottom wells’ densities are $5.41\times 10^{10}\
\text{cm}^{-2}$ and $5.82\times 10^{10}\ \text{cm}^{-2}$, respectively. The
method to extract each layers’ densities is detailed in the supplementary
material.

The asymmetric quantum well structure ensures that, in the uncoupled limit,
the ground state energies in the two wells are different. Thus, during top
gate operations, a sizable electron density can accumulate in the wider bottom
layer before electrons start to populate the top layer, drastically reducing
further electron accumulation in the bottom layer through screening. Such a
design is necessary owing to the large separation, ${\sim}100\
\upmu\textrm{m}$, between the bottom gate and the bottom quantum well,
limiting the degree of density tuning achievable from the bottom gate alone.

To enable electrical measurements on the bilayer, a heterostructure field-
effect transistor (HFET) was patterned following standard fabrication
processes described in detail in a previous letter.Laroche  _et al._ (2015)
First, heavy phosphorus ion implantation was used for ohmic contacts. Prior to
the deposition of the Hall bar-shaped Ti/Au gate, 1000 cycles of atomic-layer-
deposited $\text{Al}_{2}\text{O}_{3}$ were grown. Finally, the back Ti/Au gate
was implemented by mechanically lapping the substrate down to ${\sim}100\
\upmu\text{m}$, followed by a similar oxide and metal deposition scheme.

The device was tested in a cryo-free Blue Fors dilution refrigerator with a
base temperature of $T\,{\sim}\,10\ \text{mK}$. Standard low-frequency (37.3
Hz) lock-in measurement techniques were used in a four-terminal geometry with
a constant 50 nA excitation current. Prior to the initial operation of the
device, a large negative bias ($-7$ V) was applied to evacuate negative
charges trapped at the dielectric layer interface, improving device
reliability. The measurements were found to be consistent over several thermal
cycles.

The Hall density of the device is extracted from the slope of the Hall
resistance $(\text{R}_{\text{xy}})$ at small magnetic fields ($B\leq
0.1\,\text{T}$), and its dependence on top (bottom) gate voltage at
$V_{\text{bot}}=0$ V ($V_{\text{top}}=0.58$ V) is shown in the top (bottom)
panel of Fig. 1(b). Both curves show near-linear behavior, while the bottom
gate voltage dependence has a smaller slope owing to a larger gate-to-well
separation.

Total Hall densities were tuned from $6.8\times 10^{10}\ \text{cm}^{-2}$ to
$34.2\times 10^{10}\ \text{cm}^{-2}$, with a maximal Hall mobility of
$30.8\times 10^{4}\ \text{cm}^{2}/\text{(V}\cdot\text{s)}$ achieved at the
highest density. The crossover density (the electronic density required to
start populating the top quantum well) was observed at
$n_{\text{crossover}}=8.22\times 10^{10}\ \text{cm}^{-2}$. Identification of
$n_{\text{crossover}}$ and confirmation of bilayer behavior is assessed from
the mobility drop in the density versus mobility plot, caused by inter-layer
scattering,Störmer, Gossard, and Wiegmann (1982) as exemplified in the main
panel of Fig. 1(c).

The experimental crossover density was also replicated using an iterative,
self-consistent SP simulation. Using the nominal growth parameters and
thicknesses, the experimental crossover density was reproduced via the
simulation by shifting the lower boundary condition, corresponding to a
constant offset in the bottom gate voltage. The insets of Fig. 1(c) illustrate
the simulated band diagram and electron density distribution before and after
$n_{\text{crossover}}$. To prevent charge accumulation at either interface,
the conduction band of the upper (lower) half of the heterostructure’s top
(bottom) spacer is artificially increased to be 500 meV above that of the
quantum wells.Laroche  _et al._ (2015) We also note that it was possible to
reproduce the experimental crossover density without utilizing a bottom gate
voltage offset by considering a modified structure with a top well width of
$3.8\ \text{nm}$, an interlayer barrier of $2.2\ \text{nm}$ and a bottom well
width of $15\ \text{nm}$.

Figure 2:  (a) 2D color map exhibiting the longitudinal resistivity,
$\rho_{\text{xx}}$ as a function of $n_{\text{tot}}$ and $B$ for variable
imbalanced densities. The matched density case is shown in panel (c). (b)
Density dependence of the extracted activation gaps at variable imbalanced
densities with confidence intervals for $\nu_{\text{T}}=1$ (red squares) and
$\nu_{\text{T}}=2$ (blue circles). $n_{\text{crossover}}$ is indicated by the
blue dotted line. In addition, a black line indicates the points used to fit
the estimated disorder broadening, $\Gamma$, and the linear valley splitting
coefficient, $c_{\text{B}}$. The matched density case is shown in panel (d).
Inset (b): Relative filling fraction imbalance as a function of total density
(same scale as main graph) for the variable imbalanced density case. As the
bottom well’s density marginally increases beyond the crossover density, the
top layer’s density slowly catches up to the bottom layer’s density, thereby
reducing the filling fraction imbalance. Inset (d): Temperature dependence of
the $\rho_{\text{xx}}$ minima in an Arrhenius plot for $\nu_{\text{T}}=2$ at
$V_{\text{top}}=0.64$ V and $V_{\text{bot}}=-12.368$ V in the matched density
regime.

Two density regimes are studied in this letter. The variable imbalanced
density regime describes experimental conditions where only $V_{\text{top}}$
is increased from 0 to positive values, while $V_{\text{bot}}$ is fixed at 0
V, as exemplified in Fig. 1(c). Beyond $n_{\text{crossover}}$, the density in
the top layer steadily increases towards matched densities, while the density
in the bottom layer remains nearly constant, changing by less than 4$\%$
according to simulations. In contrast, both gates are biased in the matched
density regime, yielding densities within 6$\%$ of one another in each layer,
under the assumption of ideal gate screening by the 2DEGs. Independent
contacts to both layers or a capacitive measurement scheme Dorozhkin  _et al._
(2018) would be needed to further refine this estimate. Details on the tuning
scheme can be found in the supplementary material. An instance of this later
case is shown in Fig. 1(d); where the traces of the Hall resistance,
$R_{\text{xy}}$, and of the longitudinal resistivity, $\rho_{\text{xx}}$, at a
front gate voltage of $V_{\text{top}}=0.63$ V and a back gate voltage of
$V_{\text{bot}}=-22.51$ V are presented. In both regimes, the data clearly
exhibits quantum Hall plateaus developing at $\nu_{\text{T}}=1$ and
$\nu_{\text{T}}=2$, within $3\%$ of $h/e^{2}$.

The evolution of $\rho_{\text{xx}}$ as a function of magnetic field and total
density is presented in Fig. 2 for both variable imbalanced (a) and matched
(c) densities. Quantum Hall plateaus are clearly observed for
$\nu_{\text{T}}=1$ and $\nu_{\text{T}}=2$, indicated by dotted yellow lines.
Although less developed, the $\nu_{\text{T}}=4$ minima is also observed.

The main results of this letter are presented in Fig. 2(b) and Fig. 2(d) where
the quantum Hall excitation gaps, extracted from the temperature dependence of
the $\rho_{xx}$ minima, are shown for the variable imbalanced and matched
density regimes, respectively. The inset in Fig. 2(d) shows a typical example
of the dependence $\uprho_{xx}$ has on the temperature in an Arrhenius plot.
The linear part of the plot was systematically determined and fitted to an
activation model,Boebinger  _et al._ (1990) and the extracted slope was used
to determine the excitation gap $\Delta$:
$\rho_{\text{xx}}(T)=\rho_{0}\text{e}^{-\Delta/2T}$. The deviation from
linearity at low temperatures is attributed to interstitial hopping.Ebert  _et
al._ (1983) The shaded confidence intervals in the main panel are determined
from either adding and/or subtracting an additional point to/from the
determined linear activation interval, whenever possible. For additional
details on the fitting procedure, see the supplementary material.

In the single layer regime, strained Si contains two degeneracies per Landau
Level: spins $(\uparrow\downarrow)$ and valleys (+/-). For
$n_{\text{tot}}=10\times 10^{10}\ \text{cm}^{-2}$, the single layer valley
splitting is estimated to be ${\sim}1.35\ \text{K}$ at 4.13 T.Wuetz  _et al._
(2020) At the same magnetic field, the spin gap would be ${\sim}2.79\
\text{K}$ without considering many-body effects. This implies that the
$\nu_{\text{T}}=1$ ($\nu_{\text{T}}=2$) state is associated with the valley
(spin) degree of freedom. Upon entering the bilayer regime, an additional
degeneracy, the layer degree of freedom (S/AS), comes into play. The tunneling
gap, $\Delta_{\text{SAS}}$, Boebinger  _et al._ (1990) characterizing this
degeneracy can be estimated from the SP simulation. As we do not observe
Landau level crossing in our magneto-transport measurements (Fig. 2(a)), we
can safely attribute the bilayer $\nu_{\text{T}}=1$ state to valley splitting.
The $\nu_{\text{T}}=2$ state is attributed either to spin-splitting or to
interlayer effects. In the latter case, two mechanisms can give rise to the
quantum Hall state: a single-particle state caused by interlayer tunnelling
Boebinger  _et al._ (1990); He, Sarma, and Xie (1993) or a many-body state
resulting from spontaneous interlayer coherence (SIC). He, Sarma, and Xie
(1993); Murphy  _et al._ (1994); Sawada  _et al._ (1998)

Although valley splitting in vertically coupled Si/SiGe bilayers has not been
investigated yet, the single layer regime has attracted a lot of interest;
both experimentallyLai  _et al._ (2004); Goswami  _et al._ (2007) and
theoretically.Friesen, Eriksson, and Coppersmith (2006); Friesen  _et al._
(2007) For single layers, the conduction band valley degeneracy is lifted by a
sharp quantum well interface. Atomic-scale disorder, attributed to random SiGe
alloys and steps in the quantum well, causes destructive interference that
suppresses valley splitting. Theoretical calculations predict a linear
dependence that extrapolates to a value that coincides with sample-dependent
Landau level disorder broadening at zero density.Friesen  _et al._ (2007)

The evolution of the $\nu_{\text{T}}=1$ gap in the single layer regime,
presented in red in Fig. 2(b), follows these expectations. Linearly fitting
the data below $n_{\text{crossover}}$, we obtain an estimate for the disorder
broadening $\Gamma\,{\sim}\,0.327$ K and for a linear coefficient
$c_{\text{B}}\,{\sim}\,0.29$ $\text{K}/\text{T}$, consistent with previous
studies in Si/SiGe heterostructures reporting similar results:Wuetz  _et al._
(2020) $\Gamma\,{\sim}\,0.435$ K and $c_{\text{B}}\,{\sim}\,0.326$
$\text{K}/\text{T}$. However, this behavior changes drastically in the bilayer
regime. The single layer linear dependence strongly decreases and nearly
flattens out. This remains true in the matched density case where the linear
coefficient is nearly reduced by a factor of 2 to $c_{\text{B}}\,{\sim}\,0.15$
$\text{K}/\text{T}$. Further theoretical and experimental studies of this
unexpected behavior are required to determine whether Si bilayers exhibit
large valley splitting at low magnetic fields or even spontaneous valley
polarization.Friesen, Eriksson, and Coppersmith (2006); Hossain  _et al._
(2020)

Figure 3:  (a) Simulated energy fan diagram depicting the 0 and 1 Landau
levels for variable imbalanced densities. The matched density case is shown in
panel (c). The symmetric (S) state subband is shown as black lines while the
asymmetric (AS) state subband is shown as blue lines. Spin is denoted by
arrows. The dotted line indicates where $n_{\text{crossover}}$ occurs. Valley
splitting is omitted from the calculations. The yellow shaded regions indicate
the energy gaps considered for panel (b) and (d). (b) As a function of
density, $\text{E}_{\text{2}}$ can either result from a spin state transition
(shaded black) or from a layer state transition (shaded blue). (d)
$\Delta_{\text{SAS}}$ values for exactly matched densities (blue) and for the
slightly mismatched densities measured experimentally (red); see Table S.I. in
the supplementary material. The dotted lines in panels (b) and (d) take into
consideration disorder broadening. The experimentally measured
$\nu_{\text{T}}=2$ excitation gaps are shown as green lines in panels (b) and
(d).

To understand the evolution of the $\nu_{\text{T}}=2$ state, one has to
consider the interplay between spin-splitting and interlayer effects. One
interlayer mechanism that may induce the $\nu_{\text{T}}=2$ quantum Hall state
in coherent bilayersEisenstein  _et al._ (1992) occurs when the tunneling
between the two quantum wells is sufficiently strong, resulting in an overlap
of the wave functions in each well. As a result, the individual layer
eigenstates hybridize into symmetric and antisymmetric states with single-
particle tunnel splitting gap,
$\Delta_{\text{SAS}}=\text{E}_{\text{AS}}-\text{E}_{\text{S}}$.Hu and
MacDonald (1992) $\Delta_{\text{SAS}}$ depends strongly on the height and
width of the interwell barrier.

Spontaneous interlayer coherence (SIC) is another interlayer mechanism that
can give rise to the $\nu_{\text{T}}=2$ state. As a signature of exciton
condensation between spatially separated electrons and holes, SIC induces a
coherent state that can occur solely due to Coulomb interactions in the limit
of low interlayer tunneling; when the ratio between interlayer and intralayer
Coulomb interactions is less than a critical value: $d/l_{\text{b}}\lesssim
1.8$.Spielman  _et al._ (2000); Kellogg  _et al._ (2002) Here, $d$ is the
distance between the centers of each quantum well and
$l_{\text{b}}=(\hbar/eB)^{1/2}$ is the magnetic length. The device reported in
this letter has $d/l_{\text{b}}$ values ranging between $0.62$ and $1.24$,
well within the range where SIC can be observed. Furthermore, the ratio of the
SP simulated $\Delta_{\text{SAS}}$ gap to the Coulomb energy ($e^{2}/\epsilon
l_{\text{b}}$) is ${\sim}5\times 10^{-4}$, indicating that the device is
indeed in the Coulomb dominated regime.

Both spin-splitting and interlayer tunnelling are captured within the SP
simulations. In addition to the two different cases used to reproduce
$n_{\text{crossover}}$, two different interwell barrier heights,
$\phi_{\text{b}}$, were used. The first one corresponds to a Ge content of
28.8%, giving $\phi_{\text{b}}=141.2$ meV; while the other assumes a 25.8$\%$
Ge content, giving $\phi_{\text{b}}=122.8$ meV.Schäffler (1997) Together, we
have thus simulated four cases. We note that further increasing the Ge content
of the barrier would have negligible effect on its height. Only the case with
as-grown parameters and higher interwell barrier is shown in Fig. 3. The other
three cases are shown in the supplementary material. The general shape of the
density dependence of $\Delta_{\text{SAS}}$ is nearly identical in all four
cases. However, the magnitude of $\Delta_{\text{SAS}}$ can be up to
${\sim}96\%$ larger for the case of thinner top well width and weaker
interlayer barrier, as opposed to as-grown parameters and higher interwell
barrier. Panel (a) (Panel (c)) of Fig. 3 depicts the evolution of the energy
fan diagrams for the variable imbalanced (matched) density regime. Here, only
spin splitting and $\Delta_{\text{SAS}}$ are considered while valley splitting
is omitted from the calculations, resulting in 4 energy sub-levels per Landau
level. As such, the simulation strictly provides information about
$\text{E}_{\text{2n}}$, the energy of the even quantum Hall states. The yellow
shaded regions denote the simulated $\text{E}_{2}$. The simulated second
Landau energy level, $\text{E}_{2}$, for the variable imbalanced density case
is also shown in Fig. 3(b) while the simulated $\Delta_{\text{SAS}}$ for the
matched density case is shown in Fig. 3(d).

From theory,Wuetz  _et al._ (2020) we expect the experimental gap to be
$\Delta_{2}=\text{E}_{2}-\text{E}_{1}-\Gamma=\text{E}_{2}-\Delta_{1}-2\Gamma$.
In the variable imbalanced density case at low density, the simulated
$\text{E}_{2}$ and the experimentally measured $\Delta_{2}$ are similar.
Starting in the single layer regime, a linear dependence is observed as long
as $\text{E}_{2}$ originates from a spin state transition. A downturn in the
gap magnitude is then observed with increasing density once the state
originates from $\Delta_{\text{SAS}}$. However, the measured magnitude of
$\Delta_{2}$ is notably smaller than the simulated $\text{E}_{2}$. In
addition, the theoretical downturn in $\text{E}_{2}$ is reduced to the point
of flattening, and its experimental onset occurs at a reduced density,
$9.8\times 10^{10}\ \text{cm}^{-2}$. As shown in Fig. 3(b), disorder
broadening alone cannot explain the observed discrepancy. A non-linearity of
the valley splitting excitation gap, $\Delta_{1}$, might explain this
discrepancy, but the same non-linearity prevents us from extrapolating the
experimentally measured $\Delta_{1}$ to lower densities and magnetic fields.

At higher densities, the behavior of $\Delta_{2}$ departs more strongly from
the simulated $\text{E}_{2}$, shown in Fig. 3(b). The downturn is rapidly
inverted and the overall magnitude of $\Delta_{2}$ overtakes that of the
disorder-broadening corrected $\Delta_{\text{SAS}}$ as the layers approach
matched densities, where $\Delta_{\text{SAS}}$ is calculated to be at its
weakest. Both the qualitative departure from the $\Delta_{\text{SAS}}$ density
dependence as well as the increased $\Delta_{2}$ value are consistent with the
onset of SIC.Murphy  _et al._ (1994) This is further confirmed by the behavior
of $\Delta_{2}$ in the density matched case (Fig. 2(d)). Here, $\Delta_{2}$
has a concavity opposite to that of $\Delta_{\text{SAS}}$ (Fig. 3(d)) and a
magnitude of up to 10 times larger than that of the disorder-broadened
$\Delta_{\text{SAS}}$, at $\Delta_{2}=2.37$ K. Overall, the SP simulations
show that neither $\Delta_{\text{SAS}}$, spin-splitting nor Landau level
spacing alone reproduces the qualitative behavior obtained experimentally.
Extrapolating the density dependence of $\Delta_{1}$ to lower magnetic fields
is also insufficient to account for the differences between the measured
$\Delta_{2}$ and the simulated $\Delta_{\text{SAS}}$.

Considering the density imbalance between the layers, shown in the inset of
Fig. 2(b), the onset of $\Delta_{2}$’s magnitude increase coincides with a
relative density imbalance
$\frac{\Delta\nu}{\nu_{\text{T}}}=\frac{\nu_{1}-\nu_{2}}{2}\approx 0.5$ at
$n_{\text{tot}}=11.07\times 10^{10}\ \text{cm}^{-2}$. Prior studies in GaAs
bilayersSpielman  _et al._ (2004); Champagne  _et al._ (2008) have reported
that SIC remains intact for relative density imbalances as large as
$\frac{\Delta\nu}{\nu_{\text{T}}}=0.5$, which is in quantitative agreement
with our observations.

Three main parameters influence the strength of $\Delta_{2}$ induced by SIC:
the ratio between interlayer and intralayer Coulomb interactions
$d/l_{\text{b}}$, the tunneling strength $\Delta_{\text{SAS}}$ and the density
imbalance $\Delta\nu$. As the total density is increased in our device,
$d/l_{\text{b}}$ increases while $\Delta_{\text{SAS}}$ decreases. $\Delta\nu$
also decreases for the density imbalanced regime. While the excitation gap has
previously been shown to increase with decreasing $d/l_{\text{b}}$He, Sarma,
and Xie (1993); Murphy  _et al._ (1994) and increasing $\Delta\nu$,Tutuc  _et
al._ (2003); Clarke  _et al._ (2005); Wiersma  _et al._ (2006); Champagne  _et
al._ (2008) tunnelling has been shown to weaken the SIC induced quantum
state.He, Sarma, and Xie (1993); Lay  _et al._ (1994) Therefore, we attribute
the increase in $\Delta_{2}$ with increasing density to the reduced
tunnelling. Indeed, going towards higher density decreases tunneling as the
wave functions are pushed farther apart. The downturn of $\Delta_{2}$ at the
highest densities in the matched density case coincides with the saturation of
the simulated $\Delta_{\text{SAS}}$ (Fig. 3(d)) and can be explained by the
increased $d/l_{\text{b}}$ at constant tunnelling. These effects are likely
noticeable in our device since tunnelling ($\Delta_{\text{SAS}}\geq 500$ mK)
is over 1000 times larger than prior studies probing similar effects.Tutuc
_et al._ (2003); Clarke  _et al._ (2005); Zhang  _et al._ (2013) While our
observations are consistent with a SIC-induced $\Delta_{2}$, we cannot fully
rule out that a significant reduction of valley splitting at lower magnetic
fields is at the origin of this phenomenon since the value of $\Delta_{1}$ at
$B\sim 2$ T is neither know experimentally nor theoretically.

In conclusion, we have fabricated and performed magneto-transport measurements
on a novel back-gated Si/SiGe bilayer system. The density dependence of the
excitation gaps for $\nu_{\text{T}}=1$ and $\nu_{\text{T}}=2$ has been studied
in both the variable imbalanced and matched density regimes and has been
compared to values obtained from SP simulations. The $\nu_{\text{T}}=1$ state
is attributed to valley-splitting while the $\nu_{\text{T}}=2$ state is
attributed to spin-splitting at low density and interlayer effects at larger
density. Stark departure in the density evolution of the $\nu_{\text{T}}=1$
excitation gap between the single and the bilayer regimes has been reported,
highlighting the need for theoretical work on valley splitting in bilayer
Si/SiGe structures. Evidence for exciton condensation has been observed in the
$\nu_{\text{T}}=2$ state for $\frac{\Delta\nu}{\nu_{\text{T}}}<0.5$. Achieving
independent contacts to both layers in order to perform tunneling
conductance,Spielman  _et al._ (2000) Coulomb dragKellogg  _et al._ (2002) and
counterflow measurementsTutuc, Shayegan, and Huse (2004); Yoon  _et al._
(2010) would be required to further confirm this conclusion. Reducing the
separation between the bottom quantum well and the bottom gate as well as
utilizing larger magnetic fields would enable studies of these effects over a
wider density range. These results highlight the prospects of Si/SiGe bilayers
as a platform for valley splitting tuning and exciton condensation, which may
have potential applications towards realizing SiGe-basedXu  _et al._ (2020) or
Bose-Einstein condensate based Byrnes, Wen, and Yamamoto (2012) qubits.

See the supplementary material for details on the methods used to achieve
matched density between the layers, how to extract the excitation gaps as well
as results from the four cases used in the SP simulations.

This work was performed, in part, at the Center for Integrated
Nanotechnologies, an Office of Science User Facility operated for the U.S.
Department of Energy (DOE) Office of Science. Sandia National Laboratories is
a multimission laboratory managed and operated by National Technology &
Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell
International, Inc., for the U.S. DOE’s National Nuclear Security
Administration under contract DE-NA-0003525. The views expressed in the
article do not necessarily represent the views of the U.S. DOE or the United
States Government. This work was partially supported by the National High
Magnetic Field Laboratory. The National High Magnetic Field Laboratory is
supported by the National Science Foundation through NSF/DMR-1644779 and the
State of Florida. The National Taiwan University (NTU) group is supported by
the Ministry of Science and Technology (MOST), Taiwan under contracts
MOST-110-2634-F-009-027, MOST-109-2622-8-002-003, and MOST-110-2218-E-002-030.
The authors have no conflicts to disclose.

## Data Availability Statement

The data that support the findings of this study are openly available in
IR@UF, reference number IR00011549.

## References

## References

  * Ng (2002) K. K. Ng,  _Complete guide to semiconductor devices_ (IEEE Press, 2002).
  * Lu  _et al._ (2009) T. Lu, D. Tsui, C.-H. Lee, and C. Liu, Applied Physics Letters 94, 182102 (2009).
  * Lu  _et al._ (2012) T. Lu, W. Pan, D. Tsui, C.-H. Lee, and C. Liu, Physical Review B 85, 121307 (2012).
  * Weitz  _et al._ (1996) P. Weitz, R. Haug, K. Von Klitzing, and F. Schäffler, Surface science 361, 542 (1996).
  * Neyens  _et al._ (2018) S. F. Neyens, R. H. Foote, B. Thorgrimsson, T. Knapp, T. McJunkin, L. Vandersypen, P. Amin, N. K. Thomas, J. S. Clarke, D. Savage, _et al._ , Applied Physics Letters 112, 243107 (2018).
  * McJunkin  _et al._ (2021) T. McJunkin, E. R. MacQuarrie, L. Tom, S. F. Neyens, J. P. Dodson, B. Thorgrimsson, J. Corrigan, H. E. Ercan, D. E. Savage, M. G. Lagally, R. Joynt, S. N. Coppersmith, M. Friesen, and M. A. Eriksson, Phys. Rev. B 104, 085406 (2021).
  * Xue  _et al._ (2019) X. Xue, T. Watson, J. Helsen, D. R. Ward, D. E. Savage, M. G. Lagally, S. N. Coppersmith, M. Eriksson, S. Wehner, and L. Vandersypen, Physical Review X 9, 021011 (2019).
  * Mills  _et al._ (2019) A. Mills, D. Zajac, M. Gullans, F. Schupp, T. Hazard, and J. Petta, Nature communications 10, 1 (2019).
  * Lawrie  _et al._ (2020) W. Lawrie, H. Eenink, N. Hendrickx, J. Boter, L. Petit, S. Amitonov, M. Lodari, B. Paquelet Wuetz, C. Volk, S. Philips, _et al._ , Applied Physics Letters 116, 080501 (2020).
  * Simmons  _et al._ (2011) C. Simmons, J. Prance, B. Van Bael, T. S. Koh, Z. Shi, D. Savage, M. Lagally, R. Joynt, M. Friesen, S. Coppersmith, _et al._ , Physical review letters 106, 156804 (2011).
  * Russ  _et al._ (2018) M. Russ, D. M. Zajac, A. J. Sigillito, F. Borjans, J. M. Taylor, J. R. Petta, and G. Burkard, Physical Review B 97, 085421 (2018).
  * Petta  _et al._ (2005) J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, Science 309, 2180 (2005).
  * Sarma and Pinczuk (2008) S. D. Sarma and A. Pinczuk,  _Perspectives in quantum hall effects: Novel quantum liquids in low-dimensional semiconductor structures_ (John Wiley & Sons, 2008).
  * Bolotin  _et al._ (2009) K. I. Bolotin, F. Ghahari, M. D. Shulman, H. L. Stormer, and P. Kim, Nature 462, 196 (2009).
  * Shimazaki  _et al._ (2015) Y. Shimazaki, M. Yamamoto, I. V. Borzenets, K. Watanabe, T. Taniguchi, and S. Tarucha, Nature Physics 11, 1032 (2015).
  * Malissa  _et al._ (2006) H. Malissa, D. Gruber, D. Pachinger, F. Schäffler, W. Jantsch, and Z. Wilamowski, Superlattices and Microstructures 39, 414 (2006).
  * Laroche  _et al._ (2015) D. Laroche, S.-H. Huang, E. Nielsen, C. W. Liu, J.-Y. Li, and T.-M. Lu, Applied Physics Letters 106, 143503 (2015).
  * Eisenstein  _et al._ (1992) J. Eisenstein, G. Boebinger, L. Pfeiffer, K. West, and S. He, Physical review letters 68, 1383 (1992).
  * Tutuc, Shayegan, and Huse (2004) E. Tutuc, M. Shayegan, and D. Huse, Physical review letters 93, 036802 (2004).
  * Tiemann  _et al._ (2008) L. Tiemann, W. Dietsche, M. Hauser, and K. Von Klitzing, New Journal of Physics 10, 045018 (2008).
  * Burg  _et al._ (2017) G. W. Burg, N. Prasad, B. Fallahazad, A. Valsaraj, K. Kim, T. Taniguchi, K. Watanabe, Q. Wang, M. J. Kim, L. F. Register, _et al._ , Nano letters 17, 3919 (2017).
  * Li  _et al._ (2017) J. Li, T. Taniguchi, K. Watanabe, J. Hone, and C. Dean, Nature Physics 13, 751 (2017).
  * Goswami  _et al._ (2007) S. Goswami, K. Slinker, M. Friesen, L. McGuire, J. Truitt, C. Tahan, L. Klein, J. Chu, P. Mooney, D. W. Van Der Weide, _et al._ , Nature Physics 3, 41 (2007).
  * Störmer, Gossard, and Wiegmann (1982) H. Störmer, A. Gossard, and W. Wiegmann, Solid State Communications 41, 707 (1982).
  * Dorozhkin  _et al._ (2018) S. Dorozhkin, A. Kapustin, I. Fedorov, V. Umansky, K. Von Klitzing, and J. Smet, Journal of Applied Physics 123, 084301 (2018).
  * Boebinger  _et al._ (1990) G. Boebinger, H. Jiang, L. Pfeiffer, and K. West, Physical review letters 64, 1793 (1990).
  * Ebert  _et al._ (1983) G. Ebert, K. Von Klitzing, C. Probst, E. Schuberth, K. Ploog, and G. Weimann, Solid State Communications 45, 625 (1983).
  * Wuetz  _et al._ (2020) B. P. Wuetz, M. P. Losert, A. Tosato, M. Lodari, P. L. Bavdaz, L. Stehouwer, P. Amin, J. S. Clarke, S. N. Coppersmith, A. Sammak, _et al._ , Physical review letters 125, 186801 (2020).
  * He, Sarma, and Xie (1993) S. He, S. D. Sarma, and X. Xie, Physical Review B 47, 4394 (1993).
  * Murphy  _et al._ (1994) S. Murphy, J. Eisenstein, G. Boebinger, L. Pfeiffer, and K. West, Physical review letters 72, 728 (1994).
  * Sawada  _et al._ (1998) A. Sawada, Z. Ezawa, H. Ohno, Y. Horikoshi, Y. Ohno, S. Kishimoto, F. Matsukura, M. Yasumoto, and A. Urayama, Physical review letters 80, 4534 (1998).
  * Lai  _et al._ (2004) K. Lai, W. Pan, D. Tsui, S. Lyon, M. Mühlberger, and F. Schäffler, Physical review letters 93, 156805 (2004).
  * Friesen, Eriksson, and Coppersmith (2006) M. Friesen, M. Eriksson, and S. Coppersmith, Applied physics letters 89, 202106 (2006).
  * Friesen  _et al._ (2007) M. Friesen, S. Chutia, C. Tahan, and S. Coppersmith, Physical Review B 75, 115318 (2007).
  * Hossain  _et al._ (2020) M. S. Hossain, M. Ma, K. Rosales, Y. Chung, L. Pfeiffer, K. West, K. Baldwin, and M. Shayegan, arXiv preprint arXiv:2011.06721 (2020).
  * Hu and MacDonald (1992) J. Hu and A. MacDonald, Physical Review B 46, 12554 (1992).
  * Spielman  _et al._ (2000) I. Spielman, J. Eisenstein, L. Pfeiffer, and K. West, Physical review letters 84, 5808 (2000).
  * Kellogg  _et al._ (2002) M. Kellogg, I. Spielman, J. Eisenstein, L. Pfeiffer, and K. West, Physical review letters 88, 126804 (2002).
  * Schäffler (1997) F. Schäffler, Semiconductor Science and Technology 12, 1515 (1997).
  * Spielman  _et al._ (2004) I. Spielman, M. Kellogg, J. Eisenstein, L. Pfeiffer, and K. West, Physical Review B 70, 081303 (2004).
  * Champagne  _et al._ (2008) A. Champagne, A. Finck, J. Eisenstein, L. Pfeiffer, and K. West, Physical Review B 78, 205310 (2008).
  * Tutuc  _et al._ (2003) E. Tutuc, S. Melinte, E. De Poortere, R. Pillarisetty, and M. Shayegan, Physical review letters 91, 076802 (2003).
  * Clarke  _et al._ (2005) W. Clarke, A. Micolich, A. Hamilton, M. Simmons, C. B. Hanna, J. Rodriguez, M. Pepper, and D. Ritchie, Physical Review B 71, 081304 (2005).
  * Wiersma  _et al._ (2006) R. Wiersma, J. Lok, W. Dietsche, K. v. Klitzing, D. Schuh, and W. Wegscheider, physica status solidi (b) 243, 3643 (2006).
  * Lay  _et al._ (1994) T. Lay, Y. Suen, H. Manoharan, X. Ying, M. Santos, and M. Shayegan, Physical Review B 50, 17725 (1994).
  * Zhang  _et al._ (2013) D. Zhang, S. Schmult, V. Venkatachalam, W. Dietsche, A. Yacoby, K. von Klitzing, and J. Smet, Physical Review B 87, 205304 (2013).
  * Yoon  _et al._ (2010) Y. Yoon, L. Tiemann, S. Schmult, W. Dietsche, K. von Klitzing, and W. Wegscheider, Physical review letters 104, 116802 (2010).
  * Xu  _et al._ (2020) Y. Xu, F. Unseld, A. Corna, A. Zwerver, A. Sammak, D. Brousse, N. Samkharadze, S. Amitonov, M. Veldhorst, G. Scappucci, _et al._ , Applied Physics Letters 117, 144002 (2020).
  * Byrnes, Wen, and Yamamoto (2012) T. Byrnes, K. Wen, and Y. Yamamoto, Physical Review A 85, 040306 (2012).


# Dephasing-assisted entanglement in a system of strongly coupled qubits

I. V. Vovcenko  Moscow Institute of Physics and Technology, 141700, Moscow
region, Russia  Institute for Theoretical and Applied Electromagnetics,
125412, Moscow, Russia  V. Yu. Shishkov  Moscow Institute of Physics and
Technology, 141700, Moscow region, Russia  Institute for Theoretical and
Applied Electromagnetics, 125412, Moscow, Russia  Dukhov Research Institute of
Automatics, 127055, Moscow, Russia  E. S. Andrianov<EMAIL_ADDRESS>Moscow Institute of Physics and Technology, 141700, Moscow region, Russia
Institute for Theoretical and Applied Electromagnetics, 125412, Moscow, Russia
Dukhov Research Institute of Automatics, 127055, Moscow, Russia

(August 27, 2024)

entanglement, strong-coupling, dephasing

## abstract

Creation of entangled states of quantum systems with low decoherence rates is
a cornerstone in practical implementation of quantum computations. Processes
of separate dephasing in each qubit in experimentally feasible systems is
commonly accepted to destroy entanglement. In this work, we consider a system
of two strongly coupled qubits that interact with dephasing reservoirs. We
demonstrate that interaction with dephasing reservoirs can contribute to the
formation of a long-lived mixed entangled state with nonzero concurrence. The
weight of the subradiant state in this mixed state tends toward unity if the
dephasing rate is much larger than the radiative rate and less than the
coupling constant between qubits. The lifetime of this state is proportional
to the exponent of the ratio of the coupling constant to environmental
temperature and can be, by orders of magnitude, larger than the system’s
characteristic dephasing and dissipation times. Therefore, high dephasing,
along with strong coupling, contributes to the creation of an entangled state
with a long lifetime. This result paves the way for creation of long-lived
entangled states.

##  I Introduction

Entangled states of two or more qubits are the building blocks of quantum
computers and elements of quantum communication systems Krantz  _et al._
(2019). The concept of entanglement was first proposed for the pure state,
which can be described in terms of wave function Bennett  _et al._ (1996);
Kumar  _et al._ (2002); Wootters (1998); Hill and Wootters (1997). However,
all real systems inevitably interact with the environment. This leads to a
state of the system ceases to be pure and becomes a mixed state, which can be
described only in terms of a density matrix. The concept of entanglement can
be extended to mixed states, as well Aolita  _et al._ (2015). The system-
environment interaction leads to processes of dissipation and dephasing or
decoherence. The first leads to change in both energy and coherence, while the
second process does not change the energy and results in destruction of
coherence. This manifests as dissipation of only the non-diagonal density
matrix elements. In real systems, the dephasing rate is up to five orders of
magnitude higher than the dissipation rate Blanton  _et al._ (1997); Shim and
Guyot-Sionnest (1999) and it is responsible for the fast destruction of the
entanglement.

One basic element for entanglement creation is an ensemble of two two-level
systems (TLSs) interacting with the environment and, possibly, each other
Barenco (1995); Poyatos  _et al._ (1997); Benhelm  _et al._ (2008); Zheng and
Guo (2000); Nielsen and Chuang (2010). The problem of entanglement creation in
such a system has been investigated in many works Braun (2002); Ficek and
Tanaś (2002, 2006); Scala  _et al._ (2008); Benatti  _et al._ (2010);
Campagnano  _et al._ (2010); Scala  _et al._ (2011), and several opportunities
to create and conserve entanglement exist. One can enhance the interaction
between TLSs in such a way that the ground state of the system becomes
entangled Scala  _et al._ (2008, 2011) or use coherent external drive to move
the system into the entangled state Li and Paraoanu (2009); Cecoi  _et al._
(2018). The use of common reservoirs for two qubits is another way to create
entanglement Braun (2002); Ficek and Tanaś (2006); Benatti  _et al._ (2010);
Hor-Meyll  _et al._ (2009).

Free space modes of electromagnetic field are an example of common reservoirs
for dipole moments of atoms or molecules. When they occupy subwavelength
volume, an entangled subradiant state, which is antisymmetric with respect to
mutual permutation, is protected from interaction with a common radiative
reservoir. For this reason, the creation of a subradiant state is particularly
interesting. There are quantum algorithms and quantum memory schemes based on
manipulation with such a state Facchinetti  _et al._ (2016); Kalachev and
Samartsev (2005); Plankensteiner  _et al._ (2015); Kalachev (2007); Scully
(2015). Additionally, one can use interfaces, e.g., waveguides Gonzalez-
Ballestero  _et al._ (2014), plasmonic waveguides Gonzalez-Tudela  _et al._
(2011); Martin-Cano  _et al._ (2011) or chiral waveguides Gonzalez-Ballestero
_et al._ (2015), to move the system, more effectively, towards the subradiant
state.

Thus, a common reservoir can contribute to the formation of entanglement. This
concerns the dephasing reservoir as well. If one considers dephasing
reservoirs common for two subsystems, long-lived entanglement is possible
Campagnano  _et al._ (2010). At the same time, dephasing reservoirs usually
are separate reservoirs because they are associated with nuclear vibration
degrees of freedom for each molecule or a quantum dot Breuer and Petruccione
(2002); Carmichael (2009). Thus, such type of reservoirs are expected to
result in the entanglement destruction.

However, when the system contains interacting subsystems, relaxation of one
system inevitably leads to relaxation of the others. The last statement can be
elucidated in the Born-Markov approximation by means of a master equation in
the Lindblad form. As the Lindblad approach gives an answer consistent with
thermodynamic laws, using eigenstates of the whole system is necessary
Shishkov  _et al._ (2019); Kosloff (2013). In the case of well-separated
eigenstates, the Lindblad superoperator describes the transitions between
eigenstates. Because the eigenstate of the interacting subsystem is a
superposition of isolated subsystem eigenstates, mentioned Lindblad
superoperators will result in a cross-relaxation process when relaxation of
one subsystem leads to relaxation of the others Shishkov  _et al._ (2019,
2020). Recently, construction of the correct master equation for the system of
coupled TLSs, interacting with different types of reservoirs, was actively
investigated Santos and Semiao (2014); De Chiara  _et al._ (2018); Cattaneo
_et al._ (2019). In view of the above, revising the problem of two-qubit
entanglement, in case of strong coupling, using recent progress in physics of
open quantum systems is important Cattaneo  _et al._ (2019).

In this paper, we consider the system of two strongly coupled qubits each of
which interacts with its own dephasing reservoir. We show that dephasing not
only results in the relaxation of non-diagonal terms of density matrix, but
also leads to the transition between super- and sub-radiant states. If
coupling constant between qubits is much larger than dephasing rate then there
is exist long-lived mixed entangled state. Its lifetime is proportional to the
exponent of ratio of coupling constant to environmental temperature. If
coupling constant between qubits is much larger than the temperature of
dephasing reservoir, the lifetime of the mentioned state is much larger than
the characteristic dephasing and dissipation times of the system. It is
demonstrated that the weight of entangled subradiant state in this mixed state
tends to unity, if dephasing rate is much higher than the radiative rate,
i.e., high dephasing contributes to entanglement. As a result, concurrence of
this state is greater than zero which is a criterion of a state to be
entangled. The obtained results open the way to use dephasing as a source for
the creation of long-lived entanglement.

##  II Formation of long-lived entanglement by separate dephasing resevoirs in
the system of strongly-coupled qubits

We consider the system of two strongly coupled TLSs with transition
frequencies ${\omega_{1}}\approx{\omega_{2}}$ lying in the optical range and
with the distance between TLSs $r$ that is much smaller than the optical
wavelength, i.e. $r\ll\lambda_{1,2}=2\pi c/\omega_{1,2}$. We suppose that both
TLSs interact with a common radiative reservoir, and each TLS interacts with
its own dephasing reservoir. Each dephasing reservoir is associated with
vibrational degrees of freedom (e.g. phonons Carmichael (2009)) of qubits.
Because we are interested in the case of nearby qubit frequencies,
$\left|{{\omega_{1}}-{\omega_{2}}}\right|\ll{\omega_{1,2}}$ we will use the
rotating wave approximation Scully and Zubairy (1999). Thus, the Hamiltonian
of the system takes the form Breuer and Petruccione (2002); Carmichael (2009):

${\hat{H}_{\rm{S}}}=\hbar{\omega_{1}}\hat{\sigma}_{1}^{\dagger}{\hat{\sigma}_{1}}+\hbar{\omega_{2}}\hat{\sigma}_{2}^{\dagger}{\hat{\sigma}_{2}}+\hbar\Omega\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{2}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{1}}}\right).$
(1)

Here, ${\hat{\sigma}_{1,2}}$ are lowering operators for the first and second
TLSs, respectively. The last term in Eq. (1) describes the dipole-dipole
interaction in the rotating wave approximation when dipoles occupy
subwavelength volume Agarwal (1974),
$\Omega=\left({{{\bf{d}}_{1}}{{\bf{d}}_{2}}-3\left({{{\bf{d}}_{1}}{\bf{n}}}\right)\left({{{\bf{d}}_{2}}{\bf{n}}}\right)}\right)/\hbar{r^{3}}$
is the constant of interaction, ${{\bf{d}}_{1,2}}$ are dipole transition
matrix elements of the first and second TLSs, and ${\bf{n}}$ is the unit
vector from one TLS to another. First, we consider the case of zero detuning,
${\omega_{1}}={\omega_{2}}=\omega$, and ${\Omega}>0$. As such, the system
eigenstates are the excited state of two TLSs, $\left|{ee}\right\rangle$, the
superradiant symmetric state
$\left|s\right\rangle=\left({\left|{eg}\right\rangle+\left|{ge}\right\rangle}\right)/\sqrt{2}$,
the subradiant antisymmetric state
$\left|as\right\rangle=\left({\left|{eg}\right\rangle-\left|{ge}\right\rangle}\right)/\sqrt{2}$,
and the ground state $\left|{gg}\right\rangle$, with the eigenfrequencies
$2\omega$, $\omega+\Omega$, $\omega-\Omega$, and $0$, respectively.

The master equation for the system density matrix can be obtained by the Born-
Markov approximation via standard procedure Breuer and Petruccione (2002);
Kosloff (2013); Davies (1974), assuming that the reservoirs are in thermal
equilibrium. As a result, we obtain the following Lindblad equation (see
Appendix):

$\begin{array}[]{l}{{\dot{\hat{\rho}}}_{\rm{S}}}(t)=-i{\hbar^{-1}}[{{\hat{H}}_{\rm{S}}},{{\hat{\rho}}_{\rm{S}}}(t)]+\sum\limits_{{\rm{k}}=1,2}{\sum\limits_{{\rm{j}}=1,2,3}{\frac{{{\gamma_{{\rm{k}}{\rm{,dp}}}}\left({2\Omega{\theta_{\rm{j}}}}\right)}}{2}\mathcal{L}_{\rm{dp,kj}}\left[\hat{\rho}\right]}}+\\\
+\sum\limits_{{\rm{k}}=1,2}{\frac{{{\gamma_{{\rm{rad}}}}\left({-\left({\omega+{{(-1)}^{{\rm{k}}-1}}\Omega}\right)}\right)}}{2}\mathcal{L}_{\rm{rad,k}}\left[\hat{\rho}\right]}+\sum\limits_{{\rm{k}}=1,2}{\frac{{{\gamma_{{\rm{rad}}}}\left({\left({\omega+{{(-1)}^{{\rm{k}}-1}}\Omega}\right)}\right)}}{2}\mathcal{L}_{\rm{rad,k}}\left[\hat{\rho}\right]}{\kern
1.0pt},\end{array}$ (2)

where Lindblad superoperators have the standard form, ${{\cal
L}_{i}}\left[{\hat{\rho}}\right]=2{{\hat{L}}_{i}}{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{i}^{\dagger}-{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{i}^{\dagger}{{\hat{L}}_{i}}-\hat{L}_{i}^{\dagger}{{\hat{L}}_{i}}{{\hat{\rho}}_{\rm{S}}}(t)$,
and ${\cal
L}_{i}^{\dagger}\left[{\hat{\rho}}\right]=2\hat{L}_{i}^{\dagger}{{\hat{\rho}}_{\rm{S}}}(t){{\hat{L}}_{i}}-{{\hat{\rho}}_{\rm{S}}}(t){{\hat{L}}_{i}}\hat{L}_{i}^{\dagger}-{{\hat{L}}_{i}}\hat{L}_{i}^{\dagger}{{\hat{\rho}}_{\rm{S}}}(t)$.
Operators $\hat{L}_{i}$ are found to be

$\begin{array}[]{l}{{\hat{L}}_{{\rm{rad}}{\rm{,1}}}}={{\hat{\sigma}}_{1}}+{{\hat{\sigma}}_{2}}-\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}-{{\hat{\sigma}}_{1}}\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}},\
{{\hat{L}}_{{\rm{rad}}{\rm{,2}}}}=\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}+{{\hat{\sigma}}_{1}}\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}},\\\
{{\hat{L}}_{{\rm{dp}}{\rm{,11}}}}={{\hat{L}}_{{\rm{dp}}{\rm{,21}}}}=\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}/2+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}/2,\
{{\hat{L}}_{{\rm{dp}}{\rm{,12}}}}=-{{\hat{L}}_{{\rm{dp}}{\rm{,22}}}}=\left({\hat{\sigma}_{1}^{\dagger}+\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}-{{\hat{\sigma}}_{2}}}\right)/4,\\\
{{\hat{L}}_{{\rm{dp}}{\rm{,13}}}}=-\,{{\hat{L}}_{{\rm{dp}}{\rm{,23}}}}=\hat{L}_{{\rm{dp}}{\rm{,12}}}^{\dagger},\end{array}$
(3)

where parameter ${\theta_{\rm{j}}}=0,\,1,\,-1$ for ${\rm{j}}=1,\,2,\,3$. The
dissipation constants, $\gamma_{i}$, are determined according to
${\gamma_{i}}\left(\omega\right)=\int\limits_{-\infty}^{\infty}{d\tau\exp\left({-i\omega\tau}\right)\left\langle{{\hat{\tilde{R}}_{i}}(t+\tau){\hat{\tilde{R}}_{i}}(t)}\right\rangle}$
Shishkov  _et al._ (2019); Kosloff (2013), which, for a reservoir in thermal
equilibrium at temperature $T$, satisfies the Kubo-Martin-Shwinger (KMS)
condition
${\gamma_{i}}\left(\omega\right)=\exp\left({-\hbar\omega/k{T_{i}}}\right){\gamma_{i}}\left({-\omega}\right)$.

The main consequence of usage of correct Lindblad operators (3) is that
dephasing not only results in the relaxation of non-diagonal density matrix
terms (in the eigenstate basis) but also leads to the transition between
super- and subradiant states. Indeed, operators
${\hat{L}_{\rm{dp},11}},\;{\hat{L}_{\rm{dp},21}}$ describe dephasing of both
TLSs while operators ${\hat{L}_{\rm{dp},12}},\;{\hat{L}_{\rm{dp},22}}$ and
${\hat{L}_{\rm{dp},13}},\;{\hat{L}_{\rm{dp},23}}$ describe forward and
backward energy flows between the system and dephasing reservoirs,
respectively. Notably, these operators have nonzero matrix elements between
the system eigenstates $\left|s\right\rangle$ and $\left|as\right\rangle$.
This means that dephasing reservoirs result in transitions between super- and
subradiant states. The reason for the energy flow between the system and the
dephasing reservoir is interaction between the TLSs and expansion of the
system reservoir interaction operators over eigenstates of interacting TLSs.
These terms do not appear in the local approach, in which Lindblad
superoperators are assumed as in the case of non-interacting TLSs Shishkov
_et al._ (2019, 2020). Note that in Cecoi  _et al._ (2018) it has been
mentioned that the dephasing reservoirs in the case of coupled qubits result
in transition between super- and subradiant states. However, in Cecoi  _et
al._ (2018) it has been considered stationary entanglement which is created by
the external coherent drive. In contrast, here we are interested in the
temporal dynamics of the entanglement between two qubits without any external
coherent source.

Operator ${\hat{L}_{\rm{rad},1}}$ defines the influence of the radiative
reservoir and consist of two parts. The first part represents relaxation from
the state $\left|{ee}\right\rangle$ to the symmetric state
$\left|s\right\rangle$ and from the state $\left|s\right\rangle$ to the ground
state $\left|{gg}\right\rangle$. The second part is energy flow from one TLS
to another. Operator ${\hat{L}_{\rm{rad},2}}$ describes relaxation of
nondiagonal density matrix terms.

Note that the matrix elements between the states $\left|as\right\rangle$ and
$\left|{gg}\right\rangle$ of the Lindblad superoperators, arising from both
dephasing and radiative reservoirs, are equal to zero. Thus, neither radiative
nor dephasing reservoirs lead to transition from the subradiant state to the
ground state. The only possibility of the system to leave the state
$\left|as\right\rangle$ is a temperature induced transition to the state
$\left|s\right\rangle$ which has a higher eigenfrequency. According to the KMS
condition, the rate of this transition is smaller than the rate of transition
from the state $\left|s\right\rangle$ to the state $\left|as\right\rangle$ by
the factor $\exp\left({-2\hbar{\Omega}/kT}\right)$. From these qualitative
arguments, it follows that if the system reaches the state
$\left|as\right\rangle$, then the time during which the system remains in this
state is larger than the characteristic relaxation times by the factor
$\exp\left({2\hbar{\Omega}/kT}\right)$. Because this state is entangled, one
can expect that entanglement between two TLSs will also be conserved during
the time that is larger than the characteristic relaxation times by the factor
$\exp\left({2\hbar{\Omega}/kT}\right)$.

Figure 1:  The dependencies of density matrix elements $p_{ee}$ (black dashed
line), $p_{gg}$ (blue dash-dotted line), $p_{as}$ (orange solid line), $p_{s}$
(green dashed line) (a), and entropy (blue dashed line), concurrence (red
dash-dotted line), antisymmetric state occupancy (orange solid line) (b) on
time when the system interacts with two separate dephasing reservoirs and
common radiative reservoir.

Numerical simulations and analytical evaluations confirm these qualitative
arguments. Results of the simulation of master equation (2), with the initial
condition $\hat{\rho}_{S}(0)=|ee\rangle\langle ee|$, are shown in the Fig. 1.
Experimentally feasible parameters for molecules and quantum dots are
${\gamma_{\rm{dp}}}\simeq{10^{11}}{\ \rm{s}^{-1}}$,
${\gamma_{\rm{rad}}}\simeq{10^{8}}-{10^{9}}{\ \rm{s}^{-1}}$. The value of the
Rabi constant depends on the dipole moment matrix element. We are interested
in the case of strong coupling. Thus, we consider the case of quantum dots
with a large dipole moment $\simeq 50{\ \rm{D}}$ Blanton  _et al._ (1997);
Shim and Guyot-Sionnest (1999) and $\simeq 10{\ \rm{nm}}$ distance between
them, such that the corresponding Rabi constant is ${\Omega}\simeq{2\cdot
10^{12}}{\ \rm{s}^{-1}}$. In numerical simulations, we use the value
$\simeq{\rm{0.01\ eV}}$ as a dimensionless unit of frequency such that the
other parameters are
${\gamma_{\rm{dp}1,2}}\left({-2{\Omega}}\right)=\gamma_{\rm{dp}}={2\cdot
10^{-2}}$,
${\gamma_{\rm{rad}}}\left({-\omega\pm\Omega}\right)=\gamma_{\rm{rad}}={2\cdot
10^{-4}}$, ${T_{\rm{dp}}}={T_{\rm{rad}}}=2\cdot{10^{-2}}$ ($\simeq 5{\rm K}$),
and $\Omega={10^{-1}}$.

In Fig. 1a, the dependencies of the non-zero matrix elements of the density
matrix are shown, namely, the probabilities ${p_{ee}}$, ${p_{as}}$, ${p_{s}}$,
and ${p_{gg}}$ of the system to be in the states $\left|{ee}\right\rangle$,
$\left|as\right\rangle$, $\left|s\right\rangle$, and
$\left|{gg}\right\rangle$, respectively. There are two stages in the system
dynamics. During the first stage ($t\leq{10^{5}}$), the system starting from
the state $\left|{ee}\right\rangle$ relaxes and primarily occupies the
antisymmetric subradiant state. At the second stage, ($t\geq{10^{5}}$), the
system leaves this state and relaxes to the ground state of the system,
$\left|{gg}\right\rangle$. We calculated the entropy,
$S=-{\rm{Tr}}\left({\hat{\rho}_{\rm{S}}\log\hat{\rho}_{\rm{S}}}\right)$, and
concurrence,
$C=\max\\{0,\sqrt{{\lambda_{1}}}-\sqrt{{\lambda_{2}}}-\sqrt{{\lambda_{3}}}-\sqrt{{\lambda_{4}}}\\}$,
where ${\lambda_{i}}$ are eigenvalues of the matrix
$\hat{\rho}_{\rm{S}}\left[{\left({{{\hat{\sigma}}_{y}}\otimes{{\hat{\sigma}}_{y}}}\right)\hat{\rho}_{\rm{S}}^{*}\left({{{\hat{\sigma}}_{y}}\otimes{{\hat{\sigma}}_{y}}}\right)}\right]$
Wootters (1998). The evolution of entropy, concurrence and occupancy of the
antisymmetric state are presented in Fig. 1b.

At the first stage, entropy of the system grows (Fig. 1b, blue dashed line)
because the system transits from the state $\left|{ee}\right\rangle$ to the
subspace formed by states $\left|{eg}\right\rangle$ and
$\left|{ge}\right\rangle$. After the first stage finishes, the system
primarily occupies the antisymmetric state (Fig. 1b, orange solid line) and
stays in this state for a time which is much larger than all relaxation times
in the system. As a consequence, the entropy at the second stage decreases
while concurrence grows from zero to 0.95 (Fig. 1b, red dash-dotted line).
This indicates that the system is in the entangled state. Finally, the system
relaxes to the ground state which is also accompanied by the temporary grows
of entropy.

As has been mentioned in the introduction, in general, separate reservoirs
destruct entanglement. In our case, the main contribution to the formation of
entanglement is made by separated dephasing reservoirs. We illustrate this
fact with two examples II.1 and II.2 described below. Additionally, we examine
the system with help of local approach in the third example II.3 to show that
the local approach does not reveal the entanglement.

###  II.1 Long-lived entanglement in the case of one dephasing reservoir

In the first example, we turn off the dephasing reservoirs associated with the
first TLS, while the dephasing reservoir associated with the second TLS is not
changed (i.e., $\gamma_{\rm{1,dp}}=0$ in Eq. (2)). In other words, we break
the symmetry of interaction between the TLSs and dephasing reservoirs.

Numerical simulation shows (see Fig. 2) that no significant changes of
entropy, concurrence, or occupancy of the antisymmetric state are observed,
and time dynamics are similar to those from Fig. 1b. This is a consequence of
the cross-relaxation processes described above.

Figure 2: The dependencies of density matrix elements $p_{ee}$ (black dashed
line), $p_{gg}$ (blue dash-dotted line), $p_{as}$ (orange solid line), $p_{s}$
(green dashed line) (a), and entropy (blue dashed line), concurrence (red
dash-dotted line), antisymmetric state occupancy (orange solid line) (b) on
time with one dephasing reservoir removed.

###  II.2 Disappearance of the entanglement in the case of absence of
dephasing reservoirs

In the second example, we remove both dephasing reservoirs
($\gamma_{\rm{1,dp}}=\gamma_{\rm{2,dp}}=0$), such that only the common
radiative reservoir interacts with the TLSs. This case is analogous to the
Dicke model of superradiance for two TLSs Gross and Haroche (1982), such that
only states $\left|{ee}\right\rangle$, $\left|s\right\rangle$, and
$\left|{gg}\right\rangle$ are occupied (see Fig. 3a). Relaxation of the system
to the ground state with the rate $\simeq{2\gamma_{\rm{rad}}}$ occurs.
Concurrence, in all moments of time, is zero (Fig. 3b, red dash-dotted line),
which means that the system is in the disentangled state. Indeed, in the Dicke
model of superradiance, entanglement is not generated Yu (2016); Wolfe and
Yelin (2014).

Figure 3: The dependencies of density matrix elements $p_{ee}$ (black dashed
line), $p_{gg}$ (blue dash-dotted line), $p_{as}$ (orange solid line), $p_{s}$
(green dashed line) (a), and entropy (blue dashed line), concurrence (red
dash-dotted line), symmetric state occupancy (green dashed line) (b) on time
with both dephasing reservoirs removed.

###  II.3 Disappearance of entanglement in the local approach

In third example, we use local approach to describe the system. Under local
approach, the Lindblad superoperators, arising from dephasing reservoirs, are
assumed to coincide with the ones for noninteracting TLSs and have the form

$\begin{array}[]{l}{{\hat{L}}_{{\rm{rad}}{\rm{,1}}}}={{\hat{\sigma}}_{1}}+{{\hat{\sigma}}_{2}}-\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}-{{\hat{\sigma}}_{1}}\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}},\
{{\hat{L}}_{{\rm{rad}}{\rm{,2}}}}=\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}+{{\hat{\sigma}}_{1}}\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}},\\\
{{\hat{L}}_{{\rm{dp}}{\rm{1}}}}={\hat{\sigma}_{z}^{\left(1\right)}},\
{{\hat{L}}_{{\rm{dp}}{\rm{2}}}}={\hat{\sigma}_{z}^{\left(2\right)}}.\end{array}$
(4)

In the Fig. 4, the results of this model are presented. Both symmetric and
antisymmetric states have equal populations at all moments of time.
Concurrence is zero in all moments of time, thus, entanglement does not occur.

Figure 4: The dependencies of density matrix elements $p_{ee}$ (black dashed
line), $p_{gg}$ (blue dash-dotted line), $p_{as}$ (orange solid line), $p_{s}$
(green dashed line) (a), and entropy (blue dashed line), concurrence (red
dash-dotted line), antisymmetric state occupancy (orange solid line),
symmetric state occupancy (green dashed line) (b) on time in local approach.

##  III Analytical evaluation of the probability of subradiant state
occupation and its lifetime

To clarify the nature of the long-lived entangled state, from master equation
(2), we find the approximate solution for $p_{as}$ in the limit
$\hbar{\Omega}\gg kT$ for time scale $t\geq
2\gamma_{\rm{dp}}^{-1}{\rm{exp}}\left({{\hbar\Omega}/{kT_{\rm{dp}}}}\right)$.
We rewrite the master equation in terms of occupancies ${p_{ee}}$, ${p_{s}}$,
${p_{as}}$, ${p_{gg}}$ of the system eigenstates $\left|{ee}\right\rangle$,
$\left|s\right\rangle$, $\left|{as}\right\rangle$, $\left|{gg}\right\rangle$
and obtain the following equations:

${\dot{p}_{ee}}=-2{\gamma_{{\rm{rad}}}}{p_{ee}},\ \ \ \ \ \ \
{\dot{p}_{s}}=2{\gamma_{{\rm{rad}}}}{p_{ee}}-\left({{\gamma_{{\rm{dp}}}}/2+2{\gamma_{{\rm{rad}}}}}\right){p_{s}}+\frac{{{\gamma_{{\rm{dp}}}}}}{2}{\rm{exp}}\left({-2\Omega/{T_{{\rm{dp}}}}}\right){p_{as}},$
(5)
${\dot{p}_{as}}=\frac{{{\gamma_{{\rm{dp}}}}}}{2}{p_{s}}-\frac{{{\gamma_{{\rm{dp}}}}}}{2}{\rm{exp}}\left({-2\Omega/{T_{{\rm{dp}}}}}\right){p_{as}},\
\ \ \ \ \ \ {\dot{p}_{gg}}=2{\gamma_{{\rm{rad}}}}{p_{s}},$ (6)

with the initial condition

${p_{ee}}\left(0\right)=1,\ \ \ \ \ \ \
{p_{s}}\left(0\right)={p_{as}}\left(0\right)={p_{gg}}\left(0\right)=0.$ (7)

Because the case of strong coupling is examined, $\hbar\Omega\gg kT$, at time
scale $\gamma_{{\rm{rad}}}^{-1}\ll
t\ll\gamma_{{\rm{dp}}}^{-1}{\rm{exp}}\left({\hbar\Omega/k{T_{{\rm{dp}}}}}\right)$,
we can suppose that the rate of transition from the sub- to superradiant state
is zero. In such an assumption, we obtain the following quasi-stationary
solution:

$p_{ee}^{{\rm{qs}}}=p_{s}^{{\rm{qs}}}=0,\ \ \ \ \ \ \
p_{as}^{{\rm{qs}}}=\frac{{{\gamma_{{\rm{dp}}}}/2}}{{2{\gamma_{{\rm{rad}}}}+{\gamma_{{\rm{dp}}}}/2}},\
\ \ \ \ \ \
p_{gg}^{{\rm{qs}}}=\frac{{2{\gamma_{{\rm{rad}}}}}}{{2{\gamma_{{\rm{rad}}}}+{\gamma_{{\rm{dp}}}}/2}}.$
(8)

For the time scale $t\geq
2\gamma_{{\rm{dp}}}^{-1}{\rm{exp}}\left({\hbar\Omega/k{T_{{\rm{dp}}}}}\right)$
we can use quasi-stationary solution (8) as the initial condition for Eqs. (5)
– (6). The answer takes the form:

$p_{as}^{an}\left(t\right)\simeq\frac{{{\gamma_{\rm{dp}}/2}}}{{{2\gamma_{{\rm{rad}}}}+{\gamma_{\rm{dp}}/2}}}{\rm{exp}}\left({-\frac{{{\gamma_{\rm{dp}}}{\gamma_{{\rm{rad}}}}t}}{{{\gamma_{\rm{dp}}/2}+{2\gamma_{{\rm{rad}}}}}}{\rm{exp}}\left({-2\hbar\Omega/k{T_{{\rm{dp}}}}}\right)}\right),\
\ \ \ \ \ \ p_{gg}^{an}\left(t\right)\simeq 1-p_{as}^{an}\left(t\right).$ (9)

Analytical evaluation of antisymmetric state occupancy,
${p^{an}_{as}}\left(t\right)$, is in a good agreement with exact solution,
${p_{as}}\left(t\right)$ at times $t\geq
2\gamma_{\rm{dp}}^{-1}{\rm{exp}}\left({{\hbar\Omega}/{kT_{\rm{dp}}}}\right)$
(see Fig 5). Eq. (9) shows that the characteristic lifetime of the subradiant
state, in the limit ${\gamma_{\rm{dp}}}\gg{\gamma_{{\rm{rad}}}}$, is of the
order ${t_{{\rm{ent}}}}\simeq
0.5\gamma_{{\rm{rad}}}^{-1}{\rm{exp}}\left({2{\hbar\Omega}/{kT_{\rm{dp}}}}\right)\gg
2\gamma_{\rm{dp}}^{-1},\,0.5\gamma_{{\rm{rad}}}^{-1}$, i.e., much greater than
the characteristic time of dissipation and dephasing. Note that the long-lived
state described above is a special case of metastable states proposed in
Macieszczak  _et al._ (2016).

Figure 5: The dependence of the probability ${p_{as}}$ of the system to be in
the antisymmetric state obtained from the numerical simulation of master
equation (2) (red line) and $p_{as}^{an}$ from analytical approximation (9)
(blue dashed line).

##  IV Concurrence in the case of non-zero detuning between qubit frequencies

The long-lived entangled state is also exists in the case of non-zero detuning
between qubit frequencies, $\omega_{1}\neq\omega_{2}$. In the Appendix, we
generalize the master equation for this case (see Eqs. (26) – (29)). Modeling
this equation with initial condition
$\hat{\rho}\left(0\right)=\left|{ee}\right\rangle\left\langle{ee}\right|$, one
can get concurrence as a function of time and detuning (we fix ${\omega_{2}}$
and change ${\omega_{1}}$) (see Fig 6). Its seen that the lifetime of the
entangled state decreases when detuning grows. However, there exists detuning
range ($\left|{{\omega_{1}}-{\omega_{2}}}\right|\sim 0.1\Omega$), for which
the system stays in the entangled state (concurrence is grater than $0.9$) for
the time which is much larger than the radiative and dephasing times. If
detuning is out of this range, entangled state also is formed, but its
lifetime is much shorter and its concurrence is smaller.

##  V Conclusion

To conclude, in this work, we consider the system of two coupled two-level
qubits with optical frequencies, which interact with the electromagnetic field
of free space modes and with dephasing reservoirs. We show that when the
coupling constant between qubits is much larger than the dephasing rate and
dephasing rate is much larger than radiative rate, a long-lived mixed
entangled state exists. The weight of the subradiant state in this mixed
state, as well as concurrence, are equal to one up to the ratio
${\gamma_{\rm{rad}}}/{\gamma_{\rm{dp}}}$. The lifetime of this state is larger
than the characteristic radiative time by the factor
${\rm{exp}}\left({2{\hbar\Omega}/{kT_{\rm{dp}}}}\right)$.

The obtained result, in principle, can be extended to the case of many qubits.
If one creates such interaction between qubits that symmetry with respect to
permutation is protected at least approximately, then super- and subradiant
states may still be eigenstates of the system. For example, when qubits are
placed in the ring and approximate interaction Hamiltonian is
$\hbar\Omega\sum\limits_{k=1}^{N}{\left({\hat{\sigma}_{k}^{\dagger}{{\hat{\sigma}}_{k-1}}+\hat{\sigma}_{k-1}^{\dagger}{{\hat{\sigma}}_{k}}}\right)}$
(where, by definition, ${\hat{\sigma}_{0}}={\hat{\sigma}_{n}}$ and $N$ is
even), there are the superradiant and subradiant states which for the case of
one net excitation have the form
$\frac{1}{{\sqrt{N}}}\sum\limits_{k}{\left|{g...{e_{k}}...g}\right\rangle}$
and
$\frac{1}{{\sqrt{N}}}\sum\limits_{k}{{{(-1)}^{k}}\left|{g\ldots{e_{k}}\ldots
g}\right\rangle}$, respectively (see Scully (2015)). There is dephasing-
induced transition from the super- to subradiant state, analogously to the
case of two qubits. Careful investigation we leave on further works.

Figure 6: The dependence of concurrence on time $t$ and frequency detuning
${\omega_{1}}-{\omega_{2}}$. The frequency of the first qubit is changed while
other parameters of the system are the same as in the main text.

The obtained results pave the way for using dephasing as a resource for the
creation of long-lived entanglement at experimentally realizable conditions.
In the case of molecules, the characteristic dipole moment is $\simeq 1{\
\rm{D}}$. If one achieves the characteristic distance between them $\simeq
10{\ \rm{nm}}$, then the Rabi constant of interaction is
${\Omega}\simeq{10^{9}}{{\ \rm{s}}^{-1}}$. Thus, the factor
${\rm{exp}}\left({2\hbar{\Omega}/k{T_{\rm{dp}}}}\right)$ is much larger than
unity for temperatures $T_{\rm{dp}}\leq 5{\ \rm{mK}}$, and the lifetime of the
entangled state can be $\geq\gamma_{\rm{rad}}^{-1}\simeq{10^{-9}}\ \rm{s}$. If
one considers semiconductor quantum dots with the characteristic dipole moment
$\simeq 50{\ \rm{D}}$ Blanton  _et al._ (1997); Shim and Guyot-Sionnest (1999)
and distance between them $\simeq 10{\ \rm{nm}}$, one obtains ${\Omega}\simeq
2\cdot{10^{12}}{{\ \rm{s}}^{-1}}$ and the factor
$\hbar{\Omega}/k{T_{\rm{dp}}}$ much larger than unity for the temperatures
$T_{\rm{dp}}\leq 10{\ \rm{K}}$. For example, for temperature
$T_{\rm{dp}}\leq{5\ \rm{K}}$ the factor
$\exp\left({2\hbar{\Omega}/k{T_{\rm{dp}}}}\right)\simeq{10^{3}}$ and the
lifetime of the entangled state may reach $\simeq 1\ \mu{\rm{s}}$.

The predicted value of concurrence may be useful in contest of DLSZ protocol
Kimble (2008); Laurat  _et al._ (2007), where the value $0.9\pm 0.3$ of
concurrence is used. In this work, we obtain the value $C\approx 0.9$. Also,
this value is quite similar to maximal concurrence of two entangled photons
achieved in recent experiments on optical quantum networks Jones  _et al._
(2018).

Note that the simple model developed in this work may shed light on the
recently demonstrated robustness of the subradiant state decay time for an
ensemble of cold atoms with the increasing of temperature Weiss  _et al._
(2019). In this regard, the toy-model that describes long-lived entanglement
may be useful for qualitative analyses.

## Appendix

Here we derive the master equation (2) and Lindblad superoperators (3) from
the main text and consider the case of non-zero detuning.

To describe the interaction of two coupled qubits with the electromagnetic
field of free space and separate dephasing reservoirs, we start with the
Hamiltonian of the entire system. This Hamiltonian consists of three terms.
${\hat{H}_{\rm{R}}}$ represents the electromagnetic (EM) field of free space
modes and degrees of freedom of dephasing reservoirs. The term
${\hat{H}_{\rm{S}}}$ describes the coupled qubits, and ${\hat{H}_{\rm{SR}}}$
is responsible for the interaction between the EM field and qubit dipole
moments, and for the interaction between dephasing reservoirs and qubits. The
corresponding Hamiltonians have the form:

$\hat{H}={\hat{H}_{\rm{S}}}+{\hat{H}_{\rm{R}}}+{\hat{H}_{\rm{SR}}},\
{\hat{H}_{\rm{S}}}=\hbar{\omega_{1}}\hat{\sigma}_{1}^{\dagger}{\hat{\sigma}_{1}}+\hbar{\omega_{2}}\hat{\sigma}_{2}^{\dagger}{\hat{\sigma}_{2}}+\hbar\Omega\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{2}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{1}}}\right),$
(10)
${\hat{H}_{\rm{R}}}=\sum\limits_{\rm{k}}{\hbar{\omega_{\rm{rad},{k}}}\hat{b}_{\rm{rad},{k}}^{\dagger}{{\hat{b}}_{\rm{rad},{k}}}}+\sum\limits_{\rm{k}}{\hbar{\omega_{\rm{1dp},{k}}}\hat{b}_{\rm{1dp},{k}}^{\dagger}{{\hat{b}}_{\rm{1dp},{k}}}}+\sum\limits_{\rm{k}}{\hbar{\omega_{\rm{2dp},{k}}}\hat{b}_{\rm{2dp},{k}}^{\dagger}{{\hat{b}}_{\rm{2dp},{k}}}},$
(11)
${\hat{H}_{\rm{SR}}}={\hat{H}_{\rm{S},\rm{1dp}}}+{\hat{H}_{\rm{S},\rm{2dp}}}+{\hat{H}_{\rm{S},\rm{rad}}},\
{\hat{H}_{\rm{S},\rm{rad}}}=\sum\limits_{\rm{k}}{\hbar\gamma_{\rm{k}}^{\rm{rad}}}\left({\hat{\sigma}_{1}^{\dagger}+{{\hat{\sigma}}_{1}}+\hat{\sigma}_{2}^{\dagger}+{{\hat{\sigma}}_{2}}}\right)\left({\hat{b}_{\rm{rad},{k}}^{\dagger}+{{\hat{b}}_{\rm{rad},{k}}}}\right),$
(12)
${\hat{H}_{\rm{S},\rm{1dp}}}=\sum\limits_{\rm{k}}{\hbar\gamma_{\rm{k}}^{\rm{1dp}}\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}}\left({\hat{b}_{\rm{1dp},{k}}^{\dagger}+{{\hat{b}}_{\rm{1dp},{k}}}}\right),\
{\hat{H}_{\rm{S},\rm{2dp}}}=\sum\limits_{\rm{k}}{\hbar\gamma_{\rm{k}}^{\rm{2dp}}\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}}\left({\hat{b}_{\rm{2dp},{k}}^{\dagger}+{{\hat{b}}_{\rm{2dp},{k}}}}\right).$
(13)

Here ${\hat{\sigma}_{1}}$ and ${\hat{\sigma}_{2}}$ are lowering operators for
the first and second qubit, respectively; ${\omega_{1}}$ and ${\omega_{2}}$
are transition frequencies. The last term in Eq. (10) describes dipole-dipole
interaction in the rotating wave approximation Scully and Zubairy (1999). The
coupling constant between qubits is
$\Omega=\left({{{\bf{d}}_{1}}{{\bf{d}}_{2}}-3\left({{{\bf{d}}_{1}}{\bf{n}}}\right)\left({{{\bf{d}}_{2}}{\bf{n}}}\right)}\right)/\hbar{r^{3}}$,
${{\bf{d}}_{1,2}}$ are matrix elements of the first and second TLS dipole
transitions, and ${\bf{n}}$ is the unit vector from one TLS to another.
Operators ${\hat{b}_{\rm{1dp},k}},\;{\hat{b}_{\rm{2dp},k}}$ denote lowering
operators for the k-th degree of freedom of the first and second dephasing
reservoirs, respectively. We suppose that they obey commutation relations
$\left[{{{\hat{b}}_{\rm{1dp},k}},\;\hat{b}_{\rm{1dp},k^{\prime}}^{\dagger}}\right]={\delta_{\rm{kk^{\prime}}}}$,
$\left[{{{\hat{b}}_{\rm{2dp},k}},\;\hat{b}_{\rm{2dp},k^{\prime}}^{\dagger}}\right]={\delta_{\rm{kk^{\prime}}}}$,
and
$\left[{{{\hat{b}}_{\rm{1dp},k}},\;\hat{b}_{\rm{2dp},k^{\prime}}^{\dagger}}\right]=0$;
this means that the reservoirs are independent. Operator
${\hat{b}_{\rm{rad},k}}$ denotes the lowering operator for the k-th EM field
mode of free space with commutation relation
$\left[{{{\hat{b}}_{\rm{rad},k}},\;\hat{b}_{\rm{rad},k^{\prime}}^{\dagger}}\right]={\delta_{\rm{kk^{\prime}}}}$
The constants $\gamma_{\rm{k}}^{\rm{1dp}}$, $\gamma_{\rm{k}}^{\rm{2dp}}$, and
$\gamma_{\rm{k}}^{\rm{rad}}$ are interaction constants of the qubits with the
k-th degree of freedom of dephasing and radiative reservoirs, respectively.

To compute time dependence of the density matrix, first, we rewrite the von-
Neuman equation in interaction representation using substitutions
$\hat{\rho}(t)=U\hat{\tilde{\rho}}(t)U^{\dagger}$,
${\hat{H}_{\rm{SR}}}(t)=U{\hat{\tilde{H}}_{\rm{SR}}}(t)U^{\dagger}$, where
$U=\exp\left(-i(\hat{H}_{\rm{S}}+\hat{H}_{\rm{R}})t/\hbar\right)$.

Note that all terms describing the system-reservoir interaction have the form
$\hat{S}\hat{R}$, where $\hat{S}$ depends on system variables and $\hat{R}$
depends on the reservoir variables. Thus, each term in ${\hat{H}_{\rm{SR}}}$
can be presented in the form
${\hat{\tilde{H}}_{\rm{SR}}}(t)=\sum\limits_{m}{{{\hat{\tilde{S}}}_{m}}{{\hat{\tilde{R}}}_{m}}}$,
(subindex m denotes each of reservoirs) where
${\hat{\tilde{S}}_{m}}(t)=\exp\left({i\hat{H}_{\rm{S}}t/\hbar}\right){\hat{S}_{m}}(t)\exp\left({-i\hat{H}_{\rm{S}}t/\hbar}\right)$,
${\hat{\tilde{R}}_{m}}(t)=\exp\left({i\hat{H}_{\rm{R}}t/\hbar}\right){\hat{R}_{m}}(t)\exp\left({-i\hat{H}_{\rm{R}}t/\hbar}\right)$.
After this substitution we get the von-Neuman equation in the following form:

$\frac{{\partial\hat{\tilde{\rho}}(t)}}{{\partial
t}}=\frac{i}{\hbar}\left[{\hat{\tilde{\rho}},{{\hat{\tilde{H}}}_{\rm{SR}}}(t)}\right].$
(14)

First, we consider $\omega_{1}=\omega_{2}$. To calculate
${\hat{\tilde{S}}_{m}}\left(t\right)$, one can establish the action of the
operator
$\exp\left({i{{\hat{H}}_{\rm{S}}}t}\right)\hat{\sigma}(t)\exp\left({-i{{\hat{H}}_{\rm{S}}}t}\right)$
on an arbitrary state $\left|\psi\right\rangle$ and, then, write this action
as a linear combination of operators from the basis set, which contains
operators ${\hat{\sigma}_{1}}$, ${\hat{\sigma}_{2}}$,
$\hat{\sigma}_{1}^{\dagger}{\hat{\sigma}_{2}}$,
${\hat{\sigma}_{1}}{\hat{\sigma}_{2}}$,
$\hat{\sigma}_{1}^{\dagger}{\hat{\sigma}_{1}}{\hat{\sigma}_{2}}$, and
$\hat{\sigma}_{2}^{\dagger}{\hat{\sigma}_{2}}{\hat{\sigma}_{1}}$, their
Hermitian conjugation, and Hermitian operators
$\hat{\sigma}_{1}^{\dagger}{\hat{\sigma}_{1}}$,
$\hat{\sigma}_{2}^{\dagger}{\hat{\sigma}_{2}}$,
$\hat{\sigma}_{1}^{\dagger}{\hat{\sigma}_{1}}\hat{\sigma}_{2}^{\dagger}{\hat{\sigma}_{2}}$,
$\hat{1}$. After straightforward calculations, we obtain

$\begin{array}[]{c}{\hat{\tilde{\sigma}}_{1}}(t)+{\hat{\tilde{\sigma}}_{2}}(t)=\\\
\left({{{\hat{\sigma}}_{1}}+{{\hat{\sigma}}_{2}}-\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}-\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}{{\hat{\sigma}}_{1}}}\right){e^{-i\left({\omega+\Omega}\right)t}}+\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}{{\hat{\sigma}}_{1}}}\right){e^{-i\left({\omega-\Omega}\right)t}},\end{array}$
(15)

$\begin{array}[]{l}\hat{\tilde{\sigma}}_{1}^{\dagger}(t){{\hat{\tilde{\sigma}}}_{1}}(t)=\frac{1}{2}\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}}\right)+\frac{{\hat{\sigma}_{1}^{\dagger}+\hat{\sigma}_{2}^{\dagger}}}{2}\cdot\frac{{{{\hat{\sigma}}_{1}}-{{\hat{\sigma}}_{2}}}}{2}{e^{i2\Omega
t}}+\frac{{\hat{\sigma}_{1}^{\dagger}-\hat{\sigma}_{2}^{\dagger}}}{2}\cdot\frac{{{{\hat{\sigma}}_{1}}+{{\hat{\sigma}}_{2}}}}{2}{e^{-i2\Omega
t}}.\end{array}$ (16)

(Alternatively, one can use the Baker–Campbell–Hausdorff formula Suzuki
(1985).) The Lindblad superoperators can be obtained via standard procedure
Breuer and Petruccione (2002); Davies (1974); Kosloff (2013), assuming that
the reservoirs are in thermal equilibrium. As a result, we obtain the
following Lindblad equation

$\begin{array}[]{l}{{\dot{\hat{\rho}}}_{\rm{S}}}(t)=-i{\hbar^{-1}}[{{\hat{H}}_{\rm{S}}},{{\hat{\rho}}_{\rm{S}}}(t)]+\\\
\sum\limits_{\rm{k}=1,2}{\sum\limits_{\rm{j}=1,2,3}{\frac{{{\gamma_{{\rm{k}}{\rm{,dp}}}}\left({2\Omega{\theta_{\rm{j}}}}\right)}}{2}\left({2{{\hat{L}}_{{\rm{dp}}{\rm{,kj}}}}{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{dp}}{\rm{,kj}}}^{\dagger}-{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{dp}}{\rm{,kj}}}^{\dagger}{{\hat{L}}_{{\rm{dp}}{\rm{,kj}}}}-\hat{L}_{{\rm{dp}}{\rm{,kj}}}^{\dagger}{{\hat{L}}_{{\rm{dp}}{\rm{,kj}}}}{{\hat{\rho}}_{\rm{S}}}(t)}\right)}}+\\\
\sum\limits_{\rm{k}=1,2}{\frac{{{\gamma_{{\rm{rad}}}}\left({-\left({\omega+{{(-1)}^{{\rm{k}}-1}}\Omega}\right)}\right)}}{2}\left({2{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}-{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}-\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}{{\hat{\rho}}_{\rm{S}}}(t)}\right)}+\\\
\sum\limits_{\rm{k}=1,2}{\frac{{{\gamma_{{\rm{rad}}}}\left({\left({\omega+{{(-1)}^{{\rm{k}}-1}}\Omega}\right)}\right)}}{2}\left({2\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{\rho}}_{\rm{S}}}(t){{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}-{{\hat{\rho}}_{\rm{S}}}(t){{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}-{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{\rho}}_{\rm{S}}}(t)}\right)}{\kern
1.0pt},\end{array}$ (17)

where Lindblad superoperators have the form

$\begin{array}[]{l}{{\hat{L}}_{{\rm{rad}}{\rm{,1}}}}={{\hat{\sigma}}_{1}}+{{\hat{\sigma}}_{2}}-\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}-{{\hat{\sigma}}_{1}}\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}},\,\,{{\hat{L}}_{{\rm{rad}}{\rm{,2}}}}=\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}+{{\hat{\sigma}}_{1}}\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}},\\\
{{\hat{L}}_{{\rm{dp}}{\rm{,11}}}}=\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}/2+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}/2,\,\,{{\hat{L}}_{{\rm{dp}}{\rm{,12}}}}=\left({\hat{\sigma}_{1}^{\dagger}+\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}-{{\hat{\sigma}}_{2}}}\right)/4,\,\,{{\hat{L}}_{{\rm{dp}}{\rm{,13}}}}=\left({\hat{\sigma}_{1}^{\dagger}-\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}+{{\hat{\sigma}}_{2}}}\right)/4,\\\
{{\hat{L}}_{{\rm{dp}}{\rm{,21}}}}=\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}/2+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}/2,\,\,{{\hat{L}}_{{\rm{dp}}{\rm{,22}}}}=\left({\hat{\sigma}_{1}^{\dagger}+\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{2}}-{{\hat{\sigma}}_{1}}}\right)/4,\,\,{{\hat{L}}_{{\rm{dp}}{\rm{,23}}}}=\left({\hat{\sigma}_{2}^{\dagger}-\hat{\sigma}_{1}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}+{{\hat{\sigma}}_{2}}}\right)/4,\end{array}$
(18)

and parameter ${\theta_{\rm{j}}}=0,\,\,1,\,\,-1$ for ${\rm{j}}=1,\,\,2,\,\,3$.
The ${\gamma_{m}}\left(\omega\right)$ functions are determined according to
${\gamma_{m}}\left(\omega\right)=\int\limits_{-\infty}^{\infty}{d\tau\exp\left({-i\omega\tau}\right)\left\langle{{{\hat{\tilde{R}}}_{m}}(t+\tau){{\hat{\tilde{R}}}_{m}}(t)}\right\rangle}$,
and, in the case of reservoirs which are in thermal equilibrium at temperature
$T$, they satisfy the Kubo-Martin-Shwinger (KMS) condition
${\gamma_{m}}\left(\omega\right)=\exp\left({-\hbar\omega/kT}\right){\gamma_{m}}\left({-\omega}\right)$.

If ${\omega_{1}}\neq{\omega_{2}}$ the states $\left|{ee}\right\rangle$ and
$\left|{gg}\right\rangle$ remain the system eigenstates. Indeed,
${\hat{H}_{\rm{S}}}\left|{ee}\right\rangle=\left({{\omega_{1}}+{\omega_{2}}}\right)\left|{ee}\right\rangle$
and ${\hat{H}_{\rm{S}}}\left|{gg}\right\rangle=0\left|{gg}\right\rangle$.
Other two eigenstates have the form

$\left|+\right\rangle=\left({\left|s\right\rangle+\frac{\Delta}{{2\Omega+\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}\left|{as}\right\rangle}\right)/\sqrt{1+{{\left({\frac{\Delta}{{2\Omega+\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}}\right)}^{2}}},$
(19)
$\left|-\right\rangle=\left({\left|{as}\right\rangle+\frac{{2\Omega-\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{\Delta}\left|s\right\rangle}\right)/\sqrt{1+{{\left({\frac{{2\Omega-\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{\Delta}}\right)}^{2}}},$
(20)

with eigenvalues

${\omega_{+}}=\frac{{\theta+\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2},\ \ \ \ \ \
\ \ {\omega_{-}}=\frac{{\theta-\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2}$ (21)

respectively. Here $\theta={\omega_{1}}+{\omega_{2}}$,
$\Delta={\omega_{1}}-{\omega_{2}}$.

To obtain new Lindblad equation, we, first, need to calculate
${\hat{\tilde{\sigma}}_{1,2}}(t)=\exp\left({i{{\hat{H}}_{\rm{S}}}t}\right){\hat{\sigma}_{1,2}}(t)\exp\left({-i{{\hat{H}}_{\rm{S}}}t}\right)$.
After straightforward calculations, we obtain

$\begin{array}[]{*{20}{c}}{{{\hat{\tilde{\sigma}}}_{1}}(t)}&=&{\frac{{{{\hat{\sigma}}_{1}}\frac{{y+2\sqrt{{y^{2}}/4+1}}}{2}+{{\hat{\sigma}}_{2}}}}{{2\sqrt{{y^{2}}/4+1}}}\exp\left({i\frac{{-\theta-\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2}t}\right)-\frac{{\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}}}{{\sqrt{{y^{2}}/4+1}}}\exp\left({i\frac{{-\theta-\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2}t}\right)}\\\
{}\hfil&{}\hfil&{+\frac{{-{{\hat{\sigma}}_{1}}\frac{{y-2\sqrt{{y^{2}}/4+1}}}{2}-{{\hat{\sigma}}_{2}}}}{{2\sqrt{{y^{2}}/4+1}}}\exp\left({i\frac{{-\theta+\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2}t}\right)+\frac{{\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}}}{{\sqrt{{y^{2}}/4+1}}}\exp\left({i\frac{{-\theta+\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2}t}\right),}\end{array}$
(22)

where $y=\left({{\omega_{1}}-{\omega_{2}}}\right)/\Omega$. Analogous
expression, up to replacement $1\leftrightarrow 2$ and $y\to-y$, is valid for
${\hat{\tilde{\sigma}}_{2}}(t)$. Using these expressions, we obtain

$\begin{array}[]{*{20}{c}}{{{\hat{\tilde{\sigma}}}_{1}}(t)+{{\hat{\tilde{\sigma}}}_{2}}(t)+\hat{\tilde{\sigma}}_{1}^{\dagger}(t)+\hat{\tilde{\sigma}}_{2}^{\dagger}(t)=}\\\
{\frac{1}{{\sqrt{{y^{2}}+4}}}\left[{\left({{{\hat{\sigma}}_{1}}\left({1+\frac{{y+\sqrt{{y^{2}}+4}}}{2}}\right)+{{\hat{\sigma}}_{2}}\left({1+\frac{{-y+\sqrt{{y^{2}}+4}}}{2}}\right)-2\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}-2\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}{{\hat{\sigma}}_{1}}}\right)\exp\left({i\frac{{-\theta-\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2}t}\right)+}\right.}\\\
{+\left({{{\hat{\sigma}}_{1}}\left({-1-\frac{{y-\sqrt{{y^{2}}+4}}}{2}}\right)+{{\hat{\sigma}}_{2}}\left({-1+\frac{{y+\sqrt{{y^{2}}+4}}}{2}}\right)+2\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}+2\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}{{\hat{\sigma}}_{1}}}\right)\exp\left({i\frac{{-\theta+\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}}{2}t}\right)+}\\\
{\left.{+h.\;c.}\right],}\end{array}$ (23)

$\begin{array}[]{*{20}{c}}{\hat{\tilde{\sigma}}_{1}^{\dagger}(t){{\hat{\tilde{\sigma}}}_{1}}(t)=}{\frac{1}{{{y^{2}}+4}}\left[{2\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}\frac{{{y^{2}}+2}}{2}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}}\right)+y\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{2}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{1}}}\right)}\right.}\\\
{-\left({\hat{\sigma}_{1}^{\dagger}\frac{{y+\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}\frac{{y-\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{2}}}\right)\exp\left({i\sqrt{{\Delta^{2}}+4{\Omega^{2}}}t}\right)}\\\
{-\left.{\left({\hat{\sigma}_{1}^{\dagger}\frac{{y-\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}\frac{{y+\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{2}}}\right)\exp\left({-i\sqrt{{\Delta^{2}}+4{\Omega^{2}}}t}\right)}\right],}\end{array}$
(24)

$\begin{array}[]{*{20}{c}}{\hat{\tilde{\sigma}}_{2}^{\dagger}(t){{\hat{\tilde{\sigma}}}_{2}}(t)=}{\frac{1}{{{y^{2}}+4}}\left[{2\left({\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}\frac{{{y^{2}}+2}}{2}+\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}}\right)-y\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{2}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{1}}}\right)}\right.}\\\
{-\left({\hat{\sigma}_{2}^{\dagger}\frac{{-y+\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{1}^{\dagger}}\right)\left({{{\hat{\sigma}}_{2}}\frac{{-y-\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{1}}}\right)\exp\left({i\sqrt{{\Delta^{2}}+4{\Omega^{2}}}t}\right)}\\\
{-\left.{\left({\hat{\sigma}_{2}^{\dagger}\frac{{-y-\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{1}^{\dagger}}\right)\left({{{\hat{\sigma}}_{2}}\frac{{-y+\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{1}}}\right)\exp\left({-i\sqrt{{\Delta^{2}}+4{\Omega^{2}}}t}\right)}\right].}\end{array}$
(25)

Using standard procedure for derivation of Linblad superoperators, one can get

$\begin{array}[]{l}{{\dot{\hat{\rho}}}_{\rm{S}}}(t)=-i{\hbar^{-1}}[{{\hat{H}}_{\rm{S}}},{{\hat{\rho}}_{\rm{S}}}(t)]+\\\
+\sum\limits_{{\rm{k}}=1,2}{\sum\limits_{{\rm{j}}=1,2,3}{\frac{{{\gamma_{{\rm{k}}{\rm{,dp}}}}\left({\sqrt{{\Delta^{2}}+4{\Omega^{2}}}{\theta_{\rm{j}}}}\right)}}{2}\left({2{{\hat{L}}_{{\rm{dp}}{\rm{,kj}}}}{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{dp}}{\rm{,kj}}}^{\dagger}-{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{dp}}{\rm{,kj}}}^{\dagger}{{\hat{L}}_{{\rm{dp}}{\rm{,kj}}}}-\hat{L}_{{\rm{dp}}{\rm{,kj}}}^{\dagger}{{\hat{L}}_{{\rm{dp}}{\rm{,kj}}}}{{\hat{\rho}}_{\rm{S}}}(t)}\right)}}+\\\
+\sum\limits_{{\rm{k}}=1,2}{\frac{{{\gamma_{{\rm{rad}}}}\left({-\left({\theta+{{(-1)}^{{\rm{k}}-1}}\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}\right)/2}\right)}}{2}\left({2{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}-{{\hat{\rho}}_{\rm{S}}}(t)\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}-\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}{{\hat{\rho}}_{\rm{S}}}(t)}\right)}+\\\
+\sum\limits_{{\rm{k}}=1,2}{\frac{{{\gamma_{{\rm{rad}}}}\left({\left({\theta+{{(-1)}^{{\rm{k}}-1}}\sqrt{{\Delta^{2}}+4{\Omega^{2}}}}\right)/2}\right)}}{2}\left({2\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{\rho}}_{\rm{S}}}(t){{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}-{{\hat{\rho}}_{\rm{S}}}(t){{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}-{{\hat{L}}_{{\rm{rad}}{\rm{,k}}}}\hat{L}_{{\rm{rad}}{\rm{,k}}}^{\dagger}{{\hat{\rho}}_{\rm{S}}}(t)}\right)}{\kern
1.0pt},\end{array}$ (26)

where Lindblad superoperators have the form

$\begin{array}[]{l}{{\hat{L}}_{{\rm{rad}}{\rm{,1}}}}=\frac{1}{{\sqrt{{y^{2}}+4}}}\left({{{\hat{\sigma}}_{1}}\left({1+\frac{{y+\sqrt{{y^{2}}+4}}}{2}}\right)+{{\hat{\sigma}}_{2}}\left({1+\frac{{-y+\sqrt{{y^{2}}+4}}}{2}}\right)-2\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}-2\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}{{\hat{\sigma}}_{1}}}\right),\\\
{{\hat{L}}_{{\rm{rad}}{\rm{,2}}}}=\frac{1}{{\sqrt{{y^{2}}+4}}}\left({{{\hat{\sigma}}_{1}}\left({-1-\frac{{y-\sqrt{{y^{2}}+4}}}{2}}\right)+{{\hat{\sigma}}_{2}}\left({-1+\frac{{y+\sqrt{{y^{2}}+4}}}{2}}\right)+2\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}{{\hat{\sigma}}_{2}}+2\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}{{\hat{\sigma}}_{1}}}\right),\end{array}$
(27)

$\begin{array}[]{l}{{\hat{L}}_{{\rm{dp}}{\rm{,11}}}}=\frac{1}{{{y^{2}}+4}}\left({2\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}\frac{{{y^{2}}+2}}{2}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}}\right)+y\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{2}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{1}}}\right)}\right),\\\
\,{{\hat{L}}_{{\rm{dp}}{\rm{,12}}}}=-\frac{1}{{{y^{2}}+4}}\left({\hat{\sigma}_{1}^{\dagger}\frac{{y+\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}\frac{{y-\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{2}}}\right),\\\
{{\hat{L}}_{{\rm{dp}}{\rm{,13}}}}=-\frac{1}{{{y^{2}}+4}}\left({\hat{\sigma}_{1}^{\dagger}\frac{{y-\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{2}^{\dagger}}\right)\left({{{\hat{\sigma}}_{1}}\frac{{y+\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{2}}}\right),\end{array}$
(28)

$\begin{array}[]{l}{{\hat{L}}_{{\rm{dp}}{\rm{,21}}}}=\frac{1}{{{y^{2}}+4}}\left({2\left({\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{2}}\frac{{{y^{2}}+2}}{2}+\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{1}}}\right)-y\left({\hat{\sigma}_{1}^{\dagger}{{\hat{\sigma}}_{2}}+\hat{\sigma}_{2}^{\dagger}{{\hat{\sigma}}_{1}}}\right)}\right),\\\
\,{{\hat{L}}_{{\rm{dp}}{\rm{,22}}}}=-\frac{1}{{{y^{2}}+4}}\left({\hat{\sigma}_{2}^{\dagger}\frac{{-y+\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{1}^{\dagger}}\right)\left({{{\hat{\sigma}}_{2}}\frac{{-y-\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{1}}}\right),\\\
{{\hat{L}}_{{\rm{dp}}{\rm{,23}}}}=-\frac{1}{{{y^{2}}+4}}\left({\hat{\sigma}_{2}^{\dagger}\frac{{-y-\sqrt{{y^{2}}+4}}}{2}+\hat{\sigma}_{1}^{\dagger}}\right)\left({{{\hat{\sigma}}_{2}}\frac{{-y+\sqrt{{y^{2}}+4}}}{2}+{{\hat{\sigma}}_{1}}}\right).\end{array}$
(29)

## Funding

Russian Science Foundation (20-72-10057).

## Acknowledgments

The study was supported by a grant from Russian Science Foundation (project
No. 20-72-10057). E.S.A. thank foundation for the advancement of theoretical
physics and mathematics “Basis”.

## References

  * Krantz  _et al._ (2019) P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Applied Physics Reviews 6, 021318 (2019).
  * Bennett  _et al._ (1996) C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Physical Review A 54, 3824 (1996).
  * Kumar  _et al._ (2002) P. Kumar, G. M. D’Ariano, and O. Hirota,  _Quantum Communication, Computing and Measurement 2_ (Springer, 2002).
  * Wootters (1998) W. K. Wootters, Physical Review Letters 80, 2245 (1998).
  * Hill and Wootters (1997) S. Hill and W. K. Wootters, Physical Review Letters 78, 5022 (1997).
  * Aolita  _et al._ (2015) L. Aolita, F. de Melo, and L. Davidovich, Reports on Progress in Physics 78, 042001 (2015).
  * Blanton  _et al._ (1997) S. A. Blanton, R. L. Leheny, M. A. Hines, and P. Guyot-Sionnest, Physical Review Letters 79, 865 (1997).
  * Shim and Guyot-Sionnest (1999) M. Shim and P. Guyot-Sionnest, The Journal of Chemical Physics 111, 6955 (1999).
  * Barenco (1995) A. Barenco, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 449, 679 (1995).
  * Poyatos  _et al._ (1997) J. F. Poyatos, J. I. Cirac, and P. Zoller, Physical Review Letters 78, 390 (1997).
  * Benhelm  _et al._ (2008) J. Benhelm, G. Kirchmair, C. F. Roos, and R. Blatt, Nature Physics 4, 463 (2008).
  * Zheng and Guo (2000) S.-B. Zheng and G.-C. Guo, Physical Review Letters 85, 2392 (2000).
  * Nielsen and Chuang (2010) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, Cambridge, 2010).
  * Braun (2002) D. Braun, Physical Review Letters 89, 277901 (2002).
  * Ficek and Tanaś (2002) Z. Ficek and R. Tanaś, Physics Reports 372, 369 (2002).
  * Ficek and Tanaś (2006) Z. Ficek and R. Tanaś, Physical Review A 74, 024304 (2006).
  * Scala  _et al._ (2008) M. Scala, R. Migliore, and A. Messina, Journal of Physics A: Mathematical and Theoretical 41, 435304 (2008).
  * Benatti  _et al._ (2010) F. Benatti, R. Floreanini, and U. Marzolino, Physical Review A 81, 012105 (2010).
  * Campagnano  _et al._ (2010) G. Campagnano, A. Hamma, and U. Weiss, Physics Letters A 374, 416 (2010).
  * Scala  _et al._ (2011) M. Scala, R. Migliore, A. Messina, and L. Sánchez-Soto, The European Physical Journal D 61, 199 (2011).
  * Li and Paraoanu (2009) J. Li and G. Paraoanu, New Journal of Physics 11, 113020 (2009).
  * Cecoi  _et al._ (2018) E. Cecoi, V. Ciornea, A. Isar, and M. A. Macovei, JOSA B 35, 1127 (2018).
  * Hor-Meyll  _et al._ (2009) M. Hor-Meyll, A. Auyuanet, C. Borges, A. Aragão, J. Huguenin, A. Khoury, and L. Davidovich, Physical Review A 80, 042327 (2009).
  * Facchinetti  _et al._ (2016) G. Facchinetti, S. D. Jenkins, and J. Ruostekoski, Physical Review Letters 117, 243601 (2016).
  * Kalachev and Samartsev (2005) A. A. Kalachev and V. V. Samartsev, Quantum Electronics 35, 679 (2005).
  * Plankensteiner  _et al._ (2015) D. Plankensteiner, L. Ostermann, H. Ritsch, and C. Genes, Scientific reports 5, 16231 (2015).
  * Kalachev (2007) A. Kalachev, Physical Review A 76, 043812 (2007).
  * Scully (2015) M. O. Scully, Physical Review Letters 115, 243602 (2015).
  * Gonzalez-Ballestero  _et al._ (2014) C. Gonzalez-Ballestero, E. Moreno, and F. Garcia-Vidal, Physical Review A 89, 042328 (2014).
  * Gonzalez-Tudela  _et al._ (2011) A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, Physical Review Letters 106, 020501 (2011).
  * Martin-Cano  _et al._ (2011) D. Martin-Cano, A. González-Tudela, L. Martín-Moreno, F. Garcia-Vidal, C. Tejedor, and E. Moreno, Physical Review B 84, 235306 (2011).
  * Gonzalez-Ballestero  _et al._ (2015) C. Gonzalez-Ballestero, A. Gonzalez-Tudela, F. J. Garcia-Vidal, and E. Moreno, Physical Review B 92, 155304 (2015).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Carmichael (2009) H. Carmichael,  _An open systems approach to quantum optics: lectures presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_ , Vol. 18 (Springer Science & Business Media, 2009).
  * Shishkov  _et al._ (2019) V. Y. Shishkov, E. S. Andrianov, A. A. Pukhov, A. P. Vinogradov, and A. A. Lisyansky, Physics-Uspekhi 62, 510 (2019).
  * Kosloff (2013) R. Kosloff, Entropy 15, 2100 (2013).
  * Shishkov  _et al._ (2020) V. Y. Shishkov, E. Andrianov, A. Pukhov, A. Vinogradov, and A. Lisyansky, Physical Review A 102, 032207 (2020).
  * Santos and Semiao (2014) J. Santos and F. Semiao, Physical Review A 89, 022128 (2014).
  * De Chiara  _et al._ (2018) G. De Chiara, G. Landi, A. Hewgill, B. Reid, A. Ferraro, A. J. Roncaglia, and M. Antezza, New Journal of Physics 20, 113024 (2018).
  * Cattaneo  _et al._ (2019) M. Cattaneo, G. L. Giorgi, S. Maniscalco, and R. Zambrini, New Journal of Physics 21, 113045 (2019).
  * Scully and Zubairy (1999) M. O. Scully and M. S. Zubairy, “Quantum optics,” (1999).
  * Agarwal (1974) G. S. Agarwal, in  _Quantum Optics_ (Springer, 1974) pp. 1–128.
  * Davies (1974) E. B. Davies, Communications in mathematical Physics 39, 91 (1974).
  * Gross and Haroche (1982) M. Gross and S. Haroche, Physics Reports 93, 301 (1982).
  * Yu (2016) N. Yu, Physical Review A 94, 060101 (2016).
  * Wolfe and Yelin (2014) E. Wolfe and S. F. Yelin, Physical Review Letters 112, 140402 (2014).
  * Macieszczak  _et al._ (2016) K. Macieszczak, M. Guţă, I. Lesanovsky, and J. P. Garrahan, Physical Review Letters 116, 240404 (2016).
  * Kimble (2008) H. J. Kimble, Nature 453, 1023 (2008).
  * Laurat  _et al._ (2007) J. Laurat, K. Choi, H. Deng, C. Chou, and H. Kimble, Physical Review Letters 99, 180504 (2007).
  * Jones  _et al._ (2018) D. E. Jones, B. T. Kirby, and M. Brodsky, npj Quantum Information 4, 1 (2018).
  * Weiss  _et al._ (2019) P. Weiss, A. Cipris, M. Araújo, R. Kaiser, and W. Guerin, Physical Review A 100, 033833 (2019).
  * Suzuki (1985) M. Suzuki, Journal of mathematical physics 26, 601 (1985).