[ust]organization=Department of Mathematics and Physics, College of Science,
University of Santo Tomas, addressline=Espana Boulevard, Sampaloc,
city=Manila, postcode=1008, state=, country=Philippines
[rcnas]organization=Research Center for the Natural and Applied Sciences,
University of Santo Tomas, addressline=Espana Boulevard, Sampaloc,
city=Manila, postcode=1008, state=, country=Philippines

# Driven - Dissipative Dynamics of Ultracold Atoms Trapped in an Array of
Harmonic Potentials

Roland Cristopher F. Caballar<EMAIL_ADDRESS>
###### Abstract

We investigate the dynamics of a gas of ultracold atoms that are trapped in an
array of harmonic potentials and that interacts with a Bose - Einstein
condensate (BEC) that acts as a reservoir of Bogoliubov excitations. The
ground and excited energy levels of these trapped ultracold atoms are coupled
to each other via detuned Raman lasers with corresponding Rabi frequencies.
Once excited via the Raman lasers, these trapped ultracold atoms then return
to their ground energy levels, but not necessarily to their original trap
locations, by emitting Bogoliubov excitations into the BEC. This combination
of driving via Raman lasers to excited energy levels and dissipation via
interaction with the BEC resulting in emission of Bogoliubov excitations into
it will result in the trapped ultracold atoms approaching a steady state,
whereby the expectation value of the number of trapped ultracold atoms in each
harmonic trap of the array will attain a constant value over time. One can
then use this system to prepare states that require a definite number of atoms
in a particular energy level, such as BECs and atom lasers used for atom
interferometry and for tests of the foundations of quantum mechanics.

##  1 Introduction

In recent years, there has been an increase in the study and application of
open quantum systems in quantum information and for quantum technologies [1,
2, 3]. Advances in optics and in atom trapping and cooling in particular have
made it possible to create quantum systems that use the environment with which
they interact as resources that drive their time evolution towards a
particular quantum state [4, 5, 6, 7, 8]. A particular mechanism that enables
this is by using what is known as driven - dissipative dynamics, wherein a
quantum system is driven from a lower energy level to a higher energy state,
then decays back to its original energy level by interacting with an
environment that acts as a reservoir of excitations that carry away energy
from the system as it returns to its initial state. This continuous driving
towards a higher energy level and dissipation of excitations to return to its
original energy state will then cause the quantum system, coupled to an
environment, to evolve towards a given steady - state, which is dependent on
the choice of system, environment and interaction between the two. Driven -
dissipative dynamics is widely used in dissipative quantum state preparation,
wherein an initial quantum state evolves towards a desired final state by an
appropriate choice of environment and interaction with this chosen
environment. Examples of such dissipative quantum state preparation schemes
are given in Refs. [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24], which make use of various systems and various environments to prepare
quantum states of interest in quantum information and quantum computing. Of
particular interest among dissipative quantum state preparation schemes are
those that make use of driven - dissipative dynamics to prepare of Bose -
Einstein condensates (BECs), such as those described in Refs. [25, 26] due to
their widespread use and application in quantum computation and fundamental
tests of quantum mechanics [27, 28]. One particular dissipative BEC
preparation scheme is detailed in Refs. [29, 30], wherein a universal quantum
dissipative process is formulated using an appropriate choice of jump
operators for a Markovian master equation describing the dissipative dynamics
of spin - 1/2 bosons on a spatial lattice, with the system evolving via this
master equation towards a BEC of hard-core spin - 1/2 bosons. Another
dissipative BEC preparation scheme was described in Ref. [31], wherein it was
experimentally shown that driven - dissipative dynamics can be used in a dye -
filled microcavity with light leaking out on both sides to create a photon
BEC, albeit one with a small population (roughly of order of magnitude 1) of
photons.

In most dissipative quantum state preparation schemes, the environment is
either a thermal bath that absorbs thermal excitations emitted by trapped
atoms, or an optical bath that absorbs photons emitted from leaky optical
cavities. However, in Ref. [32], a new type of dissipative quantum preparation
scheme for encoding qubits was formulated, one which made use of a system of
trapped ultracold atoms immersed in a superfluid bath, which would play the
role of a reservoir of Bogoliubov excitations emitted by the trapped ultracold
atoms as they are driven from their initial energy states to excited energy
states. This dissipative quantum state preparation scheme was later adapted in
Ref. [33] for a gas of ultracold atoms trapped in an optical lattice and
immersed in a bath of Bogoliubov excitations, with the trapped ultracold atoms
evolving over time towards a many-body pure state via the driven - dissipative
dynamics of this quantum state preparation scheme. This scheme was also
adapted and further refined in Refs. [34, 35], wherein a gas of ultracold
atoms was trapped in a double well potential immersed in a harmonic trap and
coupled to a background BEC acting as a reservoir of Bogoliubov excitations,
with the driven - dissipative dynamics causing the trapped ultracold atoms to
evolve over time towards a many - body phase or particle number squeezed state
(which was done in Ref. [34]) or a many - body spin steady state (as was done
in Ref. [35]).

Having shown that it is possible to formulate a dissipative quantum
preparation scheme which uses a background BEC as a reservoir of Bogoliubov
excitations which is coupled to a system of trapped ultracold atoms, we then
investigate in this paper whether it is possible to use a variation of this
scheme to prepare a many - body quantum state whose particle number
expectation value remains constant over time, i. e. a trapped ultracold atom
gas with a constant particle number expectation value. In doing so, the
resulting dissipative quantum state preparation scheme can then be used as a
basis for formulating a dissipative BEC preparation scheme, since one of the
requirements for a BEC is that almost all, if not all, of the bosons
comprising the condensate occupy the same energy level, thus providing an
alternative to the dissipative BEC preparation schemes mentioned in Refs. [25,
29, 30, 31, 26]. To do that, we first describe, in section 2 of this paper,
the trapped ultracold atom system and the background BEC that we will be using
for our scheme, together with the manner of interaction between system and
background BEC, expressed in terms of the system, background BEC and
interaction Hamiltonians, respectively. In particular, the trapping potential
to be used for our ultracold atom system is an array of harmonic potentials,
wherein three of these harmonic potentials will be loaded with atoms from our
ultracold atom gas system, and will occupy the ground state energy levels of
this trap. These atoms occupying the ground state energy levels in these
harmonic potentials will then be driven, using Raman lasers with corresponding
Rabi frequencies, to excited energy levels in the adjacent harmonic potential
traps in the array. Once we have done this, we will then, in section 3, derive
the master equation that will describe the driven - dissipative dynamics of
our trapped ultracold atom system while it interacts with the background BEC.
In deriving the master equation, we will make use of certain approximations,
such as the phonon approximation for the Bogoliubov excitations emitted by the
trapped ultracold atom gas into the background BEC, and the Born - Markov
approximation for the coupling between the system and the background BEC. Once
we have obtained the master equation, we present, in section 4, the numerical
results of our time evolution of various initial states of the trapped
ultracold atom gas of the system using this derived master equation. Our
results show that, for particular values of physical parameters describing
both the trapped ultracold atom system and the background BEC, such as the
coupling strength between the system and the background BEC, the expectation
value of the number of atoms initially located in each of the harmonic
potentials in the trap array, and the particular form of the initial state of
the trapped ultracold atom gas, the trapped ultracold atom gas will, via the
driven - dissipative dynamics resulting from its interaction with the
background BEC as well as the continued excitation of the atoms in the ground
state to the excited state, evolve towards a state where the expectation value
of the particle number in each of the harmonic potentials in the trap array
will remain constant over time.

##  2 Derivation of the Interaction Hamiltonian for the System

For our trapped ultracold atom system, we assume that the trapping array
consists of five adjacent harmonic oscillator potentials, whose schematic
diagram is shown in figure 1 below.

Figure 1:  Schematic diagram of the harmonic trap array and the driven -
dissipative dynamics that occur as a result of the coupling between the ground
state energy levels $\varepsilon_{g}$ (indicated by the blue lines) at traps
1, 2, and 3 (with centers located at $x=-x_{0}$, $x=+x_{0}$ and $x=+3x_{0}$,
respectively) and the excited state energy levels $\varepsilon_{e}$ (indicated
by the orange lines) at the harmonic traps with centers located at $x=0$ and
$x=2x_{0}$, and the interaction between the trapped ultracold atoms and the
background BEC. The coupling between energy levels is achieved by Raman lasers
with Rabi frequencies $\Omega_{1}$ (green lines) for the ground state energy
levels at traps 1 and 2, and $\Omega_{2}$ (pink lines) for the ground state
energy levels at traps 2 and 3. These Raman lasers have corresponding
detunings $\Delta_{1}$ (yellow line) and $\Delta_{2}$ (red line). On the other
hand, as a result of the interaction between the trapped ultracold atoms and
the background BEC, Bogoliubov excitations (indicated by dashed grey lines for
atoms in traps 1 and 2 and by dashed lavender lines for atoms in traps 2 and
3) are emitted by atoms excited from the energy level $\varepsilon_{g}$ to the
energy level $\varepsilon_{e}$ as they decay back to $\varepsilon_{g}$.

The corresponding Hamiltonian for the system can then be written as

$\hat{H}_{S}=\varepsilon_{g,1}\hat{a}^{\dagger}_{g,1}\hat{a}_{g,1}+\varepsilon_{g,2}\hat{a}^{\dagger}_{g,2}\hat{a}_{g,2}+\varepsilon_{g,3}\hat{a}^{\dagger}_{g,3}\hat{a}_{g,3}+\varepsilon_{e,1}\hat{a}^{\dagger}_{e,1}\hat{a}_{e,1}+\varepsilon_{e,2}\hat{a}^{\dagger}_{e,2}\hat{a}_{e,2}$
(1)

We note that the ground state energies in nodes 1, 2 and 3 are all equal to
each other, i. e.
$\varepsilon_{g,1}=\varepsilon_{g,2}=\varepsilon_{g,3}=\varepsilon_{g}$, and
that the same is true for the excited state energies in the nodes between
those corresponding to the ground state energies, i. e.
$\varepsilon_{e,1}=\varepsilon_{e,2}=\varepsilon_{e}$. As shown in the figure,
a pair of Rabi lasers with Raman frequencies $\Omega_{1}$ and detuning
$\Delta_{1}$ are used to couple the ground states with energies
$\varepsilon_{g,1}$ and $\varepsilon_{g,2}$ with the excited state with energy
$\varepsilon_{e,1}$, while a second pair of Rabi lasers with Raman frequencies
$\Omega_{2}$ and detuning $\Delta_{2}$ are used to couple the ground states
with energies $\varepsilon_{g,2}$ and $\varepsilon_{g,3}$ with the excited
state with energy $\varepsilon_{e,2}$. Hence, the trap sites corresponding to
the excited states with energies $\varepsilon_{e,1}$ and $\varepsilon_{e,2}$
are coupled to the same site, which corresponds to the ground state with
energy $\varepsilon_{g,2}$, in contrast to the other two sites corresponding
to the ground states with energies $\varepsilon_{g,1}$ and
$\varepsilon_{g,3}$, to which there is only one excited state coupled to each
of these sites.

On the other hand, for the background BEC with which the trapped ultracold
atom gas interacts, its Hamiltonian is given as

$\hat{H}_{B}=\sum_{k}E_{k}\hat{b}^{\dagger}_{k}\hat{b}_{k}$ (2)

One can consider the gas of ultracold atoms trapped in the array of harmonic
potentials and the background BEC with which it interacts as a mixture of two
ultracold atomic gases, for which the interaction Hamiltonian will have the
following form:

$\hat{H}_{SB}=\frac{2\pi a_{SB}}{\mu}\int
dx\;\hat{\psi}^{\dagger}_{S}\hat{\psi}_{S}\hat{\psi}^{\dagger}_{B}\hat{\psi}_{B}$
(3)

Here, $\hat{\psi}_{S}$ and $\hat{\psi}_{B}$ are the field operators for the
system of trapped ultracold atoms and the background BEC, respectively. For
the system field operator, it can be written as

$\hat{\psi}_{S}=\phi_{g,1}(x)\hat{a}_{g,1}+\phi_{g,2}(x)\hat{a}_{g,2}+\phi_{g,3}(x)\hat{a}_{g,3}+\psi_{e,1}(x)\hat{a}_{e,1}+\phi_{e,2}(x)\hat{a}_{e,2}$
(4)

where $\hat{a}_{n,j}$ and $\hat{a}^{\dagger}_{n,j}$ are the annihilation and
creation operators, respectively, for the system corresponding to the energy
level $\varepsilon_{n,j}$, which has corresponding energy eigenstates
$\phi_{n,j}(x)$, with $n=g,e$ denoting the ground state or the excited state
of the system, and $j=1,2,3$ denoting the location of the center of the energy
eigenstate. In particular, for the ground state eigenstates,

$\displaystyle\phi_{g,1}(x)=\left(\frac{m_{s}\omega_{g}}{\pi\hbar}\right)^{1/4}\exp\left(-\frac{m_{s}\omega_{g}}{2\hbar}(x+x_{0})^{2}\right)\underset{\sigma_{g}\rightarrow
0}{\longrightarrow}\sqrt{2}\pi^{1/4}\sigma_{g}^{1/2}\delta(x+x_{0})$
$\displaystyle\phi_{g,2}(x)=\left(\frac{m_{s}\omega_{g}}{\pi\hbar}\right)^{1/4}\exp\left(-\frac{m_{s}\omega_{g}}{2\hbar}(x-x_{0})^{2}\right)\underset{\sigma_{g}\rightarrow
0}{\longrightarrow}\sqrt{2}\pi^{1/4}\sigma_{g}^{1/2}\delta(x-x_{0})$
$\displaystyle\phi_{g,3}(x)=\left(\frac{m_{s}\omega_{g}}{\pi\hbar}\right)^{1/4}\exp\left(-\frac{m_{s}\omega_{g}}{2\hbar}(x-3x_{0})^{2}\right)\underset{\sigma_{g}\rightarrow
0}{\longrightarrow}\sqrt{2}\pi^{1/4}\sigma_{g}^{1/2}\delta(x-3x_{0})$ (5)

wherein we make use of the property of Dirac delta distributions that

$\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{\sigma^{2}}(x-x_{0})^{2}\right)\underset{\sigma\rightarrow
0}{\longrightarrow}\delta(x-x_{0})$

and where these eigenstates have width
$\sigma_{g}=\sqrt{\frac{\hbar}{m_{s}\omega_{g}}}$. From these expressions, we
see that the ground - state eigenstates $\phi_{g,1}(x)$, $\phi_{g,2}(x)$ and
$\phi_{g,3}(x)$ of the system are centered at $x=-x_{0}$, $x=x_{0}$ and
$x=3x_{0}$, respectively. On the other hand, for the excited state
eigenstates,

$\displaystyle\phi_{e,1}(x)=\left(\frac{m_{s}\omega_{e}}{\pi\hbar}\right)^{1/4}\sqrt{\frac{2m_{s}\omega_{e}}{\hbar}}x\exp\left(-\frac{m_{s}\omega_{e}}{2\hbar}x^{2}\right)=\sqrt{2}\pi^{-1/4}\sigma_{e}^{-3/4}x\exp\left(-\frac{x^{2}}{2\sigma_{e}^{2}}\right)$
$\displaystyle\phi_{e,2}(x)=\left(\frac{m_{s}\omega_{e}}{\pi\hbar}\right)^{1/4}\sqrt{\frac{2m_{s}\omega_{e}}{\hbar}}(x-2x_{0})\exp\left(-\frac{m_{s}\omega_{e}}{2\hbar}(x-2x_{0})^{2}\right)$
$\displaystyle=\sqrt{2}\pi^{-1/4}\sigma_{e}^{-3/4}(x-2x_{0})\exp\left(-\frac{(x-2x_{0})^{2}}{2\sigma_{e}^{2}}\right)$
(6)

These eigenstates have width $\sigma_{e}=\sqrt{\frac{\hbar}{m_{s}\omega_{e}}}$
and are centered at $x=0$ and $x=2x_{0}$, respectively. Multiplying the field
operator of the system with its adjoint, and considering only terms up to
first order in the system creation and annihilation operators
$\hat{a}^{\dagger}_{n,j}$ and $\hat{a}_{n,j}$, we obtain

$\displaystyle\hat{\psi}^{\dagger}_{S}\hat{\psi}_{S}=\phi_{e,1}(x)\phi_{g,1}(x)(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})+\phi_{e,1}(x)\phi_{g,2}(x)(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})$
$\displaystyle+\phi_{e,2}(x)\phi_{g,2}(x)(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})+\phi_{e,2}(x)\phi_{g,3}(x)(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})$
(7)

We note that in computing for this product, we neglected the particle number
terms $\hat{a}^{\dagger}_{n,j}\hat{a}_{n,j}$ as well as terms proportional to
$\hat{a}^{\dagger}_{g,3}\hat{a}_{e,1}$ and its conjugate. The former is due to
the assumption that particle number is conserved (and hence remains constant
over time, so these terms have no contribution to the time evolution of the
system), while the latter is due to the nearest - neighbor approximation,
whereby atoms trapped in one site can only be coupled to energy levels in
adjoining sites. This approximation allows us to avoid any complications that
may arise due to the possibility of atoms tunneling between non-adjacent trap
array sites.

On the other hand, the field operator for the background BEC has the form

$\hat{\psi}_{B}(x)=\sqrt{\rho_{B}}+\delta\hat{\psi}_{B}(x)$ (8)

Here, $\sqrt{\rho_{B}}$ is the BEC density, while $\delta\hat{\psi}_{B}$ is
the condensate excitation term, and has the explicit form

$\delta\hat{\psi}_{B}=\frac{1}{\sqrt{L}}\sum_{k}\left(u_{k}e^{ikx}\hat{b}_{k}+v_{\mathbf{k}}e^{-ikx}\hat{b}^{\dagger}_{k}\right),$
(9)

In this term, $L$ is the length of the BEC, $u_{k}=(1-L^{2}_{k})^{-1/2}$,
$v_{k}=L_{k}(1-L^{2}_{k})^{-1/2}$,
$L_{k}=(E_{k}-(k^{2}/2m_{B})-m_{B}c^{2})/m_{B}c^{2}$, and $E_{k}$ is the
excitation energy given by

$E_{k}=ck\sqrt{1+\left(\frac{k}{2m_{B}c}\right)^{2}}$ (10)

with $c=\sqrt{g_{BB}\rho_{B}/m_{B}}$ being the velocity of sound in the
condensate and $g_{BB}$ the interaction strength among the condensate atoms.
Here, we make use of the assumption that the excitatons emitted into and from
the BEC are sound - like, so that $E_{k}\approx ck$. However, this
approximation implies that $k^{2}/4(m_{B}c)^{2}<<1$. As such, the terms in the
BEC field operator will then have the following form:

$\displaystyle
L_{k}=\frac{1}{m_{B}c^{2}}\left(E_{k}-\frac{k^{2}}{2m_{B}}-m_{B}c^{2}\right)\approx\frac{k}{m_{B}c}-1$
$\displaystyle
u_{k}+v_{k}=\frac{1+L_{k}}{\sqrt{1-L^{2}_{k}}}=\sqrt{\frac{1+L_{k}}{1-L_{k}}}\approx\sqrt{\frac{k/m_{B}c}{2-k/m_{B}c}}$
$\displaystyle=\sqrt{\frac{k}{2m_{B}c}}\left(1-\frac{k}{2m_{B}c}\right)^{-1/2}\approx\sqrt{\frac{k}{2m_{B}c}}$

The second of the two equations given above will be crucial in the derivation
of the explicit form of $\hat{\psi}^{\dagger}_{B}\hat{\psi}_{B}$, which, if we
are to keep only terms up to first order in $\hat{b}_{k}$ and
$\hat{b}^{\dagger}_{k}$, can be written as follows:

$\displaystyle\hat{\psi}^{\dagger}_{B}\hat{\psi}_{B}=\left(\sqrt{\rho}_{B}+\frac{1}{\sqrt{L}}\sum_{k}(u_{k}e^{-ikx}\hat{b}^{\dagger}_{k}+v_{k}e^{ikx}\hat{b}_{k})\right)\left(\sqrt{\rho}_{B}+\frac{1}{\sqrt{L}}\sum_{k}(u_{k}e^{ikx}\hat{b}_{k}+v_{k}e^{-ikx}\hat{b}^{\dagger}_{k})\right)$
$\displaystyle\approx\sqrt{\frac{\rho_{B}}{L}}\sum_{k}(u_{k}+v_{k})(e^{ikx}\hat{b}_{k}+e^{-ikx}\hat{b}^{\dagger}_{k})=\sqrt{\frac{\rho_{B}}{2m_{B}cV}}\sum_{k}\sqrt{k}(e^{ikx}\hat{b}_{k}+e^{-ikx}\hat{b}^{\dagger}_{k})$
(12)

Thus, substituting Eqs. 7 and 12 into Eq. 3, we obtain the following
expression:

$\displaystyle\hat{H}_{SB}=\frac{2\pi a_{SB}}{\mu}\int
dx\;\hat{\psi}^{\dagger}_{S}\hat{\psi}_{S}\hat{\psi}^{\dagger}_{B}\hat{\psi}_{B}$
$\displaystyle=\frac{2\pi
a_{SB}}{\mu}\sqrt{\frac{\rho_{B}}{2m_{B}cV}}\sum_{k}\sqrt{k}\int
dx\;\left(\phi_{e,1}(x)\phi_{g,1}(x)(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})+\phi_{e,1}(x)\phi_{g,2}(x)\right.$
$\displaystyle\left.\times(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})+\phi_{e,2}(x)\phi_{g,2}(x)(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})+\phi_{e,2}(x)\phi_{g,3}(x)(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)$
$\displaystyle\times(e^{ikx}\hat{b}_{k}+e^{-ikx}\hat{b}^{\dagger}_{k})$
$\displaystyle=\frac{2\pi
a_{SB}}{\mu}\sqrt{\frac{\rho_{B}}{2m_{B}cV}}\sum_{k}\sqrt{k}\left((\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})\left(\hat{b}_{k}\int
dx\;e^{ikx}\phi_{e,1}(x)\phi_{g,1}(x)\right.\right.$
$\displaystyle\left.+\hat{b}^{\dagger}_{k}\int
dx\;e^{-ikx}\phi_{e,1}(x)\phi_{g,1}(x)\right)+(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\left(\hat{b}_{k}\int
dx\;e^{ikx}\phi_{e,1}(x)\phi_{g,2}(x)\right.$
$\displaystyle\left.+\hat{b}^{\dagger}_{k}\int
dx\;e^{-ikx}\phi_{e,1}(x)\phi_{g,2}(x)\right)+(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})\left(\hat{b}_{k}\int
dx\;e^{ikx}\phi_{e,2}(x)\phi_{g,2}(x)\right.$
$\displaystyle\left.+\hat{b}^{\dagger}_{k}\int
dx\;e^{-ikx}\phi_{e,2}(x)\phi_{g,2}(x)\right)+(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\left(\hat{b}_{k}\int
dx\;e^{ikx}\phi_{e,2}(x)\phi_{g,3}(x)\right.$
$\displaystyle\left.\left.+\hat{b}^{\dagger}_{k}\int
dx\;e^{-ikx}\phi_{e,2}(x)\phi_{g,3}(x)\right)\right)$ (13)

In this expression for the interaction Hamiltonian, we see that there are a
number of overlap integrals of the form $\int dx\;e^{\pm
i\mathbf{k}\cdot\mathbf{r}}\phi_{n,j}(x)\phi_{n^{\prime},j^{\prime}}(x)$.
Substituting the expressions for the ground state and excited state
eigenstates given by Eqs. 5 and 6 into the overlap integrals, we find that

$\displaystyle\int dx\;e^{\pm
ikx}\phi_{e,1}(x)\phi_{g,1}(x)=2\sigma_{g}^{1/2}\sigma_{e}^{-3/4}\int
dx\;e^{\pm
ikx}x\exp\left(-\frac{x^{2}}{2\sigma_{e}^{2}}\right)\delta(x+x_{0})$
$\displaystyle=-\frac{2}{\sigma^{1/4}_{e}}\sqrt{\frac{\sigma_{g}}{\sigma_{e}}}e^{\mp
ikx_{0}}x_{0}\exp\left(-\frac{x_{0}^{2}}{2\sigma_{e}^{2}}\right)$
$\displaystyle\int dx\;e^{\pm
ikx}\phi_{e,1}(x)\phi_{g,2}(x)=2\sigma_{g}^{1/2}\sigma_{e}^{-3/4}\int
dx\;e^{\pm
ikx}x\exp\left(-\frac{x^{2}}{2\sigma_{e}^{2}}\right)\delta(x-x_{0})$
$\displaystyle=\frac{2}{\sigma^{1/4}_{e}}\sqrt{\frac{\sigma_{g}}{\sigma_{e}}}e^{\pm
ikx_{0}}x_{0}\exp\left(-\frac{x_{0}^{2}}{2\sigma_{e}^{2}}\right)$
$\displaystyle\int dx\;e^{\pm
ikx}\phi_{e,2}(x)\phi_{g,2}(x)=2\sigma_{g}^{1/2}\sigma_{e}^{-3/4}\int
dx\;e^{\pm
ikx}(x-2x_{0})\exp\left(-\frac{(x-2x_{0})^{2}}{2\sigma_{e}^{2}}\right)\delta(x-x_{0})$
$\displaystyle=-\frac{2}{\sigma^{1/4}_{e}}\sqrt{\frac{\sigma_{g}}{\sigma_{e}}}e^{\pm
2ikx_{0}}x_{0}\exp\left(-\frac{x_{0}^{2}}{2\sigma_{e}^{2}}\right)$
$\displaystyle\int dx\;e^{\pm
ikx}\phi_{e,2}(x)\phi_{g,3}(x)=2\sigma_{g}^{1/2}\sigma_{e}^{-3/4}\int
dx\;e^{\pm
ikx}(x-2x_{0})\exp\left(-\frac{(x-2x_{0})^{2}}{2\sigma_{e}^{2}}\right)\delta(x-3x_{0})$
$\displaystyle=\frac{2}{\sigma^{1/4}_{e}}\sqrt{\frac{\sigma_{g}}{\sigma_{e}}}e^{\pm
3ikx_{0}}x_{0}\exp\left(-\frac{x_{0}^{2}}{2\sigma_{e}^{2}}\right)$ (14)

Now before we use these evaluated overlap integrals to simplify the
interaction Hamiltonian between the trapped ultracold atom system and the
background BEC, we first use the approximation that $kx_{0}<<1$ for all $k$.
This is consistent with our assumption that the Bogoliubov excitations emitted
by the ultracold atoms as they return to the ground state energy in the trap
array are phononic, or have energies that vary linearly with $k$, i. e.
$E_{k}=ck$. For this assumption to hold, the momenta $\hbar k$ of the
excitations must be very small, implying that $k<<1$. At the same time, to
avoid tunneling between traps, we must make $x_{0}$ much larger than the de
Broglie wavelength of the ultracold atoms, so that $x_{0}>>1$. Therefore,
given these orders of magnitude for both $k$ and $x_{0}$, we then obtain the
approximation mentioned earlier, $kx_{0}<<1$. Under this approximation,
$e^{\pm ikx_{0}}\approx 1$, so that these factors which are present in the
overlap integrals as given in Eq. 14 can be neglected. Thus, following this
condition, we can then substitute the evaluated overlap integrals, sans the
factors $e^{\pm ikx_{0}}$, into Eq 13 to obtain the following form of the
interaction Hamiltonian:

$\displaystyle\hat{H}_{SB}=\frac{4\pi
a_{SB}}{\mu\sigma_{e}^{1/4}}\sqrt{\frac{\rho_{B}}{2m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}}x_{0}\exp\left(-\frac{x_{0}^{2}}{2\sigma_{e}^{2}}\right)\sum_{k}\sqrt{k}\left(-(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})\right.$
$\displaystyle+(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})-(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})$
$\displaystyle\left.+(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)(\hat{b}_{k}+\hat{b}^{\dagger}_{k})$
(15)

We then evolve this Hamiltonian over time, making use of the Baker - Campbell
- Hausdorff (BCH) identity as well as the commutation relations between the
system and the BEC creation and annihilation operators in doing so.
Explicitly, using the BCH identity, a time - evolved quantum operator
$\hat{A}$ in a Hilbert space is given as

$\hat{A}(t)=\exp\left(\frac{i}{\hbar}t\hat{H}\right)\hat{A}\exp\left(-\frac{i}{\hbar}t\hat{H}\right)=\hat{A}+\frac{i}{\hbar}t\left[\hat{H},\hat{A}\right]+\frac{1}{2!}\left(\frac{i}{\hbar}t\right)^{2}\left[\hat{H},\left[\hat{H},\hat{A}\right]\right]+...$
(16)

On the other hand, the creation and annihilation operators for the system and
the background BEC obey the following commutation relations:

$\left[\hat{a}_{n,j},\hat{a}^{\dagger}_{n^{\prime},j^{\prime}}\right]=\delta_{n,n^{\prime}}\delta_{j,j^{\prime}},\;\left[\hat{b}_{k},\hat{b}^{\dagger}_{k^{\prime}}\right]=\delta_{k,k^{\prime}}$
(17)

Thus, substituting the interaction Hamiltonian given by Eq. 15 into Eq. 16,
with the Hamiltonian term in the latter being the sum of the system and
background BEC Hamiltonians, $\hat{H}=\hat{H}_{S}+\hat{H}_{B}$, with
$\hat{H}_{S}$ and $\hat{H}_{B}$ given by Eqs. 1 and 2, respectively, and
making use of the commutation relations given by Eq. 17 to evaluate the
resulting commutators for the system and bath creation and annihilation
operators, the time - evolved interaction Hamiltonian will then have the
following form:

$\displaystyle\hat{H}_{SB}(t)=\frac{4\pi
a_{SB}}{\mu\sigma_{e}^{1/4}}\sqrt{\frac{\rho_{B}}{2m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}}x_{0}\exp\left(-\frac{x_{0}^{2}}{2\sigma_{e}^{2}}\right)\sum_{k}\sqrt{k}$
$\displaystyle\times\left(-(e^{\frac{it}{\hbar}(\varepsilon_{e,1}-\varepsilon_{g,1})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+e^{-\frac{it}{\hbar}(\varepsilon_{e,1}-\varepsilon_{g,1})}\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})\right.$
$\displaystyle+(e^{\frac{it}{\hbar}(\varepsilon_{e,1}-\varepsilon_{g,2})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+e^{-\frac{it}{\hbar}(\varepsilon_{e,1}-\varepsilon_{g,2})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})-(e^{\frac{it}{\hbar}(\varepsilon_{e,2}-\varepsilon_{g,2})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+e^{-\frac{it}{\hbar}(\varepsilon_{e,2}-\varepsilon_{g,2})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})$
$\displaystyle\left.+(e^{\frac{it}{\hbar}(\varepsilon_{e,2}-\varepsilon_{g,3})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+e^{-\frac{it}{\hbar}(\varepsilon_{e,2}-\varepsilon_{g,3})}\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)(e^{-\frac{it}{\hbar}E_{k}}\hat{b}_{k}+e^{\frac{it}{\hbar}E_{k}}\hat{b}^{\dagger}_{k})$
(18)

##  3 Derivation of the Master Equation for the Dissipative Dynamics of the
System

Having derived the interaction Hamiltonian describing the interaction between
the trapped ultracold atom system and the background BEC, we now proceed to
derive the master equation that describes the dissipative dynamics of the
system. The general form of the master equation under the Born - Markov
approximation is given as (see, for example, Ref. [1])

$\frac{d}{dt}\hat{\rho}_{S}=-\int_{0}^{\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left[\hat{H}_{SB}(t),\left[\hat{H}_{SB}(t-t^{\prime}),\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\right]\right]$
(19)

Expanding the commutation relations in the integrand, this master equation can
be written as

$\displaystyle\frac{d}{dt}\hat{\rho}_{S}=-\int_{0}^{\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\right)$
$\displaystyle+\int_{0}^{\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\right)+\int_{0}^{\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t)\right)$
$\displaystyle-\int_{0}^{\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\hat{H}_{SB}(t)\right)$
(20)

In evaluating the commutators and the trace with respect to the background BEC
variables in the integrand, we make use of the following identities, following
Ref. [32]:

$\mathrm{Tr}_{B}\left(\hat{b}_{k}\hat{b}^{\dagger}_{k^{\prime}}\hat{\rho}_{B}\right)=\delta_{k,k^{\prime}},\;\mathrm{Tr}_{B}\left(\hat{b}^{\dagger}_{k}\hat{b}_{k^{\prime}}\hat{\rho}_{B}\right)=0,\;\mathrm{Tr}_{B}\left(\hat{b}_{k}\hat{b}_{k^{\prime}}\hat{\rho}_{B}\right)=\mathrm{Tr}_{B}\left(\hat{b}^{\dagger}_{k}\hat{b}^{\dagger}_{k^{\prime}}\hat{\rho}_{B}\right)=0$
(21)

These identities are implemented to ensure that no excitations are emitted
from the BEC to the ultracold atom gas to couple the ground states to the
excited states (with excitation emission by the BEC corresponding to the
second identity), while at the same time ensuring that the ultracold atoms in
the excited energy state are able to decay back to one of the three available
ground states in the trap via emission of excitations into the BEC (with
excitation absorption by the BEC corresponding to the first identity). At the
same time, we can make use of the cyclic nature of the trace, i. e.
$\mathrm{Tr}_{B}\left(\hat{b}_{k}\hat{b}^{\dagger}_{k^{\prime}}\hat{\rho}_{B}\right)=\mathrm{Tr}_{B}\left(\hat{b}^{\dagger}_{k^{\prime}}\hat{\rho}_{B}\hat{b}_{k}\right)=\mathrm{Tr}_{B}\left(\hat{\rho}_{B}\hat{b}_{k}\hat{b}^{\dagger}_{k^{\prime}}\right)=\delta_{k,k^{\prime}}$.
Thus, evaluating the commutator and the trace over the background BEC
variables, and recalling that
$\varepsilon_{g,1}=\varepsilon_{g,2}=\varepsilon_{g,3}=\varepsilon_{g}$ and
$\varepsilon_{e,1}=\varepsilon_{e,2}=\varepsilon_{e}$ we obtain the following
terms in the integrand:

$\displaystyle\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\right)=\frac{8\pi^{2}a^{2}_{SB}}{\mu^{2}\sigma_{e}^{1/2}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}x^{2}_{0}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)\sum_{k}ke^{\frac{it^{\prime}}{\hbar}(\varepsilon_{e}-\varepsilon_{g}-ck)}$
$\displaystyle\times\left((\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)$
$\displaystyle\times\left((e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)\hat{\rho}_{S}(t)$
(22)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\right)=\frac{8\pi^{2}a^{2}_{SB}}{\mu^{2}\sigma_{e}^{1/2}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}x^{2}_{0}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)\sum_{k}ke^{\frac{it^{\prime}}{\hbar}(\varepsilon_{e}-\varepsilon_{g}+ck)}$
$\displaystyle\times\left((\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)\hat{\rho}_{S}(t)$
$\displaystyle\times\left((e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)$
(23)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t)\right)=\frac{8\pi^{2}a^{2}_{SB}}{\mu^{2}\sigma_{e}^{1/2}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}x^{2}_{0}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)\sum_{k}ke^{-\frac{it^{\prime}}{\hbar}(\epsilon_{e}-\epsilon_{g}+ck)}$
$\displaystyle\times\left((\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)\hat{\rho}_{S}(t)$
$\displaystyle\times\left((e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)$
(24)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\hat{H}_{SB}(t)\right)=\frac{8\pi^{2}a^{2}_{SB}}{\mu^{2}\sigma_{e}^{1/2}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}x^{2}_{0}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)\sum_{k}ke^{-\frac{it^{\prime}}{\hbar}(\varepsilon_{e}-\varepsilon_{g}-ck)}$
$\displaystyle\times\hat{\rho}_{S}(t)\left((\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)$
$\displaystyle\times\left((e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,1}+\hat{a}^{\dagger}_{g,1}\hat{a}_{e,1})-(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,1}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,1})\right.$
$\displaystyle\left.+(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,2}+\hat{a}^{\dagger}_{g,2}\hat{a}_{e,2})-(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{a}^{\dagger}_{e,2}\hat{a}_{g,3}+\hat{a}^{\dagger}_{g,3}\hat{a}_{e,2})\right)$
(25)

Having calculated the explicit form of the integrands in Eq. 20, we now
proceed to simplify these terms. To do so, we make use of adiabatic
elimination to express the creation and annihilation operators for the excited
states of the system as linear combinations of the creation and annihilation
operators of its ground states, as follows:

$\hat{a}_{e,1}=\frac{\Omega_{1}}{\sqrt{2}\Delta_{1}}(\hat{a}_{g,1}+\hat{a}_{g,2}),\;\hat{a}_{e,2}=\frac{\Omega_{2}}{\sqrt{2}\Delta_{2}}(\hat{a}_{g,2}+\hat{a}_{g,3})$
(26)

Substituting these operators into Eqs. 22 to 25 and simplifying the resulting
expressions, we then obtain

$\displaystyle\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\right)=\frac{4\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{\sigma_{e}}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\times\sum_{k}ke^{\frac{it^{\prime}}{\hbar}(\varepsilon_{e}-\varepsilon_{g}-ck)}\left(((\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon((\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
$\displaystyle\times\left((e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3})\right)\hat{\rho}_{S}(t)$
(27)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\right)=\frac{4\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{\sigma_{e}}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\times\sum_{k}ke^{\frac{it^{\prime}}{\hbar}(\varepsilon_{e}-\varepsilon_{g}+ck)}\left(((\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon((\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)\hat{\rho}_{S}(t)$
$\displaystyle\times\left((e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(e^{\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
(28)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t)\right)=\frac{4\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{\sigma_{e}}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\times\sum_{k}ke^{-\frac{it^{\prime}}{\hbar}(\epsilon_{e}-\epsilon_{g}+ck)}\left(((\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon((\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)\hat{\rho}_{S}(t)$
$\displaystyle\times\left((e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
(29)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\hat{H}_{SB}(t)\right)=\frac{4\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{\sigma_{e}}}\frac{\rho_{B}}{m_{B}cL}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\times\sum_{k}ke^{-\frac{it^{\prime}}{\hbar}(\varepsilon_{e}-\varepsilon_{g}-ck)}\hat{\rho}_{S}(t)\left(((\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon((\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+e^{-\frac{2i(t-t^{\prime})}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
$\displaystyle\times\left((e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
(30)

To simplify these terms further, we now proceed to integrating them over time.
At the same time, we replace the summation over $k$ with an integration over
$k$, treating the summation variable as continuous to be able to do so. In
doing so, we make use of the definition of the Dirac delta distribution in
terms of oscillatory integrals, giving us the following integration
identities:

$\displaystyle\sum_{k}\int_{0}^{+\infty}dt^{\prime}e^{\pm\frac{i}{\hbar}(\varepsilon_{e}-\varepsilon_{g}-ck)t^{\prime}}k^{\prime}=\frac{\hbar}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}k^{\prime}\delta(\varepsilon_{e}-\varepsilon_{g}-ck)=\frac{1}{\sqrt{2\pi}}\frac{\hbar}{c}(\varepsilon_{e}-\varepsilon_{g})$
$\displaystyle\sum_{k}\int_{0}^{+\infty}dt^{\prime}e^{\pm\frac{i}{\hbar}(\varepsilon_{e}-\varepsilon_{g}+ck)t^{\prime}}k^{\prime}=\frac{\hbar}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}k^{\prime}\delta(\varepsilon_{e}-\varepsilon_{g}+ck)=-\frac{1}{\sqrt{2\pi}}\frac{\hbar}{c}(\varepsilon_{e}-\varepsilon_{g})$
(31)

Thus, integrating the terms given by Eqs. 27 to 30 over time and over $k$, and
making use of the integration identities given by Eq. 31 in doing so, we find
that

$\displaystyle\int_{0}^{+\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\right)$
$\displaystyle=\frac{4\hbar\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{2\pi\sigma_{e}}}\frac{\rho_{B}}{m_{B}c^{2}L}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}(\varepsilon_{e}-\varepsilon_{g})\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\times\left(((\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon((\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
$\displaystyle\times\left((-e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(-e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3})\right)\hat{\rho}_{S}(t)$
(32)

$\displaystyle\int_{0}^{+\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\right)$
$\displaystyle=\frac{4\hbar\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{2\pi\sigma_{e}}}\frac{\rho_{B}}{m_{B}c^{2}L}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}(\varepsilon_{e}-\varepsilon_{g})\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\left(((\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon((\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)\hat{\rho}_{S}(t)$
$\displaystyle\times\left((e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})-(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})-(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
(33)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t)\right)$
$\displaystyle=\frac{4\hbar\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{2\pi\sigma_{e}}}\frac{\rho_{B}}{m_{B}c^{2}L}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}(\varepsilon_{e}-\varepsilon_{g})\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\times\left((-(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(-(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)\hat{\rho}_{S}(t)$
$\displaystyle\times\left((e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
(34)

$\displaystyle\mathrm{Tr}_{B}\left(\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\hat{H}_{SB}(t)\right)$
$\displaystyle=\frac{4\hbar\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{2\pi\sigma_{e}}}\frac{\rho_{B}}{m_{B}c^{2}L}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}(\varepsilon_{e}-\varepsilon_{g})\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
$\displaystyle\times\hat{\rho}_{S}(t)\left(((\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})-e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}-\hat{a}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon((\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})-e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
$\displaystyle\times\left((e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+(\hat{a}^{\dagger}_{g,1}-\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}+\hat{a}_{g,2}))\right.$
$\displaystyle\left.+\epsilon(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})+(\hat{a}^{\dagger}_{g,2}-\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}+\hat{a}_{g,3}))\right)$
(35)

Now let us introduce the following jump operator to simplify the expressions
from Eqs. 32 to 35:

$\hat{c}=(\hat{a}^{\dagger}_{g,1}+\hat{a}^{\dagger}_{g,2})(\hat{a}_{g,1}-\hat{a}_{g,2})+\epsilon(\hat{a}^{\dagger}_{g,2}+\hat{a}^{\dagger}_{g,3})(\hat{a}_{g,2}-\hat{a}_{g,3})$
(36)

Also, we designate $A$ as a constant that denotes the strength of coupling
between the system and the background BEC as follows:

$A=\frac{4\hbar\pi^{2}a^{2}_{SB}}{\mu^{2}\sqrt{2\pi\sigma_{e}}}\frac{\rho_{B}}{m_{B}c^{2}L}\frac{\sigma_{g}}{\sigma_{e}}\left(\frac{\Omega_{1}}{\Delta_{1}}x_{0}\right)^{2}(\varepsilon_{e}-\varepsilon_{g})\exp\left(-\frac{x_{0}^{2}}{\sigma_{e}^{2}}\right)$
(37)

Substituting this into Eqs. 32 to 35, we then obtain

$\displaystyle\int_{0}^{+\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\right)=A(\hat{c}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger})(\hat{c}^{\dagger}-e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c})\hat{\rho}_{S}(t)$
$\displaystyle=A(\hat{c}\hat{c}^{\dagger}\hat{\rho}_{S}(t)-\hat{c}^{\dagger}\hat{c}\hat{\rho}_{S}(t)+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger}\hat{c}^{\dagger}\hat{\rho}_{S}(t)-e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}\hat{c}\hat{\rho}_{S}(t))$
(38)

$\displaystyle\int_{0}^{+\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t)\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\right)=A(\hat{c}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger})\hat{\rho}_{S}(t)(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}-\hat{c}^{\dagger})$
$\displaystyle=A(\hat{c}^{\dagger}\hat{\rho}_{S}(t)\hat{c}-\hat{c}\hat{\rho}_{S}(T)\hat{c}^{\dagger}-e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger}\hat{\rho}_{S}(t)\hat{c}^{\dagger}+e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}\hat{\rho}_{S}(t)\hat{c})$
(39)

$\displaystyle\int_{0}^{+\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{H}_{SB}(t-t^{\prime})\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t)\right)=A(-\hat{c}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger})\hat{\rho}_{S}(t)(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}+\hat{c}^{\dagger})$
$\displaystyle=A(\hat{c}^{\dagger}\hat{\rho}_{S}(t)\hat{c}-\hat{c}\hat{\rho}_{S}(t)\hat{c}^{\dagger}+e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger}\hat{\rho}_{S}(t)\hat{c}^{\dagger}-e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}\hat{\rho}_{S}(t)\hat{c})$
(40)

$\displaystyle\int_{0}^{+\infty}dt^{\prime}\;\mathrm{Tr}_{B}\left(\hat{\rho}_{S}(t)\otimes\hat{\rho}_{B}\hat{H}_{SB}(t-t^{\prime})\hat{H}_{SB}(t)\right)=A\hat{\rho}_{S}(t)(\hat{c}-e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger})(e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}+hat{c}^{\dagger})$
$\displaystyle=A\hat{\rho}_{S}(t)(\hat{c}\hat{c}^{\dagger}-\hat{c}^{\dagger}\hat{c}+e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}\hat{c}-e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\hat{c}^{\dagger}\hat{c}^{\dagger})$
(41)

Finally, substituting Eqs. 38 to 41 into Eq. 20 and evaluating the resulting
expression, we obtain the following form of the master equation governing the
dissipative dynamics due to the interaction between the trapped ultracold atom
system and the background BEC:

$\displaystyle\frac{d}{dt}\hat{\rho}_{S}(t)=A\left((2\hat{c}^{\dagger}\hat{\rho}_{S}(t)\hat{c}-\left\\{\hat{c}\hat{c}^{\dagger},\hat{\rho}_{S}(t)\right\\})-(2\hat{c}\hat{\rho}_{S}(t)\hat{c}^{\dagger}-\left\\{\hat{c}^{\dagger}\hat{c},\hat{\rho}_{S}(t)\right\\})\right.$
$\displaystyle\left.+e^{\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\left[\hat{c}\hat{c},\hat{\rho}_{S}(t)\right]-e^{-\frac{2it}{\hbar}(\varepsilon_{e}-\varepsilon_{g})}\left[\hat{c}^{\dagger}\hat{c}^{\dagger},\hat{\rho}_{S}(t)\right]\right)$
(42)

##  4 Time Evolution of the Trapped Ultracold Atom System Via the Master
Equation

Having derived the master equation that describes the driven - dissipative
dynamics of the trapped ultracold atom system as it interacts with the
background BEC, we now proceed to use this equation to evolve the system over
time, and evaluate some physical properties that are of interest for this
system. Before we do so, however, let us make some notes regarding the
equation itself, in particular the coupling constant $A$, given by Eq. 37.
First, the constant $A$ is known as the coupling constant because its
magnitude specifies how strongly the trapped ultracold atom system interacts
with the background BEC. The greater the magnitude of $A$, the stronger the
interaction between the system and the BEC. At the same time, the explicit
form of the coupling constant shows that the coupling strength between the
trapped ultracold atom system and the background BEC can be adjusted by
adjusting one or more physical parameters governing this coupling, such as the
trap widths at the ground and excited states, $\sigma_{g}$ and $\sigma_{e}$,
as well as the scattering length $a_{S}B$. However, the parameter that is
easiest to access from an experimental standpoint would be the Rabi frequency
$\Omega_{1}$ of the Raman laser coupling the ground states located in the
first and third harmonic traps to the excited state in the second harmonic
trap, which can be done by replacing the Raman laser coupling these energy
levels with another of a different frequency. We note that changing this Rabi
frequency $\Omega_{1}$ also changes the detuning $\Delta_{1}$ of the Raman
laser.

This ability to vary the coupling strength between the trapped ultracold atom
system and the background BEC is crucial to the manner by which the coupled
trapped ultracold atom - background BEC system will be evolved over time. This
is because, following Refs. [34] and [35], over the course of the time
evolution of the trapped ultracold atoms, we will be carrying out a
stroboscopic coupling of the trapped ultracold atoms and the background BEC,
wherein we vary the strength of the coupling between the trapped atoms and the
background BEC at particular instants of time. In particular, we couple the
ground states and the excited states of the ultracold atoms in their
respective trap locations using the Raman lasers while they interact with the
background BEC for only a finite amount of time before we turn off the Raman
lasers and replace one of them, specifically the first Raman laser which
couples the ground states in the first and third harmonic traps to the excited
state of the second excited trap, with another Raman laser with a different
Rabi frequency, afterwhich we turn both lasers on again to couple the energy
levels and repeat the process over a number of intervals of time. This
stroboscopic coupling results in the trapped ultracold atoms evolving via a
combination of driven and dissipative dynamics towards a steady state, as
shown in the figures below.

For our simulations, due to computational power limitations, we consider the
case wherein the trapped ultracold atom gas has an order of magnitude equal to
3. Also, our initial states will be a superposition of four number states,
which has the following explicit form:

$\displaystyle\left|\psi_{0}\right\rangle=\frac{1}{\sqrt{4}}\left(\left(\left|(N_{g,1})_{1}\right\rangle+\left|(N_{g,1})_{2}\right\rangle+\left|(N_{g,1})_{3}\right\rangle+\left|(N_{g,1})_{4}\right\rangle\right)+\left(\left|(N_{g,2})_{1}\right\rangle+\left|(N_{g,2})_{2}\right\rangle\right.\right.$
$\displaystyle\left.\left.+\left|(N_{g,2})_{3}\right\rangle+\left|(N_{g,2})_{4}\right\rangle\right)+\left(\left|(N_{g,3})_{1}\right\rangle+\left|(N_{g,3})_{2}\right\rangle+\left|(N_{g,3})_{3}\right\rangle+\left|(N_{g,3})_{4}\right\rangle\right)\right)$
(43)

Here, the state $\left|(N_{g,j})_{n}\right\rangle$ is an eigenstate of the
particle number operator $\hat{N}_{g,j}=\hat{a}^{\dagger}_{g,j}\hat{a}_{g,j}$
corresponding to the ground state of the system located at the $j$th node with
particle number $(N_{g,j})_{n}$, so that
$\hat{N}_{g,j}\left|(N_{g,j})_{n}\right\rangle=(N_{g,j})_{n}\left|(N_{g,2})_{1}\right\rangle$.
The density matrix corresponding to the initial state will then have the form
$\hat{\rho}_{0}=\left|\psi_{0}\right\rangle\left\langle\psi_{0}\right|$.

###  4.1 $\left\langle N_{1}(t)\right\rangle>\left\langle
N_{2}(t)\right\rangle$ and $\left\langle N_{1}(t)\right\rangle>\left\langle
N_{3}(t)\right\rangle$

We first consider the case wherein most of the trapped ultracold atoms are
initially loaded into the first trap, i. e. the trap where the first ground
state is located, such that $\left\langle N_{1}(t)\right\rangle>\left\langle
N_{2}(t)\right\rangle$ and $\left\langle N_{1}(t)\right\rangle>\left\langle
N_{3}(t)\right\rangle$. For this case, the particle number eigenstates which
are components of the initial state $\left|\psi_{0}\right\rangle$ have
corresponding eigenvalues

$\displaystyle(N_{g,1})_{1}=325,(N_{g,1})_{2}=330,(N_{g,1})_{3}=370,(N_{g,1})_{4}=375,(N_{g,2})_{1}=120,(N_{g,2})_{2}=130,$
$\displaystyle(N_{g,2})_{3}=150,(N_{g,2})_{4}=200,(N_{g,3})_{1}=50,(N_{g,3})_{2}=75,(N_{g,3})_{3}=125,(N_{g,3})_{4}=150$
(44)

These eigenstates, in turn, will correspond to the particle number expectation
values
$\left\langle\hat{N}_{g,1}\right\rangle=350,\left\langle\hat{N}_{g,2}\right\rangle=150,\left\langle\hat{N}_{g,3}\right\rangle=100$
which all add up to a total particle number expectation value of
$\left\langle\hat{N}\right\rangle=\left\langle\hat{N}_{g,1}\right\rangle+\left\langle\hat{N}_{g,2}\right\rangle+\left\langle\hat{N}_{g,3}\right\rangle=600$.
In evolving this initial state, we set the coupling constant in the master
equation used to evolve the state to have the value $A=9.00\times 10^{-}5$,
and the characteristic timescale over which the state is evolved to have the
value $\tau_{E}=0.1$, with the state evolved over 100 timesteps.

Figure 2:  Time evolution of the expectation values of the particle number
$N(t)$ in the first (upper left), second (upper right) and third (lower left)
harmonic traps containing the ground state of the system, calculated using a
superposition of four particle number eigenstates per ground state location,
with corresponding eigenvalues given in Eq. 44.

As can be seen in Fig. 2, for the given parameters of the system, each of the
components of the ultracold atom gas trapped in the three harmonic potentials
corresponding to the ground states of the system will evolve via a combination
of their interaction with the background BEC and their stroboscopic coupling
to the excited states of the system located in neighboring harmonic traps in
such a way that the expectation value of the particle number of these
ultracold atom gas components will approach a steady - state value over time.
In particular, the expectation value of the number of ultracold atoms in the
first harmonic trap will continuously increase until it reaches its steady -
state value, while the expectation value of the number of ultracold atoms in
the second harmonic trap will continuously decrease until it reaches its
steady - state value. On the other hand, the expectation value of the number
of ultracold atoms in the third harmonic trap oscillates until its steady -
state value is attained. This dynamical behavior of the three components of
the ultracold atom gas in their respective harmonic traps can be seen for a
wide range of values of the coupling constant $A$ and for a large set of
combinations of particle number eigenstates $\left|(N_{g,j})_{n}\right\rangle$
that form the initial state. However, as can also be seen in Fig. 2, the
change in the initial number of particles in each trapped component of the
ultracold atom gas is minimal, of the order of magnitude 1 for this example.
But for a given set of particle number eigenstates
$\left|(N_{g,j})_{n}\right\rangle$ and coupling constants $A$, this change in
the particle number can be magnified to be of order of magnitude 2 or greater.

This is shown in this second case case wherein most of the trapped ultracold
atoms are again initially loaded into the first trap, i. e. the trap where the
first ground state is located. But for this case, the particle number
eigenstates which are components of the initial state
$\left|\psi_{0}\right\rangle$ have corresponding eigenvalues

$\displaystyle(N_{g,1})_{1}=425,(N_{g,1})_{2}=426,(N_{g,1})_{3}=575,(N_{g,1})_{4}=574,(N_{g,2})_{1}=20,(N_{g,2})_{2}=25,$
$\displaystyle(N_{g,2})_{3}=80,(N_{g,2})_{4}=70,(N_{g,3})_{1}=20,(N_{g,3})_{2}=25,(N_{g,3})_{3}=80,(N_{g,3})_{4}=75$
(45)

These eigenstates, in turn, will correspond to the particle number expectation
values
$\left\langle\hat{N}_{g,1}\right\rangle=500,\left\langle\hat{N}_{g,2}\right\rangle=50,\left\langle\hat{N}_{g,3}\right\rangle=50$
which all add up to a total particle number expectation value of
$\left\langle\hat{N}\right\rangle=\left\langle\hat{N}_{g,1}\right\rangle+\left\langle\hat{N}_{g,2}\right\rangle+\left\langle\hat{N}_{g,3}\right\rangle=600$.
In evolving this initial state, we set the coupling constant in the master
equation used to evolve the state to have the value $A=1.00\times 10^{-4}$,
and the characteristic timescale over which the state is evolved to have the
value $\tau_{E}=0.1$, with the state evolved over 100 timesteps.

Figure 3:  Time evolution of the expectation values of the particle number
$N(t)$ in the first (upper left), second (upper right) and third (lower left)
harmonic traps containing the ground state of the system, calculated using a
superposition of four particle number eigenstates per ground state location,
with corresponding eigenvalues given in Eq. 45.

As can be seen in Fig. 3, for this combination of parameters, we can see that
the number of particles in the first trap corresponding to the ground state of
the system has increased by an order of magnitude of 2. However, we see a
negligible change in the number of particles in the other two traps
corresponding to the system’s ground state. This then suggests that, due to
the interaction with the background BEC, atoms from the BEC have been
scattered into the trapped ultracold atom gas, mixing with them and becoming
part of the system. One possible mechanism that can facilitate such scattering
from the background BEC into the trapped ultracold atom gas is hyperfine
scattering between background BEC atoms and the trapped ultracold gas’s atoms.
This may not be a problem if both the background BEC and the trapped ultracold
atom gas are composed of atoms of the same species, but in cases wherein the
background BEC and the trapped ultracold atom gas are composed of atoms of
different species, such intermixing of atoms may present a problem especially
if the trapped ultracold atom gas needs to be homogeneous, as is the case if
the trapped ultracold atom gas will eventually undergo Bose - Einstein
condensation. Hence, great care must be exercised in selecting the particle
number eigenstates that will be used to form the initial state describing the
trapped ultracold atom gas.

To further underscore the importance of our choice of particle number
eigenstates that will form our initial state, let us consider the following
figure, which shows how the expectation value of the particle number in trap 1
varies for different combinations of particle number eigenstates used to
construct the component of the initial state corresponding to trap 1. As we
can see in this figure, it is possible to construct an initial state wherein
the particle number expectation value in trap 1 first increases, then
decreases over time until it attains its steady - state value. However, the
order of magnitude of the particle number decrease in this case is around 2,
which is a relatively large amount. On the other hand, a small variation in
the component particle number eigenstates for this initial state will result
in large changes in the dynamical behavior of the system. As shown in the
figure, if we change the initial state in trap 1 from

$\frac{1}{\sqrt{4}}\left(\left|(N_{g,1})_{1}=421\right\rangle+\left|(N_{g,1})_{2}=579\right\rangle+\left|(N_{g,1})_{3}=422\right\rangle+\left|(N_{g,1})_{4}=578\right\rangle\right)$
(46)

to

$\frac{1}{\sqrt{4}}\left(\left|(N_{g,1})_{1}=422\right\rangle+\left|(N_{g,1})_{2}=578\right\rangle+\left|(N_{g,1})_{3}=423\right\rangle+\left|(N_{g,1})_{4}=577\right\rangle\right)$
(47)

then the initial state suddenly changes from a state whose particle number
expectation value decreases by an order of magnitude of 2 as it evolves
towards its steady state value to a state whose particle number expectation
value increases by the same order of magnitude as it evolves over time towards
its steady state value. At the same time, any further changes to the initial
state will not result in a noticeable change in the manner by which it evolves
over time, as can be seen in figure 4, wherein the eigenvalues corresponding
to the component particle number eigenstates of the initial state in trap 1
are changed by a value equal to 1. As the eigenvalues corresponding to the
first and third component particle number eigenstates increase (with a
corresponding decrease in the eigenvalues corresponding to the second and
fourth component particle number eigenstates), the initial state in trap 1
will eventually evolve in such a way that the expectation value of its
particle number increases over time and attains a definite steady - state
value. However, the magnitude of this steady - state value remains constant
even as the eigenvalues corresponding to the first and third component
particle number eigenstates of the initial state in trap 1 continue to
increase in magnitude. This case then shows that it is possible to find a set
of parameters for the initial state of the system and evolve it in such a way
that not only will the expectation value of the number of particles in trap 1
attain a steady - state value irrespective of the combination of particle
number eigenstates that form the component of the system’s initial state in
trap 1, but also that this steady - state expectation value is greater than
its initial value.

Figure 4:  Time evolution of the expectation values of the particle number
$N(t)$ in the first harmonic trap containing the ground state of the system,
calculated using different superpositions of four particle number eigenstates
per ground state location, with corresponding eigenvalues given in the figure.

Now let us consider what happens to the time evolution of the initial state as
the expectation value of the number of particles of the initial state in trap
1, $\left\langle N_{1}\right\rangle$ increases, when most of the particles in
the trapped ultracold atom gas are initially in trap 1. To do this, we need to
compute the particle number fraction, which we define as follows:

$PF(t)=\frac{\left\langle N_{j}(t)\right\rangle-\left\langle
N_{j}\right\rangle}{\left\langle N_{j}\right\rangle}$ (48)

Here, $\left\langle N_{j}(t)\right\rangle$ is the expectation value of the
particle number in trap j, $j=1,2,3$, of the time - evolved state at the
instant of time $t$, while $\left\langle N_{j}\right\rangle$ is the
expectation value of the particle number in trap j for the initial state. As
shown in Fig. 5, as $\left\langle N_{1}\right\rangle$ increases, the particle
number fraction in trap 1 also increases for any given instant of time $t$,
and will eventually approach a steady - state value. This signifies that the
greater the expectation value of the initial state of the system in trap 1,
the greater the number of atoms from the background BEC that can mix with
those from the ultracold atom gas in this trap, consequently increasing the
expectation value of the trapped ultracold atom gas’s particle number over
time.

Figure 5:  Time evolution of the particle number fraction $PF(t)$ of the
ultracold atom system’s component in trap 1, with varying initial values of
the particle number expectation value as indicated in the figure.

Finally, let us see what happens when we adjust the coupling strength between
the trapped ultracold atom system and the background BEC, which is done by
varying the value of the coupling constant $A$ in the master equation. As
shown in Fig. 6, varying the magnitude of the coupling constant will have no
effect on the manner by which the initial state in Trap 1 of the system will
evolve over time; it will still evolve in such a way that $\left\langle
N(t)\right\rangle$ will attain a steady - state value. However, as $A$
decreases, so too will the steady - state value of $\left\langle
N(t)\right\rangle$. This signifies that weakening the coupling between the
trapped ultracold atom system and the background BEC will also reduce the
number of atoms from the background BEC mixing with the trapped ultracold atom
gas if most of the atoms in the trapped ultracold gas are located in trap 1.

Figure 6:  Time evolution of the expectation values of the particle number
$N(t)$ in the first harmonic trap containing the ground state of the system,
calculated with varying values of the coupling constant $A$ in the master
equation, with these values of $A$ indicated in the figure.

###  4.2 $\left\langle N_{2}(t)\right\rangle>\left\langle
N_{1}(t)\right\rangle$ and $\left\langle N_{2}(t)\right\rangle>\left\langle
N_{3}(t)\right\rangle$

Let us now consider the case where most of the particles in the ultracold atom
gas are loaded in trap 2, such that $\left\langle
N_{2}(t)\right\rangle>\left\langle N_{1}(t)\right\rangle$ and $\left\langle
N_{2}(t)\right\rangle>\left\langle N_{3}(t)\right\rangle$. It is important to
note that, from Eq. 36, creation and annihilation operators corresponding to
the location of trap 2 appear twice in the jump operator corresponding to the
driven - dissipative dynamics of the system. This can be explained by the fact
that, as shown in figure 1, there are two excited energy levels, both of the
same magnitude, that are coupled to the ground state energy level located in
trap 2. As such, we expect that if the particle number expectation value at
trap 2 of the initial state of the system is greater than the corresponding
particle number expectation values at the other two traps for the same initial
state, the dynamics of the system will be different from the case wherein the
particle number expectation value at trap 1 of the system’s initial state is
greater than for the other two. This is, indeed the case when the particle
number expectation values for the component particle number eigenstates of the
initial state of the system have the following values:

$\displaystyle(N_{g,1})_{1}=120,(N_{g,1})_{2}=130,(N_{g,1})_{3}=150,(N_{g,1})_{4}=200,(N_{g,2})_{1}=325,(N_{g,2})_{2}=330,$
$\displaystyle(N_{g,2})_{3}=370,(N_{g,2})_{4}=375,(N_{g,3})_{1}=50,(N_{g,3})_{2}=75,(N_{g,3})_{3}=150,(N_{g,3})_{4}=125$
(49)

For this case, the corresponding particle number expectation values in each
trap are
$\left\langle\hat{N}_{g,1}\right\rangle=150,\left\langle\hat{N}_{g,2}\right\rangle=350,\left\langle\hat{N}_{g,3}\right\rangle=100$
which all add up to a total particle number expectation value of
$\left\langle\hat{N}\right\rangle=\left\langle\hat{N}_{g,1}\right\rangle+\left\langle\hat{N}_{g,2}\right\rangle+\left\langle\hat{N}_{g,3}\right\rangle=600$.

Figure 7:  Time evolution of the expectation values of the particle number
$N(t)$ in the first (upper left), second (upper right) and third (lower left)
harmonic traps containing the ground state of the system, calculated using a
superposition of four particle number eigenstates per ground state location,
with corresponding eigenvalues given in Eq. 49.

Now the first thing that we should note for this case is that the coupling
constant $A$ of the master equation used to evolve this state is larger by an
order of magnitude equal to 1 compared to the coupling constant $A$ of the
master equation used to evolve the first two initial states. In particular,
for this case, the coupling constant $A$ has a value $A=1.00\times 10^{-3}$.
This in itself already signifies that the dynamics for this case is already
different from the previous two, due to the coupling being stronger for this
case. It should be noted that if this magnitude of the coupling constant were
used for the previous two cases, with all other parameters for those cases
remaining unchanged, the dynamics of the time - evolved initial state will be
very different from what we have seen in Fig. 7. At the same time, if one is
to compare the dynamics of the time - evolved initial state from Figs. 2 and
7, we find that qualitatively, they may look similar, with the particle number
expectation value in trap 1 steadily increasing until the steady - state value
is reached; the particle number expectation value in trap 2 steadily
decreasing until the steady - state value is reached; and the particle number
expectation value in trap 3 oscillating until its steady - state value is
reached. However, there is a fundamental difference in the magnitude of
increase and decrease in particle number expectation values for traps 1 and 2,
respectively, for these two different initial states. In particular, for trap
1, the magnitude of increase of the particle number expectation value is much
greater for the case where $\left\langle N_{1}(t)\right\rangle>\left\langle
N_{2}(t)\right\rangle$ than in the case where $\left\langle
N_{1}(t)\right\rangle<\left\langle N_{2}(t)\right\rangle$, while for trap 2,
the magnitude of decrease of the particle number expectation value is much
greater for the case where $\left\langle N_{1}(t)\right\rangle<\left\langle
N_{2}(t)\right\rangle$ than in the case where $\left\langle
N_{1}(t)\right\rangle>\left\langle N_{2}(t)\right\rangle$. Furthermore, the
order of magnitude of the difference between the particle number expectation
values in trap 1 for the cases $\left\langle
N_{1}(t)\right\rangle>\left\langle N_{2}(t)\right\rangle$ and $\left\langle
N_{1}(t)\right\rangle<\left\langle N_{2}(t)\right\rangle$ is much less than
the order of magnitude of the difference between the particle number
expectation values in trap 2 for the cases $\left\langle
N_{1}(t)\right\rangle>\left\langle N_{2}(t)\right\rangle$ and $\left\langle
N_{1}(t)\right\rangle<\left\langle N_{2}(t)\right\rangle$.

Now let us consider what happens when we vary the particle number eigenstates
that form the initial state of the system for this case where most of the
particles are to be found in Trap 2. For the case under consideration, the
particle number expectation values in each trap are given as $\left\langle
N_{1}\right\rangle=100,\left\langle N_{2}\right\rangle=500,\left\langle
N_{3}\right\rangle=100$. As shown in Fig. 8, the expectation value of the
number of particles in Trap 2 for this case will vary as we change the
component particle number eigenstates, and their corresponding eigenvalues,
that form the initial state of the system in Trap 2. However, despite this
variation in the component particle number eigenstates, there is no observed
qualitative change in how the trapped ultracold atom system’s component in
Trap 2 evolves over time, as can be seen in Fig.8. In particular, it can be
seen that the magnitude of the particle number expectation value in Trap 2
oscillates, eventually settling towards a central value when
$\left(N_{g,2}\right)_{1}=440,\left(N_{g,2}\right)_{2}=560,\left(N_{g,2}\right)_{3}=441,\left(N_{g,2}\right)_{4}=559$.
However, the manner in which the time - evolved initial state in Trap 2
evolves towards this steady - state remains the same, whereby the expectation
value of its particle number continuously decreases until it reaches the
steady - state value. At the same time, for this case, it can be seen that
even if the eigenvalues corresponding to the first and third component
particle number eigenstates for the initial state in Trap 2 increase in
magnitude, the steady - state attained by the time - evolved state in this
location will remain the same.

Figure 8:  Time evolution of the expectation values of the particle number
$N(t)$ in the second harmonic trap containing the ground state of the system,
calculated using different superpositions of four particle number eigenstates
per ground state location, with corresponding eigenvalues given in the figure.

Looking at how the particle number fraction $PF(t)$ in trap 2 as defined in
Eq. 48 will vary for this case as we increase the expectation number of the
particle number of the initial state in trap 2, we find that, as shown in Fig.
9, $PF(t)$ behaves differently from the previous case. In particular, for this
case, the steady - state value of $PF(t)$ first increases slightly as
$\left\langle N_{2}\right\rangle$ increases from $\left\langle
N_{2}\right\rangle=300$ to $\left\langle N_{2}\right\rangle=400$, then
decreases significantly as $\left\langle N_{2}\right\rangle$ increases from
$\left\langle N_{2}\right\rangle=400$ to $\left\langle
N_{2}\right\rangle=500$. However, for all three cases, just as in the previous
case, $PF(t)$ will again increase with time, eventually attaining a steady -
state value.

Figure 9:  Time evolution of the particle number fraction $PF(t)$ of the
ultracold atom system’s component in trap 2, with varying initial values of
the particle number expectation value as indicated in the figure.

Finally, let us now consider what happens when we vary the coupling constant
$A$ for this case, which consequently varies the coupling strength between the
trapped ultracold atom system and the background BEC. As can be seen in Fig.
10, lowering the value of the coupling constant $A$ will also lower the steady
- state value of the expectation value of the particle number of the time -
evolved component of the trapped ultracold atom state in trap 2, $\left\langle
N_{2}(t)\right\rangle$. However, just as in the previous case, the manner in
which the initial state in trap 2 evolves over time will not change. Hence,
this implies that for the case wherein the expectation value of the particle
number in trap 2 is greater than in the other two traps, weakening the
coupling strength between the background BEC and the trapped ultracold atom
gas will reduce the number of atoms in trap 2 that are lost as the state
evolves over time, while ensuring that this state will still attain a steady -
state over the course of its time evolution.

Figure 10:  Time evolution of the expectation values of the particle number
$N(t)$ in the second harmonic trap containing the ground state of the system,
calculated with varying values of the coupling constant $A$ in the master
equation, with these values of $A$ indicated in the figure.

###  4.3 $\left\langle N_{3}(t)\right\rangle>\left\langle
N_{1}(t)\right\rangle$ and $\left\langle N_{3}(t)\right\rangle>\left\langle
N_{2}(t)\right\rangle$

Having considered the cases wherein $\left\langle
N_{1}\right\rangle>\left\langle N_{2}\right\rangle and\left\langle
N_{1}\right\rangle>\left\langle N_{3}\right\rangle$ as well as $\left\langle
N_{2}\right\rangle>\left\langle N_{1}\right\rangle and\left\langle
N_{2}\right\rangle>\left\langle N_{3}\right\rangle$, we now consider the case
where $\left\langle N_{3}\right\rangle>\left\langle N_{1}\right\rangle
and\left\langle N_{3}\right\rangle>\left\langle N_{2}\right\rangle$. Let us
first look at how each component of the trapped ultracold atom gas in each of
the harmonic potentials of the trap array evolve over time for this case. The
particle number eigenvalues corresponding to the component particle number
eigenstates that are used to construct the components of the initial state of
the ultracold atom gas in traps 1, 2 and 3 of the harmonic trap array have the
following values:

$\displaystyle(N_{g,1})_{1}=20,(N_{g,1})_{2}=80,(N_{g,1})_{3}=25,(N_{g,1})_{4}=75,(N_{g,2})_{1}=20,(N_{g,2})_{2}=80,$
$\displaystyle(N_{g,2})_{3}=25,(N_{g,2})_{4}=75,(N_{g,3})_{1}=414,(N_{g,3})_{2}=586,(N_{g,3})_{3}=415,(N_{g,3})_{4}=585$
(50)

As can be seen in Fig. 11, each of the trapped ultracold atom gas’s components
in traps 1, 2 and 3 will evolve over time in a manner similar to the previous
two cases, which was shown in Figs. 3 and 7. However, unlike the previous two
cases, the change in $\left\langle N_{3}(t)\right\rangle$ at any instant of
time is not large compared to the variation in $\left\langle
N_{2}(t)\right\rangle$ and $\left\langle N_{1}(t)\right\rangle$. In fact, the
oscillatory behavior for $\left\langle N_{3}(t)\right\rangle$ is still present
for this case, just as it was present in the previous two cases, which is
characteristic of the time evolution of components of the trapped ultracold
atom gas in trap 3. Furthermore, for this case, the resulting steady - state
value for $\left\langle N_{3}(t)\right\rangle$ is very close to, if not equal
to, its initial value $\left\langle N_{3}\right\rangle$ despite fluctuating in
value before settling to its steady - state value, which is the same behavior
shown in the previous two cases for the time - evolved component of the
trapped ultracold atom gas in trap 3.

Figure 11:  Time evolution of the expectation values of the particle number
$N(t)$ in the first (upper left), second (upper right) and third (lower left)
harmonic traps containing the ground state of the system, calculated using a
superposition of four particle number eigenstates per ground state location,
with corresponding eigenvalues given in Eq. 50.

Having seen how the trapped ultracold atom gas evolves over time in each of
the harmonic potentials of the trap array for this case, let us now see what
the effect of varying the component particle number eigenstates of the initial
state corresponding to the component of the ultracold atom gas in Trap 3 will
be on its time evolution and steady - state attained. As can be seen in Fig.
12, varying the initial state of the system in Trap 3 will not change the
manner by which the initial state of the ultracold atom gas’s component in
trap 3 evolves over time, but it will change the steady - state value of the
particle number expectation value for this component of the trapped ultracold
atom gas. In particular, from a steady-state value of $\left\langle
N_{3}(t)\right\rangle\approx 800$ when the initial state in trap 3 has the
form

$\frac{1}{\sqrt{4}}\left(\left|(N_{g,3})_{1}=411\right\rangle+\left|(N_{g,3})_{2}=589\right\rangle+\left|(N_{g,3})_{3}=412\right\rangle+\left|(N_{g,3})_{4}=588\right\rangle\right),$
(51)

this steady - state value decreases to $\left\langle
N_{3}(t)\right\rangle\approx 470$ when the initial state in trap 3 has the
form

$\frac{1}{\sqrt{4}}\left(\left|(N_{g,3})_{1}=412\right\rangle+\left|(N_{g,3})_{2}=588\right\rangle+\left|(N_{g,3})_{3}=413\right\rangle+\left|(N_{g,3})_{4}=587\right\rangle\right)$
(52)

Eventually, as the value of $(N_{g,3})_{1}$ and $(N_{g,3})_{3}$ increase, the
steady - state value of $\left\langle N_{3}(t)\right\rangle$ will eventually
approach a value of $\left\langle N_{3}(t)\right\rangle=500$, which is the
particle number expectation value of the initial state of the ultracold atom
gas’s component localized in trap 3 of the array.

Figure 12:  Time evolution of the expectation values of the particle number
$N(t)$ in the first harmonic trap containing the ground state of the system,
calculated using different superpositions of four particle number eigenstates
per ground state location, with corresponding eigenvalues given in the figure.

Having seen the effect of varying the components of the initial state of the
component of the trapped ultracold atom gas in trap 3, let us now see what
happens when we vary the particle number expectation value of the initial
state of this component of the ultracold atom gas. Following our analysis of
the effects of varying $\left\langle N_{1}\right\rangle$ and $\left\langle
N_{2}\right\rangle$ when $\left\langle N_{1}\right\rangle>\left\langle
N_{2}\right\rangle and\left\langle N_{1}\right\rangle>\left\langle
N_{3}\right\rangle$ and $\left\langle N_{2}\right\rangle>\left\langle
N_{1}\right\rangle and\left\langle N_{2}\right\rangle>\left\langle
N_{3}\right\rangle$, respectively, we take a look at what happens to the
particle fraction of the time - evolved component of the trapped ultracold
atom gas in trap 3 for this case, with the particle fraction given by Eq. 48,
with $j=3$. As shown in Fig. 13, the manner by which $PF(t)$ evolves for this
case will not change as $\left\langle N_{3}\right\rangle$ increases; rather,
it is the maximum value and steady - state value of $PF(t)$ for the time -
evolved component of the trapped ultracold atom gas in trap 3 that changes,
both increasing with $\left\langle N_{3}\right\rangle$.

Figure 13:  Time evolution of the particle number fraction $PF(t)$ of the
ultracold atom system’s component in trap 3, with varying initial values of
the particle number expectation value as indicated in the figure.

Finally, let us examine the effect of weakening the coupling between the
trapped ultracold atom system and the background BEC for this case by
decreasing the magnitude of the coupling constant, $A$. As can be seen in Fig.
14, decreasing the coupling constant $A$’s magnitude will increase the
magnitude of the steady - state value of the particle number expectation value
$\left\langle N_{3}(t)\right\rangle$ of the time-evolved component of the
trapped ultracold atom system in trap 3 for this case. At the same time, this
is accompanied by a corresponding decrease in the maximum value of
$\left\langle N_{3}(t)\right\rangle$ for this system. It is of interest to
note that the manner in which the component of the trapped ultracold atom
system in trap 3 evolves over time does not change even as the magnitude of
the coupling constant $A$ changes. In particular, the instant of time when
$\left\langle N_{3}(t)\right\rangle$ attains its maximum value does not, from
Fig. 14, change even as $A$ changes, and the oscillations present in
$\left\langle N_{3}(t)\right\rangle$ still occur at the same instants of time
no matter the magnitude of $A$, though their amplitude will change as $A$’s
magnitude changes.

Figure 14:  Time evolution of the expectation values of the particle number
$N(t)$ in the second harmonic trap containing the ground state of the system,
calculated with varying values of the coupling constant $A$ in the master
equation, with these values of $A$ indicated in the figure.

##  5 Discussion

Based on the results presented, we can see that, for an ultracold atom gas
trapped in an array of 5 harmonic potentials arranged beside each other, which
is coupled to a background BEC, it will evolve in such a way that if the
trapped ultracold atom gas is trapped in 3 of these 5 potentials and occupy
the ground state energies in these traps, and are in addition coupled to the
excited state energies in the other two traps, then the combination of the
driving of these ultracold atoms from the ground state to the excited state
and the dissipation of Bogoliubov excitations into the background BEC by these
driven ultracold atoms to enable them to return to the ground state energy
will result in the ultracold atom gas evolving over time in such a way that
the expectation value of the particle number of each component of the
ultracold atom gas occupying one of the harmonic potential traps in the trap
array will approach a steady - state value. At the same time, the manner by
which the trapped ultracold atom gas evolves over time, as well as the
particle number expectation value at the system’s steady state, will depend on
the particle number expectation value of the initial state of the system at
each of the harmonic traps in the array, on the particle number eigenstates
which form the initial state of the trapped ultracold atom system, and on the
strength of coupling between the trapped ultracold atom system and the
background BEC.

However, what seems to be surprising among the results obtained is the
distinctive manner by which the ultracold atom states in each harmonic trap in
the array evolve towards their respective steady states. Notably, the results
show that mixing between the atoms of the trapped ultracold atom gas and the
background BEC can occur, increasing the expectation value of the particle
number of the trapped ultracold atom gas, so long as the initial state of the
system has more atoms in the first trap than in the other two. At the moment,
the author cannot yet explain why indeed this is the case, and doing so is
currently beyond the scope of this research work. However, the author would
like to note that this behavior is a matter of interest emerging from this
work, and being able to explain why this behavior in this system emerges in
future research work may actually tell us something new about how two
ultracold atom gases interact with each other.

Nonetheless, the emergence of steady states with respect to the expectation
value of the particle number for the trapped ultracold atom system undergoing
driven - dissipative dynamics as described above implies that this system can
possibly be used in the preparation of multiple ultracold atom states
requiring a fixed particle number, in particular BECs. This dissipative
quantum state preparation scheme becomes especially important when we consider
that novel applications of BECs continue to emerge, such as for fundamental
tests of quantum gravity [36, 37, 38] and for quantum computing and quantum
simulation [39, 40, 41, 42]. The need for BECs for these and other possible
applications that are of interest in quantum mechanics and quantum technology,
now and in the future, will ensure that dissipative quantum state preparation
schemes for constant particle number states such as the one formulated in this
paper, as well as other quantum state preparation schemes for these types of
states, will continue to be of interest among theorists and experimentalists
working in the field of ultracold atoms and quantum technology.

## Acknowledgments

The author would like to thank the Research Center for the Natural and Applied
Sciences (RCNAS) and the College of Science of the University of Santo Tomas
for financial support during the conduct of this research. Part of this work
was conducted when the author was a visiting researcher at the Institute for
Basic Science - Center for Theoretical Physics of Complex Systems (PCS - IBS),
Daejeon, Republic of Korea.

## References

  * Breuer and Petruccione [2007] H. P. Breuer, F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, 2007.
  * Breuer et al. [2016] H. P. Breuer, E. M. Laine, J. Piilo, B. Vacchini,  Rev. Mod. Phys. 88 (2016) 021002.
  * de Vega and Alonso [2017] I. de Vega, D. Alonso,  Rev. Mod. Phys. 89 (2017) 015001.
  * Verstraete et al. [2009] F. Verstraete, M. M. Wolf, J. I. Cirac,  Nat. Phys. 5 (2009) 633–636.
  * Kastoryano et al. [2013] M. J. Kastoryano, M. M. Wolf, J. Eisert,  Phys. Rev. Lett. 110 (2013) 110501.
  * Kordas et al. [2015] G. Kordas, D. Witthaut, P. Buonsante, A. Vezzani, R. Burioni, A. Karanikas, S. Wimberger,  Eur. Phys. J. Special Topics 224 (2015) 2127–2171.
  * Reiter et al. [2017] F. Reiter, A. Sørensen, P. Zoller, C. Muschik,  Nat. Comm. 8 (2017) 1822.
  * Nosov et al. [2023] P. A. Nosov, D. S. Shapiro, M. Goldstein, I. S. Burmistrov,  Phys. Rev. B 103 (2023) 174312.
  * Kastoryano et al. [2011] M. J. Kastoryano, F. Reiter, A. S. Sørensen,  Phys. Rev. Lett. 106 (2011) 090502.
  * Watanabe and Mäkelä [2012] G. Watanabe, H. Mäkelä,  Phys. Rev. A 85 (2012) 023604.
  * Stannigel et al. [2012] K. Stannigel, P. Rabl, P. Zoller,  New J. Phys. 14 (2012) 063014.
  * Gong and Yao [2013] Z. R. Gong, W. Yao,  Phys. Rev. A 87 (2013) 032314.
  * Sweke et al. [2013a] R. Sweke, I. Sinayskiy, F. Petruccione,  J. Phys. B: At. Mol. Opt. Phys. 46 (2013a) 104004\.
  * Sweke et al. [2013b] R. Sweke, I. Sinayskiy, F. Petruccione,  Phys. Rev. A 87 (2013b) 042323.
  * Su et al. [2014] S.-L. Su, X.-Q. Shao, H.-F. Wang, S. Zhang,  Sci. Rep 4 (2014) 07566.
  * Shao et al. [2017] X. Q. Shao, J. H. Wu, X. X. Yi, G.-L. Long,  Phys. Rev. A 96 (2017) 062315.
  * Li and Shao [2019] D. X. Li, X. Q. Shao,  Phys. Rev. A 99 (2019) 032348.
  * Kouzelis et al. [2020] A. Kouzelis, K. Macieszczak, J. Minář, I. Lesanovsky,  Phys. Rev. A 101 (2020) 043847.
  * Zheng et al. [2021] R.-H. Zheng, Y. Xiao, S.-L. Su, Y.-H. Chen, Z.-C. Shi, J. Song, Y. Xia, S.-B. Zheng,  Phys. Rev. A 103 (2021) 052402.
  * Cole et al. [2021] D. C. Cole, J. J. Wu, S. D. Erickson, P.-Y. Hou, A. C. Wilson, D. Leibfried, F. Reiter,  New J. Phys. 23 (2021) 073001.
  * Sharma and Mueller [2021] V. Sharma, E. J. Mueller,  Phys. Rev. A 103 (2021) 043322.
  * Seetharam et al. [2022] K. Seetharam, A. Lerose, R. Fazio, J. Marino,  Phys. Rev. Research. 4 (2022) 013089.
  * Marino [2022] J. Marino,  Phys. Rev. Lett. 129 (2022) 050603.
  * Yang et al. [2023] F. Yang, P. Molignini, E. J. Berghlotz,  Phys. Rev. Res. 5 (2023) 043229.
  * Vorberg et al. [2013] D. Vorberg, W. Wustmann, R. Ketzmerick, A. Eckardt,  Phys. Rev. Lett. 111 (2013) 240405.
  * Zhang and Bartel [2024] Y. Zhang, T. Bartel,  Phys. Rev. A 109 (2024) L021301.
  * Pethick and Smith [2008] C. J. Pethick, H. Smith, Bose - Einstein Condensation in Dilute Gases, Cambridge University Press, 2008.
  * Pitaevskii and Stringari [2016] L. Pitaevskii, S. Stringari, Bose - Einstein Condensation and Superfluidity, Oxford University Press, 2016.
  * Caspar et al. [2016a] S. Caspar, F. Hebenstreit, D. Mesterhazy, U.-J. Wiese,  Phys. Rev. A 93 (2016a) 021602(R).
  * Caspar et al. [2016b] S. Caspar, F. Hebenstreit, D. Mesterhazy, U.-J. Wiese,  New J. Phys. 18 (2016b) 073015.
  * Walker et al. [2018] B. T. Walker, L. C. Flatten, H. J. Hesten, F. Mintert, D. Hunger, A. A. P. Trichet, J. M. Smith, R. A. Nyman,  Nat. Phys. 14 (2018) 1173.
  * Daley et al. [2004] A. J. Daley, P. O. Fedichev, P. Zoller,  Phys. Rev. A 69 (2004) 022306.
  * Diehl et al. [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, P. Zoller,  Nature Phys. 4 (2008) 878.
  * Caballar et al. [2014] R. C. F. Caballar, S. Diehl, H. Mäkelä, M. Oberthaler, G. Watanabe,  Phys. Rev. A 89 (2014) 013620.
  * Caballar [2024] R. C. F. Caballar,  Ann. Phys. 462 (2024) 169609.
  * Howl et al. [2021] R. Howl, V. Vedral, D. Naik, M. Christodoulou, C. Rovelli, A. Iyer,  PRX Quantum 2 (2021) 010325.
  * Overstreet et al. [2023] C. Overstreet, J. Curti, M. Kim, P. Asenbaum, M. Kasevich, F. Giacomini,  Phys. Rev. D 108 (2023) 084038.
  * Thompson et al. [2023] R. J. Thompson, D. C. Aveline, S.-W. Chiow, E. R. Elliot, J. R. Kellogg, J. M. Kohel, M. S. Sbroscia, C. Schneider, J. R. Williams, N. Lundblad,  Quantum Sci. Technol. 8 (2023) 024004.
  * Byrnes et al. [2012] T. Byrnes, K. Wen, Y. Yamamoto,  Phys. Rev. A 85 (2012) 040306(R).
  * Ghasemian [2023] E. Ghasemian,  J. Opt. Soc. Am. B 40(2) (2023) 247–259.
  * Fujii [2023] G. Fujii,  Phys. Scr. 98 (2023) 065109.
  * Bello et al. [2022] M. Bello, G. Platero, A. González-Tudela,  PRX Quantum 3 (2022) 010336.