*Article* **Non-Markovian Quantum Dynamics in a Squeezed Reservoir**

**Valentin Link 1,†, Walter T. Strunz <sup>1</sup> and Kimmo Luoma 1,2,\*,†**


**Abstract:** We study non-Markovian dynamics of an open quantum system system interacting with a nonstationary squeezed bosonic reservoir. We derive exact and approximate descriptions for the open system dynamics. Focusing on the spin boson model, we compare exact dynamics with Redfield theory and a quantum optical master equation for both short and long time dynamics and in non-Markovian and Markov regimes. The squeezing of the bath results in asymptotic oscillations in the stationary state, which are captured faithfully by the Redfield master equation in the case of weak coupling. Furthermore, we find that the bath squeezing direction modifies the effective system–environment coupling strength and, thus, the strength of the dissipation.

**Keywords:** open quantum systems; squeezed states; non-Markovian dynamics

#### **1. Introduction**

In 1926 and 1927, two families of quantum states for the quantum harmonic oscillator were proposed by Schrödinger [1] and Kennard [2], respectively. The first family of states included coherent states [3–5] and the second family included squeezed states [6]. These can be distinguished by considering the Heisenberg uncertainty principle:

$$
\Delta x^2 \Delta p^2 \ge \frac{\hbar^2}{4},
$$

which was discovered also in 1927 [7]. The coherent states saturate this uncertainty principle with equal variances Δ*x*<sup>2</sup> = Δ*p*<sup>2</sup> = *h*¯ /2 in both quadratures. This makes them the closest to points in phase space; therefore, coherent states are often considered the most classical quantum states. In contrast, for a squeezed state, the variances of one of the quadratures can be smaller that *h*¯ /2. The state is, thus, squeezed along a certain direction in phase space. According to the Heisenberg uncertainty principle, the variance of the other quadrature must then be larger than *h*¯ /2. Often, squeezed states are regarded as nonclassical states [8,9]. In fact, so-called two mode squeezing generates entanglement between two oscillators [10]. In contrast to coherent states, squeezed states are not stationary states of the harmonic oscillator.

Ensembles of harmonic oscillators are commonly considered as quantum environments of open quantum systems. In most applications, these environments are stationary, meaning that the initial state of the reservoir oscillators is a stationary state, for instance, a thermal state (which becomes a coherent state for zero temperature). In a large thermal bath, the system is expected to reach thermal equilibrium [11]. If instead the bath can be prepared in a squeezed state, one directly violates stationarity, resulting in a a breakdown of equilibration.

This can have a drastic impact on the dissipation induced by the bath. In 1986, it was proposed that a nonstationary reservoir consisting of harmonic oscillators prepared in

**Citation:** Link, V.; Strunz, W. T.; Luoma, K. Non-Markovian Quantum Dynamics in a Squeezed Reservoir. *Entropy* **2022**, *24*, 352. https:// 10.3390/e24030352

Academic Editor: Bassano Vacchini, Andrea Smirne and Nina Megier

Received: 15 February 2022 Accepted: 22 February 2022 Published: 28 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

squeezed states can lead to an inhibition of phase decay of an atom [12]. Only in 2013, this prediction was experimentally verified using superconducting qubits [13].

Previous theoretical works on squeezed reservoirs were based on quantum master equations with constant coefficients [14–19]. The validity of these equations requires severe assumptions on the system and bath coupling strength and possibly also a separation of time scales such that a rotating wave approximation can be performed. For this reason, such master equations are unable to capture certain phenomena in an exact manner. For instance, they can predict that the very short time dynamics of an observable is linear in time, whereas the Schrödinger equation actually yields quadratic short time dynamics [20,21].

In this article, we investigate open quantum system dynamics in a non-stationary reservoir consisting of oscillators prepared in two-mode squeezed states (broad band squeezed reservoir) [18]. We derive an exact description of the reduced state dynamics using non-Markovian quantum state diffusion (NMQSD) [22–24]. Then, using NMQSD as a starting point, we derive a Hierarchy of Equations of Motion [25] (HEOM) for the density matrix of the open system based on the Hierarchy of stochastic Pure States (HOPS) [26]. For the example of a single two level system, we compare the short and long time dynamics of the numerically exact HEOM theory with a weak coupling perturbation theory, so-called Redfield theory, and the commonly used standard master equation with time independent coefficients. Our main findings are that, in a parameter regime where a quantum optical master equation would work fine for stationary reservoirs, it can not capture accurately short and long time dynamics. Remarkably, if the bath memory time is short the Redfield theory performs extremely well, matching with exact dynamics and predicting accurate long time dynamics.

#### **2. Model**

We consider the following commonly used model for open quantum systems:

$$H = H\_{\mathbb{S}} + L \sum\_{\lambda} \mathbb{g}\_{\lambda} (a\_{\lambda} + a\_{\lambda}^{\dagger}) + \sum\_{\lambda} \omega\_{\lambda} a\_{\lambda}^{\dagger} a\_{\lambda},\tag{1}$$

where the environment consists of independent harmonic oscillators (bosonic modes) [*aλ*, *a*† *<sup>μ</sup>*] = *δλμ* [19]. The system Hamiltonian *HS* and the system coupling operator are left arbitrary, except that we demand *L* = *L*†. Moving to an interaction representation with respect to the free bath Hamiltonian, one obtains the following.

$$H(t) = H\_S + LB(t), \qquad B(t) = \sum\_{\lambda} \mathbf{g}\_{\lambda} (\mathbf{e}^{-i\omega\_{\lambda}t} a\_{\lambda} + \mathbf{e}^{i\omega\_{\lambda}t} a\_{\lambda}^{\dagger}) \tag{2}$$

If the bath's initial state is Gaussian and the coupling of the system to the bath is linear, the response of the bath to the system is fully characterized by the first and second moment of the bath response operator *B*(*t*). Typically, and without loss of generality, one considers *B*(*t*) = 0. Then, the so-called bath correlation function (BCF) is the only relevant quantity describing the influence of the bath [27]:

$$
\pi(t, s) = \langle B(t)B(s) \rangle,\tag{3}
$$

where the expectation value is with respect to the bath initial state. As an initial condition, for the bath, usually a stationary state of the free bath Hamiltonian is considered. For example, a zero temperature bath is described by vacuum state |0, which is a Gaussian state; hence, it is fully characterized by the following correlations.

$$\begin{aligned} \langle 0|a\_{\lambda}|0\rangle &= 0, \\ \langle 0|a\_{\lambda}a\_{\mu}^{\dagger}|0\rangle &= \delta\_{\lambda\mu}, \\ \langle 0|a\_{\lambda}^{\dagger}a\_{\mu}|0\rangle &= 0, \\ \langle 0|a\_{\lambda}a\_{\mu}|0\rangle &= 0 \end{aligned} \tag{4}$$

This results in a stationary BCF that depends on the time difference only.

$$\varkappa(t,s) \equiv \varkappa\_0(t-s) = \sum\_{\lambda} \mathcal{g}\_{\lambda}^2 \mathbf{e}^{-i\omega\_\lambda(t-s)}\tag{5}$$

An alternative characterization of a stationary bath is given in terms of spectral density:

$$J(\omega) = \frac{1}{\pi} \sum\_{\lambda} \mathbf{g}\_{\lambda}^{2} \delta(\omega - \omega\_{\lambda}),\tag{6}$$

which is the Fourier transform of the stationary BCF.

In this article, we are interested in the situation where the initial state of the bath is a squeezed vacuum. In particular, we consider two-mode squeezing which is symmetric around a reference frequency *ω*0. We denote our squeezed vacuum state by |*φ* = *S*|0, where *<sup>S</sup>* is a unitary squeezing operator. If we order the bath modes according to *<sup>ω</sup>*2*λ*0−*<sup>λ</sup>* = 2*ω*<sup>0</sup> − *ωλ*, the state |*φ* is characterized by the following correlations.

$$\begin{aligned} \langle \phi | a\_{\lambda} | \phi \rangle &= 0, \\ \langle \phi | a\_{\lambda} a\_{\mu}^{\dagger} | \phi \rangle &= \iota^{2} \delta\_{\lambda \mu \prime} \\ \langle \phi | a\_{\lambda}^{\dagger} a\_{\mu} | \phi \rangle &= | \upsilon |^{2} \delta\_{\lambda \mu \prime} \\ \langle \phi | a\_{\lambda} a\_{\mu} | \phi \rangle &= -\upsilon \iota \delta\_{\mu 2 \lambda\_{0} - \lambda} \end{aligned} \tag{7}$$

Because *S* is a unitary operator, the squeezing parameters *u* ∈ R and *v* ∈ C satisfy *<sup>u</sup>*<sup>2</sup> <sup>−</sup> <sup>|</sup>*v*<sup>|</sup> <sup>2</sup> = 1. We also have assumed that the two mode squeezing is homogeneous, i.e., *u* and *v* are the same for all modes *λ*. This is often called broadband squeezing. In order to parametrize the squeezing, we introduce real variables *r* and *ϕ* and write

$$u = \cosh(r), \qquad v = \cosh(r)\mathbf{e}^{i\boldsymbol{\wp}},\tag{8}$$

where *r* is squeezing strength, and *φ* is the squeezing direction in phase space. *φ* = 0 corresponds to squeezing the *p*-quadrature, whereas *φ* = *π* corresponds to squeezing the *x*-quadrature. In case of no squeezing *r* = 0, we recover the vacuum correlations (4). To continue, we further assume that coupling constants *g<sup>λ</sup>* have a symmetry property *<sup>g</sup>*2*λ*0−*<sup>λ</sup>* = *<sup>g</sup>λ*, which implies that the spectral density is symmetric around *<sup>ω</sup>*0. From this, we obtain that the bath correlation function has the following structure:

$$a(t,s) = a\_0(t-s)(\mu^2 - vue^{-2i\omega\_0 s} - v^\*\mu e^{2i\omega\_0 t}) + a\_0^\*(t-s)|v|^2,\tag{9}$$

where *α*0(*t* − *s*) is the zero temperature bath correlation function from Equation (5). It is easy to check that this is a valid BCF obeying *α*(*t*,*s*) = *α*∗(*s*, *t*). Note that this function does depend explicitly on both *t* and *s* and not only on their difference. This is because the bath initial state is not stationary with respect to the bath Hamiltonian. We assume the following model for the vacuum BCF:

$$
\alpha\_0(t-s) = \frac{\gamma \Gamma}{2} e^{-\Gamma|t-s| - i\omega\_R(t-s)} \, \_\prime \tag{10}
$$

which corresponds to a continuous bath with Lorentzian spectral density.

$$J(\omega) = \frac{\gamma}{2} \frac{\Gamma^2}{\Gamma^2 + (\omega\_B - \omega)^2} \tag{11}$$

This bath satisfies the required symmetry if *ω<sup>B</sup>* = *ω*0, which we assume in the following. Notably, for such a bath, the white noise limit *α*0(*t* − *s*) → *γδ*(*t* − *s*) exits when Γ → ∞.

The non-Markovian open system dynamics of this model can be described with non-Markovian quantum state diffusion (NMQSD) [22,24,28]. NMQSD is a stochastic unraveling of reduced open system dynamics in terms of a Gaussian colored noise process *z*(*t*) with statistics *E*[*z*(*t*)*z*∗(*s*)] = *α*(*t*,*s*) and *E*[*z*(*t*)] = *E*[*z*(*t*)*z*(*s*)] = 0. The system state is obtained as the ensemble average of stochastic pure states |*ψ*(*z*∗, *t*), which depend on this noise process.

$$\rho(t) = E[|\psi(z^\*, t)\rangle\langle\psi(z^\*, t)|]\tag{12}$$

The stochastic states obey the NMQSD evolution equation.

$$\left|\partial\_{t}\vert\psi(z^{\*},t)\right\rangle = \left(-\mathrm{i}H\_{\mathrm{S}} + Lz^{\*}(t)\right)\vert\psi(z^{\*},t)\rangle - L\int\_{0}^{t}\mathrm{d}s\,a(t,s)\frac{\delta}{\delta z^{\*}(s)}\vert\psi(z^{\*},t)\rangle\tag{13}$$

Different solution strategies for this equation exits, notably the *O*-operator method [23,29] and the exact HOPS method [26,30,31]. Here, we use NMQSD as a tool to derive perturbative and exact master equations for the squeezed bath problem (9).

Although most results are very general, an example we will consider is the spin boson model. For this model, the system is a simple two level system with the following operators

$$H\_S = \frac{\Omega}{2} \sigma\_{z\prime} \qquad L = \sigma\_x. \tag{14}$$

#### **3. Perturbative Master Equations**

If the coupling to the bath is weak and/or the BCF decays rapidly, we can make a perturbative approximation to NMQSD. The lowest order perturbative equation is given by the following.

$$\begin{split} \int\_{0}^{t} \mathrm{d}s \, a(t,s) \frac{\delta}{\delta z^{\*}(s)} |\psi(z^{\*},t)\rangle &= O(z^{\*},t) |\psi(z^{\*},t)\rangle, \\ O(z^{\*},t) \approx \int\_{0}^{t} \mathrm{d}s \, a(t,s) \mathrm{e}^{-\mathrm{i}H\_{\mathrm{S}}(t-s)} \mathrm{d}\mathbf{e}^{\mathrm{i}H\_{\mathrm{S}}(t-s)}. \end{split} \tag{15}$$

Because in this approximation the operator *O*¯ does not depend on the stochastic process, a master equation can be directly derived from NMQSD as follows [23].

$$\begin{split} \partial\_t \rho(t) &= -i[H\_{\mathbb{S}}, \rho(t)] + [L, \rho(t)\mathcal{O}^\dagger(t)] + [\mathcal{O}(t)\rho(t), L] \\ &:= \mathcal{L}\_t^{\mathbb{R}}(\rho(t)) \end{split} \tag{16}$$

This is the well known Redfield master equation [32].

Under certain conditions and in a frame rotating with the frequency *ω*0, one can derive a master equation with constant coefficients starting from the Redfield equation. The resulting equation is well known from quantum optics textbooks [17,18]. For a derivation, we assume the spin boson model with *L* = *σ<sup>x</sup>* and *HS* = *ωσz*/2. In particular, we consider the case where Ω = *ω*<sup>0</sup> + *δ* and *ω*<sup>0</sup> defines the fastest timescale *ω*<sup>0</sup> Γ *δ* and *ω*<sup>0</sup> *γ*. The state in the rotating frame is given by the following.

$$\rho(t) = \mathcal{R}(t)\rho(t)\mathcal{R}^\dagger(t), \qquad \mathcal{R}(t) = \exp\left(\mathrm{i}\omega\_0 \frac{\sigma\_z}{2}t\right) \tag{17}$$

To obtain a master equation with constant coefficients, analogous to the standard quantum optical master equation [32], one explicitly computes the integral in (15) and neglects exponentially decaying terms as well as terms of order *δ*/Γ. Plugging the result into the Redfield Equation (16) in the rotating frame, one can identify counter-rotating terms that oscillate at frequency 2*ω*0. Under the above assumptions, these terms can be neglected and one obtains a master equation with time-independent coefficients.

$$\begin{split} \partial\_t \bar{\rho}(t) &= -i\frac{\delta}{2} [\sigma\_z, \bar{\rho}(t)] \\ &+ \gamma \mu^2 \Big( \sigma\_- \bar{\rho}(t) \sigma\_+ - \frac{1}{2} \{ \sigma\_+, \sigma\_-, \bar{\rho}(t) \} \Big) \\ &+ \gamma |\upsilon|^2 \Big( \sigma\_+ \bar{\rho}(t) \sigma\_- - \frac{1}{2} \{ \sigma\_-, \sigma\_+, \bar{\rho}(t) \} \Big) \\ &+ \gamma \mu \upsilon \sigma\_- \bar{\rho}(t) \sigma\_- + \gamma \mu \upsilon^\* \sigma\_+ \bar{\rho}(t) \sigma\_+ \\ &:= \mathcal{L}^M(\bar{\rho}(t)) \end{split} \tag{18}$$

Note that the coefficients of the last two terms become time dependent if one moves back to the laboratory frame.

In order to assess the validity of these perturbative master equations, we have to compare their predictions with exact results. We discuss in the following how the exact reduced dynamics can be computed using non-Markovian open system methods.

#### **4. Exact Method**

We can utilize non-Markovian open system methods to compute the exact reduced dynamics in the squeezed reservoir. In particular, due to the exponential form of the bath correlation function, we can easily generalize hierarchical methods for the squeezed bath. To this aim, we decompose the BCF in the following manner:

$$u(t,s) = u\_1(t,s) + u\_2(t,s) \tag{19}$$

$$u\_1(t,s) = \frac{\gamma \Gamma}{2} (u^2 - vu e^{-2i\omega\_0 s}) e^{-i\omega\_0 (t-s) - \Gamma |t-s|} \tag{20}$$

$$w\_2(t,s) = \frac{\gamma \Gamma}{2} (|v|^2 - v^\* \mu \mathbf{e}^{2i\omega\_0 s}) \mathbf{e}^{i\omega\_0 (t-s) - \Gamma |t-s|} \tag{21}$$

To solve the linear NMQSD Equation (13), we further define the corresponding functional differential operators:

$$\mathcal{D}\_{j} = \int\_{0}^{t} \mathrm{d}s \, a\_{j}(t,s) \frac{\delta}{\delta z\_{s}^{\*}} \, \qquad j = 1,2,\tag{22}$$

so that (13) becomes the following.

$$\left|\partial\_{\mathbf{i}}\vert\psi\rangle = -\mathbf{i}H\vert\psi\rangle + L z\_{\mathbf{i}}^{\*}\vert\psi\rangle - L \sum\_{j=1,2} \mathcal{D}\_{j}\vert\psi\rangle\tag{23}$$

Because of the exponential form of the kernels (20), we can compute the time derivative of the functional differential operators.

$$
\partial\_1 \mathcal{D}\_{\dot{\jmath}} = -\mathcal{W}\_{\dot{\jmath}} \mathcal{D}\_{\dot{\jmath}} + a\_{\dot{\jmath}}(t, t) \frac{\delta}{\delta z\_t^\*}, \qquad \mathcal{W}\_1 = \mathcal{W}\_2^\* = \mathrm{i}\omega \mathbf{\hat{u}} + \Gamma \tag{24}
$$

This allows employing the hierarchy of pure states (HOPS) scheme to solve NMQSD. In particular, one defines an 'auxiliary state' for every double-index *<sup>n</sup>* <sup>∈</sup> <sup>N</sup><sup>2</sup> 0.

$$<\langle \psi^{\mathfrak{u}} \rangle = (\mathcal{D}\_1)^{\mathfrak{u}\_1} (\mathcal{D}\_2)^{\mathfrak{u}\_2} |\psi\rangle \tag{25}$$

These states obey the hierarchy of pure states equation of motion.

$$\partial\_t |\boldsymbol{\psi}^{\boldsymbol{\mu}}\rangle = \left(-\mathrm{i}H\_S + Lz\_1^\* - \sum\_{j=1,2} n\_j \mathcal{W}\_j\right) |\boldsymbol{\psi}^{\boldsymbol{\mu}}\rangle - L \sum\_{j=1,2} |\boldsymbol{\psi}^{\boldsymbol{\mu} + \varepsilon\_j}\rangle + L \sum\_{j=1,2} a\_j(t, t)n\_j |\boldsymbol{\psi}^{\boldsymbol{\mu} - \varepsilon\_j}\rangle \tag{26}$$

In this equation, we set (*ej*)*<sup>i</sup>* = *δij*. The exact NMQSD solution is simply the 'root' state of the hierarchy <sup>|</sup>*ψ* <sup>=</sup> <sup>|</sup>*ψ*(**0**). Truncating the hierarchy at a finite depth yields a closed set of equations, which can be numerically integrated. Typically, for a weak coupling, only a few auxiliary states have to be taken into account to achieve convergence. There also exists a hierarchy of operators complementary to HOPS, which is known as the hierarchy of equations of motions (HEOM) [25,30]. We define auxilliary operators with two multiindices (*n*, *m*) as follows.

$$\rho^{(\mathfrak{u},\mathfrak{m})} = E\left[|\Psi^{(\mathfrak{u})}\rangle\langle\Psi^{(\mathfrak{m})}|\right] \tag{27}$$

Upon employing the HOPS Equation (26) with relation *E*[*zt* ...] = *E* <sup>∑</sup>*<sup>j</sup>* <sup>D</sup>*<sup>j</sup>* ... , we find the hierarchical equations of motion:

$$\begin{split} \partial\_t \rho^{(\mathfrak{n}, \mathfrak{m})} &= -\mathbf{i} \left[ H\_\prime \rho^{(\mathfrak{n}, \mathfrak{m})} \right] - \sum\_{j=1,2} \left( n\_j \mathcal{W}\_j + m\_j \mathcal{W}\_j^\* \right) \rho^{(\mathfrak{n}, \mathfrak{m})} \\ &+ \sum\_{j=1,2} \left( n\_j a\_j(t, t) L \rho^{(\mathfrak{n} - \mathfrak{e}\_j, \mathfrak{m})} + m\_j a\_j^\*(t, t) \rho^{(\mathfrak{n}, \mathfrak{m} - \mathfrak{e}\_j)} L \right) \\ &+ \sum\_{j=1,2} \left( \left[ \rho^{(\mathfrak{n} + \mathfrak{e}\_j, \mathfrak{m})} , L \right] + \left[ L \, \rho^{(\mathfrak{n}, \mathfrak{m} + \mathfrak{e}\_j)} \right] \right), \end{split} \tag{28}$$

where the reduced state of the system is simply given by *ρ*(*t*) = *ρ*(**0**,**0**)(*t*). We can use a truncated hierarchy to accurately describe the non-Markovian open system dynamics of the model, even for strong coupling (large *γ*) and long memory time (small Γ).

#### **5. Dynamics**

The fidelity *Ft* <sup>=</sup> |*ψ*0|*ρ*(*t*)|*ψ*0|<sup>2</sup> of time-evolved state *<sup>ρ</sup>*(*t*) with an initial pure system state |*ψ*0 can be used to investigate how master equations capture short time dynamics. In general, the fidelity behaves for short times as [20]:

$$F\_l^F = 1 - \left(\Delta H\_S^2 + \Gamma\_F\right) t^2,\tag{29}$$

where <sup>Γ</sup>*<sup>F</sup>* <sup>=</sup> <sup>2</sup>*α*(0, 0) and <sup>Δ</sup>*X*<sup>2</sup> <sup>=</sup> *X*2−*X*<sup>2</sup> are the variances of operator *<sup>X</sup>*.

A quantum optical master equation that has constant coefficients, such as Equation (18), instead, leads to a linear short time behavior for the fidelity of the state of the system alone [21]:

$$F\_t^M = 1 - \Gamma\_M t + \mathcal{O}(t^2),\tag{30}$$

where Γ*<sup>M</sup>* = ∑*<sup>μ</sup> cμ*(*RμLμ*−*RμLμ*). Here, operators *Lμ*, *R<sup>μ</sup>* are the operators, which act on the left or on the right, respectively, in the sandwich term of the master equation.

The quadratic short time behavior is captured exactly by the Redfield theory, which is the lowest order perturbation theory.

$$F\_t^R = 1 - (\Delta H\_S^2 + \Gamma\_F)t^2\tag{31}$$

In Figure 1, we compare the short time linear rate Γ*<sup>M</sup>* computed for master Equation (18) and the quadratic rate Γ*<sup>F</sup>* for the full system. We observe that the full rate has a maximum at *ϕ* = *π*, which corresponds to squeezing along the *x*-direction. The two-level system is coupled to the *x*-quadrature of the bath, and squeezing in this direction increases the coupling strength which should lead to faster decay. In contrast, the master equation predicts a minimum for *φ* = *π*, which is not in line with the microscopic theory.

**Figure 1.** (**Left**): Linear rate Γ*<sup>M</sup>* for short time dynamics due to the Markov master equation divided by the overall coupling strength *γ*. (**Right**): Quadratic rate for short time dynamics for the total system Γ*<sup>F</sup>* divided by *γ*Γ, where Γ is the inverse bath correlation time scale. Γ*<sup>F</sup>* is proportional to |*u* − *v*| 2 , which has a maximum when *ϕ* = *π*. Γ*<sup>M</sup>* is proportional to *u*(*v* + *v*∗), which has a minimum at *ϕ* = *π*. Here, we consider the spin boson model (14) in a squeezed bath (9) with *u* = cosh(*r*) and *<sup>v</sup>* <sup>=</sup> sinh(*r*)*ei<sup>ϕ</sup>* and initial state <sup>|</sup>*ψ*0 <sup>=</sup> <sup>|</sup>+ :<sup>=</sup> <sup>√</sup><sup>1</sup> 2 (|0 + |1), where |0 and |1 are the eigenstates of *σz*. We have chosen the following parameters: *r* = 1/2, Γ = 5*γ*, and Ω = *ω<sup>B</sup>* = *ω*<sup>0</sup> = *γ*

In Figure 2, we compare the logarithm of the fidelity computed using the HEOM (28), the Redfield theory (16), and the master equation (18) with quadratic short time estimates (29), (31) and the linear estimate (30).

**Figure 2.** Short time dynamics of the spin boson model (14) in the squeezed bath (9) with Γ = 5*γ*, Ω = *ω<sup>B</sup>* = *ω*<sup>0</sup> = *γ*, *r* = 0.5, and |*ψ*0 = |+. The logarithm of the fidelity *Ft* is displayed for short times within different approximate and exact descriptions. The quadratic short time dynamics holds up to *γt* = 0.002 for *ϕ* = *π* as the HEOM and the Redfield theory start to deviate from the short time expansion. For *ϕ* = 0, the quadratic short time dynamics, HEOM, and the Redfield agree well in the range of the plot. This is explained by the effectively strong system environment coupling for *ϕ* = *π*, which sets a different regime of validity for the short time expansion. The parameters are chosen as in Figure 1.

The quadratic estimates; the HEOM and the Redfield curves are in line with each other for both squeezing directions *ϕ* = 0 and *ϕ* = *π*. The linear prediction and the master equation curves also are on top of each other but show a different behavior for different *φ*. For *ϕ* = 0, the decay of the fidelity is faster than for *ϕ* = *π*, as expected from Figure 1.

For longer times, the fidelity is still captured correctly by the Redfield equation, while the Markovian equation becomes valid only asymptoically, as seen in Figure 3.

**Figure 3.** Dynamics of the fidelity as in Figure 2 but for long times and for two squeezing directions *ϕ* = 0 (**top**) and for *ϕ* = *π* (**bottom**). The decay of the fidelity is faster for *ϕ* = *π* for HEOM and Redfield as expected from the microscopic model. Predictions from the quantum optical master equation (Markov) show an opposite behavior and the fidelity even has a positive slope at intermediate times.

We now turn to the long time dynamics, which reveals the steady state properties. In Figure 4, we compare the time evolution computed from HEOM, Redfield, and master equation for *ϕ* = 0. We do the same for *ϕ* = *π* in Figure 5. The same conclusions hold as before: For *ϕ* = *π*, the decay is faster, as can be observed from HEOM and Redfield curves. The steady state prediction of the master equation is far-off from the Refield and HEOM computations, which both agree exptremely well for the choosen parameters. This is because the memory time of the bath 1/Γ is chosen to be very short. As expected for a spin boson model, the steady state reaches a finite *σz* value. However, due to the non-stationary bath, some residual oscillations remain indefinitely in long time dynamics.

**Figure 4.** Long time dynamics of the spin boson model (14) in the squeezed bath (9) as in Figure 2 with *ϕ* = 0. Both approximate methods describe correctly that the *x* and *y* Bloch sphere components decay to zero. The *σz* expectation value asymptotically acquires a finite value modulated by weak oscillations which are captured properly in the Redfield theory.

**Figure 5.** Long time dynamics as in Figure 4, except *ϕ* = *π*.

Next, we investigate the Markov limit where we expect the master Equation (18) to be applicable. We chose Ω = *ω*<sup>0</sup> = 10*γ*, Γ = 3*γ*, which is in a regime where the assumptions for the rotating wave approximation are satisfied. Moreover, since Γ/*γ* = 3, the BCF decays faster than the time scale set for system bath interaction by *γ*. As we can observe from Figure 6, the master equation, the Redfield theory, and HEOM are in good qualitative agreement. The rapid oscillations present in HEOM and Redfield solutions are removed by the secular approximation behind the Markovian master equation.

**Figure 6.** Dynamics of the spin boson model (14) in the squeezed bath (9) with Γ = 3*γ*, Ω = 10*γ*, *ω*<sup>0</sup> = Ω, *r* = 0.5, *ϕ* = 0, and initial state |*ψ*0 = |+. The Markovian master Equation (18) agrees well with the exact dynamics but does not capture small asymptotic oscillations.

When we change the direction of the squeezing to *ϕ* = *π*, hence effectively coupling the system and the environment more strongly, the agreement is not as good, as shown in Figure 7. This in particularly visible for the expectation value *σy*. The asymptotic values are, however, captured accurately.

**Figure 7.** Dynamics of the spin boson model as in Figure 6, except *ϕ* = *π*.

Lastly, we investigate a regime for non-Markovian dynamics. The Redfield theory is a weak system–environment coupling theory. We expect it to fail when the bath correlation function decays on a slower time scale than the system environment interactions; that is, when Γ < *γ*. For short times, the Redfield theory agrees with the exact results. We confirm these statements in Figure 8, where one observes that the Redfield and the quantum optical master equation are both far-off from the exact solution. In fact, in this regime, the solution

of the Redfield equation ceases to be positive, as can be observed from *σz* < −1. To compute the exact dynamics, we have estimated convergence of HEOM by increasing the hierarchy depth until the relative change in the solution is less than 10−2.

**Figure 8.** Highly non-Markovian dynamics in the spin boson model. We have chosen |*ψ*0 = |+, Γ = *γ*/2, Ω = *ω<sup>B</sup>* = *ω*<sup>0</sup> = *γ*/2, and *ϕ* = 0. We are in the strong system–envinronment coupling regime since *γ* > Γ. This means that system–environment dynamics occurs on a faster time scale than the bath correlation function decay time. The failure of weak coupling master equations is expected.

#### **6. Conclusions**

In this article, we have revisited an old quantum optical model, namely, the decay to a squeezed reservoir [12]. This model has been realized experimentally using superconducting qubits [13] and is an interesting example of non-stationary open system dynamics. We investigated the short and long time dynamics of the model using textbook master equations [17–19], perturbative Redfield theory [32], and exact hierarchical methods [25,26], which we generalize to treat our non-Markovian squeezed bath model. Assuming broadband two mode squeezing, a model for the bath correlation function can be proposed where the bath memory time can be easily controlled by a single parameter Γ. This corresponds to a Lorentzian spectral density, which yields an exponentially decaying bath correlation function. We showed how in the limit where Γ *γ*, where *γ* is the overall system environment coupling, the Redfield theory performs extremely well. However, the dynamics is far from a regime where a usual "quantum optical master equation" is valid. This is due to the fact that there are other time scales involved which need to satisfy a strict hierarchical order so that the Markov approximation can be justified. These timescales are determined by the free system evolution (Ω) and the bath resonance (*ω*0). The quantum optical master equation is valid only in a regime where *ω*<sup>0</sup> is by far the fastest time-scale. All in all, the Redfield theory, which is based on weak coupling approximation only, produces the correct short and long time dynamics in a wide parameter range and is conceptually and numerically much simpler then more advanced methods such as HEOM. To those concerned about the possible positivity violations when using the Redfield theory, we would like to point out that this should be seen merely as a symptom of the breakdown of the weak coupling perturbation theory [32].

**Author Contributions:** V.L. and K.L. contributed equally to this work. W.T.S. contributed to the development of the concept and initial ideas behind this project. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

NMQSD Non-Markovian quantum state diffusion; HEOM Hierarchical equations Of motion; HOPS Hierarchy of pure states; BCF Bath correlation function.

#### **References**

