A Low-Rank Hierarchical Framework for the Non-Markovian Stochastic Schrödinger Equation with Convergence Analysis 3
July 16, 2026
Abstract
We propose and analyze a novel numerical framework for the non-Markovian stochastic Schrödinger equation (NMSSE) based on a low-rank approximation of the bath correlation functions. By decomposing the memory kernel into a finite-dimensional representation, we derive a truncated system of hierarchical equations that effectively balances computational tractability with physical fidelity. A rigorous convergence analysis is established for the hierarchical framework under mild assumptions. We demonstrate that our formulation serves as a mathematical generalization of the Hierarchy of Pure States (HOPS), encompassing it as a special case while offering a more flexible representation of non-Markovian effects. Numerical experiments across several benchmark models are presented to illustrate the validity and efficacy of the proposed method.
Non-Markovian dynamics, Stochastic Schrödinger equation, Low-rank approximation, Convergence analysis, Hierarchical equations.
65C30, 81S22, 60H35, 65L03.
1 Introduction↩︎
Open quantum systems (OQS) play a central role in quantum physics and chemistry, with applications in quantum information processing [1], quantum optics [2], and chemical physics [3]. A key difficulty in their simulation is the non-Markovian memory encoded in the bath correlation function (BCF). Depending on how the reduced dynamics is represented, existing numerical methods are often organized into two broad classes.
The first class evolves the reduced density matrix (RDM). Representative approaches include path integrals [4], such as QuAPI [5] and the Dyson series iterative method [6]. To reduce the computational burden, i-QuAPI [7], [8] truncates memory, and [9] combines low-rank approximation with the frozen Gaussian approximation. Another line of work starts from the Nakajima-Zwanzig equation [10], [11], including TTM [12] and HEOM [13], [14]. Despite this progress, RDM-based methods still face a severe curse of dimensionality for large systems or complex baths.
The second class is stochastic unravelling [15], which represents the reduced state by an ensemble of stochastic trajectories. In the Markovian regime, the stochastic Schrödinger equation [2], [16] derived from the Lindblad master equation [17] can be solved by standard SDE methods [18], [19]. For general non-Markovian dynamics, however, the NMSSE [20], [21] is difficult to solve directly. Existing strategies include time-local formulations based on the \(O\)-operator [21], [22] and the hierarchy of pure states (HOPS) [23], but these methods are still constrained by the accuracy of the \(O\)-operator approximation or by the exponential decomposition of the BCF. A rigorous convergence analysis for hierarchical stochastic solvers also remains limited.
Motivated by these limitations, we propose a low-rank hierarchical method for the NMSSE. The method uses a low-rank approximation of the BCF to compress environmental memory, leading to a new infinite-dimensional hierarchy that generalizes HOPS. We then introduce a finite-order truncation and prove its convergence under the assumption that the NMSSE solution is unique. Numerical experiments on benchmark OQS models confirm the efficiency and accuracy of the proposed algorithm.
The remainder of this paper is organized as follows. Section 2 introduces the NMSSE. Section 3 develops the low-rank approximation induced hierarchical framework and its finite truncation. Section 4 establishes the convergence of the hierarchical framework, including the low-rank approximation and the finite truncation. Section 5 presents numerical schemes for the finite hierarchical system. Section 6 reports numerical experiments, and Section 7 concludes the paper.
2 Preliminaries↩︎
Let \(\mathcal{H}\) be a separable Hilbert space with inner product \(\langle\cdot,\cdot\rangle_{\mathcal{H}}\) and induced norm \(\|\cdot\|_{\mathcal{H}}\). Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space. We consider the following non-Markovian stochastic Schrödinger equation [20], [21]: \[\begin{align} \label{eq:NMSSE} \frac{\partial}{\partial t} \psi_t = \left(-iH+z_tL\right)\psi_t - L^{\dagger}\int_{0}^{t}\mathrm{d}s\,\alpha(t,s)\frac{\delta}{\delta z_s}\psi_t, \quad t\in[0, T]. \end{align}\tag{1}\] The mathematical constituents are defined as follows:
\(\alpha(\cdot,\cdot)\in \mathcal{C}([0,T]^2,\mathbb{C})\) is the two-time bath correlation function, which characterizes temporal correlations of the bath. It is Hermitian and positive semidefinite; i.e., for any \(f\in L^2([0,T],\mathbb{C})\), \[\int_0^T\int_0^T f^*(u)\alpha(u,v)f(v)\mathrm{d}u\mathrm{d}v \geq 0.\]
\(z\in L^2(\Omega; \mathcal{C}([0,T],\mathbb{C}))\) is a complex-valued centered Gaussian process over the probability space \((\Omega, \mathcal{F}, \mathbb{P})\) with covariance \(\alpha\) and pseudo-covariance \(0\), i.e., \[\mathbb{E}(z_t) = 0, \quad \mathbb{E}(z_tz_s) = 0,\quad \mathbb{E}(z_t^*z_s) = \alpha(t,s),\quad \forall t,s\in [0, T]. \label{eq:noise95condition}\tag{2}\] Let \(\mathcal{F}_t^z\) be the \(\sigma\)-algebra generated by \(z\), i.e., \(\mathcal{F}_t^z = \sigma(z_s: s\in[0,t])\).
\(H:D(H)\subset\mathcal{H}\to\mathcal{H}\), which represents the system Hamiltonian, is a self-adjoint linear operator such that \(-iH\) serves as the infinitesimal generator of a unitary group \(S(\cdot): S(t)=e^{-itH}\).
\(L\in\mathcal{L}(\mathcal{H})\) is the system-bath coupling operator, which is bounded with respect to \(\|\cdot\|_{\mathcal{H}}\), i.e., \[\|L\|_{\mathcal{H}} \mathrel{\vcenter{:}}= \sup_{\|\psi\|_{\mathcal{H}}= 1}\|L\psi\|_{\mathcal{H}} < \infty.\]
The initial value \(\psi_0\) is an \(\mathcal{H}\)-valued random variable independent of \(\mathcal{F}_T^z\). We consider the solution \(\psi_t\) adapted to the filtration \(\mathcal{F}_t^{z,\psi_0}\) generated by \(\psi_0\) and \(z_s, s\in[0,t]\).
For the last term on the right-hand side of 1 , we denote \[\mathcal{R}_t\psi_t \coloneq \int_{0}^{t} \mathrm{d}s\, \alpha(t,s)\frac{\delta}{\delta z_s}\psi_t,\] where \(\mathcal{R}_t : D(\mathcal{R}_t)\subset\mathcal{Y}_t\to\mathcal{H}\) is a linear operator, and \(\mathcal{Y}_t\) denotes the set of all \(\mathcal{F}_t^{z,\psi_0}\)-measurable \(\mathcal{H}\)-valued random variables over \((\Omega, \mathcal{F}, \mathbb{P})\). The term \(\mathcal{R}_t\psi_t\) is interpreted in the sense of the Gâteaux derivative, i.e., for any \(\omega\in\Omega\), \[\mathcal{R}_t\psi_t(\omega) = \left.\frac{\mathrm{d}}{\mathrm{d}\varepsilon}\mathscr{M}_t\left(\psi_0(\omega),z(\omega)|_{[0,t]} + \varepsilon\alpha(t,\cdot)|_{[0,t]}\right)\right|_{\varepsilon=0},\] where \(\mathscr{M}_t: \mathcal{H}\times \mathcal{C}([0,t],\mathbb{C}) \to \mathcal{H}\) is a mapping such that \[\psi_t(\omega) = \mathscr{M}_t\left(\psi_0(\omega),z(\omega)|_{[0,t]}\right)\qquad \forall\, \omega\in \Omega.\] The existence of \(\mathscr{M}_t\) is guaranteed by the Doob-Dynkin lemma [24].
An \(\mathcal{F}_t^{z,\psi_0}\)-adapted \(\mathcal{H}\)-valued process \(\psi_t\), \(t\in[0,T]\), is a mild solution [25] to 1 if \(\psi_t\) takes values in \(D(\mathcal{R}_t)\), a.s., \[\begin{align} \tag{3} & \mathbb{P}\left(\int_{0}^{T}\left(1+|z_t|\|L\|_{\mathcal{H}}\right)\|\psi_t\|_{\mathcal{H}}\mathrm{d}t < \infty\right) = 1, \\ \tag{4} & \mathbb{P}\left(\int_{0}^{T}\left\|\mathcal{R}_t\psi_t\right\|_{\mathcal{H}}\mathrm{d}t < \infty\right) = 1, \end{align}\] and for arbitrary \(t\in[0,T]\), it holds a.s. that \[\label{eq:mild95formulation} \psi_t = S(t)\psi_0 + \int_0^t S(t-s) \left(z_sL - L^{\dagger}\mathcal{R}_s\right)\psi_s \mathrm{d}s.\tag{5}\]
Let \(\widetilde{\psi}_t = S(-t)\psi_t\) and \(\widetilde{L}(s) = S(-s)L S(s)\) denote the interaction picture of the mild solution \(\psi_t\) and the operator \(L\), respectively. Then it holds that \[\label{eq:strong95solution} \widetilde{\psi}_t = \widetilde{\psi}_0 + \int_0^t \left(z_s\widetilde{L}(s) - \widetilde{L}^{\dagger}(s)\mathcal{R}_s\right)\widetilde{\psi}_s \mathrm{d}s,\tag{6}\] which follows from the fact that \(\mathcal{R}_s\) commutes with \(S(\cdot)\). This identity indicates that the mild solution of 1 is equivalent to the strong solution of the NMSSE in the interaction picture: \[\label{eq:NMSSE95interaction} \frac{\partial}{\partial t}\widetilde{\psi}_t = \left(z_t\widetilde{L}(t) - \widetilde{L}^{\dagger}(t)\mathcal{R}_t\right)\widetilde{\psi}_t.\tag{7}\] In what follows, we shall exploit this equivalence to analyze the strong solution of the interaction picture equation whenever it is notationally or analytically advantageous.
The main numerical challenge in the NMSSE is the functional-derivative integral together with colored noise. This feature prevents a direct application of standard solvers for ordinary differential equations and stochastic differential equations, and makes direct simulation of non-Markovian dynamics computationally difficult. To address this issue, existing approaches usually approximate or reformulate the functional derivative, thereby replacing the original temporally nonlocal problem with a tractable evolution system. In many cases, this corresponds to embedding the non-Markovian dynamics into an enlarged time-local, or effectively Markovian, system.
However, existing methods often impose restrictive assumptions on the bath correlation function or lack rigorous convergence guarantees. Motivated by these limitations, we develop a generalized hierarchical framework based on a low-rank approximation of the BCF, and we provide a rigorous convergence analysis to establish its theoretical validity.
3 A low-rank hierarchical framework↩︎
3.1 Infinite hierarchy↩︎
We begin by approximating the bath correlation function through a low-rank approximation. Specifically, \[\alpha(t,s) \approx \sum_{j=1}^{r} \lambda_j V_j^*(t)V_j(s),\qquad \forall t,s\in[0,T], \label{eq:low95rank95decomp95corr95func}\tag{8}\] where \(\lambda_1\geq\cdots\geq\lambda_r\geq0\) and \(V_j\in L^2([0,T],\mathbb{C})\) for \(j=1,\cdots,r\). The existence of such an approximation follows from the positive semidefiniteness of the BCF and Mercer’s theorem [26]–[29]. In practice, we compute this representation via singular value decomposition (SVD) or eigenvalue decomposition of the discretized bath correlation matrix \(A = \{\alpha(t_i,t_j)\}\), and retain the leading \(r\) components.
Substituting 8 into 1 , we obtain \[\frac{\partial}{\partial t} \psi_t = \left(-iH+z_tL\right)\psi_t - L^{\dagger}\sum_{j=1}^{r}\lambda_jV_j^*(t)\int_{0}^{t}\mathrm{d}sV_j(s)\frac{\delta}{\delta z_s}\psi_t. \label{eq:low95rank95decomp95NMSSE}\tag{9}\] To derive a hierarchy following the HOPS formalism, we introduce the operators \[D_j = \int_{0}^{\infty} \mathrm{d}sV_j(s)\frac{\delta}{\delta z_s}, \qquad j = 1,2,\cdots,r. \label{eq:operator95Dj}\tag{10}\] Although the upper limit of integration in 10 is extended to \(\infty\) in contrast to \(t\) in 9 , the following identity holds: \[D_j \psi_t = \int_{0}^{\infty} \mathrm{d}sV_j(s)\frac{\delta}{\delta z_s}\psi_t = \int_{0}^{t} \mathrm{d}sV_j(s)\frac{\delta}{\delta z_s}\psi_t.\] This equivalence follows from the fact that \(\psi_t\) is \(\mathcal{F}_t^{z,\psi_0}\)-adapted. Next, let \(\boldsymbol{k}\in\mathbb{N}^r\) be a multi-index whose magnitude is defined by \(|\boldsymbol{k}|\mathrel{\vcenter{:}}=\sum_{j=1}^rk_j\). We define the auxiliary wavefunctions recursively as \[\psi_t^{(\boldsymbol{0})}\mathrel{\vcenter{:}}= \psi_t, \qquad \psi^{(\boldsymbol{e}_j)}_t\mathrel{\vcenter{:}}= D_j\psi_t^{(\boldsymbol{0})}, \qquad \psi_t^{(\boldsymbol{k}+\boldsymbol{e}_j)} \mathrel{\vcenter{:}}= D_j\psi_t^{(\boldsymbol{k})},\] where \(\boldsymbol{e}_j\) denotes the unit vector with a \(1\) at the \(j\)th position and \(0\) elsewhere. Consequently, the evolution equation for the primary wavefunction \(\psi_t^{(\boldsymbol{0})}\) is given by \[\frac{\partial}{\partial t}\psi_t^{(\boldsymbol{0})} = (-iH + z_tL)\psi_t^{(\boldsymbol{0})} - L^{\dagger}\sum_{j=1}^{r}\lambda_jV_j^*(t)\psi_t^{(\boldsymbol{e}_j)}.\]
We now derive the evolution equation for \(\psi_t^{(\boldsymbol{e}_j)}\). Differentiating in time and applying the chain rule gives \[\begin{align} \tag{11} \frac{\partial}{\partial t}\psi_t^{(\boldsymbol{e}_j)} &= \frac{\partial}{\partial t}\left(D_j\psi_t^{(\boldsymbol{0})}\right) \\ &= \frac{\partial D_j}{\partial t}\psi_t^{(\boldsymbol{0})} + D_j\frac{\partial\psi_t^{(\boldsymbol{0})}}{\partial t} \\ &= D_j\left[\left(-iH+z_tL\right)\psi_t^{(\boldsymbol{0})} - L^{\dagger}\sum_{k=1}^{r}\lambda_kV_k^*(t)\psi_t^{(\boldsymbol{e}_k)}\right] \\ &= \left(-iH+z_tL\right)\psi_t^{(\boldsymbol{e}_j)} + V_j(t)L\psi_t^{(\boldsymbol{0})} - L^{\dagger}\sum_{k=1}^{r}\lambda_kV^*_k(t)\psi_t^{(\boldsymbol{e_j}+\boldsymbol{e}_k)}, \tag{12} \end{align}\] In 12 , we have utilized the fact that \(D_j\) commutes with both \(H\) and \(L\), along with the identity \(D_jz_t = V_j(t)\). This relation demonstrates that the evolution of \(\psi_t^{(\boldsymbol{e}_j)}\) depends on the higher-order auxiliary wavefunctions \(\psi_t^{(\boldsymbol{e}_j+\boldsymbol{e}_k)}\). By induction, we can present the coupled hierarchy for an arbitrary multi-index \(\boldsymbol{k}\in\mathbb{N}^r\): \[\frac{\partial}{\partial t} \psi_t^{(\boldsymbol{k})} = (-iH+z_tL)\psi_t^{(\boldsymbol{k})} + L\sum_{j=1}^{r}k_jV_j(t)\psi_t^{(\boldsymbol{k}-\boldsymbol{e}_j)} - L^{\dagger}\sum_{j=1}^{r}\lambda_jV_j^*(t)\psi_t^{(\boldsymbol{k}+\boldsymbol{e}_j)},\] where \(\psi_t^{(\boldsymbol{k})}\equiv 0\) for \(\boldsymbol{k}\notin\mathbb{N}^r\).
3.2 Finite truncation↩︎
To obtain a computationally tractable model, we truncate the infinite hierarchy at a finite order. In contrast to the various truncation schemes developed for HOPS [23], we adopt a straightforward zero-truncation strategy; a detailed convergence analysis of this approach is provided in the subsequent section. Let \(N\) denote the truncation order. The state of the system is then described by the set of auxiliary wavefunctions \[\{\psi_t^{(\boldsymbol{k})}:\boldsymbol{k}\in\mathcal{K}^r_N\},\] where \(\mathcal{K}^r_N\mathrel{\vcenter{:}}=\{\boldsymbol{k}\in\mathbb{N}^r:|\boldsymbol{k}|\leq N\}\). By imposing \(\psi_t^{(\boldsymbol{k})}\equiv 0\) for \(\boldsymbol{k}\notin\mathcal{K}^r_N\), we obtain the closed system \[\frac{\partial}{\partial t} \psi_t^{(\boldsymbol{k})} = (-iH+z_tL)\psi_t^{(\boldsymbol{k})} + L\sum_{\substack{j=1 \\ \boldsymbol{k}-\boldsymbol{e}_j\in\mathcal{K}^r_N}}^{r}k_jV_j(t)\psi_t^{(\boldsymbol{k}-\boldsymbol{e}_j)} - L^{\dagger}\sum_{\substack{j=1\\ \boldsymbol{k}+\boldsymbol{e}_j\in\mathcal{K}^r_N}}^{r}\lambda_jV_j^*(t)\psi_t^{(\boldsymbol{k}+\boldsymbol{e}_j)}, \label{eq:finite95hierarchy}\tag{13}\] for \(\boldsymbol{k}\in\mathcal{K}^r_N\), with the initial condition \[\psi_0^{(\boldsymbol{0})} =\psi_0,\quad \text{and} \quad \psi_0^{(\boldsymbol{k})}=0,\quad \forall\boldsymbol{k}\in\mathcal{K}^r_N, |\boldsymbol{k}|>0.\]
4 Convergence of the hierarchical framework↩︎
In this section, we establish convergence of the proposed hierarchical framework via a two-stage argument. We first introduce the following uniqueness assumption:
- For every initial condition \(\psi_0\in\mathcal{H}\), the NMSSE 1 admits at most one mild solution on \([0,T]\). More precisely, if \(\psi\) and \(\phi\) are two \(\mathcal{F}_t^{z,\psi_0}\)-adapted \(\mathcal{H}\)-valued processes with continuous sample paths satisfying 3 , 4 and 5 with the same initial condition and driven by the same noise, then \[\mathbb{P}\left(\sup_{t\in[0,T]}\left\|\psi_t-\phi_t\right\|_{\mathcal{H}}=0\right) = 1.\]
Existence of the mild solution will be established via an infinite-series construction in the following analysis. The main result is stated in the next theorem:
Theorem 1. Let \(\sum_{j=1}^{r}\lambda_jV_j^*(t)V_j(s)\) be a rank-\(r\) approximation of \(\alpha(t,s)\). Under Assumption [assum:exist95uniq], there exist \(r_0>0\) such that for all \(r>r_0\) and \(N>0\), \[\mathbb{E}\left(\sup_{t\in[0,T]}\left\|\psi_t - \psi_t^{r,N}\right\|_{\mathcal{H}}\right) \leq C_1\varepsilon_r + C_2\sqrt{\frac{(\beta r)^N}{N!}},\] where \(C_1, C_2, \beta>0\) are generic constants which only depend on \(\alpha\) and \(T\), \(\psi_t\) is the mild solution of 1 , \(\psi_t^{r,N}\) denotes the mild solution \(\psi_t^{(\boldsymbol{0})}\) of 13 with rank \(r\) and truncation order \(N\), and \[\varepsilon_r \coloneq \sqrt{\sum_{j=r+1}^{\infty}\lambda_j^2} = o\left(r^{-\frac{1}{2}}\right). \label{eq:epsilon95r}\qquad{(1)}\]
The proof of 1 is based on the decomposition \[\mathbb{E}\left(\left\|\psi_t - \psi_t^{r,N}\right\|_{\mathcal{H}}\right) \leq \mathbb{E}\left(\left\|\psi_t - \psi_t^{r}\right\|_{\mathcal{H}}\right) + \mathbb{E}\left(\left\|\psi_t^{r} - \psi_t^{r,N}\right\|_{\mathcal{H}}\right),\] where \(\psi_t^{r}\) denotes the mild solution of 9 . In the first part, we estimate the error between the exact NMSSE solution and the solution associated with the low-rank BCF approximation. Under Assumption [assum:exist95uniq], we show that this perturbation vanishes as the rank increases. In the second part, we analyze the truncation error induced by retaining only finitely many hierarchy levels, and prove convergence to the infinite-hierarchy solution as the truncation order tends to infinity. Combining the two estimates yields convergence of the proposed framework to the target dynamics. The specific details of these derivations are provided in the following two subsections.
4.1 Convergence of the low-rank approximation↩︎
To establish the convergence of the solution to the NMSSE under the low-rank BCF approximation, we first derive an infinite-series representation of the mild solution.
Theorem 2. The mild solution of 1 is given by the infinite series \[\begin{align} \psi_t =& S(t)\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\widetilde{\psi}_t^{(n,m)} \\ \coloneq& S(t)\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\int_{0<\tau_1<\cdots<\tau_n<t}\mathrm{d}^{n}\boldsymbol{\tau}z_{\tau_n}\cdots z_{\tau_1} \\ &\qquad\quad \int_{0<s_1<\cdots<s_{2m}<t}\mathrm{d}^{2m}\boldsymbol{s} (-1)^m\sum_{\boldsymbol{q}\in\mathscr{Q}(\boldsymbol{s})}\alpha(q_{2m},q_{2m-1})\cdots\alpha(q_2,q_1) \\ &\qquad\quad \mathscr{T}\left\{\widetilde{L}(\tau_n)\cdots\widetilde{L}(\tau_1)\widetilde{L}^{\dagger}(q_{2m})\widetilde{L}(q_{2m-1})\cdots\widetilde{L}^{\dagger}(q_2)\widetilde{L}(q_1)\right\}\psi_0, \label{eq:infinite95series95solution} \end{align}\qquad{(2)}\] where \(\mathscr{T}\) is the time-ordering operator, \(\widetilde{L}(\tau)\mathrel{\vcenter{:}}= S(-\tau)LS(\tau)\) is the interaction picture of \(L\), and \(\mathscr{Q}(\boldsymbol{s})=\mathscr{Q}(s_{1},\cdots,s_{2m})\) denotes the set of all possible ordered pairings of \(\{s_1,\cdots,s_{2m}\}\). For example, \[\begin{align} \mathscr{Q}(s_1,s_2) =& \big\{\{(s_1,s_2)\}\big\}, \\ \mathscr{Q}(s_1,s_2,s_3,s_4) =& \big\{\{(s_1,s_2),(s_3,s_4)\}, \{(s_1,s_3),(s_2,s_4)\}, \{(s_1,s_4),(s_2,s_3)\}\big\}. \end{align}\] The first element of each pair is less than the second. If \(\boldsymbol{q}\) is equal to some set of ordered pairings, it indicates that the elements of \(\boldsymbol{q}\) are, in order, equal to all the elements in that set. For example, \(\boldsymbol{q}=\{(s_1,s_4),(s_2,s_3)\}\) implies \(q_1=s_1, q_2=s_4, q_3=s_2, q_4=s_3\).
The infinite series ?? absolutely converges a.s. and has an upper bound \[\begin{align} \left\|\psi_t\right\|_{\mathcal{H}}\leq& \mathcal{M}_t\|\psi_0\|_{\mathcal{H}}\\ \coloneq& \exp\left\{\|L\|_{\mathcal{H}}\int_{0}^t|z_s|\mathrm{d}s + \frac{\|L\|^2_{\mathcal{H}}}{2}\int_0^t\int_0^t|\alpha(u,v)|\mathrm{d}u\mathrm{d}v\right\}\|\psi_0\|_{\mathcal{H}}. \label{eq:upper95bound95inf95series} \end{align}\qquad{(3)}\]
Proof. The verification that ?? satisfies 3 , 4 and 5 is deferred to the supplementary materials. Here we prove absolute convergence of the series in ?? , i.e., \[\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\left\|\widetilde{\psi}_t^{(n,m)}\right\|_{\mathcal{H}} < \infty\] for \(t\in[0,T]\) a.s.. Note that \(\|\widetilde{L}(t)\|_{\mathcal{H}}=\|S(-t)LS(t)\|_{\mathcal{H}}=\|L\|_{\mathcal{H}}\). Hence, \[\begin{align} \left\|\widetilde{\psi}_t^{(n,m)}\right\|_{\mathcal{H}} \leq& \int_{0<\tau_1<\cdots<\tau_n<t}\mathrm{d}^n\boldsymbol{\tau}\left|z_{\tau_n}\cdots z_{\tau_1}\right|\int_{0<s_1<\cdots<s_{2m}<t}\mathrm{d}^{2m}\boldsymbol{s} \\ & \sum_{\boldsymbol{q}\in\mathscr{Q}(\boldsymbol{s})}\left|\alpha(q_{2m},q_{2m-1})\cdots\alpha(q_2,q_1)\right|\cdot \|L\|_{\mathcal{H}}^{n+2m} \cdot\|\psi_0\|_{\mathcal{H}}. \label{eq:absolute95converg95estimate} \end{align}\tag{14}\] We now rewrite the two integrals over the cubes \([0,t]^n\) and \([0,t]^{2m}\). By symmetry, \[\int_{0<\tau_1<\cdots<\tau_n<t}\mathrm{d}^n \boldsymbol{\tau}\left|z_{\tau_n}\cdots z_{\tau_1}\right|\cdot\|L\|_{\mathcal{H}}^n = \frac{1}{n!}\left(\|L\|_{\mathcal{H}}\int_0^t|z_s|\mathrm{d}s\right)^n. \label{eq:4958}\tag{15}\] Let \(f(\boldsymbol{s})\mathrel{\vcenter{:}}=\sum_{\boldsymbol{q}\in\mathscr{Q}(\boldsymbol{s})}|\alpha(q_{2m},q_{2m-1})\cdots\alpha(q_2,q_1)|\). Since \(\alpha(t,s)=\alpha^*(s,t)\), we have \(|\alpha(t,s)|=|\alpha(s,t)|\), and thus \(f(\boldsymbol{s})\) is symmetric under permutations of \(s_1,\cdots,s_{2m}\). Because \(\mathscr{Q}(s_1,\cdots,s_{2m})\) contains \((2m-1)!!\) elements, we obtain \[\begin{align} \int_{0<s_1<\cdots<s_{2m}<t}\mathrm{d}^{2m}\boldsymbol{s}f(\boldsymbol{s}) \|L\|_{\mathcal{H}}^{2m} =& \frac{1}{(2m)!}\int_{[0,t]^{2m}}\mathrm{d}^{2m}\boldsymbol{s}f(\boldsymbol{s})\|L\|_{\mathcal{H}}^{2m} \\ =& \frac{(2m-1)!!}{(2m)!}\left(\int_0^t\int_0^t\left|\alpha(u,v)\right|\mathrm{d}u\mathrm{d}v\right)^m\|L\|_{\mathcal{H}}^{2m} \\ =& \frac{1}{m!}\left(\frac{\|L\|_{\mathcal{H}}^2}{2}\int_0^t\int_0^t|\alpha(u,v)|\mathrm{d}u\mathrm{d}v\right)^{m}. \label{eq:4959} \end{align}\tag{16}\] Substituting 15 and 16 into 14 yields \[\left\|\widetilde{\psi}_t^{(n,m)}\right\|_{\mathcal{H}} \leq \frac{1}{n!}\left(\|L\|_{\mathcal{H}}\int_0^t|z_s|\mathrm{d}s\right)^n \frac{1}{m!}\left(\frac{\|L\|_{\mathcal{H}}^2}{2}\int_0^t\int_0^t|\alpha(u,v)|\mathrm{d}u\mathrm{d}v\right)^m \|\psi_0\|_{\mathcal{H}}. \label{eq:49510}\tag{17}\] Summing over \(n,m\geq 0\) yields \[\begin{align} & \sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\left\|\widetilde{\psi}_t^{(n,m)}\right\|_{\mathcal{H}} \\ &\qquad\leq \exp\left\{\|L\|_{\mathcal{H}}\int_{0}^t|z_s|\mathrm{d}s + \frac{\|L\|^2_{\mathcal{H}}}{2}\int_0^t\int_0^t|\alpha(u,v)|\mathrm{d}u\mathrm{d}v\right\}\|\psi_0\|_{\mathcal{H}} \\ &\qquad< \infty,\quad \text{a.s.}, \end{align}\] which proves ?? . ◻
The following lemma provides an essential estimate for \(\mathbb{E}(\exp\{\int_0^t|z_s|\mathrm{d}s\})\), which will be utilized in the subsequent analysis.
Lemma 1. Let \(z\) be the Gaussian process satisfying 2 . It holds that \[\mathbb{E}\left(\exp\left\{\int_0^t|z_s|\mathrm{d}s\right\}\right) \leq \int_0^t \sqrt{\pi\alpha(s,s)}\exp\left\{\frac{\alpha(s,s)t^2}{4}\right\}\mathrm{d}s + 1.\]
Proof. By the convexity of the exponential function, Jensen’s inequality implies \[\exp\left\{\int_0^t|z_s|\mathrm{d}s\right\} = \exp\left\{\frac{1}{t}\int_0^tt|z_s|\mathrm{d}s\right\} \leq \frac{1}{t}\int_0^te^{t|z_s|}\mathrm{d}s.\] Taking expectation gives \[\mathbb{E}\left(\exp\left\{\int_0^t|z_s|\mathrm{d}s\right\}\right) \leq \frac{1}{t}\mathbb{E}\left(\int_0^te^{t|z_s|}\mathrm{d}s\right) = \frac{1}{t}\int_0^t\mathbb{E}\left(e^{t|z_s|}\right)\mathrm{d}s, \label{eq:estimate95414}\tag{18}\] where exchanging expectation and integration follows from Tonelli’s theorem. Note that \(|z_s|\) follows a Rayleigh distribution with density \[f(r) = \frac{2r}{\alpha(s,s)}e^{-r^2/\alpha(s,s)},\qquad r\geq 0.\] A direct calculation gives \[\begin{align} \label{eq:estimate95416} \mathbb{E}\left(e^{t|z_s|}\right) =& 1+\sqrt{\pi\alpha(s,s)}te^{\alpha(s,s)t^2/4}\Phi\left(\sqrt{\frac{\alpha(s,s)}{2}}t\right) \\ \leq& 1+\sqrt{\pi\alpha(s,s)}te^{\alpha(s,s)t^2/4}, \end{align}\tag{19}\] where \(\Phi\) is the cumulative distribution function of the standard normal distribution. Substituting 19 into 18 completes the proof. ◻
Let \(\widehat{\alpha}_r(t,s)=\sum_{j=1}^r\lambda_jV_j^*(t)V_j(s)\) be the low-rank approximation of \(\alpha(t,s)\) of rank \(r\). The uniform convergence is given by the following lemma (see [27]–[29]).
Lemma 2 (Mercer’s theorem). Let \(\alpha(\cdot,\cdot)\in \mathcal{C}([0,T]^2,\mathbb{C})\) be positive semi-definite. Define the integral operator \(\mathcal{I}:L^2([0,T],\mathbb{C})\to L^2([0,T],\mathbb{C})\) as \[(\mathcal{I}f)(t) \mathrel{\vcenter{:}}= \int_0^T \alpha^*(t,s)f(s)\mathrm{d}s.\] Let \(V_j\in L^2([0,T],\mathbb{C})\) be the orthonormal eigenfunctions of \(\mathcal{I}\) associated with the eigenvalues \(\lambda_j\geq 0\), sorted in non-increasing order. Then
the eigenvalues \(\{\lambda_j\}_{j=1}^{\infty}\) are absolutely summable;
\(\alpha(t,s)=\sum_{j=1}^{\infty} \lambda_jV_j^*(t)V_j(s)\) holds for almost every \((t,s)\in[0,T]^2\), where the series converges absolutely and uniformly almost everywhere. And \[\int_0^T\int_0^T\Big|\alpha(t,s) - \sum_{j=1}^{r}\lambda_jV_j^*(t)V_j(s)\Big|^2\mathrm{d}s\mathrm{d}t = \sum_{j=r+1}^{\infty}\lambda_j^2.\]
Let \(\psi_t\) and \(\psi_t^r\) be the mild solutions of 1 and 9 , respectively, with the same initial condition. In both equations, the noise \(z\) is generated by \(\alpha(t,s)\). The next theorem states convergence of the low-rank approximation.
Theorem 3. Under Assumption [assum:exist95uniq], there exists a constant \(r_0>0\) such that for all \(r>r_0\), we have \[\mathbb{E}\left(\sup_{t\in[0,T]}\left\|\psi_t-\psi_t^{r}\right\|_{\mathcal{H}}\right) \leq C\varepsilon_r, \label{eq:converg95r}\qquad{(4)}\] where \(C\) is independent of \(r\), and \(\varepsilon_r\) is defined as ?? .
Proof. It follows from 2 that \[\int_0^T\int_0^T\left|\alpha(t,s)-\widehat{\alpha}_r(t,s)\right|^2\mathrm{d}t\mathrm{d}s = \varepsilon_r^2.\] Applying the Cauchy-Schwarz inequality, we obtain \[\int_0^T\int_0^T\left|\alpha(t,s)-\widehat{\alpha}_r(t,s)\right|\mathrm{d}t\mathrm{d}s \leq T\varepsilon_r.\] According to 2, \(\psi_t\) and \(\psi_t^{r}\) admit the following series representations: \[\psi_t = S(t)\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\widetilde{\psi}_t^{(n,m)}, \qquad \psi_t^{r} = S(t)\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\widetilde{\psi}_t^{r(n,m)}.\] By exploiting the symmetry of the integrand, we have \[\begin{align} g(t,\alpha,\widehat{\alpha}_r) &\coloneq \int_{0<s_1<\cdots<s_{2m}<t}\mathrm{d}^{2m}\boldsymbol{s}\sum_{\boldsymbol{q}\in\mathscr{Q}(\boldsymbol{s})}\left|\prod_{j=1}^{m}\alpha(q_{2j},q_{2j-1})-\prod_{j=1}^{m}\widehat{\alpha}_{r}(q_{2j},q_{2j-1})\right| \\ &= \frac{(2m-1)!!}{(2m)!} \int_{[0,t]^{2m}} \mathrm{d}^{2m}\boldsymbol{s} \left|\prod_{j=1}^{m}\alpha(s_{2j},s_{2j-1})-\prod_{j=1}^{m}\widehat{\alpha}_{r}(s_{2j},s_{2j-1})\right|. \label{eq:49522} \end{align}\tag{20}\] To estimate the integrand in 20 , note that \[\begin{align} & \left|\prod_{j=1}^{m}\alpha(s_{2j},s_{2j-1})-\prod_{j=1}^{m}\widehat{\alpha}_{r}(s_{2j},s_{2j-1})\right| \\ \leq& \sum_{j=1}^{m} \prod_{k=j+1}^{m} |\widehat{\alpha}_r(s_{2k},s_{2k-1})| \cdot \left|\alpha(s_{2j},s_{2j-1})-\widehat{\alpha}_r(s_{2j},s_{2j-1})\right| \cdot \prod_{k=1}^{j-1} |\alpha(s_{2k},s_{2k-1})|. \label{eq:49523} \end{align}\tag{21}\] Let \(A\coloneq \int_0^T\int_0^T|\alpha(t,s)|\mathrm{d}t\mathrm{d}s\). Substituting 21 into 20 gives \[\begin{align} g(t,\alpha,\widehat{\alpha}_r) \leq g(T, \alpha, \widehat{\alpha}_r) \leq \frac{1}{m!2^m}\left[\left(A+T\varepsilon_r\right)^m - A^m\right]. \end{align}\] Consequently, we have \[\begin{align} & \left\|\widetilde{\psi}_t^{(n,m)} - \widetilde{\psi}_t^{r(n,m)}\right\|_{\mathcal{H}} \\ &\quad \leq \frac{\|L\|_{\mathcal{H}}^n}{n!}\left(\int_0^t|z_s|\mathrm{d}s\right)^n \frac{\|L\|_{\mathcal{H}}^{2m}}{m!}\left[\left(\frac{A+T\varepsilon_r}{2}\right)^m - \left(\frac{A}{2}\right)^m\right] \|\psi_0\|_{\mathcal{H}}. \end{align}\] Summing over \(n\) and \(m\), and invoking the absolute convergence of the series, we obtain \[\begin{align} \sup_{t\in[0,T]}\left\|\psi_t-\psi_t^{r}\right\|_{\mathcal{H}} \leq& \sup_{t\in[0,T]}\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\left\|\widetilde{\psi}_t^{(n,m)} - \widetilde{\psi}_t^{r(n,m)}\right\|_{\mathcal{H}} \\ \leq& e^{\|L\|_{\mathcal{H}}\int_{0}^T|z_s|\mathrm{d}s}\left(e^{\frac{\|L\|_{\mathcal{H}}^2}{2}(A+T\varepsilon_r)} - e^{\frac{\|L\|_{\mathcal{H}}^2}{2}A}\right) \|\psi_0\|_{\mathcal{H}}. \end{align}\] Finally, taking the expectation and applying the estimate from 1, we establish the desired result ?? . ◻
4.2 Convergence of the finite truncation↩︎
In this subsection, we establish convergence of the finitely truncated hierarchical system 13 . We equip \(\mathbb{N}^r\) with the graded lexicographic order: for \(\boldsymbol{j},\boldsymbol{k}\in\mathbb{N}^r\), we write \(\boldsymbol{j}\prec\boldsymbol{k}\) if either \(|\boldsymbol{j}|<|\boldsymbol{k}|\), or \(|\boldsymbol{j}|=|\boldsymbol{k}|\) and \(\boldsymbol{j}\) precedes \(\boldsymbol{k}\) in the usual lexicographic order. Utilizing the operator \(D_j\) defined in 10 , let \(\Psi_t^r\) denote the infinite dimensional vector: \[\Psi_t^r\mathrel{\vcenter{:}}=\left(\prod_{j=1}^r D_j^{k_j}\psi^r_t\right)_{\boldsymbol{k}\in\mathbb{N}^r},\] where \(\psi_t^r\) is the mild solution of 9 . To characterize the solution of the hierarchical system, we introduce the Banach space \[\mathcal{H}_w \mathrel{\vcenter{:}}= \left\{\Psi=\left(\psi^{(\boldsymbol{k})}\right): \boldsymbol{k}\in\mathbb{N}^r,\;\;\psi^{(\boldsymbol{k})}\in\mathcal{H},\;\; \|\Psi\|_{w}<\infty\;\;\mathrm{a.s.}\right\},\] with the norm \(\|\cdot\|_w\) defined as \[\left\|\Psi\right\|_w^2 \mathrel{\vcenter{:}}= \sum_{\boldsymbol{k}\in\mathbb{N}^r} \prod_{j=1}^{r}\frac{\lambda_{j}^{k_j}}{k_j!}\left\|\psi^{(\boldsymbol{k})}\right\|_\mathcal{H}^2,\] where \(\lambda_j\) correspond to the low-rank BCF approximation in 8 . The following lemma shows that \(\Psi_t^r\in\mathcal{H}_w\) for \(t\in[0,T]\).
Lemma 3. The norm of each auxiliary component \(\prod_{j=1}^rD_j^{k_j}\psi^r_t\) satisfies the following upper bound: \[\left\|\prod_{j=1}^r D_j^{k_j}\psi^r_t\right\|_{\mathcal{H}} \leq \|L\|_{\mathcal{H}}^{|\boldsymbol{k}|}\prod_{j=1}^r\left(\int_0^t|V_j(s)|\mathrm{d}s\right)^{k_j} \mathcal{M}_t^r \|\psi_0\|_{\mathcal{H}}, \label{eq:estimate95D95j}\qquad{(5)}\] where \[\label{eq:M95t95r} \mathcal{M}_t^r \coloneq \exp\left\{\|L\|_{\mathcal{H}}\int_{0}^t|z_s|\mathrm{d}s + \frac{\|L\|^2_{\mathcal{H}}}{2}\int_0^t\int_0^t|\widehat{\alpha}_r(u,v)|\mathrm{d}u\mathrm{d}v\right\}.\qquad{(6)}\] Here, \(\widehat{\alpha}_r\) is the rank-\(r\) approximation of \(\alpha\). It follows that \(\Psi^r_t\in\mathcal{H}_w\), i.e., \[\mathbb{P}\left(\left\|\Psi^r_t\right\|_{w}<\infty\right) = 1.\]
Proof. We consider the infinite series representation of \(\psi_t^r\) in 2. We have an estimate of \(D_j\widetilde{\psi}_t^{r(n,m)}\) following the procedure of 14 –17 : \[\begin{align} \left\|D_j\widetilde{\psi}_t^{r(n,m)}\right\|_{\mathcal{H}} \leq& \left(\|L\|_{\mathcal{H}}\int_0^t|V_j(s)|\mathrm{d}s\right) \frac{1}{(n-1)!}\left(\|L\|_{\mathcal{H}}\int_0^t|z_s|\mathrm{d}s\right)^{n-1} \\ &\frac{1}{m!}\left(\frac{\|L\|_{\mathcal{H}}^2}{2}\int_0^t\int_0^t|\widehat{\alpha}_r(u,v)|\mathrm{d}u\mathrm{d}v\right)^m \|\psi_0\|_{\mathcal{H}}. \end{align}\] Summing them up yields \[\left\|D_j\psi_t^r\right\|_{\mathcal{H}} \leq \left(\|L\|_{\mathcal{H}}\int_0^t|V_j(s)|\mathrm{d}s\right) \mathcal{M}_t^r \|\psi_0\|_{\mathcal{H}}.\] By recursively applying the operator \(D_j\) and iterating the aforementioned procedure, we obtain the bound ?? . Substituting the estimate into the norm \(\|\cdot\|_w\) implies \[\|\Psi_t^r\|_w^2 \leq \sum_{\boldsymbol{k}\in\mathbb{N}^r}\prod_{j=1}^{r}\frac{1}{k_j!}\left(\|L\|_{\mathcal{H}}\int_0^t\sqrt{\lambda_j}|V_j(s)|\mathrm{d}s\right)^{2k_j} \left(\mathcal{M}_t^r\right)^2\|\psi_0\|_{\mathcal{H}}^2. \label{eq:49529}\tag{22}\] Furthermore, an application of the Cauchy-Schwarz inequality yields \[\begin{align} \left(\int_0^t\sqrt{\lambda_j}|V_j(s)|\mathrm{d}s\right)^2 \leq \int_0^t1^2\mathrm{d}s\int_0^t\lambda_jV_j^*(s)V_j(s)\mathrm{d}s \leq t\int_0^t\alpha(s,s)\mathrm{d}s. \label{eq:49530} \end{align}\tag{23}\] Combining 22 and 23 , we arrive at \[\|\Psi_t^r\|_w^2 \leq \sum_{\boldsymbol{k}\in\mathbb{N}^r}\prod_{j=1}^{r}\frac{1}{k_j!}\left(t\|L\|_{\mathcal{H}}^2\int_0^t\alpha(s,s)\mathrm{d}s\right)^{k_j} \left(\mathcal{M}_t^r\right)^2\|\psi_0\|_{\mathcal{H}}^2 < \infty, \quad \text{a.s.},\] which completes the proof. ◻
For the following analysis, we now turn to the interaction picture of the truncated hierarchical system 13 , which reads \[\label{eq:finite95hierarchy95interaction} \frac{\partial}{\partial t} \widetilde{\psi}_t^{(\boldsymbol{k})} = z_t\widetilde{L}(t)\widetilde{\psi}_t^{(\boldsymbol{k})} + \widetilde{L}(t)\sum_{\substack{j=1 \\ \boldsymbol{k}-\boldsymbol{e}_j\in\mathcal{K}^r_N}}^{r}k_jV_j(t)\widetilde{\psi}_t^{(\boldsymbol{k}-\boldsymbol{e}_j)} - \widetilde{L}^{\dagger}(t)\sum_{\substack{j=1\\ \boldsymbol{k}+\boldsymbol{e}_j\in\mathcal{K}^r_N}}^{r}\lambda_jV_j^*(t)\widetilde{\psi}_t^{(\boldsymbol{k}+\boldsymbol{e}_j)}.\tag{24}\] Let \(\widetilde{\Psi}_t^{r,N}\) denote the strong solution to 24 subject to the initial conditions \[\widetilde{\psi}_0^{(\boldsymbol{0})} = \widetilde{\psi}_0 = \psi_0, \quad \text{and}\quad \widetilde{\psi}_0^{(\boldsymbol{k})} = 0,\quad \forall \boldsymbol{k}\in\mathcal{K}^r_N, |\boldsymbol{k}|>0.\] We embed \(\widetilde{\Psi}_t^{r,N}\) into \(\mathcal{H}_w\) by defining its components as \(\widetilde{\psi}_t^{(\boldsymbol{k})}\) for \(\boldsymbol{k} \in \mathcal{K}_N^r\) and setting \(\widetilde{\psi}_t^{(\boldsymbol{k})} = 0\) for all \(\boldsymbol{k} \notin \mathcal{K}_N^r\). Consequently, to establish the convergence of \(\psi_t^{r,N}\) to \(\psi_t^r\), it suffices to show that \(\widetilde{\Psi}_t^{r,N}\) converges to \(\widetilde{\Psi}_t^r\) in \(\mathcal{H}_w\) as \(N\to\infty\).
The interaction picture of the infinite hierarchical system can be reformulated as \[\frac{\partial}{\partial t}\widetilde{\Psi}_t^{r} = \mathcal{A}(t)\widetilde{\Psi}_t^r,\] where \(\mathcal{A}(t)\) serves as the infinitesimal generator of the underlying propagator. We further define the projection operator \(P_N \in \mathcal{L}(\mathcal{H}_w)\), which maps the infinite hierarchy onto its finite-dimensional subspace by retaining components indexed by \(\boldsymbol{k} \in \mathcal{K}_N^r\) and mapping all higher-order components to zero; specifically, \[(P_N\Psi)_{\boldsymbol{k}} = \begin{cases} \psi^{(\boldsymbol{k})}, & \boldsymbol{k}\in\mathcal{K}^r_N, \\ 0, & \boldsymbol{k}\notin\mathcal{K}^r_N. \end{cases}\] The finite system 24 is reformulated as \[\frac{\partial}{\partial t} \widetilde{\Psi}_t^{r,N} = \mathcal{A}_N(t)\widetilde{\Psi}_t^{r,N}, \label{eq:49534}\tag{25}\] where \(\mathcal{A}_N(t) \mathrel{\vcenter{:}}= P_N\mathcal{A}(t)P_N\). Consider the error function \(E_N(t)\mathrel{\vcenter{:}}= \widetilde{\Psi}_t^{r,N}-P_N\widetilde{\Psi}_t^r\) subject to \(E_N(0)=0\). By construction, \(E_N(t)\) satisfies \[\begin{align} \frac{\partial}{\partial t} E_N(t) =& \frac{\partial}{\partial t}\widetilde{\Psi}_t^{r,N} - \frac{\partial}{\partial t}\left(P_N\widetilde{\Psi}_t^{r}\right) \\ =& \mathcal{A}_N(t)\widetilde{\Psi}_t^{r,N} - P_N\mathcal{A}(t)\widetilde{\Psi}_t^r \\ =& \mathcal{A}_N(t)\left(\widetilde{\Psi}_t^{r,N}-P_N\widetilde{\Psi}_t^r\right) + \left(\mathcal{A}_N(t)P_N - P_N\mathcal{A}(t)\right)\widetilde{\Psi}_t^r \\ =& \mathcal{A}_N(t)E_N(t) + P_N\mathcal{A}(t)(P_N-I)\widetilde{\Psi}_t^r, \end{align}\] where \(I\) denotes the identity operator in \(\mathcal{L}(\mathcal{H}_w)\) and the last step follows from \(P_N^2=P_N\). To derive a bound for the error \(E_N(t)\), we first establish a stability estimate for 25 using an energy-based approach.
Lemma 4. The strong solution of 25 has an estimate \[\left\|\widetilde{\Psi}_t^{r,N}\right\|_w \leq \left\|\widetilde{\Psi}_0^{r,N}\right\|_w \exp\left\{\|L\|_{\mathcal{H}}\int_0^t|z_s|\mathrm{d}s\right\}. \label{eq:49536}\qquad{(7)}\]
Proof. Taking the time derivative of \(\frac{1}{2}\|\widetilde{\Psi}_t^{r,N}\|_w^2\) and substituting 24 yields \[\begin{align} \frac{1}{2}\frac{\partial}{\partial t}\left\|\widetilde{\Psi}_t^{r,N}\right\|_w^2 =\,& \mathrm{Re}\Bigg\{\sum_{\boldsymbol{k}\in\mathcal{K}_N^r}\prod_{j=1}^{r}\frac{\lambda_j^{k_j}}{k_j!}\left\langle\widetilde{\psi}_t^{(\boldsymbol{k})}, z_t\widetilde{L}(t)\widetilde{\psi}_t^{(\boldsymbol{k})}\right\rangle_{\mathcal{H}} \\ &+ \sum_{\boldsymbol{k}\in\mathcal{K}_N^r}\sum_{\substack{\ell=1\\ \boldsymbol{k}-\boldsymbol{e}_{\ell}\in\mathcal{K}_N^r}}^{r}\prod_{j=1}^{r}\frac{\lambda_j^{k_j}}{k_j!}\left\langle\widetilde{\psi}_t^{(\boldsymbol{k})}, k_{\ell}V_{\ell}(t)\widetilde{L}(t)\widetilde{\psi}_t^{(\boldsymbol{k}-\boldsymbol{e}_{\ell})}\right\rangle_{\mathcal{H}} \\ &- \sum_{\boldsymbol{k}\in\mathcal{K}_N^r}\sum_{\substack{\ell=1\\ \boldsymbol{k}+\boldsymbol{e}_{\ell}\in\mathcal{K}_N^r}}^{r}\prod_{j=1}^{r}\frac{\lambda_{j}^{k_j}}{k_j!}\left\langle\widetilde{\psi}_t^{(\boldsymbol{k})},\lambda_{\ell}V_{\ell}^*(t)\widetilde{L}^{\dagger}(t)\widetilde{\psi}_t^{(\boldsymbol{k}+\boldsymbol{e}_{\ell})}\right\rangle_{\mathcal{H}}\Bigg\}. \label{eq:energy95estimate} \end{align}\tag{26}\] Note that \(\forall\phi,\psi\in\mathcal{H}\), the following identity holds: \[\left\langle\phi, V_{\ell}^*(t)\widetilde{L}^{\dagger}(t)\psi\right\rangle_{\mathcal{H}} = \left\langle V_{\ell}(t)\widetilde{L}(t)\phi, \psi \right\rangle_{\mathcal{H}} = \left\langle\psi,V_{\ell}(t)\widetilde{L}(t)\phi\right\rangle_{\mathcal{H}}^*.\] By relabeling the indices, the last two terms on the right-hand side of 26 vanish. Consequently, we obtain the following inequality: \[\begin{align} \frac{1}{2}\frac{\partial}{\partial t}\left\|\widetilde{\Psi}_t^{r,N}\right\|_w^2 \leq |z_t| \|L\|_{\mathcal{H}}\left\|\widetilde{\Psi}_t^{r,N}\right\|_w^2. \end{align}\] Finally, an application of Gronwall’s inequality yields the desired estimate ?? . ◻
With the stability estimate from 4, we now proceed to establish the convergence by invoking Duhamel’s principle.
Theorem 4. Let \(\psi_t^{r,N}\) and \(\psi_t^{r}\) denote the \(\boldsymbol{0}\)-th solution of 13 and the mild solution of 9 , respectively. There exist constants \(C, \beta > 0\) independent of \(r\) and \(N\) such that \[\label{eq:truncation95converg95order} \mathbb{E}\left(\sup_{t\in[0,T]}\left\|\psi_t^{r} - \psi_t^{r,N}\right\|_{\mathcal{H}}\right) \leq C\sqrt{\frac{(\beta r)^N}{N!}}.\qquad{(8)}\]
Proof. Let \(\boldsymbol{v}^{s}\) denote the solution of the homogeneous equation \(\frac{\partial}{\partial t}\boldsymbol{v}^s = \mathcal{A}_N(t)\boldsymbol{v}^s\) with the initial condition \(\boldsymbol{v}^s(s) = P_N\mathcal{A}(s)(P_N-I)\widetilde{\Psi}_s^r\). By Duhamel’s principle, we have \(E_N(t) = \int_0^t \boldsymbol{v}^s(t)\mathrm{d}s\). By 4, we have \[\begin{align} \left\|\boldsymbol{v}^s(t)\right\|_w \leq& \exp\left\{\|L\|_{\mathcal{H}}\int_s^t|z_u|\mathrm{d}u\right\} \left\|P_N\mathcal{A}(s)(P_N-I)\widetilde{\Psi}_s^r\right\|_w. \end{align}\] Consequently, \[\begin{align} \left\|E_N(t)\right\|_w \leq& \exp\left\{\|L\|_{\mathcal{H}}\int_0^t|z_s|\mathrm{d}s\right\} \int_0^t \left\|P_N\mathcal{A}(s)(P_N-I)\widetilde{\Psi}_s^r\right\|_w \mathrm{d}s. \label{eq:49542} \end{align}\tag{27}\] It remains to bound \(\left\|P_N\mathcal{A}(s)(P_N-I)\widetilde{\Psi}_s^r\right\|_w\). We observe that \(P_N \mathcal{A}(s)(P_N - I) \widetilde{\Psi}_s^r\) vanishes except for those components associated with multi-indices \(\boldsymbol{k}\) on the truncation boundary, satisfying \(|\boldsymbol{k}| = N\). For any such multi-index, the corresponding \(\boldsymbol{k}\)-th component is given by \[\sum_{j=1}^r \lambda_jV_j^*(s)\widetilde{L}^{\dagger}(s)\widetilde{\psi}_s^{(\boldsymbol{k}+\boldsymbol{e}_j)} \coloneq \sum_{j=1}^r \lambda_jV_j^*(s)\widetilde{L}^{\dagger}(s) D_j\prod_{m=1}^r D_m^{k_m} \widetilde{\psi}_s^r.\] By definition, \[\int_0^t \left\|P_N\mathcal{A}(s)(P_N-I)\widetilde{\Psi}_s^r\right\|_w \mathrm{d}s = \int_0^t\sqrt{\sum_{|\boldsymbol{k}|=N}\prod_{j=1}^r\frac{\lambda_j^{k_j}}{k_j!}\left\|\sum_{j=1}^r\lambda_jV_j^*(s)\widetilde{L}^{\dagger}(s)\widetilde{\psi}_s^{(\boldsymbol{k}+\boldsymbol{e}_j)}\right\|_{\mathcal{H}}^2}\mathrm{d}s. \label{eq:49544}\tag{28}\] Applying the Cauchy-Schwarz inequality and 3, we obtain: \[\begin{align} \left\|\sum_{j=1}^r\lambda_jV_j^*(s)\widetilde{\psi}_s^{(\boldsymbol{k}+\boldsymbol{e}_j)}\right\|_{\mathcal{H}}^2 \leq& \Bigg(\sum_{j=1}^r \lambda_jV_j^*(s)V_j(s)\Bigg)\cdot \Bigg(\|L\|_{\mathcal{H}}^{|\boldsymbol{k}|+1}\mathcal{M}_t^r\|\psi_0\|_{\mathcal{H}} \\ & \prod_{\ell=1}^r\left(\int_0^t|V_{\ell}(u)|\mathrm{d}u\right)^{k_{\ell}}\sum_{j=1}^r\int_0^t\sqrt{\lambda_j}|V_j(u)|\mathrm{d}u \Bigg)^2, \label{eq:49545} \end{align}\tag{29}\] where \(\mathcal{M}_t^r\) is defined in ?? . Here, we have utilized the fact that \(s\leq t\) and substituted \(t\) into the final term on the right-hand side. Applying Cauchy-Schwarz again, together with \(\sum_{j=1}^r\lambda_jV_j^*(s)V_j(s) \leq \alpha(s,s)\), we obtain \[\sum_{j=1}^r \left(\int_0^t \sqrt{\lambda_j}|V_j(u)|\mathrm{d}u\right)^2 \leq t\int_0^t\alpha(s,s)\mathrm{d}s. \label{eq:49546}\tag{30}\] Combining 23 , 27 , 28 , 29 , and 30 yields \[\begin{align} \label{eq:49554} \sup_{t\in[0,T]}\left\|\psi_t^{r} - \psi_t^{r,N}\right\|_{\mathcal{H}} =& \sup_{t\in[0,T]}\left\|\widetilde{\psi}_t^{r} - \widetilde{\psi}_t^{r,N}\right\|_{\mathcal{H}} \leq \sup_{t\in[0,T]}\left\|E_N(t)\right\|_w \\ \leq& \exp\left\{\|L\|_{\mathcal{H}}\int_0^T|z_s|\mathrm{d}s\right\} \|L\|_{\mathcal{H}}^2\mathcal{M}_T^r\|\psi_0\|_{\mathcal{H}}\int_0^T\sqrt{\alpha(s,s)}\mathrm{d}s \\ &\;T\int_0^T\alpha(s,s)\mathrm{d}s \sqrt{\sum_{|\boldsymbol{k}|=N}\prod_{j=1}^r \frac{1}{k_j!}\left(T\|L\|_{\mathcal{H}}^2\int_0^T\alpha(s,s)\mathrm{d}s\right)^{k_j}}. \end{align}\tag{31}\] By Cauchy-Schwarz inequality, \[\begin{align} \label{eq:49555} \int_0^T\int_0^T\left|\widehat{\alpha}_r(t,s)\right|\mathrm{d}t\mathrm{d}s \leq& \,T\sqrt{\int_0^T\int_0^T\left|\widehat{\alpha}_r(t,s)\right|^2\mathrm{d}t\mathrm{d}s} \\ \leq& \,T\sqrt{\int_0^T\int_0^T\left|\alpha(t,s)\right|^2\mathrm{d}t\mathrm{d}s}. \end{align}\tag{32}\] Note that for any \(\beta > 0\), it holds that \[\sum_{|\boldsymbol{k}|=N}\prod_{j=1}^r\frac{\beta^{k_j}}{k_j!} = \frac{(\beta r)^N}{N!}.\] Substituting 32 into \(\mathcal{M}_T^r\) in 31 , taking expectation on both sides of 31 , and applying 1, we complete the proof. ◻
4.3 Proof of Theorem 1↩︎
The proof of 1 follows directly by synthesizing the results of 3 and 4.
5 Implementation↩︎
In this section, we describe the numerical generation of the colored noise paths and the low-rank approximation of the BCF. Furthermore, we present the numerical schemes employed to solve the finitely truncated hierarchical system 13 .
5.1 Generation of the noise path and eigenfunctions↩︎
The colored noise path \(z_t\) is generated by Cholesky factorization. We discretize \([0,T]\) using \(t_i=i\Delta t\), \(\Delta t=T/M\), and set \(\boldsymbol{z}=(z_{t_0},\ldots,z_{t_M})^{\top}\). Its covariance matrix is \(S_{ij}=\alpha^*(t_i,t_j)\). Given a factorization \(S=LL^{\dagger}\), one sample path is obtained from \[\boldsymbol{z}=L\boldsymbol{\xi}, \qquad \boldsymbol{\xi}\sim\mathcal{CN}(0,I_{M+1}).\]
The eigenpairs of \(\alpha\) are computed from the grid matrix \(A_{ij}=\alpha(t_i,t_j)\). If \(A=U\Sigma V^{\dagger}\) is its singular value decomposition, then the discrete eigenvalues are approximated by the diagonal entries of \(\Delta t\,\Sigma\), while the sampled eigenfunctions are given by the columns of \(V/\sqrt{\Delta t}\).
5.2 Numerical scheme for the finite hierarchical system↩︎
We employ a Strang splitting scheme to solve the finite hierarchical system 13 . This approach is motivated by the non-commutativity between the Hamiltonian \(H\) (representing the conservative dynamics) and the operator \(L\) (associated with the dissipative terms). Let \(\Delta t > 0\) denote the time step, and define the discrete time points \(t_n = n\Delta t\) for \(n=0, 1, \dots\). The Strang splitting update for the state \(\Psi^{r,N}_t\) over one time step \([t_n, t_{n+1}]\) is formulated as \[\label{eq:strang95splitting} \begin{array}{lll} \displaystyle \frac{\partial}{\partial t}u_t^{(\boldsymbol{k})} = -iHu_t^{(\boldsymbol{k})}, &\quad u_{t_n}^{(\boldsymbol{k})} = \psi_{t_n}^{(\boldsymbol{k})}, &\quad t\in\left[t_n,t_{n+1/2}\right], \\[2ex] \displaystyle \frac{\partial}{\partial t}v_t^{(\boldsymbol{k})} = f^{(\boldsymbol{k})}(t,\boldsymbol{v}_t), &\quad v_{t_{n}}^{(\boldsymbol{k})}=u_{t_{n+1/2}}^{(\boldsymbol{k})}, &\quad t\in\left[t_n, t_{n+1}\right], \\[2ex] \displaystyle \frac{\partial}{\partial t}\phi_t^{(\boldsymbol{k})} = -iH\phi_t^{(\boldsymbol{k})}, &\quad \phi_{t_{n+1/2}}^{(\boldsymbol{k})} = v_{t_{n+1}}^{(\boldsymbol{k})}, &\quad t\in[t_{n+1/2}, t_{n+1}], \end{array}\tag{33}\] for \(\boldsymbol{k}\in\mathcal{K}^r_N\), where \[f^{(\boldsymbol{k})}(t,\boldsymbol{v}_t) = z_tLv_t^{(\boldsymbol{k})} + L\sum_{\substack{j=1 \\ \boldsymbol{k}-\boldsymbol{e}_j\in\mathcal{K}^r_N}}^{r}k_jV_j(t)v_t^{(\boldsymbol{k}-\boldsymbol{e}_j)} - L^{\dagger}\sum_{\substack{j=1 \\ \boldsymbol{k}+\boldsymbol{e}_j\in\mathcal{K}^r_N}}^{r}\lambda_jV_j^*(t)v_t^{(\boldsymbol{k}+\boldsymbol{e}_j)}.\] The numerical solution at the \((n+1)\)-th time step is obtained as \(\psi_{t_{n+1}}^{(\boldsymbol{k})} = \phi_{t_{n+1}}^{(\boldsymbol{k})}\). Regarding the splitting steps defined in 33 , the first and third subproblems correspond to the unitary evolution under the system Hamiltonian, which is computed via the propagator \(S(\Delta t/2)\). The second subproblem, which involves the dissipative and stochastic terms, is integrated using a second-order Runge-Kutta method. In our implementation, we specifically employ Heun’s method, which is formulated as follows: \[\begin{align} w_1^{(\boldsymbol{k})} =& f^{(\boldsymbol{k})}(t_n,\boldsymbol{v}_{t_n}), \\ w_2^{(\boldsymbol{k})} =& f^{(\boldsymbol{k})}\left(t_{n+1},\boldsymbol{v}_{t_n} + \Delta t \boldsymbol{w}_1\right), \\ v_{t_{n+1}}^{(\boldsymbol{k})} =& v_{t_n}^{(\boldsymbol{k})} + \frac{\Delta t}{2}\left(w_1^{(\boldsymbol{k})} + w_2^{(\boldsymbol{k})}\right). \end{align}\] A key computational advantage of this scheme in the present context is that it circumvents the requirement for evaluating the stochastic path \(z_t\) and the eigenfunctions \(V_j(t)\) at intermediate time steps.
6 Numerical experiments↩︎
In this section, we present several numerical experiments to validate the proposed method for the NMSSE. We consider the two-level spin-boson model coupled to a bosonic bath, which is a benchmark model commonly used in the study of open quantum systems. Consequently, the system Hilbert space is identified as \(\mathcal{H}=\mathbb{C}^2\). For any vectors \(\phi, \psi \in \mathbb{C}^2\), the inner product and its associated norm are now defined as \[\langle \phi, \psi \rangle = \phi^{\dagger}\psi, \quad\mathrm{and}\quad \|\psi\| = \sqrt{\langle \psi, \psi \rangle}.\] The expectation value of an observable \(O\in\mathcal{L}(\mathbb{C}^2)\) at time \(t\) is given by \[\langle O \rangle_t \coloneq \frac{\mathbb{E}\left(\left\langle\psi_t,O\psi_t\right\rangle\right)}{\mathbb{E}\left(\left\langle\psi_t,\psi_t\right\rangle\right)},\] where \(\psi_t\) is the solution to the NMSSE. To evaluate the accuracy of the proposed method, we benchmark our results against a reference solution obtained via the Time-Evolving Matrix Product Operators (TEMPO) method [30], which serves as a high-precision numerical standard for non-Markovian dynamics.
The decay of the eigenvalues in the low-rank approximation is governed by the regularity of the BCF: if \(\alpha\in \mathcal{C}^p\), then \(\lambda_n=o(n^{-(p+1)})\) as \(n\to\infty\) [31]. Hence smoother BCFs require smaller numerical ranks. For a spectral density \(J\), the BCF is generally given by \[\alpha(t,s) = \int_0^{\infty} J(\omega) \left[ \coth\left(\frac{\beta\omega}{2}\right) \cos(\omega(t-s)) - i \sin(\omega(t-s)) \right] \mathrm{d}\omega,\] where \(\beta\) is the inverse temperature and \[J(\omega) = \eta \omega^s f(\omega/\omega_c).\] Here \(\eta>0\), \(s>0\), \(\omega_c\) is the cutoff frequency, and \(f\) specifies the high-frequency cutoff, e.g., \(f(x)=e^{-x}\) or \(f(x)=(x^2+1)^{-1}\). For the exponential cutoff with \(s\geq 1\), \(\alpha\) is smooth and the eigenvalues decay rapidly, so a small rank is sufficient. For less regular choices such as the Drude–Lorentz cutoff, one may instead use a Matsubara sum-of-exponentials representation [14]; in this setting our relaxed formulation recovers HOPS, as discussed in Section 6.2.
6.1 Exponential cutoff BCF with Ohmic spectral density↩︎
We first evaluate the performance of our method using a spin-boson model with the following configuration:
The system Hamiltonian is given by \(H=\varepsilon\sigma_z + \sigma_x\), where \(\sigma_z\) and \(\sigma_x\) are the Pauli matrices.
The Lindblad operator is given by \(L=\sigma_z\).
The bath correlation function is given by [32] \[\alpha(t,s) = \sum_{j=1}^J\frac{c_{j}^2}{2\omega_{j}}\left[\coth\left(\frac{\beta\omega_j}{2}\right)\cos\left(\omega_j(t-s)\right) - i\sin\left(\omega_j(t-s)\right)\right], \label{eq:BCF95Ohmic}\tag{34}\] where \[\begin{align} \omega_j &= -\omega_c\log\left(1-\frac{j}{J}\left[1-\exp(-\omega_{\mathrm{max}}/\omega_c)\right]\right), \\ c_j &= \omega_j\sqrt{\frac{\xi\omega_c}{J}\left[1-\exp(-\omega_{\mathrm{max}}/\omega_c)\right]}. \end{align}\]
The physical parameters are set as \(\beta=5\), \(\omega_c=2.5\), and \(\omega_{\mathrm{max}}=4\omega_c\), with the number of bath modes \(J=200\). We conduct experiments by varying the energy bias \(\varepsilon\) and the coupling strength \(\xi\). The system is initialized in the state \(\psi_0 = (1,0)^{\top}\). Numerical simulations are performed over the time interval \([0,5]\) with a uniform time step \(\Delta t=0.05\).
This BCF 34 can be written as a sum of exponentials: \[\alpha(t,s) = \sum_{j=1}^J \frac{c_j^2}{4\omega_j}\left(g_je^{i\omega_j(t-s)}+h_je^{-i\omega_j(t-s)}\right),\] where \(g_j=\coth\left(\beta\omega_j/2\right)-1\) and \(h_j=\coth\left(\beta\omega_j/2\right)+1\). However, it is formidable to solve this problem with HOPS, since there are 400 exponentials in the BCF, which leads to a multi-index of length 400 to represent the auxiliary wave functions. The hierarchical system is thus too large to be solved in practice. In contrast, the proposed method can solve this problem efficiently by using a low-rank approximation of the BCF. 1 shows that the BCF 34 can be well approximated by a low-rank approximation with rank \(r=10\).
2 illustrates the evolution of \(\langle\sigma_z\rangle_t\) for various values of \(\varepsilon=0,1,2\), \(\xi=0.2,0.4\), rank \(r=10\), and truncation order \(N=8\). Under this configuration, the total number of hierarchical equations is given by \[\sum_{n=0}^N\binom{n+r-1}{r-1} = \binom{N+r}{r} = \binom{18}{10} = 43758.\] This cardinality is significantly lower than that required by the standard HOPS method, which entails \(\binom{408}{400}\approx 10^{15}\) equations. 2 also compares the results obtained using different sample sizes for \(z\). These numerical trajectories demonstrate that the proposed method converges to the reference solution as the sample size increases.
Figure 2: Evolution of \(\langle\sigma_z\rangle_t\) for different parameter settings (top to bottom: \(\xi=0.2,0.4\))..
Since the size of the truncated hierarchy grows rapidly with the rank \(r\), we restrict the convergence-order test to the truncation order \(N\). For each fixed rank \(r=2,4,8\), we vary \(N\) and take the solution with \(N=16\) as the reference solution, denoted by \(\psi_t^{r,\mathrm{ref}}\). The error is measured by \[e_N\coloneq \mathbb{E}\left(\sup_{t\in[0,T]}\left\|\psi_t^{r,\mathrm{ref}}-\psi_t^{r,N}\right\|\right).\] According to ?? , for fixed \(r\) this error is expected to satisfy \[e_N=\mathcal{O}\left(\sqrt{\frac{(\beta r)^N}{N!}}\right).\] To test this prediction, 3 reports two normalized quantities. The first is \(\log(e_N\sqrt{N!})\): after removing the factorial factor, the remaining dependence on \(N\) should be at most exponential, so the curves in 3 (a) are expected to be approximately affine. The second diagnostic is the rescaled consecutive error ratio. From the same estimate, it should satisfy \[\frac{e_{N+1}}{e_N}\sqrt{N+1}\approx \sqrt{\beta r},\] up to constants independent of \(N\). Thus the bounded and stabilizing profiles in 3 (b) give a complementary check of the predicted order. In addition, the variation of the slopes in 3 (a) and of the limiting levels in 3 (b) with respect to \(r\) is consistent with the \(r\)-dependent factor in ?? , further supporting the convergence estimate.
Figure 3: Convergence-order tests with respect to the truncation order \(N\) for fixed ranks \(r=2,4,8\)..
6.2 Exponentially decaying BCF↩︎
Next, we examine a spin-boson configuration characterized by an exponentially decaying BCF:
The system Hamiltonian is given by \(H=\frac{1}{2}\sigma_z\).
The Lindblad operator is given by \(L=\sqrt{2}\sigma_z\).
The bath correlation function is given by \[\alpha(t,s) = \frac{\gamma}{2}e^{-\gamma|t-s|}. \label{eq:simple95bcf}\tag{35}\]
The system is initialized in the state \(\psi_0 = (\sqrt{2}/2, \sqrt{2}/2)^{\top}\). Numerical simulations are performed over the time interval \([0,2]\) with a uniform time step \(\Delta t=0.02\).
Under this configuration, the problem has an exact solution. The reduced density matrix of the system at time \(t\) is given by \[\rho(t) = \begin{pmatrix} \rho_{11}(0) & \rho_{12}(0)e^{-f(t)} \\ \rho_{21}(0)e^{-f(t)^*} & \rho_{22}(0) \end{pmatrix},\] where \[f(t) = \frac{4}{\gamma}(e^{-\gamma t}+\gamma t-1) + i t.\] For any observable \(O\in\mathcal{L}(\mathbb{C}^2)\), the expectation value can be computed as \[\langle O \rangle_t = \mathrm{Tr}\left(O\rho(t)\right).\]
The BCF 35 is a single exponential corresponding to the simplest case for the HOPS method. Notably, this exponentially decaying BCF does not possess an intrinsic low-rank structure. 4 illustrates the amplitude of 35 with \(\gamma=4\) and its low-rank approximation at ranks \(r=10, 50\). Clearly, a high rank is required to achieve a sufficient approximation. Nevertheless, our proposed framework remains applicable by relaxing the approximation requirements. Specifically, we represent the BCF 35 using two functions \(U,V\in L^2([0,T],\mathbb{C})\) such that \[\alpha(t,s) = \lambda U(t)V(s), \qquad t\geq s, \label{eq:relax95approx}\tag{36}\] with \(\lambda=\gamma/2\), \(U(t) = e^{-\gamma t}\), and \(V(t) = e^{\gamma t}\). Substituting 36 into the NMSSE 1 and following the procedure established in Section 3, we obtain the hierarchical system \[\frac{\partial}{\partial t}\psi_t^{(k)} = (-iH+z_tL)\psi_t^{(k)} + k\lambda V(t)L\psi_t^{(k-1)} - U(t)L^{\dagger}\psi_t^{(k+1)} \label{eq:relax95hierarchy}\tag{37}\] for \(k=0,1,\cdots,N\) with \(\psi_t^{(-1)} \equiv \psi_t^{(N+1)} \equiv 0\). This system is equivalent to the one derived via the standard HOPS method; see 8 for a formal proof.
5 compares the numerical results obtained from the proposed framework 37 with those from the HOPS method. The agreement between the two approaches demonstrates that for exponential or sum-of-exponential BCFs, the proposed framework recovers the HOPS results by relaxing the low-rank approximation requirements. This observation suggests that the proposed framework serves as a generalization of the HOPS method, extending its applicability to a broader class of BCFs.
7 Conclusion↩︎
We have developed a low-rank hierarchical framework for the non-Markovian stochastic Schrödinger equation. By approximating the bath correlation function through Mercer’s theorem, the non-local stochastic dynamics are reformulated as a coupled hierarchical system that generalizes HOPS without relying on a prescribed multi-exponential expansion of the memory kernel. Under the uniqueness assumption, we proved convergence of the finite-rank, finitely truncated approximation in the \(L^1\) sense as both the rank \(r\) and the truncation order \(N\) increase. The numerical experiments support the theoretical convergence results and illustrate the effectiveness of the proposed method.
Several issues remain for future work. These include extending the analysis to spatially dependent models with unbounded system–bath couplings, establishing well-posedness results that remove the imposed uniqueness assumption, and deriving sharper error estimates to guide the choice of \(r\) and \(N\) in practical computations.
8 Equivalence between the relaxed hierarchical framework and the HOPS for exponential BCFs↩︎
In this section, we prove the equivalence between the relaxed hierarchical framework 37 and HOPS for the exponential BCF \(\alpha(t,s)=g e^{-w(t-s)}\) for \(t\geq s\). Setting \(\lambda = g\), \(U(t) = e^{-wt}\), and \(V(t) = e^{wt}\) in 37 yields \[\frac{\partial}{\partial t}\psi_t^{(k)} = (-iH+z_tL)\psi_t^{(k)} + kg e^{wt}L\psi_t^{(k-1)} - e^{-wt}L^{\dagger}\psi_t^{(k+1)}.\] Denote \(\widetilde{\psi}_t^{(k)} = e^{-kwt}\psi_t^{(k)}\). It follows that \[\begin{align} \frac{\partial}{\partial t}\widetilde{\psi}_t^{(k)} =& -kw e^{-kwt}\psi_t^{(k)} + (-iH+z_tL)e^{-kwt}\psi_t^{(k)} \\ &+ kgL e^{-(k-1)wt}\psi_t^{(k-1)} - L^{\dagger}e^{-(k+1)wt}\psi_t^{(k+1)} \\ =& (-iH - kw + z_tL)\widetilde{\psi}_t^{(k)} + kgL\widetilde{\psi}_t^{(k-1)} - L^{\dagger}\widetilde{\psi}_t^{(k+1)}, \end{align}\] which is identical to the hierarchical system derived via HOPS (cf. [23]). The same transformation applies to BCFs expressed as sums of exponentials.
Acknowledgments↩︎
The authors would like to thank Quanhui Zhu, Hongfei Zhan, and Shuigen Liu for their valuable suggestions and discussions.