Non-Markovian dynamics in nonstationary Gaussian baths

TL;DR

This work addresses simulating open quantum systems coupled to nonstationary Gaussian baths, where the bath correlation function $\alpha(t,s)$ depends explicitly on both times. It generalizes the hierarchy of pure states (HOPS) by introducing a nonstationary bath decomposition $\alpha(t,s)=\sum_{j=1}^N \alpha_j(t-s) f_j(t) g_j^*(s)$ and representing the bath with a pseudo-Fock space, leading to a nonlinear stochastic Schrödinger equation with noise $Z(t)$ satisfying $\mathbb{E}[Z(t) Z^*(s)]=\alpha(t,s)$. The framework also yields a deterministic hierarchy of master equations (HME) and, for the special case $f_j=g_j$, a pseudomode master equation (PME) with Lindblad damping, along with a PSSE for stochastic sampling. Benchmarking on squeezed-bath models (single-mode and three-mode BCFs) shows rapid convergence with modest hierarchy depth and highlights regimes where PME or PSSE can outperform HOPS. These results enable efficient simulation of squeezed light–matter interactions and driven quantum materials in nonstationary environments.

Abstract

Building on the standard hierarchy of pure states (HOPS) approach, we construct a generalized formulation suitable for open quantum systems interacting with nonstationary Gaussian baths, potentially extending its applicability to nonequilibrium baths. This is achieved by extending the conventional exponential decomposition of bath correlation functions (BCF) to allow explicitly time-dependent forms. We demonstrate the method's performance on two examples of nonstationary squeezed reservoirs generated via uniform squeezing and degenerate parametric amplification. Benchmarking against the associated hierarchy of master equations shows that HOPS achieves superior efficiency under hierarchy truncation. In cases where each contribution in the BCF expansion can be associated with an independent physical bath, the formalism can be simplified in a pseudomode representation which is more efficient in a strongly non-Markovian regime. Our results highlight HOPS as a versatile and powerful tool for simulating open quantum systems in nonstationary baths, with potential applications ranging from squeezed light-matter interactions to driven quantum materials and dissipative phase transitions.

Figures (5)

• Figure 1: Expectation values of the Bloch vector components $\langle \hat{\sigma}_\alpha(t) \rangle$, where $\alpha=x,y,z$, computed for different values of $\Gamma$. The parameters are fixed as $\omega_0 = 5$, $\gamma = 1$, $r = 1.5$, and $\varphi = 0$. The atom is initially prepared in the pure state $\frac{|e\rangle + e^{-i\pi/4} |g\rangle}{\sqrt{2}}$. The dynamics is obtained by solving Eq. \ref{['eq: pseudomode master equation']}, using $100$ basis states for the pseudo-Fock space. The lower row shows the Bloch vector components in the rotating frame, defined as $\tilde{\sigma}_x(t) = \langle\hat{\sigma}_+(t) \rangle e^{-i\omega_0 t} + \text{c.c.}$ and $\tilde{\sigma}_y(t) = -i\langle\hat{\sigma}_+(t) \rangle e^{-i\omega_0 t} + \text{c.c.}$
• Figure 2: Comparison of the numerical error of the hierarchy of master equations (HME) \ref{['eq: HEOM']} and pseudomode master equation (PME) \ref{['eq: pseudomode master equation']} for different values of $\Gamma$. The other parameters are the same as in Fig. \ref{['fig: atomic dynamics']}: $\omega_0=5$, $\gamma=1$, $r=1.5$, and $\varphi = 0$. The total error includes contributions from the integration method and from hierarchy truncation [see Appendix \ref{['sec: error ME']}]. The upper row (a)-(d) shows the time dependence of the error for the hierarchy of master equations. The corresponding plot for the pseudomode master equation is not shown, as it exhibits the same qualitative behavior. The lower row (e)-(h) shows the root-mean-square error. For both approaches, the reference density operator is obtained by solving the corresponding equation (Eq.\ref{['eq: HEOM']} or Eq.\ref{['eq: pseudomode master equation']}) with a hierarchy depth $n^\text{max} = 100$.
• Figure 3: Estimated numerical error for the HOPS method and the pseudomode stochastic Schrödinger equation (PSSE). The other parameters are the same as in Figs. \ref{['fig: atomic dynamics']} and \ref{['fig: ME truncation']}. Statistical averages were obtained using $10^5$ trajectories. For the PSSE case with $\Gamma = 0.2$ shown in panel (e), a time step of $T \times 10^{-5}$ was used; in all other cases, the time step was $T \times 10^{-4}$. The upper row (a)-(d) shows the time dependence of the error for HOPS. The corresponding plot for the PSSE is not shown, as it exhibits the same qualitative behavior. The lower row (e)–(h) shows the time-averaged Euclidean norm of the error across different methods. These panels also include the mean errors for HME and PME from the corresponding panels in Fig. \ref{['fig: ME truncation']} for direct comparison. The reference density matrix was the same as in Fig. \ref{['fig: ME truncation']}.
• Figure 4: (a) Squeezed light generation in a one-sided cavity with a degenerate parametric amplifier. A nonlinear crystal inside the cavity is driven at frequency $2\omega_0$ with pump amplitude $\epsilon$, and the down-conversion process is described by $\hat{H}_\text{DPA}$. After reflection from the perfect mirror, the intracavity field exits through the transmitting mirror and interacts with a two-level atom, located far from the cavity. The vacuum input is modeled with a Lorentzian spectrum centered at $\omega_0$ and width $2\Gamma_0$. (b)–(e) Mean Bloch vector dynamics of the atom under this driving, characterized by the BCF in Eq. \ref{['eq: Gardiner BCF']}. Parameters: $\omega_0 = 5$, $\Gamma_0=2$, $\Gamma = 1$, $\epsilon = 0.5$, $\gamma = 1$ and $\varphi = \pi$. The atom is initially in the state $\frac{|e\rangle + e^{-i\pi/4} |g\rangle}{\sqrt{2}}$.
• Figure 5: Estimated numerical error of the HOPS method. Parameters are the same as in Fig. \ref{['fig: atomic dynamics Gardiner']}. Statistical averages were obtained from $10^5$ trajectories with a time step of $20\times10^{-4}$. The upper row (a)–(d) shows the time dependence of the error: (a) $n^{\text{max}}_2=n^{\text{max}}_3=20$, varying $n^{\text{max}}_1$; (b) $n^{\text{max}}_1=n^{\text{max}}_3=20$, varying $n^{\text{max}}_2$; (c) $n^{\text{max}}_1=n^{\text{max}}_2=20$, varying $n^{\text{max}}_3$; (d) triangular truncation with $n_1+n_2+n_3\leq n^{\text{sum}}$, varying $n^{\text{sum}}$. The lower row (e)–(h) shows the corresponding time-averaged errors. The reference density matrix is obtained using the same time step, averaging over $10^6$ trajectories and truncation levels $n^{\text{max}}_1=n^{\text{max}}_2=n^{\text{max}}_3=20$.