\begin{equation}gin{document}

\thetaitle{Decoherence in strongly coupled quantum oscillators}
\author{M. A. de Ponte$^{1}$, M. C. de Oliveira$^{2}$, and M. H. Y. Moussa$
^{3}$}
\address{Departamento de F\'{\i}sica, CCET, Universidade Federal de S\~{a}o\\
Carlos, Via Washington Luiz Km 235, S\~{a}o Carlos, 13565-905, SP, Brazil.}
\title{Decoherence in strongly coupled quantum oscillators}

\begin{equation}gin{abstract}
In this paper we present a comprehensive analysis of the coherence
phenomenon of two coupled dissipative oscillators. The action of a classical
driving field on one of the oscillators is also analyzed. Master equations
are derived for both regimes of weakly and strongly interacting oscillators
from which interesting results arise concerning the coherence properties of
the joint and the reduced system states. The strong coupling regime is
required to achieve a large frequency shift of the oscillator normal modes,
making it possible to explore the whole profile of the spectral density of
the reservoirs. We show how the decoherence process may be controlled by
shifting the normal mode frequencies to regions of small spectral density of
the reservoirs. Different spectral densities of the reservoirs are
considered and their effects on the decoherence process are analyzed. For
oscillators with different damping rates, we show that the worse-quality
system is improved and vice-versa, a result which could be useful for
quantum state protection. State recurrence and swap dynamics are analyzed as
well as their roles in delaying the decoherence process.
\end{abstract}

\pacs{PACS number: 42.50.Ct, 42.50.Dv, 03.65.Bz, 32.80.-t}

\section{Introduction}

The process of decoherence of quantum states has long been a central issue
in the description of quantum measurements \cite{Neumann,Zurek,CL}. In
recent years, experimental advances in the domain of cavity QED and trapped
ions have allowed the decoherence of photon \cite{Haroche} and phonon \cite
{Wineland} field states to be probed in more depth, providing insights into
the borderline between the classical and quantum descriptions of the
physical world. The decoherence time of a superposition of coherent states
in a cavity field was measured \cite{Haroche} and shown to be in full
agreement with theoretical predictions \cite{Glauber,WM}. In trapped ions
systems, the observed damping of Rabi oscillations has motivated a number of
articles on the main sources of noise leading to decoherence \cite
{Milburn,Vogel,Serra,Budini}. Such experimental achievements in matter-field
interactions have also encouraged a deep dialog between theoretical and
experimental physics, resulting in a degree of mastery of fundamental
quantum phenomena that may herald a new stage in the technology of
communication \cite{Comm} and computation \cite{Comp}.

The exploration of the borderline between quantum and classical descriptions
of nature \cite{Zurek} has impelled the generation of superposition states
of mesoscopic systems, known as \thetaextquotedblleft Schr\"{o}dinger cat
states\thetaextquotedblright\ \cite{Haroche,Wineland}. Such superpositions are
irreversibly affected by their surroundings, whose effect is to destroy
probability interference (coherence), and continuously transformed into
statistical mixtures. Thus, the environment plays a key role in the
establishment of a direct correspondence between quantum and classical
dynamics. While the decoherence time of a superposition state depends on the
amplitude of the field, the relaxation does not, since the model adopted for
the relaxation process is amplitude damping, achieved by coupling the
systems bilinearly to the degrees of freedom of the reservoir.

Decoherence and its dependence upon the amplitude of the superposition state
is the main obstacle to the implementation of a logic network based on
quantum gates \cite{nielsen,cnot}. The dream of quantum communication and
computation comes up against the nightmare of decoherence mechanisms \cite
{Haroche1}, owing not only to the inevitable action of the surrounding
environment but also to the intrinsic fluctuations in the interaction
parameters required for logic operations \cite{Milburn,Vogel,Serra}. The
need for huge superpositions of qubit states in the practical implementation
of logical operations imposes the requirements that the quantum systems be
totally isolated from the environment and that the interaction parameters
involved be tightly controlled. For this reason, investigation of the
sources of noise in such promising quantum systems is a crucial step towards
the realization of a quantum logic processor. There is also a major effort
being made in present-day research, to discover mechanisms to prevent
decoherence occurring in actual physical systems, by using parity kicks \cite
{Vitali}, stroboscopic feedback \cite{Tombesi}, engineered driving fields 
\cite{pump1} or an engineered reservoir \cite{Zoller,Wineland1,Matos}. In
this light, the main concern of the present work is to analyze the coherence
dynamics and decoherence process in a network composed of two coupled
dissipative oscillators, which may be field modes in dissipative cavities 
\cite{Raimond}, phonon modes of trapped ions \cite{Ions}, phonon modes of
surface electrons in liquid helium \cite{Nelson}, etc. Master equations are
derived for both weakly and strongly interacting oscillators, leading to
interesting results concerning the coherence properties of the joint and the
reduced system states. This work constitutes a first step towards a more
comprehensive treatment of the decoherence process in multipartite quantum
systems.

On attempting to extend the work on decoherence to interacting quantum
oscillators coupled to distinct reservoirs, one faces the problem of
deriving a master equation for different regimes of coupling between
oscillators. Leaving to one side the difficulty of engineering an arbitrary
coupling strength between the oscillators, in the present work we analyze
not only the weak, but also the strong coupling regime, where the coupling
strength between the oscillators is near the typical oscillators
frequencies. In both regimes we assume that the coupling strength between
the oscillators is considerably larger than the system damping rates.

In weak coupling, the coupling strength between the oscillators, labelled $
\ell =1,2$ from here on, is considerably smaller than the typical
frequencies of either oscillator, and the resulting master equation is as if
the two oscillators were decoupled and a decay channel, described by the
Liouville operator ${\cal L}_{\ell }\rho $, can simply be inserted into the
master equation for each oscillator considered. In that case, assuming both
oscillators have the same damping constant, the decoherence time for each
oscillator is unaffected by the interaction with the other one. However,
when the oscillators have different damping constants (field modes in
cavities with different quality factors, for example), we observe that the
``good-quality''\ oscillator gets worse, while the ``bad-quality''\
oscillator gets better, a result which can be employed for quantum state
protection.

In the strong coupling regime, we observe that a cross-decay channel ${\cal L
}_{12}\rho $ appears, besides the usual system-reservoir individual decay
channels ${\cal L}_{\ell }\rho $. This cross-decay channel modifies the
decoherence process of both the joint and the reduced system state, to an
extent depending crucially on the spectral density of the reservoirs. In
fact, in the strong coupling regime, the normal-mode frequencies are
substantially shifted from the typical oscillator frequencies, enabling us
to explore the whole profile of the spectral densities of the reservoirs. We
show how the decoherence process may be controlled by shifting the
normal-mode frequencies to regions of small spectral density of the
reservoirs. Apart from these spectral densities, the competition between the
cross-decay and the usual channels can give rise to a computed delay or
advance of the decoherence process, for eigenstates of the system normal
modes. Note that if a system pointer variable does not commute with the
operator responsible for its coupling to other system, it is clear that the
internal dynamics must interfere in its decoherence time. Thus the
derivation of master equations for strongly interacting systems is a central
task \cite{Walls,Charmicael} in the study of decoherence in quantum networks.

It is worth mentioning some previous work concerned with coupled systems. In
Ref. \cite{Raimond} the authors describe a proposal to achieve reversible
decoherence of a mesoscopic superposition of field states. This proposal is
based on the possibility of performing a reversible coupling between two
Fabry-Perot cavities. In Ref. \cite{Nemes}, a theoretical model of the
experimental proposal in Ref. \cite{Raimond} is given, but in \cite{Nemes}
the inevitable coupling of the resonators to their environment is taken into
account when the reversibility of coherence loss is analyzed. A system of
two coupled cavities is also analyzed in Ref. \cite{Zoubi}, where just one
of the cavities is interacting with a reservoir. In Ref. \cite{Zoubi}, a
master equation is derived in the case of strongly coupled cavities and it
is shown that the relaxation term is not simply the standard one, obtained
by neglecting the interaction between the cavities. It is the aim of the
present paper, in the context of cavity QED, to analyze the reversible
decoherence process of Refs. \cite{Raimond,Nemes}, where two dissipative
cavities are considered, but investigating also the regime of strongly
coupled cavities, as done in \cite{Zoubi}, in which some remarkable
coherence properties appear. A central result extracted from our discussion
is that in a strongly interacting quantum network the decoherence time may
not decrease as the number of systems considered grows. In fact, it may
happen that with many coupled sites \cite{Net}, as with the two coupled
systems analyzed here, the decoherence time increases, depending on the
spectral density of the reservoirs.

Together with the strong coupling between two oscillators we consider a
classical driving field feeding one of the oscillators continuously, which
is intended to drive each of the coupled oscillators to a stationary
coherent state \cite{pump1,MMC}, other than the vacuum state. We select a
specific coupling between the oscillators, which may be responsible for the
dynamics of local transfer of states or state swap between them, as
discussed in Ref. \cite{meu} in relation to cavity QED with weak coupling.

The outline of this paper is as follows. In Sec. II we develop a master
equation for the two coupled lossy oscillators, one of which is under the
action of a classical driving field, and we analyze the weak and strong
field coupling regimes. As the spectral density of the reservoirs plays a
crucial role in the strong coupling regime, in Sec. III we analyze
particular cases of spectral densities. In Sec. IV we develop a c-number
version of the master equation and solve it with reservoir temperatures set
to zero. The solution of the master equation derived for strongly coupled
oscillators is also analyzed in Sec. IV for special initial field states. In
Sec. V we discuss the recurrence and swapping dynamics of the system states.
The central result of the paper, namely the coherence properties of the
system states, which depend on the spectral density of the reservoir and are
strongly affected by the regime of coupling, is presented in Sec. VI. In
Sec. VII we consider the two oscillators to have different damping rates and
demonstrate that the coupling between them, assumed to be larger than these
damping rates, makes the good-quality oscillator worse and the bad-quality
oscillator better. A careful analysis of the entropy excess in our network
is developed in Sec. VIII. Finally, Sec. IX concludes the paper.

\section{The problem of coupled dissipative oscillators: derivation of the
master equation}

General results can be extracted from specific examples of quantum
oscillators, such as field modes in coupled cavities, trapped ions, or
surface electrons in liquid helium. Let us consider a system of two
interacting oscillators under the action of a driving field, as pictured in
Fig. 1. We start from a positive-defined Hamiltonian so that the energy
spectrum has a lower bound which is equal to zero in the absence of the
driving field \cite{FLO,Amir}. The system Hamiltonian is then given by
\begin{equation}gin{eqnarray}
H &=&\hbar \omegamega _{10}\lambdaeft( a_{1}^{\partialtaagger }+\frac{\lambdaambda }{2\omegamega _{10}}
a_{2}^{\partialtaagger }\right) \lambdaeft( a_{1}+\frac{\lambdaambda }{2\omegamega _{10}}
a_{2}\right) +\hbar \omegamega _{20}\lambdaeft( a_{2}^{\partialtaagger }+\frac{\lambdaambda }{
2\omegamega _{20}}a_{1}^{\partialtaagger }\right) \lambdaeft( a_{2}+\frac{\lambdaambda }{2\omegamega
_{20}}a_{1}\right)  \nonumber \\
&&+\hbar \sum_{k}\omegamega _{1k}\lambdaeft( b_{1k}^{\partialtaagger }+\frac{V_{1k}}{\omegamega
_{1k}}a_{1}^{\partialtaagger }\right) \lambdaeft( b_{1k}+\frac{V_{1k}}{\omegamega _{1k}}
a_{1}\right) +\hbar \sum_{k}\omegamega _{2k}\lambdaeft( b_{2k}^{\partialtaagger }+\frac{V_{2k}
}{\omegamega _{2k}}a_{2}^{\partialtaagger }\right) \lambdaeft( b_{2k}+\frac{V_{2k}}{\omegamega
_{2k}}a_{2}\right)  \nonumber \\
&&+\hbar F\lambdaeft( e^{i\omegamega t}+a_{2}^{\partialtaagger }\right) \lambdaeft( e^{-i\omegamega
t}+a_{2}\right) ,  \lambdaabel{Eq1}
\end{eqnarray}
where $a_{\ell }^{\partialtaagger }$ and $a_{\ell }$ are, respectively, the creation
and annihilation operators for the oscillator mode of frequency $\omegamega
_{\ell 0}$, whereas $b_{\ell k}$ and $b_{\ell k}^{\partialtaagger }$ are the
analogous operators for the $k$th bath mode of oscillator $\ell $, whose
corresponding frequency and coupling strength are $\omegamega _{\ell k}$ and $
V_{\ell k}$, respectively. The coupling strength between the oscillators is $
\lambdaambda $ and the classical driving field applied to oscillator
$2$ has intensity $F$ and frequency $\omegamega $. Assuming that the
coupling between the oscillators and their reservoirs satisfies
the condition $\sum_{k}\lambdaeft( V_{\ell k}\right) ^{2}/\omegamega _{\ell
k}\lambdal \omegamega _{\ell 0}$, \ and shifting the origin of the energy
scale to $\hbar F$, we obtain from Eq. (\ref{Eq1}) the Hamiltonian
$H=\sum_{\ell }H_{\ell }+H_{I}$, given by
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
H_{\ell } &=&\hbar \omegamega _{\ell }a_{\ell }^{\partialtaagger }a_{\ell }+\hbar
\sum_{k}\omegamega _{\ell k}b_{\ell k}^{\partialtaagger }b_{\ell k}+\hbar
\sum_{k}V_{\ell k}\lambdaeft( a_{\ell }b_{\ell k}^{\partialtaagger }+a_{\ell }^{\partialtaagger
}b_{\ell k}\right) {\rm {,}}  \lambdaabel{Eq2a} \\
&&+\hbar F\lambdaeft( a_{2}^{\partialtaagger }

\mathop{\rm e}

\nolimits^{-i\omegamega t}+a_{2}

\mathop{\rm e}

{}^{i\omegamega t}\right) \partialtaelta _{\ell 2}{\rm {,}}  \nonumber \\
H_{I} &=&\hbar \lambdaambda \lambdaeft( a_{1}a_{2}^{\partialtaagger }+a_{1}^{\partialtaagger
}a_{2}\right) {\rm {.}}  \lambdaabel{Eq2b}
\end{eqnarray}
Here, $\omegamega _{\ell }$ is related to the natural frequency $\omegamega _{\ell
0} $ by
\end{mathletters}
\begin{equation}gin{equation}
\omegamega _{\ell }=\omegamega _{\ell 0}\lambdaeft( 1+\frac{\lambdaambda ^{2}}{4\omegamega
_{10}\omegamega _{20}}+\partialtaelta _{\ell 2}\frac{F}{\omegamega _{20}}\right) {\rm {,}}
\lambdaabel{EqF}
\end{equation}
such that in the limit of weak coupling between the oscillators ($\lambdaambda
/\omegamega _{\ell 0}\lambdal 1$) and weak amplification process ($F/\omegamega _{20}\lambdal
1 $) we obtain the natural frequencies $\omegamega _{\ell }=\omegamega _{\ell 0}$.
In this limit, we could have started from Hamiltonian $H$ given by Eqs. (\ref
{Eq2a}) and (\ref{Eq2b}) instead of (\ref{Eq1}). Note that since
we are assuming weak couplings between the oscillators and their
reservoirs, it is unnecessary to write these interactions in a
positive-defined form, as done in Eq. (\ref{Eq1}). However, the
positive-defined form for the coupling between the oscillators
ensures an energy spectrum with a lower bound (equal to zero when
the driving field is switched off), whatever the value of the
coupling $\lambdaambda $. Under a unitary transformation represented by
the operator
\begin{equation}gin{equation}
U(t)=\exp \lambdaeft[ -i\omegamega t\sum_{\ell }\lambdaeft( a_{\ell }^{\partialtaagger }a_{\ell
}+\sum_{k}b_{\ell k}^{\partialtaagger }b_{\ell k}\right) \right] ,  \lambdaabel{Eq3}
\end{equation}
we obtain the transformed time-independent Hamiltonian
\begin{equation}gin{equation}
{\cal H}=U^{\partialtaagger }(t)HU(t)+i\frac{dU^{\partialtaagger }(t)}{dt}U(t)=\sum_{\ell }
{\cal H}_{\ell }+{\cal H}_{I}{\rm {,}}  \lambdaabel{Eq4}
\end{equation}
where
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
{\cal H}_{\ell } &=&\hbar \omegamega _{\ell }^{^\primeime }a_{\ell }^{\partialtaagger
}a_{\ell }+\hbar \sum_{k}\omegamega _{\ell k}^{^\primeime }b_{\ell k}^{\partialtaagger
}b_{\ell k}+\hbar \lambdaambda \lambdaeft( a_{1}a_{2}^{\partialtaagger }+a_{1}^{\partialtaagger
}a_{2}\right) {\rm {,}}  \lambdaabel{Eq5a} \\
{\cal H}_{I} &=&\hbar \sum_{k}V_{\ell k}\lambdaeft( a_{\ell }b_{\ell k}^{\partialtaagger
}+a_{\ell }^{\partialtaagger }b_{\ell k}\right) +\hbar F\lambdaeft( a_{2}^{\partialtaagger
}+a_{2}\right) \partialtaelta _{\ell 2}{\rm {,}}  \lambdaabel{Eq5b}
\end{eqnarray}
and the shifted frequencies are given by
\end{mathletters}
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\omegamega _{\ell }^{^\primeime } &=&\omegamega _{\ell }-\omegamega {\rm {,}}  \lambdaabel{Eq6a}
\\
\omegamega _{\ell k}^{^\primeime } &=&\omegamega _{\ell k}-\omegamega {\rm {.}}  \lambdaabel{Eq6b}
\end{eqnarray}

\qquad\ From here on we consider the specific case where $\omegamega
_{1}^{^\primeime }=\omegamega _{2}^{^\primeime }=\Omegamega $ (or $\omegamega _{1}=\omegamega _{2}$
), which links the amplitude of the driving field with the natural
frequencies, as $F=(\omegamega _{10}-\omegamega _{20})\lambdaeft[ 1+\lambdaambda ^{2}/\lambdaeft(
4\omegamega _{10}\omegamega _{20}\right) \right] $. Therefore, in the absence of the
driving field, the condition $\omegamega _{1}^{^\primeime }=\omegamega _{2}^{^\primeime
}=\Omegamega $, implies that $\omegamega _{10}=\omegamega _{20}$. With the condition $
\omegamega _{1}^{^\primeime }=\omegamega _{2}^{^\primeime }=\Omegamega $, the Hamiltonian in Eq.
(\ref{Eq5a}) can be diagonalized through the canonical transformation
\end{mathletters}
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
A_{1} &=&\frac{1}{\sqrt{2}}\lambdaeft( a_{1}+a_{2}\right) ,  \lambdaabel{Eq7a} \\
A_{2} &=&\frac{1}{\sqrt{2}}\lambdaeft( a_{1}-a_{2}\right) {\rm {,}}  \lambdaabel{Eq7b}
\end{eqnarray}
where $A_{1}$ and $A_{2}$ satisfy the same commutation relation as $a_{1}$
and $a_{2}$: $[A_{i},A_{j}]=0$ and$\lambdaeft[ A_{i},A_{j}^{\partialtaagger }\right]
=\partialtaelta _{ij}$. The purpose of these new operators is to decouple the direct
interaction between oscillators $1$ and $2$, described by $\hbar \lambdaambda
\lambdaeft( a_{1}a_{2}^{\partialtaagger }+a_{1}^{\partialtaagger }a_{2}\right) $. Consequently,
indirect interactions between oscillators $1$ and $2$ will be created
through their respective reservoirs, as described by Hamiltonian ${\bf H}=
{\bf H}_{0}+{\bf H}_{I}$, where

\end{mathletters}
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
{\bf H}_{0} &=&\hbar \sum_{\ell }\lambdaeft[ \Omegamega _{\ell }A_{\ell }^{\partialtaagger
}A_{\ell }-(-1)^{\ell }\frac{F}{\sqrt{2}}\lambdaeft( A_{\ell }^{\partialtaagger }+A_{\ell
}\right) +\sum_{k}\omegamega _{\ell k}^{^\primeime }b_{\ell k}^{\partialtaagger }b_{\ell k}
\right] ,  \lambdaabel{Eq8a} \\
{\bf H}_{I} &=&\frac{\hbar }{\sqrt{2}}\sum_{\ell ,k}\lambdaeft[ V_{1k}\lambdaeft(
A_{\ell }b_{1k}^{\partialtaagger }+A_{\ell }^{\partialtaagger }b_{1k}\right) -(-1)^{\ell
}V_{2k}\lambdaeft( A_{\ell }b_{2k}^{\partialtaagger }+A_{\ell }^{\partialtaagger }b_{2k}\right) 
\right] ,  \lambdaabel{Eq8b}
\end{eqnarray}
and the frequencies are given by $\Omegamega _{\ell }=\Omegamega -(-1)^{\ell
}\lambdaambda $. Without direct coupling between oscillators $1$ and $2$, as
modelled by Hamiltonian ${\bf H}$, it becomes simpler to derive the master
equation, following the reasoning in Ref. \cite{Walls}. In the interaction
picture, to the second order of perturbation, the evolution of the density
matrix of the coupled oscillators is given by
\end{mathletters}
\begin{equation}gin{equation}
\frac{d\rho _{12}(t)}{dt}=-\frac{1}{\hbar ^{2}}\int_{0}^{t}dt^{^\primeime }{\rm {
Tr}}_{R}\lambdaeft[ V(t),\lambdaeft[ V(t^{^\primeime }),\rho _{R}(0)\omegatimes \rho _{12}(t)
\right] \right] ,  \lambdaabel{Eq9}
\end{equation}
where $V(t)=\exp \lambdaeft( i{\bf H}_{0}t/\hbar \right) {\bf H}_{I}\exp \lambdaeft( -i
{\bf H}_{0}t/\hbar \right) $. Note that the density matrix in the
interaction picture, $\rho _{12}(t)$, follows from the state vector
transformed by both unitary operators: $U^{\partialtaagger }(t)$ defined by Eq. (\ref
{Eq3}) and $\exp \lambdaeft( i{\bf H}_{0}t/\hbar \right) $. Defining the operator
${\cal O}_{\ell }^{\partialtaagger }(t)=\sum_{k}V_{\ell k}b_{\ell k}^{\partialtaagger }\exp
\lambdaeft( i\omegamega _{\ell k}^{^\primeime }t\right) /\sqrt{2}$, we proceed to obtain
the master equation, assuming that the reservoir frequencies are very
closely spaced, to allow a continuum summation. We have to solve the
integrals appearing in Eq. (\ref{Eq9}), related to correlation functions of
the form
\begin{equation}gin{eqnarray}
\int_{0}^{t}dt^{^\primeime }\lambdaeft\lambdaangle {\cal O}_{\ell }^{\partialtaagger }(t){\cal O}
_{\ell }(t^{^\primeime })\right\rangle

\mathop{\rm e}

\nolimits^{-i\Omegamega _{m}t+i\Omegamega _{n}t^{^\primeime }} &=&\int_{0}^{t}dt^{^\primeime
}\int_{0}^{\infty }\frac{d\omegamega _{\ell k}}{4\pi }V_{\ell }^{2}\lambdaeft( \omegamega
_{\ell k}\right) \sigmama _{\ell }^{2}\lambdaeft( \omegamega _{\ell k}\right) (N_{\ell
}\lambdaeft( \omegamega _{\ell k}\right) +1)  \nonumber \\
&&\thetaimes

\mathop{\rm e}

\nolimits^{i(\omegamega _{\ell k}-\Omegamega _{n}-\omegamega )(t-t^{^\primeime })}

\mathop{\rm e}

\nolimits^{-i(\Omegamega _{m}-\Omegamega _{n})t}{\rm {,}}  \lambdaabel{Eq10}
\end{eqnarray}
where, from here on, $m,n=1,2$, the function $N_{\ell }\lambdaeft( \omegamega _{\ell
k}\right) $ is defined by

\begin{equation}gin{equation}
\lambdaeft\lambdaangle b_{\ell }^{\partialtaagger }(\omegamega _{\ell k})b_{\ell }(\omegamega _{\ell
k^{^\primeime }})\right\rangle =2\pi N_{\ell }\lambdaeft( \omegamega _{\ell k}\right)
\partialtaelta \lambdaeft( \omegamega _{\ell k}-\omegamega _{\ell k^{^\primeime }}\right)
\lambdaabel{Eq11}
\end{equation}
and $\sigmama _{\ell }\lambdaeft( \omegamega _{\ell k}\right) $ is the density of
states of reservoir $\ell $. Performing the variable transformations $\thetaau
=t-t^{^\primeime }$ and $\varepsilon =\omegamega _{\ell k}-\Omegamega _{n}-\omegamega $, we
obtain, for $n=1$ and $2$, respectively,
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\int_{0}^{t}dt^{^\primeime }\lambdaeft\lambdaangle {\cal O}_{\ell }^{\partialtaagger }(t){\cal O}
_{\ell }(t^{^\primeime })\right\rangle

\mathop{\rm e}

\nolimits^{-i\Omegamega _{m}t+i\Omegamega _{1}t^{^\primeime }} &=&\frac{1}{4}

\mathop{\rm e}

{}^{-i(\Omegamega _{m}-\Omegamega _{1})t}\int_{-\omegamega _{\ell }^{+}}^{\infty
}d\varepsilon V_{\ell }^{2}\lambdaeft( \varepsilon +\omegamega _{\ell }^{+}\right)
\sigmama _{\ell }^{2}\lambdaeft( \varepsilon +\omegamega _{\ell }^{+}\right)  \nonumber
\\
&&\thetaimes \lambdaeft[ N_{\ell }\lambdaeft( \varepsilon +\omegamega _{\ell }^{+}\right) +1
\right] \int_{0}^{t}d\thetaau

\mathop{\rm e}

\nolimits^{i\varepsilon \thetaau },  \lambdaabel{Eq12a} \\
\int_{0}^{t}dt^{^\primeime }\lambdaeft\lambdaangle {\cal O}_{\ell }^{\partialtaagger }(t){\cal O}
_{\ell }(t^{^\primeime })\right\rangle

\mathop{\rm e}

\nolimits^{-i\Omegamega _{m}t+i\Omegamega _{2}t^{^\primeime }} &=&\frac{1}{4}

\mathop{\rm e}

{}^{-i(\Omegamega _{m}-\Omegamega _{2})t}\int_{-\omegamega _{\ell }^{-}}^{\infty
}d\varepsilon V_{\ell }^{2}\lambdaeft( \varepsilon +\omegamega _{\ell }^{-}\right)
\sigmama _{\ell }^{2}\lambdaeft( \varepsilon +\omegamega _{\ell }^{-}\right)  \nonumber
\\
&&\thetaimes \lambdaeft[ N_{\ell }\lambdaeft( \varepsilon +\omegamega _{\ell }^{-}\right) +1
\right] \int_{0}^{t}d\thetaau

\mathop{\rm e}

\nolimits^{i\varepsilon \thetaau },  \lambdaabel{Eq12b}
\end{eqnarray}
where the frequency $\omegamega _{\ell }$ has been split into two effective
frequencies corresponding to the normal modes of the coupled oscillators
(note that $\omegamega _{1}=\omegamega _{2}$):
\end{mathletters}
\begin{equation}gin{equation}
\omegamega _{\ell }^{\pm }=\omegamega _{\ell }\pm \lambdaambda .  \lambdaabel{Eq13}
\end{equation}
We have suppressed the contribution of the Cauchy principal value since it
represents only a small shift in the frequency $\omegamega _{\ell }^{\pm }$. We
note that the minimum value of $\omegamega _{\ell }^{-}$ is $F\lambdaeft[ \omegamega
_{10}/\lambdaeft( \omegamega _{10}+\omegamega _{20}\right) \right] $ which follows from $
\lambdaambda =2\omegamega _{20}$. Next, we discuss both regimes: $(i)$ the weak
coupling regime $\lambdaambda /\omegamega _{\ell 0}\lambdal 1$, and $(ii)$ the strong
coupling regime $\lambdaambda /\omegamega _{\ell 0}\approx 1$. For the strong
coupling regime we will assume $\lambdaambda /\omegamega _{\ell 0}=2$ in order to
minimize $\omegamega _{\ell }^{-}$, which becomes zero when the driving field is
switched off.

As usual, we consider that $V_{\ell }$, $\sigmama _{\ell }$,\ and $N_{\ell }$\
are functions that vary slowly around the frequency $\omegamega _{\ell }^{\pm }$
, an assumption which does not apply to the function $N_{\ell }\lambdaeft( \omegamega
_{\ell }^{-}\right) =\lambdaeft[ \exp \lambdaeft( \hbar \omegamega _{\ell }^{-}/kT\right)
-1\right] ^{-1}$ (taking the reservoir to be in thermal equilibrium at
temperature $T$) in the strong coupling regime when $F=0$, since in this
case $\omegamega _{\ell }^{-}\approx 0$. However, this regime (even with $F=0$)
can safely be applied to a reservoir at absolute zero, the situation we
analyze in the present work. We observe that, in practice, $N_{\ell }\lambdaeft(
\omegamega _{\ell }^{-}\right) \approx 0$ whenever the shift in the frequency $
\omegamega _{\ell }^{-}$, arising from the contribution of the Cauchy principal
value, becomes sufficiently greater than $kT/\hbar $. Note that the last
integrals in Eqs. (\ref{Eq12a}) and (\ref{Eq12b}) contribute significantly
only when $|\varepsilon t|\lambdaessapprox 1$, so that if we extend the upper
limit of the time integration to infinity, the expressions for the
correlation functions become
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\int_{0}^{t}dt^{^\primeime }\lambdaeft\lambdaangle {\cal O}_{\ell }^{\partialtaagger }(t){\cal O}
_{\ell }(t^{^\primeime })\right\rangle

\mathop{\rm e}

\nolimits^{-i\Omegamega _{m}t+i\Omegamega _{1}t^{^\primeime }} &=&\frac{\gammaamma _{\ell
}(\omegamega _{\ell }^{+})}{2}\lambdaeft[ N_{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right)
+1\right]

\mathop{\rm e}

{}^{-i(\Omegamega _{m}-\Omegamega _{1})t},  \lambdaabel{Eq14a} \\
\int_{0}^{t}dt^{^\primeime }\lambdaeft\lambdaangle {\cal O}_{\ell }^{\partialtaagger }(t){\cal O}
_{\ell }(t^{^\primeime })\right\rangle

\mathop{\rm e}

\nolimits^{-i\Omegamega _{m}t+i\Omegamega _{2}t^{^\primeime }} &=&\frac{\gammaamma _{\ell
}(\omegamega _{\ell }^{-})}{2}\lambdaeft[ N_{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right)
+1\right]

\mathop{\rm e}

{}^{-i(\Omegamega _{m}-\Omegamega _{2})t},  \lambdaabel{Eq14b}
\end{eqnarray}
where the damping rates are defined as
\end{mathletters}
\begin{equation}gin{equation}
\gammaamma _{\ell }(\omegamega _{\ell }^{\pm })=\frac{1}{2}V_{\ell }^{2}\lambdaeft(
\omegamega _{\ell }^{\pm }\right) \sigmama _{\ell }^{2}\lambdaeft( \omegamega _{\ell }^{\pm
}\right) \int_{-\omegamega _{\ell }^{\pm }}^{\infty }d\varepsilon \partialtaelta
(\varepsilon ).  \lambdaabel{Eq15}
\end{equation}

Defining $\Gammaamma _{\ell }=$ $V_{\ell }^{2}\lambdaeft( \omegamega _{\ell 0}\right)
\sigmama _{\ell }^{2}\lambdaeft( \omegamega _{\ell 0}\right) $, in the weak coupling
regime, where $\omegamega _{\ell }^{\pm }\approx \omegamega _{\ell }\approx \omegamega
_{\ell 0}$, we obtain from Eq. (\ref{Eq15}) the result
\begin{equation}gin{equation}
\gammaamma _{\ell }(\omegamega _{\ell 0})\approx \frac{1}{2}V_{\ell }^{2}\lambdaeft(
\omegamega _{\ell 0}\right) \sigmama _{\ell }^{2}\lambdaeft( \omegamega _{\ell 0}\right) =
\frac{\Gammaamma _{\ell }}{2},  \lambdaabel{Eq16}
\end{equation}
while in the strong coupling regime, where $\omegamega _{\ell }^{+}\gammag $ $\omegamega
_{\ell }^{-}$, we have
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\gammaamma _{\ell }(\omegamega _{\ell }^{+}) &\approx &\frac{1}{2}V_{\ell
}^{2}\lambdaeft( \omegamega _{\ell }^{+}\right) \sigmama _{\ell }^{2}\lambdaeft( \omegamega
_{\ell }^{+}\right) ,  \lambdaabel{Eq17a} \\
\gammaamma _{\ell }(\omegamega _{\ell }^{-}) &\approx &\frac{1}{4}V_{\ell
}^{2}\lambdaeft( \omegamega _{\ell }^{-}\right) \sigmama _{\ell }^{2}\lambdaeft( \omegamega
_{\ell }^{-}\right) {\rm {\quad }}{\rm {when\;}}F=0  \lambdaabel{Eq17b} \\
\gammaamma _{\ell }(\omegamega _{\ell }^{-}) &\approx &\frac{1}{2}V_{\ell
}^{2}\lambdaeft( \omegamega _{\ell }^{-}\right) \sigmama _{\ell }^{2}\lambdaeft( \omegamega
_{\ell }^{-}\right) {\rm {\quad }}{\rm {when\;}}F\neq 0{\rm {.}}
\lambdaabel{Eq17c}
\end{eqnarray}

>From the above results for the correlation functions, we observe that the
master equation for the strong coupling regime includes that for the weak
coupling regime. In fact, with $\gammaamma _{\ell }(\omegamega _{\ell }^{\pm
})\approx \Gammaamma _{\ell }/2$ and $N_{\ell }(\omegamega _{\ell }^{\pm })\approx
N_{\ell }(\omegamega _{\ell 0})$ we get the master equation for the weak
coupling regime from that for the strong coupling regime which, described in
the Schr\"{o}dinger picture via the mode operators $a_{\ell }$ and $a_{\ell
}^{\partialtaagger }$, reads
\end{mathletters}
\begin{equation}gin{eqnarray}
\frac{d\rho _{12}}{dt} &=&\sum_{\ell }\lambdaeft\{ i\lambdaeft[ \rho _{12},\Omegamega
a_{\ell }^{\partialtaagger }a_{\ell }+\lambdaambda \sum_{m\neq \ell }a_{\ell }^{\partialtaagger
}a_{m}+F\lambdaeft( a_{2}^{\partialtaagger }+a_{2}\right) \partialtaelta _{\ell 2}\right] \right.
\nonumber \\
&&+\frac{F}{2\lambdaeft( \Omegamega ^{2}-\lambdaambda ^{2}\right) }\lambdaeft\{ \lambdaeft( \Omegamega
\partialtaelta _{\ell 1}-\lambdaambda \partialtaelta _{\ell 2}\right) \lambdaeft[ \gammaamma _{\ell
}\lambdaeft( \omegamega _{\ell }^{-}\right) -\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell
}^{+}\right) \right] +\lambdaeft( \Omegamega \partialtaelta _{\ell 2}-\lambdaambda \partialtaelta _{\ell
1}\right) \right.  \nonumber \\
&&\thetaimes \lambdaeft. \lambdaeft[ 2\gammaamma _{\ell }\lambdaeft( \omegamega \right) -\gammaamma _{\ell
}\lambdaeft( \omegamega _{\ell }^{-}\right) -\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell
}^{+}\right) \right] \right\} \lambdaeft[ \rho _{12},a_{\ell }-a_{\ell }^{\partialtaagger
}\right]  \nonumber \\
&&+\frac{1}{2}\lambdaeft[ \gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right)
N_{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right) +\gammaamma _{\ell }\lambdaeft( \omegamega
_{\ell }^{-}\right) N_{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right) \right]
\lambdaeft( \lambdaeft[ \lambdaeft[ a_{\ell }^{\partialtaagger },\rho _{12}\right] ,a_{\ell }\right]
+\lambdaeft[ a_{\ell }^{\partialtaagger },\lambdaeft[ \rho _{12},a_{\ell }\right] \right]
\right)  \nonumber \\
&&+\frac{1}{2}\lambdaeft[ \gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right)
+\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right) \right] \lambdaeft( \lambdaeft[
a_{\ell }\rho _{12},a_{\ell }^{\partialtaagger }\right] +\lambdaeft[ a_{\ell },\rho
_{12}a_{\ell }^{\partialtaagger }\right] \right)  \nonumber \\
&&+\frac{1}{2}\sum_{m\neq \ell }\lambdaeft[ \gammaamma _{\ell }\lambdaeft( \omegamega _{\ell
}^{+}\right) -\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right) \right]
\lambdaeft( \lambdaeft[ a_{m}\rho _{12},a_{\ell }^{\partialtaagger }\right] +\lambdaeft[ a_{\ell
},\rho _{12}a_{m}^{\partialtaagger }\right] \right)  \nonumber \\
&&\lambdaeft. +\frac{1}{2}\lambdaeft[ \gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right)
N_{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right) -\gammaamma _{\ell }\lambdaeft( \omegamega
_{\ell }^{-}\right) N_{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right) \right]
\lambdaeft( \lambdaeft[ \lambdaeft[ a_{2}^{\partialtaagger },\rho _{12}\right] ,a_{1}\right] +\lambdaeft[
a_{1}^{\partialtaagger },\lambdaeft[ \rho _{12},a_{2}\right] \right] \right) \right\} .
\lambdaabel{Eq18}
\end{eqnarray}
We stress that when the driving field is switched off, the second term in
Eq. (\ref{Eq18}) under the summation on $\ell $ disappears. Otherwise,
noting that $\Omegamega ^{2}-\lambdaambda ^{2}=(\omegamega _{\ell }^{-}-\omegamega )(\omegamega
_{\ell }^{+}-\omegamega )$, when $F\neq 0$ and $\omegamega =\omegamega _{\ell }^{\pm }$,
this term becomes
\begin{equation}gin{equation}
\frac{F}{4\lambdaambda }\lambdaeft[ \gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right)
-\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right) \right] \lambdaeft( \partialtaelta
_{\ell 2}\mp \partialtaelta _{\ell 1}\right) \lambdaeft[ \rho _{12},a_{\ell }-a_{\ell
}^{\partialtaagger }\right] {\rm {.}}  \lambdaabel{Eq18l}
\end{equation}

\qquad \qquad Next, we derive the master equations in the weak and strong
coupling regime. For simplicity, we define for the strong coupling regime $
\gammaamma _{\ell }(\omegamega _{\ell }^{\pm })\equiv $ $\gammaamma _{\ell }^{\pm }$.
Besides, we assume from here on the resonance condition for the driving
field $\omegamega =\omegamega _{20}$.

\subsection{Weak coupling regime}

In the weak coupling regime, where $\omegamega _{\ell }^{\pm }\approx \omegamega
_{\ell }\approx \omegamega _{\ell 0}$, the resonant condition for the driving
field implies $\gammaamma _{\ell }\lambdaeft( \omegamega \right) \approx \gammaamma _{\ell
}^{\pm }\approx \Gammaamma _{\ell }/2$, so that the master equation becomes
\begin{equation}gin{equation}
\frac{d\rho _{12}}{dt}=\frac{i}{\hbar }\lambdaeft[ \rho _{12},{\sf H}_{0}\right]
+\sum_{\ell }{\cal L}_{\ell }\rho _{12}{\rm {,}}  \lambdaabel{Eq19}
\end{equation}
where
\begin{equation}gin{equation}
{\sf H}_{0}=\hbar \sum_{\ell }\lambdaeft[ \Omegamega a_{\ell }^{\partialtaagger }a_{\ell
}+\lambdaambda \sum_{m\neq \ell }a_{\ell }^{\partialtaagger }a_{m}+F\lambdaeft( a_{2}^{\partialtaagger
}+a_{2}\right) \partialtaelta _{\ell 2}\right] {\rm {,}}  \lambdaabel{Eq20}
\end{equation}
and the Liouville operator ${\cal L}_{\ell }\rho _{12}$ is given by the
usual operator structure
\begin{equation}gin{eqnarray}
{\cal L}_{\ell }\rho _{12} &\equiv &\frac{1}{2}\Gammaamma _{\ell }\lambdaeft\{
N_{\ell }\lambdaeft( \omegamega _{\ell 0}\right) \lambdaeft( \lambdaeft[ a_{\ell }^{\partialtaagger
}\rho _{12},a_{\ell }\right] +\lambdaeft[ a_{\ell }^{\partialtaagger },\rho _{12}a_{\ell }
\right] \right) \right.  \nonumber \\
&&+\lambdaeft. \lambdaeft( N_{\ell }\lambdaeft( \omegamega _{\ell 0}\right) +1\right) \lambdaeft(
\lambdaeft[ a_{\ell }\rho _{12},a_{\ell }^{\partialtaagger }\right] +\lambdaeft[ a_{\ell
},\rho _{12}a_{\ell }^{\partialtaagger }\right] \right) \right\} {\rm {.}}
\lambdaabel{Eq21}
\end{eqnarray}
Therefore, in the weak coupling regime and assuming that the driving field
is resonant with oscillator $2$, the interaction between the field modes and
the amplification process appears only in the von Neumann term of the master
equation, and does not affect the dissipative mechanism of the individual
cavities. However, when the driving field is out of resonance with
oscillator $2$, a correction term is added to the Liouville operator ${\cal L
}_{\ell }\rho _{12}$, given by
\begin{equation}gin{equation}
\frac{F}{\lambdaeft( \Omegamega ^{2}-\lambdaambda ^{2}\right) }\lambdaeft( \Omegamega \partialtaelta _{\ell
2}-\lambdaambda \partialtaelta _{\ell 1}\right) \lambdaeft[ \gammaamma _{\ell }\lambdaeft( \omegamega
\right) -\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell 0}\right) \right] \lambdaeft[ \rho
_{12},a_{\ell }-a_{\ell }^{\partialtaagger }\right] .  \lambdaabel{Eq21l}
\end{equation}
For a strong amplification process ($F/\omegamega _{20}\approx 1$), this
correction becomes $F\Omegamega /\lambdaeft( \Omegamega ^{2}-\lambdaambda ^{2}\right) \approx
1/2$ for mode $\ell =2$, being ignored for mode $\ell =1$, since $F\lambdaambda
/\lambdaeft( \Omegamega ^{2}-\lambdaambda ^{2}\right) \lambdal 1$.

\subsection{Strong coupling regime}

Below we present the master equation for the strong coupling regime
considering the general situation where both driving field and reservoir
temperature are present. However, as discussed above, when switching off the
driving field we must consider reservoirs at absolute zero for our results
to be valid. In this regime, where $\lambdaambda /\omegamega _{\ell }\approx 1$, the
master equation is written
\begin{equation}gin{equation}
\frac{d\rho _{12}}{dt}=\frac{i}{\hbar }\lambdaeft[ \rho _{12},{\sf H}_{0}\right]
+\sum_{\ell }{\cal L}_{\ell }\rho _{12}+{\cal L}_{12}\rho _{12},
\lambdaabel{Eq22}
\end{equation}
where
\begin{equation}gin{eqnarray}
{\cal L}_{\ell }\rho _{12} &=&\frac{1}{2}\lambdaeft( \gammaamma _{\ell }^{+}N_{\ell
}\lambdaeft( \omegamega _{\ell }^{+}\right) +\gammaamma _{\ell }^{-}N_{\ell }\lambdaeft(
\omegamega _{\ell }^{-}\right) \right) \lambdaeft( \lambdaeft[ \lambdaeft[ a_{\ell }^{\partialtaagger
},\rho _{12}\right] ,a_{\ell }\right] +\lambdaeft[ a_{\ell }^{\partialtaagger },\lambdaeft[
\rho _{12},a_{\ell }\right] \right] \right)  \nonumber \\
&&+\frac{1}{2}\lambdaeft( \gammaamma _{\ell }^{+}+\gammaamma _{\ell }^{-}\right) \lambdaeft( 
\lambdaeft[ a_{\ell }\rho _{12},a_{\ell }^{\partialtaagger }\right] +\lambdaeft[ a_{\ell
},\rho _{12}a_{\ell }^{\partialtaagger }\right] \right)  \nonumber \\
&&+\frac{F}{2\lambdaeft( \Omegamega ^{2}-\lambdaambda ^{2}\right) }\lambdaeft\{ \lambdaeft( \Omegamega
\partialtaelta _{\ell 1}-\lambdaambda \partialtaelta _{\ell 2}\right) \lambdaeft( \gammaamma _{\ell
}^{-}-\gammaamma _{\ell }^{+}\right) +\lambdaeft( \Omegamega \partialtaelta _{\ell 2}-\lambdaambda
\partialtaelta _{\ell 1}\right) \right.  \nonumber \\
&&\thetaimes \lambdaeft. \lambdaeft[ 2\gammaamma _{\ell }\lambdaeft( \omegamega _{20}\right) -\gammaamma
_{\ell }^{-}-\gammaamma _{\ell }^{+}\right] \right\} \lambdaeft[ \rho _{12},a_{\ell
}-a_{\ell }^{\partialtaagger }\right] ,  \lambdaabel{Eq23}
\end{eqnarray}
and a cross-decay channel is included via the Liouville operator
\begin{equation}gin{eqnarray}
{\cal L}_{12}\rho _{12} &=&\frac{1}{2}\sum_{\ell }\lambdaeft\{ \sum_{m\neq \ell
}\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) \lambdaeft( \lambdaeft[
a_{m}\rho _{12},a_{\ell }^{\partialtaagger }\right] +\lambdaeft[ a_{\ell },\rho
_{12}a_{m}^{\partialtaagger }\right] \right) \right.  \nonumber \\
&&\lambdaeft. +\lambdaeft[ \gammaamma _{\ell }^{+}N_{\ell }\lambdaeft( \omegamega _{\ell
}^{+}\right) -\gammaamma _{\ell }^{-}N_{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right) 
\right] \lambdaeft( \lambdaeft[ \lambdaeft[ a_{2}^{\partialtaagger },\rho _{12}\right] ,a_{1}\right]
+\lambdaeft[ a_{1}^{\partialtaagger },\lambdaeft[ \rho _{12},a_{2}\right] \right] \right)
\right\} .  \lambdaabel{Eq24}
\end{eqnarray}
Note that in the weak coupling regime, where $\gammaamma _{\ell }^{\pm }\approx
\Gammaamma _{\ell }/2$, we obtain Eq.(\ref{Eq19}) from Eq. (\ref{Eq22}).
Evidently, the Liouville operator accounting for the cross-decay channel, $
{\cal L}_{12}\rho _{12}$, owing to the strong coupling between the
oscillators, can be of the same order of magnitude as the direct-decay
channels ${\cal L}_{1}\rho _{12}$ and ${\cal L}_{2}\rho _{12}$. In this
regime, in contrast to the weak coupling case, a strong driving field ($
F/\omegamega _{20}\approx 1$) modifies the Liouville operator ${\cal L}_{\ell
}\rho _{12}$, such that $F\Omegamega /\lambdaeft( \Omegamega ^{2}-\lambdaambda ^{2}\right)
,F\lambdaambda /\lambdaeft( \Omegamega ^{2}-\lambdaambda ^{2}\right) \approx 1$, independently
of the resonance consition. We observe that for both reservoirs at absolute
zero, the cross-decay channel is lost for the case where $\gammaamma _{\ell
}^{+}=\gammaamma _{\ell }^{-}$, which may occur, as discussed below, depending
on the spectral density of the reservoirs. In what follows it will become
clear that the cross-decay channel represented by the Liouville operator (
\ref{Eq24}) leads to interesting results concerning the decoherence process
in strongly coupled oscillators.

\subsection{The split of the damping rate}

It is interesting to note that for the coupled dissipative oscillators the
damping rate $\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell 0}\right) =\Gammaamma _{\ell
}/2$ for mode $\ell $ splits into $\gammaamma _{\ell }^{+}$ and $\gammaamma _{\ell
}^{-}$. To illustrate this mechanism, we assume a Lorentzian coupling $
V_{\ell }$ between the oscillators and their respective reservoirs, such
that the damping function $\Gammaamma _{\ell }\lambdaeft( \chi \right) $, centered on
frequency $\chi _{0}$, is given by

\begin{equation}gin{equation}
\Gammaamma _{\ell }\lambdaeft( \chi \right) =\sigmama _{\ell }^{2}\frac{\Upsilon }{
\lambdaeft( \chi -\chi _{0}\right) ^{2}+\Upsilon ^{2}}{\rm {,}}  \lambdaabel{Eq25}
\end{equation}
with the parameter $\Upsilon $ accounting for the spectral sharpness around
the mode frequency.{\lambdaarge \ }From the above expression and remembering,
from Eq. (\ref{Eq16}), that the frequency $\omegamega _{\ell }$ splits into two
shifted frequencies $\omegamega _{\ell }^{\pm }$, we obtain the double
Lorentzian function depicted in Fig. 2
\begin{equation}gin{eqnarray}
\Gammaamma _{\ell }\lambdaeft( \chi \right) &=&\frac{\Upsilon \sigmama _{\ell }^{2}}{2}
\lambdaeft( \frac{1}{\lambdaeft( \chi -\omegamega _{\ell }^{+}\right) ^{2}+\Upsilon ^{2}}+
\frac{1}{\lambdaeft( \chi -\omegamega _{\ell }^{-}\right) ^{2}+\Upsilon ^{2}}\right)
\nonumber \\
&=&\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right) +\gammaamma _{\ell }\lambdaeft(
\omegamega _{\ell }^{-}\right) {\rm {,}}  \lambdaabel{Eq26}
\end{eqnarray}
with maxima on $\omegamega _{\ell }^{\pm }$. In fact, from master equation (\ref
{Eq18}), we observe that in the weak coupling regime, when $\gammaamma _{\ell
}\lambdaeft( \omegamega _{\ell }^{+}\right) +\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell
}^{-}\right) =\Gammaamma _{\ell }(\omegamega _{\ell 0})=$ $\Gammaamma _{\ell }$ (since $
\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right) =\gammaamma _{\ell }\lambdaeft(
\omegamega _{\ell }^{-}\right) =\Gammaamma _{\ell }/2$) , we obtain the expected
Liouville form for two independent dissipative oscillators. From Eq. (\ref
{Eq26}) it is immediately obvious that in the weak coupling regime, where $
\omegamega _{\ell }^{\pm }\approx \omegamega _{\ell 0}$, the damping function
presents only one peak, shown by the dotted line in Fig. 2. In this regime,
the damping rate, assumed to be the maximum of a sharp-peaked damping
function, i.e., $\Gammaamma _{\ell }$ (for a small value of $\Upsilon $),
becomes twice the value designated for $\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell
}^{\pm }\right) $. As $\lambdaambda =\pm (\omegamega _{\ell }^{\pm }-\omegamega _{\ell })$
increases, the damping function splits into two Lorentzian functions whose
peak heights are half the original value $\Gammaamma _{\ell }$, as dictated by
the master equation (\ref{Eq22}) and shown by the solid line in Fig. 2. The
dashed line shows the situation where the two peaks can be clearly
distinguished, on the way to the strong coupling regime, $\lambdaambda /\omegamega
_{\ell 0}\approx 1$, where the peak centered on $\omegamega _{\ell }^{-}$ shifts
to around the value $F\lambdaeft[ \omegamega _{10}/\lambdaeft( \omegamega _{10}+\omegamega
_{20}\right) \right] $, which can be made as smaller as we wish by
decreasing the amplitude of the driving field. In practice, the effect of
the strong coupling between the oscillators is essentially to shift the
normal-mode frequency $\omegamega _{\ell }^{\pm }$ to regions far way from the
natural frequency of the oscillator $\omegamega _{\ell 0}$, where the spectral
density of the reservoir may be significantly different from that around $
\omegamega _{\ell 0}$. In this connection, the spectral density of the reservoir
$\sigmama _{\ell }\lambdaeft( \omegamega _{\ell k}\right) $ plays a crucial role in the
dissipative dynamics of strongly coupled oscillators, since the magnitude of
\ the damping rate $\gammaamma _{\ell }\lambdaeft( \omegamega _{\ell }\right) $ depends
on $\sigmama _{\ell }\lambdaeft( \omegamega _{\ell }\right) $. For this reason, we next
analyze reservoirs with different spectral densities in order to illustrate
the interesting features arising from the strong coupling regime. Evidently,
the physical systems under consideration and their dissipative mechanisms
(i.e., the nature of the reservoirs and their spectral densities), will be
decisive for our analysis.

\section{Spectral densities of the reservoirs}

It is possible that specific spectral densities could be achieved through
engineered reservoirs, a program which has recently attracted considerable
attention also in attempts to control the decoherence process of quantum
states \cite{Zoller,Wineland1,Matos}. Therefore, the results we present
below, depending crucially on the spectral density of the reservoir, might
provide a motivation for future theoretical proposals on engineered
reservoirs. For the following discussion we consider the strong coupling
regime and set both driving field and reservoir temperatures to zero ($F,T=0$
) and remember, from the definition $\gammaamma _{\ell }(\omegamega _{\ell }^{\pm
})\equiv $ $\gammaamma _{\ell }^{\pm }$, that the parameter $\gammaamma _{\ell
}^{\pm }$\ depends on the reservoir spectral density around $\omegamega _{\ell
}^{\pm }$.

\subsection{Markovian white noise}

We start with the simplest case of Markovian white noise, where the spectral
density of the reservoir is invariant over translation in frequency space,
as depicted in Fig 3(a). In this case, assuming a Lorentzian coupling
between the oscillators and their respective reservoirs, centered around the
effective frequencies $\omegamega _{\ell }^{\pm }$, as in Eq. (\ref{Eq26}), we
get $\gammaamma _{\ell }^{-}=\gammaamma _{\ell }^{+}/2=\Gammaamma _{\ell }/4$. In fact,
the effective frequency $\omegamega _{\ell }^{-}$ shifts to around zero (when $
F=0)$, and so, the system-reservoir coupling $\gammaamma _{\ell }^{-}$ becomes
half the value $\gammaamma _{\ell }^{+}=\Gammaamma _{\ell }/2$, as can be concluded
immediately from Eqs. (\ref{Eq17a}), (\ref{Eq17b}) and \ref{Eq17c}. The
system-reservoir couplings around $\omegamega _{\ell }^{\pm }$ are represented
by shaded regions in Fig. 3(a). Therefore, for a Markovian white noise
reservoir, strong coupling between the oscillators delays the decoherence
time of a joint state which is an eigenstate of normal mode $\omegamega _{\ell
}^{-}$, as will be discussed latter. Next, we analyze two cases of
non-Markovian colored noise, still assuming that $V_{\ell }$, $\sigmama _{\ell
}$,\ and $N_{\ell }$\ are functions that vary slowly around the frequency $
\omegamega _{\ell }^{\pm }$, as discussed above.

\subsection{A Lorentzian spectral density}

Let us consider a Lorentzian spectral density of the reservoir which goes to
zero at both effective frequencies $\omegamega _{\ell }^{-}$ and $\omegamega _{\ell
}^{+}$ (Fig. 3(b)). The achievement of a reservoir spectral density for
which $\sigmama _{\ell }\lambdaeft( \omegamega _{\ell }^{\pm }\right) \approx 0$ would
result in a damping function $\gammaamma _{\ell }(\omegamega _{\ell }^{\pm })$
arising from both terms on the right hand side of Eq. (\ref{Eq26}) (i.e., $
\Upsilon \sigmama _{\ell }^{2}/\lambdaeft[ \lambdaeft( \chi -\omegamega _{\ell }^{\pm
}\right) ^{2}+\Upsilon ^{2}\right] \approx 0$), neither of which would
contribute significantly to the relaxation process. Therefore, for such a
Lorentzian coupling between the field modes and their respective reservoirs,
we get from Eqs. (\ref{Eq17a}), (\ref{Eq17b}) and \ref{Eq17c} the rates $
\gammaamma _{\ell }^{\pm }\lambdal \Gammaamma _{\ell }$, since both effective frequencies
$\omegamega _{\ell }^{+}$ and $\omegamega _{\ell }^{-}$ shift to regions where the
spectral densities of the reservoir are considerable smaller than that
around $\omegamega _{\ell 0}$. In view of the difficulty of obtaining a
Lorentzian spectral density sufficiently small around $\omegamega _{\ell }^{+}$,
such that $\sigmama _{\ell }\lambdaeft( \omegamega _{\ell }^{+}\right) \approx 0$, we
next turn to a more realistic case, which we call a wide Lorentzian spectral
density.

\subsection{A wide Lorentzian spectral density}

Finally, we consider a wide Lorentzian spectral density, which goes to zero
at the effective frequency $\omegamega _{\ell }^{-}$, reaches a flat maximum
around $\omegamega _{\ell 0}$, continues on it up to around $\omegamega _{\ell }^{+}$
and only reaches zero again at frequencies far beyond $\omegamega _{\ell }^{+}$,
as depicted in Fig. 3(c). Differently from the Lorentzian spectral density,
in this case only the damping function $\gammaamma _{\ell }^{-}$ arising from
the second term on the right hand side of Eq. (\ref{Eq26}) (i.e., $\Upsilon
\sigmama _{\ell }^{2}/\lambdaeft[ \lambdaeft( \chi -\omegamega _{\ell }^{-}\right)
^{2}+\Upsilon ^{2}\right] \approx 0$) does not contribute significantly to
the relaxation process; i.e., $\gammaamma _{\ell }^{-}\lambdal \Gammaamma _{\ell }$,
since the resonances around the effective frequency $\omegamega _{\ell }^{-}$
can be disregarded. Only the system-reservoir coupling around the effective
frequency $\omegamega _{\ell }^{+}$ contributes to the relaxation process of the
cavity mode ($\gammaamma _{\ell }^{+}=\Gammaamma _{\ell }/2$). Such a model can be
considered as similar to subhomic dissipation \cite{CL,Leggett} and can be
reasonably applied for a variety of phonon-like spectral densities where the
de Debye model applies very well in the domain of small frequencies.

\section{The Fokker-Planck equation}

Using the standard procedures, we derive a c-number version of the master
equation (\ref{Eq18}) for the Glauber P-function \cite{MLW,WM} given by
\begin{equation}gin{equation}
\frac{dP}{dt}=\frac{1}{2}\sum_{\ell ,m}\lambdaeft[ \Pi _{\ell }+C_{\ell }(\alpha
_{1},\alpha _{2})\frac{\partial }{\partial \alpha _{\ell }}+D_{\ell m}\frac{
\partial ^{2}}{\partial \alpha _{\ell }\partial \alpha _{m}^{\ast }}+H.c.
\right] P{\rm {,}}  \lambdaabel{Eq27}
\end{equation}
where the function $C_{\ell }(\alpha _{1},\alpha _{2})$ and the matrix
elements $D_{\ell m}$ satisfy
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
C_{\ell }(\alpha _{1},\alpha _{2}) &=&B_{\ell }+\alpha _{\ell }E_{\ell
}^{+}+\alpha _{\ell -(-1)^{\ell }}E_{\ell }^{-},  \lambdaabel{Eq28a} \\
D_{\ell \ell } &=&\gammaamma _{\ell }^{+}N_{\ell }\lambdaeft( \omegamega _{\ell
}^{+}\right) +\gammaamma _{\ell }^{-}N_{\ell }\lambdaeft( \omegamega _{\ell }^{-}\right) ,
\lambdaabel{Eq28b} \\
D_{12} &=&D_{21}=\frac{1}{2}\sum_{\ell }\lambdaeft[ \gammaamma _{\ell }^{+}N_{\ell
}\lambdaeft( \omegamega _{\ell }^{+}\right) -\gammaamma _{\ell }^{-}N_{\ell }\lambdaeft(
\omegamega _{\ell }^{-}\right) \right] ,  \lambdaabel{Eq28c}
\end{eqnarray}
while the parameters $\Pi _{\ell }$, $E_{\ell }^{\pm }$\ and $B_{\ell }$\
are defined by{\LambdaARGE \ }
\end{mathletters}
\begin{equation}gin{eqnarray}
\Pi _{\ell } &=&\frac{1}{2}\lambdaeft( \gammaamma _{\ell }^{-}+\gammaamma _{\ell
}^{+}\right) ,  \lambdaabel{Eq29a} \\
E_{\ell }^{+} &=&\frac{1}{2}\lambdaeft( \gammaamma _{\ell }^{+}+\gammaamma _{\ell
}^{-}+2i\Omegamega \right) ,  \lambdaabel{Eq29b} \\
E_{\ell }^{-} &=&\frac{1}{2}\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell
}^{-}+2i\lambdaambda \right) ,  \lambdaabel{Eq29c} \\
B_{\ell } &=&iF\partialtaelta _{\ell 2}-\frac{F}{2\lambdaeft( \Omegamega ^{2}-\lambdaambda
^{2}\right) }\lambdaeft[ \lambdaeft( \Omegamega \partialtaelta _{\ell 2}-\lambdaambda \partialtaelta _{\ell
1}\right) \lambdaeft( 2\gammaamma _{\ell }(\omegamega _{20})-\gammaamma _{\ell }^{-}-\gammaamma
_{\ell }^{+}\right) \right.  \nonumber \\
&&+\lambdaeft. \lambdaeft( \Omegamega \partialtaelta _{\ell 1}-\lambdaambda \partialtaelta _{\ell 2}\right)
\lambdaeft( \gammaamma _{\ell }^{-}-\gammaamma _{\ell }^{+}\right) \right] {.}
\lambdaabel{Eq29d}
\end{eqnarray}
Note that the drift coefficient $C_{\ell }$ is linearly dependent upon both
variables $\alpha _{1}$ and $\alpha _{2}$, while the diffusion coefficient $
D_{mn}$ is constant, determining an Orstein-Uhlenbeck process \cite{Gardiner}
.

In an extreme case of the weak coupling regime, which we do not consider in
the present work, where $\lambdaambda \lambdaesssim $ $\Gammaamma _{\ell }$, the drift
coefficient $C_{\ell }$ depends mainly on $\alpha _{\ell }$ and the
diffusion coefficient is given by $D_{mn}=\partialtaelta _{mn}D_{mm}$, which is
strictly positive. In this particular case, the dynamics of the two
oscillators decouple, since the two-mode state relaxes, due to the damping
process, before they have time to interact. Thus
\begin{equation}gin{equation}
\frac{dP}{dt}=\sum_{\ell }\lambdaeft[ \Pi _{\ell }+C_{\ell }(\alpha _{\ell })
\frac{\partial }{\partial \alpha _{\ell }}+\frac{D_{\ell \ell }}{2}\frac{
\partial ^{2}}{\partial \alpha _{\ell }\partial \alpha _{\ell }^{\ast }}+H.c.
\right] P,  \lambdaabel{Eq30}
\end{equation}
and the general solution is simply the summation of the independent
solutions for each mode.

\subsection{Solution of the Fokker-Planck equation at absolute zero}

>From now on we assume both reservoirs set to absolute zero, such that $
N_{\ell }(\omegamega _{\ell }^{\pm })=0$. This assumption is equivalent to
setting $D_{mn}=0$, and thus the Fokker-Planck equation (\ref{Eq27}) reduces
to a simple drift equation
\begin{equation}gin{equation}
\frac{dP}{dt}=\sum_{\ell }\lambdaeft[ \Pi _{\ell }+C_{\ell }(\alpha _{1},\alpha
_{2})\frac{\partial }{\partial \alpha _{\ell }}+H.c.\right] P.  \lambdaabel{Eq31}
\end{equation}

With the substitution $P=P^{^\primeime }\exp \lambdaeft[ 2\lambdaeft( \Pi _{1}+\Pi
_{2}\right) t\right] $, we simplify the above equation to the form
\begin{equation}gin{equation}
\frac{dP^{^\primeime }}{dt}=\sum_{\ell }\lambdaeft[ C_{\ell }(\alpha _{1},\alpha _{2})
\frac{\partial }{\partial \alpha _{\ell }}+H.c.\right] P^{^\primeime }{\rm {,}}
\lambdaabel{Eq32}
\end{equation}
and assuming that $P^{^\primeime }(\alpha _{1},\alpha _{2},t)=P^{^\primeime }(\alpha
_{1}(t),\alpha _{2}(t))$,{\Lambdaarge \ }we get
\begin{equation}gin{equation}
\frac{d}{dt}P^{^\primeime }(\alpha _{1}(t),\alpha _{2}(t))=\sum_{\ell }\lambdaeft(
\frac{\partial \alpha _{\ell }}{\partial t}\frac{\partial }{\partial \alpha
_{\ell }}+H.c\right) P^{^\primeime }(\alpha _{1}(t),\alpha _{2}(t)){\rm {.}}
\lambdaabel{Eq33}
\end{equation}
Therefore, from Eqs. (\ref{Eq32}) and (\ref{Eq33}) we obtain the system of
coupled equations
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\frac{\partial \alpha _{1}}{\partial t} &=&B_{1}+\alpha _{1}E_{1}^{+}+\alpha
_{2}E_{1}^{-},  \lambdaabel{Eq34a} \\
\frac{\partial \alpha _{2}}{\partial t} &=&B_{2}+\alpha _{2}E_{2}^{+}+\alpha
_{1}E_{2}^{-},  \lambdaabel{Eq34b}
\end{eqnarray}
which leads to
\end{mathletters}
\begin{equation}gin{equation}
\frac{\partial ^{2}\alpha _{1}}{\partial t^{2}}-\lambdaeft(
E_{1}^{+}+E_{2}^{+}\right) \frac{\partial \alpha _{1}}{\partial t}+\alpha
_{1}\lambdaeft( E_{1}^{+}E_{2}^{+}-E_{1}^{-}E_{2}^{-}\right) +\lambdaeft(
B_{1}E_{2}^{+}-B_{2}E_{1}^{-}\right) =0.  \lambdaabel{Eq34l}
\end{equation}

Next, we define the parameters (where we have disregarded terms of order O$
^{2}$($\gammaamma _{\ell }^{\pm }/\lambdaambda $)){\LambdaARGE \ }$\Lambdaambda =\lambdaeft( \gammaamma
_{1}^{+}+\gammaamma _{2}^{+}+\gammaamma _{1}^{-}+\gammaamma _{2}^{-}\right) /4+i\Omegamega $
, $\Deltaelta =\lambdaeft( \gammaamma _{1}^{+}-\gammaamma _{2}^{+}+\gammaamma _{1}^{-}-\gammaamma
_{2}^{-}\right) /4$, $\Phi =\lambdaeft( \gammaamma _{1}^{+}+\gammaamma _{2}^{+}-\gammaamma
_{1}^{-}-\gammaamma _{2}^{-}\right) /4$, $\Thetaheta =\lambdaeft( \gammaamma _{1}^{+}-\gammaamma
_{2}^{+}-\gammaamma _{1}^{-}+\gammaamma _{2}^{-}\right) /4$, and the time-dependent
functions
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
W^{\pm } &=&\cosh \Phi t\lambdaeft( \cos \lambdaambda t\pm \frac{\Deltaelta }{\lambdaambda }
\sin \lambdaambda t\right) +i\sinh \Phi t\lambdaeft( \sin \lambdaambda t\mp \frac{\Deltaelta }{
\lambdaambda }\cos \lambdaambda t\right) ,  \lambdaabel{Eq35a} \\
Z^{\pm } &=&\lambdaeft( \sinh \Phi t\cos \lambdaambda t\pm \frac{\Thetaheta }{\lambdaambda }
\cosh \Phi t\sin \lambdaambda t\right) +i\lambdaeft( \cosh \Phi t\sin \lambdaambda t\mp
\frac{\Thetaheta }{\lambdaambda }\sinh \Phi t\cos \lambdaambda t\right) .  \lambdaabel{Eq35b}
\end{eqnarray}
Defining also the parameters
\end{mathletters}
\begin{equation}gin{eqnarray}
G_{1} &=&\frac{B_{2}E_{1}^{-}-B_{1}E_{2}^{+}}{
E_{1}^{+}E_{2}^{+}-E_{1}^{-}E_{2}^{-}},  \lambdaabel{Eq36a} \\
G_{2} &=&\frac{B_{1}E_{2}^{-}-B_{2}E_{1}^{+}}{
E_{1}^{+}E_{2}^{+}-E_{1}^{-}E_{2}^{-}},  \lambdaabel{Eq36b}
\end{eqnarray}
we finally obtain the solution of Eqs. (\ref{Eq34a}, \ref{Eq34b}), written
in the compact form
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\alpha _{1}(t) &=&

\mathop{\rm e}

{}^{\Lambdaambda t}\lambdaeft[ \lambdaeft( \alpha _{1}^{0}-G_{1}\right) W^{+}+\lambdaeft( \alpha
_{2}^{0}-G_{2}\right) Z^{+}\right] +G_{1}{\rm {,}}  \lambdaabel{Eq37a} \\
\alpha _{2}(t) &=&

\mathop{\rm e}

{}^{\Lambdaambda t}\lambdaeft[ \lambdaeft( \alpha _{2}^{0}-G_{2}\right) W^{-}+\lambdaeft( \alpha
_{1}^{0}-G_{1}\right) Z^{-}\right] +G_{2}{\rm {,}}  \lambdaabel{Eq37b}
\end{eqnarray}
where $\alpha _{\ell }^{0}$ indicates $\alpha _{\ell }(t=0)$. From the
solution of Eqs. (\ref{Eq37a}) and (\ref{Eq37b}) it follows that the
Fokker-Planck equation (\ref{Eq31}) can be solved to give
\end{mathletters}
\begin{equation}gin{equation}
P(\alpha _{1},\alpha _{2},t)=

\mathop{\rm e}

{}^{2\lambdaeft( \Pi _{1}+\Pi _{2}\right) t}P(\alpha _{1},\alpha _{2},0)|_{\alpha
_{\ell }\rightarrow \alpha _{\ell }(t)}{\rm {,}}  \lambdaabel{Eq38}
\end{equation}
where $P(\alpha _{1},\alpha _{2},0)$ is the P-function for the initial
state. The evolved P-function is then obtained by simply replacing $\alpha
_{\ell }$ by the evolved parameters given by Eqs. (\ref{Eq37a}) and (\ref
{Eq37b}). Finally, from Eq. (\ref{Eq38}) the evolved density operator
related to (\ref{Eq18}) follows from
\begin{equation}gin{equation}
\rho _{12}(t)=\int d^{2}\alpha _{1}d^{2}\alpha _{2}P(\alpha _{1},\alpha
_{2},t)|\alpha _{1},\alpha _{2}\rangle \lambdaangle \alpha _{1},\alpha _{2}|.
\lambdaabel{Eq39}
\end{equation}

\subsection{Initial joint states}

Next, we obtain the density operator $\rho _{12}(t)$ supposing that the
modes $1$ and $2$ are prepared in a superposition of coherent states of the
form
\begin{equation}gin{equation}
\lambdaeft| \varphi \right\rangle ={\cal N}_{\pm }\lambdaeft( \lambdaeft| \begin{equation}ta
_{I}^{1},\begin{equation}ta _{I}^{2}\right\rangle \pm \lambdaeft| \begin{equation}ta _{II}^{1},\begin{equation}ta
_{II}^{2}\right\rangle \right) ,  \lambdaabel{Eq40}
\end{equation}
where ${\cal N}_{\pm }$ stands for the normalization factor and the
subscripts $I,II$ are related to different coherent states participating in
the superposition. From Eqs. (\ref{Eq39}) and (\ref{Eq40}), the
corresponding evolved density operator is given by
\begin{equation}gin{equation}
\rho _{12}(t)=\sum_{m,n=I,II}C_{mn}(t)\lambdaeft| \varsigma _{m},\zetaeta
_{m}\right\rangle \lambdaeft\lambdaangle \varsigma _{n},\zetaeta _{n}\right| ,
\lambdaabel{Eq41}
\end{equation}
where the labels $\varsigma $ and $\zetaeta $ represent the states of modes $1$
and $2$, respectively. The expansion coefficients read
\begin{equation}gin{eqnarray}
C_{mn}(t) &=&{\cal N}_{\pm }^{2}\lambdaeft( \pm 1\right) ^{1-\partialtaelta
_{mn}}\lambdaeft\lambdaangle \begin{equation}ta _{n}^{1}\right. \lambdaeft| \begin{equation}ta _{m}^{1}\right\rangle
^{1-\exp \lambdaeft( -2{\rm Re}\lambdaeft\{ \Lambdaambda \right\} t\right) \lambdaeft( \lambdaeft|
W^{-}\right| ^{2}+\lambdaeft| Z^{-}\right| ^{2}\right) }  \nonumber \\
&&\thetaimes \lambdaeft\lambdaangle \begin{equation}ta _{n}^{2}\right. \lambdaeft| \begin{equation}ta
_{m}^{2}\right\rangle ^{1-\exp \lambdaeft( -2{\rm Re}\lambdaeft\{ \Lambdaambda \right\}
t\right) \lambdaeft( \lambdaeft| W^{+}\right| ^{2}+\lambdaeft| Z^{+}\right| ^{2}\right) }

\mathop{\rm e}

{}^{i(\thetaheta _{mn}^{(1)}+\thetaheta _{mn}^{(3)})+\thetaheta _{mn}^{(2)}+\thetaheta
_{mn}^{(4)}},  \lambdaabel{Eq42}
\end{eqnarray}
while the states $\varsigma $ and $\zetaeta $ are given by
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\varsigma _{\ell } &=&

\mathop{\rm e}

{}^{-\Lambdaambda t}\lambdaeft[ \lambdaeft( \begin{equation}ta _{\ell }^{1}-G_{1}\right) W^{-}-\lambdaeft(
\begin{equation}ta _{\ell }^{2}-G_{2}\right) Z^{+}\right] +G_{1},  \lambdaabel{Eq43a} \\
\zetaeta _{\ell } &=&

\mathop{\rm e}

{}^{-\Lambdaambda t}\lambdaeft[ \lambdaeft( \begin{equation}ta _{\ell }^{2}-G_{2}\right) W^{+}-\lambdaeft(
\begin{equation}ta _{\ell }^{1}-G_{1}\right) Z^{-}\right] +G_{2}  \lambdaabel{Eq43b}
\end{eqnarray}
We define the time-dependent functions
\end{mathletters}
\begin{equation}gin{mathletters}
\begin{equation}gin{eqnarray}
\thetaheta _{mn}^{(1)} &=&{\rm Im}\lambdaeft\{

\mathop{\rm e}

{}^{-\Lambdaambda ^{\ast }t}\lambdaeft[ \lambdaeft( W^{+}\right) ^{\ast }\lambdaeft( \begin{equation}ta
_{m}^{2\ast }-\begin{equation}ta _{n}^{2\ast }\right) -\lambdaeft( Z^{-}\right) ^{\ast }\lambdaeft(
\begin{equation}ta _{m}^{1\ast }-\begin{equation}ta _{n}^{1\ast }\right) \right] \right.  \nonumber \\
&&\thetaimes \lambdaeft. \lambdaeft[

\mathop{\rm e}

{}^{-\Lambdaambda t}\lambdaeft( G_{1}Z^{-}-G_{2}W^{+}\right) +G_{2}\right] \right\} ,
\lambdaabel{Eq44a} \\
2\thetaheta _{mn}^{(2)} &=&

\mathop{\rm e}

{}^{-2{\rm Re}\lambdaeft\{ \Lambdaambda \right\} t}\lambdaeft\{ W^{+}\lambdaeft( Z^{-}\right)
^{\ast }\lambdaeft[ \begin{equation}ta _{n}^{1\ast }\lambdaeft( \begin{equation}ta _{m}^{2}-\begin{equation}ta
_{n}^{2}\right) -\lambdaeft( \begin{equation}ta _{m}^{1\ast }-\begin{equation}ta _{n}^{1\ast }\right) \begin{equation}ta
_{m}^{2}\right] \right.  \nonumber \\
&&-\lambdaeft. \lambdaeft( W^{+}\right) ^{\ast }Z^{-}\lambdaeft[ \begin{equation}ta _{m}^{1}\lambdaeft( \begin{equation}ta
_{m}^{2\ast }-\begin{equation}ta _{n}^{2\ast }\right) -\lambdaeft( \begin{equation}ta _{m}^{1}-\begin{equation}ta
_{n}^{1}\right) \begin{equation}ta _{n}^{2\ast }\right] \right\} ,  \lambdaabel{Eq44b} \\
\thetaheta _{mn}^{(3)} &=&{\rm Im}\lambdaeft\{

\mathop{\rm e}

{}^{-\Lambdaambda ^{\ast }t}\lambdaeft[ \lambdaeft( W^{-}\right) ^{\ast }\lambdaeft( \begin{equation}ta
_{m}^{1\ast }-\begin{equation}ta _{n}^{1\ast }\right) -\lambdaeft( Z^{+}\right) ^{\ast }\lambdaeft(
\begin{equation}ta _{m}^{2\ast }-\begin{equation}ta _{n}^{2\ast }\right) \right] \right.  \nonumber \\
&&\thetaimes \lambdaeft. \lambdaeft[

\mathop{\rm e}

{}^{-\Lambdaambda t}\lambdaeft( G_{2}Z^{+}-G_{1}W^{-}\right) +G_{1}\right] \right\} ,
\lambdaabel{Eq44c} \\
2\thetaheta _{mn}^{(4)} &=&

\mathop{\rm e}

{}^{-2{\rm Re}\lambdaeft\{ \Lambdaambda \right\} t}\lambdaeft\{ W^{-}\lambdaeft( Z^{+}\right)
^{\ast }\lambdaeft[ \lambdaeft( \begin{equation}ta _{m}^{1}-\begin{equation}ta _{n}^{1}\right) \begin{equation}ta _{n}^{2\ast
}-\begin{equation}ta _{m}^{1}\lambdaeft( \begin{equation}ta _{m}^{2\ast }-\begin{equation}ta _{n}^{2\ast }\right) \right]
\right.  \nonumber \\
&&-\lambdaeft. \lambdaeft( W^{-}\right) ^{\ast }Z^{+}\lambdaeft[ \lambdaeft( \begin{equation}ta _{m}^{1\ast
}-\begin{equation}ta _{n}^{1\ast }\right) \begin{equation}ta _{m}^{2}-\begin{equation}ta _{n}^{1\ast }\lambdaeft( \begin{equation}ta
_{m}^{2}-\begin{equation}ta _{n}^{2}\right) \right] \right\} .  \lambdaabel{Eq44d}
\end{eqnarray}
Finally, the reduced density operator can be obtained easily from Eq. (\ref
{Eq41}), being
\end{mathletters}
\begin{equation}gin{eqnarray}
\rho _{1}(t) &=&{}{\cal N}_{\pm }^{2}\sum_{m,n=I,II}\lambdaeft( \pm 1\right)
^{1-\partialtaelta _{mn}}\lambdaeft\lambdaangle \begin{equation}ta _{n}^{1}\right. \lambdaeft| \begin{equation}ta
_{m}^{1}\right\rangle ^{1-\exp \lambdaeft( -2{\rm Re}\lambdaeft\{ \Lambdaambda \right\}
t\right) \lambdaeft| W^{-}\right| ^{2}}  \nonumber \\
&&\thetaimes \lambdaeft\lambdaangle \begin{equation}ta _{n}^{2}\right| \lambdaeft. \begin{equation}ta
_{m}^{2}\right\rangle ^{1-\exp \lambdaeft( -2{\rm Re}\lambdaeft\{ \Lambdaambda \right\}
t\right) \lambdaeft| Z^{+}\right| ^{2}}

\mathop{\rm e}

{}^{i\thetaheta _{mn}^{(3)}+\thetaheta _{mn}^{(4)}}\lambdaeft| \varsigma
_{m}\right\rangle \lambdaeft\lambdaangle \varsigma _{n}\right| {\rm {,}}  \lambdaabel{Eq45a}
\\
\rho _{2}(t) &=&{\cal N}_{\pm }^{2}\sum_{m,n=I,II}\lambdaeft( \pm 1\right)
^{1-\partialtaelta _{mn}}\lambdaeft\lambdaangle \begin{equation}ta _{n}^{1}\right| \lambdaeft. \begin{equation}ta
_{m}^{1}\right\rangle ^{1-\exp \lambdaeft( -2{\rm Re}\lambdaeft\{ \Lambdaambda \right\}
t\right) \lambdaeft| Z^{-}\right| ^{2}}  \nonumber \\
&&\thetaimes \lambdaeft\lambdaangle \begin{equation}ta _{n}^{2}\right. \lambdaeft| \begin{equation}ta
_{m}^{2}\right\rangle ^{1-\exp \lambdaeft( -2{\rm Re}\lambdaeft\{ \Lambdaambda \right\}
t\right) \lambdaeft| W^{+}\right| ^{2}}

\mathop{\rm e}

{}^{i\thetaheta _{mn}^{(1)}+\thetaheta _{mn}^{(2)}}\lambdaeft| \zetaeta _{m}\right\rangle
\lambdaeft\lambdaangle \zetaeta _{n}\right| .  \lambdaabel{Eq45b}
\end{eqnarray}

\section{State recurrence and swap dynamics}

Let us suppose that the joint system $1+2$ is prepared in the superposition
state
\begin{equation}gin{equation}
|\psi _{12}\rangle ={}{\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha
\rangle \right) _{1}\omegatimes |\eta \rangle _{2}{\rm {,}}  \lambdaabel{Eq46}
\end{equation}
which represents a particular case of Eq. (\ref{Eq40}), where $\begin{equation}ta
_{I}^{1}=-\begin{equation}ta _{II}^{1}=\alpha $ and $\begin{equation}ta _{I}^{2}=\begin{equation}ta _{II}^{2}=\eta
. $

In this section we analyze, in the weak and strong coupling regimes,\ the
effect of dissipation on two phenomena: $i)$ the joint-system state
recurrence and $ii)$ the state swap between the systems. Considering the
dynamics of the coupled systems, governed by the Fokker-Planck equation (\ref
{Eq31}), we calculate $i)$ the probability that in each system the initial
state recurs -- the recurrence time for both systems being the same -- and $
ii)$ the probability of a swapping of states between the systems, i.e., the
probability\ of oscillator $1$ being in the initial state of oscillator $2$,
and vice-versa. Recurrence is guaranteed whenever a dynamical system is
closed and its evolution is limited \cite{Bocchieri}, while state swap is
only possible for specific interactions. Here the two oscillators are
coupled by a bilinear Hamiltonian which, in principle, allows state swap to
occur between the systems \cite{meu}.

Considering the initial superposition (\ref{Eq46}), the time-evolved joint
state following from the density operator (\ref{Eq41}) reads
\begin{equation}gin{equation}
\rho _{12}(t)=\sum_{m,n=I,II}{\cal C}_{mn}(t)\lambdaeft| \varsigma _{m},\zetaeta
_{m}\right\rangle \lambdaeft\lambdaangle \varsigma _{n},\zetaeta _{n}\right| {\rm {,}}
\lambdaabel{Eq47}
\end{equation}
where
\begin{equation}gin{eqnarray}
{\cal C}_{mn}(t) &=&{\cal N}_{\pm }^{2}\lambdaeft( \pm 1\right) ^{1-\partialtaelta
_{mn}}\exp \lambdaeft\{ -2\lambdaeft| \alpha \right| ^{2}\lambdaeft[ 1-\lambdaeft( \lambdaeft|
W^{-}\right| ^{2}+\lambdaeft| Z^{-}\right| ^{2}\right) \exp \lambdaeft( -2{\rm Re}
\lambdaeft\{ \Lambdaambda \right\} t\right) \right] \right\}   \nonumber \\
&&\thetaimes \exp \lambdaeft\{ i\lambdaeft[ \thetaheta _{mn}^{(1)}+\thetaheta _{mn}^{(3)}+

\mathop{\rm Im}

\lambdaeft( \thetaheta _{mn}^{(2)}+\thetaheta _{mn}^{(4)}\right) \right] \right\} {}{\rm {
,}}  \lambdaabel{Eq48}
\end{eqnarray}
and the probability of recurrence is given by
\begin{equation}gin{eqnarray}
P_{R}(t) &\equiv &Tr_{12}\lambdaeft[ \rho _{12}(t)\rho _{12}(0)\right]   \nonumber
\\
&=&{\cal N}_{\pm }^{2}\sum_{m,n=I,II}{\cal C}_{mn}(t)\lambdaeft\lambdaangle \varsigma
_{n}\right| \lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{1}\lambdaeft(
\lambdaeft\lambdaangle \alpha \right| \pm \lambdaeft\lambdaangle -\alpha \right| \right) \lambdaeft|
\varsigma _{m}\right\rangle \lambdaeft\lambdaangle \zetaeta _{m}\right| \lambdaeft. \eta
\right\rangle _{2}\lambdaeft\lambdaangle \eta \right| \lambdaeft. \zetaeta _{m}\right\rangle
{\rm {.}}  \lambdaabel{Eq49}
\end{eqnarray}
Analogously to the recurrence probability, the state-swap probability is
calculated by swapping the labels of each initial state of the two fields,
as defined in \cite{meu}:
\begin{equation}gin{eqnarray}
P_{S}(t) &\equiv &Tr\lambdaeft[ \rho _{12}(t)\rho _{12}(0)|_{1(2)\rightarrow 2(1)}
\right]   \nonumber \\
&=&{\cal N}_{\pm }^{2}\sum_{m,n=I,II}{\cal C}_{mn}(t)\lambdaeft\lambdaangle \varsigma
_{n}\right| \lambdaeft. \eta \right\rangle _{2}\lambdaeft\lambdaangle \eta \right. \lambdaeft|
\varsigma _{m}\right\rangle \lambdaeft\lambdaangle \zetaeta _{m}\right| \lambdaeft( |\alpha
\rangle \pm |-\alpha \rangle \right) _{1}\lambdaeft( \lambdaeft\lambdaangle \alpha \right|
\pm \lambdaeft\lambdaangle -\alpha \right| \right) \lambdaeft| \zetaeta _{m}\right\rangle {\rm
{.}}  \lambdaabel{Eq50}
\end{eqnarray}

In Fig. 4 we plot the state-swap probability $P_{S}(t)$ (dotted line) and
the recurrence probability $P_{R}(t)$ (solid line) against the scaled time $
\lambdaambda t$, taking $\alpha =\eta =1$ as real parameters. In Figs. 3(a, b and
c) we set the driving field strength $F$ to zero (so that $\omegamega
_{10}=\omegamega _{20}$ and the minimum of $\omegamega _{\ell }^{-}$, occurring for $
\lambdaambda =2\omegamega _{20}$, also becomes zero) and assume the absence of
dissipation, so as to take these figures as references. In Fig. 4(a) we
consider the weak coupling regime, assuming $\lambdaambda /\omegamega _{10}=2\thetaimes
10^{-2}$ (a somewhat exaggerated ratio used to visualize better the strong
oscillations of probabilities $P_{R}(t)$ and $P_{S}(t)$). We observe (dotted
line), as can be deduced from Eq. (\ref{Eq50}), that the modes will swap
their states whenever
\begin{equation}gin{equation}
\lambdaambda t_{S}=(2n+1)\frac{\pi }{2},{\rm {\ }}n=0,1,2,...,  \lambdaabel{Eq51}
\end{equation}
which means that in the swap time $t_{S}$, the state of mode $1$ becomes $
|\begin{equation}ta \rangle _{1}$, while the state of mode $2$ becomes the superposition $
\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{2}$. From the solid
line and Eq. (\ref{Eq49}), we observe that the joint-system state recurs
whenever
\begin{equation}gin{equation}
\lambdaambda t_{R}=n\pi ,{\rm {\ }}n=0,1,2,...,  \lambdaabel{Eq52}
\end{equation}
i.e., in the recurrence time $t_{R}$ the mode $1$($2$) returns to its
initial state $\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{1}$ ($
|\begin{equation}ta \rangle _{2}$). It is evident that the shape of Fig. 4(a) results
mainly from the small value of $\lambdaambda $\ (compared to $\Omegamega $) which
defines an envelope function.

In Fig. 4(b) an intermediate coupling is assumed, such that $\lambdaambda /\omegamega
_{10}=1$. As can be observed, the recurrence dynamics remains the same
(solid line), while the swap dynamics begins to be affected by the coupling
strength. In Fig. 4(c) we consider the strong coupling regime where $\lambdaambda
=2\omegamega _{10}$ ($\omegamega _{10}=\omegamega _{20}$), such that $\Omegamega /\lambdaambda =1$
(note that $\Omegamega /\lambdaambda =\omegamega _{10}/\lambdaambda +\lambdaambda /4\omegamega _{20}$).
As in Fig. 4(a), in Fig. 4(c) the recurrence process remains unchanged,
still obeying Eq. (\ref{Eq52}), so that the recurrence time becomes smaller
due to the strong coupling parameter $\lambdaambda $. However, the swap dynamics
is almost completely lost, the remaining oscillations arising from the
nonorthogonality between the states $\lambdaeft( |\alpha \rangle \pm |-\alpha
\rangle \right) _{1}$ and $|\begin{equation}ta \rangle _{2}$.

To understand the behavior of $P_{R}(t)$ and $P_{S}(t)$ in Figs. 4(b and c),
we note that in the absence of dissipation the recurrence and state-swap
probabilities are computed as $P_{R}(t)=\lambdaeft| \lambdaeft\lambdaangle \psi
_{12}(0)\right. \lambdaeft| \psi _{12}(t)\right\rangle \right| ^{2}$ and $
P_{S}(t)=\lambdaeft| \lambdaeft\lambdaangle \psi _{12}(0)\right. \lambdaeft| \lambdaeft. \psi
_{12}(t)\right| _{1(2)\rightarrow 2(1)}\right\rangle \right| ^{2}$,
respectively. For $P_{R}(t)$, the probability amplitude
\begin{equation}gin{eqnarray}
\lambdaeft\lambdaangle \psi _{12}(0)\right. \lambdaeft| \psi _{12}(t)\right\rangle &=&{\cal 
N}_{\pm }^{2}\lambdaeft( \lambdaeft\lambdaangle \alpha \right. \lambdaeft| \varsigma
_{I}\right\rangle \lambdaeft\lambdaangle \eta \right. |\zetaeta _{I}\rangle +\exp \lambdaeft(
2i\alpha G_{2}\sin \lambdaeft( \lambdaambda t\right) \right) \lambdaeft\lambdaangle -\alpha
\right. \lambdaeft| \varsigma _{II}\right\rangle \lambdaeft\lambdaangle \eta \right. |\zetaeta
_{II}\rangle \right.  \nonumber \\
&&\lambdaeft. \pm \lambdaeft\lambdaangle -\alpha \right. \lambdaeft| \varsigma _{I}\right\rangle
\lambdaeft\lambdaangle \eta \right. |\zetaeta _{I}\rangle \pm \exp \lambdaeft( 2i\alpha
G_{2}\sin \lambdaeft( \lambdaambda t\right) \right) \lambdaeft\lambdaangle \alpha \right. \lambdaeft|
\varsigma _{II}\right\rangle \lambdaeft\lambdaangle \eta \right. |\zetaeta _{II}\rangle
\right) ,  \lambdaabel{Eq53}
\end{eqnarray}
tends to unity under the conditions $\zetaeta _{I},\zetaeta _{II}\rightarrow \eta $
, $\varsigma _{I}\rightarrow \alpha $ ($-\alpha $), and $\varsigma
_{II}\rightarrow -\alpha $ ($\alpha $), which are satisfied when $\cos
\Omegamega t_{R}\cos \lambdaambda t_{R}=1$. This relation implies that
\begin{equation}gin{equation}
t_{R}=

{\partialtaisplaystyle{n\pi  \omegaver \lambdaambda }}

=\frac{m\pi }{\Omegamega },  \lambdaabel{Eq54}
\end{equation}
where $n,m$ are integers, both being even or odd. In Fig. 4(b), the
intermediate coupling $\lambdaambda /\omegamega _{10}=1$ follows from $n=(4/5)m$,
such that the joint-system state recurs whenever
\begin{equation}gin{equation}
\lambdaambda t_{R}=8n\pi ,n=0,1,2,...,  \lambdaabel{Eq55}
\end{equation}
a result which explains the shift in the scaled time $\lambdaambda t$ observed in
Fig. 4(b).

Concerning the state-swap probability $P_{S}(t)$, we obtain the expression
\begin{equation}gin{eqnarray}
\lambdaeft\lambdaangle \psi _{12}(0)|_{1(2)\rightarrow 2(1)}\right. \lambdaeft| \psi
_{12}(t)\right\rangle &=&{\cal N}_{\pm }^{2}\lambdaeft( \lambdaeft\lambdaangle \alpha
\right. \lambdaeft| \zetaeta _{I}\right\rangle \lambdaeft\lambdaangle \eta \right. |\varsigma
_{I}\rangle +\exp \lambdaeft( 2iG_{2}\alpha \sin \lambdaeft( \lambdaambda t\right) \right)
\lambdaeft\lambdaangle -\alpha \right. \lambdaeft| \zetaeta _{II}\right\rangle \lambdaeft\lambdaangle
\eta \right. |\varsigma _{II}\rangle \right.  \nonumber \\
&&\lambdaeft. \pm \lambdaeft\lambdaangle -\alpha \right. \lambdaeft| \zetaeta _{I}\right\rangle
\lambdaeft\lambdaangle \eta \right. |\varsigma _{I}\rangle \pm \exp \lambdaeft(
2iG_{2}\alpha \sin \lambdaeft( \lambdaambda t\right) \right) \lambdaeft\lambdaangle \alpha
\right. \lambdaeft| \zetaeta _{II}\right\rangle \lambdaeft\lambdaangle \eta \right. |\varsigma
_{II}\rangle \right) ,  \lambdaabel{Eq56}
\end{eqnarray}
which tends to unity when $\varsigma _{I},\varsigma _{II}\rightarrow \eta $,
$\zetaeta _{I}\rightarrow \alpha $ ($-\alpha $), and $\zetaeta _{II}\rightarrow
-\alpha $ ($\alpha $). These conditions are satisfied only when $\sin \Omegamega
t_{S}\sin \lambdaambda t_{S}=-1${\LambdaARGE \ }and{\LambdaARGE \ }$G_{1}$,$G_{2}\approx 0$
. For the special case of $F=0$, as in Fig. 4(b and c), the condition $G_{1}$
,$G_{2}\approx 0$ is automatically satisfied (as can be deduced from Eqs. (
\ref{Eq36a}) and (\ref{Eq36b})) and the modes will swap their states
whenever
\begin{equation}gin{equation}
t_{S}=\frac{\lambdaeft( 2n-1\right) \pi }{2\lambdaambda }=\frac{\lambdaeft( 2m+1\right) \pi
}{2\Omegamega }{\rm {,}}  \lambdaabel{Eq57}
\end{equation}
with the additional condition that $\sin \lambdaeft[ \lambdaeft( n-1/2\right) \pi 
\right] =-\sin \lambdaeft[ \lambdaeft( m+1/2\right) \pi \right] $. Therefore, the
{\lambdaarge \ }maxima of the swap probability are eliminated, as observed in
Fig. 4(b) (for the intermediate coupling $\lambdaambda /\omegamega _{10}=1$) due to
the fact that the relation $2n-1=4\lambdaeft( 2m+1\right) /5$ (and so $\sin
\Omegamega t_{S}\sin \lambdaambda t_{S}=-1$) cannot be satisfied for any pair ($n,m$
). However, owing to the oscillations of $P_{S}(t)$ within the envelop
function defined by $\lambdaambda $($=4\Omegamega /5$), the neighborhood of these
maxima still survives. Differently, in Fig. 4(c) (for the strong coupling
regime $\Omegamega /\lambdaambda =1$), the relation $\sin \Omegamega t_{S}\sin \lambdaambda
t_{S}=\sin ^{2}\lambdaambda t_{S}\neq -1$ implies that the single maximum of the
swap probability is eliminated (since $\Omegamega /\lambdaambda =1$, there are no
oscillations of $P_{S}(t)$ but the envelope function). A heuristic
explanation of the behaviors of Figs. 4(b and c) will be provided below. In
fact, as discussed below, while the superposition state ${\cal N}_{\pm
}\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{1}$ swaps to
oscillator $2$, the coherent state $|\eta \rangle _{2}$ does not swap to
oscillator $1$, even though both return to their respective systems.
Evidently, the phase mismatching between the coupling parameter $\lambdaambda $
and the field-shifted frequencies $\omegamega _{\ell }=\Omegamega $ (when $F=0$),
represented by the relation $\sin \Omegamega t_{S}\sin \lambdaambda t_{S}=\sin
^{2}\lambdaambda t_{S}\neq -1$, is the core of this result.

In Fig. 4(d), also setting $F=0$, we include dissipation, taking both
systems with the same damping rate $\gammaamma _{\ell }^{\pm }=\Gammaamma /2$, where
$\Gammaamma /\omegamega _{10}=2\thetaimes 10^{-3}$, and $\lambdaambda /\omegamega _{10}=2\thetaimes
10^{-2}$ (as in the weak coupling regime) and observing, as expected, the
relaxation of both probabilities $P_{R}(t)$ (solid line) and $P_{S}(t)$
(dotted line). To compare the relaxation of probabilities $P_{R}(t)$ and $
P_{S}(t)$ in both regimes, in Fig. 4(e) we assume the parameters $F=0$, $
\gammaamma _{\ell }^{+}=2\gammaamma _{\ell }^{-}=\Gammaamma /2$ (assuming Markovian
white noise$)$, $\Gammaamma /\omegamega _{10}=2\thetaimes 10^{-3}$, and $\lambdaambda /\omegamega
_{10}=2$ (as in the strong coupling regime). We observe, comparing Figs.
4(d) and 4(e), that in the strong coupling regime the fields recur more
frequently, within the relaxation time, than in the weak coupling regime.
Owing to the strong coupling, in Fig. 4(e) the swap dynamics is almost
completely lost, as in Fig. 4(c). In Fig. 4(e), a dashed-dotted line has
been drawn at unity in order to display the slow decay of $P_{R}(t)$.

Finally, in Fig. 4(f) we again disregard dissipation, but turn on the
driving field, taking $F/\omegamega _{20}=1$ (i.e., $\omegamega _{10}=2\omegamega _{20}$
), and $\lambdaambda /\omegamega _{10}=2\thetaimes 10^{-2}$ (as in the weak coupling
regime), and observe a reduction of the swap probability compared to the
recurrence probability. (Note that the driving field occurs in the von
Neumann term of Eq. (\ref{Eq18}) even though it does not influence the
Liouville operator.) This behavior can be explained by the result
\begin{equation}gin{equation}
G_{1}-G_{2}=\frac{F}{\Omegamega -\lambdaambda }{\rm {,}}  \lambdaabel{Eq58}
\end{equation}
which becomes $1/2$\ for the parameters considered above, preventing the
swap probability from being unity, as required by the conditions $\sin
\Omegamega t_{S}\sin \lambdaambda t_{S}=-1${\LambdaARGE \ }and{\LambdaARGE \ }$
G_{1}-G_{2}\approx 0$. We have assumed in this figure the ratio $\omegamega
/\omegamega _{10}=10^{-2}$, instead of the resonance condition $\omegamega /\omegamega
_{10}=1/2$ adopted above, in order to make clear the amplification effects.
>From Eq. (\ref{Eq56}) it follows immediately that the probability $P_{S}(t)$
is reduced by the factor $\exp \lambdaeft\{ -2\lambdaeft[ F/(\Omegamega -\lambdaambda )\right]
^{2}\right\} $, which explains the maximum value around $1/2$ for the
probability $P_{S}(t)$. It is worth noting that the expression (\ref{Eq58})
does not diverge for $F\neq 0$, since in the case where $\Omegamega =\lambdaambda $
(and switching off the dissipation), Eq. (\ref{Eq34l}) becomes

\begin{equation}gin{equation}
\frac{\partial ^{2}\alpha _{1}}{\partial t^{2}}-2\Omegamega \frac{\partial
\alpha _{1}}{\partial t}-iF\lambdaambda =0,  \lambdaabel{Eq58l}
\end{equation}
whose solution differs from that leading to Eq. (\ref{Eq58}).

It can be argued that the state-swap and recurrence dynamics are
consequences of energy transfer between modes, which in the presence of
dissipation is severely reduced, as the two modes tend to reach energy
equilibrium with the reservoirs and the driving field. Actually, state-swap
and recurrence are properties of information transfer rather than energy
transfer between systems, even though these quantities are generally
correlated. However, it can be shown \cite{meu} that for the coupling
between the modes selected above, even if each mode is kept to constant
energy, their states can be swapped, showing the independence of the two
processes. In the presence of the relaxation process, the information
transfer between modes is reduced, becouse of the absorption of information
by the reservoirs. When the systems are strongly interacting, however, the
field states recur many times on the scale of $\lambdaambda $ before the
relaxation takes place. While not preventing the modes from relaxing to
equilibrium with the reservoirs, the recurrence has striking consequences
for short-time-scale dynamics, such as the dynamics of decoherence discussed
below.

\section{coherence properties}

So far we have analyzed the dynamics of strongly-interacting modes $1$ and $
2 $ to discuss the recurrence and state swap processes. Now we analyze, also
in the strong coupling regime,\ with the driving field switched off ($F=0$)
and the reservoirs at absolute zero ($T=0$), the decoherence dynamics of the
joint state described by the density operator $\rho _{12}$, and of the state
of system $1$ ($2$), described by $\rho _{1}=

\mathop{\rm Tr}

_{2}\rho _{12}$ ($\rho _{2}=

\mathop{\rm Tr}

_{1}\rho _{12}$). In this section we also consider the case of identical
dissipative systems, $\gammaamma _{\ell }^{\pm }=\gammaamma ^{\pm }$. The case of
different decay rates will be analyzed subsequently. We consider three
different initial joint states: first, the disentangled state, given by Eq. (
\ref{Eq46}), and then the entanglements which are eigenstates associated
with the normal modes $\omegamega _{\ell }^{\pm }$, derived from Eq. (\ref{Eq40}
): with $\begin{equation}ta _{I}^{1}=\begin{equation}ta _{II}^{2}=\alpha $\ and $\begin{equation}ta _{I}^{2}=\begin{equation}ta
_{II}^{1}=-\alpha $ we obtain the eigenstate
\begin{equation}gin{equation}
\lambdaeft| \varphi _{12}^{-}\right\rangle ={\cal N}_{\pm }\lambdaeft( \lambdaeft| \alpha
,-\alpha \right\rangle _{12}\pm \lambdaeft| -\alpha ,\alpha \right\rangle
_{12}\right) {,}  \lambdaabel{Eq59}
\end{equation}
associated with the normal mode $\omegamega _{\ell }^{-}$, and with $\begin{equation}ta
_{I}^{1}=\begin{equation}ta _{I}^{2}=\alpha $\ and $\begin{equation}ta _{II}^{2}=\begin{equation}ta
_{II}^{1}=-\alpha $, we obtain
\begin{equation}gin{equation}
\lambdaeft| \varphi _{12}^{+}\right\rangle ={\cal N}_{\pm }\lambdaeft( \lambdaeft| \alpha
,\alpha \right\rangle _{12}\pm \lambdaeft| -\alpha ,-\alpha \right\rangle
_{12}\right) {,}  \lambdaabel{Eq60}
\end{equation}
which is the eigenstate associated with $\omegamega _{\ell }^{+}$.

\subsection{Decoherence time of the state{\it \ }$|^\primeotect\psi _{12}\rangle
={}{\cal N}_{\pm }\lambdaeft( |^\primeotect\alpha \rangle \pm |-^\primeotect\alpha 
\rangle \right) _{1}\omegatimes |^\primeotect\eta \rangle _{2}$}

The coherence of the joint state (\ref{Eq46}) is given, essentially, by the
term
\begin{equation}gin{equation}
\exp \lambdaeft\{ -2|\alpha |^{2}\lambdaeft[ 1-\lambdaeft( \lambdaeft| W^{-}\right| ^{2}+\lambdaeft|
Z^{-}\right| ^{2}\right) \exp \lambdaeft( -2{\rm Re}\lambdaeft\{ \Lambdaambda \right\}
t\right) \right] \right\} ,  \lambdaabel{Eq61}
\end{equation}
coming from the off-diagonal coefficients of the density operator in Eq. (
\ref{Eq47}). In the strong coupling regime, the exponential decay in Eq. (
\ref{Eq61}), computed from Eqs. (\ref{Eq35a}), (\ref{Eq35b}), reduces to the
form
\begin{equation}gin{equation}
\exp \lambdaeft[ -|\alpha |^{2}\lambdaeft( 2-

\mathop{\rm e}

\nolimits^{-\lambdaeft( \gammaamma _{1}^{+}+\gammaamma _{2}^{+}\right) t}-

\mathop{\rm e}

\nolimits^{-\lambdaeft( \gammaamma _{1}^{-}+\gammaamma _{2}^{-}\right) t}\right) \right] {
.}  \lambdaabel{Eq62}
\end{equation}
Assuming identical dissipative systems, $\gammaamma _{\ell }^{\pm }=\gammaamma ^{\pm
}$, Eq. (\ref{Eq62}) simplifies to
\begin{equation}gin{equation}
\exp \lambdaeft[ -|\alpha |^{2}\lambdaeft( 2-

\mathop{\rm e}

\nolimits^{-2\gammaamma ^{+}t}-

\mathop{\rm e}

\nolimits^{-2\gammaamma ^{-}t}\right) \right] {,}  \lambdaabel{Eq63}
\end{equation}
resulting in a decoherence time for the joint system given by
\begin{equation}gin{equation}
\thetaau _{D}=\lambdaeft[ 2|\alpha |^{2}\lambdaeft( \gammaamma ^{+}+\gammaamma ^{-}\right) \right]
^{-1}{\rm {.}}  \lambdaabel{Eq63l}
\end{equation}
This decoherence time has to be compared with that for an isolated mode in
the superposition state ${\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha
\rangle \right) $ (also obtained from the weak coupling regime: $\gammaamma
^{\pm }=\Gammaamma /2$), following from the exponential decay
\begin{equation}gin{equation}
\exp \lambdaeft[ -2|\alpha |^{2}\lambdaeft( 1-

\mathop{\rm e}

\nolimits^{-\Gammaamma t}\right) \right]  \lambdaabel{Eq64}
\end{equation}
and given by the well-known expression $\thetaau _{D}\approx \thetaau _{R}/2|\alpha
|^{2}\equiv {\cal T}_{D}$, where $\thetaau _{R}=\Gammaamma ^{-1}$ stands for the
relaxation time of the system. Analyzing the decoherence process in the
light of the spectral densities considered in Section III, we observe that
for Markovian white noise ($M$), where $\gammaamma ^{-}=\gammaamma ^{+}/2=\Gammaamma /4$
, the decoherence time in Eq. (\ref{Eq63l}) rises to
\begin{equation}gin{equation}
\thetaau _{D}^{M}\approx 4{\cal T}_{D}/3.  \lambdaabel{EqTD1}
\end{equation}
For the Lorentzian spectral density ($L$), where $\gammaamma ^{\pm }=\varepsilon
^{\pm }\Gammaamma /2\lambdal \Gammaamma $, we obtain the result
\begin{equation}gin{equation}
\thetaau _{D}^{L}\approx 2{\cal T}_{D}/(\varepsilon ^{+}+\varepsilon ^{-}),
\lambdaabel{EqTD2}
\end{equation}
which becomes large as the reservoir spectral density decreases around the
effective frequency $\omegamega _{\ell }^{\pm }=\omegamega ^{\pm }$. We stress that
for the case where $\varepsilon ^{+}=\varepsilon ^{-}$\ (i.e., $\gammaamma
^{+}=\gammaamma ^{-}$) the cross-decay channel is null.

Finally, for a wide Lorentzian spectral density ($WL$), where $\gammaamma
^{-}=\varepsilon ^{-}\Gammaamma /2\lambdal \Gammaamma $ and $\gammaamma ^{+}=\Gammaamma /2$, we
obtain from Eq. (\ref{Eq63l}), the value
\begin{equation}gin{equation}
\thetaau _{D}^{WL}\approx 2{\cal T}_{D}.  \lambdaabel{EqTD3}
\end{equation}
Therefore, for a wide Lorentzian spectral density we obtain a decoherence
time for strongly coupled systems which is twice as long as in the
weak-coupling regime. The mechanism behind these improved decoherence times
is that in the strong coupling regime (where the\ natural frequency $\omegamega
_{\ell 0}$ of the systems splits into two effective frequencies $\omegamega
_{\ell }^{\pm }=\omegamega _{\ell }\pm \lambdaambda $) the spectral density of the
reservoir plays a decisive role in the damping rate, which also splits into
two Lorentzian functions. When the oscillator effective frequency $\omegamega
_{\ell }^{\pm }$ shifts to regions where the spectral density of the
reservoir is significantly smaller than that around $\omegamega _{\ell 0}$, the
damping rate becomes smaller than its value in the weak coupling regime.
Despite the spectral densities, the competition between the cross-decay and
the usual channels is the reason for the computed delay of the decoherence
process.

Still regarding state (\ref{Eq46}) and the strong coupling regime, focusing
on the reduced systems, we obtain from Eqs. (\ref{Eq45a}) and (\ref{Eq45b})
the density operators

\begin{equation}gin{eqnarray}
\rho _{1}(t) &=&{\cal N}_{\pm }^{2}\sum_{m,n=I,II}\lambdaeft( \pm 1\right)
^{1-\partialtaelta _{mn}}\lambdaeft\lambdaangle \begin{equation}ta _{n}^{1}\right. \lambdaeft| \begin{equation}ta
_{m}^{1}\right\rangle ^{1-{\cal F}(\lambdaambda t)/4}

\mathop{\rm e}

{}^{i\thetaheta _{mn}^{(3)}+\thetaheta _{mn}^{(4)}}\lambdaeft| \varsigma
_{m}\right\rangle \lambdaeft\lambdaangle \varsigma _{n}\right| ,  \lambdaabel{Eq65a} \\
\rho _{2}(t) &=&{\cal N}_{\pm }^{2}\sum_{m,n=I,II}\lambdaeft( \pm 1\right)
^{1-\partialtaelta _{mn}}\lambdaeft\lambdaangle \begin{equation}ta _{n}^{1}\right. \lambdaeft| \begin{equation}ta
_{m}^{1}\right\rangle ^{1-{\cal F}(\lambdaambda t+\pi /2)/4}

\mathop{\rm e}

{}^{i\thetaheta _{mn}^{(1)}+\thetaheta _{mn}^{(2)}}\lambdaeft| \zetaeta _{m}\right\rangle
\lambdaeft\lambdaangle \zetaeta _{n}\right| ,  \lambdaabel{Eq65b}
\end{eqnarray}
where
\begin{equation}gin{equation}
{\cal F}(\lambdaambda t)=\exp \lambdaeft[ -(\gammaamma _{1}^{+}+\gammaamma _{2}^{+})t\right]
+\exp \lambdaeft[ -(\gammaamma _{1}^{-}+\gammaamma _{2}^{-})t\right] +2\cos (2\lambdaambda
t)\exp \lambdaeft[ -(\gammaamma _{1}^{+}+\gammaamma _{2}^{+}+\gammaamma _{1}^{-}+\gammaamma
_{2}^{-})t/2\right] .  \lambdaabel{Eq66}
\end{equation}
Observe that for identical systems, $\gammaamma _{\ell }^{\pm }=\gammaamma ^{\pm }$,
the decoherence of the reduced state in oscillator$1$ is given by
\begin{equation}gin{equation}
\exp \lambdaeft\{ -2|\alpha |^{2}\lambdaeft[ 1-\frac{1}{4}\lambdaeft( {\rm {

\mathop{\rm e}

}}^{-2\gammaamma ^{+}t}+{\rm {

\mathop{\rm e}

}}^{-2\gammaamma ^{-}t}+2\cos (2\lambdaambda t){\rm {

\mathop{\rm e}

}}^{-(\gammaamma ^{+}+\gammaamma ^{-})t}\right) \right] \right\}  \lambdaabel{Eq67a}
\end{equation}
while for the state in oscillator $2$ it is
\begin{equation}gin{equation}
\exp \lambdaeft\{ -2|\alpha |^{2}\lambdaeft[ 1-\frac{1}{4}\lambdaeft( {\rm {

\mathop{\rm e}

}}^{-2\gammaamma ^{+}t}+{\rm {

\mathop{\rm e}

}}^{-2\gammaamma ^{-}t}+2\cos (2\lambdaambda t+\pi ){\rm {

\mathop{\rm e}

}}^{-(\gammaamma ^{+}+\gammaamma ^{-})t}\right) \right] \right\} .  \lambdaabel{Eq67b}
\end{equation}
It is easily shown that the expressions in Eqs. (\ref{Eq67a}) and (\ref
{Eq67b}), associated with the decoherence times of the states of systems $1$
and $2$, respectively, oscillate below the curve for the coherence decay of
the joint system given by expression (\ref{Eq63}). In Fig. 5, assuming
Markovian white noise and setting the fictitious ratio $\lambdaambda /\Gammaamma =5$
to make the oscillations clear, the dashed and dotted lines refer to the
decoherence dynamics of systems $1$ and $2$, dictated by Eqs. (\ref{Eq67a})
and (\ref{Eq67b}), respectively. The solid line represents the coherence
decay derived from Eq. (\ref{Eq63}) and the dashed-dotted line indicates the
coherence decay for an isolated mode, computed from Eq. (\ref{Eq64}).
Therefore, for the special case of the initially disentangled state given by
Eq. (\ref{Eq46}), the decoherence times of systems $1$ and $2$ coincide with
that of the joint system.

\subsection{Decoherence time of the state{\it \ }$\lambdaeft| ^\primeotect\varphi 
_{12}^{-}\right\rangle ={\cal N}_{\pm }\lambdaeft( \lambdaeft| ^\primeotect\alpha ,-
^\primeotect\alpha \right\rangle _{12}\pm \lambdaeft| -^\primeotect\alpha ,^\primeotect\alpha 
\right\rangle _{12}\right) $}

Next, we analyze the decoherence dynamics of the entangled state (\ref{Eq59}
) for the joint system. The decoherence process of this joint state, in the
strong coupling regime, is given by the exponential decay
\begin{equation}gin{equation}
\exp \lambdaeft[ -4\lambdaeft| \alpha \right| ^{2}\lambdaeft( 1-

\mathop{\rm e}

\nolimits^{-\lambdaeft( \gammaamma _{1}^{-}+\gammaamma _{2}^{-}\right) t}\right) \right] ,
\lambdaabel{Eq69}
\end{equation}
which takes into account only the decay rate of the system-reservoir
coupling around the effective frequency $\omegamega _{\ell }^{-}$. At this point
it is interesting to identify the contribution of the cross-decay channel to
the decoherence process, rewriting Eq. (\ref{Eq69}) as
\begin{equation}gin{equation}
\exp \lambdaeft\{ -4\lambdaeft| \alpha \right| ^{2}\lambdaeft( 1-\exp \lambdaeft[ -

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}+\gammaamma _{\ell }^{-}\right) t/2
\right] \exp \lambdaeft[

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) t/2
\right] \right) \right\} .  \lambdaabel{Eq(69l)}
\end{equation}
In this expression, the term $\exp \lambdaeft[

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) t/2
\right] $ comes from the cross-decay channel and it is evident that for $

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) >0$,
its contribution makes the exponential decay slower. Assuming identical
dissipative systems, $\gammaamma _{\ell }^{-}={\gammaamma }^{-}$, the expression (
\ref{Eq69}) simplifies to
\begin{equation}gin{equation}
\exp \lambdaeft[ -4\lambdaeft| \alpha \right| ^{2}\lambdaeft( 1-

\mathop{\rm e}

\nolimits^{-2\gammaamma ^{-}t}\right) \right] ,  \lambdaabel{Eq70}
\end{equation}
leading to the decoherence time $\lambdaeft( 8\lambdaeft| \alpha \right| ^{2}{\gammaamma }
^{-}\right) ^{-1}$, which is ${\cal T}_{D}/2$ in the weak coupling regime,
as expected. However, in the strong coupling regime, this decoherence time
for an entangled state is equal to ${\cal T}_{D}$ for Markovian white noise
and is significantly improved for both Lorentzian spectral densities, where
we obtain ${\cal T}_{D}/2\varepsilon ^{-}$.

\subsection{Decoherence time of the state{\it \ }$\lambdaeft| ^\primeotect\varphi 
_{12}^{+}\right\rangle ={\cal N}_{\pm }\lambdaeft( \lambdaeft| ^\primeotect\alpha ,^\primeotect
\alpha \right\rangle _{12}\pm \lambdaeft| -^\primeotect\alpha ,-^\primeotect\alpha 
\right\rangle _{12}\right) $}

The decoherence process of the joint state $\lambdaeft| \varphi
_{12}^{+}\right\rangle $, in the strong coupling regime, is given by the
exponential decay
\begin{equation}gin{equation}
\exp \lambdaeft[ -4\lambdaeft| \alpha \right| ^{2}\lambdaeft( 1-

\mathop{\rm e}

\nolimits^{-\lambdaeft( \gammaamma _{1}^{+}+\gammaamma _{2}^{+}\right) t}\right) \right] ,
\lambdaabel{Eq71}
\end{equation}
depending only on the decay rate of the system-reservoir coupling around the
effective frequency $\omegamega _{\ell }^{+}$. Rewriting Eq. (\ref{Eq71}) as
\begin{equation}gin{equation}
\exp \lambdaeft\{ -4\lambdaeft| \alpha \right| ^{2}\lambdaeft( 1-\exp \lambdaeft[ -

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}+\gammaamma _{\ell }^{-}\right) t/2
\right] \exp \lambdaeft[ -

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) t/2
\right] \right) \right\} ,  \lambdaabel{Eq71l}
\end{equation}
we identify the term $\exp \lambdaeft[ -

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) t/2
\right] $ as coming from the cross-decay channel. It is evident that \ for $

\mathop{\thetaextstyle\sum}

\nolimits_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) >0$,
its contribution speeds up the exponential decay, in contrast to the above
situation where the eigenstate of the normal mode $\omegamega _{\ell }^{-}$ is
considered. For $\gammaamma _{\ell }^{-}={\gammaamma }^{-}$, Eq. (\ref{Eq71})
simplifies to
\begin{equation}gin{equation}
\exp \lambdaeft[ -4\lambdaeft| \alpha \right| ^{2}\lambdaeft( 1-

\mathop{\rm e}

\nolimits^{-2\gammaamma ^{+}t}\right) \right] ,  \lambdaabel{Eq72}
\end{equation}
and the decoherence time $\lambdaeft( 8\lambdaeft| \alpha \right| ^{2}{\gammaamma }
^{+}\right) ^{-1}$ becomes ${\cal T}_{D}/2$ for the weak coupling regime and
also for both spectral densities, Markovian white noise and the wide
Lorentzian spectral noise, in the strong coupling regime. For the Lorentzian
spectral density, we obtain ${\cal T}_{D}/2\varepsilon ^{+}$. In the light
of the above results we conclude, apart from the dependence of the
decoherence time upon the reservoir spectral densities, that the eigenstate
associated with the normal mode $\omegamega _{\ell }^{-}$ is less susceptible to
decoherence than that associated with $\omegamega _{\ell }^{+}$.

\section{Improving the quality factor of a dissipative system}

Let us now turn to a central result, which is the effect of the coupling $
\lambdaambda $ on the decoherence process for states of oscillators with
different damping rates $\Gammaamma _{1}$ and $\Gammaamma _{2}$. We assume that the
state of the whole system is given by Eq. (\ref{Eq46}) and that oscillator $
2 $ has a better quality factor than oscillator $1$ (i.e., $\Gammaamma
_{1}=\Gammaamma \gammag \Gammaamma _{2}$), where we have the superposition ${\cal N}
_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{1}$ which we
want to protect. The decoherence of state $1$, following from Eq. (\ref
{Eq65a}), is then governed by the exponential decay
\begin{equation}gin{equation}
\exp \lambdaeft\{ -2|\alpha |^{2}\lambdaeft[ 1-\frac{1}{4}\lambdaeft(

\mathop{\rm e}

\nolimits^{-\gammaamma ^{-}t}+

\mathop{\rm e}

\nolimits^{-\gammaamma ^{+}t}+2\cos (2\lambdaambda t)

\mathop{\rm e}

\nolimits^{-\lambdaeft( \gammaamma ^{+}+\gammaamma ^{-}\right) t/2}\right) \right]
\right\}  \lambdaabel{Eq74a}
\end{equation}
while for the field state $2$, following from (\ref{Eq65b}), it is governed
by
\begin{equation}gin{equation}
\exp \lambdaeft\{ -2|\alpha |^{2}\lambdaeft[ 1-\frac{1}{4}\lambdaeft(

\mathop{\rm e}

\nolimits^{-\gammaamma ^{-}t}+

\mathop{\rm e}

\nolimits^{-\gammaamma ^{+}t}+2\cos (2\lambdaambda t+\pi )

\mathop{\rm e}

\nolimits^{-\lambdaeft( \gammaamma ^{+}+\gammaamma ^{-}\right) t/2}\right) \right]
\right\} .  \lambdaabel{Eq74b}
\end{equation}
These expressions are to be compared with Eqs. (\ref{Eq67a}) and (\ref{Eq67b}
), respectively, where we have assumed the same quality factors for both
oscillators. The decoherence times obtained from Eqs. (\ref{Eq74a}) and (\ref
{Eq74b}) have an upper limit given by the exponential decay for the joint
system
\begin{equation}gin{equation}
\exp \lambdaeft[ -|\alpha |^{2}\lambdaeft( 2-

\mathop{\rm e}

\nolimits^{-\gammaamma ^{+}t}-

\mathop{\rm e}

\nolimits^{-\gammaamma ^{-}t}\right) \right] ,  \lambdaabel{Eq75}
\end{equation}
which is derived from Eq. (\ref{Eq62}). This exponential decay leads to a
decoherence time given by
\begin{equation}gin{equation}
\thetaau _{D}\approx \lambdaeft[ |\alpha |^{2}\lambdaeft( \gammaamma ^{+}+\gammaamma ^{-}\right) 
\right] ^{-1},  \lambdaabel{Eq76}
\end{equation}
which is twice the result found in Eq. (\ref{Eq63l}), for the case of
identical dissipative systems ($\gammaamma _{\ell }^{\pm }=\gammaamma ^{\pm }$),
whatever the spectral density of the reservoir. Therefore, even in the weak
coupling regime, the decoherence time of a system is improved when it is
coupled to another system with a better quality factor. We stress that this
conclusion holds only when $\lambdaambda \gammag \gammaamma _{\ell }^{\pm }$, a situation
we have assumed even for the weak coupling regime. Evidently, when $\lambdaambda
\lambdaessapprox \gammaamma _{\ell }^{\pm }$, so that the recurrence-swap dynamics
does not take place effectively before the relaxation time, the quality
factor of a system cannot be improved by coupling it to another system of
better quality.

In Fig. 6, assuming the weak coupling regime, $\gammaamma _{\ell }^{\pm }=\Gammaamma
_{\ell }/2$, with $\Gammaamma _{1}/\Gammaamma _{2}=10^{2}$ and $\lambdaambda /\Gammaamma
_{1}=5$ (a ratio chosen to show the dissipative dynamics clearly), the
dashed and dotted lines refer to the decoherence dynamics of systems $1$ and
$2$, derived from Eqs. (\ref{Eq74a}) and (\ref{Eq74b}), respectively. The
solid line represents the coherence decay for the joint system derived from
Eq. (\ref{Eq75}) and the dashed-dotted line indicates the coherence decay
for the isolated mode $1$, where the superpositon state ${\cal N}_{\pm
}\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{1}$ is prepared,
computed from Eq. (\ref{Eq64}). However, while the bad-quality system $1$
gets better when coupled to a good-quality system $2$, the latter, in turn,
gets worse. In fact, the thick solid line in Fig. 6, representing the
decoherence process of the superposition ${\cal N}_{\pm }\lambdaeft( |\alpha
\rangle \pm |-\alpha \rangle \right) _{2}$ prepared in an isolated
good-quality system $2$, displays a slower decay rate than the solid line
when system $2$ is coupled to a bad-quality system $1$. In conclusion, in a
network made of two oscillators with different quality factors, the
good-quality oscillator gets worse, while the bad-quality oscillator gets
better. We hope that this process may be even improved by extending the
coupling of a bad-quality oscillator to a higher number of good-quality
oscillators.

Finally, we note that the above result can be theoretically employed to
protect quantum superposition states generated, via atom-field interaction,
in open cavities. By coupling an open bad-quality cavity, where we have
prepared a quantum superposition, to a closed good-quality cavity, we can
protect the superposition through the coupling between the cavities and,
further, in an appropriate recurrence time, we can recover that
superposition state in the open cavity. The process works as if the original
superposition coherence is protected from the dissipative mechanism by the
recurrence-swap dynamics, much like a system which is put in contact with a
high-temperature reservoir, but is intermittently brought into contact with
a low-temperature reservoir. If the interval between contacts is fast
enough, the system will take longer to relax to the hot reservoir
temperature than it would in the absence of the cold reservoir. Of course,
this analogy has to be considered cautiously, since temperature and
coherence are very different features.

In the next section, we discuss how the original superposition state can be
recovered in the bad-quality system after its coupling with the good-quality
system.

\section{Entropy Excess}

To clarify the role the state-swap processes play in the coherence dynamics,
we consider the superposition state in Eq. (\ref{Eq46}), $|\psi _{12}\rangle
={}{\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right)
_{1}\omegatimes |\eta \rangle _{2}$, and plot, in Fig. 7(a and b), the
probabilities
\begin{equation}gin{equation}
{\cal P}_{\ell }(t)=\lambdaeft\lambdaangle \psi _{m}\lambdaeft| \rho _{\ell }(t)\right|
\psi _{m}\right\rangle ,  \lambdaabel{Eq77}
\end{equation}
of finding the superposition state $|\psi _{1}\rangle ={}{\cal N}_{\pm
}\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{1}$ and the coherent
state $|\psi _{2}\rangle =|\eta \rangle _{2}$ in oscillator $\ell $. In
these figures we assume the strong coupling regime ($\lambdaambda /\omegamega _{10}=2$
), identical dissipative systems ($\gammaamma _{\ell }^{\pm }=\gammaamma ^{\pm }$),
and Markovian white noise ($\gammaamma ^{+}=2\gammaamma ^{-}=\Gammaamma /2$). The ratio $
\lambdaambda /\Gammaamma =20$ is set to a fictitious scale to make clear the
state-swap and state-recurrence dynamics, as the strong oscillations
obtained with a realistic $\lambdaambda $ would be blurred by the slow
(dissipative) dynamics.

>From the solid line in Fig. 7(a), we observe that the superposition $|\psi
_{1}\rangle $ recurs to oscillator $1$, despite the dissipative process.
However, the coherent state $|\eta \rangle _{2}$ does not swap to oscillator
$1,$ as indicated by the dotted line. In fact, the swapping of \ the
coherent state $|\eta \rangle _{2}$ to oscillator $1$ would be indicated by
the occurrence of maxima of the dotted line between those of the solid line.
In Fig. (b) we observe that the coherent state $|\eta \rangle _{2}$ recurs
to oscillator $2$, as indicated by the solid line. Moreover, the state $
|\psi _{1}\rangle $ also swaps to oscillator $2$ (dotted line), indicating
that the superposition ${\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha
\rangle \right) $ is completely interchanged between the systems,
differently from the coherent state which does not swap to oscillator $1$.

Therefore, from Fig. 7(a) we conclude that the superposition state ${\cal N}
_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) $ prepared in a
bad-quality system and protected from decoherence by coupling this system to
a good-quality one, can easily be recovered in system $1$ by switching off
the coupling at the recurrence time $t_{R}=n\pi /\lambdaambda ,{\rm {\ }}
n=0,1,2,....$ After this time, the superposition state will be in oscillator
$1$ with a fidelity less than unity due to the dissipative process.

In Figs. 8(a and b) we plot the linear entropy for the joint state in Eq. (
\ref{Eq46}) (${\cal S}_{12}=1-$ {\rm Tr}$\rho _{12}^{2}$) and the reduced
states of oscillators $1$ and $2$ (${\cal S}_{\ell }=1-$ {\rm Tr}$\rho
_{\ell }^{2}$). In these figures we employ the same parameters considered in
Fig. 7, except that in Fig. 8(b) we use the ratio $\lambdaambda /\Gammaamma =2$,
instead of $\lambdaambda /\Gammaamma =20$, to show clearly the dissipative dynamics.
In Fig. 8(a) we analyze the recurrence-swap dynamics until around the
correlation time $\thetaau _{C}$, defined as the time when the entropy ${\cal S}
_{1}$(${\cal S}_{2}$) goes to about $0.1$ at the swap (recurrence) time. In
fact, as shown in Fig. 8(a), the minima of the entropy ${\cal S}_{\ell }$
move away from zero due to the development of an inevitable correlation
between the oscillators (due to the cross-decay channel) which thus become
permanently entangled. This correlation time, estimated as the time when the
minima of ${\cal S}_{\ell }$ approach $0.1$, is given by
\begin{equation}gin{equation}
\thetaau _{C}\approx \frac{1}{5\lambdaeft| \alpha \right| \sum_{\ell }\lambdaeft( \gammaamma
_{\ell }^{+}-\gammaamma _{\ell }^{-}\right) }.  \lambdaabel{Eq78}
\end{equation}
>From Eq. (\ref{Eq78}) we conclude that for the weak coupling regime, where $
\gammaamma _{\ell }^{+}=\gammaamma _{\ell }^{-}$, the correlation time goes to
infinity, i.e., the entropy ${\cal S}_{1}$(${\cal S}_{2}$) always returns to
zero in the swap (recurrence) time. Therefore, in the weak coupling regime
the oscillators do not get permanently entangled (due to the absence of the
cross-decay channel) and, since ${\cal S}_{1}$(${\cal S}_{2}$) always
returns to zero in the swap (recurrence) time, one is always able to\
recover a superposition state of a bad-quality oscillator coupled to a
good-quality one. However, even in the strong coupling regime the
correlation developed between the oscillators could not affect the process
of recovering a superposition state of a bad-quality oscillator coupled to a
good-quality one. Since the ratio $\thetaau _{C}/\thetaau _{D}$ obeys
\begin{equation}gin{equation}
\frac{\thetaau _{C}}{\thetaau _{D}}\approx \frac{\lambdaeft| \alpha \right| }{5}\frac{
\sum_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}+\gammaamma _{\ell }^{-}\right) }{
\sum_{\ell }\lambdaeft( \gammaamma _{\ell }^{+}-\gammaamma _{\ell }^{-}\right) },
\lambdaabel{Eq79}
\end{equation}
for $\thetaau _{C}/\thetaau _{D}\gammatrsim 1$ one can always recover the superposition
state ${\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) $,
with a considerable fidelity, in spite of the process of entanglement
between the oscillators.

The overall picture coming from Fig. 8(a) is that the recurrence dynamics
due to the strong coupling between the oscillators tends to restore the
coherence of the initial states while the dissipative dynamics promotes the
decoherence process. In Fig. 8(b) the linear entropy of the joint state $
{\cal S}_{12}$, represented by the thick solid line, starts from zero, goes
to a maximum due to the decoherence process and then returns to zero, since
in the asymptotic limit both oscillators reach a pure state: the vacuum or
some coherent state whose excitation depends on the amplification parameter $
F$. Meanwhile, as shown in Fig. 8(a), the linear entropy ${\cal S}_{\ell }$
of the reduced state of oscillator $\ell $, oscillates between zero and $0.5$
. The linear entropy ${\cal S}_{1}$(${\cal S}_{2}$), solid (dashed) line,
becomes zero when oscillator $1$($2$)\ assumes the state $|-\eta \rangle
_{1} $($|\eta \rangle _{2}$), as can be computed from Eqs. (\ref{Eq43a}) and
(\ref{Eq43b}).{\LambdaARGE \ }At the same time, ${\cal S}_{2}$(${\cal S}_{1}$)
bump into the thick solid line representing the linear entropy for the joint
state of the system ${\cal S}_{12}$, from above, indicating that the
superposition $|\psi _{1}\rangle $ has swapped (recurred) to oscillator $2$($
1$) on its way to decoherence. The maximal correlation between fields occurs
at the points where the recurrence and swap curves cross, as illustrated by
the dotted line representing the excess entropy, defined as
\begin{equation}gin{equation}
{\cal I}\equiv {\cal S}_{1}+{\cal S}_{2}-{\cal S}_{12}.  \lambdaabel{Eq80}
\end{equation}
We also observe from the dotted line that the minima of the excess entropy $
{\cal I}$ move away from zero, due to the development of an inevitable
correlation between the oscillators. Returning to Fig. 8(b), after reaching
its maximum the correlation ${\cal S}_{12}$ decays exponentially as a result
of the dissipation and the driving field, attaining zero in the asymptotic
limit when there is no correlation between the fields described by
stationary factorized coherent states. During the decay of correlation $
{\cal S}_{12}$, the linear entropy ${\cal S}_{\ell }$ of the reduced state
of oscillator $\ell $ does not attain the value ${\cal S}_{12}$, due to the
correlation described by the excess entropy ${\cal I}$.

\section{Conclusion}

We have presented a comprehensive treatment of the coherence dynamics in a
network composed of two coupled dissipative oscillators. First, we have
derived a master equation for both regimes of weak and strong coupling
between the oscillators. In the weak coupling regime the dissipative
mechanism of the individual oscillators is not significantly affected by
their interaction, which appears only in the von Neumann term of the master
equation. However, in the strong coupling regime the time evolution of the
density operator of the joint system is modified by a cross-decay channel
represented by a Liouville operator ${\cal L}_{12}\rho _{12}$ accounting for
the coupling between the oscillators. The appearance of this cross-decay
channel leads to interesting properties of the coherence dynamics of
strongly coupled oscillators.

After the mathematical development, we first analyzed the state-swap and the
state-recurrence dynamics, i.e., the probability that each oscillator
returns to its initial state and the probability of state-swapping between
the oscillators, respectively. In particular, we have analyzed these
processes for the case where both oscillators, prepared in the joint state $
{\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right)
_{1}\omegatimes |\eta \rangle _{2}$, present the same quality factor $\Gammaamma $.
In the weak coupling regime the recurrence and swap processes are fully
accomplished apart from the relaxation of the field states due to
dissipation. On the way to the strong coupling regime the recurrence process
remains unchanged while the swap dynamics is gradually lost as the coupling
gets stronger. In the strong coupling regime the swap dynamics is lost due
to phase mismatching between the coupling parameter $\lambdaambda $ and the
field-shifted frequencies $\omegamega _{\ell }$, as explained in Sec. V.
Evidently, in the strong coupling regime the field states recur much more
often than in the weak coupling regime before the relaxation takes place.

Next, aware of the decoherence dynamics in the weak coupling regime,
governed by the usual master equation (\ref{Eq19}) where a Liouville
operator ${\cal L}_{\ell }\rho _{12}$ accounts for the effect of the
reservoir on oscillator $\ell $, we turn to the decoherence process in the
strong coupling regime. In this regime, the normal-mode frequencies ($\omegamega
_{\ell }^{\pm }=\omegamega _{\ell }\pm \lambdaambda $) of the coupled systems are
strongly shifted, away from the oscillator frequencies $\omegamega _{\ell 0}$,
to regions of the reservoir frequency spaces where the spectral densities
may be significantly different from that around $\omegamega _{\ell 0}$. As the
spectral densities of the reservoirs play an important role in this regime,
three different spectral functions were considered for our analysis of the
decoherence process. When the normal-mode frequency $\omegamega _{\ell }^{-}$ is
shifted to regions around the origin of the frequency space, even for
Markovian white noise the coupling parameter of oscillator $\ell $ with its
reservoir, around $\omegamega _{\ell }^{-}$, becomes half its value around $
\omegamega _{\ell }^{+}$ (i.e., $\gammaamma _{\ell }^{-}=\Gammaamma /4$). When
considering, instead, a Lorentzian spectral function, which approaches zero
around the origin of the reservoir frequency space, the coupling constant $
\gammaamma _{\ell }^{-}$ becomes even smaller, resulting in decoherence times
for the joint and the reduced system states significantly longer than the
result computed for a single system plus reservoir. In Sec. VI the
decoherence times for three different states of the composed system are
computed, considering the three distinct reservoir spectral functions. Apart
from the dependence of the decoherence time upon the reservoir spectral
density, we stress that the coherence decay of the eigenstate of the normal
mode $\omegamega _{\ell }^{-}$ becomes slower than that for the eigenstate of
the normal mode $\omegamega _{\ell }^{+}$.

As discussed in Sec. III, apart from the possibility of considering
particular physical systems with appropriate reservoir spectral densities,
it is possible that specific spectral functions could be achieved through
engineered reservoirs. The difficult task of engineering strong interactions
between the oscillators, together with the achievement of specific reservoir
spectral functions, are the most sensitive problems in the way of the
physical implementation of the network here proposed. However, this proposal
might provide a motivation for future theoretical and experimental
investigations. In particular, the proposal could be applied to test the
Markovian white noise approximation or to probe the reservoir spectral
functions.

When considering the two systems to have different damping rates $\Gammaamma
_{1} $ and $\Gammaamma _{2}$ (with $\lambdaambda \gammag \Gammaamma _{1}$,$\Gammaamma _{2}$) we
demonstrated that the coupling between the systems, independently of its
strength, makes the good-quality system worse and the bad-quality system
better. This result can be employed to improve the quality factor of a
cavity, and thus to protect quantum superposition states generated, via
atom-field interaction, in open bad-quality cavities coupled to closed
good-quality ones. Evidently, this results holds for the weak coupling
regime because we have assumed that the coupling strength $\lambdaambda $ is
significantly larger than the system damping rates even for that regime. In
fact, when considering $\lambdaambda \gammag \Gammaamma _{\ell }$, the computed
improvement of the quality factor of a system follows from the
recurrence-swap dynamics which take place many times before the relaxation
time, protecting in a good-quality system a field state originally prepared
in a bad-quality one.

Finally, we have developed a careful analysis of the entropy excess in our
network. Supposing that the joint system is prepared in the factorized state
$|\psi _{12}\rangle ={}{\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha
\rangle \right) _{1}\omegatimes |\eta \rangle _{2}${,} we observed that in the
weak-coupling regime the recurrence-swap dynamics take place uninterruptedly
until the relaxation process rules it out, i.e., the entropy excess always
returns to zero in the recurrence and swap times. Therefore, in the
weak-coupling regime both oscillators always get disentangled in the
recurrence and swap times. Differently, in the strong coupling regime a
correlation is developed between the states ${\cal N}_{\pm }\lambdaeft( |\alpha
\rangle \pm |-\alpha \rangle \right) $ and $|\eta \rangle $, in a such way
that they get permanently entangled after a time interval we have called the
correlation time. In this way, one cannot recuperate - with a fidelity equal
to unity - a superposition state ${\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm
|-\alpha \rangle \right) $ prepared in a bad-quality system and protected
from decoherence through a strong coupling of such system with a
good-quality one. However, we computed the amplitude $\lambdaeft| \alpha \right| $
where the correlation time becomes longer than the decoherence time,
allowing the recuperation of the superposition state ${\cal N}_{\pm }\lambdaeft(
|\alpha \rangle \pm |-\alpha \rangle \right) $ with a good fidelity.

One of the main restrictions argued against quantum computation is that even
if each individual logic unit were only slightly affected by the decoherence
process, the coupling of a large number of such logic cells, for an actual
implementation such as the factorization of large numbers \cite{nielsen,shor}
, would decrease the decoherence time in such a way that the whole computing
process would be seriously compromised \cite{Haroche1}. Here we have shown,
at least for two strongly coupled sites, that the decoherence process is not
dictated simply by the excitation of the state involved in the logic
operation and the damping decay rate of the logic cells. In a strongly
interacting quantum network the decoherence process is delayed by the
cross-decay channel arising from the strong coupling between the logic
cells. In fact it may occur that for many coupled sites, as for the two
coupled cavities presented in this paper, the decoherence time increases.
Therefore, the discussion about strongly interacting oscillators we have
presented here is central to understand the decoherence process in a quantum
network; specifically, our conclusion is that the decoherence time depends
not only on the excitation of the entanglement involved in the logic
operation but also on the coupling strength between the logic units.

$^{1}$maponte@df.ufscar.br

$^{2}$marcos@df.ufscar.br

$^{3}$miled@df.ufscar.br

{\bf Acknowledgments}

We wish to express thanks for the support from FAPESP (under contracts
\#99/11859-3, \#00/15084-5, and \#02/02633-6) and CNPq (Intituto do Mil\^{e}
nio de Informa\c{c}\~{a}o Qu\^{a}ntica), Brazilian agencies. We also thank
R. M. Serra, C. J. Villas-B\^{o}as, V. V. Dodonov, and A. F. R. de Toledo
Piza for helpful discussions.

\begin{equation}gin{references}
\bibitemitem{Neumann} J. von Neumann, {\it Mathematical Foundations of Quantum
Mechanics}, Translated by R. T. Beyer, Princeton University Press (19995).

\bibitemitem{Zurek} W. H. Zurek, Phys. Rev. D {\bf 24}, 1516 (1981); {\bf 26},
1862 (1982); Physics Today {\bf 44}(10), 36 (1991).

\bibitemitem{CL} A. O. Caldeira and A. J. Leggett, Annals of Physics {\bf 149},
374 (1983); A. O. Caldeira and A. J. Laggett, Physica A {\bf 121}, 587
(1983).

\bibitemitem{Haroche} M. Brune, E. Hagley, J. Dreyer, X. Maitre, A. Maali, C.
Wunderlich, J. M. Raimond, and S. Haroche, Phys. Rev. Lett. {\bf 77}, 4887
(1996).

\bibitemitem{Wineland} D. M. Meekhof, C. Monroe, B. E. King, W. M. Itano, and D.
J. Wineland, Phys. Rev. Lett. {\bf 76}, 1796 (1996).

\bibitemitem{Glauber} B. R. Mollow and R. J. Glauber, Phys. Rev. {\bf 160}, 1076
(1967); {\it ibid}. {\bf 160}, 1097 (1967).

\bibitemitem{WM} D. F. Walls and J. Milburn, {\it Quantum Optics}
(Springer-Verlag, Berlin, 1994); M. O. Scully, and M. S. Zubairy, {\it 
Quantum Optics} (Cambridge University Press, Cambridge,1997).

\bibitemitem{Milburn} S. Schneider and G. J. Milburn, Phys. Rev. A {\bf 57},
3748 (1998); Phys. Rev. A {\bf 59}, 3766 (1999).

\bibitemitem{Vogel} C. Di Fidio and W. Vogel, Phys. Rev. A {\bf 62, }031802(R)
(2000).

\bibitemitem{Serra} R. M. Serra, N. G. de Almeida, W. B. da Costa, and M. H. Y.
Moussa, Phys. Rev. A {\bf 64}, 033419 (2001).

\bibitemitem{Budini} A. A. Budini, R. L. de Matos, and N. Zagury, J. Opt. B:
Quantum Semiclass. Opt. {\bf 4}, S462 (2002).

\bibitemitem{Comm} J. I. Cirac, P. Zoller, H. J. Kimble, and H. Mabuchi, Phys.
Rev. Lett. {\bf 78}, 3221 (1997); T. Pellizzari, {\it ibid}. {\bf 79}, 5242
(1997).H. J. Briegel, W. D\"{u}r, J. I. Cirac, and P. Zoller, {\it ibid}.
{\bf 81}, 5932 (1998); S. J. van Enk, H. J. Kimble, J. I. Cirac, and P.
Zoller, Phys. Rev. A {\bf 59}, 2659 (1999).

\bibitemitem{Comp} P. Shor, in: {\it Proceedings of the 35th Annual Symposium on
the Theory of Computer Science, }ed. S. Goldwasser (IEEE Computer Society
Press, Los Alamitos, CA, 1994) , p. 124; {\it ibid. }quant-ph/9508027; J. I.
Cirac and P. Zoller, Phys. Rev. Lett. {\bf 74}, 4091 (1995); Q. A.
Turchette, C. J. Hood, W. Lange, H. Mabuchi, and H. J. Kimble, ibid. {\bf 75}
, 4710 (1995); I. L. Chuang, L. M. K. Vandersypen, X. Zhou, D. W. Leung, and
S. Lloyd, Nature, {\bf 393}, 143 (1998); B. E. Kane, {\it ibid.} {\bf 393},
143 (1998); L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, M.
H. Sherwood, I. L. Chuang, {\it ibid.} {\bf 414}, 883 (2001); L. M. K.
Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, M. H. Sherwood, and I. L.
Chuang, Nature {\bf 414}, 883 (2001).

\bibitemitem{nielsen} M. A. Nielsen and I. L. Chuang, {\it Quantum Computation
and Quantum Information}, (Cambridge University Press, UK, 2000).

\bibitemitem{cnot} M. C. de Oliveira, and W. J. Munro, Phys. Rev. A {\bf 61},
42309 (2000).

\bibitemitem{Haroche1} S. Haroche, Physics Today {\bf 51}(7), 36 (1998).

\bibitemitem{Vitali} D. Vitali, and P. Tombesi, Phys. Rev. A {\bf 59}, 4178
(1999).

\bibitemitem{Tombesi} D. Vitali, P. Tombesi, and G. J. Milburn, Phys. Rev. Lett.
{\bf 79}, 2442 (1997); {\it ibid.} J. Mod. Opt. {\bf 44}, 2033 (1997).

\bibitemitem{pump1} M. C. de Oliveira, M. H. Y. Moussa, and S. S. Mizrahi, Phys.
Rev. A {\bf 61}, 63809 (2000); J. Opt B: Quantum Semiclass. Opt. {\bf 3}, 57
(2001).

\bibitemitem{Zoller} J. F. Poyatos, J. I. Cirac, and P. Zoller, Phys. Rev. Lett.
{\bf 77}, 4728 (1996).

\bibitemitem{Wineland1} C. J. Myatt, B. E. King, Q. A. Turchette, C. A. Sackett,
D. Kielpinski, W. M. Itano, C. Monroe, and D. J. Wineland, Nature {\bf 403},
269 (2000).

\bibitemitem{Matos} A. R. R. Carvalho, P. Milman, R. L. deMatosFilho, and L.
Davidovich, Phys. Rev. Lett. {\bf 86}, 4988 (2001).

\bibitemitem{Raimond} J. M. Raimond, M. Brune and S. Haroche, Phys. Rev. Lett.
{\bf 79}, 1964 (1997).

\bibitemitem{Ions} D. J. Wineland, C. Monroe, W. M. Itano, D. Leibfried, B. E.
King, and D. M. Meekhof, J. Res. NIST {\bf 103}, 259 (1998).

\bibitemitem{Nelson} N. Studart and S. S. Sokolov, in {\it Two-Dimensional
Electron System}, ed. by E. Y. Andrei, Kluwer Acabemic Publishers,
Netherlands (1997).

\bibitemitem{Walls} D. F. Walls, Z. Physik {\bf 234}, 231 (1970).

\bibitemitem{Charmicael} H. J Carmichael and D. F. Walls, Phys. A {\bf 6}, 1552
(1973).

\bibitemitem{Nemes} S. G. Mokarzel, A. N. Salgueiro, and M. C. Nemes, Phys. Rev.
A {\bf 65}, 044101 (2002).

\bibitemitem{Zoubi} H. Zoubi, M. Orenstien, and A. Ron, Phys. Rev. A {\bf 62},
033801 (2000).

\bibitemitem{Net} M. A. de Ponte, M. C. de Oliveira, and M. H. Y. Moussa, in
progress.

\bibitemitem{MMC} M. H. Y. Moussa, S. S. Mizrahi, and A. O. Caldeira, Phys.
Lett. A {\bf 221}, 145 (1996).

\bibitemitem{meu} M.~C.~de Oliveira, S.~S.~Mizrahi, and V.~V.~Dodonov, J. Opt. B
{\bf 1}, 610 (1999).

\bibitemitem{FLO} G. W. Ford, J. T. Lewis, and R. F. O'Connell, Phys. Rev. A
{\bf 37}, 4419 (1988); G. W. Ford and R. F. O'Connell, Physica A {\bf 37},
377 (1997).

\bibitemitem{Amir} M. Rosenau da Costa, A. O. Caldeira, S. M. Dutra, and H.
Westfahl, Jr., Phys. Rev. A {\bf 61}, 022107 (2000).

\bibitemitem{Leggett} A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A.
Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. {\bf 59}, 1 (1987).

\bibitemitem{MLW} G. J. Milburn, A. S. Lane, and D. F. Walls, Phys. Rev. A 27,
2804 (1983).

\bibitemitem{Gardiner} C. W. Gardiner, P. Zoller, Quantum noise : a handbook of
Markovian and non-Markovian quantum stochastic methods with applications to
quantum optics (Springer, Berlin, 2000).

\bibitemitem{Bocchieri} P. Bocchieri and A. Loinger, Phys. Rev. {\bf 107}, 337
(1957).

\bibitemitem{shor} P. W. Shor, in {\it Proceedings of the 35th Annual Symposium
on the Foundations of Computer Science, Los Amitos, CA} (IEEE Computer
Society press, NY (1994).
\end{references}

 {\bf figures caption}

Fig. 1 Sketch of the coupled dissipative oscillators, with oscillator $2$
submitted to a classical driving field.

Fig. 2 Damping function $\gammaamma _{\ell }\lambdaeft( \chi \right) $ assuming a
Lorentzian coupling $V_{\ell }$ between oscillator $\ell $ and its
respective reservoir. In the weak coupling regime the function $\gammaamma
_{\ell }\lambdaeft( \chi \right) $ is centered around $\omegamega _{\ell }$ (dotted
line). As $\lambdaambda $ increases, the damping function splits into two
Lorentzian functions whose peak heights are half the original value $\Gammaamma
_{\ell }$ (solid line). On the way to the strong coupling regime the two
peaks can be clearly distinguished as shown by the dashed line.

Fig. 3 Spectral density of the reservoir $\sigmama _{\ell }\lambdaeft( \chi \right)
$ for (a) Markovian white noise, (b) a Lorentzian spectral density, and (c)
a wide Lorentzian spectral density. The system-reservoir couplings around $
\omegamega _{\ell }^{\pm }$ are represent by the shaded regions.

Fig. 4{\bf \ }State-swap probability $P_{S}(t)$ (dotted line) and recurrence
probability $P_{R}(t)$ (solid line) as a function of the scaled time $
\lambdaambda t$, for the factorized state\ $|\psi _{12}\rangle ={}{\cal N}_{\pm
}\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right) _{1}\omegatimes |\eta
\rangle _{2}$, taking $\alpha =\eta =1$ as real parameters. In (a)-(c)
relaxation and driving field are disregarded and we assume (a) weak
coupling, (b) an intermediate coupling, and (c) strong coupling regimes. In
(d) and (e) we set $F=0$ but dissipation is included, considering (d) weak
coupling and (e) strong coupling regimes. Finally, in (f) dissipation is
disregarded and the driving field turned on.

Fig. 5 Decoherence dynamics for the factorized state $|\psi _{12}\rangle ={}
{\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right)
_{1}\omegatimes |\eta \rangle _{2}$, assuming both oscillators with the same
damping factor. Coherence decays of oscillator $1$ (dashed line), oscillator
$2$ (dotted line), joint system (solid line) and an isolated oscillator with
damping rate $\Gammaamma _{1}$ (dashed-dotted line) are indicated.

Fig. 6 Decoherence dynamics for the factorized state $|\psi _{12}\rangle ={}
{\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha \rangle \right)
_{1}\omegatimes |\eta \rangle _{2}$, assuming that the damping rate of
oscillator $1$ is larger than that of oscillator $2$. The dashed (dotted)
line refers to the decoherence dynamics of oscillator $1$($2$), \ while the
solid (dashed-dotted) line represents the coherence decay of the joint
system (an isolated system with the damping rate of oscillator $1$),
respectively. The thick solid line represents the coherence decay of a
superposition state ${\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm |-\alpha
\rangle \right) _{2}$ prepared in an isolated oscillator $2$.

Fig. 7 Probability ${\cal P}_{\ell }(t)=\lambdaeft\lambdaangle \psi _{m}\lambdaeft| \rho
_{\ell }(t)\right| \psi _{m}\right\rangle $, of finding the superposition
state $|\psi _{1}\rangle ={}{\cal N}_{\pm }\lambdaeft( |\alpha \rangle \pm
|-\alpha \rangle \right) _{1}$ and the coherent state $|\psi _{2}\rangle
=|\eta \rangle _{2}$ in oscillator $\ell $. We consider (a) $\ell =1$ and
(b) $\ell =2$. In these figures the strong coupling regime is assumed
together with identical dissipative systems and Markovian white noise.

Fig. 8 Linear entropy for the joint state (${\cal S}_{12}=1-$ {\rm Tr}$\rho
_{12}$), represented by the thick solid line, and the reduced states of
oscillators $1$ and $2$ (${\cal S}_{\ell }=1-$ {\rm Tr}$\rho _{\ell }$),
represented by solid and dashed lines, respectively. The dotted line
indicates the excess entropy. In Fig. 8(a) the recurrence-swap dynamics is
plotted until around the correlation time $\thetaau _{C}$ when a permanent
correlation is developed between the oscillators. In Fig. 8(b) we are
concerned with the relaxation of the network.

\end{document}