# Dynamics of two central spins immersed in spin baths

Devvrat Tiwari<EMAIL_ADDRESS>Indian Institute of Technology Jodhpur,
Rajasthan, India  Shounak Datta  S. N. Bose National Centre for Basic
Sciences, Block JD, Sector III, Salt Lake, Kolkata 700 098, India  Samyadeb
Bhattacharya<EMAIL_ADDRESS>Center for Security Theory and Algorithmic
Research, International Institute of Information Technology, Gachibowli,
Hyderabad, India  Subhashish Banerjee<EMAIL_ADDRESS>Indian Institute
of Technology Jodhpur, Rajasthan, India

###### Abstract

In this article we derive the exact dynamics of a two qubit (spin $1/2$)
system interacting centrally with separate spin baths composed of qubits in
thermal state. Further, each spin of a bath is coupled to every other spin of
the same bath. The corresponding dynamical map is constructed. It is used to
analyse the non-Markovian nature of the two qubit central spin dynamics. We
further observe the evolution of quantum correlations like entanglement and
discord under the influence of the environmental interaction. Moreover, we
demonstrate the comparison between this exact two qubit dynamics and the
locally acting central spin model in a spin bath. This work is a stepping
stone towards the realization of non-Markovian heat engines and other quantum
thermal devices.

###### pacs:

03.65.Ud, 03.67.Mn, 03.65.Ta

##  I Introduction

In the interactive world of quantum particles, it is almost impossible to
create a physical system devoid of any external noisy influence. Quantum
systems conducive for information theoretic tasks, i.e. trapped ions (Cirac
and Zoller, 1995), quantum dots (Loss and DiVincenzo, 1998), NMR qubits (Cory
_et al._ , 1997), Josephson junctions (Rauch, 1998) and many more, are
subjected to environmental interactions. It is therefore a necessary task to
investigate the dynamics of such quantum systems, which are under the
influence of environmental interactions. The theory of open quantum systems
(Breuer and Petruccione, 2002; Banerjee, 2018), has found numerous
applications, in recent times, in quantum information and its interface with
various facets of quantum physics (Leggett  _et al._ , 1987a; Gorini  _et al._
, 1976; Lindblad, 1976; Caldeira and Leggett, 1983; Feynman and Vernon, 1963;
Weiss, 1999; Leggett  _et al._ , 1987b; Nakajima, 1958; Zwanzig, 1960;
Chaturvedi and Shibata, 1979; Banerjee and Ghosh, 2000, 2003a; Chandrashekar
_et al._ , 2007; Banerjee  _et al._ , 2008; Thomas  _et al._ , 2018; Naikoo
_et al._ , 2019; Alok  _et al._ , 2016; Naikoo  _et al._ , 2018; Shrikant  _et
al._ , 2018; Dixit  _et al._ , 2019). Over the past few decades, our
understanding of such systems has stretched from the limitations of Markovian
dynamics to the more intriguing and challenging domain of non-Markovian
quantum systems (Breuer  _et al._ , 2004; Laine  _et al._ , 2010; Rivas  _et
al._ , 2010; Chruściński  _et al._ , 2011; Vasile  _et al._ , 2011; Lu  _et
al._ , 2010; Luo  _et al._ , 2012; Fanchini, 2014; Chanda and Bhattacharya,
2016; Haseli, 2014; Hall  _et al._ , 2014; Andersson  _et al._ , 2007;
Bhattacharya  _et al._ , 2017; Mukhopadhyay  _et al._ , 2017; Awasthi  _et
al._ , 2018; Bhattacharya  _et al._ , 2020a, b; Maity  _et al._ , 2020;
Bhattacharya and Bhattacharya, 2021; Bhattacharya  _et al._ , 2018; Das  _et
al._ , 2021; Banerjee and Ghosh, 2003b; Kumar  _et al._ , 2018a, b; Utagi  _et
al._ , 2020; Krovi  _et al._ , 2007). Even now, it is a challenging task to
construct the reduced dynamics of such a system without the Born-Markov and
stationary bath approximations (Breuer and Petruccione, 2002). One often
associates a deviation from quantum dynamical semigroup evolution of a system
to a non-Markovian process (Utagi  _et al._ , 2020). In a non-Markovian
evolution, the time scales of the system and the environment are often not
well separated, which can result in information backflow from the environment
to the system Breuer  _et al._ (2009); Rivas  _et al._ (2010). This generally
leads to recurrences of quantum properties which is important for a
fundamental understanding of system dynamics.

It would be pertinent to add here that a number of techniques have been
developed in recent times to tackle this problem (de Vega and Alonso, 2017).
Thus, for example, there are the embedding methods, such as the pseudo-modes
technique wherein the decay of an atom strongly coupled to a reservoir can be
studied by considering an enlarged system that includes a set of pseudo-modes
(Garraway and Knight, 1996; Dalton and Garraway, 2003). Another relevant
technique is the reaction-coordinate mapping (Iles-Smith  _et al._ , 2014;
Garg  _et al._ , 1985; Hughes  _et al._ , 2009). Further, a numerically exact
hierarchical equations of motion (HEOM) method has been developed (Tanimura
and Kubo, 1989).

Quantum baths are generally categorized in two broad classes; (a) Bosonic and
(b) spin bath. Archetypal examples of Bosonic baths include Caldeira-Leggett
model (Caldeira and Leggett, 1983) or spin boson model (Leggett  _et al._ ,
1987b). Exact quantum master equations for these type of models are common in
the literature (Breuer and Petruccione, 2002). On the contrary, in the case of
spin baths, we often have to rely on perturbative techniques or time non-local
master equations (Breuer  _et al._ , 2004; Fischer and Breuer, 2007).
Although, the study of these systems are extremely important in physical
systems of paramount importance such as magnetic systems, quantum spin glasses
and superconducting systems (Breuer  _et al._ , 2004), the theoretical
modelling of such systems is still lacking in many different perspectives,
especially in bi-partite or multi-partite systems relevant in various quantum
device modelling. In this work, we attempt to lay the bedrock of such a
construction from the perspective of modelling various quantum devices. Here
we develop an exact reduced dynamics of a two qubit system immersed in spin
baths, each of which is interacting centrally with the system of interest. Our
model Hamiltonian of the two qubit system is inspired from a model of quantum
thermal diode (Ordonez-Miranda  _et al._ , 2017). In recent times, motivated
by the goal of building quantum computers, a lot of effort has been given to
develop quantum versions of various thermodynamic and electronic devices like
refrigerators and heat engines (Linden  _et al._ , 2010; Wang  _et al._ ,
2022a; Ar ısoy  _et al._ , 2022; Mohanta  _et al._ , 2022), thermal diodes
(Ordonez-Miranda  _et al._ , 2017; Karg ı  _et al._ , 2019), transistors
(Joulain  _et al._ , 2016; Wijesekara  _et al._ , 2021; Majland  _et al._ ,
2020; Ghosh  _et al._ , 2021), quantum batteries and so on (Wang  _et al._ ,
2022b; Bhandari  _et al._ , 2021). Drawing the motivation from the long term
goal of realizing experimentally feasible models of quantum thermal devices,
we develop the aforementioned exact reduced dynamics of a two qubit spin
system. We further analyse various thermodynamic and information theoretic
properties of the system undergoing the specific open quantum evolution. We
also observe the fundamentally non-Markovian behaviour of the dynamical map.

The flow of the paper is as follows: In Sec. II we introduce the model. Its
reduced dynamics is developed in Sec. III. The corresponding dynamical map is
discussed next, along with its operator sum representation. The dynamical map
is then put to use in Sec. V. This includes analyzing a witness to identify
inherent non-Markovianity in the dynamics, the corresponding local and global
dynamical maps and finally the quantum correlations, including entanglement
and discord, generated between the two spins by their open system dynamics.
This is followed by the Conclusion.

##  II The Model

We present a model of two coupled qubits where each qubit is centrally coupled
to different thermal spin baths (see Fig. 1). Further, each spin of a bath is
coupled to every other spin of the same bath.

Figure 1: (Colour online) Schematic diagram of coupled central spin model
where each central spin interacts with individual spin baths.

Initially, $\rho_{SB}(0)=\rho_{S}(0)\otimes\rho_{B}(0)$, where,
$\rho_{B}(0)=\frac{e^{-H_{B}/K_{B}T}}{Z}$ for each central spin immersed in a
spin bath. The evolution of the whole system is governed by the following
Hamiltonian in our spin-bath model:

$\displaystyle H$
$\displaystyle=H_{S_{1}}+H_{S_{2}}+H_{S_{1}S_{2}}+H_{B_{1}}+H_{B_{2}}+H_{S_{1}B_{1}}+H_{S_{2}B_{2}},$
$\displaystyle=\frac{\hslash\omega_{1}}{2}\sigma_{1z}^{0}+\frac{\hslash\omega_{2}}{2}\sigma_{2z}^{0}+\frac{\hslash\delta}{2}(\sigma_{1z}^{0}\otimes\sigma_{2z}^{0})$
$\displaystyle+\frac{\hslash\omega_{a}}{2M}\sum_{i=1}^{M}\Big{\\{}\frac{1}{2}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{M}(\sigma_{1x}^{i}\sigma_{1x}^{j}+\sigma_{1y}^{i}\sigma_{1y}^{j})+\sigma_{1z}^{i}\Big{\\}}$
$\displaystyle+\frac{\hslash\omega_{b}}{2N}\sum_{i=1}^{N}\Big{\\{}\frac{1}{2}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{N}(\sigma_{2x}^{i}\sigma_{2x}^{j}+\sigma_{2y}^{i}\sigma_{2y}^{j})+\sigma_{2z}^{i}\Big{\\}}$
$\displaystyle+\frac{\hslash\epsilon_{1}}{2\sqrt{M}}\sum_{i=1}^{M}(\sigma_{1x}^{0}\sigma_{1x}^{i}+\sigma_{1y}^{0}\sigma_{1y}^{i})+\frac{\hslash\epsilon_{2}}{2\sqrt{N}}\sum_{i=1}^{N}(\sigma_{2x}^{0}\sigma_{2x}^{i}+\sigma_{2y}^{0}\sigma_{2y}^{i}),$
(1)

where $\sigma_{lk}^{i}$ or $\sigma_{lk}^{j}$ (k=x,y,z; $l$=1,2) are Pauli
matrices corresponding to $i$-th or $j$-th spin of the $l$-th bath and
$\sigma_{lk}^{0}$ (k=x,y,z; $l$=1,2) are same for $l$-th central spin.
$\omega_{1}$ and $\omega_{2}$ are the two central spin frequencies and
$\delta$ corresponds to their coupling strength. Also, $\omega_{a}$ and
$\omega_{b}$ are bath frequencies of two spin baths and $\epsilon_{l}$
($l$=1,2) are interaction parameters of $l$-th system-bath. $M$ and $N$ are
the respective number of atoms in two baths.

Using total spin angular momentum operator,
$J_{lk^{\prime}}=\frac{1}{2}\sum_{i=1}^{MorN}\sigma_{lk^{\prime}}^{i}$
($k^{\prime}$=x,y,z,+,-; $l$=1 or 2 corresponding to summation upper limit M
or N), we may rewrite bath Hamiltonians as,

$\displaystyle
H_{B_{1}}=\hslash\omega_{a}(\frac{J_{1+}J_{1-}}{M}-\frac{\openone}{2}),$
$\displaystyle
H_{B_{2}}=\hslash\omega_{b}(\frac{J_{2+}J_{2-}}{N}-\frac{\openone}{2}),$ (2)

and system-bath interaction Hamiltonians as,

$\displaystyle
H_{S_{1}B_{1}}=\frac{\hslash\epsilon_{1}}{\sqrt{M}}(\sigma_{1x}^{0}J_{1x}+\sigma_{1y}^{0}J_{1y}),$
$\displaystyle
H_{S_{2}B_{2}}=\frac{\hslash\epsilon_{2}}{\sqrt{N}}(\sigma_{2x}^{0}J_{2x}+\sigma_{2y}^{0}J_{2y}).$
(3)

Following Bhattacharya and Banerjee (2021), we then use the Holstein-Primakoff
transformation (Holstein and Primakoff, 1940) to redefine collective angular
momentum operators as,

$\displaystyle
J_{1+}=\sqrt{M}a^{\dagger}\Big{(}1-\frac{a^{\dagger}a}{M}\Big{)}^{1/2},$
$\displaystyle J_{1-}=\sqrt{M}\Big{(}1-\frac{a^{\dagger}a}{M}\Big{)}^{1/2}a,$
(4)

for first bath and the following for the second bath

$\displaystyle
J_{2+}=\sqrt{N}b^{\dagger}\Big{(}1-\frac{b^{\dagger}b}{N}\Big{)}^{1/2},$
$\displaystyle J_{2-}=\sqrt{N}\Big{(}1-\frac{b^{\dagger}b}{N}\Big{)}^{1/2}b.$
(5)

Here $a$ and $a^{\dagger}$ are the bosonic annihilation and creation operators
for the first spin bath having the property, $[a,a^{\dagger}]=1$ and $b$ and
$b^{\dagger}$ represents the same for the second spin bath. After this
transformation, the bath Hamiltonians appear as,

$\displaystyle
H_{B_{1}}=\hslash\omega_{a}\Big{\\{}a^{\dagger}a\Big{(}1-\frac{a^{\dagger}a-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}},$
$\displaystyle
H_{B_{2}}=\hslash\omega_{b}\Big{\\{}b^{\dagger}b\Big{(}1-\frac{b^{\dagger}b-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}},$
(6)

and the interaction Hamiltonians of respective spin-baths as,

$\displaystyle
H_{S_{1}B_{1}}=\hslash\epsilon_{1}\Big{\\{}\sigma_{1+}^{0}\Big{(}1-\frac{a^{\dagger}a}{M}\Big{)}^{1/2}a+\sigma_{1-}^{0}a^{\dagger}\Big{(}1-\frac{a^{\dagger}a}{M}\Big{)}^{1/2}\Big{\\}},$
$\displaystyle
H_{S_{2}B_{2}}=\hslash\epsilon_{2}\Big{\\{}\sigma_{2+}^{0}\Big{(}1-\frac{b^{\dagger}b}{N}\Big{)}^{1/2}b+\sigma_{2-}^{0}b^{\dagger}\Big{(}1-\frac{b^{\dagger}b}{N}\Big{)}^{1/2}\Big{\\}}.$
(7)

It is also important to mention the limitations and the essential
approximations. Here we are taking a central spin half system interacting
homogeneously with each of the bath spins with all the characteristic
frequencies to be of constant value. We have taken homogeneous interactions
for the sake of analytical clarity. This is not an oversimplified assumption,
as these kinds of interactions exist in physical situations, of which some
examples are quantum spin glasses, superconducting systems and NMR. The more
general system with in-homogeneous interaction parameters can only be handled
numerically. A use is then made of the Holstein-Primakoff transformation. A
Hamiltonian describing the collective behaviour of $N$ interacting spins can
be mapped to a bosonic one employing this transformation, at the expense of
having an infinite series in powers of the bosonic creation and annihilation
operators. Truncating this series up to quadratic terms allows for obtaining
analytic solutions, which become exact in the limit $N\rightarrow\infty$. In
the literature, works on similar spin environments exists which make use of
different methods, cf. (Breuer  _et al._ , 2004). The reason we choose
Holstein-Primakoff transformation here is because of its technical advantage,
both from analytical and numerical perspectives. We can work with more number
of bath spins with lesser numerical limitations. Homogeneous interactions
further enable us to modify the total Hamiltonian into a form of non-linear
Jaynes-Cummings model. Other than considering homogeneous interactions, no
further assumptions are needed for the type of system-bath model studied here.
This opens up an opportunity to study more complex quantum devices and
networks.

##  III Reduced dynamics of the two qubit central spin model

We derive the reduced dynamics of the system of two coupled central spins
evolved under the Hamiltonian $H$ along with the two baths by tracing over the
bath degrees of freedom.

Consider the evolution of the state $|\psi(0)\rangle=|11\rangle|xy\rangle$,
where two central spins are in the excited state $|1\rangle$ and $|x\rangle$
is an arbitrary state for the 1st bath while $|y\rangle$ is a state belonging
to the 2nd bath. Global unitary operator corresponding to the evolution under
Hamiltonian $H$ can be written as, $U(t)=\exp(-\frac{iHt}{\hslash})$. After
the evolution let the state at time $t$ be
$|\psi(t)\rangle=\varsigma_{1}(t)|11\rangle|x^{\prime}y^{\prime}\rangle+\varsigma_{2}(t)|10\rangle|x^{\prime\prime}y^{\prime\prime}\rangle+\varsigma_{3}(t)|01\rangle|x^{\prime\prime\prime}y^{\prime\prime\prime}\rangle$.
We exclude the case when both spins of the system are flipped simultaneously;
i.e., $|11\rangle\rightarrow|00\rangle$ transition. The transitions we have
considered convey the message sufficiently and hence we choose to exclude the
aforementioned transition to avoid unnecessary complications in the
calculations. Now we introduce three time-dependent operators
$\hat{A}(t),\hat{B}(t)$ and $\hat{C}(t)$ corresponding to the joint Hilbert
space of the two baths such that,
$\hat{A}(t)|xy\rangle=\varsigma_{1}(t)|x^{\prime}y^{\prime}\rangle$,
$\hat{B}(t)|xy\rangle=\varsigma_{2}(t)|x^{\prime\prime}y^{\prime\prime}\rangle$
and
$\hat{C}(t)|xy\rangle=\varsigma_{3}(t)|x^{\prime\prime\prime}y^{\prime\prime\prime}\rangle$.
Then we have,
$|\psi(t)\rangle=\hat{A}(t)|11\rangle|xy\rangle+\hat{B}(t)|10\rangle|xy\rangle+\hat{C}(t)|01\rangle|xy\rangle$.

From the time-dependent Schrödinger equation
$\frac{d}{dt}|\psi(t)\rangle=-\frac{iH}{\hslash}|\psi(t)\rangle$, we can have,

$\displaystyle\frac{d\hat{A}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}+\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}a^{\dagger}a\Big{(}1-\frac{a^{\dagger}a-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}b^{\dagger}b\Big{(}1-\frac{b^{\dagger}b-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{A}(t)$
$\displaystyle-i\epsilon_{2}\Big{(}1-\frac{b^{\dagger}b}{N}\Big{)}^{1/2}b\hat{B}(t)-i\epsilon_{1}\Big{(}1-\frac{a^{\dagger}a}{M}\Big{)}^{1/2}a\hat{C}(t),$
$\displaystyle\frac{d\hat{B}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}-\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}a^{\dagger}a\Big{(}1-\frac{a^{\dagger}a-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}b^{\dagger}b\Big{(}1-\frac{b^{\dagger}b-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{B}(t)$
$\displaystyle-i\epsilon_{2}b^{\dagger}\Big{(}1-\frac{b^{\dagger}b}{N}\Big{)}^{1/2}\hat{A}(t),$
$\displaystyle\frac{d\hat{C}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{-\omega_{1}+\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}a^{\dagger}a\Big{(}1-\frac{a^{\dagger}a-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}b^{\dagger}b\Big{(}1-\frac{b^{\dagger}b-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{C}(t)$
$\displaystyle-i\epsilon_{1}a^{\dagger}\Big{(}1-\frac{a^{\dagger}a}{M}\Big{)}^{1/2}\hat{A}(t).$
(8)

Now we substitute $\hat{A}(t)=\hat{A_{1}}(t)$,
$\hat{B}(t)=b^{\dagger}\hat{B_{1}}(t)$ and
$\hat{C}(t)=a^{\dagger}\hat{C_{1}}(t)$ and have,

$\displaystyle\frac{d\hat{A_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}+\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{A_{1}}(t)$
$\displaystyle-i\epsilon_{2}\Big{(}1-\frac{\hat{n}}{N}\Big{)}^{1/2}(\hat{n}+1)\hat{B_{1}}(t)$
$\displaystyle-i\epsilon_{1}\Big{(}1-\frac{\hat{m}}{M}\Big{)}^{1/2}(\hat{m}+1)\hat{C_{1}}(t),$
$\displaystyle\frac{d\hat{B_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}-\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(\hat{n}+1)\Big{(}1-\frac{\hat{n}}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{B_{1}}(t)$
$\displaystyle-i\epsilon_{2}\Big{(}1-\frac{\hat{n}}{N}\Big{)}^{1/2}\hat{A_{1}}(t),$
$\displaystyle\frac{d\hat{C_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{-\omega_{1}+\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}(\hat{m}+1)\Big{(}1-\frac{\hat{m}}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{C_{1}}(t)$
$\displaystyle-i\epsilon_{1}\Big{(}1-\frac{\hat{m}}{M}\Big{)}^{1/2}\hat{A_{1}}(t).$
(9)

Here $\hat{m}=a^{\dagger}a$ and $\hat{n}=b^{\dagger}b$ are number operators
corresponding to bosonic operators of 1st and 2nd bath, respectively, which
have the properties $\hat{m}|m\rangle=m|m\rangle$ and
$\hat{n}|n\rangle=n|n\rangle$. Therefore, we define,
$\hat{A_{1}}(t)|mn\rangle=A_{1}(m,n,t)|mn\rangle$,
$\hat{B_{1}}(t)|mn\rangle=B_{1}(m,n,t)|mn\rangle$ and
$\hat{C_{1}}(t)|mn\rangle=C_{1}(m,n,t)|mn\rangle$.

By tracing out bath modes, reduced state of the system ($|11\rangle\langle
11|$) becomes

$\displaystyle\phi(|11\rangle\langle 11|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\psi(t)\rangle\langle\psi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{\\{}|A_{1}(m,n,t)|^{2}|11\rangle\langle
11|$ $\displaystyle+(n+1)|B_{1}(m,n,t)|^{2}|10\rangle\langle 10|$
$\displaystyle+(m+1)|C_{1}(m,n,t)|^{2}|01\rangle\langle 01|\Big{\\}}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}]$
(10)

From the solution of Eq. (9), one can derive the values of
$|A_{1}(m,n,t)|^{2}$, $|B_{1}(m,n,t)|^{2}$ and $|C_{1}(m,n,t)|^{2}$ following
the steps as mentioned in the Appendix B. The partition function in the above
equation is,
$Z=\sum_{m}^{M}\sum_{n}^{N}\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]~{}\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}]$.

In a similar way, we can derive the evolution of the other elements of the
system’s reduced state. The details are given in Appendix A. The reduced state
of the system of two central spins after the global unitary evolution of the
pair of joint system-bath state can be denoted as

$\displaystyle\rho_{S_{1}S_{2}}$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}\Big{[}e^{-iHt/\hslash}\rho_{SB}(0)e^{iHt/\hslash}\Big{]}$
$\displaystyle=\begin{pmatrix}\rho_{11}(t)&\rho_{12}(t)&\rho_{13}(t)&\rho_{14}(t)\\\
\rho_{21}(t)&\rho_{22}(t)&\rho_{23}(t)&\rho_{24}(t)\\\
\rho_{31}(t)&\rho_{32}(t)&\rho_{33}(t)&\rho_{34}(t)\\\
\rho_{41}(t)&\rho_{42}(t)&\rho_{43}(t)&\rho_{44}(t)\\\ \end{pmatrix},$ (11)

where $\rho_{SB}(0)$ is the joint system-bath initial state. The elements of
the density matrix are given in Appendix C.

##  IV Construction of the dynamical map

Having constructed the density matrix of the reduced two-qubit central spin
system, we now find the dynamical map of the system. To this effect, we derive
the Kraus operators for the evolution of the reduced system.

###  IV.1 Operator sum representation

An important facet of general quantum evolution represented by a CPTP
(completely positive trace preserving) operation is the Kraus operator sum
representation, given as $\rho(t)=\sum_{i}K_{i}(t)\rho(0)K_{i}^{\dagger}(t)$.
The Kraus operators, $K_{i}^{\prime}s$ can be constructed from the eigenvalues
and eigenvectors of the Choi-Jamiołkowski (CJ) state of the corresponding
dynamical map. The CJ state of the dynamic map $\phi(\cdot)$ acting on a
$d-$dimensional system is given by
$({\mathbb{I}}_{d}\otimes\phi)|\psi\rangle\langle\psi|$, with $\ket{\psi}$
being the maximally entangled state in $d^{2}$ dimension. For our particular
case, the CJ matrix is given by,

$\displaystyle\mathcal{C}(t)=\setcounter{MaxMatrixCols}{16}\begin{pmatrix}|A_{1}|^{2}&0&0&0&0&A_{1}J_{1}^{*}&0&0&0&0&A_{1}G_{1}^{*}&0&0&0&0&A_{1}D_{1}^{*}\\\
0&(n+1)|B_{1}|^{2}&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&(m+1)|C_{1}|^{2}&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&n|K_{1}|^{2}&0&0&0&0&0&0&0&0&0&0&0\\\
A_{1}^{*}J_{1}&0&0&0&0&|J_{1}|^{2}&0&0&0&0&J_{1}G_{1}^{*}&0&0&0&0&J_{1}D_{1}^{*}\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&(m+1)|L_{1}|^{2}&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&m|I_{1}|^{2}&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
A_{1}^{*}G_{1}&0&0&0&0&J_{1}^{*}G_{1}&0&0&0&0&|G_{1}|^{2}&0&0&0&0&G_{1}D_{1}^{*}\\\
0&0&0&0&0&0&0&0&0&0&0&(n+1)|H_{1}|^{2}&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&m|F_{1}|^{2}&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0&0&0&n|E_{1}|^{2}&0\\\
A_{1}^{*}D_{1}&0&0&0&0&J_{1}^{*}D_{1}&0&0&0&0&G_{1}^{*}D_{1}&0&0&0&0&|D_{1}|^{2}\end{pmatrix}.$
(12)

The symbols and abbreviations used in Eq. (12) are explained in appendices A
and D. One can find the eigenspectrum of the matrix in Eq. (12) numerically,
and from there the Kraus operators for the evolution of the state can be
derived. These Kraus operators should satisfy the relation
$\sum_{i}K_{i}^{\dagger}K_{i}=\mathbb{I}$, which is indeed the case here.

##  V Analysis of the dynamical map

In this section, we use the dynamical map of the reduced state of the system
to shed light on the dynamics of the two-qubit central spin model.
Particularly we are interested in the study of its non-Markovian behaviour
from the backdrop of information backflow (Breuer  _et al._ , 2009). The study
of information backflow is done by observing the time evolution of certain
distance functions like trace distance or fidelity between two quantum states,
which is a monotonically decreasing function under Markovian evolution. Any
deviation from their monotonic behavior is thus considered as a signature of
non-Markovianity. Here we study the evolution of trace distance to investigate
the non-Markovian nature of the dynamics. In this context, it is interesting
to observe the difference between the action of the dynamical map in their
local and global avatars. Along with this, we also investigate the evolution
of non-classical correlations such as quantum entanglement and quantum
discord, to determine the sustenance of such resourceful quantum properties
under the dynamics in question.

###  V.1 Trace distance as a witness of non-Markovianity

The trace distance, which gives a measure of distinguishability between two
quantum states, is given as
${\mathcal{D}}(\rho_{1},\rho_{2})=\frac{1}{2}||\rho_{1}-\rho_{2}||_{1}$, where
$||(\cdot)||_{1}={\rm Tr}\sqrt{(\cdot)^{\dagger}(\cdot)}$. To this end, we
calculate the trace distance between the state evolved with time through the
dynamical map given above and the initial state, and is given by,

$\displaystyle{\mathcal{D}}(\rho_{S_{1}S_{2}}(t),\rho_{S_{1}S_{2}}(0))=\frac{1}{2}||\rho_{S_{1}S_{2}}(t)-\rho_{S_{1}S_{2}}(0)||_{1}.$
(13)

Figure 2: (Colour online) Variation of trace distance, between time evolved
state and initial state, with time for the reduced state of the system. In
(a), the initial state is taken to be $\ket{11}$, and in (b), the initial
state is taken to be $\ket{10}$. The parameters have following values:
$\omega_{1}=2.0$, $\omega_{2}=1.9$, $\delta=2.5$, $\omega_{a}=1.1$,
$\omega_{b}=1.2$, $M=N=100$, $T=1$, $\epsilon_{1}=2.6$, $\epsilon_{2}=2.5$.

In Fig. 2, we can see the variation of
${\mathcal{D}}(\rho_{S_{1}S_{2}}(t),\rho_{S_{1}S_{2}}(0))$ with time for
different initial states. The non monotonic behaviour of the trace distance
between time evolved state and the initial state is a witness of non-
Markovianity in the system.

###  V.2 Difference between local and global dynamical maps

In this work, we have derived the global dynamical map of two central spins
$(\Lambda_{12})$. Let $\Lambda_{1}$, $\Lambda_{2}$ be the local dynamical maps
derived by solving the local Lindblad equations for each central spin
(Mukhopadhyay  _et al._ , 2017), with the added proviso that the bath spins
are interacting with each other. $\Lambda_{12}$ is the global map constructed
here. We take the same parameter values for both the global and local maps, so
that we can observe the difference between them from a common footing. For a
particular initial state $\rho_{s}(0)$ we can calculate the trace distance as,

${\mathcal{D}}(\rho_{\text{global}}(t),\rho_{\text{local}}(t))=\frac{1}{2}||\Lambda_{12}(\rho_{s}(0))-\Lambda_{1}\otimes\Lambda_{2}(\rho_{s}(0))||_{1}.$
(14)

Figure 3: (Colour online) Variation of trace distance for difference between
local and global maps with time for the reduced state of the system. In (a),
the initial state is taken to be $\ket{11}$, and in (b), the initial state is
taken to be $\ket{10}$. The parameters have following values:
$\omega_{1}=2.0$, $\omega_{2}=1.9$, $\delta=5$, $\omega_{a}=1.1$,
$\omega_{b}=1.2$, $M=N=100$, $T=1$, $\epsilon_{1}=2.6$, $\epsilon_{2}=2.5$.

In Fig. 3, we can see the difference between local and global dynamical maps
by setting the same values for all the parameters. It is evident from the plot
that though there is no interaction between the two baths in the dynamics we
are evaluating, and the two baths are acting separately with individual
qubits, there is a distinct difference in the dynamical behaviour of this map
with the local maps acting separately with each of the qubits. Through the
interaction between the qubits, bath information is passed from one
environment to another and hence, despite being mutually non-interacting, the
action of the bath exhibits a global trait (Banerjee  _et al._ , 2010).
Therefore, it is our understanding that, in those situations where baths are
acting locally on a bipartite system, as presented in Eq. (II), applying local
Kraus operators for each separate system-bath interaction, does not give us
the complete picture.

###  V.3 Quantum Correlations

Quantum correlations are a very useful resource in quantum information
processing. To this end, we investigate the quantum correlations in the
reduced state of the system given in Eq. (11). Concurrence (Wootters, 1998) is
a measure of quantum entanglement in the system given by,

$\mathcal{Q_{C}}=\max\\{0,\lambda_{1}-\lambda_{2}-\lambda_{3}-\lambda_{4}\\},$
(15)

where $\lambda_{i}^{\prime}s$ are the eigenvalue of the matrix
$\sqrt{\sqrt{\rho_{S_{1}S_{2}}}\tilde{\rho}_{S_{1}S_{2}}\sqrt{\rho_{S_{1}S_{2}}}}$
in decreasing order, and
$\tilde{\rho}_{S_{1}S_{2}}=(\sigma_{y}\otimes\sigma_{y})\rho^{*}(\sigma_{y}\otimes\sigma_{y})$.

Figure 4: (Colour online) Variation of concurrence ($\mathcal{Q_{C}}$) and
quantum discord ($\mathcal{Q_{D}}$) with time for the reduced state of the
system. The parameters have following values: $\omega_{1}=2.0$,
$\omega_{2}=1.9$, $\delta=3$, $\omega_{a}=1.1$, $\omega_{b}=1.2$, $M=N=100$,
$T=1$, $\epsilon_{1}=1.3$, $\epsilon_{2}=1.25$.

Another popular candidate for the measurement of quantum correlations is
quantum discord (Ollivier and Zurek, 2001). It includes the quantum
correlations due to quantum effects in the system which may not necessarily be
due to quantum entanglement. For the reduced state given in Eq. (11), quantum
discord can be given as,

$\mathcal{Q_{D}}=S(\rho_{S_{2}})-S(\rho_{S_{1}S_{2}})+S(\rho_{S_{1}|S_{2}}),$
(16)

where $S(\rho_{S_{2}})$ and $S(\rho_{S_{1}S_{2}})$ are the von Neumann entropy
of the reduced subsystem $Tr_{S_{1}}[\rho_{S_{1}S_{2}}]$ and joint von Neumann
entropy of the reduced system $\rho_{S_{1}S_{2}}$, respectively.
$S(\rho_{S_{1}|S_{2}})$ is the quantum conditional entropy given by,

$S(\rho_{S_{1}|S_{2}})=\min_{\\{\Pi_{k}\\}}\sum_{k=1}^{2}p_{k}S(\rho_{S_{1}|\Pi_{k}}),$
(17)

where $\rho_{S_{1}|\Pi_{k}}$ is the state of the reduced system when
measurement operator $\Pi_{k}$ is operated on subsystem $S_{2}$, such that
$\rho_{S_{1}|\Pi_{k}}=\frac{1}{p_{k}}Tr_{S_{2}}(\Pi_{k}\rho_{S_{1}S_{2}}\Pi_{k}^{\dagger})$.
$p_{k}$ is the probability associated with the measurement operators $\Pi_{k}$
given by $p_{k}=Tr(\Pi_{k}\rho_{S_{1}S_{2}}\Pi_{k}^{\dagger})$. The
generalized measurement operators $\Pi_{k}$ in 2 qubits, are given by
$\Pi_{1}=\mathbb{I}_{S_{1}}\otimes|u\rangle\langle u|_{S_{2}}$ and
$\Pi_{2}=\mathbb{I}_{S_{1}}\otimes|v\rangle\langle v|_{S_{2}}$, where
$\ket{u}=\cos(\theta)\ket{1}+e^{i\phi}\sin(\theta)\ket{0}$ and
$\ket{v}=\sin(\theta)\ket{1}-e^{i\phi}\cos(\theta)\ket{0}$. The parameters
$\theta$ and $\phi$ vary in the range $0\leq\theta\leq\pi/2$ and
$0\leq\phi\leq 2\pi$. We consider the two-qubit maximally entangled state,
$\ket{\psi(0)}_{S}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$ as the initial state
of the central spin system and study the variation of quantum correlations
with time, depicted in Fig. 4. The profile of both the quantum discord and the
concurrence is similar. There exist time intervals where entanglement is zero,
but discord is non-zero. Revivals in the quantum correlations depicting the
information backflow nature of non-Markovianity of the system are observed. It
is also important to mention that, here we are considering finite bath spins
in our numerical analysis. If we extend this to the thermodynamic limit of
$M,~{}N\rightarrow\infty$, the dynamical behavior will fall into the Markovian
regime. In fact, in a previous work (Mukhopadhyay  _et al._ , 2017), one of
the present authors has shown that with increasing bath spins the dynamical
behavior asymptotically reaches the Markovian situation. Therefore, a natural
conclusion is that non-Markovianity is a result of the finite-ness of the spin
environment. On the other hand, to eliminate the possible conclusion that non-
Markovianity is an artifact of the truncated Holstein-Primakoff
transformation, we refer to a previous work on similar spin environment
treated by different methods, but also exhibiting typical non-Markovian
behaviour Breuer  _et al._ (2004). This shows that non-Markovianity is not an
artifact of the method.

##  VI Conclusion

In this article, we have derived the exact reduced dynamics of two coupled
central spins where each spin is centrally coupled to different thermal baths.
We develop the corresponding dynamical map by constructing the relevant Kraus
operators. This helps to shed light into the reduced dynamics of the coupled
central spin model. In particular, we show evidence of the non-Markovian
evolution of the system in the form of trace distance. We also calculate the
difference between the local and global dynamical maps of evolution and show
that there is a distinct difference in the exact dynamics of the two qubit
system interacting separately with two non-interacting environments and that
of the phenomenological application of local Kraus operations of the same
physical picture. It is our assertion that the later does not give us the
complete picture of the exact dynamical behaviour of such interactions, where
we overlook the information flow between the baths via the interaction between
the qubit systems. Moreover, we have also seen the effect on the quantum
correlations, in particular entanglement and quantum discord, between the
central spins as the dynamics of the system evolves in time. Revival of the
quantum correlations benchmarks the non-Markovian behaviour in the central
spin system.

## Acknowledgements

SB acknowledges support from Interdisciplinary Cyber Physical Systems (ICPS)
programme of the Department of Science and Technology (DST), India, Grant No.:
DST/ICPS/QuST/Theme-1/2019/13. SB also acknowledges support from the
Interdisciplinary Research Platform (IDRP) on Quantum Information and
Computation (QIC) at IIT Jodhpur.

##  Appendix A Derivation of the elements of the reduced state of the system.

Following the procedure in Sec. III, we now define,
$|\chi(0)\rangle=|00\rangle|xy\rangle$ and
$|\chi(t)\rangle=\hat{D}(t)|00\rangle|xy\rangle+\hat{E}(t)|01\rangle|xy\rangle+\hat{F}(t)|10\rangle|xy\rangle$.
Substituting $\hat{D}(t)=\hat{D_{1}}(t)$, $\hat{E}(t)=b\hat{E_{1}}(t)$ and
$\hat{F}(t)=a\hat{F_{1}}(t)$, we find,

$\displaystyle\frac{d\hat{D_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{-\omega_{1}-\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{D_{1}}(t)$
$\displaystyle-i\epsilon_{2}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}^{1/2}\hat{E_{1}}(t)$
$\displaystyle-i\epsilon_{1}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}^{1/2}\hat{F_{1}}(t),$
$\displaystyle\frac{d\hat{E_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{-\omega_{1}+\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(\hat{n}-1)\Big{(}1-\frac{\hat{n}-2}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{E_{1}}(t)$
$\displaystyle-i\epsilon_{2}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}^{1/2}\hat{D_{1}}(t),$
$\displaystyle\frac{d\hat{F_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}-\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}(\hat{m}-1)\Big{(}1-\frac{\hat{m}-2}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{F_{1}}(t)$
$\displaystyle-i\epsilon_{1}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}^{1/2}\hat{D_{1}}(t).$
(18)

The evolution of the reduced state ($|00\rangle\langle 00|$) by tracing over
bath modes can be found as,

$\displaystyle\phi(|00\rangle\langle 00|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\chi(t)\rangle\langle\chi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{\\{}|D_{1}(m,n,t)|^{2}|00\rangle\langle
00|$ $\displaystyle+n|E_{1}(m,n,t)|^{2}|01\rangle\langle 01|$
$\displaystyle+m|F_{1}(m,n,t)|^{2}|10\rangle\langle 10|\Big{\\}}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}]$
(19)

One may find the components $D_{1}(m,n,t)$, $E_{1}(m,n,t)$ and $F_{1}(m,n,t)$,
which are the eigenvalues of the operators, $\hat{D_{1}}(t)$, $\hat{E_{1}}(t)$
and $\hat{F_{1}}(t)$ with eigenvectors $|mn\rangle$, by following the steps
described in the Appendix B.

In a similar manner, we specify $|\xi(0)\rangle=|01\rangle|xy\rangle$ and
$|\xi(t)\rangle=\hat{G}(t)|01\rangle|xy\rangle+\hat{H}(t)|00\rangle|xy\rangle+\hat{I}(t)|11\rangle|xy\rangle$.
Now substituting with, $\hat{G}(t)=\hat{G_{1}}(t)$,
$\hat{H}(t)=b^{\dagger}\hat{H_{1}}(t)$ and $\hat{I}(t)=a\hat{I_{1}}(t)$, we
can have,

$\displaystyle\frac{d\hat{G_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{-\omega_{1}+\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{G_{1}}(t)$
$\displaystyle-i\epsilon_{2}\Big{(}1-\frac{\hat{n}}{N}\Big{)}^{1/2}(\hat{n}+1)\hat{H_{1}}(t)$
$\displaystyle-i\epsilon_{1}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}^{1/2}\hat{I_{1}}(t),$
$\displaystyle\frac{d\hat{H_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{-\omega_{1}-\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(\hat{n}+1)\Big{(}1-\frac{\hat{n}}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{H_{1}}(t)$
$\displaystyle-i\epsilon_{2}\Big{(}1-\frac{\hat{n}}{N}\Big{)}^{1/2}\hat{G_{1}}(t),$
$\displaystyle\frac{d\hat{I_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}+\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}(\hat{m}-1)\Big{(}1-\frac{\hat{m}-2}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{I_{1}}(t)$
$\displaystyle-i\epsilon_{1}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}^{1/2}\hat{G_{1}}(t).$
(20)

We may now express the evolution of the reduced state ($|01\rangle\langle
01|$) in the following way,

$\displaystyle\phi(|01\rangle\langle 01|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\xi(t)\rangle\langle\xi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{\\{}|G_{1}(m,n,t)|^{2}|01\rangle\langle
01|$ $\displaystyle+(n+1)|H_{1}(m,n,t)|^{2}|00\rangle\langle 00|$
$\displaystyle+m|I_{1}(m,n,t)|^{2}|11\rangle\langle 11|\Big{\\}}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}]$
(21)

Following the procedure written in the Appendix B, one may calculate the
eigenvalues of $\hat{G_{1}}(t)$, $\hat{H_{1}}(t)$ and $\hat{I_{1}}(t)$ i.e.
$G_{1}(m,n,t)$, $H_{1}(m,n,t)$ and $I_{1}(m,n,t)$ for eigenvectors
$|mn\rangle$.

To determine the evolution of the fourth diagonal element ($|10\rangle\langle
10|$) of the density matrix of the system of two central spins, we assume,
$|\varrho(0)\rangle=|10\rangle|xy\rangle$ and
$|\varrho(t)\rangle=\hat{J}(t)|10\rangle|xy\rangle+\hat{K}(t)|11\rangle|xy\rangle+\hat{L}(t)|00\rangle|xy\rangle$.
Then we replace the time-dependent operators as, $\hat{J}(t)=\hat{J_{1}}(t)$,
$\hat{K}(t)=b\hat{K_{1}}(t)$ and $\hat{L}(t)=a^{\dagger}\hat{L_{1}}(t)$ and we
obtain,

$\displaystyle\frac{d\hat{J_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}-\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{J_{1}}(t)$
$\displaystyle-i\epsilon_{2}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}^{1/2}\hat{K_{1}}(t)$
$\displaystyle-i\epsilon_{1}\Big{(}1-\frac{\hat{m}}{M}\Big{)}^{1/2}(\hat{m}+1)\hat{L_{1}}(t),$
$\displaystyle\frac{d\hat{K_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{\omega_{1}+\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}\hat{m}\Big{(}1-\frac{\hat{m}-1}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(\hat{n}-1)\Big{(}1-\frac{\hat{n}-2}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{K_{1}}(t)$
$\displaystyle-i\epsilon_{2}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}^{1/2}\hat{J_{1}}(t),$
$\displaystyle\frac{d\hat{L_{1}}(t)}{dt}=$
$\displaystyle-i\Big{[}\frac{-\omega_{1}-\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}(\hat{m}+1)\Big{(}1-\frac{\hat{m}}{M}\Big{)}-\frac{\openone}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}\hat{n}\Big{(}1-\frac{\hat{n}-1}{N}\Big{)}-\frac{\openone}{2}\Big{\\}}\Big{]}\hat{L_{1}}(t)$
$\displaystyle-i\epsilon_{1}\Big{(}1-\frac{\hat{m}}{M}\Big{)}^{1/2}\hat{J_{1}}(t).$
(22)

The reduced state now changes as,

$\displaystyle\phi(|10\rangle\langle 10|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\varrho(t)\rangle\langle\varrho(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{\\{}|J_{1}(m,n,t)|^{2}|10\rangle\langle
10|$ $\displaystyle+n|K_{1}(m,n,t)|^{2}|11\rangle\langle 11|$
$\displaystyle+(m+1)|L_{1}(m,n,t)|^{2}|00\rangle\langle 00|\Big{\\}}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(23)

where, by definition, $\hat{J_{1}}(t)|mn\rangle=J_{1}(m,n,t)|mn\rangle$,
$\hat{K_{1}}(t)|mn\rangle=K_{1}(m,n,t)|mn\rangle$ and
$\hat{L_{1}}(t)|mn\rangle=L_{1}(m,n,t)|mn\rangle$. The corresponding
eigenvalues can be found from Eq. (22) by following the approach mentioned in
the Appendix B.

Now the off-diagonal components of the reduced density matrix will take the
following forms:

$\displaystyle\phi(|11\rangle\langle 00|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\psi(t)\rangle\langle\chi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}A_{1}(m,n,t)D_{1}^{*}(m,n,t)|11\rangle\langle
00|\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(24) $\displaystyle\phi(|11\rangle\langle 01|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\psi(t)\rangle\langle\xi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}A_{1}(m,n,t)G_{1}^{*}(m,n,t)|11\rangle\langle
01|\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(25) $\displaystyle\phi(|11\rangle\langle 10|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\psi(t)\rangle\langle\varrho(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}A_{1}(m,n,t)J_{1}^{*}(m,n,t)|11\rangle\langle
10|\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(26) $\displaystyle\phi(|10\rangle\langle 00|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\varrho(t)\rangle\langle\chi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}J_{1}(m,n,t)D_{1}^{*}(m,n,t)|10\rangle\langle
00|\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(27) $\displaystyle\phi(|10\rangle\langle 01|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\varrho(t)\rangle\langle\xi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}J_{1}(m,n,t)G_{1}^{*}(m,n,t)|10\rangle\langle
01|\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(28) $\displaystyle\phi(|01\rangle\langle 00|)$
$\displaystyle=\operatorname{Tr}_{B_{1}B_{2}}[|\xi(t)\rangle\langle\chi(t)|]$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}G_{1}(m,n,t)D_{1}^{*}(m,n,t)|01\rangle\langle
00|\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}].$
(29)

From the hermiticity of the reduced density matrix, the other off-diagonal
elements can be expressed as,

$\displaystyle\phi(|00\rangle\langle 11|)=\Big{(}\phi(|11\rangle\langle
00|)\Big{)}^{\dagger},\phi(|01\rangle\langle
11|)=\Big{(}\phi(|11\rangle\langle 01|)\Big{)}^{\dagger},$
$\displaystyle\phi(|10\rangle\langle 11|)=\Big{(}\phi(|11\rangle\langle
10|)\Big{)}^{\dagger},\phi(|00\rangle\langle
10|)=\Big{(}\phi(|10\rangle\langle 00|)\Big{)}^{\dagger},$
$\displaystyle\phi(|01\rangle\langle 10|)=\Big{(}\phi(|10\rangle\langle
01|)\Big{)}^{\dagger},\phi(|00\rangle\langle
01|)=\Big{(}\phi(|01\rangle\langle 00|)\Big{)}^{\dagger}.$

##  Appendix B Solution of simultaneous linear differential equations from
boundary conditions: Matrix method

Suppose, three simultaneous linear differential equations are written as:

$\displaystyle\frac{dX(t)}{dt}=-iaX(t)-idY(t)-ieZ(t),$
$\displaystyle\frac{dY(t)}{dt}=-ifX(t)-ibY(t),$
$\displaystyle\frac{dZ(t)}{dt}=-igX(t)-icZ(t).$ (30)

And we have to find the solution using the given boundary conditions:
$X(0)=1$, $Y(0)=0$ and $Z(0)=0$.

The problem can be solved by plugging in the well-known Matrix method in this
context. According to this, one can write Eq. (30) using column and square
matrices of dimension 3 in a following way,

$\displaystyle\frac{d}{dt}\begin{pmatrix}X(t)\\\ Y(t)\\\ Z(t)\\\
\end{pmatrix}=-i\begin{pmatrix}a&d&e\\\ f&b&0\\\ g&0&c\\\
\end{pmatrix}\begin{pmatrix}X(t)\\\ Y(t)\\\ Z(t)\\\ \end{pmatrix}$ (31)

For simplicity, we denote, $M$=$\begin{pmatrix}a&d&e\\\ f&b&0\\\ g&0&c\\\
\end{pmatrix}$. Any of the three eigenvalues, say $\lambda_{j}$($j$=1,2,3) of
matrix, $M$ can be found by solving the characteristic Eq,

$\displaystyle(a-\lambda_{j})(b-\lambda_{j})(c-\lambda_{j})-fd(c-\lambda_{j})-eg(b-\lambda_{j})=0,~{}\forall
j.$ (32)

And corresponding to each $\lambda_{j}$, eigenvectors can be expressed as,
$\begin{pmatrix}(\lambda_{j}-b)(\lambda_{j}-c)\\\ f(\lambda_{j}-c)\\\
g(\lambda_{j}-b)\\\ \end{pmatrix}\equiv\begin{pmatrix}\alpha_{j}\\\
\beta_{j}\\\ \gamma_{j}\\\ \end{pmatrix}$, where,
$\alpha_{j}=(\lambda_{j}-b)(\lambda_{j}-c)$, $\beta_{j}=f(\lambda_{j}-c)$ and
$\gamma_{j}=g(\lambda_{j}-b)$, except for the case when $\lambda_{j}=b,\;c$.
Then according to the Matrix method, the general solution of Eq. (30) in terms
of column matrix will take the form,

$\displaystyle\begin{pmatrix}X(t)\\\ Y(t)\\\ Z(t)\\\
\end{pmatrix}=\sum_{j=1}^{3}v_{j}~{}e^{-i\lambda_{j}t}~{}\begin{pmatrix}\alpha_{j}\\\
\beta_{j}\\\ \gamma_{j}\\\ \end{pmatrix},$ (33)

where the co-efficients $v_{j}$ can be computed from three boundary conditions
given as,

$\displaystyle\sum_{j=1}^{3}v_{j}\alpha_{j}=1;~{}~{}\sum_{j=1}^{3}v_{j}\beta_{j}=0;~{}~{}\sum_{j=1}^{3}v_{j}\gamma_{j}=0.$
(34)

Now one can write from Eq. (34),

$\displaystyle
v_{1}=\frac{\beta_{3}\gamma_{2}-\beta_{2}\gamma_{3}}{\alpha_{1}(\beta_{3}\gamma_{2}-\beta_{2}\gamma_{3})+\alpha_{2}(\beta_{1}\gamma_{3}-\beta_{3}\gamma_{1})+\alpha_{3}(\beta_{2}\gamma_{1}-\beta_{1}\gamma_{2})},$
$\displaystyle
v_{2}=\frac{\beta_{1}\gamma_{3}-\beta_{3}\gamma_{1}}{\alpha_{1}(\beta_{3}\gamma_{2}-\beta_{2}\gamma_{3})+\alpha_{2}(\beta_{1}\gamma_{3}-\beta_{3}\gamma_{1})+\alpha_{3}(\beta_{2}\gamma_{1}-\beta_{1}\gamma_{2})},$
$\displaystyle
v_{3}=\frac{\beta_{2}\gamma_{1}-\beta_{1}\gamma_{2}}{\alpha_{1}(\beta_{3}\gamma_{2}-\beta_{2}\gamma_{3})+\alpha_{2}(\beta_{1}\gamma_{3}-\beta_{3}\gamma_{1})+\alpha_{3}(\beta_{2}\gamma_{1}-\beta_{1}\gamma_{2})}.$
(35)

We see that, we are able to come up with the solution (33) only when
$\lambda_{j}$ is calculated from the Eq. (32). For solving the three roots of
such cubic equation, we follow the procedure given below.

Let us say, a cubic equation is in the form:

$\displaystyle a_{1}x^{3}+b_{1}x^{2}+c_{1}x+d_{1}=0$ (36)

Then we define the following quantities,

$\displaystyle\Delta=18a_{1}b_{1}c_{1}d_{1}-4b_{1}^{3}d_{1}+b_{1}^{2}c_{1}^{2}-4a_{1}c_{1}^{3}-27a_{1}^{2}d_{1}^{2},$
$\displaystyle\Delta_{0}=b_{1}^{2}-3a_{1}c_{1},$
$\displaystyle\Delta_{1}=2b_{1}^{3}-9a_{1}b_{1}c_{1}+27a_{1}^{2}d_{1},$
$\displaystyle
Q=\sqrt[3]{\frac{\Delta_{1}+\sqrt{\Delta_{1}^{2}-4\Delta_{0}^{3}}}{2}}=\sqrt[3]{\frac{\Delta_{1}+\sqrt{-27a_{1}^{2}\Delta}}{2}}.$

where, $\Delta$ is called the discriminant of the cubic Eq. (36). Next,
denoting a complex number, $-\frac{1}{2}+\frac{\sqrt{3}}{2}i=\varphi$, we may
write the roots of Eq. (36) as,

$\displaystyle
x_{k}=-\frac{1}{3a_{1}}(b_{1}+\varphi^{k}Q+\frac{\Delta_{0}}{\varphi^{k}Q}),~{}k\in\\{1,2,3\\}.$
(37)

Note that, the root $x_{1}$ is always real and the roots $x_{2}$ and $x_{3}$
are complex and conjugate to each other only when $\Delta<0$. Otherwise, all
the roots are real and they become equal when $\Delta=0$.

In our case, expanding the Left Hand Side of Eq. (32) in a power series of
$\lambda_{j}$ for all $j$, we identify, $a_{1}=1$, $b_{1}=-(a+b+c)$,
$c_{1}=ab+bc+ac-fd-eg$ and $d_{1}=fdc+egb-abc$. As described above,
$\lambda_{j}$’s($j$=1,2,3) can be determined and the solution (33) can be
attained accordingly.

All the Eqs. (9,18,20,22) are in the form of Eq. (30). Hence, all the co-
efficients can be obtained in terms of $\omega_{1}$, $\omega_{2}$, $\delta$,
$\omega_{a}$, $\omega_{b}$, $\epsilon_{1}$, $\epsilon_{2}$, m, n, M and N.

Explicitly, the co-efficients, $A_{1}(m,n,t)$, $B_{1}(m,n,t)$ and
$C_{1}(m,n,t)$ can be obtained by solving Eq. (9) which in comparison to Eq.
(30) shows that,

$\displaystyle X(t)$ $\displaystyle=A_{1}(m,n,t),$ $\displaystyle Y(t)$
$\displaystyle=B_{1}(m,n,t),$ $\displaystyle Z(t)$
$\displaystyle=C_{1}(m,n,t),$ $\displaystyle a=$
$\displaystyle\frac{\omega_{1}+\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle b=$
$\displaystyle\frac{\omega_{1}-\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(n+1)\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle c=$
$\displaystyle\frac{-\omega_{1}+\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}(m+1)\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle d=$
$\displaystyle\epsilon_{2}\Big{(}1-\frac{n}{N}\Big{)}^{1/2}(n+1),$
$\displaystyle e=$
$\displaystyle\epsilon_{1}\Big{(}1-\frac{m}{M}\Big{)}^{1/2}(m+1),$
$\displaystyle f=$
$\displaystyle\epsilon_{2}\Big{(}1-\frac{n}{N}\Big{)}^{1/2},$ $\displaystyle
g=$ $\displaystyle\epsilon_{1}\Big{(}1-\frac{m}{M}\Big{)}^{1/2}.$

Likewise, comparing Eq. (18) with Eq. (30), we get,

$\displaystyle X(t)$ $\displaystyle=D_{1}(m,n,t),$ $\displaystyle Y(t)$
$\displaystyle=E_{1}(m,n,t),$ $\displaystyle Z(t)$
$\displaystyle=F_{1}(m,n,t),$ $\displaystyle a=$
$\displaystyle\frac{-\omega_{1}-\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle b=$
$\displaystyle\frac{-\omega_{1}+\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(n-1)\Big{(}1-\frac{n-2}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle c=$
$\displaystyle\frac{\omega_{1}-\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}(m-1)\Big{(}1-\frac{m-2}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle d=$
$\displaystyle\epsilon_{2}n\Big{(}1-\frac{n-1}{N}\Big{)}^{1/2},$
$\displaystyle e=$
$\displaystyle\epsilon_{1}m\Big{(}1-\frac{m-1}{M}\Big{)}^{1/2},$
$\displaystyle f=$
$\displaystyle\epsilon_{2}\Big{(}1-\frac{n-1}{N}\Big{)}^{1/2},$ $\displaystyle
g=$ $\displaystyle\epsilon_{1}\Big{(}1-\frac{m-1}{M}\Big{)}^{1/2}.$

Correspondingly, Eq. (20) exhibits,

$\displaystyle X(t)$ $\displaystyle=G_{1}(m,n,t),$ $\displaystyle Y(t)$
$\displaystyle=H_{1}(m,n,t),$ $\displaystyle Z(t)$
$\displaystyle=I_{1}(m,n,t),$ $\displaystyle a=$
$\displaystyle\frac{-\omega_{1}+\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle b=$
$\displaystyle\frac{-\omega_{1}-\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(n+1)\Big{(}1-\frac{n}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle c=$
$\displaystyle\frac{\omega_{1}+\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}(m-1)\Big{(}1-\frac{m-2}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle d=$
$\displaystyle\epsilon_{2}\Big{(}1-\frac{n}{N}\Big{)}^{1/2}(n+1),$
$\displaystyle e=$
$\displaystyle\epsilon_{1}m\Big{(}1-\frac{m-1}{M}\Big{)}^{1/2},$
$\displaystyle f=$
$\displaystyle\epsilon_{2}\Big{(}1-\frac{n}{N}\Big{)}^{1/2},$ $\displaystyle
g=$ $\displaystyle\epsilon_{1}\Big{(}1-\frac{m-1}{M}\Big{)}^{1/2}.$

And finally the resemblance between Eq. (22) and Eq. (30) manifests,

$\displaystyle X(t)$ $\displaystyle=J_{1}(m,n,t),$ $\displaystyle Y(t)$
$\displaystyle=K_{1}(m,n,t),$ $\displaystyle Z(t)$
$\displaystyle=L_{1}(m,n,t),$ $\displaystyle a=$
$\displaystyle\frac{\omega_{1}-\omega_{2}-\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle b=$
$\displaystyle\frac{\omega_{1}+\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}(n-1)\Big{(}1-\frac{n-2}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle c=$
$\displaystyle\frac{-\omega_{1}-\omega_{2}+\delta}{2}+\omega_{a}\Big{\\{}(m+1)\Big{(}1-\frac{m}{M}\Big{)}-\frac{1}{2}\Big{\\}}$
$\displaystyle+\omega_{b}\Big{\\{}n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\Big{\\}},$
$\displaystyle d=$
$\displaystyle\epsilon_{2}n\Big{(}1-\frac{n-1}{N}\Big{)}^{1/2},$
$\displaystyle e=$
$\displaystyle\epsilon_{1}\Big{(}1-\frac{m}{M}\Big{)}^{1/2}(m+1),$
$\displaystyle f=$
$\displaystyle\epsilon_{2}\Big{(}1-\frac{n-1}{N}\Big{)}^{1/2},$ $\displaystyle
g=$ $\displaystyle\epsilon_{1}\Big{(}1-\frac{m}{M}\Big{)}^{1/2}.$

In this way, all the co-efficients can be obtained in order to reveal the
complete dynamical map of the two-spin system.

##  Appendix C Matrix elements of the reduced state of the system

The elements of the density matrix given in Eq. (11) are given by,

$\displaystyle\rho_{11}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}|A_{1}(m,n,t)|^{2}\rho_{11}(0)+n|K_{1}(m,n,t)|^{2}\rho_{22}(0)$
$\displaystyle+m|I_{1}(m,n,t)|^{2}\rho_{33}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(38) $\displaystyle\rho_{22}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}[(n+1)\Big{(}B_{1}(m,n,t)|^{2}\rho_{11}(0)+|J_{1}(m,n,t)|^{2}\rho_{22}(0)$
$\displaystyle+m|F_{1}(m,n,t)|^{2}\rho_{44}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(39) $\displaystyle\rho_{33}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}(m+1)|C_{1}(m,n,t)|^{2}\rho_{11}(0)+|G_{1}(m,n,t)|^{2}\rho_{33}(0)$
$\displaystyle+n|E_{1}(m,n,t)|^{2}\rho_{44}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(40) $\displaystyle\rho_{44}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}(m+1)|B_{1}(m,n,t)|^{2}\rho_{22}(0)+|D_{1}(m,n,t)|^{2}\rho_{44}(0)$
$\displaystyle+(n+1)|H_{1}(m,n,t)|^{2}\rho_{33}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(41) $\displaystyle\rho_{12}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}A_{1}(m,n,t)J_{1}^{*}(m,n,t)\rho_{12}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(42) $\displaystyle\rho_{13}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}A_{1}(m,n,t)G_{1}^{*}(m,n,t)\rho_{13}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(43) $\displaystyle\rho_{14}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}A_{1}(m,n,t)D_{1}^{*}(m,n,t)\rho_{14}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(44) $\displaystyle\rho_{23}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}J_{1}(m,n,t)G_{1}^{*}(m,n,t)\rho_{23}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(45) $\displaystyle\rho_{24}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}J_{1}(m,n,t)D_{1}^{*}(m,n,t)\rho_{24}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(46) $\displaystyle\rho_{34}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}\Big{(}G_{1}(m,n,t)D_{1}^{*}(m,n,t)\rho_{34}(0)\Big{)}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}].$
(47)

##  Appendix D Elements of the Choi-Jamiołkowski matrix

We have used following abbreviations in the CJ matrix given in Eq. (12),

$\displaystyle|A_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}|A_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(48) $\displaystyle(n+1)|B_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}(n+1)|B_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(49) $\displaystyle(m+1)|C_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}(m+1)|C_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(50) $\displaystyle n|K_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}n|K_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(51) $\displaystyle|J_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}|J_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(52) $\displaystyle(m+1)|L_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}(m+1)|L_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(53) $\displaystyle m|I_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}m|I_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(54) $\displaystyle|G_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}|G_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(55) $\displaystyle(n+1)|H_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}(n+1)|H_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(56) $\displaystyle m|F_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}m|F_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(57) $\displaystyle n|E_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}n|E_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(58) $\displaystyle|D_{1}|^{2}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}|D_{1}(m,n,t)|^{2}$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(59) $\displaystyle A_{1}J_{1}^{*}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}A_{1}(m,n,t)J_{1}^{*}(m,n,t)$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(60) $\displaystyle A_{1}G_{1}^{*}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}A_{1}(m,n,t)G_{1}^{*}(m,n,t)$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(61) $\displaystyle A_{1}D_{1}^{*}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}A_{1}(m,n,t)D_{1}^{*}(m,n,t)$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(62) $\displaystyle J_{1}G_{1}^{*}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}J_{1}(m,n,t)G_{1}^{*}(m,n,t)$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(63) $\displaystyle J_{1}D_{1}^{*}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}J_{1}(m,n,t)D_{1}^{*}(m,n,t)$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}],$
(64) $\displaystyle G_{1}D_{1}^{*}$
$\displaystyle=\frac{1}{Z}\sum_{m}^{M}\sum_{n}^{N}G_{1}(m,n,t)D_{1}^{*}(m,n,t)$
$\displaystyle\times\exp[-\frac{\hslash\omega_{a}}{K_{B}T}\\{m\Big{(}1-\frac{m-1}{M}\Big{)}-\frac{1}{2}\\}]$
$\displaystyle\times\exp[-\frac{\hslash\omega_{b}}{K_{B}T}\\{n\Big{(}1-\frac{n-1}{N}\Big{)}-\frac{1}{2}\\}].$
(65)

The rest of the elements, that is, lower diagonal elements in Eq. (12) are
hermitian conjugates of the upper diagonal terms.

## References

  * Cirac and Zoller (1995) J. I. Cirac and P. Zoller, Phys. Rev. Lett. 74, 4091 (1995).
  * Loss and DiVincenzo (1998) D. Loss and D. P. DiVincenzo, Phys. Rev. A 57, 120 (1998).
  * Cory  _et al._ (1997) D. G. Cory, A. F. Fahmy, and T. F. Havel, Proceedings of the National Academy of Sciences 94, 1634 (1997), https://www.pnas.org/doi/pdf/10.1073/pnas.94.5.1634 .
  * Rauch (1998) H. Rauch, Physica Scripta T76, 24 (1998).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Great Clarendon Street, 2002).
  * Banerjee (2018) S. Banerjee,  _Open Quantum Systems_ (Springer Singapore, 152 Beach Rd, 21-01 Gateway East, Singapore 189721, 2018).
  * Leggett  _et al._ (1987a) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987a).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Lindblad (1976) G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Caldeira and Leggett (1983) A. O. Caldeira and A. J. Leggett, Annals of Physics 149, 374 (1983).
  * Feynman and Vernon (1963) R. Feynman and F. Vernon, Annals of Physics 24, 118 (1963).
  * Weiss (1999) U. Weiss,  _Quantum Dissipative Systems_ , 2nd ed., Series in Modern Condensed Matter Physics, Vol. 10 (World Scientific, 1999).
  * Leggett  _et al._ (1987b) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987b).
  * Nakajima (1958) S. Nakajima, Progress of Theoretical Physics 20, 948 (1958).
  * Zwanzig (1960) R. Zwanzig, The Journal of Chemical Physics 33, 1338 (1960).
  * Chaturvedi and Shibata (1979) S. Chaturvedi and F. Shibata, Zeitschrift für Physik B Condensed Matter 35, 297 (1979).
  * Banerjee and Ghosh (2000) S. Banerjee and R. Ghosh, Phys. Rev. A 62, 042105 (2000).
  * Banerjee and Ghosh (2003a) S. Banerjee and R. Ghosh, Phys. Rev. E 67, 056120 (2003a).
  * Chandrashekar  _et al._ (2007) C. M. Chandrashekar, R. Srikanth, and S. Banerjee, Phys. Rev. A 76, 022316 (2007).
  * Banerjee  _et al._ (2008) S. Banerjee, R. Srikanth, C. M. Chandrashekar, and P. Rungta, Phys. Rev. A 78, 052316 (2008).
  * Thomas  _et al._ (2018) G. Thomas, N. Siddharth, S. Banerjee, and S. Ghosh, Phys. Rev. E 97, 062108 (2018).
  * Naikoo  _et al._ (2019) J. Naikoo, A. K. Alok, S. Banerjee, and S. U. Sankar, Phys. Rev. D 99, 095001 (2019).
  * Alok  _et al._ (2016) A. K. Alok, S. Banerjee, and S. Uma Sankar, Nuclear Physics B 909, 65 (2016).
  * Naikoo  _et al._ (2018) J. Naikoo, A. K. Alok, and S. Banerjee, Phys. Rev. D 97, 053008 (2018).
  * Shrikant  _et al._ (2018) U. Shrikant, R. Srikanth, and S. Banerjee, Phys. Rev. A 98, 032328 (2018).
  * Dixit  _et al._ (2019) K. Dixit, J. Naikoo, S. Banerjee, and A. K. Alok, The European Physical Journal C 79, 1 (2019).
  * Breuer  _et al._ (2004) H.-P. Breuer, D. Burgarth, and F. Petruccione, Phys. Rev. B 70, 045323 (2004).
  * Laine  _et al._ (2010) E.-M. Laine, J. Piilo, and H.-P. Breuer, Phys. Rev. A 81, 062115 (2010).
  * Rivas  _et al._ (2010) A. Rivas, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050403 (2010).
  * Chruściński  _et al._ (2011) D. Chruściński, A. Kossakowski, and A. Rivas, Phys. Rev. A 83, 052128 (2011).
  * Vasile  _et al._ (2011) R. Vasile, S. Maniscalco, M. G. A. Paris, H.-P. Breuer, and J. Piilo, Phys. Rev. A 84, 052118 (2011).
  * Lu  _et al._ (2010) X.-M. Lu, X. Wang, and C. P. Sun, Phys. Rev. A 82, 042103 (2010).
  * Luo  _et al._ (2012) S. Luo, S. Fu, and H. Song, Phys. Rev. A 86, 044101 (2012).
  * Fanchini (2014) F. F. e. Fanchini, Phys. Rev. Lett. 112, 210402 (2014).
  * Chanda and Bhattacharya (2016) T. Chanda and S. Bhattacharya, Annals of Physics 366, 1 (2016).
  * Haseli (2014) S. e. Haseli, Phys. Rev. A 90, 052118 (2014).
  * Hall  _et al._ (2014) M. J. W. Hall, J. D. Cresser, L. Li, and E. Andersson, Phys. Rev. A 89, 042120 (2014).
  * Andersson  _et al._ (2007) E. Andersson, J. D. Cresser, and M. J. W. Hall, Journal of Modern Optics 54, 1695 (2007), https://doi.org/10.1080/09500340701352581 .
  * Bhattacharya  _et al._ (2017) S. Bhattacharya, A. Misra, C. Mukhopadhyay, and A. K. Pati, Phys. Rev. A 95, 012122 (2017).
  * Mukhopadhyay  _et al._ (2017) C. Mukhopadhyay, S. Bhattacharya, A. Misra, and A. K. Pati, Phys. Rev. A 96, 052125 (2017).
  * Awasthi  _et al._ (2018) N. Awasthi, S. Bhattacharya, A. Sen(De), and U. Sen, Phys. Rev. A 97, 032103 (2018).
  * Bhattacharya  _et al._ (2020a) S. Bhattacharya, B. Bhattacharya, and A. S. Majumdar, Journal of Physics A: Mathematical and Theoretical 54, 035302 (2020a).
  * Bhattacharya  _et al._ (2020b) S. Bhattacharya, B. Bhattacharya, and A. S. Majumdar, Journal of Physics A: Mathematical and Theoretical 53, 335301 (2020b).
  * Maity  _et al._ (2020) A. G. Maity, S. Bhattacharya, and A. S. Majumdar, Journal of Physics A: Mathematical and Theoretical 53, 175301 (2020).
  * Bhattacharya and Bhattacharya (2021) B. Bhattacharya and S. Bhattacharya, Quantum Information Processing 20, 253 (2021).
  * Bhattacharya  _et al._ (2018) S. Bhattacharya, S. Banerjee, and A. K. Pati, Quantum Information Processing 17, 236 (2018).
  * Das  _et al._ (2021) S. Das, S. S. Roy, S. Bhattacharya, and U. Sen, Journal of Physics A: Mathematical and Theoretical 54, 395301 (2021).
  * Banerjee and Ghosh (2003b) S. Banerjee and R. Ghosh, Phys. Rev. E 67, 056120 (2003b).
  * Kumar  _et al._ (2018a) N. P. Kumar, S. Banerjee, and C. M. Chandrashekar, Scientific reports 8, 1 (2018a).
  * Kumar  _et al._ (2018b) N. P. Kumar, S. Banerjee, R. Srikanth, V. Jagadish, and F. Petruccione, Open Systems & Information Dynamics 25, 1850014 (2018b).
  * Utagi  _et al._ (2020) S. Utagi, R. Srikanth, and S. Banerjee, Scientific Reports 10, 1 (2020).
  * Krovi  _et al._ (2007) H. Krovi, O. Oreshkov, M. Ryazanov, and D. A. Lidar, Phys. Rev. A 76, 052117 (2007).
  * Breuer  _et al._ (2009) H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. 103, 210401 (2009).
  * de Vega and Alonso (2017) I. de Vega and D. Alonso, Rev. Mod. Phys. 89, 015001 (2017).
  * Garraway and Knight (1996) B. M. Garraway and P. L. Knight, Phys. Rev. A 54, 3592 (1996).
  * Dalton and Garraway (2003) B. J. Dalton and B. M. Garraway, Phys. Rev. A 68, 033809 (2003).
  * Iles-Smith  _et al._ (2014) J. Iles-Smith, N. Lambert, and A. Nazir, Phys. Rev. A 90, 032114 (2014).
  * Garg  _et al._ (1985) A. Garg, J. N. Onuchic, and V. Ambegaokar, The Journal of Chemical Physics 83, 4491 (1985), https://doi.org/10.1063/1.449017 .
  * Hughes  _et al._ (2009) K. H. Hughes, C. D. Christ, and I. Burghardt, The Journal of Chemical Physics 131, 124108 (2009), https://doi.org/10.1063/1.3226343 .
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, Journal of the Physical Society of Japan 58, 101 (1989), https://doi.org/10.1143/JPSJ.58.101 .
  * Fischer and Breuer (2007) J. Fischer and H.-P. Breuer, Phys. Rev. A 76, 052119 (2007).
  * Ordonez-Miranda  _et al._ (2017) J. Ordonez-Miranda, Y. Ezzahri, and K. Joulain, Phys. Rev. E 95, 022128 (2017).
  * Linden  _et al._ (2010) N. Linden, S. Popescu, and P. Skrzypczyk, Phys. Rev. Lett. 105, 130401 (2010).
  * Wang  _et al._ (2022a) L. Wang, Z. Wang, C. Wang, and J. Ren, Phys. Rev. Lett. 128, 067701 (2022a).
  * Ar ısoy  _et al._ (2022) O. Ar ısoy, J.-T. Hsiang, and B.-L. Hu, Phys. Rev. E 105, 014108 (2022).
  * Mohanta  _et al._ (2022) S. Mohanta, S. Saryal, and B. K. Agarwalla, Phys. Rev. E 105, 034127 (2022).
  * Karg ı  _et al._ (2019) C. Karg ı, M. T. Naseem, T. c. v. Opatrný, O. E. Müstecaplıoğlu, and G. Kurizki, Phys. Rev. E 99, 042121 (2019).
  * Joulain  _et al._ (2016) K. Joulain, J. Drevillon, Y. Ezzahri, and J. Ordonez-Miranda, Phys. Rev. Lett. 116, 200601 (2016).
  * Wijesekara  _et al._ (2021) R. T. Wijesekara, S. D. Gunapala, and M. Premaratne, Phys. Rev. B 104, 045405 (2021).
  * Majland  _et al._ (2020) M. Majland, K. S. Christensen, and N. T. Zinner, Phys. Rev. B 101, 184510 (2020).
  * Ghosh  _et al._ (2021) R. Ghosh, A. Ghoshal, and U. Sen, Phys. Rev. A 103, 052613 (2021).
  * Wang  _et al._ (2022b) L. Wang, Z. Wang, C. Wang, and J. Ren, Phys. Rev. Lett. 128, 067701 (2022b).
  * Bhandari  _et al._ (2021) B. Bhandari, P. A. Erdman, R. Fazio, E. Paladino, and F. Taddei, Phys. Rev. B 103, 155434 (2021).
  * Bhattacharya and Banerjee (2021) S. Bhattacharya and S. Banerjee, Quanta 10, 55 (2021).
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
  * Banerjee  _et al._ (2010) S. Banerjee, V. Ravishankar, and R. Srikanth, Annals of Physics 325, 816 (2010).
  * Wootters (1998) W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
  * Ollivier and Zurek (2001) H. Ollivier and W. H. Zurek, Phys. Rev. Lett. 88, 017901 (2001).