# Extraction of autonomous quantum coherences

Artur Slobodeniuk<EMAIL_ADDRESS>Department of Optics,
Palacký University, 17. listopadu 12, 77146 Olomouc, Czech Republic
Department of Condensed Matter Physics, Faculty of Mathematics and Physics,
Charles University,
Ke Karlovu 5, CZ-121 16 Prague, Czech Republic  Tomáš Novotný
<EMAIL_ADDRESS>Department of Condensed Matter Physics, Faculty of
Mathematics and Physics, Charles University,
Ke Karlovu 5, CZ-121 16 Prague, Czech Republic  Radim Filip
<EMAIL_ADDRESS>Department of Optics, Palacký University, 17. listopadu
12, 77146 Olomouc, Czech Republic

###### Abstract

Quantum coherence is an essential resource to gain advantage over classical
physics and technology. Recently, it has been proposed that a low-temperature
environment can induce quantum coherence of a spin without an external
coherent pump. We address a critical question if such coherence is extractable
by a weak coupling to an output system dynamically affecting back the spin-
environment coupling. Describing the entire mechanism, we prove that such
extraction is generically possible for output spins (also oscillators or
fields) and, as well, in a fermionic analogue of such a process. We compare
the internal spin coherence and output coherence over temperature and
characteristic frequencies. The proposed optimal coherence extraction opens
paths for the upcoming experimental tests with atomic and solid-state systems.

##  1 Introduction

Quantum coherence of atomic and solid-state spins is at the core of many
quantum technology applications [1] and necessary to gain quantum advantage
for modern applications over classical technology. It is highly required for
quantum sensing [2, 3], quantum communication [4, 5], quantum simulators [6,
7, 8], energy harvesting [9, 10], quantum thermodynamics [11, 12, 13] and, of
course, in quantum computing [14, 15]. It also stimulates a broad recent
theoretical analysis [16, 17, 18, 19, 20, 21, 22] and new discussions [23, 24,
25, 26, 27, 28, 29, 30, 31, 32]. The simplest quantum coherence can appear in
a single spin as superposition of free-Hamiltonian eigenstates. A classical
coherent external force can be used to induce spin quantum coherence. On the
other hand, coupling of spin to the surrounding environment typically destroys
quantum coherence [33].

Can such a single spin, and later a larger spin ensemble, receive a steady-
state quantum coherence without any external coherent drive, purely from a
coherent coupling with its environment? If that is possible, quantum
technology can be, in principle, more autonomous. Maybe, as broadly discussed,
even Nature can already use such coherence to gain quantum advantages without
external drive [34, 35]. It might be hard to imagine that a low-temperature
environment coherently builds up quantum coherence. However, the recent
theoretical considerations predict such an open door. Indeed, single spin
coherence can generically appear by cooling the phononic environment to a low
temperature [36]. Quantum coherence of the spin rises there from the coherence
of engineered interaction with the environment.

Subsequent studies confirmed similar phenomena [37, 38, 39] and also bring the
first experimental proposal in the solid-state systems [40]. The path to the
first experiments seems to be open. However, an essential aspect remains still
unsolved. Is such autonomous coherence extractable, i.e. can an output system
continuously and even weakly connected to the spin extract the steady-state
coherence? Would such weak extraction not disturb the coherence generation?
These essential questions are already critical for a fundamental understanding
and detectability of quantum coherence as a resource. Simultaneously, to
answer them will be necessary for more autonomous quantum devices and their
applications. Apparently, our current picture about an extraction of coherence
caused by an external drive is not sufficient. It assumes that strong coherent
external drive is largely immune to coupling the output system. As the
autonomous coherences appear by an interaction with environment, their
extraction is still an open issue.

The existing experimental proposal [40] has opened this key question of
extraction of autonomous coherence into an intracavity field. The extraction
assumed a weak linear coupling of the cavity field to the spin in a simple
static approximation without any back action on the spin and coupling with the
environment. For the cavity-field frequency much below the spin frequency and
mean spin population almost unchanged (see equation (70) in Supplementary
Notes [40]), the spin coherence drives the mean coherent displacement of the
high-quality cavity field. It scales only by the spin-cavity coupling strength
(see equations (70)-(73) in Supplementary Notes [40]) and is inversely
proportional to the field frequency.

Here, we build upon that analysis and propose a minimal extraction mechanism
suitable for fields, oscillators, spins, and fermions. We compare that static
approximation to our full dynamical description with an essential back action
of extraction on spin-environment interaction. In this way, we thoroughly
answer the question about autonomous coherence extractability. To provide an
in-depth analysis, we adopt the spin as the output system, as illustrated in
Fig. 1 (a). However, as the spin is weakly excited, the results remain valid
for the oscillator of field mode as the output system. For the sake of
completeness, we present an adequate extension to fermionic systems.
Importantly, we confirm that the autonomous coherence is generically
extractable to all these output systems. It paves a road to the first
experimental tests. In the full dynamical description, extracted coherence
nontrivially depends on the environment temperature and characteristic system
frequencies. For the spin system, we reach coherence maximum asymptotically
for a low temperature. Interestingly, the coherence maximum appears at finite
temperature for a fermionic analogue. It is stimulating for further proposals
to obtain autonomous quantum coherences.

##  2 Spin-boson model

Figure 1:  Extraction of autonomous quantum coherence. (a) Model of the
system. The internal spin $1$ with a coherence $\langle\sigma_{1}^{x}\rangle$
interacts with the bosonic degrees of freedom of the environmental bath,
represented by the red two-headed arrows. The intensity of interaction with
the $k$-th bosonic mode is characterized by the weight parameter
$\lambda_{k}$. The coupling of the internal spin with the bath is defined by
the pair of real parameters $f_{1},f_{2}$. The internal spin interacts with
the output spin $2$ with a coherence $\langle\sigma_{2}^{x}\rangle$ via spin-
spin coupling term. The strength of coupling is defined by the parameter $t$.
The bath temperature is $T$. (b) Scheme of derivation the coherence of the
output spin $\langle\sigma_{2}^{x}\rangle$ in the static response
approximation. The internal spin with the non-zero thermal average
$\langle\sigma_{1}^{z}\rangle$ but vanishing $\langle\sigma_{1}^{x}\rangle$
(first panel) induces the non-zero shift of the boson degrees of freedom
$\langle b_{k}+b_{k}^{\dagger}\rangle$ (second panel) via the $f_{1}$ term.
The bosonic shift acts back on to the output spin via the $f_{2}$ term and
generates a non-zero coherence of the internal spin
$\langle\sigma_{1}^{x}\rangle$ (third panel). The part of this coherence is
transferred to the output spin via the $t$ term and generates non-zero
coherence $\langle\sigma_{2}^{x}\rangle$ (fourth panel). (c) The sketch of the
temperature dependence of the extracted coherence
$\langle\sigma_{2}^{x}\rangle$ derived in the static response approximation
(blue curve) and in the quantum full dynamical response (gray curve).

To analyze the extraction of autonomous quantum coherence, we minimally extend
the model [36, 40]. Fig. 1 (a) schematically illustrates the model we use to
describe the extraction of quantum coherence. The Hamiltonian of the model
$H=H_{12}+H_{B}+H_{B1}$ consists of three parts, described below. We consider
a single (internal) spin $1$ coupled weakly with a strength $t$ to an output
spin $2$. The output spin $2$ will extract the quantum coherence of the
initial spin $1$. The overall two-spin Hamiltonian

$H_{12}=\frac{\hbar\omega_{1}}{2}\sigma_{1}^{z}+\frac{\hbar\omega_{2}}{2}\sigma_{2}^{z}+\hbar
t(\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{2}^{+}\sigma_{1}^{-}),$ (1)

describes the spin-spin coupling, where $\hbar\omega_{1}$ and
$\hbar\omega_{2}$ are the energies of the first and second spins,
respectively. Here the Pauli spin operators $\sigma_{1}^{j}$ and
$\sigma_{2}^{j}$ (with $j=x,y,z$) correspond to the spins $1$ and $2$. The
raising $\sigma_{n}^{+}$ and lowering $\sigma_{n}^{-}$ operators of
$n(=1,2)$-th spin are $\sigma_{n}^{\pm}=(\sigma_{n}^{x}\pm
i\sigma_{n}^{y})/2$. Note that, the coupling term in Eq. (1) is not able to
generate the output coherence if it is not present in the internal spin. As
the coupling is weak $t\ll\omega_{1},\omega_{2}$ and the internal spin $1$ is
weakly excited, any excitation of the output spin $2$ will be small too. In
this limit, we can equivalently consider the output spin $2$ as the two lowest
levels of an oscillator or field of radiation. It generalizes the description
and covers detailed analysis also for the proposal [40]. Numerical simulations
check the validity of weak coupling approximation for small
$t<\omega_{1},\omega_{2}$, see Appendices H and I for details.

As in the previous works [36, 40], the internal spin $1$ is coupled to a
multimode phononic bath (generally, bosonic bath) described by the Hamiltonian

$H_{B}=\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k},$ (2)

through the engineered spin-boson coupling

$H_{B1}=(f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x})\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}).$
(3)

The boson system is described by creation $b_{k}^{\dagger}$ and annihilation
$b_{k}$ operators and spectrum $\Omega_{k}>0$. Subscript $k$ parameterizes all
the boson degrees of freedom. The interaction between the internal spin $1$
and bosons are defined by real-valued coupling constants $f_{1}$ and $f_{2}$
[36, 40]. Steady-state equilibrium quantum coherence autonomously appears at a
low temperature if both $f_{1},f_{2}$ are non-vanishing [36, 40]. In Ref.
[36], the Hamiltonian without the output subsystem was considered, and hence
the extraction of autonomous quantum coherence remained an open problem. In
Ref. [40], a perturbative static analysis of the extraction predicted a
coherence of the output cavity-field coupled to the spin $1$. As is described
in Fig. 1 (b), such a static approximation does not consider the crucial back-
action of the output system on the spin-environment. Although, the static
approach is an inspiring approximation, the complete dynamical description is
needed to address extractability convincingly. In such a dynamical theory,
extractable quantum coherence can not only change, as depicted schematically
at Fig. 1 (c), but it could also potentially vanish, which would be critical.

In the full treatment of the perturbative analysis, quantum coherence of the
output spin $2$ in our spin model (see Appendices A, B, C, D, and E for
details) can reach

$\displaystyle\langle\sigma_{2}^{x}\rangle_{T}=\frac{4f_{1}f_{2}t}{\omega_{2}}\tanh\Big{(}\frac{\beta\hbar\omega_{2}}{2}\Big{)}\tanh\Big{(}\frac{\beta\hbar\omega_{1}}{2}\Big{)}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\Big{(}\frac{\beta\xi}{2}\Big{)}-\hbar\omega_{1}\coth\Big{(}\frac{\beta\hbar\omega_{1}}{2})}{\xi(\xi^{2}-\hbar^{2}\omega_{1}^{2})},$
(4)

where $\beta\equiv(k_{B}T)^{-1}$ is the inverse temperature in energy units,
$k_{B}$ is the Boltzmann constant, $T$ is the temperature of the system in
Kelvin scale and
$\mathcal{I}(\xi)\equiv\sum_{k}\lambda_{k}^{2}\delta(\xi-\Omega_{k})$ is the
bath spectral density function. Simultaneously, the mean value
$\langle\sigma_{2}^{y}\rangle_{T}$ is zero. The extracted quantum coherence is
proportional to $f_{1}$ and $f_{2}$ as in the previous analysis [40].
Similarly, it linearly scales with the coupling strength $t$ of the
interaction term in Eq. (1). The frequency dependence on $\omega_{1}$ and
$\omega_{2}$ is factorized. Note that the coherence of the second spin scales
linearly with the coherence of the first spin

$\langle\sigma_{1}^{x}\rangle_{T}=-4f_{1}f_{2}\tanh\Big{(}\frac{\beta\hbar\omega_{1}}{2}\Big{)}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\Big{(}\frac{\beta\xi}{2}\Big{)}-\hbar\omega_{1}\coth\Big{(}\frac{\beta\hbar\omega_{1}}{2}\Big{)}}{\xi\left(\xi^{2}-\hbar^{2}\omega_{1}^{2}\right)},$
(5)

derived in Appendix F and satisfies the relation

$\langle\sigma_{2}^{x}\rangle_{T}=-\frac{t}{\omega_{2}}\tanh\Big{(}\frac{\beta\hbar\omega_{2}}{2}\Big{)}\langle\sigma_{1}^{x}\rangle_{T},$
(6)

valid for any spectral density function $\mathcal{I}(\xi)$. This property can
be easily observed in the static approximation (see details in Appendix G),
where the coherence of the second spin is derived in three steps, depicted in
Fig. 1 (b). In the first step, supposing the smallness of the parameter $t$,
one considers the spins as non-interacting and estimates the non-zero thermal
average of the first spin
$\langle\sigma_{1}^{z}\rangle_{T}^{st}=-\tanh(\beta\hbar\omega_{1}/2)$ (first
panel in Fig. 1 (b)). Here, the superscript “st” means “static”. Then,
replacing the operator $\sigma_{1}^{z}$ by its average value
$\langle\sigma_{1}^{z}\rangle_{T}^{st}$, one gets the linear in
$(b_{k}+b_{k}^{\dagger})$ correction to the bosonic Hamiltonian. This
correction causes the non-zero average value $\langle
b_{k}+b_{k}^{\dagger}\rangle_{T}^{st}=-2f_{1}\langle\sigma_{1}^{z}\rangle_{T}^{st}\lambda_{k}/\Omega_{k}$
for the $k$-th bosonic mode (second panel in Fig. 1 (b)). In the second step,
the average $\langle b_{k}+b_{k}^{\dagger}\rangle_{T}^{st}$ generates the non-
zero coherence of the first spin
$\langle\sigma_{1}^{x}\rangle^{st}_{T}\approx-4f_{1}f_{2}[\langle\sigma_{1}^{z}\rangle_{T}^{st}]^{2}\Omega/(\hbar\omega_{1})$
due to the $f_{2}\sigma_{1}^{x}$ term in Hamiltonian (3) (third panel in Fig.
1 (b)). Here
$\Omega\equiv\sum_{k}\lambda_{k}^{2}/\Omega_{k}=\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)/\xi$.
Finally, $\langle\sigma_{1}^{x}\rangle_{T}^{st}$ induces the non-zero
coherence of the second spin
$\langle\sigma_{2}^{x}\rangle_{T}^{st}\approx-(t/\omega_{2})\tanh(\beta\hbar\omega_{2}/2)\langle\sigma_{1}^{x}\rangle_{T}^{st}$
due to the spin-spin interaction term in Hamiltonian (1) (fourth panel in Fig.
1 (b)). Note that the obtained result is valid only if
$t\ll|\omega_{1}-\omega_{2}|,\omega_{1},\omega_{2}$, since under these
conditions both spins can be considered separately and weakly interacting,
while the results (4) and (5) are valid also for $\omega_{1}=\omega_{2}$ case
(see details in Appendices H and I).

Both temperature and frequency dependency, however, needs more detailed
analysis beyond the static approximation. For the generic case of the ohmic
spectral density function $\mathcal{I}(\xi)=A\xi\exp(-\xi/\hbar\omega_{c})$
[41, 42], the coherence expression reads

$\langle\sigma_{2}^{x}\rangle_{T}=\frac{4Af_{1}f_{2}t}{\omega_{2}}\tanh\left(\frac{\beta\hbar\omega_{2}}{2}\right)\tanh\left(\frac{\beta\hbar\omega_{1}}{2}\right)\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1}\Big{)},$
(7)

where

$\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1}\Big{)}=\int_{0}^{\infty}dx\,e^{-\frac{\omega_{1}}{\omega_{c}}x}\,\frac{x\coth\left(\frac{\beta\hbar\omega_{1}x}{2}\right)-\coth\left(\frac{\beta\hbar\omega_{1}}{2}\right)}{x^{2}-1},$
(8)

is a dimensionless function of parameters $\omega_{1}/\omega_{c}$ and
$\beta\hbar\omega_{1}$. In the low-temperature limit $\beta\hbar\omega_{1}\gg
1$, quantum coherence reaches

$\langle\sigma_{2}^{x}\rangle_{T}\approx\frac{4Af_{1}f_{2}t}{\omega_{2}}\tanh\left(\frac{\beta\hbar\omega_{2}}{2}\right)\Big{[}-e^{\frac{\omega_{1}}{\omega_{c}}}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}-\frac{\pi^{2}}{3\beta^{2}\hbar^{2}\omega_{1}^{2}}\Big{]},$
(9)

where $\textrm{Ei}(x)=-\int_{-x}^{\infty}dz\,e^{-z}/z$ is the exponential
integral function [43].

Figure 2:  Normalized autonomous coherences in the spin-boson model. (a) Log-
log scale plot of normalized coherences
$\langle\sigma_{2}^{x}\rangle_{T}/\langle\sigma_{2}^{x}\rangle_{T=0}$ (solid
lines) and
$\langle\sigma_{1}^{x}\rangle_{T}/\langle\sigma_{1}^{x}\rangle_{T=0}$ (dashed
lines) for different values of $\omega_{1}/\omega_{c}=0.01,1,100$ (blue,
yellow and green curves, respectively) for the particular case of
$\omega_{2}/\omega_{1}=3$ as a function of temperature $T$. Here $\omega_{c}$
is the cut-off energy of the bosonic spectral function of the ohmic type. The
inset is the log-linear scale plot of the same functions with the same color
notation and parameters. It represents the low-temperature behavior of
normalized coherences. The extractable coherence of the spin $2$ starts to
increase fast for $k_{B}T<\hbar\omega_{1}$ and saturates sooner at the maximum
for smaller $\omega_{1}/\omega_{c}$. The backaction of the extraction process
does not slow coherence increase in the low-temperature region. (b) Log-log
scale plot of absolute value of dimensionless ratio
$\langle\sigma_{2}^{x}\rangle_{T}/\langle\sigma_{1}^{x}\rangle_{T=0}$
multiplied by the normalization factor $\omega_{1}/t$ for
$\omega_{1}/\omega_{c}=0.01$ and different values of
$\omega_{2}/\omega_{1}=0.5,2,4,8,16$ (blue, yellow, green, red and purple
curves, respectively) as a function of temperature $T$. It represents the
dependence of the coherence $\langle\sigma_{2}^{x}\rangle_{T}$ on
$\omega_{2}$. The gray dashed line depicts the asymptote of all curves, which
scales as $T^{-3}$. The inset shows the linear-log plot of the same functions
with the same color notation and parameters at small temperatures. It is
advantageous to reach $\omega_{2}/\omega_{1}<1$, to extract quantum coherence
to the spin $2$ more efficiently.

Both the internal spin $1$ and output spin $2$ coherences, normalized to their
asymptotic values at zero temperature, are compared in Fig. 2 (a). Both the
coherences increase monotonously with temperature reduction as in the static
approximation. They demonstrate nearly linear behavior at log-log scale for
the high temperature region $k_{B}T/\hbar\omega_{1}\gg 1$. Therefore, the
normalized coherences have a power-law dependence,

$\displaystyle\langle\sigma_{1}^{x}\rangle_{T\rightarrow\infty}$
$\displaystyle\approx-
f_{1}f_{2}\frac{\Omega}{3\hbar\omega_{1}}(\beta\hbar\omega_{1})^{2}=-f_{1}f_{2}A\frac{\omega_{c}}{3\omega_{1}}(\beta\hbar\omega_{1})^{2},$
(10) $\displaystyle\langle\sigma_{2}^{x}\rangle_{T\rightarrow\infty}$
$\displaystyle\approx\frac{f_{1}f_{2}t}{\omega_{1}}\frac{\Omega}{6\hbar\omega_{1}}(\beta\hbar\omega_{1})^{3}=\frac{Af_{1}f_{2}t}{\omega_{1}}\frac{\omega_{c}}{6\omega_{1}}(\beta\hbar\omega_{1})^{3},$
(11)

which, however, remain valid only for $k_{B}T\gg\hbar\omega_{1}$. Note that
the aforementioned power-law behavior remains valid for the case of an
arbitrary spectral density function $\mathcal{I}(\xi)$ if the condition
$\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)/\xi<+\infty$ is satisfied (see Eq.
(103) in Appendix E). In particular, it is valid for the generalized Ohmic
spectral density function
$\mathcal{I}(\xi,s)\sim\xi^{s}\exp(-\xi/\hbar\omega_{c})$ with $s>0$ (sub-
Ohmic for $0<s<1$ and super-Ohmic for $s>1$, respectively). Extracted spin $2$
coherence increases with a larger slope at the log-log scale. At the high-
temperature limit, autonomous coherence is relatively negligible. As the
temperature reduces to $k_{B}T\approx\hbar\omega_{1}$, both coherences
increase almost up to their maximum values and the normalized output coherence
gets closer to the internal one. The latter is evident from the inset of Fig.
2 (a). Note that a wider band of the spectral function with smaller
$\omega_{1}/\omega_{c}$ can already stimulate a faster increase of extractable
coherence. Reducing temperature further, the coherence saturates when
$k_{B}T/\hbar\omega_{1}\ll 1$ much earlier for $\omega_{1}/\omega_{c}\ll 1$.
It reaches maximum asymptotically at zero temperature. It is advantageous to
have $\omega_{1}/\omega_{c}$ as small as possible to reach a plateau near the
maximum of extractable coherence sooner by lowering the temperature. Indeed,
the Eq. (6) proves quantum coherence is extractable, expectably, scaled by the
coupling strength $t$. The extraction depends both on temperature $T$ and
exclusively on $\omega_{2}$. Frequency $\omega_{1}$ does not advantageously
affect the extraction process in this approximation. For large temperature,
the transfer coefficient in (6) turns to
$\langle\sigma_{2}^{x}\rangle_{T}=-\frac{\hbar
t}{2k_{B}T}\langle\sigma_{1}^{x}\rangle$. High temperatures thus limit the
yielding of the extraction process; on the other hand, for small temperatures,
the extraction follows
$\langle\sigma_{2}^{x}\rangle_{T}=-\frac{t}{\omega_{2}}\langle\sigma_{1}^{x}\rangle$
and extraction is not limited by the temperature. Lower $\omega_{2}$ is
preferable as it increases the extraction. It approaches the result for the
static response qualitatively [40].

To visualize dependence of the output spin $2$ coherence as a function of
temperature for different values of $\omega_{2}$ we consider the following
dimensionless ratio

$\frac{\omega_{1}}{t}\frac{\langle\sigma_{2}^{x}\rangle_{T}}{\langle\sigma_{1}^{x}\rangle_{T=0}}=\frac{\omega_{1}}{\omega_{2}}\tanh\left(\frac{\beta\hbar\omega_{2}}{2}\right)\tanh\left(\frac{\beta\hbar\omega_{1}}{2}\right)\frac{\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1}\Big{)}}{e^{\frac{\omega_{1}}{\omega_{c}}}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}}\propto\langle\sigma_{2}^{x}\rangle_{T}$
(12)

Note that for a fixed ratio $\omega_{1}/\omega_{c}$ the coherence of the
second spin at large temperatures becomes independent of $\omega_{2}$ and
demonstrates a universal behavior. In particular, this feature can be observed
from the plots of the function

$\Big{|}\frac{\omega_{1}}{t}\frac{\langle\sigma_{2}^{x}\rangle_{T}}{\langle\sigma_{1}^{x}\rangle_{T=0}}\Big{|}_{T\rightarrow\infty}\approx\frac{\omega_{c}}{24\omega_{1}}\frac{(\beta\hbar\omega_{1})^{3}}{e^{\frac{\omega_{1}}{\omega_{c}}}\Big{|}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}\Big{|}}\propto
T^{-3},$ (13)

for different $\omega_{2}$. The transfer of coherence from spin $1$ to spin
$2$ is visualized in Fig. 2 (b) for an optimal choice of small
$\omega_{1}/\omega_{c}$ and different $\omega_{2}/\omega_{1}$. For high
temperatures, $k_{B}T/\hbar\omega_{1}\gg 1$, the increase is linear in the
log-log plot with the scaling $T^{-3}$ as predicted above. The extraction is
inefficient, and the temperature reduction down to $k_{B}T/\hbar\omega_{1}\ll
1$ is needed to reach efficient transfer of the quantum coherence. A smaller
frequency $\omega_{2}$ of the output spin $2$ is very advantageous in this
low-temperature regime as longer as the perturbation approach holds. Transfer
to the output spin $2$ coherence is substantially more efficient when
$\omega_{2}<\omega_{1}$. Together with the previous condition, the optimal
extraction to the output spin (and oscillator or field) needs
$\omega_{2}<\omega_{1}<\omega_{c}$.

Note that the coherence of the spin 2 in the low-temperature limit
$\beta\hbar\omega_{1}\gg 1$ can be approximated as

$\displaystyle\langle\sigma_{2}^{x}\rangle_{T\rightarrow
0}\approx\frac{4f_{1}f_{2}t}{\hbar\omega_{1}\omega_{2}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{[}\frac{1}{\xi}-\frac{1}{\xi+\hbar\omega_{1}}\Big{]}=\frac{4f_{1}f_{2}t}{\hbar\omega_{1}\omega_{2}}\Omega-\frac{4f_{1}f_{2}t}{\hbar\omega_{1}\omega_{2}}\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)}{\xi+\hbar\omega_{1}}.$
(14)

The first term of this expression is nothing but the coherence of spin 2 in
the static approximation $\langle\sigma_{2}^{x}\rangle_{T}^{st}$ at zero
temperature. In this term the contribution of the bosonic bath is fully
factorized and represented only by the parameter $\Omega$. The factorization
reflects the fact that the back-action effects are not taken into account in
the static approximation. The result is not sensitive to the shape of the
spectral density function $\mathcal{I}(\xi)$, since two different spectral
density functions with the same $\Omega$ give the same value of the coherence
$\langle\sigma_{2}^{x}\rangle_{T}^{st}$. The second term can be called
dynamical since it provides the back-action correction to the static
contribution. Note that the frequency parameter $\omega_{1}$ of the spin 1 and
energy dependent spectral density function $\mathcal{I}(\xi)$ are included in
the expression inseparably in contrast to the static result. Therefore, the
dynamical term is sensitive to the shape of the spectral density function
$\mathcal{I}(\xi)$. Finally, the second term represents the deviation of the
full coherence $\langle\sigma_{2}^{x}\rangle_{T}$ form the static result
$\langle\sigma_{2}^{x}\rangle_{T}^{st}$ and hence its value as a function of
the temperature answers the question depicted in Fig. 1 (c).

Since $\mathcal{I}(\xi)>0$, the static approximation provides only the upper
bound of the extracted coherence of spin 2 at low temperatures. Therefore, the
correct estimation of such a coherence needs a better understanding of the
dynamical term. However, it seems that this term doesn’t have an intuitive
explanation in terms of the spin-boson model. Moreover, the perturbation
analysis of the dynamical contribution in the spin language is bulky due to
the non-commutativity of spin operators and the absence of the Wick theorem
for them. This makes the analysis of the static and dynamical terms in higher
orders of perturbation theory hard.

Therefore, beyond a general interest, it can be also mathematically valuable
to rewrite the spin Hamiltonian, for example, in terms of fermions where the
Wick theorem is applicable. It allows the use of the Green function technique
and, potentially, makes the final result clearer for understanding. This
analysis goes beyond the previous theory for the spins [36, 40], and helps to
compare the spin and fermionic systems coupled to the same bosonic baths.
Additionally, such an analysis can stimulate a broader class of experimental
tests and discussions about a coherence in fermionic systems.

##  3 Fermion-boson model

In order to implement the aforementioned fermionization program, we apply the
Jordan-Wigner transformation (see [44] and references therein) to the
Hamiltonians (1) and (3). First, we introduce auxiliary (zeroth) subsystem,
which is described by the operators
$\\{\sigma_{0}^{x},\sigma_{0}^{y},\sigma_{0}^{z}\\}$ and has zero energy. Then
we replace the spin operators by the fermionic ones $a_{n},n=0,1,2$, with the
standard anti-commutation rules $\\{a_{n},a_{m}^{\dagger}\\}=\delta_{nm}$

$\displaystyle\sigma_{0}^{-}$ $\displaystyle\rightarrow
a_{0},\quad\sigma_{0}^{+}\rightarrow
a_{0}^{\dagger},\quad\sigma_{0}^{z}\rightarrow 2a_{0}^{\dagger}a_{0}-1,$ (15)
$\displaystyle\sigma_{n}^{-}$
$\displaystyle\rightarrow\prod_{m=0}^{n-1}(2a_{m}^{\dagger}a_{m}-1)a_{n},\quad\sigma_{n}^{+}\rightarrow\prod_{m=0}^{n-1}(2a_{m}^{\dagger}a_{m}-1)a^{\dagger}_{n},\quad\sigma_{n}^{z}\rightarrow
2a_{n}^{\dagger}a_{n}-1,$ (16)

and make the additional substitutions

$\displaystyle\sigma_{1}^{z}=$ $\displaystyle
2a_{1}^{\dagger}a_{1}-1\rightarrow
2a_{1}^{\dagger}a_{1},\quad\sigma_{2}^{z}=2a_{2}^{\dagger}a_{2}-1\rightarrow
2a_{2}^{\dagger}a_{2},$ (17) $\displaystyle\sigma_{1}^{x}=$
$\displaystyle(\sigma_{1}^{+}+\sigma_{1}^{-})\rightarrow(\sigma_{1}^{+}\sigma_{0}^{-}+\sigma_{0}^{+}\sigma_{1}^{-})=-(a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1}),$
(18)

to remove the constant terms, which do not affect the quantum dynamics
significantly. Then the overall fermion-boson Hamiltonian takes the following
form

$\displaystyle H_{F}=$
$\displaystyle\hbar\omega_{1}a_{1}^{\dagger}a_{1}+\hbar\omega_{2}a_{2}^{\dagger}a_{2}-\hbar
t(a_{1}^{\dagger}a_{2}+a_{2}^{\dagger}a_{1})+\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}+$
(19)
$\displaystyle+[2f_{1}a_{1}^{\dagger}a_{1}-f_{2}(a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1})]\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})$

from the sum of spin-boson Hamiltonians (1), (2) and (3). Note that the
fermionization procedure changes a sign of the spin-spin $t\rightarrow-t$ and
spin-boson $f_{2}\rightarrow-f_{2}$ interaction terms. In the following, we
will use the subscript $F$ to indicate that we derived the coherences from the
fermionized Hamiltonian $H_{F}$.

Considering same ohmic spectral characteristics of bosonic bath as for the
previous two-spin analysis, we can define and derive the fermionic coherence
of the first spin using the Eq. (18) (see Appendices J and K for details)

$\langle\sigma_{1}^{x}\rangle_{F,T}\equiv-\langle
a_{1}^{\dagger}a_{0}+a_{1}^{\dagger}a_{0}\rangle_{F}=-4f_{1}f_{2}A\tanh\Big{(}\frac{\beta\hbar\omega_{1}}{2}\Big{)}\Big{[}-\frac{\omega_{c}}{2\omega_{1}}+\frac{1}{2}\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1}\Big{)}\Big{]}.$
(20)

Here $\langle\dots\rangle_{F}$ is the thermal average with the Hamiltonian
$H_{F}$. The final expression represents the leading contribution to the
coherence. We define the fermion coherence of the spin 2 by analogy with the
expression for the coherence of the spin 1 and obtain

$\langle\sigma_{2}^{x}\rangle_{F,T}\equiv-\langle
a_{2}^{\dagger}a_{0}+a_{0}^{\dagger}a_{2}\rangle_{F}=-\frac{4tf_{1}f_{2}A}{\hbar\Delta\omega}\mathcal{F}_{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1},\beta\hbar\omega_{2}\Big{)},$
(21)

where $\Delta\omega=\omega_{2}-\omega_{1}$ and

$\displaystyle\mathcal{F}_{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1},\beta\hbar\omega_{2}\Big{)}=\frac{\omega_{c}}{\omega_{2}}n(\hbar\omega_{1})[1-2n(\hbar\omega_{2})]+\int_{-\infty}^{\infty}dxe^{-\frac{\omega_{1}}{\omega_{c}}|x|}\mathcal{K}(x),$
(22)

with frequency-dependent kernel

$\displaystyle\mathcal{K}(x)=[N(\hbar\omega_{1}x)+n(\hbar\omega_{1})]\Big{\\{}\frac{x}{2(x-1)}\frac{\Delta\omega}{\omega_{2}}-\frac{n[\hbar\omega_{1}(1-x)]x}{(x-1)(x+\frac{\Delta\omega}{\omega_{1}})}\frac{\Delta\omega}{\omega_{1}}+\frac{n(\hbar\omega_{2})\omega_{1}x}{\omega_{2}(x+\frac{\Delta\omega}{\omega_{1}})}\Big{\\}}.$
(23)

We introduced the notations for the Bose-Einstein
$N(\xi)=(e^{\beta\xi}-1)^{-1}$ and the Fermi-Dirac
$n(\xi)=(e^{\beta\xi}+1)^{-1}$ distributions for brevity. Note that the
coherence of each spin is a sum of the Hartree and Fock contributions (see
Appendices L, M, N and O for details). It turns out that the Hartree
contribution of $j(=1,2)$-th spin coincides with the fermionic coherence of
the $j$-th spin obtained in the static approximation (see Appendix P)

$\langle\sigma_{1}^{x}\rangle_{F,T}^{st}\approx\frac{4f_{1}f_{2}\Omega}{\hbar\omega_{1}}n(\hbar\omega_{1})\tanh\Big{(}\frac{\beta\hbar\omega_{1}}{2}\Big{)},$
(24) $\langle\sigma_{2}^{x}\rangle_{F,T}^{st}\approx\frac{4f_{1}f_{2}\Omega
t}{\hbar\Delta\omega}n(\hbar\omega_{1})\left[\frac{\tanh\Big{(}\frac{\beta\hbar\omega_{1}}{2}\Big{)}}{\omega_{1}}-\frac{\tanh\Big{(}\frac{\beta\hbar\omega_{2}}{2}\Big{)}}{\omega_{2}}\right].$
(25)

Like in the previous spin-boson case, the influence of the boson bath is fully
factorized and represented by the same parameter $\Omega$ for these static
terms. The factorization has a natural explanation in the fermion picture.
Namely, it comes from the structure of the expression for the Hartree term
(see Eq. (168) in Appendix K). Then, the dynamical contribution to the fermion
coherences is defined by the Fock terms. Again, similarly to the spin-boson
model, the mutual influence of the bosonic bath and the fermion subsystem are
not separated in such a case. This property can be seen from the expression
for the Fock contribution (see Eq. (167) in Appendix K). The clear physical
meaning of the static and dynamical terms in the fermion system provides the
potential opportunity to understand better the corresponding terms in spin
systems via their fermionization.

The final relation between the internal and output coherence reads

$\displaystyle\langle\sigma_{2}^{x}\rangle_{F,T}=\frac{t}{\omega_{2}}\langle\sigma_{1}^{x}\rangle_{F,T}\left[-\frac{2\omega_{2}}{\Delta\omega}\coth\Big{(}\frac{\beta\hbar\omega_{1}}{2}\Big{)}\frac{\mathcal{F}_{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1},\beta\hbar\omega_{2}\Big{)}}{\frac{\omega_{c}}{\omega_{1}}-\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1}\Big{)}}\right]=\frac{t}{\omega_{2}}\langle\sigma_{1}^{x}\rangle_{F,T}\mathcal{R}(T).$
(26)

Note that the ratio $\mathcal{R}(T)$ is a function of both $\omega_{1}$ and
$\omega_{2}$ parameters, while the same ratio for the spin-boson system (see
Eq. (6)) is only a function of $\omega_{2}$. Similarly as for the spin
coherence, we normalize (20) and (21) to their values at zero temperature

$\displaystyle\langle\sigma_{1}^{x}\rangle_{F,T=0}$
$\displaystyle=2f_{1}f_{2}A\Big{[}\frac{\omega_{c}}{\omega_{1}}+e^{\frac{\omega_{1}}{\omega_{c}}}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}\Big{]},\quad\langle\sigma_{2}^{x}\rangle_{F,T=0}=\frac{t}{\omega_{2}}\langle\sigma_{1}^{x}\rangle_{F,T=0},$
(27)

which simplifies their analysis to

$\displaystyle\frac{\langle\sigma_{1}^{x}\rangle_{F,T}}{\langle\sigma_{1}^{x}\rangle_{F,T=0}}$
$\displaystyle=\tanh\left(\frac{\beta\hbar\omega_{1}}{2}\right)\frac{\frac{\omega_{c}}{\omega_{1}}-\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1}\Big{)}}{\frac{\omega_{c}}{\omega_{1}}+e^{\frac{\omega_{1}}{\omega_{c}}}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}},$
(28)
$\displaystyle\frac{\langle\sigma_{2}^{x}\rangle_{F,T}}{\langle\sigma_{2}^{x}\rangle_{F,T=0}}$
$\displaystyle=-\frac{2\omega_{2}}{\Delta\omega}\frac{\mathcal{F}_{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1},\beta\hbar\omega_{2}\Big{)}}{\frac{\omega_{c}}{\omega_{1}}+e^{\frac{\omega_{1}}{\omega_{c}}}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}}.$
(29)

The low-temperature limit relation
$\langle\sigma_{2}^{x}\rangle_{F,T=0}=\frac{t}{\omega_{2}}\langle\sigma_{1}^{x}\rangle_{F,T=0}$
replicates up to the sign the previous result for the spins.

Figure 3:  Normalized autonomous coherences in the fermion-boson model. (a)
Log-log scale plot of the absolute values of the dimensionless coherences for
the fermionic system
$\langle\sigma_{2}^{x}\rangle_{F,T}/\langle\sigma_{2}^{x}\rangle_{F,T=0}$
(solid lines) and
$\langle\sigma_{1}^{x}\rangle_{F,T}/\langle\sigma_{1}^{x}\rangle_{F,T=0}$
(dashed lines) for different values of $\omega_{1}/\omega_{c}=0.01,1,100$
(blue, yellow and green curves, respectively) in the particular case of
$\omega_{2}/\omega_{1}=3$ as a function of temperature $T$. The inset shows an
advantageous maximum of both internal and output coherences at a finite
temperature for $\omega_{1}>\omega_{c}$. (b) Log-log scale plot of the
absolute value of the dimensionless ratio
$\langle\sigma_{2}^{x}\rangle_{F,T}/\langle\sigma_{1}^{x}\rangle_{F,T=0}$,
multiplied by the normalization factor $\omega_{1}/t$, for the fixed value of
$\omega_{1}/\omega_{c}=0.01$ and different values of
$\omega_{2}/\omega_{1}=0.5,2,4,8,16$ (blue, yellow, green, red and purple
curves, respectively) as a function of temperature $T$. It represents the
dependence of the coherence $\langle\sigma_{2}^{x}\rangle_{F,T}$ on
$\omega_{2}$. The gray dashed line depicts the asymptote of all curves, which
scales as $T^{-2}$. The inset shows the coherence transfer at the low-
temperatures similarly to the spin system.

The normalized internal and output coherences for the fermionic system are
visualized in Fig. 3 (a). In the large temperature limit,
$k_{B}T/\hbar\omega_{1}\gg 1$, the increase of both coherences with decreasing
temperature is linear in the log-log plot, similarly as for the spins. It
means that also for fermionic systems the increase is only by a power law for
large temperatures, and the coherence is negligible. On the other hand,
temperature dependence of the internal coherence
$\langle\sigma_{1}^{x}\rangle_{F,T}/\langle\sigma_{1}^{x}\rangle_{F,T=0}$ is
non-monotonous for a large $\omega_{1}/\omega_{c}$. It has a maximum near
$k_{B}T/\hbar\omega_{1}\approx 1$ with a much larger value than at zero
temperature. Remarkably, this fermionic model obtains larger internal quantum
coherence for a finite temperature away from the zero-temperature limit.
Advantageously, the output coherence follows the same non-monotonous tendency
for $\omega_{1}>\omega_{c}$. The origin of such nontrivial behavior results
from interesting competition between Hartree and Fock terms of the inner
fermionic coherence. The sudden drop of the output coherence due to zero-
crossing is practically irrelevant. The inset shows that the maximum output
coherence shifts to the lower temperature with lowering $\omega_{c}$. However,
relative output coherence increase is similar to the internal one. The
fermionic system brought new positive phenomena into the autonomous quantum
coherence being largest at a finite temperature.

In order to compare the spin and fermion model we consider the high
temperature limit

$-\frac{\omega_{1}}{t}\frac{\langle\sigma_{2}^{x}\rangle_{F,T\rightarrow\infty}}{\langle\sigma_{1}^{x}\rangle_{F,T=0}}=\frac{2\omega_{1}}{\Delta\omega}\frac{\mathcal{F}_{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1},\beta\hbar\omega_{2}\Big{)}}{\frac{\omega_{c}}{\omega_{1}}+e^{\frac{\omega_{1}}{\omega_{c}}}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}}\Big{|}_{\beta\rightarrow
0}\approx\frac{\omega_{c}}{12\omega_{1}}\frac{(\beta\hbar\omega_{1})^{2}}{\frac{\omega_{c}}{\omega_{1}}+e^{\frac{\omega_{1}}{\omega_{c}}}\textrm{Ei}\Big{(}-\frac{\omega_{1}}{\omega_{c}}\Big{)}}\propto(\beta\hbar\omega_{1})^{2},$
(30)

where we have taken into account

$\mathcal{F}_{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\hbar\omega_{1},\beta\hbar\omega_{2}\Big{)}\Big{|}_{\beta\rightarrow
0}\approx\frac{\Delta\omega\omega_{c}}{24\omega_{1}^{2}}(\beta\hbar\omega_{1})^{2}.$
(31)

Note that the high temperature behavior of the ratio
$-\frac{\omega_{1}}{t}\frac{\langle\sigma_{2}^{x}\rangle_{F,T\rightarrow\infty}}{\langle\sigma_{1}^{x}\rangle_{F,T=0}}$
is independent of $\omega_{2}$ and has another power-low dependence $\propto
T^{-2}$ in comparison with the same ratio for the spin system, where it
behaves as $\propto T^{-3}$. Note that the obtained power-law behavior remains
the same for the case of an arbitrary spectral density function
$\mathcal{I}(\xi)$ which satisfies the condition
$\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)/\xi<+\infty$ (see Eq. (217) in
Appendix N). In particular, it is valid for the case of the generalized Ohmic
spectral density function
$\mathcal{I}(\xi,s)\sim\xi^{s}\exp(-\xi/\hbar\omega_{c})$ for $s>0$ (sub-Ohmic
for $0<s<1$ and super-Ohmic for $s>1$, respectively).

The transfer of coherence from internal to output part in the fermionic system
is analogical to the spin system. Only at high temperatures, it scales with
$T^{-2}$ instead of $T^{-3}$ as for the spins. The sudden drops caused by the
zero-crossings are practically irrelevant. In the relevant low-temperature
regime $k_{B}T/\hbar\omega_{1}\leq 1$, the transfer improves for
$\omega_{2}<\omega_{1}$ as in the previous spin analysis. For the fermions,
the optimal regime to obtain and extract autonomous quantum coherence
efficiently at a finite temperature needs $\omega_{1}\gg\omega_{c}$ and
$\omega_{2}<\omega_{1}$. It is complementary with the spins, where the large
band was optimal. Optimal extraction still profits from the small $\omega_{2}$
until the approximation fails.

##  4 Summary

We have hereby confirmed the extractibility of autonomous quantum coherence
for all basic quantum systems coupled to a generic bosonic environment. The
results justified a qualitative prediction from static approximation [40].
They proved that the back action of a weak extraction is fortunately not
critical for the coherence generation and extraction in the low-temperature
limit. Moreover, we obtained the required frequency and temperature
dependences that are not reliably predictable from the static approximation.
We optimized the extraction for spins (oscillators and fields) and fermions
and found different roles of the environment bandwidth to obtain extractable
coherence. Importantly, all those effects are numerically confirmed beyond the
asymptotic limit for a finite strength of the extraction (detailed analysis is
presented in Appendix). Extractability principally opens a way to experimental
verifications of proposals with quantum dots [40] and also to others under
development in the superconducting circuits [45], electro-mechanical systems
[46], and trapped ions [47].

We can extend the systems used to generate and extract quantum coherence. The
fermionic system generating a maximum of coherence at a finite temperature
$k_{B}T\leq\hbar\omega_{1}$ inspires further extensions. We can investigate
three-level and four-level models to optimize the proposed methods. It is a
viable investigation area as such systems appear for both atomic and solid-
state systems. A possible maximum of quantum coherence at the finite
temperature in the region $k_{B}T>\hbar\omega_{1}$ remains an open problem.
However, further engineering of the systems and interaction may overcome such
a limit. Another open problem is an accumulation of coherence in an output
system from many such parallel interactions. Recently, it has been
demonstrated in a collision model that multiple environments can enhance
quantum coherence substantially [39]. Further research can be focused on
coherence distillation and use a possible extension of the protocols to more
spins or fermions giving an even more considerable amount of the extractable
quantum coherence. Thermodynamical analysis of autonomous coherence has been
already initiated too [37]. These directions for immediate exploration
illustrate the novelty of this subject in theory and stimulate current
experiments. It is very timely as quantum coherence from a low-temperature
environment will directly transform quantum sensors and quantum engines to be
more autonomous.

##  5 Acknowledgments

A.S. and R.F. acknowledge the project 20-16577S of the Czech Science
Foundation. R.F. acknowledges also LTAUSA19099 of MEYS CR.

## References

  * [1] Alexander Streltsov, Gerardo Adesso, and Martin B. Plenio. Colloquium: Quantum coherence as a resource. Rev. Mod. Phys. 89, 041003 (2017).
  * [2] Simon Schmitt, Tuvia Gefen, Felix M. Stürner, Thomas Unden, Gerhard Wolff, Christoph M’́uller, Jochen Scheuer, Boris Naydenov, Matthew Markham, Sebastien Pezzagna, Jan Meijer, Ilai Schwarz, Martin Plenio, Alex Retzker, Liam P. McGuinness, Fedor Jelezko. Submillihertz magnetic spectroscopy performed with a nanoscale quantum sensor. Science 356, 6340, 832-837 (2017).
  * [3] Hengyun Zhou, Joonhee Choi, Soonwon Choi, Renate Landig, Alexander M. Douglas, Junichi Isoya, Fedor Jelezko, Shinobu Onoda, Hitoshi Sumiya, Paola Cappellaro, Helena S. Knowles, Hongkun Park, and Mikhail D. Lukin. Quantum Metrology with Strongly Interacting Spin Systems.  Phys. Rev. X 10, 031003 (2020).
  * [4] Andreas Reiserer, Norbert Kalb, Machiel S. Blok, Koen J. M. van Bemmelen, Tim H. Taminiau, Ronald Hanson, Daniel J. Twitchen, and Matthew Markham. Robust Quantum-Network Memory Using Decoherence-Protected Subspaces of Nuclear Spins.  Phys. Rev. X 6, 021040 (2016).
  * [5] David D. Awschalom, Ronald Hanson, Jörg Wrachtrup & Brian B. Zhou. Quantum technologies with optically interfaced solid-state spins. Nature Photonics 12, 516–527 (2018).
  * [6] T. Hensgens, T. Fujita, L. Janssen, Xiao Li, C. J. Van Diepen, C. Reichl, W. Wegscheider, S. Das Sarma & L. M. K. Vandersypen. Quantum simulation of a Fermi–Hubbard model using a semiconductor quantum dot array. Nature 548, 70–73 (2017).
  * [7] Robert Drost, Teemu Ojanen, Ari Harju & Peter Liljeroth. & Liljeroth, P. Topological states in engineered atomic lattices. Nature Physics 13, 668–671 (2017).
  * [8] Marlou R. Slot, Thomas S. Gardenier, Peter H. Jacobse, Guido C. P. van Miert, Sander N. Kempkes, Stephan J. M. Zevenhuizen, Cristiane Morais Smith, Daniel Vanmaekelbergh & Ingmar Swart. Experimental realization and characterization of an electronic Lieb lattice. Nature Physics 13, 672-676 (2017).
  * [9] Gregory D. Scholes, Graham R. Fleming, Lin X. Chen, Alán Aspuru-Guzik, Andreas Buchleitner, David F. Coker, Gregory S. Engel, Rienk van Grondelle, Akihito Ishizaki, David M. Jonas, Jeff S. Lundeen, James K. McCusker, Shaul Mukamel, Jennifer P. Ogilvie, Alexandra Olaya-Castro, Mark A. Ratner, Frank C. Spano, K. Birgitta Whaley & Xiaoyang Zhu. Using coherence to enhance function in chemical and biophysical systems. Nature 543, 647–656 (2017).
  * [10] Elisabet Romero, Vladimir I. Novoderezhkin & Rienk van Grondelle. Quantum design of photosynthesis for bio-inspired solar-energy conversion. Nature 543, 647–656 (2017).
  * [11] James Klatzow, Jonas N. Becker, Patrick M. Ledingham, Christian Weinzetl, Krzysztof T. Kaczmarek, Dylan J. Saunders, Joshua Nunn, Ian A. Walmsley, Raam Uzdin, and Eilon Poem. Experimental demonstration of quantum effects in the operation of microscopic heat engines. Phys. Rev. Lett. 122, 110601 (2019).
  * [12] K. Ono, S. N. Shevchenko, T. Mori, S. Moriyama, and Franco Nori. F. Analog of a Quantum Heat Engine Using a Single-Spin Qubit. Phys. Rev. Lett. 125, 166802 (2020).
  * [13] Camille L. Latune, Ilya Sinayskiy & Francesco Petruccione. Roles of quantum coherences in thermal machines. Eur. Phys. J. Spec. Top. (2021).
  * [14] C. E. Bradley, J. Randall, M. H. Abobeih, R. C. Berrevoets, M. J. Degen, M. A. Bakker, M. Markham, D. J. Twitchen, and T. H. Taminiau. A Ten-Qubit Solid-State Spin Register with Quantum Memory up to One Minute. Phys. Rev. X 9, 031045 (2019).
  * [15] C. J. Stephen, B. L. Green, Y. N. D. Lekhai, L. Weng, P. Hill, S. Johnson, A. C. Frangeskou, P. L. Diggle, Y.-C. Chen, M. J. Strain, E. Gu, M. E. Newton, J. M. Smith, P. S. Salter, and G. W. Morley. Deep Three-Dimensional Solid-State Qubit Arrays with Long-Lived Spin Coherence. Phys. Rev. Applied 12, 064005 (2019).
  * [16] Flavio Del Santo and Borivoje Dakić. Coherence Equality and Communication in a Quantum Superposition. Phys. Rev. Lett. 124, 190501 (2020).
  * [17] Flavio Del Santo and Borivoje Dakić. Two-Way Communication with a Single Quantum Particle. Phys. Rev. Lett. 120, 060503 (2018).
  * [18] Kok Chuan Tan, Tyler Volkoff, Hyukjoon Kwon, and Hyunseok Jeong. Quantifying the Coherence between Coherent States. Phys. Rev. Lett. 119, 190405 (2017).
  * [19] Kaonan Micadei, John P. S. Peterson, Alexandre M. Souza, Roberto S. Sarthour, Ivan S. Oliveira, Gabriel T. Landi, Roberto M. Serra, and Eric Lutz. Experimental Validation of Fully Quantum Fluctuation Theorems Using Dynamic Bayesian Networks Phys.Rev.Lett. 127, 180603 (2021).
  * [20] Jose Joaquin Alonso, Eric Lutz, and Alessandro Romito. Thermodynamics of Weakly Measured Quantum Systems. Phys. Rev. Lett. 116, 080403 (2016).
  * [21] J. F. Haase, A. Smirne, J. Kołodyński, R. Demkowicz-Dobrzański, and S. F. Huelga, Fundamental limits to frequency estimation: a comprehensive microscopic perspective. New J. Phys. 20, 053009 (2018).
  * [22] Jan Czajkowski, Krzysztof Pawłowski and Rafał Demkowicz-Dobrzański. Many-body effects in quantum metrology. New J. Phys. 21 053031 (2019).
  * [23] Leonardo Novo, Masoud Mohseni & Yasser Omar. Disorder-assisted quantum transport in suboptimal decoherence regimes. Scientific Reports 6, 18142 (2016).
  * [24] Kaonan Micadei, Gabriel T. Landi, and Eric Lutz. Quantum Fluctuation Theorems beyond Two-Point Measurements. Phys. Rev. Lett. 124, 090602 (2020).
  * [25] María García Díaz, Giacomo Guarnieri, Mauro Paternostro. Quantum work statistics with initial coherence. Entropy 22(11), 1223 (2020).
  * [26] Rafał Demkowicz-Dobrzański, Jan Czajkowski, and Pavel Sekatski. Adaptive Quantum Metrology under General Markovian Noise. Phys. Rev. X 7, 041009 (2017).
  * [27] Stella Seah, Stefan Nimmrichter, and Valerio Scarani. Maxwell’s Lesser Demon: A Quantum Engine Driven by Pointer Measurements. Phys. Rev. Lett. 124, 100603 (2020).
  * [28] Harry J. D. Miller, Giacomo Guarnieri, Mark T. Mitchison, and John Goold. Quantum fluctuations hinder finite-time information erasure near the Landauer limit. Phys. Rev. Lett. 125, 160602 (2020).
  * [29] G. Francica, F. C. Binder, G. Guarnieri, M. T. Mitchison, J. Goold, and F. Plastina. Quantum Coherence and Ergotropy. Phys. Rev. Lett. 125, 180603 (2020).
  * [30] Stella Seah, Stefan Nimmrichter, Daniel Grimmer, Jader P. Santos, Valerio Scarani, and Gabriel T. Landi. Collisional Quantum Thermometry. Phys. Rev. Lett. 123, 180602 (2019).
  * [31] C. L. Latune, I. Sinayskiy & F. Petruccione. Quantum coherence, many-body correlations, and non-thermal effects for autonomous thermal machines. Scientific Reports 9, 3191 (2019).
  * [32] Devashish Tupkary, Abhishek Dhar, Manas Kulkarni, Archak Purkayastha. Fundamental limitations in Lindblad descriptions of systems weakly coupled to baths. Phys. Rev. A 105, 032208 (2022)
  * [33] Wojciech Hubert Zurek. Decoherence, einselection, and the quantum origins of the classical. Rev. Mod. Phys. 75, 715 (2003).
  * [34] Mohseni, M., Omar, Y., Engel, G. S., Plenio, M. Quantum Effects in Biology (Cambridge Univ Press, Cambridge, MA, 2014).
  * [35] Hong-Guang Duan, Valentyn I. Prokhorenko, Richard J. Cogdell, Khuram Ashraf, Amy L. Stevens, Michael Thorwart, and R. J. Dwayne Miller. Nature does not rely on long-lived electronic quantum coherence for photosynthetic energy transfer. PNAS 114 (32) 8493-8498 (2017).
  * [36] Giacomo Guarnieri, Michal Kolář, and Radim Filip. Steady-State Coherences by Composite System-Bath Interactions. Phys. Rev. Lett. 121, 070401 (2018).
  * [37] Giacomo Guarnieri, Daniele Morrone, Barış Çakmak, Francesco Plastina, Steve Campbelle. Non-equilibrium steady-states of memoryless quantum collision models. Phys. Lett. A 384, 24, 126576 (2020).
  * [38] Mike Reppert, Deborah Reppert, Leonardo A. Pachon, and Paul Brumer. Equilibrium stationary coherence in the multilevel spin-boson model. Phys. Rev. A 102, 012211 (2020).
  * [39] Román-Ancheyta, R., Kolář, M., Guarnieri, G., Filip, R. Enhanced steady-state coherences via repeated system-bath interactions. Phys. Rev. A 104, 062209 (2021).
  * [40] Archak Purkayastha, Giacomo Guarnieri, Mark T. Mitchison, Radim Filip & John Goold. Tunable phonon-induced steady-state coherence in a double-quantum-dot charge qubit. npj Quantum Information 6, 27 (2020).
  * [41] H.-P. Breuer, F.Petruccione. The Theory of Open Quantum Systems (Oxford University Press, 2007)
  * [42] U. Weiss. Quantum Dissipative Systems (World Scientific Publishing, 2012)
  * [43] H. Bateman, A. Erdélyi. Higher Transcendental Functions. Vol.2. Bessel Functions, Parabolic Cylinder Functions, Orthogonal Polynomials (McGraw-Hill, 1953).
  * [44] Elliott Lieb,Theodore Schultz, Daniel Mattis. Two soluble models of an antiferromagnetic chain. Annals of Physics 16, 407 (1961).
  * [45] S. Touzard, A. Kou, N. E. Frattini, V. V. Sivak, S. Puri, A. Grimm, L. Frunzio, S. Shankar, and M. H. Devoret. Gated Conditional Displacement Readout of Superconducting Qubits. Phys. Rev. Lett. 122, 080502 (2019).
  * [46] X. Ma, J. J. Viennot, S. Kotler, J. D. Teufel & K. W. Lehnert. Non-classical energy squeezing of a macroscopic mechanical oscillator.  Nature Physics 17, 322-326 (2021).
  * [47] M.-L. Cai, Z.-D. Liu, W.-D. Zhao, Y.-K. Wu, Q.-X. Mei, Y. Jiang, L. He, X. Zhang, Z.-C. Zhou & L.-M. Duan. Observation of a quantum phase transition in the quantum Rabi model with a single trapped ion. Nature Communications 12, 1126 (2021).

##  Appendix A Spin-boson Hamiltonian

We consider the Hamiltonian of two interacting spins coupled to the boson
system

$H=\frac{\omega_{1}}{2}\sigma_{1}^{z}+\frac{\omega_{2}}{2}\sigma_{2}^{z}+t(\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{2}^{+}\sigma_{1}^{-})+\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}+(f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x})\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}).$
(32)

Here $\sigma_{1}^{j}$ and $\sigma_{2}^{j}$ (with $j=x,y,z$) are the Pauli
matrices corresponding to the $1$-st and $2$-nd spins with the energies
$\omega_{1}$ and $\omega_{2}$ respectively.The real valued parameter $t>0$
defines the coupling between the spins. Note that in this section and in
further we consider $\omega_{1}$, $\omega_{2}$ and $t$ in energy units for
brevity. In order to restore the original frequency units, used in the main
text, one needs to make the following substitutions
$\omega_{1}\rightarrow\hbar\omega_{1}$, $\omega_{2}\rightarrow\hbar\omega_{2}$
and $t\rightarrow\hbar t$. The raising $\sigma_{n}^{+}$ and lowering
$\sigma_{n}^{-}$ operators of $n(=1,2)$-th spin are
$\sigma_{n}^{\pm}=(\sigma_{n}^{x}\pm i\sigma_{n}^{y})/2$. The boson system is
described by creation $b_{k}^{\dagger}$ and annihilation $b_{k}$ operators and
spectrum $\Omega_{k}>0$. Subscript $k$ parameterizes all the boson degrees of
freedom. The interaction between spins and bosons are defined by real valued
coupling constants $f_{1}$ and $f_{2}$. It is convenient to introduce the
following notations for spin Hamiltonian

$H_{12}=\frac{\omega_{1}}{2}\sigma_{1}^{z}+\frac{\omega_{2}}{2}\sigma_{2}^{z}+t(\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{2}^{+}\sigma_{1}^{-}),$
(33)

boson Hamiltonian

$H_{B}=\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k},$ (34)

and the spin-boson coupling term

$H_{B1}=(f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x})\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}).$
(35)

As a first step we write the spin Hamiltonian in the basis of triplet and
singlet states of two spins

$\displaystyle|1,1\rangle=$
$\displaystyle|\uparrow_{1}\rangle|\uparrow_{2}\rangle,\quad|1,-1\rangle=|\downarrow_{1}\rangle|\downarrow_{2}\rangle,\quad|1,0\rangle=\frac{1}{\sqrt{2}}(|\uparrow_{1}\rangle|\downarrow_{2}\rangle+|\downarrow_{1}\rangle|\uparrow_{2}\rangle),$
(36) $\displaystyle|0,0\rangle=$
$\displaystyle\frac{1}{\sqrt{2}}(|\uparrow_{1}\rangle|\downarrow_{2}\rangle-|\downarrow_{1}\rangle|\uparrow_{2}\rangle),$

where $|\uparrow_{n}\rangle$ and $|\downarrow_{n}\rangle$ are the eigenstates
of the $\sigma_{n}^{z}$ operators with eigenvalues $\pm 1$ respectively. In
this basis the Hamiltonian has the following matrix form

$H_{12}=\frac{1}{2}\left[\begin{array}[]{cccc}\omega_{1}+\omega_{2}&0&0&0\\\
0&-\omega_{1}-\omega_{2}&0&0\\\ 0&0&2t&\omega_{1}-\omega_{2}\\\
0&0&\omega_{1}-\omega_{2}&-2t\end{array}\right].$ (37)

The non-zero difference $\delta\omega=\omega_{1}-\omega_{2}\neq 0$ mixes the
spin states with zero angular momentum. We consider their linear combination
in order to diagonalize the Hamiltonian matrix

$|+,0\rangle=\cos(\theta/2)|1,0\rangle+\sin(\theta/2)|0,0\rangle,\quad|-,0\rangle=-\sin(\theta/2)|1,0\rangle+\cos(\theta/2)|0,0\rangle.$
(38)

The matrix becomes diagonal for the case
$\cos\theta=2t/\sqrt{4t^{2}+\delta\omega^{2}}$,
$\sin\theta=\delta\omega/\sqrt{4t^{2}+\delta\omega^{2}}$. The eigenvalues
become $E_{1}=(\omega_{1}+\omega_{2})/2=\omega$,
$E_{2}=-(\omega_{1}+\omega_{2})/2=-\omega$,
$E_{3}=\sqrt{t^{2}+\delta\omega^{2}/4}=\epsilon$ and
$E_{4}=-\sqrt{t^{2}+\delta\omega^{2}/4}=-\epsilon$ for the states
$|1,1\rangle$, $|1,-1\rangle$, $|+,0\rangle$ and $|-,0\rangle$ respectively.
The spin operators $\sigma_{1}^{x}$ and $\sigma_{1}^{z}$ have the following
form in the basis $|1,1\rangle,|1,-1\rangle,|+,0\rangle,|-,0\rangle$

$\sigma_{1}^{x}=\left[\begin{array}[]{cccc}0&0&\sin\phi&-\cos\phi\\\
0&0&\cos\phi&\sin\phi\\\ \sin\phi&\cos\phi&0&0\\\
-\cos\phi&\sin\phi&0&0\end{array}\right],\quad\sigma_{1}^{z}=\left[\begin{array}[]{cccc}1&0&0&0\\\
0&-1&0&0\\\ 0&0&\sin\theta&\cos\theta\\\
0&0&\cos\theta&-\sin\theta\end{array}\right],$ (39)

where $\phi=\pi/4-\theta/2$.

##  Appendix B Perturbation analysis

We treat the Hamiltonian $H_{0}=H_{12}+H_{B}$ exactly and $H_{B1}$ as a
perturbation [40]. We present the density matrix operator in the form
($\beta=1/T$, $T$ is the temperature in energetic units)

$e^{-\beta H}=e^{-\beta(H_{12}+H_{B})}X(\beta).$ (40)

Here the operator $X(\beta)$ is presented as a $\mathrm{T}$-exponent

$X(\beta)=\mathrm{T}\Big{(}-\int_{0}^{\beta}d\tau\widetilde{H}_{B1}(\tau)\Big{)}\approx
1-\int_{0}^{\beta}d\tau\widetilde{H}_{B1}(\tau)+\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\widetilde{H}_{B1}(\tau)\widetilde{H}_{B1}(\tau^{\prime}),$
(41)

where

$\widetilde{H}_{B1}(\tau)=e^{\tau(H_{12}+H_{B})}H_{B1}e^{-\tau(H_{12}+H_{B})}.$
(42)

Then the perturbation term
$H_{B1}=SB=[f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x}][\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})]$
takes the form

$\displaystyle\widetilde{H}_{B1}(\tau)=$
$\displaystyle\widetilde{S}(\tau)\widetilde{B}(\tau)=$ (43) $\displaystyle=$
$\displaystyle[f_{1}e^{\tau H_{12}}\sigma_{1}^{z}e^{-\tau H_{12}}+f_{2}e^{\tau
H_{12}}\sigma_{1}^{x}e^{-\tau H_{12}}\Big{]}\Big{[}\sum_{k}\lambda_{k}(e^{\tau
H_{B}}b_{k}^{\dagger}e^{-\tau H_{B}}+e^{\tau H_{B}}b_{k}e^{-\tau
H_{B}})\Big{]}.$

Taking into account the structure of the $H_{12}$ and $H_{B}$ Hamiltonians and
using the notation $\Lambda_{1}(\tau)=e^{\tau\omega}$,
$\Lambda_{2}(\tau)=e^{\tau\epsilon}$ we obtain

$e^{\tau H_{B}}b_{k}^{\dagger}e^{-\tau
H_{B}}=e^{\tau\Omega_{k}}b_{k}^{\dagger},\quad e^{\tau H_{B}}b_{k}e^{-\tau
H_{B}}=e^{-\tau\Omega_{k}}b_{k},$ (44)
$\displaystyle\sigma_{1}^{z}(\tau)=e^{\tau H_{12}}\sigma_{1}^{z}e^{-\tau
H_{12}}=\left[\begin{array}[]{cccc}1&0&0&0\\\ 0&-1&0&0\\\
0&0&\sin\theta&\Lambda_{2}^{2}(\tau)\cos\theta\\\
0&0&\Lambda_{2}^{-2}(\tau)\cos\theta&-\sin\theta\end{array}\right],$ (49)
$\displaystyle\sigma_{1}^{x}(\tau)=e^{\tau H_{12}}\sigma_{1}^{x}e^{-\tau
H_{12}}=$ (50)
$\displaystyle=\left[\begin{array}[]{cccc}0&0&\Lambda_{1}(\tau)\Lambda_{2}^{-1}(\tau)\sin\phi&-\Lambda_{1}(\tau)\Lambda_{2}(\tau)\cos\phi\\\
0&0&\Lambda_{1}^{-1}(\tau)\Lambda_{2}^{-1}(\tau)\cos\phi&\Lambda_{1}^{-1}(\tau)\Lambda_{2}(\tau)\sin\phi\\\
\Lambda_{1}^{-1}(\tau)\Lambda_{2}(\tau)\sin\phi&\Lambda_{1}(\tau)\Lambda_{2}(\tau)\cos\phi&0&0\\\
-\Lambda_{1}^{-1}(\tau)\Lambda_{2}^{-1}(\tau)\cos\phi&\Lambda_{1}(\tau)\Lambda_{2}^{-1}(\tau)\sin\phi&0&0\end{array}\right].$
(55)

For calculation of the averages of any operator we use the symmetrization of
the density operator

$e^{-\beta
H}=\frac{1}{2}\Big{[}e^{-\beta(H_{12}+H_{B})}X(\beta)+X^{\dagger}(\beta)e^{-\beta(H_{12}+H_{B})}\Big{]}.$
(56)

Then the partition function $Z$ can be expressed as

$Z=\mathrm{Tr}[e^{-\beta H}]=\frac{Z_{0}}{2}\langle
X(\beta)+X^{\dagger}(\beta)\rangle_{0}.$ (57)

Here
$\langle\star\rangle_{0}=Z_{0}^{-1}\mathrm{Tr}[e^{-\beta(H_{12}+H_{B})}\star]$
is the averaging by the non-interaction density function. The partition
function up to second order correction takes the following form

$Z=Z_{0}\Big{\\{}1+\frac{1}{2}\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}[\langle\widetilde{H}_{B1}(\tau)\widetilde{H}_{B1}(\tau^{\prime})\rangle_{0}+\langle\widetilde{H}_{B1}^{\dagger}(\tau^{\prime})\widetilde{H}_{B1}^{\dagger}(\tau)\rangle_{0}]\Big{\\}}.$
(58)

Taking the average over the boson degrees of freedom we get

$Z=Z_{0}\Big{\\{}1+\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\tau-\tau^{\prime})\mathrm{Re}\langle\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})\rangle_{0}\Big{\\}},$
(59)

where we introduced the bosonic correlation function

$\phi(\tau-\tau^{\prime})=\langle\widetilde{B}(\tau)\widetilde{B}(\tau^{\prime})\rangle_{0}=\langle\widetilde{B}^{\dagger}(\tau^{\prime})\widetilde{B}^{\dagger}(\tau)\rangle_{0}=\sum_{k}\lambda_{k}^{2}[e^{\Omega_{k}(\tau-\tau^{\prime})}\langle
b_{k}^{\dagger}b_{k}\rangle_{0}+e^{-\Omega_{k}(\tau-\tau^{\prime})}\langle
b_{k}b_{k}^{\dagger}\rangle_{0}].$ (60)

The spin dependent part takes the form

$\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})=f_{1}^{2}\sigma_{1}^{z}(\tau)\sigma_{1}^{z}(\tau^{\prime})+f_{2}^{2}\sigma_{1}^{x}(\tau)\sigma_{1}^{x}(\tau^{\prime})+f_{1}f_{2}\sigma_{1}^{z}(\tau)\sigma_{1}^{x}(\tau^{\prime})+f_{1}f_{2}\sigma_{1}^{x}(\tau)\sigma_{1}^{z}(\tau^{\prime}).$
(61)

We are interested in the spin sector only, therefore we integrate out the
boson degrees of freedom in total density operator and keep only the spin
dependent part $\varrho_{S}=\mathrm{Tr}_{B}[e^{-\beta H}]/Z$

$\displaystyle\varrho_{S}=$ $\displaystyle\frac{e^{-\beta
H_{12}}}{Z_{S}}\Big{[}1-\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\tau-\tau^{\prime})\mathrm{Re}\langle\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})\rangle_{0}\Big{]}+$
$\displaystyle+$
$\displaystyle\frac{1}{2}\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\tau-\tau^{\prime})\Big{[}\frac{e^{-\beta
H_{12}}}{Z_{S}}\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})+\textrm{h.c.}\Big{]}.$
(62)

Here we introduced $Z_{S}=\mathrm{Tr}[e^{-\beta
H_{12}}]=2\cosh(\beta\omega)+2\cosh(\beta\epsilon)$.

##  Appendix C Averaging procedure for the second spin operators

We are interested in the averages
$\langle\sigma_{2}^{j}\rangle=\mathrm{Tr}[\varrho_{S}\sigma_{2}^{j}]$, where
$j=x,y,z$. The matrices of corresponding spin operators in the basis
$\\{|1,1\rangle,|1,-1\rangle,|+,0\rangle,|-,0\rangle\\}$ are

$\displaystyle\sigma_{2}^{z}=$
$\displaystyle\left[\begin{array}[]{cccc}1&0&0&0\\\ 0&-1&0&0\\\
0&0&-\sin\theta&-\cos\theta\\\
0&0&-\cos\theta&\sin\theta\end{array}\right],\quad\sigma_{2}^{x}=\left[\begin{array}[]{cccc}0&0&\cos\phi&\sin\phi\\\
0&0&\sin\phi&-\cos\phi\\\ \cos\phi&\sin\phi&0&0\\\
\sin\phi&-\cos\phi&0&0\end{array}\right],$ (71) $\displaystyle\sigma_{2}^{y}=$
$\displaystyle\left[\begin{array}[]{cccc}0&0&-i\cos\phi&-i\sin\phi\\\
0&0&i\sin\phi&-i\cos\phi\\\ i\cos\phi&-i\sin\phi&0&0\\\
i\sin\phi&i\cos\phi&0&0\end{array}\right].$ (76)

where $\phi=\pi/4-\theta/2$. In this basis the procedure of taking trace
becomes matrix trace procedure. First we calculate the leading terms

$Z_{S}\langle\sigma_{2}^{z}\rangle_{0}=\mathrm{Tr}[e^{-\beta
H_{12}}\sigma_{2}^{z}]=2\sin\theta\sinh(\beta\epsilon)-2\sinh[\beta\omega],\,Z_{S}\langle\sigma_{2}^{x(y)}\rangle_{0}=\mathrm{Tr}[e^{-\beta
H_{12}}\sigma_{2}^{x(y)}]=0.$ (77)

Then we introduce the traces of the product of the spin operators

$\displaystyle F(\tau,\tau^{\prime})$
$\displaystyle=Z_{S}\mathrm{Re}\langle\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})\rangle_{0}=$
$\displaystyle
2f_{1}^{2}\Big{[}\cosh(\beta\omega)+\sin^{2}\theta\cosh(\beta\epsilon)+\cos^{2}\theta\cosh((\beta-2\tau+2\tau^{\prime})\epsilon)\Big{]}+$
$\displaystyle
2f_{2}^{2}\Big{[}\cos^{2}\phi\\{\cosh(\beta\epsilon-(\epsilon+\omega)(\tau-\tau^{\prime}))+\cosh(\beta\omega-(\epsilon+\omega)(\tau-\tau^{\prime}))\\}\Big{]}+$
$\displaystyle
2f_{2}^{2}\Big{[}\sin^{2}\phi\\{\cosh(\beta\epsilon-(\epsilon-\omega)(\tau-\tau^{\prime}))+\cosh(\beta\omega+(\epsilon-\omega)(\tau-\tau^{\prime}))\\}\Big{]},$
(78) $\displaystyle F_{1}(\tau,\tau^{\prime})=\mathrm{Tr}$
$\displaystyle[e^{-\beta
H_{12}}\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})\sigma_{2}^{z}+\widetilde{S}^{\dagger}(\tau^{\prime})\widetilde{S}^{\dagger}(\tau)e^{-\beta
H_{12}}\sigma_{2}^{z}]=$ $\displaystyle 4f_{1}^{2}$
$\displaystyle\Big{[}\sin^{3}\theta\sinh(\beta\epsilon)-\sinh(\beta\omega)+\sin\theta\cos^{2}\theta\\{\sinh(\epsilon(\beta+2\tau^{\prime}-2\tau))+$
$\displaystyle+\sinh(\epsilon(\beta-2\tau^{\prime}))-\sinh(\epsilon(\beta-2\tau))\\}\Big{]}+$
$\displaystyle 4f_{2}^{2}$
$\displaystyle\Big{[}\sin\theta\sin^{2}\phi\sinh(\beta\epsilon+(\tau^{\prime}-\tau)(\epsilon-\omega))+$
$\displaystyle\sin\theta\cos^{2}\phi\sinh(\beta\epsilon+(\tau^{\prime}-\tau)(\omega+\epsilon))+$
$\displaystyle\cos\theta\sin\phi\cos\phi\\{\sinh(\beta\epsilon+\tau^{\prime}(\omega-\epsilon)-\tau(\omega+\epsilon))-$
$\displaystyle\sinh(\beta\epsilon-\tau^{\prime}(\omega+\epsilon)+\tau(\omega-\epsilon))\\}-$
$\displaystyle\sin^{2}\phi\sinh(\omega\beta+(\tau^{\prime}-\tau)(\omega-\epsilon))-\cos^{2}\phi\sinh(\beta\omega+(\tau^{\prime}-\tau)(\omega+\epsilon))\Big{]},$
(79) $\displaystyle F_{2}(\tau,\tau^{\prime})$
$\displaystyle=\mathrm{Tr}[e^{-\beta
H_{12}}\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})\sigma_{2}^{x}+\widetilde{S}^{\dagger}(\tau^{\prime})\widetilde{S}^{\dagger}(\tau)e^{-\beta
H_{12}}\sigma_{2}^{x}]=$ $\displaystyle
4f_{1}f_{2}\Big{[}-\cos\theta\cos^{2}\phi\\{\cosh(\beta\omega+2\epsilon\tau^{\prime}-\tau(\omega+\epsilon))+\cosh(\epsilon(\beta-2\tau)+\tau^{\prime}(\omega+\epsilon))\\}+$
$\displaystyle+\cos\theta\sin^{2}\phi\\{\cosh(\beta\omega-2\epsilon\tau^{\prime}-\tau(\omega-\epsilon))+\cosh(\epsilon(\beta-2\tau)+\tau^{\prime}(\epsilon-\omega))\\}+$
$\displaystyle+\sin(2\phi)\sin\theta\\{(\cosh(\omega\tau^{\prime})\cosh(\epsilon(\beta-\tau^{\prime}))+\cosh(\tau\epsilon)\cosh(\omega(\beta-\tau))\\}+$
$\displaystyle+\sin(2\phi)\\{\sinh(\epsilon\tau^{\prime})\sinh(\omega(\beta-\tau^{\prime}))+\sinh(\tau\omega)\sinh(\epsilon(\beta-\tau)))\\}\Big{]},$
(80) $\mathrm{Tr}[e^{-\beta
H_{12}}\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})\sigma_{2}^{y}+\widetilde{S}^{\dagger}(\tau^{\prime})\widetilde{S}^{\dagger}(\tau)e^{-\beta
H_{12}}\sigma_{2}^{y}]=0.$ (81)

##  Appendix D Integral representation for the spin average values

The average values $\langle\sigma_{2}^{x}\rangle$ and
$\langle\sigma_{2}^{z}\rangle$ are defined by the convolution integrals of
$F(\tau,\tau^{\prime})$, $F_{1}(\tau,\tau^{\prime})$ and
$F_{2}(\tau,\tau^{\prime})$ with $\phi(\tau-\tau^{\prime})$ in
$\tau,\tau^{\prime}$ domain. Such the integration can be done in the following
way. We introduce the boson spectral function

$\mathcal{I}(\xi)=\sum_{k}\lambda_{k}^{2}\delta(\xi-\Omega_{k})$ (82)

and present $\phi(\tau-\tau^{\prime})$ in the form

$\phi(\tau-\tau^{\prime})=\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\left\\{\frac{e^{\xi(\tau-\tau^{\prime})}}{e^{\beta\xi}-1}+\frac{e^{-\xi(\tau-\tau^{\prime})}e^{\beta\xi}}{e^{\beta\xi}-1}\right\\},$
(83)

taking into account that $\Omega_{k}>0$. In this notation we can express the
non-zero corrections of average values of spin operators in the form

$\Delta\langle\sigma_{2}^{z}\rangle=\langle\sigma_{2}^{z}\rangle-\langle\sigma_{2}^{z}\rangle_{0}=\frac{1}{2Z_{S}}\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\tau-\tau^{\prime})\Big{[}F_{1}(\tau,\tau^{\prime})-2\langle\sigma_{2}^{z}\rangle_{0}F(\tau,\tau^{\prime})\Big{]},$
(84)
$\langle\sigma_{2}^{x}\rangle=\langle\sigma_{2}^{x}\rangle-\langle\sigma_{2}^{x}\rangle_{0}=\frac{1}{2Z_{S}}\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\tau-\tau^{\prime})F_{2}(\tau,\tau^{\prime}),$
(85)

since $\langle\sigma_{2}^{x}\rangle_{0}=0$. Taking into account the structure
of $\phi(\tau-\tau^{\prime})$ function and introducing the transformation

$\mathcal{F}(\xi)=\frac{1}{e^{\beta\xi}-1}\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}e^{\xi(\tau-\tau^{\prime})}\mathcal{F}(\tau,\tau^{\prime}),$
(86)

we can write

$\Delta\langle\sigma_{2}^{z}\rangle=\langle\sigma_{2}^{z}\rangle-\langle\sigma_{2}^{z}\rangle_{0}=\frac{1}{2Z_{S}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{\\{}F_{1}(\xi)-F_{1}(-\xi)-2\langle\sigma_{2}^{z}\rangle_{0}[F(\xi)-F(-\xi)]\Big{\\}},$
(87)
$\langle\sigma_{2}^{x}\rangle=\frac{1}{2Z_{S}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{\\{}F_{2}(\xi)-F_{2}(-\xi)\Big{\\}}.$
(88)

##  Appendix E Detuned spins: the second spin coherence

Let us consider the case when $|\omega_{1}-\omega_{2}|\gg t$. It corresponds
to the situation when the energy levels of considered spins are weakly
modified by interaction terms. For this case the calculation of the the
averages of spin operators can be simplified. Indeed, we consider the leading
terms of $\cos(\theta)$ and $\sin(\theta)$ by small parametr
$t/|\omega_{1}-\omega_{2}|$

$\displaystyle\sin\theta=$
$\displaystyle\frac{(\omega_{1}-\omega_{2})}{\sqrt{4t^{2}+(\omega_{1}-\omega_{2})^{2}}}\approx\text{sgn}(\omega_{1}-\omega_{2})\Big{[}1-\frac{2t^{2}}{(\omega_{1}-\omega_{2})^{2}}\Big{]},$
$\displaystyle\cos\theta=$
$\displaystyle\frac{2t}{\sqrt{4t^{2}+(\omega_{1}-\omega_{2})^{2}}}\approx\frac{2t}{|\omega_{1}-\omega_{2}|},$
(89)

substitute them in the expressions for $F(\tau,\tau^{\prime})$,
$F_{1}(\tau,\tau^{\prime})$ and $F_{2}(\tau,\tau^{\prime})$ and calculate the
corresponding integrals over $\tau,\tau^{\prime}$ and $\xi$. We introduce the
notations $\eta=\text{sgn}(\omega_{1}-\omega_{2})$ and
$\delta=|\omega_{1}-\omega_{2}|$ in further calculations for brevity. We focus
on the evaluation of the coherence of the second spin operator
$\langle\sigma_{2}^{x}\rangle$. In order to do it we substitute

$\displaystyle\sin^{2}\theta\approx$ $\displaystyle
1-\frac{4t^{2}}{\delta^{2}},\quad\cos^{2}\theta\approx\frac{4t^{2}}{\delta^{2}},$
$\displaystyle\cos^{2}\phi=$
$\displaystyle\frac{1+\sin\theta}{2}\approx\frac{1+\eta}{2}-\frac{\eta
t^{2}}{\delta^{2}},\quad\sin^{2}\phi=\frac{1-\sin\theta}{2}\approx\frac{1-\eta}{2}+\frac{\eta
t^{2}}{\delta^{2}},$ (90)

in Eq. (C) and keep the leading term in the series expansion of this
expression by the small parameter $t/\delta$. We have

$\displaystyle
F_{2}(\tau,\tau^{\prime})=\frac{8f_{1}f_{2}t}{\omega_{1}-\omega_{2}}$
$\displaystyle\Big{[}\sinh\left(\frac{1}{2}\left(\omega_{1}-\omega_{2}\right)\tau^{\prime}\right)\sinh\left(\frac{1}{2}\left(\omega_{1}+\omega_{2}\right)\left(\beta-\tau^{\prime}\right)\right)+$
$\displaystyle+\cosh\left(\frac{1}{2}\left(\omega_{1}+\omega_{2}\right)\tau^{\prime}\right)\cosh\left(\frac{1}{2}\left(\omega_{1}-\omega_{2}\right)\left(\beta-\tau^{\prime}\right)\right)-$
$\displaystyle-\cosh\left(\frac{1}{2}\left(\omega_{1}-\omega_{2}\right)(\beta-2\tau)+\omega_{1}\tau^{\prime}\right)-$
$\displaystyle-\cosh\left(-\frac{1}{2}\beta\left(\omega_{1}+\omega_{2}\right)+\left(\omega_{2}-\omega_{1}\right)\tau^{\prime}+\tau\omega_{1}\right)+$
$\displaystyle+\sinh\left(\frac{1}{2}\tau\left(\omega_{1}+\omega_{2}\right)\right)\sinh\left(\frac{1}{2}\left(\omega_{1}-\omega_{2}\right)(\beta-\tau)\right)+$
$\displaystyle+\cosh\left(\frac{1}{2}\tau\left(\omega_{1}-\omega_{2}\right)\right)\cosh\left(\frac{1}{2}\left(\omega_{1}+\omega_{2}\right)(\beta-\tau)\right)\Big{]}.$
(91)

After integration over $\tau,\tau^{\prime}$ we get

$F_{2}(\xi)-F_{2}(-\xi)=\frac{32f_{1}f_{2}t\sinh\left(\frac{\beta\omega_{2}}{2}\right)\sinh\left(\frac{\beta\omega_{1}}{2}\right)\left(\xi\coth\left(\frac{\beta\xi}{2}\right)-\omega_{1}\coth\left(\frac{\beta\omega_{1}}{2}\right)\right)}{\xi\omega_{2}\left(\xi^{2}-\omega_{1}^{2}\right)}$
(92)

Therefore the coherence $\langle\sigma_{2}^{x}\rangle$ takes the form

$\displaystyle\langle\sigma_{2}^{x}\rangle=$
$\displaystyle\frac{1}{2Z_{S}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{\\{}F_{2}(\xi)-F_{2}(-\xi)\Big{\\}}=$
$\displaystyle=$
$\displaystyle\frac{4f_{1}f_{2}t}{\omega_{2}}\tanh\Big{(}\frac{\beta\omega_{2}}{2}\Big{)}\tanh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\Big{(}\frac{\beta\xi}{2}\Big{)}-\omega_{1}\coth\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}}{\xi(\xi^{2}-\omega_{1}^{2})},$
(93)

where we use $Z_{S}$ in the leading order by $t/\delta$:
$Z_{S}=4\cosh(\beta\omega_{1}/2)\cosh(\beta\omega_{2}/2)$.

For the particular case of spectral density function
$\mathcal{I}(\xi)=A\xi\exp(-\xi/\omega_{c})$ this expression reads

$\langle\sigma_{2}^{x}\rangle=\frac{4Af_{1}f_{2}t}{\omega_{2}}\tanh\left(\frac{\beta\omega_{2}}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1}),$
(94)

where

$\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})=\int_{0}^{\infty}dx\,\exp(-x\omega_{1}/\omega_{c})\frac{x\coth\left(\frac{\beta\omega_{1}x}{2}\right)-\coth\left(\frac{\beta\omega_{1}}{2}\right)}{x^{2}-1},$
(95)

is a dimensionless function of parameters $\omega_{1}/\omega_{c}$ and
$\beta\omega_{1}$. Here $\omega_{c}$ is the cut-off energy of the bosonic
spectral function. The function
$\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})$ for different values of
$\omega_{1}/\omega_{c}=0.1,1,10$ is presented on the Fig. 4

Figure 4:  Dimensionless function
$\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})$ for
$\omega_{1}/\omega_{c}=0.1$ (blue curve), $\omega_{1}/\omega_{c}=1$ (yellow
curve), and $\omega_{1}/\omega_{c}=10$ (green curve) as a function of
dimensionless parameter $\beta\omega_{1}\in[0,20]$.

Let us estimate the value of the integral
$\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})$ in the low-temperature
limit $\beta\omega_{1}\gg 1$. In order to do it we present the integral in the
form

$\displaystyle\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})=$
$\displaystyle\int_{0}^{\infty}dx\frac{\exp(-x\omega_{1}/\omega_{c})}{x+1}-\frac{2}{e^{\beta\omega_{1}}-1}\text{v.p.}\int_{0}^{\infty}dx\frac{\exp(-x\omega_{1}/\omega_{c})}{x^{2}-1}+$
$\displaystyle+$
$\displaystyle\text{v.p.}\int_{0}^{\infty}dx\frac{\exp(-x\omega_{1}/\omega_{c})}{x^{2}-1}\frac{2x}{e^{\beta\omega_{1}x}-1}.$
(96)

The first term is

$\int_{0}^{\infty}dx\frac{\exp(-x\omega_{1}/\omega_{c})}{x+1}=-e^{\omega_{1}/\omega_{c}}\mathrm{Ei}(-\omega_{1}/\omega_{c})=\left\\{\begin{array}[]{c}\omega_{c}/\omega_{1},\quad\omega_{1}/\omega_{c}\gg
1\\\ -\log(\omega_{1}/\omega_{c}),\quad\omega_{1}/\omega_{c}\ll
1,\end{array}\right.$ (97)

where $\mathrm{Ei}(x)=-\int_{-x}^{\infty}dte^{-t}/t$ is the exponential
integral function [43]. This result defines the value of the function
$\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})$ at zero temperature
$\beta\omega_{1}\rightarrow\infty$. Therefore, the second and the third terms
describe the temperature dependence of
$\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})$. The second term is
asymptotically small $\sim e^{-\beta\omega_{1}}$ at $\beta\omega_{1}\gg 1$,
since the integral

$\text{v.p.}\int_{0}^{\infty}dx\frac{\exp(-x\omega_{1}/\omega_{c})}{x^{2}-1}=\frac{1}{2}\Big{[}e^{\omega_{1}/\omega_{c}}\mathrm{Ei}(-\omega_{1}/\omega_{c})-e^{-\omega_{1}/\omega_{c}}\mathrm{Ei}(\omega_{1}/\omega_{c})\Big{]}$
(98)

is finite for any $\omega_{1}/\omega_{c}$. The third term gives the leading
correction to $\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})$ at small
temperatures

$\displaystyle\text{v.p.}\int_{0}^{\infty}dx\frac{\exp(-x\omega_{1}/\omega_{c})}{x^{2}-1}\frac{2x}{e^{\beta\omega_{1}x}-1}=$
$\displaystyle-\sum_{n=1}^{\infty}\sum_{\zeta=\pm}e^{\zeta(\omega_{c}^{-1}+\beta
n)\omega_{1}}\mathrm{Ei}(-\zeta(\omega_{c}^{-1}+\beta n)\omega_{1})=$
$\displaystyle=$
$\displaystyle-\frac{\pi^{2}}{3}\frac{1}{\beta^{2}\omega_{1}^{2}}+O\Big{(}\frac{\omega_{1}/\omega_{c}}{\beta^{3}\omega_{1}^{3}}\Big{)}.$
(99)

Therefore, for the case $\omega_{c}\beta\gg 1$ we have the following
approximation formula

$\mathcal{F}(\omega_{1}/\omega_{c},\beta\omega_{1})\approx-e^{\omega_{1}/\omega_{c}}\mathrm{Ei}(-\omega_{1}/\omega_{c})-\frac{\pi^{2}}{3\beta^{2}\omega_{1}^{2}}.$
(100)

The exact and approximation functions are presented in Fig. 5 for the case of
$\omega_{1}/\omega_{c}=0.1$.

Figure 5:  Dimensionless function $\mathcal{F}(\omega_{1}/\omega_{c},1/T)$
(blue curve), and its approximation
$-e^{\omega_{1}/\omega_{c}}\mathrm{Ei}(-\omega_{1}/\omega_{c})-\pi^{2}T^{2}/3$
(yellow curve) as a function of dimensionless temperature
$T=1/\beta\omega_{1}\in[0,0.4]$ for $\omega_{1}/\omega_{c}=0.1$. The relative
deviation of the approximation function from the exact one is less than $5\%$
for considered domain.

Summarizing, one can write the expression for the coherence of the second spin
at small temperature

$\langle\sigma_{2}^{x}\rangle=\frac{4Af_{1}f_{2}t}{\omega_{2}}\tanh\left(\frac{\beta\omega_{2}}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)\Big{[}-e^{\omega_{1}/\omega_{c}}\mathrm{Ei}(-\omega_{1}/\omega_{c})-\frac{\pi^{2}}{3\beta^{2}\omega_{1}^{2}}\Big{]}.$
(101)

Note that there is the linear dependence between the the second
$\langle\sigma_{2}^{x}\rangle$ and first spin coherences
$\langle\sigma_{1}^{x}\rangle$ in the system, valid for any spectral density
function $\mathcal{I}(\xi)$

$\langle\sigma_{2}^{x}\rangle=-\frac{t}{\omega_{2}}\tanh\left(\frac{\beta\omega_{2}}{2}\right)\langle\sigma_{1}^{x}\rangle.$
(102)

The derivation of the first spin coherence $\langle\sigma_{1}^{x}\rangle$ is
presented in the Appendix F.

One can point out that at small temperatures $\beta\omega_{2}\gg 1$ the
coherence $\langle\sigma_{2}^{x}\rangle\propto 1/\omega_{2}$ formally diverges
at small $\omega_{2}$. The qualitative analysis of this situation and the
restrictions for values of $\omega_{1}$ and $\omega_{2}$ in perturbative
approach are performed in the Appendix G.

The second spin coherence in the high-temperature limit $\beta\rightarrow 0$
for an arbitrary spectral density function $\mathcal{I}(\xi)$ is

$\displaystyle\langle\sigma_{2}^{x}\rangle|_{\beta\rightarrow
0}=\frac{f_{1}f_{2}t\omega_{1}\beta^{3}}{6}\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)}{\xi}+O(\beta^{5}).$
(103)

The expressions for the coherences for arbitrary values of
$\omega_{1},\omega_{2}$ and $t$ are presented in the Appendices H and I.

##  Appendix F Detuned spins: the first spin coherence

We use the density operator $\varrho_{S}$ defined before for calculation of
the coherence of the first spin
$\langle\sigma_{1}^{x}\rangle=\langle\sigma_{1}^{x}\rangle-\langle\sigma_{1}^{x}\rangle_{0}$.
Repeating all the steps of the calculation performed for the second spin we
obtain that such a coherence is defined by the function
$G_{2}(\tau,\tau^{\prime})=\mathrm{Tr}[e^{-\beta
H_{12}}\widetilde{S}(\tau)\widetilde{S}(\tau^{\prime})\sigma_{1}^{x}+\widetilde{S}^{\dagger}(\tau^{\prime})\widetilde{S}^{\dagger}(\tau)e^{-\beta
H_{12}}\sigma_{1}^{x}]$, that is an analog of the function
$F_{2}(\tau,\tau^{\prime})$ considered before. For arbitrary angle $\theta$
this function is

$\displaystyle G_{2}(\tau,\tau^{\prime})=$ $\displaystyle
4f_{1}f_{2}\Big{[}-\sin(\theta)\sin^{2}(\phi)\\{\sinh\left(\beta\epsilon+\tau^{\prime}(\omega-\epsilon)\right)+\sinh(\beta\omega-\tau(\omega-\epsilon))\\}-$
$\displaystyle-\sin(\theta)\cos^{2}(\phi)\\{\sinh\left(\beta\epsilon-\tau^{\prime}(\omega+\epsilon)\right)-\sinh(\beta\omega-\tau(\omega+\epsilon))\\}-$
$\displaystyle-\sin^{2}(\phi)\\{\sinh\left(\beta\omega+\tau^{\prime}(\epsilon-\omega)\right)+\sinh(\beta\epsilon+\tau(\omega-\epsilon))\\}+$
$\displaystyle+\cos^{2}(\phi)\\{\sinh(\beta\epsilon-\tau(\omega+\epsilon))-\sinh\left(\beta\omega-\tau^{\prime}(\omega+\epsilon)\right)\\}+$
$\displaystyle+\cos(\theta)\sin(2\phi)\\{\sinh(\omega(\beta-\tau))\cosh\left(\epsilon\left(\tau-2\tau^{\prime}\right)\right)+$
$\displaystyle\qquad\qquad\qquad\quad+\sinh\left(\omega\tau^{\prime}\right)\cosh\left(\epsilon\left(\beta+\tau^{\prime}-2\tau\right)\right)\\}\Big{]}.$
(104)

Substituting the decomposition

$\displaystyle\sin^{2}\theta\approx$ $\displaystyle
1-\frac{4t^{2}}{\delta^{2}},\quad\cos^{2}\theta\approx\frac{4t^{2}}{\delta^{2}},$
$\displaystyle\cos^{2}\phi=$
$\displaystyle\frac{1+\sin\theta}{2}\approx\frac{1+\eta}{2}-\frac{\eta
t^{2}}{\delta^{2}},\quad\sin^{2}\phi=\frac{1-\sin\theta}{2}\approx\frac{1-\eta}{2}+\frac{\eta
t^{2}}{\delta^{2}},$ (105)

and keeping the leading order of this function by small parameter $t/\delta$
we get

$\displaystyle
G_{2}(\tau,\tau^{\prime})=8f_{1}f_{2}\sinh\Big{(}\frac{1}{2}\omega_{1}\big{(}\tau^{\prime}-\tau\big{)}\Big{)}$
$\displaystyle\Big{[}\cosh\left(\frac{1}{2}\left(\beta\left(\omega_{2}-\omega_{1}\right)+\omega_{1}\left(\tau^{\prime}+\tau\right)\right)\right)+$
$\displaystyle+\cosh\left(\frac{1}{2}\left(\beta\left(\omega_{1}+\omega_{2}\right)-\omega_{1}\left(\tau^{\prime}+\tau\right)\right)\right)\Big{]}.$
(106)

Integrating over $\tau,\tau^{\prime}$ we get

$\displaystyle G_{2}(\xi)-G_{2}(-\xi)\approx$
$\displaystyle\frac{32f_{1}f_{2}}{\xi\left(\xi^{2}-\omega_{1}^{2}\right)}\cosh\left(\frac{\beta\omega_{2}}{2}\right)\cosh\left(\frac{\beta\omega_{1}}{2}\right)\left[\omega_{1}-\xi\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)\right].$
(107)

Therefore, the coherence of the first spin in the leading order is

$\langle\sigma_{1}^{x}\rangle=-4f_{1}f_{2}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth(\beta\xi/2)\tanh(\beta\omega_{1}/2)-\omega_{1}}{\xi\left(\xi^{2}-\omega_{1}^{2}\right)}.$
(108)

This is nothing but the coherence of the single spin coupled to the bosonic
system, derived in [40].

Note that this result is valid only at small $t/\delta\ll 1$, since it doesn’t
contain the contribution from the second spin. In order to check the validity
of the obtained formula for the case when the ratio $t/\delta$ is not small we
analyze the first spin coherence as well as the second spin coherence for
arbitrary $t,\omega_{1},\omega_{2}$ in the Appendices H and I.

##  Appendix G Qualitative analysis of the spin-boson system coherences

In order to understand the proportionality between the coherences of the first
and second spin

$\langle\sigma_{2}^{x}\rangle=-\frac{t}{\omega_{2}}\tanh\left(\frac{\beta\omega_{2}}{2}\right)\langle\sigma_{1}^{x}\rangle$
(109)

and physical meaning of such proportionality we provide the following
qualitative analysis. Let us consider the bath with the Hamiltonian

$H_{B}=\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k},$ (110)

and system of two weakly interacting spins

$H_{12}=\frac{\omega_{1}}{2}\sigma_{1}^{z}+\frac{\omega_{2}}{2}\sigma_{2}^{z}+t(\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{2}^{+}\sigma_{1}^{-}),$
(111)

at temperature $T=1/\beta$. When two systems are far from each other, both of
them have corresponding thermal distributions. Taking into account that the
interaction between the spins are week ($t\ll\omega_{1},\omega_{2}$), we can
estimate their averages of
$\langle\sigma_{j}^{x}\rangle=\langle\sigma_{j}^{y}\rangle=0$ and
$\langle\sigma_{j}^{z}\rangle=-\tanh(\beta\omega_{j}/2)$, for $j=1,2$. When
system of two spins is close to the bath, the interaction term start to play
an important role

$H_{B1}=(f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x})\sum_{k}\lambda_{k}(b_{k}+b_{k}^{\dagger}).$
(112)

Due to this term the first spin start to “polarize” the bath system. This
process can be observed from the Hamiltonian

$H_{BS_{1}}=\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}+f_{1}\langle\sigma_{1}^{z}\rangle\sum_{k}\lambda_{k}(b_{k}+b_{k}^{\dagger})+\frac{\omega_{1}}{2}\langle\sigma_{1}^{z}\rangle,$
(113)

where we replace the operator $\sigma_{1}^{z}$ by its average value
$\langle\sigma_{1}^{z}\rangle$. Due to the interaction of the first spin with
the bath the boson operators get the non-zero average values $\langle
b_{k}+b_{k}^{\dagger}\rangle\neq 0$. This values can be calculated by
diagonalization of the boson part of the obtained Hamiltonian. It can be done
by substitution
$b_{k}=\beta_{k}-f_{1}\langle\sigma_{1}^{z}\rangle\lambda_{k}/\Omega_{k}$,
where $\beta_{k}$ is the set of boson annihilation operators with the same
commutation rules as $b_{k}$ operators, and zero average values. Hence, we
have

$\langle
b_{k}+b_{k}^{\dagger}\rangle=\langle\beta_{k}+\beta_{k}^{\dagger}\rangle-2f_{1}\langle\sigma_{1}^{z}\rangle\lambda_{k}/\Omega_{k}=-2f_{1}\langle\sigma_{1}^{z}\rangle\lambda_{k}/\Omega_{k}.$
(114)

The non-zero $\langle b_{k}+b_{k}^{\dagger}\rangle$ affects the term

$f_{2}\sigma_{1}^{x}\sum_{k}\lambda_{k}(b_{k}+b_{k}^{\dagger})\rightarrow-2f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\sigma_{1}^{x}\sum_{k}\lambda_{k}^{2}/\Omega_{k}=-2f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\Omega\sigma_{1}^{x}.$
(115)

Therefore, the back-action of the bath modifies the Hamiltonian of the first
spin

$H_{S_{1}}^{mod}=\frac{\omega_{1}}{2}\sigma_{1}^{z}-2f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\Omega\sigma_{1}^{x}.$
(116)

Due to non-zero $\sigma_{1}^{x}$ term, the first spin gets non-zero coherence

$\displaystyle\langle\sigma_{1}^{x}\rangle$ $\displaystyle=Tr[e^{-\beta
H_{S_{1}}^{mod}}\sigma_{1}^{x}]=$ $\displaystyle-$
$\displaystyle\tanh\Big{(}\beta\sqrt{(2f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\Omega)^{2}+(\omega_{1}/2)^{2}}\Big{)}\Big{[}\frac{-2f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\Omega}{\sqrt{(2f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\Omega)^{2}+(\omega_{1}/2)^{2}}}\Big{]}\approx-\frac{4f_{1}f_{2}\langle\sigma_{1}^{z}\rangle^{2}\Omega}{\omega_{1}}.$
(117)

The final approximate result is obtained for the case
$4f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\Omega\ll\omega_{1}$. This means that
the interaction with the bath gives negligibly small correction to the energy
and coherence of the first spin. In other words the latter result corresponds
to weak interaction limit.

The non-zero $\langle\sigma_{1}^{x}\rangle$ induces the effective
$\sigma_{2}^{x}$ term in the second spin system

$t(\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{2}^{+}\sigma_{1}^{-})\rightarrow
t(\langle\sigma_{1}^{+}\rangle\sigma_{2}^{-}+\sigma_{2}^{+}\langle\sigma_{1}^{-}\rangle)=t\frac{\langle\sigma_{1}^{x}\rangle}{2}\sigma_{2}^{x}.$
(118)

Therefore the Hamiltonian of the second spin is modified

$H_{S_{2}}^{mod}=\frac{\omega_{2}}{2}\sigma_{2}^{z}+t\frac{\langle\sigma_{1}^{x}\rangle}{2}\sigma_{2}^{x},$
(119)

which leads to the non-zero coherence

$\displaystyle\langle\sigma_{2}^{x}\rangle=$ $\displaystyle Tr[e^{-\beta
H_{S_{2}}^{mod}}\sigma_{2}^{x}]=$ $\displaystyle-$
$\displaystyle\tanh\Big{(}\beta\sqrt{(t\langle\sigma_{1}^{x}\rangle/2)^{2}+(\omega_{2}/2)^{2}}\Big{)}\frac{t\langle\sigma_{1}^{x}\rangle}{2\sqrt{(t\langle\sigma_{1}^{x}\rangle/2)^{2}+(\omega_{2}/2)^{2}}}\approx-\frac{t}{\omega_{2}}\tanh\Big{(}\frac{\beta\omega_{2}}{2}\Big{)}\langle\sigma_{1}^{x}\rangle.$
(120)

The latter approximate expression is valid only for the case
$t\langle\sigma_{1}^{x}\rangle\ll\omega_{2}$. Therefore the linear dependence
between the coherences of the second and first spins appears if two additional
restrictions are satisfied a) the weak coupling regime
$4f_{1}f_{2}\langle\sigma_{1}^{z}\rangle\Omega\ll\omega_{1}$ and b) the weak
coupling between spins $t\langle\sigma_{1}^{x}\rangle\ll\omega_{2}$. Note that
both inequalities are natural for the derivations based on the thermodynamical
perturbation approach presented before. The first inequality allows to take
only the leading term in perturbation series to evaluate the result. If this
inequality is not satisfied one needs to use the next terms of perturbation
series to get the correct answer. The second condition defines the domain of
the parameters where the linear dependence between both coherences exists.

##  Appendix H General case: the first spin coherence

Let us write the expression for $G_{2}(\tau,\tau^{\prime})$ in terms of
$\theta$ variable only

$\displaystyle G_{2}(\tau,\tau^{\prime})=4f_{1}f_{2}$
$\displaystyle\Big{[}G^{(0)}_{2}(\tau,\tau^{\prime})+\sin\theta
G^{(1)}_{2}(\tau,\tau^{\prime})+\cos(2\theta)2G^{(2)}_{2}(\tau,\tau^{\prime})\Big{]}=$
$\displaystyle=4f_{1}f_{2}$
$\displaystyle\Big{[}\frac{\cos(2\theta)}{2}\Big{\\{}\sinh(\tau^{\prime}\omega)(\cosh(\epsilon(\beta-2\tau+\tau^{\prime}))-\cosh(\epsilon(\tau^{\prime}-\beta)))+$
$\displaystyle+$
$\displaystyle\sinh(\omega(\beta-\tau))(\cosh(\epsilon(\tau-2\tau^{\prime}))-\cosh(\tau\epsilon))\Big{\\}}+$
$\displaystyle+\sin\theta$
$\displaystyle\Big{\\{}\cosh(\tau\omega)\sinh(\epsilon(\beta-\tau))-\sinh(\tau\epsilon)\cosh(\omega(\beta-\tau))+$
$\displaystyle+$
$\displaystyle\sinh(\tau^{\prime}\epsilon)\cosh(\omega(\tau^{\prime}-\beta))+\cosh(\tau^{\prime}\omega)\sinh(\epsilon(\tau^{\prime}-\beta))\Big{\\}}+$
$\displaystyle+\frac{1}{2}\Big{\\{}$
$\displaystyle\sinh(\tau^{\prime}\omega)(\cosh(\epsilon(\beta-2\tau+\tau^{\prime}))+\cosh(\epsilon(\tau^{\prime}-\beta)))+$
$\displaystyle+$
$\displaystyle\sinh(\omega(\beta-\tau))(\cosh(\epsilon(\tau-2\tau^{\prime}))+\cosh(\tau\epsilon))+$
$\displaystyle+$ $\displaystyle
2\cosh(\tau^{\prime}\epsilon)\sinh(\omega(\tau^{\prime}-\beta))-2\sinh(\tau\omega)\cosh(\epsilon(\beta-\tau))\Big{\\}}\Big{]}.$
(121)

After integration over $\tau$ and $\tau^{\prime}$ we get the following
expression

$G_{2}(\xi)-G_{2}(-\xi)=4f_{1}f_{2}\Big{[}G^{(0)}_{2}(\xi)+\sin\theta
G^{(1)}_{2}(\xi)+\cos(2\theta)G^{(2)}_{2}(\xi)\Big{]},$ (122)

where

$\displaystyle G^{(0)}_{2}(\xi)=$
$\displaystyle\frac{-1}{(\omega^{2}-\epsilon^{2})(\xi^{2}-[\omega+\epsilon]^{2})(\xi^{2}-[\omega-\epsilon]^{2})}\Big{[}\frac{4\omega\cosh(\beta\omega)\left(\left(\epsilon^{2}-\omega^{2}\right)^{2}-\xi^{2}\omega^{2}\right)}{\xi}-$
$\displaystyle-$
$\displaystyle\frac{4\omega\cosh(\beta\epsilon)\left(\xi^{2}\omega^{2}\left(\xi^{2}-\omega^{2}\right)+\epsilon^{4}\left(3\xi^{2}-4\omega^{2}\right)+2\epsilon^{2}\left(\xi^{2}\omega^{2}-\xi^{4}+\omega^{4}\right)+2\epsilon^{6}\right)}{\xi\left(\xi^{2}-4\epsilon^{2}\right)}-$
$\displaystyle-$ $\displaystyle
2\coth\left(\frac{\beta\xi}{2}\right)\Big{(}\frac{\omega\epsilon\sinh(\beta\epsilon)\left(-7\xi^{2}\omega^{2}+\xi^{4}+2\omega^{4}+\epsilon^{2}\left(3\xi^{2}+8\omega^{2}\right)-10\epsilon^{4}\right)}{\xi^{2}-4\epsilon^{2}}-$
$\displaystyle-$
$\displaystyle\sinh(\beta\omega)\left(2\omega^{2}\left(\xi^{2}-\omega^{2}\right)+\epsilon^{2}\left(\omega^{2}-\xi^{2}\right)+\epsilon^{4}\right)\Big{)}\Big{]}$
(123)

$\displaystyle G^{(1)}_{2}(\xi)=$
$\displaystyle\frac{4\epsilon\left(\xi^{2}+\omega^{2}-\epsilon^{2}\right)(\cosh(\beta\omega)+\cosh(\beta\epsilon))}{\xi(\xi^{2}-[\omega+\epsilon]^{2})(\xi^{2}-[\omega-\epsilon]^{2})}-$
$\displaystyle-$
$\displaystyle\frac{4\xi\coth\left(\frac{\beta\xi}{2}\right)\left(2\omega\epsilon\sinh(\beta\omega)+\sinh(\beta\epsilon)\left(\xi^{2}-\omega^{2}-\epsilon^{2}\right)\right)}{\xi(\xi^{2}-[\omega+\epsilon]^{2})(\xi^{2}-[\omega-\epsilon]^{2})}$
(124)

$\displaystyle G^{(2)}_{2}(\xi)=$
$\displaystyle\frac{-2\epsilon}{(\omega^{2}-\epsilon^{2})(\xi^{2}-[\omega+\epsilon]^{2})(\xi^{2}-[\omega-\epsilon]^{2})}\Big{[}-2\xi\omega\epsilon\cosh(\beta\omega)+$
$\displaystyle+$
$\displaystyle\frac{2\cosh(\beta\epsilon)\omega\epsilon\left(\xi^{4}+2\left(\omega^{4}-2\omega^{2}\left(\xi^{2}+\epsilon^{2}\right)+\epsilon^{4}\right)\right)}{\xi\left(\xi^{2}-4\epsilon^{2}\right)}+\coth\left(\frac{\beta\xi}{2}\right)\Big{(}\epsilon\sinh(\beta\omega)\times$
$\displaystyle\times$
$\displaystyle\left(\xi^{2}+\omega^{2}-\epsilon^{2}\right)+\frac{\omega\sinh(\beta\epsilon)\left(3\xi^{2}\omega^{2}-\xi^{4}-2\omega^{4}+\epsilon^{2}\left(\xi^{2}+8\omega^{2}\right)-6\epsilon^{4}\right)}{\xi^{2}-4\epsilon^{2}}\Big{)}\Big{]}$
(125)

In order to check this result we consider the limit $t\rightarrow 0$ for
finite value of $|\omega_{1}-\omega_{2}|>0$. In this limit
$\sin\theta\rightarrow\text{sign}(\omega_{1}-\omega_{2})=\eta$ and
$\cos\theta\sim t\rightarrow 0$,
$\cos(2\theta)=2\cos^{2}(\theta)-1\rightarrow-1$. The leading order of this
limit is

$G_{2}(\xi)-G_{2}(-\xi)=4f_{1}f_{2}\Big{[}G^{(0)}_{2}(\xi)+\eta
G^{(1)}_{2}(\xi)-G^{(2)}_{2}(\xi)\Big{]}_{t=0}$ (126)

The $\eta$-independent terms can be simplified

$\displaystyle G^{(0)}_{2}(\xi)-G^{(2)}_{2}(\xi)=$
$\displaystyle\frac{4\omega\left(\xi^{2}-\omega^{2}+\epsilon^{2}\right)(\cosh(\beta\omega)+\cosh(\beta\epsilon))}{\xi\left(\xi^{2}-(\omega+\epsilon)^{2}\right)\left(\xi^{2}-(\omega-\epsilon)^{2}\right)}+$
$\displaystyle+$
$\displaystyle\frac{4\xi\coth\left(\frac{\beta\xi}{2}\right)\left(\sinh(\beta\omega)\left(-\xi^{2}+\omega^{2}+\epsilon^{2}\right)-2\omega\epsilon\sinh(\beta\epsilon)\right)}{\xi\left(\xi^{2}-(\omega+\epsilon)^{2}\right)\left(\xi^{2}-(\omega-\epsilon)^{2}\right)}.$
(127)

Then taking into account that
$\eta\epsilon\rightarrow(\omega_{1}-\omega_{2})/2$, we get

$[G^{(0)}_{2}(\xi)+\eta
G^{(1)}_{2}(\xi)-G^{(2)}_{2}(\xi)]_{t=0}=\frac{8\cosh\left(\frac{\beta\omega_{2}}{2}\right)\cosh\left(\frac{\beta\omega_{1}}{2}\right)\left(\omega_{1}-\xi\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)\right)}{\xi\left(\xi^{2}-\omega_{1}^{2}\right)}$
(128)

Therefore the coherence of the first spin in the leading order by $t$ is

$\displaystyle\langle\sigma_{1}^{x}\rangle_{t=0}=$
$\displaystyle\frac{4f_{1}f_{2}}{2Z_{S}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)[G^{(0)}_{2}(\xi)+\eta
G^{(1)}_{2}(\xi)-G^{(2)}_{2}(\xi)]\Big{|}_{t=0}=$ $\displaystyle=$
$\displaystyle-4f_{1}f_{2}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)-\omega_{1}}{\xi(\xi^{2}-\omega_{1}^{2})},$
(129)

as it was obtained in the previous section. Note that this expression is valid
only in the case $t\ll|\omega_{1}-\omega_{2}|$. Let us consider the opposite
limit. To this end we put $\omega_{1}=\omega_{2}=\omega$ in general
expression, find the $G^{(0)}_{2}(\xi)-G^{(0)}_{2}(-\xi)$ as a function of
$t$, and then consider the limit $t\rightarrow 0$. In this limit
$\cos\theta=1$ and $\sin\theta=0$. Therefore we get

$G_{2}(\xi)-G_{2}(-\xi)=4f_{1}f_{2}\Big{[}G^{(0)}_{2}(\xi)+G^{(2)}_{2}(\xi)\Big{]}_{\omega_{1}=\omega_{2}}.$
(130)

Then

$\displaystyle\Big{[}G^{(0)}_{2}(\xi)+G^{(2)}_{2}(\xi)\Big{]}_{\omega_{1}=\omega_{2}}=$
$\displaystyle\frac{4\omega}{\left(\omega^{2}-t^{2}\right)\left(\xi^{2}-(t+\omega)^{2}\right)\left(\xi^{2}-(\omega-t)^{2}\right)}\times$
$\displaystyle\times$
$\displaystyle\Big{[}\frac{\cosh(\beta\omega)\left(\omega^{2}\left(\xi^{2}-\omega^{2}\right)+t^{2}\left(\xi^{2}+2\omega^{2}\right)-t^{4}\right)}{\xi}+$
$\displaystyle+$ $\displaystyle\frac{\xi\cosh(\beta
t)\left(\omega^{2}\left(\xi^{2}+6t^{2}\right)-3t^{2}\left(\xi^{2}-t^{2}\right)-\omega^{4}\right)}{\xi^{2}-4t^{2}}+$
$\displaystyle+$
$\displaystyle\coth\left(\frac{\beta\xi}{2}\right)\Big{(}\frac{t\sinh(\beta
t)\left(-5\xi^{2}\omega^{2}+\xi^{4}+\xi^{2}t^{2}-2t^{4}+2\omega^{4}\right)}{\xi^{2}-4t^{2}}-$
$\displaystyle-$
$\displaystyle\omega\sinh(\beta\omega)\left(\xi^{2}+t^{2}-\omega^{2}\right)\Big{)}\Big{]}$
(131)

This expression in the limit $t\rightarrow 0$ is

$\displaystyle\Big{[}G^{(0)}_{2}(\xi)+G^{(2)}_{2}(\xi)\Big{]}_{\omega_{1}=\omega_{2},t=0}=(1+\cosh(\beta\omega))\frac{4\left(\omega-\xi\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega}{2}\right)\right)}{\xi\left(\xi^{2}-\omega^{2}\right)}.$
(132)

Taking into account that $Z_{S}=2(1+\cosh(\beta\omega))$ we get the expression
for the coherence

$\displaystyle\langle\sigma_{1}^{x}\rangle_{t=0}=$
$\displaystyle\frac{4f_{1}f_{2}}{2Z_{S}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)[G^{(0)}_{2}(\xi)+G^{(2)}_{2}(\xi)]_{\omega_{1}=\omega_{2},t=0}=$
$\displaystyle=-4f_{1}f_{2}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega}{2}\right)-\omega}{\xi(\xi^{2}-\omega^{2})},$
(133)

which coincides with the one we got previously. Therefore this result is
relevant for the case $t\ll|\omega_{1}-\omega_{2}|$ as well as for the case
$t\gg|\omega_{1}-\omega_{2}|$, for relatively small value of parameter $t$.

In order to check the applicability of the obtained expression for for the
domains of parameters where $t\sim|\omega_{1}-\omega_{2}|$, we calculate the
ratio
$R_{1}(t,\omega_{1},\omega_{2})=\langle\sigma_{1}^{x}\rangle/\langle\sigma_{1}^{x}\rangle_{t=0}$
as a function of pair of dimensionless parameters
$(\omega_{2}/\omega_{1},t/\omega_{1})$ and find the region where this ratio is
close to unity.

$R_{1}(t,\omega_{1},\omega_{2})=-\frac{\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{\\{}G^{(0)}_{2}(\xi)+\sin\theta
G^{(1)}_{2}(\xi)+\cos(2\theta)G^{(2)}_{2}(\xi)\Big{\\}}}{4[\cosh(\beta\omega)+\cosh(\beta\epsilon)]\tanh\left(\frac{\beta\omega_{1}}{2}\right)\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)-\omega_{1}\coth\left(\frac{\beta\omega_{1}}{2}\right)}{\xi(\xi^{2}-\omega_{1}^{2})}}$
(134)

Then taking the particular case of spectral function
$\mathcal{I}(\xi)=A\xi\exp(-\xi/\omega_{c})$ and introducing the dimensionless
variable $x=\xi/\omega_{1}$ together with the function

$\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\omega_{1}\Big{)}=\int_{0}^{\infty}dxe^{-\frac{\omega_{1}}{\omega_{c}}x}\frac{x\coth\left(\frac{\beta\omega_{1}x}{2}\right)-\coth\left(\frac{\beta\omega_{1}}{2}\right)}{x^{2}-1}$
(135)

we get

$R_{1}(t,\omega_{1},\omega_{2})=-\frac{\omega_{1}^{2}\int_{0}^{\infty}dxxe^{-\frac{\omega_{1}}{\omega_{c}}x}\Big{\\{}G^{(0)}_{2}(\omega_{1}x)+\sin\theta
G^{(1)}_{2}(\omega_{1}x)+\cos(2\theta)G^{(2)}_{2}(\omega_{1}x)\Big{\\}}}{4[\cosh(\beta\omega)+\cosh(\beta\epsilon)]\tanh\left(\frac{\beta\omega_{1}}{2}\right)\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\omega_{1}\Big{)}}.$
(136)

The plots of function $R_{1}(t,\omega_{1},\omega_{2})$ for different
temperatures and ratio $\omega_{1}/\omega_{c}$ are presented on Figs. 6,7 and
8.

Figure 6:  Contour plot of $R_{1}(t,\omega_{1},\omega_{2})$ for
$\omega_{1}/\omega_{c}=0.01$ and different temperatures
$\beta\omega_{1}=0.5,1,2$ (left, middle and right figures respectively) for
$\omega_{2}/\omega_{1}\in[0,10]$ and $t/\omega_{1}\in[0,5]$.  Figure 7:
Contour plot of $R_{1}(t,\omega_{1},\omega_{2})$ for $\omega_{1}/\omega_{c}=1$
and different temperatures $\beta\omega_{1}=0.5,1,2$ (left, middle and right
figures respectively) for $\omega_{2}/\omega_{1}\in[0,10]$ and
$t/\omega_{1}\in[0,5]$.  Figure 8:  Contour plot of
$R_{1}(t,\omega_{1},\omega_{2})$ for $\omega_{1}/\omega_{c}=100$ and different
temperatures $\beta\omega_{1}=0.5,1,2$ (left, middle and right figures
respectively) for $\omega_{2}/\omega_{1}\in[0,10]$ and $t/\omega_{1}\in[0,5]$.

##  Appendix I General case: the second spin coherence

We rewrite the expression for $F_{2}(\tau,\tau^{\prime})$ in terms of $\theta$
variable only

$\displaystyle F_{2}(\tau,\tau^{\prime})=$ $\displaystyle
2f_{1}f_{2}[\cos\theta F^{(1)}_{2}(\tau,\tau^{\prime})-\cos\theta\sin\theta
F^{(2)}_{2}(\tau,\tau^{\prime})]=$ $\displaystyle=$ $\displaystyle
2f_{1}f_{2}\Big{[}\cos\theta\Big{\\{}\cosh(\epsilon(\tau^{\prime}+\beta-2\tau)-\tau^{\prime}\omega)-\cosh(\epsilon(\tau^{\prime}+\beta-2\tau)+\tau^{\prime}\omega)+$
$\displaystyle+\cosh(2\epsilon\tau^{\prime}-\epsilon\tau-\beta\omega+\tau\omega)-\cosh(2\epsilon\tau^{\prime}-\epsilon\tau+\beta\omega-\tau\omega)+$
$\displaystyle+2\sinh(\epsilon\tau^{\prime})\sinh(\omega(\beta-\tau^{\prime}))+2\sinh(\tau\omega)\sinh(\epsilon(\beta-\tau))\Big{\\}}-$
$\displaystyle-\cos\theta\sin\theta\Big{\\{}\cosh(\epsilon(\tau^{\prime}+\beta-2\tau)-\tau^{\prime}\omega)+\cosh(\epsilon(\tau^{\prime}+\beta-2\tau)+\tau^{\prime}\omega)+$
$\displaystyle+\cosh(2\epsilon\tau^{\prime}-\epsilon\tau-\beta\omega+\tau\omega)+\cosh(2\epsilon\tau^{\prime}-\epsilon\tau+\beta\omega-\tau\omega)-$
$\displaystyle-2\cosh(\omega\tau^{\prime})\cosh(\epsilon(\beta-\tau^{\prime}))-2\cosh(\tau\epsilon)\cosh(\omega(\beta-\tau))\Big{\\}}\Big{]}.$
(137)

After integration over $\tau$ and $\tau^{\prime}$ we get the following
expression

$F_{2}(\xi)-F_{2}(-\xi)=2f_{1}f_{2}\cos\theta[F^{(1)}_{2}(\xi)-\sin\theta
F^{(2)}_{2}(\xi)],$ (138)

where

$\displaystyle F^{(1)}_{2}(\xi)=$
$\displaystyle\frac{8\epsilon\frac{\sinh(\beta\omega)}{\xi}\left[\left(\omega^{2}-\epsilon^{2}\right)^{2}-\xi^{2}\left(\omega^{2}+\epsilon^{2}\right)\right]+16\xi\omega\epsilon^{2}\sinh(\beta\epsilon)}{(\omega^{2}-\epsilon^{2})(\xi^{2}-[\omega-\epsilon]^{2})(\xi^{2}-[\omega+\epsilon]^{2})}+$
$\displaystyle+\frac{8\epsilon\omega\coth\left(\frac{\beta\xi}{2}\right)\left(\xi^{2}-\omega^{2}+\epsilon^{2}\right)(\cosh(\beta\omega)-\cosh(\beta\epsilon))}{(\omega^{2}-\epsilon^{2})(\xi^{2}-[\omega-\epsilon]^{2})(\xi^{2}-[\omega+\epsilon]^{2})},$
(139) $\displaystyle F^{(2)}_{2}(\xi)=$
$\displaystyle\frac{8\epsilon\frac{\sinh(\beta\epsilon)}{\xi}\left[\left(\omega^{2}-\epsilon^{2}\right)^{2}-\xi^{2}\left(\omega^{2}+\epsilon^{2}\right)\right]+16\xi\omega\epsilon^{2}\sinh(\beta\omega)}{(\omega^{2}-\epsilon^{2})(\xi^{2}-[\omega-\epsilon]^{2})(\xi^{2}-[\omega+\epsilon]^{2})}+$
$\displaystyle+$
$\displaystyle\frac{8\epsilon^{2}\coth\left(\frac{\beta\xi}{2}\right)\left(\xi^{2}+\omega^{2}-\epsilon^{2}\right)(\cosh(\beta\epsilon)-\cosh(\beta\omega))}{(\omega^{2}-\epsilon^{2})(\xi^{2}-[\omega-\epsilon]^{2})(\xi^{2}-[\omega+\epsilon]^{2})}.$
(140)

Then taking into account that $\epsilon\cos\theta=t$ and
$\eta\epsilon\rightarrow(\omega_{1}-\omega_{2})/2$ in $t\rightarrow 0$ limit
and finite value of $|\omega_{1}-\omega_{2}|>0$ we get the answer for the
coherence of the second spin

$\displaystyle\langle\sigma_{2}^{x}\rangle=$
$\displaystyle\frac{2f_{1}f_{2}\cos(\theta)}{2Z_{S}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)[F^{(1)}_{2}(\xi)-\eta
F^{(2)}_{2}(\xi)]\Big{|}_{t=0}=$ $\displaystyle=$
$\displaystyle\frac{4f_{1}f_{2}t}{\omega_{2}}\tanh\left(\frac{\beta\omega_{2}}{2}\right)\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)-\omega_{1}}{\xi(\xi^{2}-\omega_{1}^{2})},$
(141)

as it was obtained before.

Let us consider the opposite limit $\omega_{1}\rightarrow\omega_{2}=\omega$,
for which $\cos\theta=1$ and $\sin\theta=0$. Then

$[F_{2}(\xi)-F_{2}(-\xi)]_{\omega_{1}=\omega_{2}}=2f_{1}f_{2}F^{(1)}_{2}(\xi)|_{\omega_{1}=\omega_{2}}.$
(142)

Taking the series by $t$ and keeping the leading order of the expression we
get

$F^{(1)}_{2}(\xi)|_{\omega_{1}=\omega_{2}}=\frac{8t}{\xi\omega(\xi^{2}-\omega^{2})}\Big{[}\xi\coth\left(\frac{\beta\xi}{2}\right)(\cosh(\beta\omega)-1)-\omega\sinh(\beta\omega)\Big{]},$
(143)

which gives us

$\langle\sigma_{2}^{x}\rangle_{\omega_{1,2}=\omega}=\frac{4f_{1}f_{2}t}{\omega}\tanh\left(\frac{\beta\omega}{2}\right)\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega}{2}\right)-\omega}{\xi(\xi^{2}-\omega^{2})}.$
(144)

This expression coincides with the previous expression, obtained for the case
of $t\ll|\omega_{1}-\omega_{2}|$.

In order to verify the universality of the latter result we estimate the ratio
of the coherence for arbitrary values $t/\omega_{1}\in[0,10]$ and
$\omega_{2}/\omega_{1}\in[0,9]$ to the the leading term of the coherence in
$t\rightarrow 0$ limit
$R_{2}(t,\omega_{1},\omega_{2})=\langle\sigma_{2}^{x}\rangle/\langle\sigma_{2}^{x}\rangle_{t\rightarrow
0}$

$R_{2}(t,\omega_{1},\omega_{2})=\frac{\big{(}\omega_{2}\cos\theta/[\cosh(\beta\omega)+\cosh(\beta\epsilon)]\big{)}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{\\{}F^{(1)}_{2}(\xi)-\sin\theta
F^{(2)}_{2}(\xi)\Big{\\}}}{8t\tanh\left(\frac{\beta\omega_{2}}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)-\omega_{1}\coth\left(\frac{\beta\omega_{1}}{2}\right)}{\xi(\xi^{2}-\omega_{1}^{2})}},$
(145)

which for the spectral density function
$\mathcal{I}(\xi)=A\xi\exp(-\xi/\omega_{c})$ takes the form

$R_{2}(t,\omega_{1},\omega_{2})=\frac{\omega_{1}^{2}\omega_{2}\cos\theta\int_{0}^{\infty}dx\,xe^{-\frac{\omega_{1}}{\omega_{c}}x}\Big{\\{}F^{(1)}_{2}(\omega_{1}x)-\sin\theta
F^{(2)}_{2}(\omega_{2}x)\Big{\\}}}{8t[\cosh(\beta\omega)+\cosh(\beta\epsilon)]\tanh\left(\frac{\beta\omega_{2}}{2}\right)\tanh\left(\frac{\beta\omega_{1}}{2}\right)\mathcal{F}\Big{(}\frac{\omega_{1}}{\omega_{c}},\beta\omega_{1}\Big{)}},$
(146)

where we introduced dimensionless parameter $x=\xi/\omega_{1}$. The plots of
function $R_{2}(t,\omega_{1},\omega_{2})$ for different temperatures and ratio
$\omega_{1}/\omega_{c}$ are presented on Figs. 9, 10 and 11.

Figure 9:  Contour plot of $R_{2}(t,\omega_{1},\omega_{2})$ for
$\omega_{1}/\omega_{c}=0.01$ and different temperatures
$\beta\omega_{1}=0.5,1,2$ (left, middle and right figures respectively) for
$\omega_{2}/\omega_{1}\in[0,10]$ and $t/\omega_{1}\in[0,5]$.  Figure 10:
Contour plot of $R_{2}(t,\omega_{1},\omega_{2})$ for $\omega_{1}/\omega_{c}=1$
and different temperatures $\beta\omega_{1}=0.5,1,2$ (left, middle and right
figures respectively) for $\omega_{2}/\omega_{1}\in[0,10]$ and
$t/\omega_{1}\in[0,5]$.  Figure 11:  Contour plot of
$R_{2}(t,\omega_{1},\omega_{2})$ for $\omega_{1}/\omega_{c}=100$ and different
temperatures $\beta\omega_{1}=0.5,1,2$ (left, middle and right figures
respectively) for $\omega_{2}/\omega_{1}\in[0,10]$ and $t/\omega_{1}\in[0,5]$.

##  Appendix J Fermionization of the spin-boson model

Let us consider the Hamiltonian of two spins interacting with the bosons

$H=\frac{\omega_{2}}{2}\sigma_{2}^{z}+\frac{\omega_{1}}{2}\sigma_{1}^{z}-t(\sigma_{2}^{+}\sigma_{1}^{-}+\sigma_{1}^{+}\sigma_{2}^{-})+\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}+[f_{1}\sigma_{1}^{z}-f_{2}\sigma_{1}^{x}]\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}).$
(147)

Note that we change the sign in spin-spin interaction $t\rightarrow-t$ and
spin-boson interaction $f_{2}\rightarrow-f_{2}$ parameters for simplicity in
further calculations. It was shown before that the perturbation theory
calculations with this Hamiltonian are bulky due to non-commutativity of spin
operators and absence of the Wick theorem for them. Therefore, it looks
interesting to rewrite the spin Hamiltonian in terms of bosons or fermions
where the Wick theorem is applicable. This allows us to use the Green’s
function technique and, potentially, make the final result more clear for
understanding. As an example we fermionize the Hamiltonian with the help of
the Jordan-Wigner transformation [44]. In order to do it we introduce the
additional (“zeroth”) spin subsystem, which is described by the operators
$\\{\sigma_{0}^{x},\sigma_{0}^{y},\sigma_{0}^{z}\\}$ and has zero energy. Then
we replace the spin operators by the fermion ones

$\displaystyle\sigma_{0}^{-}$
$\displaystyle=a_{0},\quad\sigma_{0}^{+}=a_{0}^{\dagger},\quad\sigma_{0}^{z}=2a_{0}^{\dagger}a_{0}-1,$
(148) $\displaystyle\quad\sigma_{1}^{-}$
$\displaystyle=(2a_{0}^{\dagger}a_{0}-1)a_{1},\quad\sigma_{1}^{+}=(2a_{0}^{\dagger}a_{0}-1)a_{1}^{\dagger},\quad\sigma_{1}^{z}=2a_{1}^{\dagger}a_{1}-1,$
(149) $\displaystyle\sigma_{2}^{-}$
$\displaystyle=(2a_{1}^{\dagger}a_{1}-1)(2a_{0}^{\dagger}a_{0}-1)a_{2},\quad\sigma_{2}^{+}=(2a_{1}^{\dagger}a_{1}-1)(2a_{0}^{\dagger}a_{0}-1)a_{2}^{\dagger},\quad\sigma_{2}^{z}=2a_{2}^{\dagger}a_{2}-1,$
(150)

and make the additional substitution

$\displaystyle\frac{\omega_{1}}{2}\sigma_{1}^{z}=$
$\displaystyle\omega_{1}(a_{1}^{\dagger}a_{1}-\frac{1}{2})\rightarrow\omega_{1}a_{1}^{\dagger}a_{1},\quad\frac{\omega_{2}}{2}\sigma_{2}^{z}=\omega_{2}(a_{2}^{\dagger}a_{2}-\frac{1}{2})\rightarrow\omega_{2}a_{2}^{\dagger}a_{2},$
(151) $\displaystyle\sigma_{1}^{x}=$
$\displaystyle(\sigma_{1}^{+}+\sigma_{1}^{-})\rightarrow(\sigma_{1}^{+}\sigma_{0}^{-}+\sigma_{0}^{+}\sigma_{1}^{-})=$
$\displaystyle=$
$\displaystyle(2a_{0}^{\dagger}a_{0}-1)a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}(2a_{0}^{\dagger}a_{0}-1)a_{1}=-(a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1}).$
(152)

Then the Hamiltonian takes the following form

$H\rightarrow\omega_{2}a_{2}^{\dagger}a_{2}+\omega_{1}a_{1}^{\dagger}a_{1}+t(a_{2}^{\dagger}a_{1}+a_{1}^{\dagger}a_{2})+\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}+[2f_{1}a_{1}^{\dagger}a_{1}+f_{2}(a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1})]\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}).$
(153)

Note that the fermionization procedure changes the sign of the spin exchange
terms. The coherence of the first spin can be written in terms of
$\mathrm{T}_{\tau}$-ordered fermion Green’s function
$G_{\alpha\beta}(\tau)=-\langle\mathrm{T}_{\tau}a_{\alpha}(\tau)a_{\beta}^{\dagger}(0)\rangle$
as

$\langle\sigma_{1}^{x}\rangle|_{t,f_{2}}\rightarrow-\langle
a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1}\rangle|_{-t,-f_{2}}=-[G_{01}(0_{-})+G_{10}(0_{-})]|_{-t,-f_{2}}.$
(154)

Here subscripts define the correspondence rule between spin and fermion models
— the coherence in the spin system with parameters $t,f_{2}$ corresponds to
the coherence in the fermion system with parameters $-t,-f_{2}$ and vice
versa. We introduce the second spin coherence by analogy

$\langle\sigma_{2}^{x}\rangle|_{t,f_{2}}\rightarrow-\langle
a_{2}^{\dagger}a_{0}+a_{0}^{\dagger}a_{2}\rangle|_{-t,-f_{2}}=-[G_{02}(0_{-})+G_{20}(0_{-})]|_{-t,-f_{2}}.$
(155)

Therefore we need to calculate the fermion Green’s function for the latter
(fermionized) Hamiltonian.

##  Appendix K Green’s function method

It is convenient to present the Hamiltonian

$H=\omega_{1}a_{1}^{\dagger}a_{1}+\omega_{2}a_{2}^{\dagger}a_{2}+t(a_{1}^{\dagger}a_{2}+a_{2}^{\dagger}a_{1})+\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}+[2f_{1}a_{1}^{\dagger}a_{1}+f_{2}(a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1})]\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})$
(156)

as a sum of three parts: fermion Hamiltonian

$H_{F}=\left[\begin{array}[]{ccc}a_{2}^{\dagger}&a_{1}^{\dagger}&a_{0}^{\dagger}\end{array}\right]\left[\begin{array}[]{ccc}\omega_{2}&t&0\\\
t&\omega_{1}&0\\\ 0&0&0\end{array}\right]\left[\begin{array}[]{ccc}a_{2}\\\
a_{1}\\\ a_{0}\\\
\end{array}\right]=a_{\mu}^{\dagger}[\mathcal{H}_{F}]_{\mu\nu}a_{\nu},$ (157)

where we introduce $3\times 3$ matrix $\mathcal{H}_{F}$ and a summation over
repeated (greek) indices; boson Hamiltonian

$H_{B}=\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k};$ (158)

and fermion-boson interaction

$H_{BF}=[2f_{1}a_{1}^{\dagger}a_{1}+f_{2}(a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1})]\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}).$
(159)

The latter term is considered as a perturbation. It gives the contribution to
the full fermion Green’s function
$G_{\alpha\beta}(\tau-\tau^{\prime})=-\langle\mathrm{T}_{\tau}a_{\alpha}(\tau)a_{\beta}^{\dagger}(\tau^{\prime})\rangle$,
where the averaging is defined as $\langle\star\rangle=\mathrm{Tr}[e^{-\beta
H}\star]/\mathrm{Tr}[e^{-\beta H}]$. We evaluate the contribution to the full
Green’s function in terms of the non-interaction fermion
$G^{(0)}_{\alpha\beta}(\tau-\tau^{\prime})$ and boson
$D^{(0)}(\tau-\tau^{\prime})$ Greens’ functions.

Let us consider the fermion subsystem first. The non-interacting Green’s
function reads

$G^{(0)}_{\alpha\beta}(\tau-\tau^{\prime})=-\langle\mathrm{T}_{\tau}a_{\alpha}(\tau)a_{\beta}^{\dagger}(\tau^{\prime})\rangle_{0},$
(160)

where the averaging is defined as
$\langle\star\rangle_{0}=\mathrm{Tr}[e^{-\beta
H_{F}}\star]/\mathrm{Tr}[e^{-\beta H_{F}}]$. This Green’s function has a
compact form in the Matsubara representation

$G^{(0)}_{\alpha\beta}(ip_{n})=\Big{[}\frac{1}{ip_{n}-\mathcal{H}_{F}}\Big{]}_{\alpha\beta}=U_{\alpha\mu}\frac{\delta_{\mu\nu}}{ip_{n}-\epsilon_{\mu}}[U^{\dagger}]_{\nu\beta}=U_{\alpha\mu}\widetilde{G}^{(0)}_{\alpha\beta}(ip_{n})[U^{\dagger}]_{\nu\beta}.$
(161)

Here $p_{n}=(2n+1)\pi/\beta$ with integer $n$,
$\epsilon_{2}=(\omega_{1}+\omega_{2})/2+\sqrt{(\omega_{1}-\omega_{2})^{2}/4+t^{2}}$,
$\epsilon_{1}=(\omega_{1}+\omega_{2})/2-\sqrt{(\omega_{1}-\omega_{2})^{2}/4+t^{2}}$,
$\epsilon_{0}=0$, $\sin\theta=2t/\sqrt{(\omega_{1}-\omega_{2})^{2}+4t^{2}}$
and
$\cos\theta=(\omega_{2}-\omega_{1})/\sqrt{(\omega_{1}-\omega_{2})^{2}+4t^{2}}$.
The $3\times 3$ matrix $U$ diagonalizes the fermion Hamiltonian matrix
$\mathcal{H}_{F}$

$\mathcal{H}_{F}=U\left[\begin{array}[]{ccc}\epsilon_{2}&0&0\\\
0&\epsilon_{1}&0\\\ 0&0&\epsilon_{0}\end{array}\right]U^{\dagger},\quad
U=\left[\begin{array}[]{ccc}\cos(\theta/2)&-\sin(\theta/2)&0\\\
\sin(\theta/2)&\cos(\theta/2)&0\\\ 0&0&1\end{array}\right].$ (162)

The non-interacting boson Green’s function is

$D^{(0)}(\tau-\tau^{\prime})=-\langle\mathrm{T}_{\tau}\phi(\tau)\phi(\tau^{\prime})\rangle_{0},$
(163)

where we introduced the boson field
$\phi=\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})$ and statistical averaging
$\langle\star\rangle_{0}=\mathrm{Tr}[e^{-\beta
H_{B}}\star]/\mathrm{Tr}[e^{-\beta H_{B}}]$. Then

$\displaystyle D^{(0)}(\tau)=$
$\displaystyle-\sum_{k,l}\lambda_{k}\lambda_{l}\langle\mathrm{T}_{\tau}(b_{k}e^{-\tau\Omega_{k}}+b_{k}^{\dagger}e^{\tau\Omega_{k}})(b_{l}+b_{l}^{\dagger})\rangle_{0}=\sum_{k}\lambda_{k}^{2}D^{(0)}(k,\tau)=$
$\displaystyle=$
$\displaystyle-\theta(\tau)\sum_{k}\lambda_{k}^{2}\Big{[}N_{k}e^{\tau\Omega_{k}}+(N_{k}+1)e^{-\tau\Omega_{k}}\Big{]}-\theta(-\tau)\sum_{k}\lambda_{k}^{2}\Big{[}N_{k}e^{-\tau\Omega_{k}}+(N_{k}+1)e^{\tau\Omega_{k}}\Big{]},$
(164)

where $N_{k}=N(\Omega_{k})=(e^{\beta\Omega_{k}}-1)^{-1}$. The corresponding
Matsubara’s Green’s function with $\omega_{m}=2\pi m/\beta$ and integer $m$
reads

$D^{(0)}(k,i\omega_{m})=\int_{0}^{\beta}d\tau
e^{i\omega_{m}\tau}D^{(0)}(k,\tau)=-\Big{[}\frac{1}{i\omega_{m}+\Omega_{k}}-\frac{1}{i\omega_{m}-\Omega_{k}}\Big{]}=-\frac{2\Omega_{k}}{\omega_{m}^{2}+\Omega_{k}^{2}}.$
(165)

The leading non-zero correction to the fermion Green’s function is

$G_{\alpha\beta}(\tau)\approx G^{(0)}_{\alpha\beta}(\tau)+\Delta
G^{F}_{\alpha\beta}(\tau)+\Delta G^{H}_{\alpha\beta}(\tau).$ (166)

It consists of Fock

$\Delta
G^{F}_{\alpha\beta}(\tau)=-\int_{0}^{\beta}d\tau_{1}\int_{0}^{\beta}d\tau_{2}D^{(0)}(\tau_{1}-\tau_{2})G^{(0)}_{\alpha\alpha^{\prime}}(\tau-\tau_{1})M_{\alpha^{\prime}\beta^{\prime}}G^{(0)}_{\beta^{\prime}\alpha^{\prime\prime}}(\tau_{1}-\tau_{2})M_{\alpha^{\prime\prime}\beta^{\prime\prime}}G^{(0)}_{\beta^{\prime\prime}\beta}(\tau_{2}).$
(167)

and Hartree terms

$\Delta
G_{\alpha\beta}^{H}(\tau)=M_{\alpha^{\prime\prime}\beta^{\prime\prime}}G^{(0)}_{\beta^{\prime\prime}\alpha^{\prime\prime}}(0)\int_{0}^{\beta}d\tau_{1}G^{(0)}_{\alpha\alpha^{\prime}}(\tau-\tau_{1})M_{\alpha^{\prime}\beta^{\prime}}G^{(0)}_{\beta^{\prime}\beta}(\tau_{1})\int_{0}^{\beta}d\tau_{2}D^{(0)}(\tau_{1}-\tau_{2}).$
(168)

Note that there is a summation over the repeated indices. The matrix
$M_{\alpha\beta}$ is

$M=\left[\begin{array}[]{ccc}0&0&0\\\ 0&2f_{1}&f_{2}\\\
0&f_{2}&0\end{array}\right].$ (169)

##  Appendix L Fock contribution

We transform the correction into the Matsubara representation using the
substitutions

$\displaystyle
G^{(0)}_{\alpha\beta}(\tau)=\frac{1}{\beta}\sum_{p_{n}}e^{-ip_{n}\tau}G^{(0)}_{\alpha\beta}(ip_{n}),\quad
D^{(0)}(\tau)=\frac{1}{\beta}\sum_{k}\lambda_{k}^{2}\sum_{\omega_{m}}e^{-i\omega_{m}\tau}D^{(0)}(k,i\omega_{m}),$
(170)

where $G^{(0)}_{\alpha\beta}(ip_{n})$ and $D^{(0)}(k,i\omega_{m})$ have been
derived before. For the Fock term we have

$\Delta
G_{\alpha\beta}^{F}(ip_{n})=G^{(0)}_{\alpha\alpha^{\prime}}(ip_{n})M_{\alpha^{\prime}\beta^{\prime}}\sum_{k}\lambda_{k}^{2}\Big{[}-\frac{1}{\beta}\sum_{\omega_{m}}D^{(0)}(k,i\omega_{m})G^{(0)}_{\beta^{\prime}\alpha^{\prime\prime}}(ip_{n}-i\omega_{m})\Big{]}M_{\alpha^{\prime\prime}\beta^{\prime\prime}}G^{(0)}_{\beta^{\prime\prime}\beta}(ip_{n}).$
(171)

We evaluate the sum in the square brackets

$\displaystyle-\frac{1}{\beta}\sum_{\omega_{m}}D^{(0)}(k,i\omega_{m})$
$\displaystyle
G^{(0)}_{\beta^{\prime}\alpha^{\prime\prime}}(ip_{n}-i\omega_{m})=$
$\displaystyle=$ $\displaystyle
U_{\beta^{\prime}\mu}\Big{\\{}-\frac{1}{\beta}\sum_{\omega_{m}}D^{(0)}(k,i\omega_{m})\widetilde{G}^{(0)}_{\mu\nu}(ip_{n}-i\omega_{m})\Big{\\}}[U^{\dagger}]_{\nu\alpha^{\prime\prime}}=$
$\displaystyle=$ $\displaystyle
U_{\beta^{\prime}\mu}\Big{\\{}\frac{\delta_{\mu\nu}}{\beta}\sum_{\omega_{m}}\frac{2\Omega_{k}}{\omega_{m}^{2}+\Omega_{k}^{2}}\frac{1}{ip_{n}-i\omega_{m}-\epsilon_{\mu}}\Big{\\}}[U^{\dagger}]_{\nu\alpha^{\prime\prime}}=$
$\displaystyle=$ $\displaystyle
U_{\beta^{\prime}\mu}\delta_{\mu\nu}\Big{[}\frac{N_{k}+n(\epsilon_{\mu})}{ip_{n}-\epsilon_{\mu}+\Omega_{k}}+\frac{N_{k}+1-n(\epsilon_{\mu})}{ip_{n}-\epsilon_{\mu}-\Omega_{k}}\Big{]}[U^{\dagger}]_{\nu\alpha^{\prime\prime}}.$
(172)

and introduce the function

$\displaystyle F(\epsilon_{\mu})=$
$\displaystyle\sum_{k}\lambda_{k}^{2}\Big{[}\frac{N_{k}+n(\epsilon_{\mu})}{ip_{n}-\epsilon_{\mu}+\Omega_{k}}+\frac{N_{k}+1-n(\epsilon_{\mu})}{ip_{n}-\epsilon_{\mu}-\Omega_{k}}\Big{]}=$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{[}\frac{N(\xi)+n(\epsilon_{\mu})}{ip_{n}-\epsilon_{\mu}+\xi}+\frac{N(\xi)+1-n(\epsilon_{\mu})}{ip_{n}-\epsilon_{\mu}-\xi}\Big{]},$
(173)

with $n(\epsilon_{\mu})=(e^{\beta\epsilon_{\mu}}+1)^{-1}$ and
$\mathcal{I}(\xi)=\sum_{k}\lambda_{k}^{2}\delta(\xi-\Omega_{k})$. Then the
Fock correction takes the form

$\Delta
G^{F}_{\alpha\beta}(ip_{n})=U_{\alpha\mu^{\prime}}\widetilde{G}^{(0)}_{\mu^{\prime}\nu^{\prime}}(ip_{n})[U^{\dagger}MU]_{\nu^{\prime}\mu}\delta_{\mu\nu}F(\epsilon_{\mu})[U^{\dagger}MU]_{\nu\mu^{\prime\prime}}\widetilde{G}^{(0)}_{\mu^{\prime\prime}\nu^{\prime\prime}}(ip_{n})[U^{\dagger}]_{\nu^{\prime\prime}\beta}.$
(174)

Then we evaluate the matrix sum

$[\widetilde{M}]_{\alpha\beta}=[U^{\dagger}MU]_{\alpha\beta}=[U^{\dagger}]_{\alpha\mu}M_{\mu\nu}U_{\nu\beta}=\left[\begin{array}[]{ccc}2f_{1}\sin^{2}(\theta/2)&f_{1}\sin\theta&f_{2}\sin(\theta/2)\\\
f_{1}\sin\theta&2f_{1}\cos^{2}(\theta/2)&f_{2}\cos(\theta/2)\\\
f_{2}\sin(\theta/2)&f_{2}\cos(\theta/2)&0\end{array}\right],$ (175)

and calculate the matrix elements

$F_{\nu^{\prime}\mu^{\prime\prime}}=[\widetilde{M}]_{\nu^{\prime}\mu}\delta_{\mu\nu}F(\epsilon_{\mu})[\widetilde{M}]_{\nu\mu^{\prime\prime}}.$
(176)

The matrix $F_{\nu^{\prime}\mu^{\prime\prime}}$ is real valued
$\text{Im}[F_{\nu^{\prime}\mu^{\prime\prime}}]=0$ and symmetric
$F_{\nu^{\prime}\mu^{\prime\prime}}=F_{\mu^{\prime\prime}\nu^{\prime}}$.
Therefore we have

$\displaystyle F_{22}=$ $\displaystyle
f_{1}^{2}\Big{[}4F(\epsilon_{2})\sin^{4}(\theta/2)+F(\epsilon_{1})\sin^{2}(\theta)\Big{]}+f_{2}^{2}F(\epsilon_{0})\sin^{2}(\theta/2),$
(177) $\displaystyle F_{11}=$ $\displaystyle
f_{1}^{2}\Big{[}F(\epsilon_{2})\sin^{2}(\theta)+4F(\epsilon_{1})\cos^{4}(\theta/2)\Big{]}+f_{2}^{2}F(\epsilon_{0})\cos^{2}(\theta/2),$
(178) $\displaystyle F_{00}=$ $\displaystyle
f_{2}^{2}\Big{[}F(\epsilon_{2})\sin^{2}(\theta/2)+F(\epsilon_{1})\cos^{2}(\theta/2)\Big{]},$
(179) $\displaystyle F_{21}=$ $\displaystyle
2f_{1}^{2}\sin\theta\Big{[}F(\epsilon_{2})\sin^{2}(\theta/2)+F(\epsilon_{1})\cos^{2}(\theta/2)\Big{]}+f_{2}^{2}F(\epsilon_{0})\sin(\theta)/2,$
(180) $\displaystyle F_{20}=$ $\displaystyle
2f_{1}f_{2}\sin(\theta/2)\Big{[}F(\epsilon_{2})\sin^{2}(\theta/2)+F(\epsilon_{1})\cos^{2}(\theta/2)\Big{]},$
(181) $\displaystyle F_{10}=$ $\displaystyle
2f_{1}f_{2}\cos(\theta/2)\Big{[}F(\epsilon_{2})\sin^{2}(\theta/2)+F(\epsilon_{1})\cos^{2}(\theta/2)\Big{]},$
(182)

Finally the Fock correction takes the simple form

$\Delta
G^{F}_{\alpha\beta}(ip_{n})=U_{\alpha\mu^{\prime}}\frac{1}{ip_{n}-\epsilon_{\mu^{\prime}}}F_{\mu^{\prime}\mu^{\prime\prime}}\frac{1}{ip_{n}-\epsilon_{\mu^{\prime\prime}}}[U^{\dagger}]_{\mu^{\prime\prime}\beta}.$
(183)

We are interested in the matrix elements $\Delta G^{F}_{01}(\tau=0_{-})$ and
$\Delta G^{F}_{02}(\tau=0_{-})$ which correspond to Fock contribution in
$\langle a_{1}^{\dagger}a_{0}\rangle$ and $\langle
a_{2}^{\dagger}a_{0}\rangle$ respectively. To do it we evaluate the sum

$\Delta
G^{F}_{\alpha\beta}(\tau)=\frac{1}{\beta}\sum_{p_{n}}e^{-ip_{n}\tau}\Delta
G^{F}_{\alpha\beta}(ip_{n}),$ (184)

and take the limit $\tau\rightarrow 0_{-}$ in the final result. We have

$\displaystyle\Delta
G^{F}_{01}(ip_{n})=2f_{1}f_{2}\frac{F(\epsilon_{2})\sin^{2}(\theta/2)+F(\epsilon_{1})\cos^{2}(\theta/2)}{ip_{n}-\epsilon_{0}}\Big{[}\frac{\sin^{2}(\theta/2)}{(ip_{n}-\epsilon_{2})}+\frac{\cos^{2}(\theta/2)}{(ip_{n}-\epsilon_{1})}\Big{]},$
(185) $\displaystyle\Delta
G^{F}_{02}(ip_{n})=f_{1}f_{2}\sin\theta\frac{F(\epsilon_{2})\sin^{2}(\theta/2)+F(\epsilon_{1})\cos^{2}(\theta/2)}{ip_{n}-\epsilon_{0}}\Big{[}\frac{1}{(ip_{n}-\epsilon_{2})}-\frac{1}{(ip_{n}-\epsilon_{1})}\Big{]}.$
(186)

These expressions are simplified the limit $t\rightarrow 0$, which equals to
the case $\sin\theta\rightarrow 0$ and $\cos\theta\rightarrow 1$. Here and
further in the text we consider the case $\omega_{2}>\omega_{1}$ for brevity.
Then in the leading order by small parameter $\sin\theta$ we have

$\displaystyle\Delta
G^{F}_{01}(ip_{n})\approx\frac{2f_{1}f_{2}F(\epsilon_{1})}{(ip_{n}-\epsilon_{0})(ip_{n}-\epsilon_{1})},\quad\Delta
G^{F}_{02}(ip_{n})\approx\sin\theta\frac{f_{1}f_{2}F(\epsilon_{1})}{(ip_{n}-\epsilon_{0})}\Big{[}\frac{1}{(ip_{n}-\epsilon_{2})}-\frac{1}{(ip_{n}-\epsilon_{1})}\Big{]}.$
(187)

One can note the correspondence between the obtained corrections

$\Delta G^{F}_{02}(ip_{n})=\frac{t}{ip_{n}-\epsilon_{2}}\Delta
G^{F}_{01}(ip_{n}).$ (188)

Let us evaluate now $\Delta G^{F}_{01}(0_{-})$ and $\Delta G^{F}_{02}(0_{-})$
in the limit $t\rightarrow 0$

$\displaystyle\Delta
G^{F}_{01}(0_{-})=2f_{1}f_{2}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{[}$
$\displaystyle-\frac{(2N(\xi)+1)\epsilon_{1}+(2n(\epsilon_{1})-1)\xi}{2\epsilon_{1}(\xi^{2}-\epsilon_{1}^{2})}+\frac{(2n(\epsilon_{1})-1)n(\epsilon_{1})}{\epsilon_{1}\xi}+$
$\displaystyle+\frac{[N(\xi)+1]n(\epsilon_{1})}{\xi(\xi-\epsilon_{1})}+\frac{N(\xi)n(\epsilon_{1})}{\xi(\xi+\epsilon_{1})}\Big{]},$
(189) $\displaystyle\Delta G^{F}_{02}(0_{-})=\sin\theta
f_{1}f_{2}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{[}W(\epsilon_{2})-W(\epsilon_{1})\Big{]},$
(190)

where we introduced the function

$\displaystyle W(w)=$
$\displaystyle\frac{1}{2w}\Big{[}\frac{N(\xi)+1-n(\epsilon_{1})}{\epsilon_{1}+\xi}-\frac{N(\xi)+n(\epsilon_{1})}{\xi-\epsilon_{1}}\Big{]}+\frac{N(\xi)+1-n(\epsilon_{1})}{(\epsilon_{1}+\xi)(\epsilon_{1}+\xi-w)}n(\epsilon_{1}+\xi)$
$\displaystyle+$
$\displaystyle\frac{n(w)}{w}\Big{[}\frac{N(\xi)+n(\epsilon_{1})}{w-\epsilon_{1}+\xi}+\frac{N(\xi)+1-n(\epsilon_{1})}{w-\epsilon_{1}-\xi}\Big{]}+\frac{N(\xi)+n(\epsilon_{1})}{(\epsilon_{1}-\xi)(\epsilon_{1}-\xi-w)}n(\epsilon_{1}-\xi).$
(191)

##  Appendix M Hartree contribution

Let us consider the correction which appears due to the Hartree term

$\Delta
G_{\alpha\beta}^{H}(\tau)=M_{\alpha^{\prime\prime}\beta^{\prime\prime}}G^{(0)}_{\beta^{\prime\prime}\alpha^{\prime\prime}}(0)\int_{0}^{\beta}d\tau_{1}G^{(0)}_{\alpha\alpha^{\prime}}(\tau-\tau_{1})M_{\alpha^{\prime}\beta^{\prime}}G^{(0)}_{\beta^{\prime}\beta}(\tau_{1})\int_{0}^{\beta}d\tau_{2}D^{(0)}(\tau_{1}-\tau_{2}).$
(192)

The first multiplier can be evaluated in the following way. We calculate this
term in Matsubara representation

$M_{\alpha^{\prime\prime}\beta^{\prime\prime}}G^{(0)}_{\beta^{\prime\prime}\alpha^{\prime\prime}}(ip_{n})=M_{\alpha^{\prime\prime}\beta^{\prime\prime}}U_{\beta^{\prime\prime}\mu}\frac{\delta_{\mu\nu}}{ip_{n}-\epsilon_{\mu}}[U^{\dagger}]_{\nu\alpha^{\prime\prime}}=\sum_{\mu=0}^{2}\frac{[U^{\dagger}MU]_{\mu\mu}}{ip_{n}-\epsilon_{\mu}}$
(193)

We evaluated the matrix $\widetilde{M}=U^{\dagger}MU$ before. Its diagonal
matrix elements are

$\widetilde{M}_{22}=2f_{1}\sin^{2}(\theta/2),\quad\widetilde{M}_{11}=2f_{1}\cos^{2}(\theta/2),\quad\widetilde{M}_{00}=0.$
(194)

It results in the following

$\displaystyle
M_{\alpha^{\prime\prime}\beta^{\prime\prime}}G^{(0)}_{\beta^{\prime\prime}\alpha^{\prime\prime}}(0)=$
$\displaystyle\frac{1}{\beta}\sum_{p_{n}}e^{-ip_{n}\tau}M_{\alpha^{\prime\prime}\beta^{\prime\prime}}G^{(0)}_{\beta^{\prime\prime}\alpha^{\prime\prime}}(ip_{n})\Big{|}_{\tau=0_{-}}=$
$\displaystyle=$ $\displaystyle
2f_{1}\sin^{2}(\theta/2)n(\epsilon_{2})+2f_{1}\cos^{2}(\theta/2)n(\epsilon_{1}).$
(195)

Then we estimate

$\displaystyle\int_{0}^{\beta}d\tau_{1}G^{(0)}_{\alpha\alpha^{\prime}}(\tau-\tau_{1})M_{\alpha^{\prime}\beta^{\prime}}G^{(0)}_{\beta^{\prime}\beta}(\tau_{1})=$
$\displaystyle\frac{1}{\beta}\sum_{p_{n}}e^{-ip_{n}\tau}U_{\alpha\mu}\frac{\delta_{\mu\nu}}{ip_{n}-\epsilon_{\nu}}\widetilde{M}_{\nu\mu^{\prime}}\frac{\delta_{\mu^{\prime}\nu^{\prime}}}{ip_{n}-\epsilon_{\mu^{\prime}}}[U^{\dagger}]_{\nu^{\prime}\beta}=$
$\displaystyle=$
$\displaystyle\frac{1}{\beta}\sum_{p_{n}}e^{-ip_{n}\tau}R_{\alpha\beta}.$
(196)

We are interested in $\Delta G_{01}^{H}(\tau)$ and $\Delta G_{02}^{H}(\tau)$.
Therefore we calculate the corresponding matrix elements

$R_{01}=\frac{f_{2}}{(ip_{n}-\epsilon_{0})}\Big{[}\frac{\sin^{2}(\theta/2)}{ip_{n}-\epsilon_{2}}+\frac{\cos^{2}(\theta/2)}{ip_{n}-\epsilon_{1}}\Big{]},\quad
R_{02}=\frac{f_{2}\sin\theta}{2(ip_{n}-\epsilon_{0})}\Big{[}\frac{1}{ip_{n}-\epsilon_{2}}-\frac{1}{ip_{n}-\epsilon_{1}}\Big{]}.$
(197)

The latter term reads

$\displaystyle\int_{0}^{\beta}d\tau_{2}D^{(0)}(\tau_{1}-\tau_{2})=$
$\displaystyle\frac{1}{\beta}\sum_{k}\lambda_{k}^{2}\sum_{\omega_{m}}D^{(0)}(k,i\omega_{m})e^{-i\omega_{m}\tau_{1}}\int_{0}^{\beta}d\tau_{2}e^{i\omega_{m}\tau_{2}}=$
$\displaystyle=$
$\displaystyle-2\int_{0}^{\infty}d\xi\frac{\mathcal{I}(\xi)}{\xi}=-2\Omega.$
(198)

As a result we have

$\Delta
G^{H}_{01}(ip_{n})=-\frac{4f_{1}f_{2}\Omega}{(ip_{n}-\epsilon_{0})}\Big{[}\sin^{2}(\theta/2)n(\epsilon_{2})+\cos^{2}(\theta/2)n(\epsilon_{1})\Big{]}\Big{[}\frac{\sin^{2}(\theta/2)}{ip_{n}-\epsilon_{2}}+\frac{\cos^{2}(\theta/2)}{ip_{n}-\epsilon_{1}}\Big{]}.$
(199) $\Delta
G^{H}_{02}(ip_{n})=-\frac{2f_{1}f_{2}\Omega\sin\theta}{ip_{n}-\epsilon_{0}}\Big{[}\sin^{2}(\theta/2)n(\epsilon_{2})+\cos^{2}(\theta/2)n(\epsilon_{1})\Big{]}\Big{[}\frac{1}{ip_{n}-\epsilon_{2}}-\frac{1}{ip_{n}-\epsilon_{1}}\Big{]}.$
(200)

Considering $\epsilon_{0}=0$ and taking the limit $\sin\theta\rightarrow
0,\cos\theta\rightarrow 1$ we get

$\Delta G^{H}_{01}(0_{-})=-\frac{4f_{1}f_{2}\Omega
n(\epsilon_{1})}{\beta}\sum_{p_{n}}\frac{e^{-ip_{n}\tau}}{ip_{n}(ip_{n}-\epsilon_{1})}\Big{|}_{\tau=0_{-}}\\!\\!\\!=\frac{2f_{1}f_{2}}{\epsilon_{1}}\Omega
n(\epsilon_{1})\tanh\Big{(}\frac{\beta\epsilon_{1}}{2}\Big{)},$ (201)
$\displaystyle\Delta G^{H}_{02}(0_{-})=$
$\displaystyle-\frac{2f_{1}f_{2}\sin\theta\Omega
n(\epsilon_{1})}{\beta}\sum_{p_{n}}\frac{e^{-ip_{n}\tau}}{ip_{n}}\Big{[}\frac{1}{ip_{n}-\epsilon_{2}}-\frac{1}{ip_{n}-\epsilon_{1}}\Big{]}\Big{|}_{\tau=0_{-}}\\!\\!\\!=$
$\displaystyle=$ $\displaystyle f_{1}f_{2}\sin\theta
n(\epsilon_{1})\Omega\Big{[}\frac{\tanh\Big{(}\frac{\beta\epsilon_{2}}{2}\Big{)}}{\epsilon_{2}}-\frac{\tanh\Big{(}\frac{\beta\epsilon_{1}}{2}\Big{)}}{\epsilon_{1}}\Big{]},$
(202)

where we took into account the sum

$\frac{1}{\beta}\sum_{p_{n}}\frac{e^{-ip_{n}\tau}}{ip_{n}(ip_{n}-w)}=-\frac{1}{2w}+\frac{e^{-\tau
w}}{w}n(w).$ (203)

##  Appendix N Analysis of the coherences

In this section we analyze the coherence of the first
$\langle\sigma_{1}^{x}\rangle_{F}$ and second
$\langle\sigma_{1}^{x}\rangle_{F}$ spins. Here the subscript $F$ indicates
that we derived the coherences from the fermionized Hamiltonian. In accordance
with the aforementioned fermion-boson correspondence we have

$\displaystyle\langle\sigma_{1}^{x}\rangle_{F}=-\langle
a_{1}^{\dagger}a_{0}+a_{0}^{\dagger}a_{1}\rangle|_{-t,-f_{2}}=$
$\displaystyle-[G_{01}(0_{-})+G_{10}(0_{-})]|_{-t,-f_{2}}=$ $\displaystyle=$
$\displaystyle-2[\Delta G^{F}_{01}(0_{-})+\Delta
G^{H}_{01}(0_{-})]|_{-t,-f_{2}},$ (204)
$\displaystyle\langle\sigma_{2}^{x}\rangle_{F}=-\langle
a_{2}^{\dagger}a_{0}+a_{0}^{\dagger}a_{2}\rangle|_{-t,-f_{2}}=$
$\displaystyle-[G_{02}(0_{-})+G_{20}(0_{-})]|_{-t,-f_{2}}=$ $\displaystyle=$
$\displaystyle-2[\Delta G^{F}_{02}(0_{-})+\Delta
G^{H}_{02}(0_{-})]|_{-t,-f_{2}}.$ (205)

Here we took into account that the non-interaction fermion Green’s function
has zero off-diagonal matrix elements and hence $G_{01}(0_{-})=\Delta
G^{F}_{01}(0_{-})+\Delta G^{H}_{01}(0_{-})$, $G_{02}(0_{-})=\Delta
G^{F}_{02}(0_{-})+\Delta G^{H}_{02}(0_{-})$.

Taking the sum of Fock $\Delta G^{F}_{01}(0_{-})$ and Hartree $\Delta
G^{F}_{01}(0_{-})$ contributions and we get

$\displaystyle\langle\sigma_{1}^{x}\rangle_{F}=$ $\displaystyle-\langle
a_{1}^{\dagger}a_{0}+a_{1}^{\dagger}a_{0}\rangle|_{-t,-f_{2}}=$
$\displaystyle=$
$\displaystyle-4f_{1}f_{2}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{[}-\frac{1}{2}\frac{\tanh(\beta\omega_{1}/2)}{\xi\omega_{1}}+\frac{\xi\coth(\beta\xi/2)\tanh(\beta\omega_{1}/2)-\omega_{1}}{2\xi(\xi^{2}-\omega_{1}^{2})}\Big{]},$
(206)

where we took into account that $\epsilon_{1}\approx\omega_{1}$ in the limit
$t\rightarrow 0$. This expression has a similar integrand as the coherence
derived in the spin representation

$\langle\sigma_{1}^{x}\rangle=-4f_{1}f_{2}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi\coth(\beta\xi/2)\tanh(\beta\omega_{1}/2)-\omega_{1}}{\xi(\xi^{2}-\omega_{1}^{2})}.$
(207)

One can observe that the coherences satisfy the equation

$\langle\sigma_{1}^{x}\rangle=2\langle\sigma_{1}^{x}\rangle_{F}-4f_{1}f_{2}\frac{\tanh(\beta\omega_{1}/2)}{\omega_{1}}\int_{0}^{\infty}d\xi\frac{\mathcal{I}(\xi)}{\xi}=2\langle\sigma_{1}^{x}\rangle_{F}-4f_{1}f_{2}\Omega\frac{\tanh(\beta\omega_{1}/2)}{\omega_{1}},$
(208)

which is valid for any spectral density function $\mathcal{I}(\xi)$. Note that