1 Introduction↩︎

Continuous-variable open quantum systems are used to model the dynamics of a wide variety of processes, including quantum Brownian motion [1], [2], surface diffusion [3]–[5] and vibrational spectroscopy [6], [7]. The total Hamiltonian is usually written in the form [1], [8] \[\begin{align} \label{1} \hat{H} = \hat{H}_\text{s} + \hat{H}_\text{b} + \hat{H}_\text{sb} \end{align}\tag{1}\] where the system Hamiltonian \[\begin{align} \label{2} \hat{H}_\text{s}={\hat{p}^2\over2m} + V(\hat{q}) \end{align}\tag{2}\] is a continuous function of \(\hat{p}\) and \(\hat{q}\), the bath Hamiltonian \(H_\text{b}\) is a collection of harmonic oscillators, and the coupling term \(\hat{H}_\text{sb}=\hat{q}\hat{F}\), where \(\hat{F}\) is the collective bath/solvent coordinate.

For an open quantum system, a non-equilibrium reduced density matrix (obtained by tracing the full density matrix over the locally equilibrated bath states) will equilibrate to the fully entangled equilibrium state after propagating for long enough with the exact quantum time-evolution operator. For a continuous-variable system, this equilibrium state can be written out explicitly in terms of imaginary-time path-integral phase-space variables [8]–[10]. The entanglement appears as a Gaussian term reducing the amplitude of thermal fluctuations around the path centroids in response to the bath coupling. Crucially, the system momenta and positions are correlated through a phase term.

In this article we investigate whether it is possible to equilibrate to this state by propagating stochastic classical trajectories in path-integral phase space. It is well known that this can be done for the position marginal, using the well-established technique of path-integral molecular dynamics, whereby a thermostat is attached to a fictitious ‘ring-polymer’ Hamiltonian [11], [12]. However, generating the full momentum–position distribution from classical trajectories would appear to be impossible on account of the phase, which stochastic classical dynamics cannot be expected to generate. Surprisingly, we show in this article that this is not the case: stochastic classical trajectories can, at least in principle, and sometimes in practice (depending on numerical stability), equilibrate to the exact quantum equilibrium state.

The key is to generate the trajectories using a generalised Langevin equation (GLE) in path-integral phase-space, which was derived recently in Ref. [10] using Matsubara dynamics. ‘Matsubara dynamics’ in this context refers to a semiclassical approximation to the exact quantum dynamics for a continuous system, which arises when initially smooth imaginary-time Feynman paths are assumed to remain smooth for all time [9], [13]. This sole approximation is sufficient to collapse the exact quantum dynamics to an ensemble of classical trajectories in the extended path-integral phase space. Ref. [10] showed that, when applied to a system–bath Hamiltonian of the form of Eq. (1 ), Matsubara dynamics can be used to derive a GLE in the path-integral phase space of the system.

In Ref. [10] it was assumed that, although this GLE conserves the quantum Boltzmann distribution, it would not be capable of taking an arbitrary initial distribution and equilibrating it, on account of the momentum–position phase ¹. Here, we show that it does precisely this, by evolving the system into the complex plane such that the resulting distribution naturally includes the phase and correctly describes the imaginary correlations between the momenta and positions. To prove that the Matsubara GLE has this property, we map it onto an auxiliary space in which the stochastic dynamics becomes Markovian, from which we construct a Fokker–Planck equation, the stationary solution of which is then shown to marginalise to the exact quantum equilibrium state.

The article is structured as follows. Sections 2–4 summarise key results from prior work on classical open systems, and path-integral descriptions of closed and open quantum systems respectively, with Sec. 4.2 reporting a slightly modified version of the Matsubara GLE of Ref. [10]. Most of the new work, including mapping the Matsubara GLE onto an auxiliary space, deriving the Matsubara Fokker–Planck equations, demonstrating that the stationary solutions marginalise to the thermal equilibrium state, and reporting a simple numerical confirmation of these theoretical predictions, is reported in Sec. 5. Sec. 6 concludes the article.

2 Classical background↩︎

2.1 Generalised Langevin equations↩︎

Consider an open classical system with the usual system–bath Hamiltonian [8], [14] \[\begin{gather} \label{eq:caldeira95leggett95classical} H(p,q, \boldsymbol{p},\boldsymbol{x}) = H_\text{s}(p,q)+ \sum_{\alpha=1}^{n_{b}}\Bigg[ \frac{p_{\alpha}^{2}}{2m_{\alpha}} \\ + \frac{1}{2} m_{\alpha} \omega_{\alpha}^{2} \left( x_{\alpha} - \frac{c_{\alpha}}{m_{\alpha}\omega_{\alpha}^{2}} q \right)^{2}\Bigg] \end{gather}\tag{3}\] where \[H_\text{s}(p,q) = \frac{p^{2}}{2m} + V(q)\] and the coefficients \(c_{\alpha}\) are generated from a spectral density \[\label{eq:spectral95density95definition} J(\omega) = \frac{\pi}{2}\sum_{\alpha=1}^{n_{b}}\frac{c_{\alpha}^{2}}{m_{\alpha}\omega_{\alpha}}\delta(\omega-\omega_{\alpha}).\tag{4}\] The dynamics of the system coordinates \((p,q)\) can be propagated using a generalised Langevin equation (GLE) \[\label{eq:GLE95classical} m\ddot{q}(t) = -V'[q(t)] - \int_{0}^{t}\!\mathop{}\!\mathrm{d}s\, \zeta(t-s)\dot{q}(s) + R(t)\tag{5}\] in which \[\label{eq:random95force95classical95randomsampling} R(t) = \sum_{\alpha=1}^{n_{b}}\sqrt{\frac{c_{\alpha}^{2}}{\beta m_{\alpha}\omega_{\alpha}^{2}}}\left[ \lambda_{\alpha}\cos{\omega_{\alpha}t} + \xi_{\alpha}\sin{\omega_{\alpha}t} \right]\tag{6}\] where \(\beta=1/k_\text{B}T\), \(\lambda_{\alpha}\) and \(\xi_{\alpha}\) are Gaussian random variates with zero mean and unit variance, and \[\begin{align} \label{eq:kernelfromJintegral} \zeta(t-s) &= \beta\langle R(s)R(t)\rangle \nonumber\\ &= \sum_{\alpha=1}^{n_{b}}\frac{c_{\alpha}^{2}}{m_{\alpha}\omega_{\alpha}^{2}}\cos{\left(\omega_{\alpha}(t-s)\right)}\nonumber\\ &= \frac{2}{\pi}\int_{0}^{\infty}\!\mathop{}\!\mathrm{d}\omega\,\frac{J(\omega)}{\omega}\cos\left(\omega (t-s)\right). \end{align}\tag{7}\] This equation is the fluctuation–dissipation relation [15] that ensures that an arbitrary, normalised, system probability density evolves to the equilibrium distribution \(Z^{-1}\exp[-\beta H_\text{s}(p,q)]\) as \(t\to\infty\).

In the special case of an ohmic spectral density \[\label{eq:white95spectral95density} J(\omega)=m\gamma\omega\tag{8}\] the noise \(R(t)\) becomes white (and will then be written \(W(t)\)), such that \[\label{eq:kernel95white} \zeta(t-s) = \beta\langle W(s)W(t)\rangle = 2m\gamma \, \delta(t-s)\tag{9}\] and the dynamics becomes Markovian, with Eq. (5 ) reducing to the Langevin equation \[\label{eq:LE95classical} m\ddot q(t) = -V'[q(t)] - m\gamma \dot{q}(t) + W(t).\tag{10}\]

2.2 Fokker–Planck equations for Markovian dynamics↩︎

A stochastic, Markovian, equation of motion such as Eq. (10 ) is equivalent to a deterministic Fokker–Planck equation (FPE) describing the time-evolution of the probability density \(\rho(p,q, t)\). There is a standard procedure for generating the FPE [16], [17]. The stochastic equations of motion are put into the form \[\label{eq:general95eom95continuous} \dot{\mathbf{x}}(t)=\mathbf{G}(\mathbf{x}(t)) + \mathbf{F}(t)\tag{11}\] where \(\mathbf{x}(t)\) is a vector composed of the \(n\) dynamical variables, and \[F_{i}(t)=\sum_{j=1}^{n}\alpha_{ij}\eta_{j}(t)\] in which \(\eta_j(t)\) are white noise-variables with zero mean and unit variance. The FPE corresponding to Eq. (11 ) can then be shown to be [16] \[\label{eq:FP95general} \frac{\partial \rho(\mathbf{x}, t)}{\partial t} = \left[-\frac{\partial}{\partial \mathbf{x}}\cdot\mathbf{G}(\mathbf{x}) + \frac{\partial}{\partial \mathbf{x}}\cdot\left(\mathbf{\Theta}\frac{\partial}{\partial \mathbf{x}}\right)\right]\rho(\mathbf{x}, t)\tag{12}\] where \[\label{eq:theta95matrix95elements} \Theta_{ij} = \frac{1}{2}\sum_{k=1}^{n}\alpha_{ik}\alpha_{jk}.\tag{13}\]

Applying this procedure to Eq. (10 ) yields the well-known Klein–Kramers equation \[\begin{gather} \label{eq:FP95white95classical} \frac{\partial \rho(p,q,t)}{\partial t} = \bigg[-\frac{\partial}{\partial q}\frac{p}{m} + \frac{\partial}{\partial{p}}\frac{\mathop{}\!\mathrm{d}V}{\mathop{}\!\mathrm{d}{q}} \\ + \gamma\frac{\partial}{\partial p}\left(p+\frac{m}{\beta}\frac{\partial}{\partial p}\right)\bigg]\rho(p,q,t). \end{gather}\tag{14}\]

2.3 Fokker–Planck equations for non-Markovian dynamics↩︎

To obtain the FPE corresponding to a non-Markovian stochastic dynamics, we need to map onto an auxiliary space in which the dynamics is Markovian [18]. To illustrate this procedure, consider the GLE of Eq. (5 ) with a Debye–Drude spectral density \[\label{eq:general95spectral95density} J(\omega) = m\omega\frac{\gamma\omega_{c}^{2}}{\omega^{2}+\omega_{c}^{2}}\tag{15}\] for which the kernel is \[\label{eq:ddnoise} \zeta(t-s) = \beta \langle R(s)R(t)\rangle = m\gamma\omega_{c} \mathrm{e}^{-\omega_{c}\left|t-s\right|}.\tag{16}\] The first step is to introduce a stochastic variable \(y(t)\), defined to be the solution to the Ornstein–Uhlenbeck equation \[\label{eq:OU95process95continuous95time} \dot{y}(t) = -\omega_{c} y(t) + \omega_{c} W(t),\tag{17}\] in which \(W(t)\) is the white-noise random variate defined in Eq. (9 ). The solution to this equation for \(t>0\) is \[y(t) = y(0) e^{-\omega_{c} t} + \omega_{c}\int_0^t\!\mathop{}\!\mathrm{d}s\,e^{-\omega_{c}(t-s)} W(s)\] which gives \[\label{eq:OU95process95correlations} \langle y(s)y(t)\rangle = {m\gamma\omega_{c}\over\beta}\mathrm{e}^{-\omega_{c}|t-s|} + \left[ \langle y^{2}(0) \rangle - {m\gamma\omega_{c}\over\beta}\right]\mathrm{e}^{-\omega_{c}(t+s)}.\tag{18}\] When \(y(0)\) is sampled from a Gaussian distribution of variance \(m\gamma\omega_{c}/\beta\), the resulting \(y(t)\) thus gives \(R(t)\) of Eq. (16 ). The second step is to introduce a variable \(v(t)\), defined to be the solution to \[\label{eq:v95eom} \dot{v}(t) = -\omega_{c}v(t) + m\gamma \omega_{c}\dot{q}(t)\tag{19}\] which is \[v(t) = v(0)\mathrm{e}^{-\omega_{c}t} + m\gamma \omega_{c}\int_{0}^{t}\!\mathop{}\!\mathrm{d}s\,\mathrm{e}^{-\omega_{c}(t-s)}\dot{q}(s).\] When \(v(0)\) is set to zero, \(v(t)\) gives the negative of the friction term in Eq. (5 ). Combining these two variables into \[z(t)=y(t)-v(t)\] allows us to write \[\begin{align} \label{eq:marko} \dot{q}(t)&=\frac{p(t)}{m} \nonumber\\ \dot{p}(t) &= -\frac{\mathop{}\!\mathrm{d}V(q(t))}{\mathop{}\!\mathrm{d}q} + z(t) \nonumber\\ \dot{z}(t)&=-\omega_{c}z(t) - \gamma\omega_{c} p(t) + \omega_{c}W(t) \end{align}\tag{20}\] where the dynamics of \((p,q)\) is identical to that generated by the GLE of Eq. (5 ).

Since the dynamics of Eq. (20 ) is Markovian, it is straightforward to obtain \({\boldsymbol{G}}\) and \({\boldsymbol{\Theta}}\) of Eq. (12 ), and thus to construct the FPE in the auxiliary space, which is \[\begin{gather} \label{eq:FP95debye95classical} \frac{\partial \rho(p,q,z,t)}{\partial t} = \bigg[-\frac{\partial}{\partial q} \frac{p}{m}+ \frac{\partial}{\partial{p}}\left(\frac{\mathop{}\!\mathrm{d}V}{\mathop{}\!\mathrm{d}{q}}-z\right) \\ + \gamma\omega_{c}\frac{\partial}{\partial z}p + \omega_{c}\frac{\partial}{\partial z}\left(z+\frac{m\gamma\omega_{c}}{\beta}\frac{\partial}{\partial z}\right)\bigg]\rho(p,q,z,t). \end{gather}\tag{21}\] The stationary solution \[\rho_\text{stat}(p,q,z)={1\over \cal{N}} \exp\left[-\beta\left({p^2\over 2m} +V(q) + {z^2\over 2m\gamma\omega_{c}}\right) \right]\] (where \(\mathcal{N}\) normalises the distribution here and throughout) marginalises to give \(Z^{-1}\exp\left[-\beta\left({p^2/ 2m} +V(q)\right)\right]\) as required.

With the above procedure in hand, it is straightforward to construct the FPE for any spectral density \(J(\omega)\), provided it can be written as a Meier–Tannor expansion [19] over \(N_\text{L}\) Lorentzians. Since each Lorentzian term is equivalent to a Debye–Drude spectral density with complex poles, the auxiliary variables resemble a straightforward \(N_\text{L}\)-dimensional generalisation of the variable \(z(t)\) in Eq. (20 ). For this reason, we will not need to generalise the derivation of Sec. 5 beyond a Debye–Drude spectral density.

3 Imaginary-time Path integrals↩︎

3.1 Equilibrium distributions↩︎

For an isolated system, with Hamiltonian \(\hat{H}_\text{s}\), the position marginal of the path-integral equilibrium distribution is [11], [20] \[\begin{align} \label{eq:pimd} \rho_\text{eq}({\boldsymbol{q}})&={1\over Z}\prod_{i=l}^N \langle q_l|e^{-\beta_N\hat{H}_\text{s}}| q_{l+1} \rangle \nonumber\\ &={1\over \cal{N}}\mathrm{e}^{-\beta_N\left[U_N({\boldsymbol{q}}) + S_N({\boldsymbol{q}})\right]} \end{align}\tag{22}\] where \(\beta_N=\beta/N\), \(q_{N+1}\equiv q_1\), and \[\begin{align} U_N({\boldsymbol{q}}) &=\sum_{l=1}^NV(q_l)\\ S_N({\boldsymbol{q}}) &=\sum_{l=1}^N {(q_{l+1}-q_l)^2m\over 2 \beta_N\hbar^2}. \end{align}\] The imaginary-time Feynman paths \({\boldsymbol{q}}\equiv\left\{q_i\right\}_{i=1}^N\) are Wiener paths in imaginary time \(\tau=0\to\beta\hbar\) which resemble jagged loops, often nicknamed ‘ring polymers’.

The paths can be Fourier-smoothed into differentiable functions of \(\tau\) by approximating them as [21] \[\label{eq:four} q\left(\tau\right)=Q_{0}+\sqrt{2}\sum_{n=1}^{\widetilde{M}}\left[Q_{n}\sin\omega_{n}\tau + Q_{\bar{n}}\cos\omega_{n}\tau\right]\tag{23}\] with \(\widetilde{M} = (M-1)/2\) ², where \(M\) is the number of Fourier modes, and \[\begin{align} \omega_n = {2 n \pi\over \beta\hbar} \end{align}\] are the Matsubara frequencies. Formally, this smoothing is produced by making \(M\) finite, then taking the limit \(N\to\infty\) in the distribution of Eq. (22 ). The mode \(Q_0\) is often referred to as the ‘centroid’ and the modes \(Q_n,n\ne0\) as the ‘Matsubara modes’ [9], [13]. The latter describe the quantum thermal fluctuations around the centroid and shrink to zero in the limit \(\beta\to0\). [Note that in the above and in what follows we use the convention that positive/negative values of \(n\) denote sin/cos modes (where \(\bar n\equiv -n\)) and that these signs are included in the definition of \(\omega_n\); e.g.\(\omega_{\bar{2}}\equiv \omega_{-2}=-4\pi/\beta\hbar\).]

On smoothing the paths, Eq. (22 ) becomes \[\label{eq:rhoQ} \rho_\text{eq}({\boldsymbol{Q}})={1\over \cal{N}}\mathrm{e}^{-\beta\left[U_M({\boldsymbol{Q}}) + S_M({\boldsymbol{Q}})\right]}\tag{24}\] where \[\begin{align}\label{eq:mats95potential} U_{M}(\mathbf{Q}) &= {1\over\beta\hbar}\int_0^{\beta\hbar}\!\mathop{}\!\mathrm{d}\tau\, V[q(\tau)]\\ S_M({\boldsymbol{Q}}) &={m\over 2\beta\hbar} \int_0^{\beta\hbar}\!\mathop{}\!\mathrm{d}\tau\, \left({\partial q(\tau)\over \partial \tau}\right)^2 \nonumber\\ &={m\over 2}\sum_{n=-\widetilde{M}}^{\widetilde{M}}\omega_n^2 Q_n^2. \end{align}\tag{25}\]

One can obtain the fully correlated position–momentum distribution from Eq. (24 ) by inverting the Fourier transforms over momentum that produced the spring term \(S_M({\boldsymbol{Q}})\), to obtain [9] \[\label{eq:rhoco} \rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})={1\over \cal{N}}\mathrm{e}^{-\beta\left[H_M({\boldsymbol{P}},{\boldsymbol{Q}}) -i\theta_{M}(\mathbf{P}, \mathbf{Q}) \right]}\tag{26}\] where \[\begin{align} \tag{27} H_M({\boldsymbol{P}},{\boldsymbol{Q}}) &={1\over\beta\hbar}\int_0^{\beta\hbar}\!d\tau\, \left({p^2(\tau)\over 2m} +V[q(\tau)]\right)\nonumber\\ &={{\boldsymbol{P}}^2\over 2m} + U_M({\boldsymbol{Q}})\\ \tag{28} \theta_{M}(\mathbf{P}, \mathbf{Q}) &= -{1\over\beta\hbar}\int_0^{\beta\hbar}\!d\tau\, {\partial q(\tau)\over \partial \tau}p(\tau) \nonumber\\ &= \sum_{n=-\widetilde{M}}^{\widetilde{M}}\omega_{n}Q_{\bar{n}}P_{n} \end{align}\] and \(p(\tau)\) and \({\boldsymbol{P}}\equiv\{P_n \}\) are defined by analogy with Eq. (23 ). We will refer to \(H_M({\boldsymbol{P}},{\boldsymbol{Q}})\) as the ‘Matsubara Hamiltonian’ and \(\theta_{M}(\mathbf{P}, \mathbf{Q})\) as the ‘Matsubara phase’. It is well known [21], [22] that the distribution \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) yields the exact quantum static expectation values of any operator function of \((\hat{p},\hat{q})\).

3.2 Matsubara dynamics↩︎

Under real-time quantum evolution, a function which initially depends only on the smooth Matsubara paths \(p(\tau)\) and \(q(\tau)\), or equivalently, the Matsubara modes \(({\boldsymbol{P}},{\boldsymbol{Q}})\), rapidly becomes dependent also on the jagged, discontinuous modes that have been excluded from Eq. (23 )—unless \(V(q)\) is harmonic. This is because the exact time evolution is generated by applying \(N\) replicas of the von-Neumann–Liouville operator, one at each of the discrete imaginary time-slices in Eq. (22 ). However, it is possible to project the jagged modes out of this operator, to ensure that an initially smooth function remains smooth for all time [9], [13]. The application of this constraint makes the dynamics classical in the extended phase space \(({\boldsymbol{P}},{\boldsymbol{Q}})\), where the trajectories follow the equations of motion \[\begin{align} \label{eq:matty} \dot{Q}_n & = {\partial H_M({\boldsymbol{P}},{\boldsymbol{Q}})\over \partial P_n}={P_n \over m} \nonumber\\ \dot{P}_n & =-{\partial H_M({\boldsymbol{P}},{\boldsymbol{Q}})\over \partial Q_n}= -{\partial U_M({\boldsymbol{Q}})\over \partial Q_n} \end{align}\tag{29}\] in which \(H_M({\boldsymbol{P}},{\boldsymbol{Q}})\) is the Matsubara Hamiltonian of Eq. (27 ). We will refer to this approximate dynamics as ‘Matsubara dynamics’.

Despite being classical, Matsubara dynamics conserves the quantum Boltzmann distribution of Eq. (26 ). This is because \(H_M({\boldsymbol{P}},{\boldsymbol{Q}})\) is conserved from Eq. (29 ), and \(\theta_M({\boldsymbol{P}},{\boldsymbol{Q}})\) is conserved because it is the momentum conjugate to a continuous symmetry transformation \[\tau\to\tau+\tau_0\] under which \(H_M({\boldsymbol{P}},{\boldsymbol{Q}})\) is invariant. This last property follows from the periodicity of the imaginary-time paths, which allows one to write \[\frac{\partial H_{M}}{\partial\tau_{0}} = 0.\] Additional properties of this imaginary-time-translation symmetry are given in Appendix 7. Two key relations that will be used later in Sec. 5 are \[\begin{align} \label{vint} \frac{\partial U_{M}}{\partial\tau_{0}} &=- \sum_{n=-\widetilde{M}}^{\widetilde{M}} \omega_{n}Q_{\bar{n}} \frac{\partial U_{M}}{\partial Q_{n}}=0 \end{align}\tag{30}\] and \[\begin{align} \label{eq:mats95mode95symmetry951} R_{n}\left(\tau_{0} + \frac{\beta\hbar}{4n}\right) &= -R_{\bar{n}}\left(\tau_{0}\right)\nonumber\\ R_{\bar{n}}\left(\tau_{0} + \frac{\beta\hbar}{4n}\right) &= R_{n}\left(\tau_{0}\right) \end{align}\tag{31}\] where \(R\) denotes \(P\) or \(Q\).

4 Path-integral description of open quantum systems↩︎

The results given in the previous section generalise straightforwardly to the full Caldeira–Leggett Hamiltonian \(\hat{H}\) of Eq. (1 ). Here we summarise the reduced picture in the system coordinates, obtained by tracing over the bath modes.

4.1 Reduced equilibrium distribution↩︎

The reduced phase-space equilibrium distribution is [10], [23] \[\label{eq:redeq} \rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})={1\over \cal{N}}\mathrm{e}^{-\beta\left[F_M({\boldsymbol{P}},{\boldsymbol{Q}}) -i\theta_{M}(\mathbf{P}, \mathbf{Q}) \right]}\tag{32}\] where \(({\boldsymbol{P}},{\boldsymbol{Q}})\) are the Matsubara modes of the system, \[\label{eq:fm} F_M({\boldsymbol{P}},{\boldsymbol{Q}}) =H_M({\boldsymbol{P}},{\boldsymbol{Q}}) + \frac{m}{2}\sum_{n=-\widetilde{M}}^{\widetilde{M}}\left|\omega_{n}\right| \widetilde{\zeta}(\left|\omega_{n}\right|)Q_{n}^{2}\tag{33}\] and \(H_M({\boldsymbol{P}},{\boldsymbol{Q}})\) is given by Eq. (27 ), with \(V(q)\) taken to be the system potential of Eq. (2 ). The phase \(\theta_{M}(\mathbf{P}, \mathbf{Q})\) is now the system contribution to the total Matsubara phase, the bath contribution having been integrated out to give the second term in Eq. (33 ) which describes the system–bath entanglement. The term \(\widetilde{\zeta}(s)\) is given by \[\begin{align} \label{eq:zeta95hat} \widetilde{\zeta}(s) &= \sum_{\alpha=1}^{n_{b}}\frac{c_{\alpha}^{2}}{ m_{\alpha}\omega_{\alpha}^{2}}\frac{s}{s^{2}+\omega_{\alpha}^{2}}\nonumber\\ &= \int_0^\infty\!\mathop{}\!\mathrm{d}t\,\zeta(t) e^{-st}. \end{align}\tag{34}\] For the white-noise spectral density of Eq. (8 ), \(\widetilde{\zeta}(|\omega_n|)=\gamma|\omega_n|\), which illustrates clearly that increasing the system–bath coupling strength \(\gamma\) decreases the variance of the quantum thermal fluctuations around the centroid; this trend also holds for other spectral densities.

4.2 Matsubara generalised Langevin equation↩︎

Reference [10] reported a GLE which was derived by starting from the Matsubara dynamics of the system plus bath, then eliminating the bath modes. This derivation assumed a direct product initial condition. Appendix 8 of this article reports a modified derivation, which assumes an initial distribution in which the bath is locally equilibrated around the system coordinate. The resulting GLE is \[\begin{gather} \label{eq:GLE95matsubara} m\ddot{Q}_{n}(t) = f^{(n)}_{M}(\mathbf{Q}(t)) - \int_{0}^{t}\!\mathop{}\!\mathrm{d}s\, \zeta(t-s)\dot{Q}_{n}(s) + R_{n}^{\mathbb{C}}(t) \\ - \left[Q_{n}(0)K_{n}(t) - iQ_{\bar{n}}(0)L_{n}(t) \right] \end{gather}\tag{35}\] where \[\begin{align} \label{eq:force} f^{(n)}_{M}(\mathbf{Q}) = -{\partial U_M(\mathbf{Q})\over \partial Q_n} \end{align}\tag{36}\] and the complex-valued noise is \[\begin{gather} \label{eq:complex95noise} R_{n}^{\mathbb{C}}(t) = \sum_{\alpha=1}^{n_{b}}\sqrt{\frac{c_{\alpha}^{2}}{\beta m_{\alpha}\omega_{\alpha}^{2}}}\bigg[ \frac{\omega_{\alpha}}{\omega_{\alpha n}}\lambda_{\alpha n}\cos{\omega_{\alpha}t} \\ + \left(\xi_{\alpha n} + i\frac{\omega_{n}}{\omega_{\alpha n}}\lambda_{\alpha \bar{n}} \right) \sin{\omega_{\alpha}t} \bigg] \end{gather}\tag{37}\] where \(\lambda_{\alpha n}\) and \(\xi_{\alpha n}\) are Gaussian random variates with zero mean and unit variance and \(\omega_{\alpha n}^{2}=\omega_{\alpha}^{2}+\omega_{n}^{2}\). The memory kernels of \(R_{n}^{\mathbb{C}}(t)\) are \[\begin{align} \nonumber \langle R_{n}^{\mathbb{C}}(s)R_{n}^{\mathbb{C}}(t)\rangle &= \frac{\zeta(t-s)}{\beta}-\frac{K_{n}(t-s)}{\beta}\\ \nonumber \langle R_{n}^{\mathbb{C}}(s)R_{\bar{n}}^{\mathbb{C}}(t)\rangle &= -\frac{iL_{n}(t-s)}{\beta} \\ \langle R_{n}^{\mathbb{C}}(s)R_{n^{\prime}}^{\mathbb{C}}(t)\rangle &= 0,\quad\quad |n|\ne|n'| \end{align}\] where \[\tag{38} \begin{align} K_{n}(t) &= \frac{2\omega_{n}^{2}}{\pi}\int_{0}^{\infty}\!\mathop{}\!\mathrm{d}\omega\,\frac{J(\omega)}{\omega}\frac{\cos\omega t}{\omega^{2}+\omega_{n}^{2}} \tag{39}\\ L_{n}(t) &= \frac{2\omega_{n}}{\pi}\int_{0}^{\infty}\!\mathop{}\!\mathrm{d}\omega\,J(\omega)\frac{\sin\omega t}{\omega^{2}+\omega_{n}^{2}}. \end{align}\]

Equation (35 ) differs from the GLE of Ref. [10] only in the form of the transient driving terms in the second line. These terms result from locally equilibrating the bath modes around the system coordinate at time \(t=0\) and cannot be removed by adjusting the initial distribution of the bath modes. They arise because the non-Markovian \(R_{n}^{\mathbb{C}}(t)\) lacks history at \(t=0\). As time advances, the noise gradually builds up a history, and these terms consequently decay to zero. They therefore have no effect on the ability of Eq. (35 ) to equilibrate the system in the long time limit, and in what follows shall be neglected.

5 Equilibration properties of the Matsubara GLE↩︎

This Section focuses on the main goal of the article, which is to show that the GLE of Eq. (35 ) is capable of equilibrating an arbitrary initial system phase-space distribution \(\rho_0({\boldsymbol{P}},{\boldsymbol{Q}})\). We first consider the special case of a white-noise spectral density in Sec. 5.1, then extend to the Debye–Drude spectral density in Sec. 5.2. As mentioned above, the latter is sufficiently general to imply that the GLE will also equilibrate any other spectral density that can be expanded as a sum of Lorentzians.

5.1 White-noise spectral density↩︎

For a white-noise spectral density, neglecting the transient terms, Eq. (35 ) simplifies to \[\begin{align} \label{eq:LE95matsubara} m\ddot{Q}_{n}(t) &= f^{(n)}_{M}(\mathbf{Q}(t)) - m\gamma\dot{Q}_{n}(s) + R_{n}^{\mathbb{C}}(t) \end{align}\tag{40}\] and the non-zero components of the noise kernel become \[\begin{align} \label{eq:mats95fd95relations95white} \langle R_{n}^{\mathbb{C}}(s)R_{n}^{\mathbb{C}}(t)\rangle &= \frac{2m\gamma}{\beta}\,\delta(t-s)-\frac{m\gamma|\omega_{n}|}{\beta}\mathrm{e}^{-|\omega_{n}(t-s)|}\nonumber\\ \langle R_{n}^{\mathbb{C}}(s)R_{\bar{n}}^{\mathbb{C}}(t)\rangle &= -\frac{im\gamma\omega_{n}}{\beta}\,\text{sgn}(t-s)\mathrm{e}^{-|\omega_{n}(t-s)|}. \end{align}\tag{41}\] Superficially Eq. (40 ) resembles a Langevin equation. However, the dynamics it describes is non-Markovian owing to the correlations in the noise kernel. Nevertheless, by introducing a collection of auxiliary variables each undergoing an Ornstein–Uhlenbeck process to treat the noise, the dynamics can be mapped onto a Markovian stochastic process.

5.1.1 Markovian auxiliary space and Fokker–Planck equation↩︎

We introduce the auxiliary variables by expanding \(R_{n}^{\mathbb{C}}(t)\) as \[\label{eq:white95undetermined95noise} R_{n}^{\mathbb{C}}(t) = a_{n}W_{n}(t) + b_{n}Z_{n}(t) + c_{n}W_{\bar{n}}(t) + d_{n}Z_{\bar{n}}(t)\tag{42}\] where the \(W_{n}(t)\) are independent white noise variates with the same kernels as Eq. (9 ), and the \(Z_n(t)\) are solutions to the Ornstein–Uhlenbeck equation \[\label{eq:Z95eom} \dot{Z}_{n}(t)=-|\omega_{n}|Z_{n}(t)+|\omega_{n}|W_{n}(t).\tag{43}\] The complex-valued coefficients \(\left\{a_n,b_n,c_n,d_n \right\}\) satisfy a set of constraints. First, since \(R_{n}^{\mathbb{C}}\), \(W_n\) and \(Z_n\) are Matsubara variables they must satisfy the imaginary-time translation relations of Eq. (31 ), which requires that \[\label{eq:noise95coeffs95negative95relations} a_{n}=a_{\bar{n}},\quad b_{n}=b_{\bar{n}},\quad c_{n}=-c_{\bar{n}},\quad d_{n}=-d_{\bar{n}}.\tag{44}\] Second, the kernels of Eq. (42 ) ³, \[\begin{align} \langle R_{n}^{\mathbb{C}}(s)R_{n}^{\mathbb{C}}(t)\rangle &= \frac{2m\gamma}{\beta}\,\delta(t-s)\left[a_{n}^{2}+c_{n}^{2}\right] \nonumber\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! +\frac{m\gamma|\omega_{n}|}{\beta}\mathrm{e}^{|\omega_{n}(t-s)|}\left[b_{n}^{2} + d_{n}^{2} + 2\left(a_{n}b_{n}+c_{n}d_{n}\right)\right]\nonumber\\ \langle R_{n}^{\mathbb{C}}(s)R_{\bar{n}}^{\mathbb{C}}(t)\rangle &= - \frac{m\gamma|\omega_{n}|}{\beta}\,\text{sgn}(t-s)\mathrm{e}^{|\omega_{n}(t-s)|} \nonumber\\ &\quad\quad\quad\quad\times\left[2\left(a_{n}d_{n}-b_{n}c_{n}\right)\right] \end{align}\] must match those of Eq. (41 ), which requires that \[\begin{align} \label{eq:noise95coeffs95relations95white} a_{n}^{2}+c_{n}^{2}&=1\nonumber\\ b_{n}^{2} + d_{n}^{2} + 2\left(a_{n}b_{n}+c_{n}d_{n}\right) &= -1\nonumber\\ a_{n}d_{n}-b_{n}c_{n} &= \frac{i\,\text{sgn}(n)}{2}. \end{align}\tag{45}\] These constraints, together with Eq. (44 ), are not sufficient to specify the coefficients uniquely, so we have some freedom in how they are chosen. A convenient choice, used to obtain the numerical results in Sec. 5.1.4, is \[\begin{align} a_n = 1,\quad b_n =-\frac{3}{2},\quad c_n = 0,\quad d_n = \frac{i\,\text{sgn}(n)}{2}. \end{align}\]

Using Eq. (42 ), we can now rewrite the dynamics of Eq. (40 ) as \[\begin{align} \label{eq:markovian95LE} \dot{Q}_{n}(t) &= \frac{P_{n}(t)}{m}\nonumber\\ \dot{P}_{n}(t) &= -\frac{\partial U_{M}(\mathbf{Q}(t))}{\partial Q_{n}}-\gamma P_{n}(t) + \chi_{n}(t) + a_{n}W_{n}(t) \nonumber\\ &\quad\quad\quad\quad\quad+ c_{n}W_{\bar{n}}(t)\nonumber\\ \dot{\chi}_{n}(t) &= -|\omega_{n}|\chi_{n}(t) + |\omega_{n}|\left[b_{n}W_{n}(t)+d_{n}W_{\bar{n}}(t)\right] \end{align}\tag{46}\] where we have defined \[\label{eq:combining95Zs} \chi_{n}(t) = b_{n}Z_{n}(t)+d_{n}Z_{\bar{n}}(t).\tag{47}\] These equations of motion are Markovian and of the form of Eq. (11 ). We can therefore construct the matrix \[\label{eq:white95noise95matrix} \mathbf{\Theta} = \frac{m\gamma}{\beta}\left(\begin{matrix} \mathbf{I} && \mathbf{0} &&\boldsymbol{\kappa} \\ \mathbf{0} && \mathbf{0} && \mathbf{0} \\ \boldsymbol{\kappa}^{T} && \mathbf{0} && \boldsymbol{\iota} \end{matrix}\right)\tag{48}\] indexed by \(({\boldsymbol{P}},{\boldsymbol{Q}} ,{\boldsymbol{\chi}})\), in which the \(M\)-dimensional submatrix elements are \[\begin{align} {{\iota}}_{nn'} &= \left(b_{n}^{2}+d_{n}^{2}\right)\omega_{n}^{2} \delta_{nn'}\nonumber\\ {{\kappa}}_{nn'} &=(a_{n}b_{n}+c_{n}d_{n})|\omega_{n}|\delta_{nn'}-i{\omega_{n}\over2}\delta_{\bar{n} n'} \end{align}\] from which we obtain the Fokker–Planck equation \[\begin{gather} \label{eq:fp95mats95white} \frac{\partial \rho(\mathbf{P}, \mathbf{Q}, \boldsymbol{\chi}, t)}{\partial t} = \sum_{n=-\widetilde{M}}^{\widetilde{M}}\bigg[-\frac{\partial}{\partial Q_{n}}\frac{P_{n}}{m} \\ + \frac{\partial}{\partial P_{n}}\left(\frac{\partial U_{M}(\mathbf{Q})}{\partial Q_{n}} - \chi_{n}\right) + \gamma \frac{\partial}{\partial P_{n}}\left(P_{n}+\frac{m}{\beta}\frac{\partial}{\partial P_{n}}\right) \\ + |\omega_{n}| \frac{\partial}{\partial \chi_{n}}\left(\chi_{n}+\frac{m\gamma|\omega_{n}|}{\beta}\left(b_{n}^{2}+d_{n}^{2}\right)\frac{\partial}{\partial \chi_{n}}\right) \\ + \frac{2m\gamma|\omega_{n}|}{\beta}\left(a_{n}b_{n}+c_{n}d_{n}\right) \frac{\partial^{2}}{\partial P_{n}\partial \chi_{n}} \\- \frac{i m \gamma\omega_{n}}{\beta}\frac{\partial^{2}}{\partial P_{n}\partial\chi_{\bar{n}}} \bigg]\rho(\mathbf{P}, \mathbf{Q}, \boldsymbol{\chi}, t). \end{gather}\tag{49}\]

5.1.2 Form of the FPE stationary solution↩︎

To demonstrate that Eq. (40 ) will equilibrate any initial distribution \(\rho_0({\boldsymbol{P}},{\boldsymbol{Q}})\) given sufficient time, we need to find the stationary solution of the FPE in Eq. (49 ) and show that it marginalises to \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) of Eq. (32 ).

Finding stationary solutions to FPEs is in general difficult [17]. However, Eq. (49 ) turns out to be readily solvable because one can first compute the solution in the special case that \(V(q)\) is harmonic. This is easily done using the method of moments [17], as described in Ref. [23]. One then finds that the harmonic \(V(q)\) reappears in the stationary solution and replacing it by the general case ⁴ gives \[\label{eq:mats95dist95harmonic95white95full} \rho_\text{stat}(\mathbf{P},\mathbf{Q},\boldsymbol{\chi}) ={1\over \cal{N}} \mathrm{e}^{-\beta\Phi(\mathbf{P},\mathbf{Q},\boldsymbol{\chi})} \delta\left(\chi_{0}\right)\tag{50}\] in which \[\begin{gather} \label{eq:mats95dist95exponent95general95white95full} \Phi(\mathbf{P},\mathbf{Q},\boldsymbol{\chi}) = F_{M}(\mathbf{P},\mathbf{Q}) -i\theta_{M}(\mathbf{P},\mathbf{Q}) \\ +\sum_{n=-\widetilde{M}}^{\widetilde{M}}\!\!\!\!\!\!\!\!{\phantom{\frac{x}{y}}}^{\prime}\,\,\, {1\over b_{n}^{2}+d_{n}^{2}}\left[\frac{\chi_{n}^{2}}{2m\gamma|\omega_{n}|} + \chi_{n}\left(Q_{n}-i\,\text{sgn}(n)Q_{\bar{n}}\right) \right] \end{gather}\tag{51}\] where \(F_{M}(\mathbf{P},\mathbf{Q})\) is defined in Eq. (33 ) and the prime indicates that \(n=0\) is omitted from the sum. Substituting into Eq. (49 ), we obtain \[\begin{align} \label{eq:drho95dt} \frac{\partial \rho_\text{stat}}{\partial t} &= i\beta\rho_\text{stat}\sum_{n=-\widetilde{M}}^{\widetilde{M}}\omega_{n}Q_{\bar{n}}\frac{\partial U_{M}}{\partial Q_{n}}\nonumber\\ &=0 \end{align}\tag{52}\] where the second equality follows from the imaginary-time translation symmetry of \(U_{M}\)—see Eq. (30 ). Thus \(\rho_{\mathrm{stat}}\) is the general form of the stationary solution to the FPE of Eq. (49 ). On marginalising over \(\boldsymbol{\chi}\), we obtain \[\int\!d\boldsymbol{\chi}\, \rho_\text{stat}(\mathbf{P},\mathbf{Q},\boldsymbol{\chi})=\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}}).\]

We have therefore shown that the white-noise Matsubara GLE of Eq. (40 ) is capable of equilibrating an arbitrary initial distribution of system phase-space points to the distribution \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) of Eq. (32 ).

5.1.3 The imaginary momentum–position correlation↩︎

The Matsubara Langevin equation of Eq. (40 ) produces the equilibrium distribution \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) without explicitly invoking the phase \(\exp[-i\theta_M({\boldsymbol{P}},{\boldsymbol{Q}})]\). In other words, an expectation value of some property \(A({\boldsymbol{P}},{\boldsymbol{Q}})\) can be evaluated using \[\begin{gather} \label{eq:langap} \int\! d{\boldsymbol{P}}\int\! d{\boldsymbol{Q}}\, \rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}}) A({\boldsymbol{P}},{\boldsymbol{Q}})=\\\lim_{N_\text{t}\to\infty}{1\over N_\text{t}}\sum_{i=1}^{N_\text{t}} A({\boldsymbol{P}}(t_i),{\boldsymbol{Q}}(t_i)) \end{gather}\tag{53}\] where the phase-space points \(\left\{{\boldsymbol{P}}(t_i),{\boldsymbol{Q}}(t_i) \right\}\) are taken at \(N_\text{t}\) times \(t_i\) from a trajectory propagated using Eq. (40 ). The phase is included naturally in the stationary distribution as a result of the complex-valued random force \(R_{n}^{\mathbb{C}}\) which pushes the phase-space variables away from the real axis, such that \(P_n\) and \(Q_{\bar{n}}\) maintain the purely imaginary correlation given by \(\theta_M({\boldsymbol{P}},{\boldsymbol{Q}})\).

Another way of looking at this is that the propagation of Eq. (40 ) has the effect of analytically continuing the distribution \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) into the complex plane, such that the phase no longer appears. It therefore eliminates the dreaded ‘phase problem’ but has replaced it with another numerical difficulty, since classical dynamics in the complex plane is well known to be numerically unstable. From Eq. (46 ), we can see that the instability is likely to be tamed by increasing \(\gamma\) and worsened by increasing \(|n|\). This suggests that any attempt to propagate trajectories using Eq. (40 ) will be numerically stable only if \(M\) is sufficiently small.

5.1.4 A simple numerical test↩︎

We have investigated these numerical properties for the simple case of a quartic system potential \[\label{eq:quartic95potential} V(q)=\frac{q^{4}}{4}\tag{54}\] with a moderate damping strength of \(\gamma=1\). At a temperature of \(\beta=50\), we find that the dynamics is stable for \(M=13\), which is far lower than the value of \(M\simeq 100\) needed to converge \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) at this temperature. Nonetheless, this small value of \(M\) is sufficient to confirm numerically that Eq. (40 ) is capable of generating the imaginary momentum–position correlations in \(\rho_\text{eq}(\mathbf{P},\mathbf{Q})\). A total of \(2\times 10^{5}\) phase-space points were sampled from an initial distribution \[\begin{gather} \label{eq:mats95rp95dist95white} \rho_{0}(\mathbf{P},\mathbf{Q}) ={1\over\cal{N}} \exp\Bigg(-\beta\Bigg[ U_{M}(\mathbf{Q}) + \sum_{n=-\widetilde{M}}^{\widetilde{M}} \bigg[ \frac{P_{n}^{2}}{2m} \\ + \frac{1}{2}m\left(\omega_{n}^{2} + \gamma|\omega_{n}|\right)Q_{n}^{2}\bigg]\Bigg]\Bigg) \end{gather}\tag{55}\] where \(m=1\) and \(\hbar=1\). The points were propagated until \(t=20\) (reduced units) by integrating Eq. (46 ) using a modified velocity Verlet algorithm with time-step \(\Delta t= 0.005\). Figure 1 plots various moments of \(\rho(\mathbf{Q},\mathbf{P},t)\) against time, with \(\langle P_nQ_{\bar{n}} \rangle\) illustrating clearly the growth of the imaginary momentum–position correlation between \(P_n\) and \(Q_{\bar{n}}\).

5.2 Debye–Drude spectral density↩︎

For the Debye–Drude spectral density, the derivation follows the same steps as those just given for white noise, with the addition of extra auxiliary variables to account for the non-Markovian decay of \(\zeta(t)\).

The Debye–Drude noise kernels, obtained by substituting \(J(\omega)\) of Eq. (15 ) into Eq. (38 ), are \[\begin{align} \label{eq:mats95fd95relations95nn95general} \langle R_{n}^{\mathbb{C}}(s)R_{n}^{\mathbb{C}}(t)\rangle &= \frac{m}{\beta}\frac{\gamma\omega_{c}^{2}}{\omega_{c}^{2}-\omega_{n}^{2}}\nonumber\\ &\quad\quad\times\left(\omega_{c}\mathrm{e}^{-\omega_{c}|t-s|}-|\omega_{n}|\mathrm{e}^{-|\omega_{n}(t-s)|}\right)\nonumber\\ \langle R_{n}^{\mathbb{C}}(s)R_{\bar{n}}^{\mathbb{C}}(t)\rangle &= \frac{im\omega_{n}}{\beta}\,\text{sgn}(t-s)\frac{\gamma\omega_{c}^{2}}{\omega_{c}^{2}-\omega_{n}^{2}}\nonumber\\ &\quad\quad\times\left(\mathrm{e}^{-\omega_{c}|t-s|}-\mathrm{e}^{-|\omega_{n}(t-s)|}\right) \end{align}\tag{56}\] and can be reproduced by expanding the noise as \[\begin{gather} \label{eq:general95undetermined95noise} R_{n}^{\mathbb{C}}(t) = a_{n}Y_{n}(t) + b_{n}Z_{n}(t) + c_{n}Y_{\bar{n}}(t) + d_{n}Z_{\bar{n}}(t) \end{gather}\tag{57}\] where the \(Z_n(t)\) are the solutions of Eq. (43 ), and the \(Y_n(t)\) (analogous to \(y(t)\) of Sec. 2.3) are the solutions of \[\label{eq:Yi95eom} \dot{Y}_{n}(t)=-\omega_{c}Y_{n}(t)+\omega_{c}W_{n}(t).\tag{58}\] The coefficients \(\left\{a_n,b_n,c_n,d_n \right\}\) satisfy the symmetry constraints given in Eq. (44 ). Additional constraints are found by matching the kernels to Eq. (56 ) to obtain ⁵ \[\begin{align} a_{n}^{2}+c_{n}^{2} + \frac{2|\omega_{n}|}{\omega_{c}+|\omega_{n}|}\left(a_{n}b_{n} + c_{n}d_{n}\right)&=\frac{\omega_{c}^{2}}{\omega_{c}^{2}-\omega_{n}^{2}}\nonumber\\ b_{n}^{2} + d_{n}^{2} + \frac{2\omega_{c}}{\omega_{c}+|\omega_{n}|}\left(a_{n}b_{n}+c_{n}d_{n}\right) &= -\frac{\omega_{c}^{2}}{\omega_{c}^{2}-\omega_{n}^{2}}\nonumber\\ a_{n}d_{n}-b_{n}c_{n} = \frac{i\,\omega_{c}\,\text{sgn}(n)}{2\left(\omega_{c}-|\omega_{n}|\right)}&. \end{align}\] To deal with the non-Markovian friction term in the GLE, we introduce the variables \(V_{n}(t)\) (analogous to \(v(t)\) of Sec. 2.3) which satisfy \[\dot{V}_{n}(t)=-\omega_{c}V_{n}(t) + m\gamma \omega_{c} \dot{Q}_{n}(t)\] with initial conditions \(V_{n}(0)=0\). As a result, \[V_{n}(t) = m\gamma\omega_{c}\int_{0}^{t}\!\mathop{}\!\mathrm{d}s\,\mathrm{e}^{-\omega_{c}(t-s)}\dot{Q}_{n}(s)\] which regenerates the corresponding friction term in Eq. (35 ) (for the Debye–Drude spectral density). Taking the linear combination \[\label{eq:combining95V95Ys} \psi_{n}(t)=-V_{n}(t)+a_{n}Y_{n}(t)+c_{n}Y_{\bar{n}}(t)\tag{59}\]

Substituting \(\rho_{\mathrm{stat}}\) into Eq. (60 ) then yields the expression given on the right-hand side of Eq. (52 ), thus confirming that \(\rho_{\mathrm{stat}}\) is the stationary solution. Finally, marginalising over \((\boldsymbol{\psi}, \boldsymbol{\chi})\) gives \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) of Eq. (32 ), with \[\tilde{\zeta}(|\omega_{n}|) = m\frac{ \gamma\omega_{c}}{\omega_{c} + |\omega_{n}|}\] which is the expression for \(\tilde{\zeta}(|\omega_{n}|)\) obtained by substituting the Debye–Drude spectral density [Eq. (15 )] into Eq. (34 ).

We have thus shown that the Matsubara GLE of Eq. (35 ) equilibrates to the thermal equlibrium state \(\rho_\text{eq}({\boldsymbol{P}},{\boldsymbol{Q}})\) for a Debye–Drude spectral density. As mentioned above, this is sufficiently general to cover any spectral density that can be expanded as a sum of Lorentzians.

6 Conclusions↩︎

It is surprising that stochastic classical trajectories in an extended space can equilibrate to the exact quantum equilibrium state, correctly including the purely imaginary correlations between the momentum and position variables. Nonetheless, there are no free lunches, and the stochastic classical dynamics produces these correlations by evolving the trajectories into the complex plane. Inevitably, this leads to numerical instability, so we are not advocating the use of this approach as a practical method. It may, however, be useful as a starting point from which to derive more approximate and practical methods for simulating continuous-variable open quantum dynamics.

Data availability↩︎

The data that support the findings of this article are publicly available in Ref. [24].

7 Imaginary-time symmetries↩︎

The effect of an imaginary-time translation on a variable \(r(\tau)\) (where \(r\) is \(p\) or \(q\) and \(R_{n}\equiv R_{n}\left(\tau_{0}=0\right)\)) is \[\begin{gather} r(\tau_{0}+\tau) = R_{0}+\sqrt{2}\sum_{n=1}^{\tilde{M}}\big[R_{n}\sin\left(\omega_{n}\left(\tau_{0}+\tau\right)\right) \\ + R_{\bar{n}}\cos\left(\omega_{n}\left(\tau_{0}+\tau\right)\right)\big] \end{gather}\] which can be rearranged into \[\begin{gather} r(\tau_{0}+\tau) = R_{0}(\tau_{0})+\sqrt{2}\sum_{n=1}^{\tilde{M}}\big[R_{n}(\tau_{0})\sin\omega_{n}\tau \\ + R_{\bar{n}}(\tau_{0})\cos\omega_{n}\tau\big] \end{gather}\] where \[\begin{align} R_{n}(\tau_{0}) &= R_{n}\cos\omega_{n}\tau_{0} - R_{\bar{n}}\sin\omega_{n}\tau_{0}\label{eq:Qn95imag95time}\\ R_{\bar{n}}(\tau_{0}) &= R_{n}\sin\omega_{n}\tau_{0} + R_{\bar{n}}\cos\omega_{n}\tau_{0}. \end{align}\tag{63}\] Since \(\omega_{n}=2n\pi/\beta\hbar\), it follows that \[\begin{align} R_{n}(\tau_{0} + \beta\hbar/4n) &= -R_{n}\sin\omega_{n}\tau_{0} - R_{\bar{n}}\cos\omega_{n}\tau_{0} \nonumber\\ R_{\bar{n}}(\tau_{0} + \beta\hbar/4n) &= R_{n}\cos\omega_{n}\tau_{0} - R_{\bar{n}}\sin\omega_{n}\tau_{0} \end{align}\] from which we obtain Eq. (31 ).

Equation (63 ) also gives \[\label{eq:Xn95imag95time95derivative} \frac{\partial R_{n}(\tau_{0})}{\partial\tau_{0}} = -\omega_{n}R_{n}\sin\omega_{n}\tau_{0} - \omega_{n}R_{\bar{n}}\cos\omega_{n}\tau_{0}\tag{64}\] so that \[\frac{\partial R_{n}}{\partial\tau_{0}} \equiv \frac{\partial R_{n}(0)}{\partial\tau_{0}} = - \omega_{n}R_{\bar{n}}.\] Subsituting this last expression (with \(R\to Q\)) into \[{\partial U_M({\boldsymbol{Q}})\over\partial \tau_0}= 0 = \sum_{n=-\widetilde{M}}^{\widetilde{M}} {\partial Q_n\over\partial \tau_0}{\partial U_M({\boldsymbol{Q}})\over\partial Q_n}\] gives Eq. (30 ).

8 Derivation of the Matsubara GLE↩︎

The Matsubara Hamiltonian corresponding to Eq. (3 ) is \[\begin{gather} H_M^\text{tot}(\mathbf{P},\mathbf{Q},{\boldsymbol{P}}_\text{b},{\boldsymbol{X}})=H_M(\mathbf{P},\mathbf{Q}) + \sum_{\alpha=1}^{n_{b}}\sum_{n=-\widetilde{M}}^{\widetilde{M}}\Bigg[ \frac{P_{\alpha,n}^{2}}{2m_{\alpha}} \\ + \frac{1}{2} m_{\alpha} \omega_{\alpha}^{2} \left( X_{\alpha,n} - \frac{c_{\alpha}}{m_{\alpha}\omega_{\alpha}^{2}} Q_n \right)^{2}\Bigg] \end{gather}\] where \((P_{\alpha,n} ,X_{\alpha,n})\) are the Matsubara modes corresponding to bath variables \((p_\alpha,x_\alpha)\). To obtain an initial distribution in which the bath is locally equilibrated around the system variables, and the system is prepared in some arbitrary initial density \(\rho_0\), we start with the full system–bath equilibrium distribution \[\rho_M({\boldsymbol{P}},{\boldsymbol{Q}},{\boldsymbol{P}}_\text{b},\mathbf{X})={1\over \cal{N}}\mathrm{e}^{-\beta\left[H_M^\text{tot}({\boldsymbol{P}},{\boldsymbol{Q}},{\boldsymbol{P}}_\text{b},\mathbf{X}) -i\theta^\text{tot}_{M}(\mathbf{P}, \mathbf{Q},{\boldsymbol{P}}_\text{b},\mathbf{X}) \right]}\] where \[\begin{gather} \theta^\text{tot}_{M}(\mathbf{P}, \mathbf{Q},{\boldsymbol{P}}_\text{b},\mathbf{X})= \theta_{M}(\mathbf{P}, \mathbf{Q}) \\+ \sum_{\alpha=1}^{n_{b}}\sum_{n=-\widetilde{M}}^{\widetilde{M}}\omega_{n}X_{\alpha,{\bar{n}}}P_{\alpha,n}. \end{gather}\] We then make the substitution \[P_{\alpha,n}=\Pi_{\alpha, n} + im_{\alpha}\omega_{n}X_{\alpha \bar{n}}\] which is equivalent to shifting the momentum integration contour so as to eliminate the bath contribution to the phase, giving \[\label{eq:equil95dist95bath95AC} \rho_M({\boldsymbol{P}},{\boldsymbol{Q}},\mathbf{\Pi}_{b},\mathbf{X})={1\over \cal{N}}\mathrm{e}^{-\beta\left[F_M({\boldsymbol{P}},{\boldsymbol{Q}})+ B_M({\boldsymbol{\Pi}}_{b},{\boldsymbol{X}};{\boldsymbol{Q}}) -i\theta_{M}(\mathbf{P}, \mathbf{Q}) \right]}\tag{65}\] where \[B_M({\boldsymbol{\Pi}}_{b},{\boldsymbol{X}};{\boldsymbol{Q}}) = \sum_{\alpha=1}^{n_{b}}\sum_{n=-\widetilde{M}}^{\widetilde{M}} \left(\frac{\Pi_{\alpha,n}^{2}}{2m_{\alpha}} + \frac{1}{2}m_{\alpha}\omega_{\alpha n}^{2}{\overline{X}}_{\alpha,n}^2 \right)\] and \[{\overline{X}}_{\alpha,n}\equiv X_{\alpha,n} - \frac{c_{\alpha}}{m_{\alpha}\omega_{\alpha n}^{2}}Q_{n}.\] Replacing the system-dependent part of Eq. (65 ) by \(\rho_0({\boldsymbol{P}},{\boldsymbol{Q}})\), we then obtain \[\label{eq:app95qb95dist95bath95AC} \rho_M({\boldsymbol{P}},{\boldsymbol{Q}},\mathbf{\Pi}_{b},\mathbf{X})={\rho_0({\boldsymbol{P}},{\boldsymbol{Q}})\over \cal{N}}\mathrm{e}^{-\beta B_M({\boldsymbol{\Pi}}_{b},{\boldsymbol{X}};{\boldsymbol{Q}})}\tag{66}\] which describes an initial distribution in which the bath variables \(X_{\alpha,n}\) are equilibrated about the system coordinate \(Q_n\), and the system distribution is arbitrary.

From Eq. (29 ), \(H_M^\text{tot}(\mathbf{P},\mathbf{Q},{\boldsymbol{P}}_\text{b},{\boldsymbol{X}})\) generates the Matsubara-dynamics equations of motion \[\begin{align} \label{eq:app95eoms} m\ddot{Q}_{n} &= f^{(n)}_{M}(\mathbf{Q}) + \sum_{\alpha=1}^{n_{b}}c_{\alpha}\left(X_{\alpha n} - \frac{c_{\alpha}}{m_{\alpha}\omega_{\alpha}^{2}}Q_{n}\right)\nonumber\\ \ddot{X}_{\alpha n} &= -\omega_{\alpha}^{2}\left(X_{\alpha n} - \frac{c_{\alpha}}{m_{\alpha}\omega_{\alpha}^{2}}Q_{n}\right) \end{align}\tag{67}\] with \(f^{(n)}_{M}(\mathbf{Q})\) defined in Eq. (36 ). Solving for the bath dynamics up to time \(t\), we obtain \[\begin{gather} \label{eq:app95mats95GLE95derivation} m\ddot{Q}_{n}(t) = f^{(n)}_{M}(\mathbf{Q}(t)) - \int_{0}^{t}\!\mathop{}\!\mathrm{d}s\, \zeta(t-s)\dot{Q}_{n}(s) \\ + \sum_{\alpha=1}^{n_{b}}c_{\alpha}\left[\left(X_{\alpha, n} - \frac{c_{\alpha}}{m_{\alpha}\omega_{\alpha}^{2}} Q_{n} \right)\cos{\omega_{\alpha}t} + \frac{P_{\alpha, n}}{m_{\alpha}\omega_{\alpha}}\sin{\omega_{\alpha}t}\right] \end{gather}\tag{68}\] with \(\zeta(t)\) given in Eq. (7 ), \(X_{\alpha,n}\equiv X_{\alpha,n}(t=0)\) and similarly for \(Q_n\) and \(P_{\alpha,n}\). Converting from \({\boldsymbol{P}}_\text{b}\) to \({\boldsymbol{\Pi}}_\text{b}\), we then obtain \[\begin{gather} m\ddot{Q}_{n}(t) = f^{(n)}_{M}(\mathbf{Q}(t)) - \int_{0}^{t}\!\mathop{}\!\mathrm{d}s\, \zeta(t-s)\dot{Q}_{n}(s) \\+ \sum_{\alpha=1}^{n_{b}}c_{\alpha}\bigg[ {\overline{X}}_{\alpha,n} \cos{\omega_{\alpha}t} + \left(\frac{\Pi_{\alpha, n}}{m_{\alpha}\omega_{\alpha}} + i \frac{\omega_{n}}{\omega_{\alpha}}{\overline{X}}_{\alpha, \bar{n}}\right)\sin{\omega_{\alpha}t} \\- \frac{c_{\alpha}\omega_{n}}{m_{\alpha}\omega_{\alpha}\omega_{\alpha n}^{2}}\left(\frac{\omega_{n}}{\omega_{\alpha}}Q_{n}\cos\omega_{\alpha}t - iQ_{\bar{n}}\sin\omega_{\alpha}t\right) \bigg]. \end{gather}\] Comparison of the variances of \({\boldsymbol{\Pi}}_\text{b}\) and \({\boldsymbol{X}}\) in Eq. (66 ) with Eq. (37 ) shows that the first three terms in the sum correspond to \(R_{n}^{\mathbb{C}}(t)\) in the GLE [Eq. (35 )]; comparison with Eq. (38 ) shows that the last two terms correspond to \(- \left[Q_{n}K_{n}(t) - iQ_{\bar{n}}L_{n}(t) \right]\).

References↩︎