# Synthesizing and multiplexing 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 a crucial prerequisite for quantum technologies.
Therefore, the robust generation, as autonomous as possible, of quantum
coherence remains the essential problem for developing this field. We consider
a method of synthesizing and multiplexing quantum coherence from a spin
systems without any direct drives only coupled to a bosonic baths. The
previous studies in this field have demonstrated that a back-action of the
bath to the spin subsystem is important to generate it, however, it
simultaneously gives a significant limits to the generated coherence. We
propose a viable approach with the bosonic bath that allows overcoming these
limits by avoiding the destructive effect of the back-action processes. Using
this approach, we suggest advanced synthesis of the quantum coherence non-
perturbatively in the spin-boson coupling parameters of multiple bosonic baths
to increase and multiplex it for upcoming proof-of-principle experiments.

##  1 Introduction

Quantum coherence [1] is a significant and diverse subject of modern quantum
physics, phase estimation and thermodynamics [2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16, 17] and a crucial resource of contemporary quantum
technology, specifically, quantum metrology [18, 19], quantum communication
[20, 21], quantum simulators [22, 23, 24], energy harvesting [25, 26], quantum
thermodynamics [27, 28, 29], and quantum computing [30, 31]. A classical,
external and strong coherent drive typically generates such quantum coherence
as a superposition of energy states. Recently, it has been proposed that there
might be a more autonomous alternative, the quantum coherence, from the
coupling between a basic system, like a two-level system, and a thermal bath
[32]. It used a composite interaction between the two-level system and thermal
bath; in one direction, the incoherent energy of the system pushed the bath
coherently, but simultaneously it could receive that coherence of the bath.
Both interactions must be present to obtain quantum coherence in a single two-
level system without any external drive, just from a coherent interaction with
a bath. It is therefore conceptually different from a coherence for a pair of
two-level systems from thermal baths [33]. It triggered further analysis [34,
35, 36, 37, 38]; however, it is still not a fully explored phenomenon, without
a direct experimental test. For more extensive feasibility, other system-bath
topologies generating and detecting more autonomous coherence have to be
found.

Here, we present two crucial steps toward such experimental verifications,
considering many separate two-level systems to push the bath coherently, many
to receive quantum coherence in parallel and also more baths assisting the
process in parallel. Advantageously, we split the single two-level system used
in Ref. [32] to two separate (drive and output) ones. Using these allowed
topologies providing autonomous quantum coherences without a back-action, we
propose and study an autonomous synthetisation of coherence from many systems,
multiplexing it to many systems and employing jointly different baths to
generate the coherences. From the detailed analysis of these cases, we proved
a significant result: many systems and baths could be used in parallel to
obtain and broadcast autonomous quantum coherences in the experiments.

The paper is organized as follows. In Sec. 2, we propose a general method of
the calculation of spin coherences in systems which contain a bosonic bath
interacting with many spins and verify this method with the previously
obtained results for single spin in Refs. [36, 37, 38]. In Sec. 3, we apply
the new method to the case of two spins, input and output, interacting
separately with the bosonic bath. In Sec. 4, we extend the previous system to
the $M$ input spins and analyze the expression for coherence of single output
spin as a function of $M$. Then, in Sec. 5 we examine the case of $M$ input
and two output spins and consider the correlation effects between the output
spins. In Sec. 6, we explore the case of two bosonic baths coupled separately
to $M$ and $N$ input spins, while the output spin is coupled to both baths. We
develop a systematic scheme and calculate the coherence of the output spin in
a general form, and then discuss the generalization of this problem to larger
number of baths and output spins. Finally, in Sec. 7, we further generalize
our method of calculation by substituting the output spin by an oscillator. In
Sec. 8, we summarize our results, discuss the advantages and limits of the
proposed mechanisms of the generation of the cohherence in the spin-bath
systems. Technical details are presented in 4 Appendices.

##  2 Multi-spin interaction with thermal bath

First, we modify the model considered in the Refs. [36, 37, 38] to open more
possibilities for the synthesization and multiplexing by separating a single
two-level system into driving and receiving two-level systems.

We consider the Hamiltonian of the system $H=H_{B}+H_{S}+H_{SB}=H_{0}+H_{SB}$.
Here

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

is the Hamiltonian of bosonic excitations with the spectrum $\Omega_{k}>0$.
The operators $b_{k},b_{k}^{\dagger}$ satisfy the canonical commutation
relations $[b_{k},b_{q}^{\dagger}]=\delta_{kq}$, where $\delta_{kq}$ is the
Kronecker symbol. Furthermore,

$H_{S}=\sum_{j=1}^{M+1}\frac{\omega_{j}}{2}\sigma_{j}^{z},$ (2)

with $\omega_{j}>0$ is the Hamiltonian of $M+1$ spins, $M$ driving the bath to
get coherence there and the $(M+1)$th receiving the coherence. Finally,

$\displaystyle H_{SB}=$
$\displaystyle\sum_{j=1}^{M+1}\boldsymbol{\sigma}_{j}\cdot\boldsymbol{n}_{j}\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}),$
(3)

describes the interaction of the spins with the bosonic system. Here, we have
introduced the notation
$\boldsymbol{\sigma}\cdot\boldsymbol{n}\equiv\sigma^{x}n^{x}+\sigma^{y}n^{y}+\sigma^{z}n^{z}$,
with the Pauli matrices $\sigma^{x},\sigma^{y},\sigma^{z}$ and vector of the
coupling strength parameters $\boldsymbol{n}=(n^{x},n^{y},n^{z})$.

We are interested in the reduced density matrix of the spin system, which can
be obtained from the full canonical density matrix $\rho\equiv e^{-\beta H}/Z$
by tracing out over the bosonic degrees of freedom

$\rho_{S}=Z^{-1}\text{Tr}_{B}\Big{[}e^{-\beta H}\Big{]},$ (4)

where $Z$ is the partition function of the full system

$Z=\text{Tr}_{S}\Big{[}\text{Tr}_{B}\Big{[}e^{-\beta H}\Big{]}\Big{]}.$ (5)

We evaluate the reduced density matrix using the following method. First we
present the operator exponent in the form

$e^{-\beta H}=e^{-\beta H_{0}}\Big{(}e^{\beta H_{0}}e^{-\beta
H}\Big{)}=e^{-\beta H_{0}}U(\beta).$ (6)

The operator $U(\tau)$ satisfies the differential equation in the domain
$\tau\in[0,\beta]$ with the initial condition $U(0)=1$

$\frac{\partial U(\tau)}{\partial\tau}=-\widetilde{H}_{SB}(\tau)U(\tau),$ (7)

where $\widetilde{H}_{SB}(\tau)\equiv e^{\tau H_{0}}H_{SB}e^{-\tau H_{0}}$.
The solution of the equation can be presented in the form of the
chronologically ordered (in the imaginary Matsubara time $\tau$) exponent

$\displaystyle U(\beta)=$ $\displaystyle
T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}=1-\int_{0}^{\beta}d\tau_{1}\widetilde{H}_{SB}(\tau_{1})+$
$\displaystyle+$
$\displaystyle\int_{0}^{\beta}d\tau_{1}\widetilde{H}_{SB}(\tau_{1})\int_{0}^{\tau_{1}}d\tau_{2}\widetilde{H}_{SB}(\tau_{2})+\dots.$
(8)

Using this result we rewrite the reduced density matrix as

$\rho_{S}=\frac{e^{-\beta
H_{S}}}{Z_{S}}\frac{\Big{\langle}T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}\Big{\rangle}_{B}}{\Big{\langle}\Big{\langle}T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}\Big{\rangle}_{B}\Big{\rangle}_{S}}.$
(9)

Here the averaging procedures over the spin $(S)$ and bosonic $(B)$ degrees of
freedom read

$\displaystyle\langle\star\rangle_{S}\equiv$ $\displaystyle
Z_{S}^{-1}\text{Tr}_{S}\Big{[}e^{-\beta H_{S}}\star\Big{]},$ (10)
$\displaystyle\langle\star\rangle_{B}\equiv$ $\displaystyle
Z_{B}^{-1}\text{Tr}_{B}\Big{[}e^{-\beta H_{B}}\star\Big{]}.$ (11)

with

$Z_{S}\equiv\text{Tr}_{S}\Big{[}e^{-\beta H_{S}}\Big{]},\quad
Z_{B}\equiv\text{Tr}_{B}\Big{[}e^{-\beta H_{B}}\Big{]},$ (12)

being the partition functions of the free spin $(Z_{S})$ and boson $(Z_{B})$
subsystems, respectively. Note that the numerator of Eq. (9) can be presented
as $e^{-H_{eff}}$, the denominator then transforms into
$\text{Tr}_{S}[e^{-H_{eff}}]$, see details in Ref. [39]. However, in the
current study we use another way by evaluating perturbatively the
$T_{\tau}$-ordered exponent.

The reduced density operator $\rho_{S}$ depends on the $M+1$ spin degrees of
freedom. Our idea is based on the observation that the bosonic bath
effectively couples the spin degrees of freedom. Part of these coupling terms
is responsible for the generation of a non-zero value of the
$\sigma^{x}_{M+1}$-operator (coherence), which can be calculated as

$\langle\sigma^{x}_{M+1}\rangle\equiv\text{Tr}_{S}[\rho_{S}\sigma^{x}_{M+1}].$
(13)

The aforementioned spin-spin interaction terms correspond to the leading terms
of the perturbation series for the reduced density matrix $\rho_{S}$ in the
spin-boson coupling parameters and can be deduced from the following
expression in the weak-coupling regime (see derivation in Appendix A)

$\displaystyle\Big{\langle}T_{\tau}\Big{\\{}$ $\displaystyle
e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}\Big{\rangle}_{B}\approx
1+\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\times$ $\displaystyle\times$
$\displaystyle\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\xi,\tau-\tau^{\prime})F(\tau)F(\tau^{\prime}).$
(14)

Here, we have introduced the bosonic spectral density function

$\mathcal{I}(\xi)\equiv\sum_{k}\lambda_{k}^{2}\delta(\xi-\Omega_{k}),$ (15)

the function

$\phi(\xi,\tau-\tau^{\prime})=\frac{e^{(\tau-\tau^{\prime})\xi}}{e^{\beta\xi}-1}+\frac{e^{-(\tau-\tau^{\prime})\xi}}{1-e^{-\beta\xi}},$
(16)

and $\tau$-dependent multi-spin operator

$F(\tau)=e^{\tau
H_{S}}\Big{[}\sum_{j=1}^{M+1}\boldsymbol{\sigma}_{j}\cdot\boldsymbol{n}_{j}\Big{]}e^{-\tau
H_{S}}.$ (17)

Figure 1:  Generation of the spin coherence in two different schemes of the
spin-bath interaction. (a) Self-induced spin coherence method from Refs. [38,
37, 36]. The single spin “polarizes” the bosonic subsystem with spectral
density function $\mathcal{I}(\xi)$ at temperature $T$ via the
$f_{1}\sigma_{1}^{z}$ term. The polarized bosonic bath generates indirectly
the coherence in the same spin $\langle\sigma_{1}^{x}\rangle$ via the
$f_{2}\sigma_{1}^{x}$ term, i.e., as the back reaction to the spin system from
the bosonic bath. (b) Externally induced spin coherence method. The first
(input) spin polarizes the bosonic bath via the $f_{1}\sigma_{1}^{z}$ term.
The polarized bosonic bath generates the coherence in the output spin
$\langle\sigma_{2}^{x}\rangle$ via the $f_{2}\sigma_{2}^{x}$ term, i.e., as
the transfer of spin coherence through the bath to the second spin system.

For a comparison, we first reconsider the case with the following single-spin
Hamiltonian [38, 37, 36]

$H_{BS}=[f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x}]\sum_{k}\lambda_{k}(b_{k}+b_{k}^{\dagger}).$
(18)

It corresponds to the situation of one spin coupled to the bosonic thermal
bath simultaneously via the $\sigma^{z}$ and $\sigma^{x}$ spin operators,
i.e., $\boldsymbol{n}=(f_{2},0,f_{1})$. The corresponding system is depicted
in Fig. 1 a).

Using the result of Appendix A, we can write the leading spin-spin correlation
term of the reduced spin density operator in this case

$\displaystyle F(\tau)F(\tau^{\prime})=\,$
$\displaystyle[f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x}(\tau)][f_{1}\sigma_{1}^{z}+f_{2}\sigma_{1}^{x}(\tau^{\prime})]=$
$\displaystyle=$ $\displaystyle
f_{1}^{2}+f_{2}^{2}\cosh(\omega_{1}(\tau-\tau^{\prime}))+f_{2}^{2}\sinh(\omega_{1}(\tau-\tau^{\prime}))\sigma_{1}^{z}+$
$\displaystyle+$ $\displaystyle
f_{1}f_{2}[\sinh(\omega_{1}\tau^{\prime})-\sinh(\omega_{1}\tau)]\sigma_{1}^{x}+if_{1}f_{2}[\cosh(\omega_{1}\tau^{\prime})-\cosh(\omega_{1}\tau)]\sigma_{1}^{y},$
(19)

In the first line of Eq. (2) we introduced the $\tau$-dependent operators

$\sigma^{j}(\tau)\equiv
e^{\tau\frac{\omega}{2}\sigma^{z}}\sigma^{j}e^{-\tau\frac{\omega}{2}\sigma^{z}},$
(20)

with $j=x,y,z$. Using the algebra of Pauli matrices
$\sigma_{j}\sigma_{k}=\delta_{jk}\mathcal{I}+i\varepsilon_{jkl}\sigma_{l}$,
where $\mathcal{I}$ is a $2\times 2$ unit matrix and $\varepsilon_{jkl}$ is
the Levi-Civita symbol, we evaluated
$\sigma^{x}(\tau)=\cosh(\omega\tau)\sigma^{x}+i\sinh(\omega\tau)\sigma^{y}$,
$\sigma^{y}(\tau)=\cosh(\omega\tau)\sigma^{y}-i\sinh(\omega\tau)\sigma^{x}$,
$\sigma^{z}(\tau)=\sigma^{z}$, and then calculated the coherence
$\langle\sigma_{1}^{x}\rangle$ in the leading order of the perturbation theory

$\displaystyle\langle\sigma_{1}^{x}\rangle\equiv$
$\displaystyle\text{Tr}_{S}\Big{[}\rho_{S}\sigma_{1}^{x}\Big{]}\approx\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\xi,\tau-\tau^{\prime})\langle
F(\tau)F(\tau^{\prime})\sigma_{1}^{x}\rangle_{S}=$ $\displaystyle=$
$\displaystyle-4f_{1}f_{2}\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\left(\xi^{2}-\omega_{1}^{2}\right)}.$
(21)

This answer coincides with the previously obtained one [36, 37, 38]. Note that
the complex structure of the integrand is a result of the dynamical back
reaction of the bosonic bath onto the spin system. This is a consequence of
the special coupling of the spin system to the bosonic bath, containing both
the $\sigma^{x}$ and $\sigma^{z}$ coupling terms. It allows self-induced
coherence through the bath by the spin itself but also creates dynamical terms
that limit the amount of coherence. In order to eliminate these dynamical
terms we introduce two groups of spin systems, where the spins of the first
(second) group interact with the bosonic bath only via the $\sigma^{z}$
($\sigma^{x}$) coupling term. In this configuration the first group of spins
influences the bosonic system as driving spins and then the affected bosonic
bath generates the nonzero coherence in the second group of output spins. The
simplest case of such systems is considered in the next section.

##  3 New two-spin method $(M=1)$

For further development and comparison, we propose the basic case
$\boldsymbol{n}_{1}=(0,0,f_{1})$, $\boldsymbol{n}_{2}=(f_{2},0,0)$, with
$f_{1},f_{2}\in\textrm{R}$. The corresponding system is depicted in Fig. 1 b).
For this case the spin-spin correlation term reads

$\displaystyle F(\tau)F(\tau^{\prime})=$
$\displaystyle[f_{1}\sigma_{1}^{z}+f_{2}\sigma_{2}^{x}(\tau)][f_{1}\sigma_{1}^{z}+f_{2}\sigma_{2}^{x}(\tau^{\prime})]=$
$\displaystyle=$ $\displaystyle
f_{1}^{2}+f_{2}^{2}\cosh(\omega_{2}(\tau-\tau^{\prime}))+f_{2}^{2}\sinh(\omega_{2}(\tau-\tau^{\prime}))\sigma_{2}^{z}+$
$\displaystyle+$ $\displaystyle
f_{1}f_{2}\Big{(}[\cosh(\omega_{2}\tau^{\prime})+\cosh(\omega_{2}\tau)]\sigma_{1}^{z}\sigma_{2}^{x}+i[\sinh(\omega_{2}\tau^{\prime})+\sinh(\omega_{2}\tau)]\sigma_{1}^{z}\sigma_{2}^{y}\Big{)}.$
(22)

Using this result we obtain

$\displaystyle\langle\sigma_{2}^{x}\rangle\equiv\text{Tr}_{S}\Big{[}\rho_{S}\sigma_{2}^{x}\Big{]}\approx-4f_{1}f_{2}\tanh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}\frac{\tanh\left(\frac{\beta\omega_{2}}{2}\right)}{\omega_{2}}\Omega,$
(23)

where we have introduced the quantity
$\Omega\equiv\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)/\xi$, with meaning of the
reorganization energy of the bosonic bath. Note that this result coincides
with the mean field (static) result for the original case in the low-
temperature limit [38].

Let us compare the results: the original self-induced method [38, 36], where a
single spin was both the driving and output system, and the new method
splitting these roles to separate spins. To do it we first put
$\omega_{2}=\omega_{1}=\omega$ in the new method to discuss a resonant case
and rewrite both expressions in the following form

$\displaystyle\langle\sigma_{1}^{x}\rangle=$
$\displaystyle-4f_{1}f_{2}\tanh\left(\frac{\beta\omega}{2}\right)\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)}{\xi}\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)-\omega\coth\left(\frac{\beta\omega}{2}\right)}{\left(\xi^{2}-\omega^{2}\right)}=$
$\displaystyle=-\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)}{\xi}\mathcal{F}_{1}(\beta,\omega,\xi),$
(24) $\displaystyle\langle\sigma_{2}^{x}\rangle=$
$\displaystyle-4f_{1}f_{2}\tanh\Big{(}\frac{\beta\omega}{2}\Big{)}\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)}{\xi}\frac{\tanh\left(\frac{\beta\omega}{2}\right)}{\omega}=-\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)}{\xi}\mathcal{F}_{2}(\beta,\omega,\xi).$
(25)

Taking into account the positivity of $\mathcal{I}(\xi)>0$ and the fact that
$\mathcal{F}_{2}(\beta,\omega,\xi)>\mathcal{F}_{1}(\beta,\omega,\xi)>0$ for
any value of $\beta$ and $\omega$ we obtain that
$|\langle\sigma_{2}^{x}\rangle|>|\langle\sigma_{1}^{x}\rangle|$ for the
resonant case. It demonstrates that the new method of generation of the
coherence using two spins with distributed roles is more effective than the
previously proposed method [36], where a single spin had to play double role,
both to coherently displace the bath by its thermal population and receive
that coherence back.

In order to know the effectiveness of the generation of the coherence with the
new method for a couple of spins we consider a general case of different
frequencies $\omega_{1}\neq\omega_{2}$ with the canonical spectral density
function

$\mathcal{I}(\xi)=\lambda\frac{\xi^{s}}{\omega_{c}^{s-1}}e^{-\xi/\omega_{c}}.$
(26)

The dimensionless parameter $\lambda$ describes the strength of the spectral
density function, while $\omega_{c}$ represents the energy cut-off [32, 40,
41]. Note that for this spectral density the coherence
$\langle\sigma_{2}^{x}\rangle$ (23) can be calculated analytically

$\displaystyle\langle\sigma_{2}^{x}\rangle=$
$\displaystyle-4f_{1}f_{2}\tanh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}\frac{\tanh\left(\frac{\beta\omega_{2}}{2}\right)}{\omega_{2}}\lambda\omega_{c}\Gamma(s),$
(27)

where $\Gamma(s)$ is the Gamma function. The normalized coherence
$\langle\sigma_{2}^{x}\rangle/(-4f_{1}f_{2}\lambda)$ for different values of
the parameter $s=0.5;1;2$ as a function of the dimensionless parameters
$\beta\omega_{1}$ and $\omega_{2}/\omega_{1}$ is presented in Fig. 2.

Figure 2: Density plot of normalized coherence
$\langle\sigma_{2}^{x}\rangle/(-4f_{1}f_{2}\lambda)$ for (a) sub-Ohmic with
$s=0.5$, (b) Ohmic with $s=1$, and (c) super-Ohmic with $s=2$ cases in the
domain $\omega_{2}/\omega_{1}\in[0,10]$ and $\beta\omega_{1}\in[0,5]$. One can
observe that the results for $s=1$ and $s=2$ coincide, which reflects the fact
that $\langle\sigma_{2}^{x}\rangle\propto\Gamma(s)$.

The coherence linearly growths with the strength $\lambda$ of the spectral
density function and its cut-off energy $\omega_{c}$. As a function of the
temperature $T=1/\beta$ the coherence is maximized in the $T\rightarrow 0$
limit

$\displaystyle\langle\sigma_{2}^{x}\rangle|_{T=0}=$
$\displaystyle-4f_{1}f_{2}\lambda\frac{\omega_{c}}{\omega_{2}}\Gamma(s).$ (28)

Finally, the generated coherence for the fixed temperature $T=1/\beta$ can be
increased by increasing the frequency of the first spin $\omega_{1}$ and
decreasing the frequency of the second spin $\omega_{2}$. However, the limit
$\omega_{2}\rightarrow 0$ can’t be taken since it will violate the conditions
of the perturbation theory. Namely, the aforementioned perturbation analysis
is applicable as long as the condition $\omega_{2}\gg 4f_{1}f_{2}\Omega$ is
satisfied, see details in Appendix D. The generalized non-perturbative
analysis of the case with arbitrary $\omega_{2}$ is presented later in Sec. 6.

Figure 3: Density plot of
$\langle\sigma_{2}^{x}\rangle/\langle\sigma_{1}^{x}\rangle$ for (a) sub-Ohmic
with $s=0.5$, (b) Ohmic with $s=1$, and (c) super-Ohmic with $s=2$ cases in
the domain $\omega_{2}/\omega_{1}\in[0,10]$ and $\beta\omega_{1}\in[0,5]$. One
can observe that the new method of generation of the coherence becomes
ineffective comparing to the previous one, i.e.,
$\langle\sigma_{2}^{x}\rangle<\langle\sigma_{1}^{x}\rangle$, only for
relatively large values of the parameter $\omega_{2}/\omega_{1}$.

Note that $\langle\sigma_{1}^{x}\rangle$ and $\langle\sigma_{2}^{x}\rangle$
are non-monotonic functions of the parametr $s$. Therefore it is not obvious
for which parameters $\omega_{2}$, $\beta$ and $s$ the coherence
$\langle\sigma_{2}^{x}\rangle$ is larger than the coherence
$\langle\sigma_{1}^{x}\rangle$. In order to clarify this issue we calculate
the ratio $\langle\sigma_{2}^{x}\rangle/\langle\sigma_{1}^{x}\rangle$ for the
dimensionless parameters $\omega_{2}/\omega_{1}$ and $\beta\omega_{1}$

$\frac{\langle\sigma_{2}^{x}\rangle}{\langle\sigma_{1}^{x}\rangle}=\Big{(}\frac{\omega_{c}}{\omega_{1}}\Big{)}^{s}\frac{\Gamma(s)\tanh\Big{(}\frac{\beta\omega_{1}}{2}\frac{\omega_{2}}{\omega_{1}}\Big{)}}{(\omega_{2}/\omega_{1})}\left[\int_{0}^{\infty}dxx^{s-1}e^{-\omega_{1}x/\omega_{c}}\frac{x\coth\Big{(}\frac{\beta\omega_{1}}{2}\frac{\omega_{2}}{\omega_{1}}\Big{)}-\coth\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}}{x^{2}-1}\right]^{-1},$
(29)

and for three values of the parameter $s=0.5;1;2$. The corresponding plots for
the case $\omega_{1}/\omega_{c}=0.1$ are presented in Fig. 3. The plots
demonstrate that the new method of generation becomes effective for small
ratio $\omega_{2}/\omega_{1}$. It can be understood from analysis of the
expression for the coherence of the second spin Eq. (23)

$\displaystyle\langle\sigma_{2}^{x}\rangle=$
$\displaystyle-4f_{1}f_{2}\tanh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}\frac{\tanh\left(\frac{\beta\omega_{2}}{2}\right)}{\omega_{2}}\Omega.$
(30)

As one can see the larger $\omega_{1}$ corresponds to the larger value of
$0<\tanh(\beta\omega_{1}/2)<1$. From other side the multiplier
$\tanh(\beta\omega_{2}/2)/\omega_{2}$ as a function of $\omega_{2}$ is a
decreasing function. Therefore, the larger value of this multiplier can be
reached at small values of $\omega_{2}$.

##  4 Synthesizing coherence from $M$ spins through a single bath

Due to the new method we consider we can now address the first general
problem, if the coherence can be synthesized from many driving spins, in
general, with specific different couplings and frequencies.

Let us consider the case $\boldsymbol{n}_{j}=(0,0,f_{1})$, with $j=1,2,\dots
M$, $\boldsymbol{n}_{M+1}=(f_{2},0,0)$. Then we have

$\displaystyle F(\tau)F(\tau^{\prime})=$ $\displaystyle
f_{1}f_{2}\sum_{j=1}^{M}[\cosh(\omega_{M+1}\tau^{\prime})+\cosh(\omega_{M+1}\tau)]\sigma_{j}^{z}\sigma_{M+1}^{x}+$
$\displaystyle+$ $\displaystyle
if_{1}f_{2}\sum_{j=1}^{M}[\sinh(\omega_{M+1}\tau^{\prime})+\sinh(\omega_{M+1}\tau)]\sigma_{j}^{z}\sigma_{M+1}^{y}+$
$\displaystyle+$ $\displaystyle
f_{1}^{2}\sum_{j,l=1}^{M}\sigma_{j}^{z}\sigma_{l}^{z}+f_{2}^{2}\cosh(\omega_{M+1}(\tau-\tau^{\prime}))+f_{2}^{2}\sinh(\omega_{M+1}(\tau-\tau^{\prime}))\sigma_{M+1}^{z}.$
(31)

Using this result we calculate
$\langle\sigma_{M+1}^{x}\rangle=\text{Tr}_{S}[\rho_{S}\sigma_{M+1}^{x}]$

$\displaystyle\langle\sigma_{M+1}^{x}\rangle\approx-4f_{1}f_{2}\Big{[}\sum_{j=1}^{M}\tanh\Big{(}\frac{\beta\omega_{j}}{2}\Big{)}\Big{]}\frac{\tanh\left(\frac{\beta\omega_{M+1}}{2}\right)}{\omega_{M+1}}\Omega.$
(32)

Comparing this result with the previous case (23) one concludes that each of
$M$ spins contributes cumulatively to the $(M+1)$th spin coherence
$\langle\sigma_{M+1}^{x}\rangle\propto\sum_{j=1}^{M}\langle\sigma_{j}^{z}\rangle_{S}$.
It means, the bath can equally accumulate the coherence from the same driving
spins and then uses it to make the output spin equivalently coherent even if
the coupling $f_{1}f_{2}$ in Eq. (32) decreases $M$ times.

This result can be generalized for the case of $M$ spins with different
couplings to the bath system, i.e., for the case

$\displaystyle
H_{SB}=\Big{[}\sum_{j=1}^{M}\boldsymbol{\sigma}_{j}\cdot\boldsymbol{n}_{j}+\boldsymbol{\sigma}_{M+1}\cdot\boldsymbol{n}_{M+1}\Big{]}\sigma_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}),$
(33)

with different $\boldsymbol{n}_{j}=(0,0,f_{1}^{(j)})$, for $j=1,2,\dots M$ and
$\boldsymbol{n}_{M+1}=(f_{2},0,0)$. Repeating the previous calculations for
this case we obtain

$\displaystyle\langle\sigma_{M+1}^{x}\rangle\approx-4f_{2}\Big{[}\sum_{j=1}^{M}f_{1}^{(j)}\tanh\Big{(}\frac{\beta\omega_{j}}{2}\Big{)}\Big{]}\frac{\tanh\left(\frac{\beta\omega_{M+1}}{2}\right)}{\omega_{M+1}}\Omega.$
(34)

It is convenient to introduce the density function of the coupling parameters

$\displaystyle
f_{1}(\omega)\equiv\sum_{j=1}^{M}\delta(\omega-\omega_{j})f_{1}^{(j)},$ (35)

and rewrite the expression for the coherence in the form

$\displaystyle\langle\sigma_{M+1}^{x}\rangle\approx-4f_{2}\Big{[}\int_{-\infty}^{\infty}\\!\\!\\!d\omega
f_{1}(\omega)\tanh\Big{(}\frac{\beta\omega}{2}\Big{)}\Big{]}\frac{\tanh\left(\frac{\beta\omega_{M+1}}{2}\right)}{\omega_{M+1}}\Omega.$
(36)

This expression easily demonstrates that the main contribution to the
generated coherence comes from the domain $\omega\gg 2/\beta$, where
$\tanh(\beta\omega/2)\approx 1$. Therefore one can use the simplified formula
for the coherence

$\displaystyle\langle\sigma_{M+1}^{x}\rangle\approx-4f_{2}\Big{[}\int_{2/\beta}^{\infty}d\omega
f_{1}(\omega)\Big{]}\frac{\tanh\left(\frac{\beta\omega_{M+1}}{2}\right)}{\omega_{M+1}}\Omega.$
(37)

Advantageously, even a broad distribution of coupling parameters
$f_{1}(\omega)$ in the frequency domain can be sufficient to induce nearly $M$
times higher coherence if the function $f_{1}(\omega)$ is localized in the
region well above $2/\beta$.

Note that the calculation of the average value of $\sigma_{M+1}^{z}$ is more
complicated, since it needs to take into account the denominator terms in the
reduced density operator $\rho_{S}$ from Eq. (9). Such a derivation is
presented in Appendix B.

##  5 Coherence multiplexing to two spins ($M+2$ case)

Now we are in position to analyse another problem concerning detection of
autonomous coherences, when more spins can be used as output systems. If they
all independently receive the same coherence from a single bath, we can
observe it more easily without need to repeat the experiment in time. For
analysis of coherence multiplexing, we extend the previous study to the case
of $(M+2)$ spins coupled to the bosonic bath, see Fig. 4.

Figure 4:  Synthesizing and multiplexing of autonomous coherence. The first
group of $M$ driving spins polarizes the bosonic bath via the
$\sum_{j=1}^{M}f_{1}\sigma_{j}^{z}$ term. The polarized bosonic bath
synthetically generates the coherence in the second group of the output spins
via the $f_{2}\sigma_{M+p}^{x}$ terms, where index $p=1,2,\dots$ marks the
target spins which get the coherence from the bath. For the coherence
multiplexing, the figure represents the case of two output spins. The
generated coherences of these spins
$\langle\sigma_{M+1}^{x}\rangle=\langle\sigma_{M+2}^{x}\rangle$ are $M$ times
larger than the coherence $\langle\sigma_{2}^{x}\rangle$ (colored in red)
generated in the direct two-spins scheme, see Fig. 1(b) and Eq. (23).

The Hamiltonians of the spin system and spin-boson coupling term are

$H_{S}=\sum_{j=1}^{M}\frac{\omega_{j}}{2}\sigma_{j}^{z}+\sum_{p=1}^{2}\frac{\omega_{M+p}}{2}\sigma_{M+p}^{z}$
(38)

and

$H_{SB}=\Big{[}\sum_{j=1}^{M}\boldsymbol{\sigma}_{j}\cdot\boldsymbol{n}_{j}+\sum_{p=1}^{2}\boldsymbol{\sigma}_{M+p}\cdot\boldsymbol{n}_{M+p}\Big{]}\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}),$
(39)

respectively. We again consider the special case
$\boldsymbol{n}_{j}=(0,0,f_{1})$, with $j=1,2,\dots M$,
$\boldsymbol{n}_{M+1}=\boldsymbol{n}_{M+2}=(f_{2},0,0)$. Repeating the
averaging procedure over the bosonic degrees of freedom, see Appendix A, we
obtain the following spin-spin correlation term

$\displaystyle F(\tau)F(\tau^{\prime})=$ $\displaystyle
f_{1}f_{2}\sum_{j=1}^{M}\sigma_{j}^{z}\sum_{p=1}^{2}[\cosh(\omega_{M+p}\tau^{\prime})+\cosh(\omega_{M+p}\tau)]\sigma_{M+p}^{x}+$
$\displaystyle+$ $\displaystyle
if_{1}f_{2}\sum_{j=1}^{M}\sigma_{j}^{z}\sum_{p=1}^{2}[\sinh(\omega_{M+p}\tau^{\prime})+\sinh(\omega_{M+p}\tau)]\sigma_{M+p}^{y}+$
$\displaystyle+$ $\displaystyle
f_{1}^{2}\sum_{j,l=1}^{M}\sigma_{j}^{z}\sigma_{l}^{z}+f_{2}^{2}\sum_{p=1}^{2}\Big{[}\cosh(\omega_{M+p}(\tau-\tau^{\prime}))+\sinh(\omega_{M+p}(\tau-\tau^{\prime}))\sigma_{M+p}^{z}\Big{]}+$
$\displaystyle+$ $\displaystyle
f_{2}^{2}\sigma_{M+1}^{x}\sigma_{M+2}^{x}[\cosh(\omega_{M+1}\tau)\cosh(\omega_{M+2}\tau^{\prime})+\cosh(\omega_{M+2}\tau)\cosh(\omega_{M+1}\tau^{\prime})]-$
$\displaystyle-$ $\displaystyle
f_{2}^{2}\sigma_{M+1}^{y}\sigma_{M+2}^{y}[\sinh(\omega_{M+1}\tau)\sinh(\omega_{M+2}\tau^{\prime})+\sinh(\omega_{M+2}\tau)\sinh(\omega_{M+1}\tau^{\prime})]+$
$\displaystyle+$ $\displaystyle
if_{2}^{2}\sigma_{M+1}^{x}\sigma_{M+2}^{y}[\cosh(\omega_{M+1}\tau)\sinh(\omega_{M+2}\tau^{\prime})+\sinh(\omega_{M+2}\tau)\cosh(\omega_{M+1}\tau^{\prime})]+$
$\displaystyle+$ $\displaystyle
if_{2}^{2}\sigma_{M+1}^{y}\sigma_{M+2}^{x}[\sinh(\omega_{M+1}\tau)\cosh(\omega_{M+2}\tau^{\prime})+\cosh(\omega_{M+2}\tau)\sinh(\omega_{M+1}\tau^{\prime})].$
(40)

Consequently,

$\displaystyle\langle F(\tau)$ $\displaystyle
F(\tau^{\prime})\sigma_{M+p}^{x}\rangle_{S}=-f_{1}f_{2}\sum_{j=1}^{M}\tanh\Big{(}\frac{\beta\omega_{j}}{2}\Big{)}\times$
$\displaystyle\times$
$\displaystyle\Big{\\{}[\cosh(\omega_{M+p}\tau^{\prime})+\cosh(\omega_{M+p}\tau)]-[\sinh(\omega_{M+p}\tau^{\prime})+\sinh(\omega_{M+p}\tau)]\tanh\Big{(}\frac{\beta\omega_{M+p}}{2}\Big{)}\Big{\\}}.$
(41)

After the substituting this result into the expression for the coherence of
the $(M+p)$th spin one can observe that it has the same structure as the
coherence, obtained in Sec. 4 for the $(M+1)$th spin. Therefore, the level of
the coherence of the $(M+1)$th spin doesn’t decrease the level of the
coherence of $(M+2)$th spin and vice versa. However, the pair of output spins
can be still correlated which will reduce their application as independent
resources. In order to check the correlation between these spins we calculate
the correlation parameter

$\sigma\equiv\sqrt{\langle\sigma_{M+1}^{x}\sigma_{M+2}^{x}\rangle-\langle\sigma_{M+1}^{x}\rangle\langle\sigma_{M+2}^{x}\rangle}.$
(42)

In the leading order in coupling parameters $f_{1},f_{2}$ we get the following
expression (see intermediate results of the calculation in Appendix C)

$\sigma^{2}\approx
4f_{2}^{2}\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)}{\xi}\mathcal{G}(\xi,\omega_{M+1},\omega_{M+2}),$
(43)

where

$\displaystyle\mathcal{G}(\xi,x,y)=$
$\displaystyle\frac{xy\xi\tanh\left(\frac{\beta
x}{2}\right)\tanh\left(\frac{\beta
y}{2}\right)}{\left(\xi^{2}-x^{2}\right)\left(\xi^{2}-y^{2}\right)}\coth\Big{(}\frac{\beta\xi}{2}\Big{)}+\frac{xy\xi^{2}}{(x^{2}-y^{2})}\Big{[}\frac{\tanh\left(\frac{\beta
x}{2}\right)}{y(\xi^{2}-y^{2})}-\frac{\tanh\left(\frac{\beta
y}{2}\right)}{x(\xi^{2}-x^{2})}\Big{]}.$ (44)

One can observe that the integrand $\mathcal{G}(\xi,x,y)$ is a positive
regular function on $\xi\in[0,\infty)$ domain. It is a symmetric function of
the parameters $x,y$ and its values belong to the domain
$\mathcal{G}(\xi,x,y)\in[\mathcal{G}_{min}(x,y),\mathcal{G}_{max}(x,y)]$ with

$\displaystyle\mathcal{G}_{min}(x,y)=$
$\displaystyle\frac{2\tanh\left(\frac{\beta x}{2}\right)\tanh\left(\frac{\beta
y}{2}\right)}{\beta xy},$ (45) $\displaystyle\mathcal{G}_{max}(x,y)=$
$\displaystyle\frac{x\tanh\left(\frac{\beta
x}{2}\right)-y\tanh\left(\frac{\beta y}{2}\right)}{x^{2}-y^{2}}.$ (46)

Therefore we have the following bands on the correlation value

$\mathcal{G}_{min}(\omega_{M+1},\omega_{M+2})\leq\frac{\sigma^{2}}{4f_{2}^{2}\Omega}\leq\mathcal{G}_{max}(\omega_{M+1},\omega_{M+2}).$
(47)

In the high temperature limit
$\mathcal{G}_{min}(\omega_{M+1},\omega_{M+2}),\mathcal{G}_{max}(\omega_{M+1},\omega_{M+2})\rightarrow\beta/2$.
Therefore in this case we have the following answer

$\sigma^{2}=2\beta f_{2}^{2}\Omega.$ (48)

Note that the result has the universal character, i.e., it doesn’t depend on
the parameters of the spin systems $\omega_{M+1},\omega_{M+2}$. In the low-
temperature limit $\beta\omega_{M+1},\beta\omega_{M+2}\gg 1$ the $\sigma^{2}$
has a non-vanishing value

$\sigma^{2}=\frac{4f_{2}^{2}}{\omega_{M+1}+\omega_{M+2}}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\frac{\xi+\omega_{M+1}+\omega_{M+2}}{(\xi+\omega_{M+1})(\xi+\omega_{M+2})}.$
(49)

Note that the non-zero correlation between both spins appears in the leading
order of the perturbation theory and it can be interpreted as the reaction of
the second spin to the first one (and vice versa) using the states of the
bosonic bath as intermediate virtual states. Indeed, the leading term of the
correlation value contains only the coupling parameter $f_{2}$, but not the
coupling parameters $f_{1}$, so the non-zero correlation value will be also at
$f_{1}=0$.

In order to understand the effectiveness of such a scheme, we estimate the
signal-to-noise ratio $|\langle\sigma_{M+1}^{x}\rangle|/\sigma$. We consider
the case of the input and output spins with the same frequencies, namely,
$\omega_{j}=\omega_{1}$ for $j=1,2,\dots M$ and
$\omega_{M+1}=\omega_{M+2}=\omega$, as a typical example

$\displaystyle\frac{|\langle\sigma_{M+1}^{x}\rangle|}{\sigma}=\frac{2f_{1}\Omega
M\tanh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}\tanh\Big{(}\frac{\beta\omega}{2}\Big{)}}{\omega\sqrt{\int_{0}^{\infty}d\xi\frac{\mathcal{I}(\xi)}{\xi}\mathcal{G}(\xi,\omega,\omega)}}.$
(50)

Using the special case of the power-law (generalized Ohmic) spectral density
function (26) we obtain the following expression

$\displaystyle\frac{\langle\sigma_{M+1}^{x}\rangle}{\sigma}=\eta\frac{\Big{(}\frac{\omega_{c}}{\omega_{1}}\Big{)}^{\frac{s+1}{2}}\Gamma(s)\tanh\Big{(}\frac{b}{2}\Big{)}\tanh\Big{(}\frac{bw}{2}\Big{)}}{\sqrt{w}\sqrt{\int_{0}^{\infty}dxx^{s}e^{-x\frac{\omega_{1}}{\omega_{c}}}g(x)}}.$
(51)

Here we have introduced the dimensionless parameter
$\eta=4f_{1}M\sqrt{\lambda}$, which doesn’t depend on the input and output
spin frequencies $\omega_{1},\omega$ and inverse temperature $\beta$,

$\displaystyle g(x)=$
$\displaystyle\frac{\tanh^{2}\left(\frac{bw}{2}\right)\left[bwx\left(w^{2}-x^{2}\right)+4w^{3}\coth\left(\frac{bx}{2}\right)\right]}{\left(w^{2}-x^{2}\right)^{2}}+$
$\displaystyle+$
$\displaystyle\frac{bwx\left(x^{2}-w^{2}\right)+\left(2x^{3}-6w^{2}x\right)\tanh\left(\frac{bw}{2}\right)}{\left(w^{2}-x^{2}\right)^{2}},$
(52)

and short notations $w=\omega/\omega_{1}$, $b=\beta\omega_{1}$ for brevity.

Using the obtained formula we plot in Fig. 5 the normalized signal-to-noise
ratio $\langle\sigma_{M+1}^{x}\rangle/\eta\sigma$ as a function of the
dimensionless parameters $\beta\omega_{1}$ and $\omega/\omega_{1}$ for
different values of $s=0.5;1;2$.

Figure 5: Density plot of normalized signal-to-noise ratio
$\langle\sigma_{M+1}^{x}\rangle/\eta\sigma$ for (a) sub-Ohmic with $s=0.5$,
(b) Ohmic with $s=1$, and (c) super-Ohmic with $s=2$ cases in the domain
$\omega/\omega_{1}\in[0,10]$ and $\beta\omega_{1}\in[0,5]$.

##  6 Synthesizing coherence from two independent baths ($M+N+1$ case)

Now, the last intriguing problem remains for our new method: do we really need
a single bath that is coherently manipulated by $M$ driving spins? Cannot we
use two different and independent baths coherently pushed by $M$ and $N$ spins
if they both transfer coherence through these baths to the same output spin?

We consider an extended version of the previous system, see Fig. 6. We
introduce two reservoirs of bosons, with spectra $\Omega_{k},\Upsilon_{k}>0$
and creation (annihilation) operators
$b_{k}^{\dagger},d_{k}^{\dagger}(b_{k},d_{k})$, respectively. First (second)
reservoir is coupled to $M$ ($N$) spins via
$f_{1}\sigma_{j}^{z}$($g_{1}\sigma_{j}^{z}$) term. Both reservoirs are coupled
to the $(M+N+1)$th spin via $f_{2}\sigma^{x}$ and $g_{2}\sigma^{x}$ terms.

Figure 6:  Autonomous quantum coherence from several independent baths. The
group of first $M$ driving spins polarizes the first bosonic bath via
$\sum_{j=1}^{M}f_{1}\sigma_{j}^{z}$ term, the group of another $N$ driving
spins polarizes the second bosonic bath via
$\sum_{j=M+1}^{M+N}g_{1}\sigma_{j}^{z}$ term. The polarized baths
simultaneously generate the coherence $\langle\sigma^{x}\rangle$ of the
$(M+N+1)$th output spin via $f_{2}\sigma^{x}$ and $g_{2}\sigma^{x}$ terms. The
coherence of this spin in the leading order in the coupling parameters
$f_{1},f_{2},g_{1},g_{2}$ is $(M+N)$ times larger than the coherence
$\langle\sigma_{2}^{x}\rangle$ (colored in red) generated in the direct two-
spins scheme, see Fig. 1(b) and Eq. (23).

The Hamiltonian of the system

$\displaystyle H=H_{S}+H_{B}+H_{SB}$ (53)

consists of the spin Hamiltonian

$\displaystyle H_{S}=$
$\displaystyle\sum_{j=1}^{M}\frac{\omega_{1}}{2}\sigma_{j}^{z}+\sum_{j=M+1}^{M+N}\frac{\omega_{2}}{2}\sigma_{j}^{z}+\frac{\omega}{2}\sigma^{z},$
(54)

the bosnic Hamiltonian

$\displaystyle H_{B}=$
$\displaystyle\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}+\sum_{k}\Upsilon_{k}d_{k}^{\dagger}d_{k},$
(55)

and the spin-bath interaction term

$\displaystyle H_{SB}=$ $\displaystyle
f_{1}\sum_{j=1}^{M}\sigma_{j}^{z}\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})+g_{1}\sum_{j=M+1}^{M+N}\sigma_{j}^{z}\sum_{k}\kappa_{k}(d_{k}^{\dagger}+d_{k})+$
$\displaystyle+$
$\displaystyle\sigma^{x}\sum_{k}\Big{[}f_{2}\lambda_{k}(b^{\dagger}_{k}+b_{k})+g_{2}\kappa_{k}(d^{\dagger}_{k}+d_{k})\Big{]}.$
(56)

The first line of the $H_{SB}$ describes the interaction of the $M$ input
spins with the left bath, the second line describes the interaction of the
additional $N$ input spins with the right bath, and finally the last term in
the $H_{SB}$ describes the coupling of the remaining $(M+N+1)$th spin with
both baths. The first and second baths are characterized by spectral density
functions

$\mathcal{I}_{1}(\xi)=\sum_{k}\lambda_{k}^{2}\delta(\xi-\Omega_{k}),\quad\mathcal{I}_{2}(\xi)=\sum_{k}\kappa_{k}^{2}\delta(\xi-\Upsilon_{k}),$
(57)

the reorganization energies

$\displaystyle\Omega=\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)/\xi,\quad\Upsilon=\int_{0}^{\infty}d\xi\,\mathcal{I}_{2}(\xi)/\xi).$
(58)

We are interested in the evaluation of the average

$\displaystyle\langle\sigma^{x}\rangle=\frac{1}{Z}\text{Tr}\Big{[}e^{-\beta
H}\sigma^{x}\Big{]},$ (59)

where $Z\equiv\text{Tr}[e^{-\beta H}]$. The coherence
$\langle\sigma^{x}\rangle$ reads (details of the calculation are presented in
the Appendix D)

$\displaystyle\langle\sigma^{x}\rangle=$
$\displaystyle\frac{Z_{B}}{Z}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\text{Tr}_{S}\Big{[}\exp\Big{\\{}-\beta\frac{R_{mn}}{2}\sigma^{z}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}f_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}g_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{2}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$
$\displaystyle(\cos\theta_{mn}\sigma^{x}-\sin\theta_{mn}\sigma^{z})\Big{]}.$
(60)

Here we introduced the partition function of the bosonic system $Z_{B}$, and
the notations $R_{mn}=(\omega^{2}+\omega^{2}_{mn})^{1/2}$, with
$\omega_{mn}=4f_{1}f_{2}\Omega(M-2m)+4g_{1}g_{2}\Upsilon(N-2n)$,
$\sin\theta_{mn}=\omega_{mn}/R_{mn}$ and $\cos\theta_{mn}=\omega/R_{mn}$.
Furthermore,
$\Sigma^{x}(\tau)=\cosh(R_{mn}\tau)[\cos\theta_{mn}\sigma^{x}+i\cos\theta_{mn}\sigma^{y}]-\sin\theta_{mn}\sigma^{z}$,
and $G(\xi,\tau-\tau^{\prime})$ is defined in Eq. (81). Finally, the partition
function $Z$ of the whole system is

$\displaystyle Z=$ $\displaystyle Z_{B}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose
m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\text{Tr}_{S}\Big{[}\exp\Big{\\{}-\beta\frac{R_{mn}}{2}\sigma^{z}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}f_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}g_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{2}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\Big{]}.$
(61)

Note that the perturbative regime corresponds to the situation
$\omega\gg[\omega_{mn}]_{\text{max}}=4f_{1}f_{2}\Omega M+4g_{1}g_{2}\Upsilon
N$. In this limit $R_{mn}\approx\omega$,
$\sin\theta_{mn}\approx[4f_{1}f_{2}\Omega(M-2m)+4g_{1}g_{2}\Upsilon(N-2n)]/\omega$,
$\cos\theta_{mn}\approx 1$.

Let’s check this formula for the previously considered case of $M+1$ spins,
where $4f_{1}f_{2}\Omega M\ll\omega$, $g_{2}=0$. For this case the expression
for the statistical sum takes the form

$\displaystyle Z=$ $\displaystyle Z_{B}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose
m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\text{Tr}_{S}\Big{[}\exp\Big{\\{}-\beta\frac{R_{m}}{2}\sigma^{z}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}f_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\Big{]},$
(62)

where we introduced $R_{m}=\sqrt{\omega^{2}+[4f_{1}f_{2}\Omega(M-2m)]^{2}}$.
We keep the leading term in $T_{\tau}$ exponent, i.e., replace it by $1$. Then
we use the approximation $R_{m}\approx\omega$, and remove the small terms
$\propto f_{1}^{2}\Omega,g_{1}^{2}\Upsilon$ in the expression for the
exponent. Then we get

$\displaystyle Z\approx$ $\displaystyle
Z_{B}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\text{Tr}_{S}\Big{[}\exp\Big{\\{}-\beta\frac{\omega}{2}\sigma^{z}\Big{\\}}\Big{]}$
$\displaystyle=$ $\displaystyle
Z_{B}\Big{[}2\cosh\Big{(}\frac{\beta\omega}{2}\Big{)}\Big{]}\Big{[}2\cosh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}\Big{]}^{M}\Big{[}2\cosh\Big{(}\frac{\beta\omega_{2}}{2}\Big{)}\Big{]}^{N},$
(63)

which is nothing but the statistical sum of the non-interacting spin and boson
subsystems. Using the same approximation and
$\sin\theta_{m}=4f_{1}f_{2}\Omega(M-2m)/R_{m}\approx
4f_{1}f_{2}\Omega(M-2m)/\omega$ we obtain the following expression for the
coherence in the system

$\displaystyle\langle\sigma^{x}\rangle\approx$
$\displaystyle\frac{Z_{B}}{Z}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$
$\displaystyle\text{Tr}_{S}\Big{[}\exp\Big{\\{}-\beta\frac{\omega}{2}\sigma^{z}\Big{\\}}(-\sin\theta_{m}\sigma^{z})\Big{]}=\frac{Z_{B}}{Z}\frac{4f_{1}f_{2}\Omega}{\omega}\Big{[}2\sinh\Big{(}\frac{\beta\omega}{2}\Big{)}\Big{]}\Big{(}-\frac{2}{\beta}\Big{)}\times$
$\displaystyle\times$
$\displaystyle\frac{\partial}{\partial\omega_{1}}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose
m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}=$
$\displaystyle=$
$\displaystyle-\frac{4f_{1}f_{2}\Omega}{\omega}\tanh\Big{(}\frac{\beta\omega}{2}\Big{)}M\tanh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)},$
(64)

which coincides with the previously obtained result (32). Considering the
general case with $g_{2}\neq 0$ in the same approximation we get

$\displaystyle\langle\sigma^{x}\rangle\approx-4\frac{\tanh\Big{(}\frac{\beta\omega}{2}\Big{)}}{\omega}\Big{[}f_{1}f_{2}\Omega
M\tanh\Big{(}\frac{\beta\omega_{1}}{2}\Big{)}+g_{1}g_{2}\Upsilon
N\tanh\Big{(}\frac{\beta\omega_{2}}{2}\Big{)}\Big{]}.$ (65)

Therefore, the impact of several baths on the coherence of the $\sigma$’s spin
has an additive character. It means that independent baths can constructively
join to generate coherence of the output spins.

Note the several advantages of the results (6), (6), and (D): i) it can be
easily generalized for the arbitrary number of baths; ii) it is non-
perturbative in the coupling parameters $f_{,}f_{2},g_{1},g_{2}$; iii) it
allows to obtain the answers beyond the weak coupling limit $4f_{1}f_{2}\Omega
M+4g_{1}g_{2}\Upsilon N\ll\omega$, $2f_{1}^{2}\Omega M\ll\omega_{1}$,
$2g_{1}g_{2}\Upsilon N\ll\omega_{2}$. However, the evaluation of the average
of the coherence beyond the weak coupling limit requires either more
sophisticated analytical methods or direct numerical calculations.

##  7 Generalization for the oscillator coherences

To extend the experimental possibilities we finally replace the output spin
system of the previous case by an oscillator. The simplest model Hamiltonian
of the corresponding system analogously to Sec. 2 is
$H=H_{A}+H_{S}+H_{B}+H_{int}$, where

$\displaystyle H_{A}=E\Big{(}a^{\dagger}a+\frac{1}{2}\Big{)},$ (66)

is the oscillator Hamiltonian,

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

the Hamiltoinan of the bosonic excitations of the bath,

$\displaystyle H_{S}=\frac{\omega}{2}\sigma^{z},$ (68)

is the spin Hamiltonian and finally the term

$\displaystyle
H_{int}=\Big{[}f\sigma^{z}+g(a^{\dagger}+a)\Big{]}\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k}),$
(69)

describes an interaction of the spin and the oscillator with the bosonic bath.

Applying the general formula for this case we obtain

$\rho_{A}=\frac{e^{-\beta
H_{A}}}{Z_{A}}\frac{\Big{\langle}T_{\tau}\exp\Big{\\{}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})F(\tau)F(\tau^{\prime})\Big{\\}}\Big{\rangle}_{S}}{\Big{\langle}\Big{\langle}T_{\tau}\exp\Big{\\{}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})F(\tau)F(\tau^{\prime})\Big{\\}}\Big{\rangle}_{S}\Big{\rangle}_{A}},$
(70)

where we have introduced

$\displaystyle F(\tau)=f\sigma^{z}+g(a^{\dagger}e^{\tau E}+ae^{-\tau E})$ (71)

The reduced density matrix of the oscillator can be calculated as a power
series in the parameters $f,g$. We calculate the leading order of the
$T_{\tau}$-exponent in the numerator of Eq. (70)

$\displaystyle\Big{\langle}T_{\tau}\exp\\{\dots\\}\Big{\rangle}_{S}\approx$
$\displaystyle
1+\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{\\{}\frac{f^{2}\beta}{\xi}-\frac{2fg}{E\xi}\tanh\Big{(}\frac{\beta\omega}{2}\Big{)}\Big{[}a^{\dagger}(e^{\beta
E}-1)+a(1-e^{-\beta E})\Big{]}+$ $\displaystyle+$ $\displaystyle
g^{2}\Big{[}\frac{E\coth\left(\frac{\beta\xi}{2}\right)-\xi\coth\left(\frac{\beta
E}{2}\right)}{2E\left(E^{2}-\xi^{2}\right)}\Big{\\{}(a^{\dagger})^{2}\left(e^{\beta
E}-1\right)^{2}+a^{2}\left(1-e^{-\beta E}\right)^{2}\Big{\\}}+$
$\displaystyle+$
$\displaystyle\Big{(}\frac{\beta(\xi-E)+e^{\beta(E-\xi)}-1}{\left(1-e^{-\beta\xi}\right)(E-\xi)^{2}}+\frac{-\beta(\xi+E)+e^{\beta(\xi+E)}-1}{\left(e^{\beta\xi}-1\right)(\xi+E)^{2}}\Big{)}a^{\dagger}a+$
$\displaystyle+$
$\displaystyle\Big{(}\frac{\beta(E-\xi)+e^{\beta(\xi-E)}-1}{\left(e^{\beta\xi}-1\right)(E-\xi)^{2}}+\frac{\beta(\xi+E)+e^{-\beta(\xi+E)}-1}{\left(1-e^{-\beta\xi}\right)(\xi+E)^{2}}\Big{)}aa^{\dagger}\Big{]}\Big{\\}}.$
(72)

The denominator of the reduced density matrix then reads

$\displaystyle\Big{\langle}T_{\tau}\exp\\{\dots\\}\Big{\rangle}_{A,S}\approx
1+\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\Big{\\{}\frac{f^{2}\beta}{\xi}+g^{2}\frac{\beta\left[E\coth\left(\frac{\beta\xi}{2}\right)-\xi\coth\left(\frac{\beta
E}{2}\right)\right]}{E^{2}-\xi^{2}}\Big{\\}}$ (73)

Let us use the obtained result for $\rho_{A}$ and calculate the average values
for the dimensionless coordinate operator $x\equiv(a^{\dagger}+a)/\sqrt{2}$

$\langle
x\rangle\approx-\frac{4fg}{\sqrt{2}E}\Omega\tanh\Big{(}\frac{\beta\omega}{2}\Big{)}.$
(74)

Note that the non-zero value of $x$ operator exists if both values $g$ and $f$
are non-zero. This result can be understood in terms of a two-step model.
First, the spin polarizes the bosonic bath via $f\sigma^{z}$ term and produces
the non-zero value $\langle b_{k}+b_{k}^{\dagger}\rangle\propto f$ in the
leading order of the perturbation theory. Then the polarized bath generates
the non-zero coordinate shift $\langle a^{\dagger}+a\rangle\propto g\langle
b_{k}+b_{k}^{\dagger}\rangle\propto gf$ via the $g(a^{\dagger}+a)$ term. The
logic of the method advantageously follows the previous case of the output
spin. The square of the coordinate operator takes the form

$\displaystyle\langle x^{2}\rangle\approx$
$\displaystyle\frac{1}{2}\coth\Big{(}\frac{\beta
E}{2}\Big{)}+2Eg^{2}\int_{0}^{\infty}\\!\\!d\xi\mathcal{I}(\xi)\frac{E\coth\left(\frac{\beta\xi}{2}\right)-\xi\coth\left(\frac{\beta
E}{2}\right)}{\left(E^{2}-\xi^{2}\right)^{2}}-$ $\displaystyle-$
$\displaystyle\frac{g^{2}}{2E}\frac{\beta E+\sinh(\beta E)}{\sinh^{2}(\beta
E/2)}\int_{0}^{\infty}\\!\\!d\xi\frac{\mathcal{I}(\xi)\xi}{\left(E^{2}-\xi^{2}\right)}.$
(75)

As a result for the variance $\sigma^{2}\equiv\langle(x-\langle
x\rangle)^{2}\rangle$ we have $\sigma^{2}\approx\langle x^{2}\rangle$ in the
leading order. Note that the spin subsystem doesn’t have an impact onto
$\sigma^{2}$ in the leading order of the perturbation theory in coupling
parameters $f,g$.

Using the latter results one can estimate the signal-to-noise ratio $|\langle
x\rangle|/\sigma$, which, in particular, has a simple form in the low-
temperature limit $\beta(=1/T)\rightarrow\infty$

$\displaystyle\left.\frac{|\langle x\rangle|}{\sigma}\right|_{T\rightarrow
0}\approx 4fg\frac{\Omega}{E}.$ (76)

Hence, the signal-to-noise ratio in such a system can be increased by
increaing the ratio of the reorganization energy of the bosonic bath $\Omega$
to the oscillator’s energy $E$.

##  8 Summary

To conclude, using a new, more efficient method then in previous case, we
proved that autonomous coherences could be synthesised from many independent
driving spins through separate low-temperature baths and multiplxed to many
output spins. It is a crucial step to allow their experiment at investigation
for a much more diverse class of experimental platforms and final verification
of such intriguing phenomena. The driving coupling is present as pure
dephasing, for example, in quantum dots. Our analysis shows that another
output spin system embedded in the same environment can exhibit such
autonomous coherence in principle, although it cannot generate it. Moreover,
more driving spins will be advantageous, and the same bath can supply many
spins with the same autonomous coherence. Additionally, it is not required
that driving spins have to generate coherence through the same bath; many
baths coupled to the output spin will be equally good.

Our approach qualitatively overcomes the first theoretical proposal suggesting
the existence of these new autonomous spin coherences in Refs. [32, 36, 38].
The groundbreaking idea there is limited by a double role of spin in this
method and, therefore, unavoidable disruptive back-action effects. We removed
that limitation and opened a road towards further much broader investigations.
This new approach and extension are critical to observing autonomous quantum
coherences experimentally and exploiting them in diverse quantum technology
applications.

##  9 Acknowledgments

R.F. acknowledges the grant of 22-27431S of Czech Science Foundation. A.S.
acknowledges the grant LTAUSA19099 of Czech Ministry of Education, Youth and
Sport and the grant by Czech Science Foundation (project GA23-06369S).

## References

  * [1] Alexander Streltsov, Gerardo Adesso, and Martin B. Plenio. Colloquium: Quantum coherence as a resource. Rev. Mod. Phys. 89, 041003 (2017).
  * [2] Leonardo Novo, Masoud Mohseni & Yasser Omar. Disorder-assisted quantum transport in suboptimal decoherence regimes. Scientific Reports 6, 18142 (2016).
  * [3] Jose Joaquin Alonso, Eric Lutz, and Alessandro Romito. Thermodynamics of Weakly Measured Quantum Systems. Phys. Rev. Lett. 116, 080403 (2016).
  * [4] Kok Chuan Tan, Tyler Volkoff, Hyukjoon Kwon, and Hyunseok Jeong. Quantifying the Coherence between Coherent States. Phys. Rev. Lett. 119, 190405 (2017).
  * [5] Rafał Demkowicz-Dobrzański, Jan Czajkowski, and Pavel Sekatski. Adaptive Quantum Metrology under General Markovian Noise. Phys. Rev. X 7, 041009 (2017).
  * [6] 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).
  * [7] Jan Czajkowski, Krzysztof Pawłowski and Rafał Demkowicz-Dobrzański. Many-body effects in quantum metrology. New J. Phys. 21 053031 (2019).
  * [8] Stella Seah, Stefan Nimmrichter, Daniel Grimmer, Jader P. Santos, Valerio Scarani, and Gabriel T. Landi. Collisional Quantum Thermometry. Phys. Rev. Lett. 123, 180602 (2019).
  * [9] 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).
  * [10] 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).
  * [11] Flavio Del Santo and Borivoje Dakić. Coherence Equality and Communication in a Quantum Superposition. Phys. Rev. Lett. 124, 190501 (2020).
  * [12] Kaonan Micadei, Gabriel T. Landi, and Eric Lutz. Quantum Fluctuation Theorems beyond Two-Point Measurements. Phys. Rev. Lett. 124, 090602 (2020).
  * [13] María García Díaz, Giacomo Guarnieri, Mauro Paternostro. Quantum work statistics with initial coherence. Entropy 22(11), 1223 (2020).
  * [14] Stella Seah, Stefan Nimmrichter, and Valerio Scarani. Maxwell’s Lesser Demon: A Quantum Engine Driven by Pointer Measurements. Phys. Rev. Lett. 124, 100603 (2020).
  * [15] 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).
  * [16] 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).
  * [17] 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)
  * [18] 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).
  * [19] 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).
  * [20] 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).
  * [21] 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).
  * [22] 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).
  * [23] Robert Drost, Teemu Ojanen, Ari Harju & Peter Liljeroth. & Liljeroth, P. Topological states in engineered atomic lattices. Nature Physics 13, 668–671 (2017).
  * [24] 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).
  * [25] 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).
  * [26] Elisabet Romero, Vladimir I. Novoderezhkin & Rienk van Grondelle. Quantum design of photosynthesis for bio-inspired solar-energy conversion. Nature 543, 647–656 (2017).
  * [27] 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).
  * [28] 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).
  * [29] Camille L. Latune, Ilya Sinayskiy & Francesco Petruccione. Roles of quantum coherences in thermal machines. Eur. Phys. J. Spec. Top. (2021).
  * [30] 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).
  * [31] 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).
  * [32] Giacomo Guarnieri, Michal Kolář, and Radim Filip. Steady-State Coherences by Composite System-Bath Interactions. Phys. Rev. Lett. 121, 070401 (2018).
  * [33] Jonatan Bohr Brask, Géraldine Haack, Nicolas Brunner, and Marcus Huber. Autonomous quantum thermal machine for generating steady-state entanglement. New J. Phys. 17, 113029 (2015).
  * [34] 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).
  * [35] 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).
  * [36] 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).
  * [37] 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).
  * [38] Artur Slobodeniuk, Tomáš Novotný, and Radim Filip. Extraction of autonomous quantum coherences. Quantum 6, 689 (2022).
  * [39] M. Łobejko, M. Winczewski, G. Suárez, R. Alicki, and M. Horodecki. Towards reconciliation of completely positive open system dynamics with the equilibration postulate. arXiv:2204.00643v1
  * [40] Ulrich Weiss. Quantum Dissipative Systems (World Scientific Publishing, 2012)
  * [41] A. J. Leggett, S. Chakravarty, A. T. Dorsey, Matthew P. A. Fisher, Anupam Garg, and W. Zwerger. Dynamics of the dissipative two-state system. Rev. Mod. Phys. 59, 1 (1987).

##  Appendix A Evaluation of the $T_{\tau}$-exponent

To evaluate the reduced density operator in Eq. (9) we need to perform the
averaging over bosonic degrees of freedom of the $T_{\tau}$ ordered exponent
(2). To do it we present the spin-boson interaction term (3) as

$H_{SB}\equiv\Big{(}\sum_{j=1}^{M+1}\boldsymbol{\sigma}_{j}\cdot\boldsymbol{n}_{j}\Big{)}\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})\equiv
F\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})\equiv\sum_{k}(f_{k}b_{k}^{\dagger}+f_{k}b_{k}),$
(77)

and use the commutativity of the operators under the $T_{\tau}$ ordering
operation we present Eq. (2) in the form

$\displaystyle
T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}=$
$\displaystyle
T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\sum_{k}[f_{k}(\tau)b^{\dagger}_{k}(\tau)+f_{k}(\tau)b_{k}(\tau)]}\Big{\\}}=$
$\displaystyle=$ $\displaystyle
T_{\tau}\prod_{k}\Big{\\{}e^{-\int_{0}^{\beta}d\tau
f_{k}(\tau)b^{\dagger}_{k}(\tau)-\int_{0}^{\beta}d\tau^{\prime}f_{k}(\tau^{\prime})b_{k}(\tau^{\prime})}\Big{\\}}=$
$\displaystyle=$ $\displaystyle
T_{\tau}\prod_{k}\Big{[}\Big{\\{}e^{-\int_{0}^{\beta}d\tau
f_{k}(\tau)b^{\dagger}_{k}(\tau)}\Big{\\}}\Big{\\{}e^{-\int_{0}^{\beta}d\tau^{\prime}f_{k}(\tau^{\prime})b_{k}(\tau^{\prime})}\Big{\\}}\Big{]}.$
(78)

Then we the averaging over bosonic degrees of freedom of Eq. (2) takes the
form

$\displaystyle\Big{\langle}T_{\tau}$
$\displaystyle\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}\Big{\rangle}_{B}=\prod_{k}\sum_{n,m=0}^{\infty}\frac{(-1)^{n+m}}{n!m!}\int_{0}^{\beta}\dots\int_{0}^{\beta}d\tau_{1}\dots
d\tau_{n}\int_{0}^{\beta}\dots\int_{0}^{\beta}d\tau^{\prime}_{l}\dots
d\tau^{\prime}_{m}\times$ $\displaystyle\times$ $\displaystyle
T_{\tau}\Big{\\{}f_{k}(\tau_{1})\dots
f_{k}(\tau_{n})f_{k}(\tau^{\prime}_{1})\dots
f_{k}(\tau^{\prime}_{m})\Big{\\}}\Big{\langle}T_{\tau}\Big{\\{}b_{k}^{\dagger}(\tau_{1})\dots
b_{k}^{\dagger}(\tau_{n})b_{k}(\tau^{\prime}_{1})\dots
b_{k}(\tau^{\prime}_{m})\Big{\\}}\Big{\rangle}_{B}=$ $\displaystyle=$
$\displaystyle\prod_{k}\sum_{n=0}^{\infty}\frac{1}{(n!)^{2}}\prod_{j=1}^{n}\int_{0}^{\beta}\int_{0}^{\beta}d\tau_{j}d\tau^{\prime}_{j}\,T_{\tau}\Big{\\{}f_{k}(\tau_{1})\dots
f_{k}(\tau_{n})f_{k}(\tau^{\prime}_{1})\dots
f_{k}(\tau^{\prime}_{n})\Big{\\}}\times$ $\displaystyle\times$
$\displaystyle\sum_{P}\Big{\langle}T_{\tau}b_{k}^{\dagger}(\tau_{1})b_{k}(\tau^{\prime}_{j_{1}})\Big{\rangle}_{B}\dots\Big{\langle}T_{\tau}b_{k}^{\dagger}(\tau_{n})b_{k}(\tau^{\prime}_{j_{n}})\Big{\rangle}_{B}.$
(79)

The sum $\sum_{P}$ represents the summation over all the transpositions
$P(1,\dots,n)=j_{1},\dots,j_{n}$. Replacing the integration indices
$\tau_{j_{1}},\dots,\tau_{j_{n}}\rightarrow\tau_{1},\dots,\tau_{n}$ for each
separate term (the total number of such terms is $n!$) in the expression we
obtain

$\displaystyle\Big{\langle}T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}\Big{\rangle}_{B}=$
$\displaystyle\prod_{k}\sum_{n=0}^{\infty}\frac{1}{n!}\prod_{j=1}^{n}\int_{0}^{\beta}\int_{0}^{\beta}d\tau_{j}d\tau^{\prime}_{j}\Big{\langle}T_{\tau}b_{k}^{\dagger}(\tau_{j})b_{k}(\tau^{\prime}_{j})\Big{\rangle}_{B}\times$
$\displaystyle\times$ $\displaystyle T_{\tau}\Big{\\{}f_{k}(\tau_{1})\dots
f_{k}(\tau_{n})f_{k}(\tau^{\prime}_{1})\dots
f_{k}(\tau^{\prime}_{m})\Big{\\}}=$ $\displaystyle=$ $\displaystyle
T_{\tau}\Big{\\{}\exp\Big{[}\sum_{k}\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\Omega_{k},\tau-\tau^{\prime})f_{k}(\tau)f_{k}(\tau^{\prime})\Big{]}\Big{\\}}.$
(80)

Here we introduced a thermal boson Green’s function $\langle
T_{\tau}b_{k}^{\dagger}(\tau)b_{k}(\tau^{\prime})\rangle_{B}=G(\Omega_{k},\tau-\tau^{\prime})$,
where

$G(\xi,\tau-\tau^{\prime})=\frac{e^{(\tau-\tau^{\prime})\xi}}{e^{\beta\xi}-1}\theta(\tau-\tau^{\prime})+\frac{e^{(\tau-\tau^{\prime})\xi}}{1-e^{-\beta\xi}}\theta(\tau^{\prime}-\tau).$
(81)

Introducing the spectral density function
$\mathcal{I}(\xi)=\sum_{k}\lambda_{k}^{2}\delta(\xi-\Omega_{k})$ one can write

$\Big{\langle}T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}\Big{\rangle}_{B}=T_{\tau}\Big{\\{}\exp\Big{[}\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})F(\tau)F(\tau^{\prime})\Big{]}\Big{\\}},$
(82)

with
$f_{k}(\tau)=\lambda_{k}F(\tau)=\lambda_{k}\sum_{j=1}^{M+1}\boldsymbol{\sigma}_{j}(\tau)\cdot\boldsymbol{n}_{j}$,
see Eq. (17) for details. Therefore the leading term of the aforementioned
expression is

$\displaystyle\Big{\langle}T_{\tau}\Big{\\{}e^{-\int_{0}^{\beta}d\tau\widetilde{H}_{SB}(\tau)}\Big{\\}}\Big{\rangle}_{B}\approx$
$\displaystyle
1+\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\,T_{\tau}\Big{\\{}\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})F(\tau)F(\tau^{\prime})\Big{\\}}=$
$\displaystyle=$ $\displaystyle
1+\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\xi,\tau-\tau^{\prime})F(\tau)F(\tau^{\prime}),$
(83)

where we introduced the notation

$\phi(\xi,\tau-\tau^{\prime})=\frac{e^{(\tau-\tau^{\prime})\xi}}{e^{\beta\xi}-1}+\frac{e^{-(\tau-\tau^{\prime})\xi}}{1-e^{-\beta\xi}}.$
(84)

Remembering that the term $F(\tau)$ is a linear combination of the spin
operators, we conclude that the leading term of the decomposition
$F(\tau)F(\tau^{\prime})$ describes the spin-spin correlation between the spin
degrees of freedom in spin subsystem induced by its interaction with bosonic
bath.

##  Appendix B Average value of $\sigma_{M+1}^{z}$ operator for the case of
$M+1$ spins

The expression for the average value of the $\sigma_{M+1}^{z}$ operator up to
the second order of the perturbation theory in parameters $f_{1},f_{2}$ has
the following form

$\displaystyle\langle\sigma_{M+1}^{z}\rangle\approx$
$\displaystyle\langle\sigma_{M+1}^{z}\rangle_{S}\Big{(}1-\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\xi,\tau-\tau^{\prime})\langle
F(\tau)F(\tau^{\prime})\rangle_{S}\Big{)}+$ $\displaystyle+$
$\displaystyle\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\xi,\tau-\tau^{\prime})\langle
F(\tau)F(\tau^{\prime})\sigma_{M+1}^{z}\rangle_{S},$ (85)

where $\phi(\xi,\tau-\tau^{\prime})$ is defined by Eq. (84) and
$F(\tau)=\sum_{j=1}^{M+1}\boldsymbol{\sigma}_{j}(\tau)\cdot\boldsymbol{n}_{j}$,
see Eq. (17). Evaluating the integrands

$\displaystyle\langle
F(\tau)F(\tau^{\prime})\rangle_{S}\langle\sigma_{M+1}^{z}\rangle_{S}=$
$\displaystyle\Big{[}f_{1}^{2}\sum_{j,l=1}^{M}\langle\sigma_{j}^{z}\sigma_{l}^{z}\rangle_{S}+f_{2}^{2}\cosh(\omega_{M+1}(\tau-\tau^{\prime}))\Big{]}\langle\sigma_{M+1}^{z}\rangle_{S}+$
$\displaystyle+f_{2}^{2}\sinh(\omega_{M+1}(\tau-\tau^{\prime}))\langle\sigma_{M+1}^{z}\rangle_{S}^{2},$
(86) $\displaystyle\langle
F(\tau)F(\tau^{\prime})\sigma_{M+1}^{z}\rangle_{S}=$ $\displaystyle
f_{1}^{2}\sum_{j,l=1}^{M}\langle\sigma_{j}^{z}\sigma_{l}^{z}\rangle_{S}\langle\sigma_{M+1}^{z}\rangle_{S}+f_{2}^{2}\cosh(\omega_{M+1}(\tau-\tau^{\prime}))\langle\sigma_{M+1}^{z}\rangle_{S}+$
$\displaystyle+f_{2}^{2}\sinh(\omega_{M+1}(\tau-\tau^{\prime})),$ (87)

and substituting them into the Eq. (B) we get

$\displaystyle\langle\sigma_{M+1}^{z}\rangle\approx$
$\displaystyle\langle\sigma_{M+1}^{z}\rangle_{S}+f_{2}^{2}\Big{[}1-\langle\sigma_{M+1}^{z}\rangle_{S}^{2}\Big{]}\times$
$\displaystyle\times\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\xi,\tau-\tau^{\prime})\sinh(\omega_{M+1}(\tau-\tau^{\prime})),$
(88)

where $\langle\sigma_{M+1}^{z}\rangle_{S}=-\tanh(\beta\omega_{M+1}/2)$ is the
average value of the $\sigma_{M+1}^{z}$ spin operator in the absence of the
spin-bath coupling. Therefore the value
$\delta\langle\sigma_{M+1}^{z}\rangle\equiv\langle\sigma_{M+1}^{z}\rangle-\langle\sigma_{M+1}^{z}\rangle_{S}$,
describes the correction to the $z$ component of the $(M+1)$th spin induced by
the spin-bath interaction. The final result reads

$\displaystyle\delta\langle$
$\displaystyle\sigma_{M+1}^{z}\rangle=f_{2}^{2}\int_{0}^{\infty}d\xi\frac{\mathcal{I}(\xi)}{\left(\xi^{2}-\omega_{M+1}^{2}\right)^{2}}\times$
$\displaystyle\times\Big{\\{}\frac{\coth\left(\frac{\beta\xi}{2}\right)}{\cosh^{2}\Big{(}\frac{\beta\omega_{M+1}}{2}\Big{)}}\left[\beta\omega_{M+1}\left(\xi^{2}-\omega_{M+1}^{2}\right)+\left(\xi^{2}+\omega_{M+1}^{2}\right)\sinh(\beta\omega_{M+1})\right]-4\xi\omega_{M+1}\Big{\\}}.$
(89)

##  Appendix C Details of the calculation of $\sigma^{2}$, Eq. (42)

In the leading order of the perturbation theory the Eq. (42) reads
$\sigma^{2}\approx\langle\sigma_{M+1}^{x}\sigma_{M+2}^{x}\rangle$. Then one
obtains

$\displaystyle\sigma^{2}\approx$
$\displaystyle\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\tau}d\tau^{\prime}\phi(\xi,\tau-\tau^{\prime})\langle
F(\tau)F(\tau^{\prime})\sigma_{M+1}^{x}\sigma_{M+2}^{x}\rangle_{S}=$
$\displaystyle=$ $\displaystyle
f_{2}^{2}\int_{0}^{\infty}d\xi\,\frac{\mathcal{I}(\xi)\mathcal{K}(\xi)}{\left(\xi^{2}-\omega_{M+1}^{2}\right)\left(\xi^{2}-\omega_{M+2}^{2}\right)},$
(90)

where $\phi(\xi,\tau-\tau^{\prime})$ is defined by Eq. (84),
$F(\tau)=\sum_{j=1}^{M+2}\boldsymbol{\sigma}_{j}(\tau)\cdot\boldsymbol{n}_{j}$.
The function $\mathcal{K}(\xi)$ reads

$\displaystyle\mathcal{K}(\xi)=$ $\displaystyle
a(\xi)+\tanh\Big{(}\frac{\beta\omega_{M+1}}{2}\Big{)}\tanh\Big{(}\frac{\beta\omega_{M+2}}{2}\Big{)}b(\xi)-$
$\displaystyle-$
$\displaystyle\tanh\Big{(}\frac{\beta\omega_{M+2}}{2}\Big{)}c(\xi)-\tanh\Big{(}\frac{\beta\omega_{M+1}}{2}\Big{)}d(\xi),$
(91)

where we introduced

$\displaystyle a(\xi)=$
$\displaystyle\xi\omega_{M+1}\sinh(\beta\omega_{M+1})\left(\cosh(\beta\omega_{M+2})\frac{\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}-1\right)-$
$\displaystyle-$
$\displaystyle\xi\omega_{M+2}\sinh(\beta\omega_{M+2})\left(\cosh(\beta\omega_{M+1})\frac{\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}+1\right)+$
$\displaystyle+$
$\displaystyle\omega_{M+1}\omega_{M+2}\coth\left(\frac{\beta\xi}{2}\right)\sinh(\beta\omega_{M+1})\sinh(\beta\omega_{M+2})-$
$\displaystyle-$ $\displaystyle
4\xi^{2}\coth\left(\frac{\beta\xi}{2}\right)\sinh^{2}\left(\frac{\beta\omega_{M+1}}{2}\right)\sinh^{2}\left(\frac{\beta\omega_{M+2}}{2}\right),$
(92) $\displaystyle b(\xi)=$
$\displaystyle\xi\omega_{M+2}\sinh(\beta\omega_{M+1})\left(1-\frac{\cosh(\beta\omega_{M+2})\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}\right)+$
$\displaystyle+$
$\displaystyle\xi\omega_{M+1}\sinh(\beta\omega_{M+2})\left(\frac{\cosh(\beta\omega_{M+1})\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}+1\right)+$
$\displaystyle+$ $\displaystyle
4\omega_{M+1}\omega_{M+2}\coth\left(\frac{\beta\xi}{2}\right)\sinh^{2}\left(\frac{\beta\omega_{M+1}}{2}\right)\sinh^{2}\left(\frac{\beta\omega_{M+2}}{2}\right)-$
$\displaystyle-$
$\displaystyle\xi^{2}\coth\left(\frac{\beta\xi}{2}\right)\sinh(\beta\omega_{M+1})\sinh(\beta\omega_{M+2}).$
(93) $\displaystyle c(\xi)=$
$\displaystyle\xi\omega_{M+1}\sinh(\beta\omega_{M+1})\sinh(\beta\omega_{M+2})\frac{\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}-$
$\displaystyle-$
$\displaystyle\xi\omega_{M+2}\Big{[}\cosh(\beta\omega_{M+1})\cosh(\beta\omega_{M+2})-1\Big{]}\frac{\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}+$
$\displaystyle+$ $\displaystyle
2\omega_{M+1}\omega_{M+2}\coth\left(\frac{\beta\xi}{2}\right)\sinh(\beta\omega_{M+1})\sinh^{2}\left(\frac{\beta\omega_{M+2}}{2}\right)-$
$\displaystyle-$ $\displaystyle
2\xi^{2}\coth\left(\frac{\beta\xi}{2}\right)\sinh^{2}\left(\frac{\beta\omega_{M+1}}{2}\right)\sinh(\beta\omega_{M+2})+$
$\displaystyle+$
$\displaystyle\xi\omega_{M+2}\Big{[}\cosh(\beta\omega_{M+1})-\cosh(\beta\omega_{M+2})\Big{]},$
(94) $\displaystyle d(\xi)=$
$\displaystyle-\xi\omega_{M+2}\sinh(\beta\omega_{M+1})\sinh(\beta\omega_{M+2})\frac{\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}+$
$\displaystyle+$
$\displaystyle\xi\omega_{M+1}\Big{[}\cosh(\beta\omega_{M+1})\cosh(\beta\omega_{M+2})-1\Big{]}\frac{\left(2\xi^{2}-\omega_{M+1}^{2}-\omega_{M+2}^{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}+$
$\displaystyle+$ $\displaystyle
2\omega_{M+1}\omega_{M+2}\coth\left(\frac{\beta\xi}{2}\right)\sinh^{2}\left(\frac{\beta\omega_{M+1}}{2}\right)\sinh(\beta\omega_{M+2})-$
$\displaystyle-$ $\displaystyle
2\xi^{2}\coth\left(\frac{\beta\xi}{2}\right)\sinh(\beta\omega_{M+1})\sinh^{2}\left(\frac{\beta\omega_{M+2}}{2}\right)+$
$\displaystyle+$
$\displaystyle\xi\omega_{M+1}\Big{[}\cosh(\beta\omega_{M+2})-\cosh(\beta\omega_{M+1})\Big{]}.$
(95)

After simplification one gets the following expression

$\displaystyle\mathcal{K}(\xi)=4\Big{[}$
$\displaystyle\omega_{M+1}\omega_{M+2}\coth\left(\frac{\beta\xi}{2}\right)\tanh\left(\frac{\beta\omega_{M+1}}{2}\right)\tanh\left(\frac{\beta\omega_{M+2}}{2}\right)+$
$\displaystyle+\frac{\xi\omega_{M+1}(\xi^{2}-\omega_{M+1}^{2})\tanh\left(\frac{\beta\omega_{M+1}}{2}\right)-\xi\omega_{M+2}(\xi^{2}-\omega_{M+2}^{2})\tanh\left(\frac{\beta\omega_{M+2}}{2}\right)}{\omega_{M+1}^{2}-\omega_{M+2}^{2}}\Big{]},$
(96)

which after substituting it in the expression (C) finally leads to Eqs. (43)
and (44) for $\sigma^{2}$.

##  Appendix D Average value of $\sigma^{x}$ for the case of $M+N+1$ spins and
two independent baths

We calculate the average $\langle\sigma^{x}\rangle$ using the formula

$\displaystyle\langle\sigma^{x}\rangle=$
$\displaystyle\frac{1}{Z}\text{Tr}\Big{[}e^{-\beta H}\sigma^{x}\Big{]},$ (97)

with the Hamiltonian (53). We rewrite this expression in the form

$\displaystyle\langle\sigma^{x}\rangle=$
$\displaystyle\frac{1}{Z}\text{Tr}\Big{[}e^{-\beta
H}\sigma^{x}\Big{]}=\frac{1}{Z}\text{Tr}\Big{[}e^{-\beta
U^{\dagger}HU}U^{\dagger}\sigma^{x}U\Big{]},$ (98)

where $U$ is a unitary operator and $Z=\text{Tr}[e^{-\beta
H}]=\text{Tr}[e^{-\beta U^{\dagger}HU}]$. We choose the unitary operator in
the form of the Lang-Firsov type

$\displaystyle U=$
$\displaystyle\exp\Big{[}-f_{1}\sum_{j=1}^{M}\sum_{k}\frac{\lambda_{k}}{\Omega_{k}}(b_{k}^{\dagger}-b_{k})\sigma_{j}^{z}\Big{]}\times\exp\Big{[}-g_{1}\sum_{j=M+1}^{M+N}\sum_{k}\frac{\kappa_{k}}{\Upsilon_{k}}(d_{k}^{\dagger}-d_{k})\sigma_{j}^{z}\Big{]}.$
(99)

Therefore $U^{\dagger}\sigma^{x}U=\sigma^{x}$ and the transformed Hamiltonian
$\mathcal{H}=U^{\dagger}HU$

$\mathcal{H}=\mathcal{H}_{S1}+\mathcal{H}_{S2}+\mathcal{H}_{B}+\mathcal{H}_{S}+\mathcal{H}_{\text{int},S}+\mathcal{H}_{\text{int},S12},$
(100)

where

$\displaystyle\mathcal{H}_{S1}=\frac{\omega_{1}}{2}\sum_{j=1}^{M}\sigma_{j}^{z}-\Omega
f_{1}^{2}\sum_{j,l=1}^{M}\sigma_{j}^{z}\sigma_{l}^{z},\quad\mathcal{H}_{S2}=\frac{\omega_{2}}{2}\sum_{j=M+1}^{M+N}\sigma_{j}^{z}-\Upsilon
g_{1}^{2}\sum_{j,l=M+1}^{M+N}\sigma_{j}^{z}\sigma_{l}^{z},$ (101)

are the new spin Hamiltonians, corresponding to the group of spins coupled to
the first ($\mathcal{H}_{S1}$) and the second ($\mathcal{H}_{S2}$) bosonic
baths. Here we introduced the reorganization energies
$\Omega=\sum_{k}\lambda_{k}^{2}/\Omega_{k}$,
$\Upsilon=\sum_{k}\kappa_{k}^{2}/\Upsilon_{k}$ of the first and the second
bath, respectively. The Hamiltonian of the corresponding baths is

$\displaystyle\mathcal{H}_{B}=\sum_{k}(\Omega_{k}b_{k}^{\dagger}b_{k}+\Upsilon_{k}d_{k}^{\dagger}d_{k}).$
(102)

Then,

$\mathcal{H}_{S}=\frac{\omega}{2}\sigma^{z},$ (103)

is the Hamiltonian of the remaining (output) spin coupled to both bathes by
term

$\displaystyle\mathcal{H}_{\text{int},S}=$ $\displaystyle
f_{2}\sigma^{x}\sum_{k}\lambda_{k}(b_{k}^{\dagger}+b_{k})+g_{2}\sigma^{x}\sum_{k}\kappa_{k}(d_{k}^{\dagger}+d_{k}).$
(104)

Finally, the term

$\displaystyle\mathcal{H}_{\text{int},S12}=$
$\displaystyle-2f_{1}f_{2}\Omega\sigma^{x}\sum_{j=1}^{M}\sigma_{j}^{z}-2g_{1}g_{2}\Upsilon\sigma^{x}\sum_{j=M+1}^{M+N}\sigma_{j}^{z},$
(105)

describes effective coupling between $M+N$ input spins $\sigma_{j}^{z}$ and
the output spin $\sigma^{x}$.

Using the structure the spin Hamiltonians $\mathcal{H}_{S1}$ and
$\mathcal{H}_{S2}$ one can calculate their eigenstates and eigenvalues.
Namely, the eigenstates can be written as a product of all possible
combinations of the spin-up $|\uparrow\rangle$ and spin-down
$|\downarrow\rangle$ states. For example, for the Hamiltonian
$\mathcal{H}_{S1}$ has $2^{M}$ eigenstates

$|s_{1},s_{2},\dots
s_{M}\rangle=|s_{1}\rangle|s_{2}\rangle\dots|s_{M}\rangle,$ (106)

where $s_{j}=\uparrow,\downarrow$ for each $j=1,2,\dots M$. The eigenvalues of
the Hamiltonian don’t depend on the position of the spin-up and spin-down
states in the eigenstate in $|s_{1},s_{2},\dots s_{M}\rangle$, but only on the
total number of spin-down and spin-up in the corresponding eigenstate.
Therefore all the eigenstates with $m$ spin-down and $M-m$ spin-up states are
characterized by the eigenvalue

$E_{1,m}=\frac{\omega_{1}}{2}(M-2m)-\Omega f_{1}^{2}(M-2m)^{2},$ (107)

which represents ${M\choose m}=M!/m!(M-m)!$ times degenerated energy level.
The analogous consideration can be performed for the Hamiltonian
$\mathcal{H}_{S2}$. Note that the Hamiltonians $\mathcal{H}_{S1}$ and
$\mathcal{H}_{S2}$ commute with the rest of the terms in the $\mathcal{H}$ and
therefore, we can take a trace over the input spin states

$\displaystyle\text{Tr}_{S1}\Big{[}\text{Tr}_{S2}\Big{[}e^{-\beta\mathcal{H}}\Big{]}\Big{]}=$
$\displaystyle\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\times$
$\displaystyle\times$
$\displaystyle\exp\Big{\\{}-\beta\big{[}\mathcal{H}_{B}+\mathcal{H}_{\text{int},S}+\big{(}\mathcal{H}_{S}-\frac{\omega_{mn}}{2}\sigma^{x}\big{)}\big{]}\Big{\\}},$
(108)

where $\omega_{mn}\equiv 4f_{1}f_{2}\Omega(M-2m)+4g_{1}g_{2}\Upsilon(N-2n)$.
We transform the last exponent in the expression to the form

$\displaystyle\exp\Big{\\{}-\beta$
$\displaystyle\big{[}\mathcal{H}_{mn}+\mathcal{H}_{\text{int},S}\big{]}\Big{\\}}=\exp\Big{\\{}-\beta\mathcal{H}_{mn}\Big{\\}}T_{\tau}\Big{\\{}\exp\Big{(}-\int_{0}^{\beta}d\tau\widetilde{H}_{\text{int},S}(\tau)\Big{)}\Big{\\}}$
(109)

using the notations
$\mathcal{H}_{mn}=\mathcal{H}_{B}+\big{(}\mathcal{H}_{S}-\frac{\omega_{mn}}{2}\sigma^{x}\big{)}$
and

$\displaystyle\widetilde{H}_{\text{int},S}(\tau)=$
$\displaystyle\exp\Big{\\{}\tau\mathcal{H}_{mn}\Big{\\}}H_{\text{int},S}\exp\Big{\\{}-\tau\mathcal{H}_{mn}\Big{\\}}=$
$\displaystyle=$
$\displaystyle\sum_{k}f_{k}(\tau)[b_{k}^{\dagger}(\tau)+b_{k}(\tau)]+\sum_{k}g_{k}(\tau)[d_{k}^{\dagger}(\tau)+d_{k}(\tau)].$
(110)

Here $f_{k}(\tau)\equiv f_{2}\lambda_{k}\sigma^{x}(\tau)$, $g_{k}(\tau)\equiv
g_{2}\kappa_{k}\sigma^{x}(\tau)$, $b_{k}(\tau)=e^{-\Omega_{k}\tau}b_{k}$,
$b_{k}^{\dagger}(\tau)=e^{\Omega_{k}\tau}b_{k}^{\dagger}$,
$d_{k}(\tau)=e^{-\Upsilon_{k}\tau}d_{k}$,
$d_{k}^{\dagger}(\tau)=e^{\Upsilon_{k}\tau}d_{k}^{\dagger}$,

$\displaystyle\sigma^{x}(\tau)=$
$\displaystyle[\sin^{2}\theta_{mn}+\cos^{2}\theta_{mn}\cosh(\tau
R_{mn})]\sigma^{x}+i\cos\theta_{mn}\sinh(\tau R_{mn})\sigma^{y}+$
$\displaystyle+$ $\displaystyle\sin\theta_{mn}\cos\theta_{mn}[\cosh(\tau
R_{mn})-1]\sigma^{z},$ (111)

with $\theta_{mn}\equiv\arccos(\omega/R_{mn})$,
$R_{mn}=\sqrt{\omega^{2}+\omega_{mn}^{2}}$. Now we trace out the bosonic
degrees of freedom from the reduced density operator. To do it we introduce
the average

$\displaystyle\langle\star\rangle_{B}=\frac{1}{Z_{B}}\text{Tr}_{B}\Big{[}e^{-\beta\mathcal{H}_{B}}\star\Big{]},$
(112)

where $Z_{B}=\text{Tr}_{B}\Big{[}e^{-\beta\mathcal{H}_{B}}\Big{]}$. Then, the
$\sigma$-dependent reduced density matrix is

$\rho_{S}=\frac{Z_{B}}{Z}\Big{\langle}\text{Tr}_{S1}\Big{[}\text{Tr}_{S2}\Big{[}e^{-\beta\mathcal{H}}\Big{]}\Big{]}\Big{\rangle}_{B},$
(113)

where $Z$ is the partition function of the full system. After the averaging
procedure we obtain

$\displaystyle\rho_{S}=$
$\displaystyle\frac{Z_{B}}{Z}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\exp\Big{\\{}-\beta\big{(}\mathcal{H}_{S}-\frac{\omega_{mn}}{2}\sigma^{x}\big{)}\Big{\\}}\times$
$\displaystyle\times$
$\displaystyle\Big{\langle}T_{\tau}\Big{\\{}\exp\Big{(}-\int_{0}^{\beta}d\tau\widetilde{H}_{\text{int},S}(\tau)\Big{)}\Big{\\}}\Big{\rangle}_{B}.$
(114)

Then we calculate the average of the $T_{\tau}$-exponents using the results of
Appendix A

$\displaystyle\rho_{S}=$
$\displaystyle\frac{Z_{B}}{Z}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\exp\Big{\\{}-\beta\big{(}\mathcal{H}_{S}-\frac{\omega_{mn}}{2}\sigma^{x}\big{)}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}f_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\sigma^{x}(\tau)\sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}g_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{2}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\sigma^{x}(\tau)\sigma^{x}(\tau^{\prime})\Big{\\}},$
(115)

with the Green’s function (81) and and spectral density functions (57). In
order to reduce this expression for the most convenient form for further
calculations we use the identity

$\displaystyle\exp\Big{\\{}$
$\displaystyle-\beta\big{(}\frac{\omega}{2}\sigma^{z}-\frac{\omega_{mn}}{2}\sigma^{x}\big{)}\Big{\\}}=\exp\Big{\\{}i\frac{\theta_{mn}}{2}\sigma^{y}\Big{\\}}\exp\Big{\\{}-\beta\frac{R_{mn}}{2}\sigma^{z}\Big{\\}}\exp\Big{\\{}-i\frac{\theta_{mn}}{2}\sigma^{y}\Big{\\}}$
(116)

and we rewrite $\rho_{S}$ in the form

$\displaystyle\rho_{S}=$
$\displaystyle\exp\Big{\\{}i\frac{\theta_{mn}}{2}\sigma^{y}\Big{\\}}\Big{[}\frac{Z_{B}}{Z}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose
m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\exp\Big{\\{}-\beta\frac{R_{mn}}{2}\sigma^{z}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}f_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}g_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{2}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\Big{]}\exp\Big{\\{}-i\frac{\theta_{mn}}{2}\sigma^{y}\Big{\\}}.$
(117)

Here we introduced

$\displaystyle\Sigma^{x}(\tau)=$ $\displaystyle
e^{-i\frac{\theta_{mn}}{2}\sigma^{y}}\sigma^{x}(\tau)e^{i\frac{\theta_{mn}}{2}\sigma^{y}}=\cos\theta_{mn}[\cosh(R_{mn}\tau)\sigma^{x}+i\sinh(R_{mn}\tau)\sigma^{y}]-\sin\theta_{mn}\sigma^{z}$
$\displaystyle\equiv$ $\displaystyle\mathbf{N}(\tau)\cdot\boldsymbol{\sigma},$
(118)

which gives the following expression

$\displaystyle\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})=$
$\displaystyle\mathbf{N}(\tau)\cdot\mathbf{N}(\tau^{\prime})+i\boldsymbol{\sigma}\cdot(\mathbf{N}(\tau)\times\mathbf{N}(\tau^{\prime}))=$
$\displaystyle=$
$\displaystyle(\cos^{2}\theta_{mn}\cosh[R_{mn}(\tau-\tau^{\prime})]+\sin^{2}\theta_{mn})\mathcal{I}+\sigma^{z}\cos^{2}\theta_{mn}\sinh[R_{nm}(\tau-\tau^{\prime})]+$
$\displaystyle+$
$\displaystyle\sigma^{x}\sin\theta_{mn}\cos\theta_{mn}[\sinh(R_{nm}\tau)-\sinh(R_{nm}\tau^{\prime})]+$
$\displaystyle+$ $\displaystyle
i\sigma^{y}\sin\theta_{mn}\cos\theta_{mn}[\cosh(R_{nm}\tau)-\cosh(R_{nm}\tau^{\prime})]$
(119)

Taking into account the identity $\text{Tr}_{S}[\rho_{S}]=1$ and cyclic
properties of trace operation we obtain

$\displaystyle Z=$ $\displaystyle Z_{B}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose
m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\text{Tr}_{S}\Big{[}\exp\Big{\\{}-\beta\frac{R_{mn}}{2}\sigma^{z}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}f_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}g_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{2}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\Big{]}.$
(120)

Therefore the expression for the coherence takes the form

$\displaystyle\langle\sigma^{x}\rangle=$
$\displaystyle\frac{Z_{B}}{Z}\sum_{m=0}^{M}\sum_{n=0}^{N}{M\choose m}{N\choose
n}\exp\Big{\\{}-\beta\big{[}\frac{\omega_{1}}{2}(M-2m)+\frac{\omega_{2}}{2}(N-2n)\big{]}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle\exp\Big{\\{}\beta\big{[}\Omega
f_{1}^{2}(M-2m)^{2}+\Upsilon
g_{1}^{2}(N-2n)^{2}\big{]}\Big{\\}}\text{Tr}_{S}\Big{[}\exp\Big{\\{}-\beta\frac{R_{mn}}{2}\sigma^{z}\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}f_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{1}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\exp\Big{\\{}g_{2}^{2}\int_{0}^{\infty}d\xi\,\mathcal{I}_{2}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}\times$
$\displaystyle\times$
$\displaystyle(\cos\theta_{mn}\sigma^{x}-\sin\theta_{mn}\sigma^{z})\Big{]}.$
(121)

Note that this is the exact result which takes into account non-perturbatively
the coupling constants $f_{1},f_{2},g_{1},g_{2}$ of the Hamiltonian. Moreover,
the perturbative calculation of $T_{\tau}$ exponents in this expression
requires $f_{2},g_{2}\ll 1$, which is much better assumption than one used in
Secs. 4 and 5, namely $f_{2}\ll 1/M$ with the large number $M$.

In order to complete the new scheme of the derivation of
$\langle\sigma^{x}\rangle$, developed in this section, we provide the
algorithm of calculation of the $T_{\tau}$-exponent in Eqs. (D) and (D)
perturbatively in the small coupling parameter
$\kappa(=f_{2}\,\text{or}\,g_{2})$

$\displaystyle T_{\tau}\exp\Big{\\{}\kappa^{2}$
$\displaystyle\int_{0}^{\infty}d\xi\,\mathcal{I}(\xi)\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}G(\xi,\tau-\tau^{\prime})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})\Big{\\}}=$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{\kappa^{2n}}{n!}\left[\prod_{j=1}^{n}\int_{0}^{\infty}d\xi_{j}\,\mathcal{I}(\xi_{j})\int_{0}^{\beta}d\tau_{2j-1}\int_{0}^{\beta}d\tau_{2j}\,G(\xi_{j},\tau_{2j-1}-\tau_{2j})\right]\times$
$\displaystyle\times$ $\displaystyle
T_{\tau}\Big{[}\Sigma^{x}(\tau_{1})\Sigma^{x}(\tau_{2})\dots\Sigma^{x}(\tau_{2n})\Big{]}=$
$\displaystyle=$ $\displaystyle
1+\kappa^{2}\int_{0}^{\infty}d\xi_{1}\,\mathcal{I}(\xi_{1})\int_{0}^{\beta}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}\,\phi(\xi_{1},\tau_{1}-\tau_{2})\Sigma^{x}(\tau_{1})\Sigma^{x}(\tau_{2})+$
$\displaystyle+$
$\displaystyle\kappa^{4}\int_{0}^{\infty}d\xi_{1}\,\mathcal{I}(\xi_{1})\int_{0}^{\infty}d\xi_{2}\,\mathcal{I}(\xi_{2})\int_{0}^{\beta}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}\int_{0}^{\tau_{2}}d\tau_{3}\int_{0}^{\tau_{3}}d\tau_{4}\times$
$\displaystyle\times$
$\displaystyle\Big{[}\phi(\xi_{1},\tau_{1}-\tau_{2})\phi(\xi_{2},\tau_{3}-\tau_{4})+\phi(\xi_{1},\tau_{1}-\tau_{3})\phi(\xi_{2},\tau_{2}-\tau_{4})+$
$\displaystyle+\phi(\xi_{1},\tau_{1}-\tau_{4})\phi(\xi_{2},\tau_{2}-\tau_{3})\Big{]}\,\,\Sigma^{x}(\tau_{1})\Sigma^{x}(\tau_{2})\Sigma^{x}(\tau_{3})\Sigma^{x}(\tau_{4})+O(\kappa^{6}).$
(122)

Here $\mathcal{I}(\xi)$ represents the bosonic spectral density function,
which in our case is $\mathcal{I}_{1}(\xi)$/$\mathcal{I}_{2}(\xi)$ for the
left/right bath. Using the previously obtained result for the product of
operators $\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})$ one can derive the
leading correction term of the $T_{\tau}$-exponent, which is needed for the
calculation of the coherence $\langle\sigma^{x}\rangle$ in higher orders in
the coupling parameters $f_{2},g_{2}$

$\displaystyle\int_{0}^{\beta}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}\,\phi(\xi,\tau_{1}-\tau_{2})\Sigma^{x}(\tau)\Sigma^{x}(\tau^{\prime})=A\mathcal{I}+B\sigma^{z}+C\sigma^{x}+D\sigma^{y},$
(123)

where

$\displaystyle A=$
$\displaystyle\cos^{2}\theta_{mn}\frac{\beta\xi\left(\xi^{2}-R_{mn}^{2}\right)+2\left(\xi^{2}+R_{mn}^{2}\right)\coth\left(\frac{\beta\xi}{2}\right)\sinh^{2}\left(\frac{\beta
R_{mn}}{2}\right)}{(\xi^{2}-R_{mn}^{2})^{2}}-$ $\displaystyle-$
$\displaystyle\cos^{2}\theta_{mn}\frac{2\xi R_{mn}\sinh(\beta
R_{mn})}{(\xi^{2}-R_{mn}^{2})^{2}}+\sin^{2}\theta_{mn}\frac{\beta}{\xi},$
(124) $\displaystyle B=$
$\displaystyle\cos^{2}\theta_{mn}\frac{\left(\xi^{2}+R_{mn}^{2}\right)\coth\left(\frac{\beta\xi}{2}\right)\sinh(\beta
R_{mn})-\beta
R_{mn}^{3}\coth\left(\frac{\beta\xi}{2}\right)}{\left(\xi^{2}-R_{mn}^{2}\right)^{2}}+$
$\displaystyle+$ $\displaystyle\cos^{2}\theta_{mn}\frac{\xi
R_{mn}\left[\beta\xi\coth\left(\frac{\beta\xi}{2}\right)-4\cosh^{2}\left(\frac{\beta
R_{mn}}{2}\right)\right]}{\left(\xi^{2}-R_{mn}^{2}\right)^{2}},$ (125)
$\displaystyle C=$ $\displaystyle
2\sin\theta_{mn}\cos\theta_{mn}\frac{\xi\coth\left(\frac{\beta\xi}{2}\right)\sinh(\beta
R_{mn})-2R_{mn}\cosh^{2}\left(\frac{\beta
R_{mn}}{2}\right)}{\xi(\xi^{2}-R^{2}_{mn})},$ (126) $\displaystyle D=$
$\displaystyle
2i\sin\theta_{mn}\cos\theta_{mn}\frac{2\xi\coth\left(\frac{\beta\xi}{2}\right)\sinh^{2}\left(\frac{\beta
R_{mn}}{2}\right)-R_{mn}\sinh(\beta
R_{mn})}{\xi\left(\xi^{2}-R_{mn}^{2}\right)}.$ (127)

Integrating the result (123) with the spectral density function
$\mathcal{I}(\xi)$ over the parameter $\xi$ one gets the expression for
$T_{\tau}$-exponent up to quadratic terms of the parameter
$\kappa(=f_{2}\,\text{and}\,g_{2})$ including. Substituting this result into
Eqs. (D) and (D), and then taking the trace $\text{Tr}_{S}[\cdots]$ over
$\sigma$’s spin degrees of freedom one obtains next-to-leading-order
expression in coupling parameters $f_{2},g_{2}$ for the coherence.