Modelling a Bistable System Strongly Coupled to a Debye Bath: A   Quasiclassical Approach Based on the Generalised Langevin Equation

L. Stella; H. Ness; C.D. Lorenz; and L. Kantorovich

arXiv:1704.03918·cond-mat.stat-mech·July 23, 2018

Modelling a Bistable System Strongly Coupled to a Debye Bath: A Quasiclassical Approach Based on the Generalised Langevin Equation

L. Stella, H. Ness, C.D. Lorenz, and L. Kantorovich

PDF

TL;DR

This paper extends the Generalised Langevin Equation (GLE) to model quantum fluctuations in a strongly coupled bistable system, capturing key quantum effects and transition rate behaviors with a quasiclassical approach.

Contribution

The authors develop a quasiclassical GLE method that incorporates quantum fluctuations, enabling more accurate modeling of transition rates in strongly coupled bistable systems.

Findings

01

Isotopic effects observed in transition rates.

02

Quantum and classical rates converge at strong coupling.

03

Saturation of rates at low temperature qualitatively matches predictions.

Abstract

Bistable systems present two degenerate metastable configurations separated by an energy barrier. Thermal or quantum fluctuations can promote the transition between the configurations at a rate which depends on the dynamical properties of the local environment (i.e., a thermal bath). In the case of classical systems, strong system-bath interaction has been successfully modelled by the Generalised Langevin Equation (GLE) formalism. Here we show that the efficient GLE algorithm introduced in Phys. Rev. B 89, 134303 (2014) can be extended to include some crucial aspects of the quantum fluctuations. In particular, the expected isotopic effect is observed along with the convergence of the quantum and classical transition rates in the strong coupling limit. Saturation of the transition rates at low temperature is also retrieved, in qualitative, yet not quantitative, agreement with the…

Tables2

Table 1. Table 1: Numerical values of the parameters, A 𝐴 A , B 𝐵 B , and ω ‡ superscript 𝜔 ‡ \omega^{\ddagger} (errors on the last digit) appearing in Eq. ( 41 ) as a function of the dimensionless coupling strength, γ 𝛾 \gamma , along with the estimates of the quantum temperature, T q subscript 𝑇 𝑞 T_{q} , the critical temperature, T c subscript 𝑇 𝑐 T_{c} , and the Grote-Hynes-Pollak frequency, ω g h p subscript 𝜔 𝑔 ℎ 𝑝 \omega_{ghp} , in the case of the hydrogen mass.

$γ$	$A$ [fs^-1]	$B$ [eV]	$ω^{‡}$ [fs^-1]	$T_{q}$ [K]	$T_{c}$ [K]	$ω_{g h p}$ [fs^-1]
$0.1$	$0.00157$	$0.0359$	$0.0144$	$67.1$	$15.1$	$0.0124$
$0.2$	$0.00503$	$0.0371$	$0.0139$	$62.9$	$14.6$	$0.0120$
$0.5$	$0.00916$	$0.0403$	$0.0097$	$40.5$	$11.2$	$0.0092$
$1.0$	$0.00439$	$0.0426$	$0.0040$	$15.9$	$4.6$	$0.0038$

Table 2. Table 2: Numerical values of the parameters, A 𝐴 A , B 𝐵 B , and ω ‡ superscript 𝜔 ‡ \omega^{\ddagger} (errors on the last digit) appearing in Eq. ( 41 ) as a function of the dimensionless coupling strength, γ 𝛾 \gamma , along with the estimates of the quantum temperature, T q subscript 𝑇 𝑞 T_{q} , the critical temperature, T c subscript 𝑇 𝑐 T_{c} , and the Grote-Hynes-Pollak frequency, ω g h p subscript 𝜔 𝑔 ℎ 𝑝 \omega_{ghp} , in the case of the deuterium mass.

$γ$	$A$ [fs^-1]	$B$ [eV]	$ω^{‡}$ [fs^-1]	$T_{q}$ [K]	$T_{c}$ [K]	$ω_{g h p}$ [fs^-1]
$0.1$	$0.00157$	$0.0360$	$0.0105$	$48.8$	$10.6$	$0.0087$
$0.2$	$0.00502$	$0.0386$	$0.0101$	$43.8$	$10.0$	$0.0083$
$0.5$	$0.00699$	$0.0431$	$0.0064$	$25.0$	$6.6$	$0.0055$
$1.0$	$0.00236$	$0.0434$	$0.0029$	$11.1$	$2.3$	$0.0019$

Equations88

\overset{r}{˙}

\overset{r}{˙}

\overset{p}{˙}

\overset{s}{˙}_{1}^{(k)}

\overset{s}{˙}_{2}^{(k)}

⟨ ξ_{1}^{(k)} (t) ⟩ = ⟨ ξ_{2}^{(k)} (t) ⟩ = ⟨ ξ_{1}^{(k)} (t) ξ_{2}^{(k^{'})} (t^{'}) ⟩ = ⟨ ξ_{2}^{(k)} (t) ξ_{1}^{(k^{'})} (t^{'}) ⟩ = 0,

⟨ ξ_{1}^{(k)} (t) ⟩ = ⟨ ξ_{2}^{(k)} (t) ⟩ = ⟨ ξ_{1}^{(k)} (t) ξ_{2}^{(k^{'})} (t^{'}) ⟩ = ⟨ ξ_{2}^{(k)} (t) ξ_{1}^{(k^{'})} (t^{'}) ⟩ = 0,

⟨ ξ_{2}^{(k)} (t) ξ_{2}^{(k^{'})} (t^{'}) ⟩ = ⟨ ξ_{1}^{(k)} (t) ξ_{1}^{(k^{'})} (t^{'}) ⟩ = δ_{k k^{'}} δ (t - t^{'}) .

⟨ ξ_{2}^{(k)} (t) ξ_{2}^{(k^{'})} (t^{'}) ⟩ = ⟨ ξ_{1}^{(k)} (t) ξ_{1}^{(k^{'})} (t^{'}) ⟩ = δ_{k k^{'}} δ (t - t^{'}) .

\overset{r}{˙}

\overset{r}{˙}

\overset{p}{˙}

K (t - t^{'}; r, r^{'}) = μ k = 0 \sum K [G_{k} (r) G_{k} (r^{'}) e^{- \frac{1}{τ _{k}} (t - t^{'})} cos (ω_{k} (t - t^{'}))] θ (t - t^{'}),

K (t - t^{'}; r, r^{'}) = μ k = 0 \sum K [G_{k} (r) G_{k} (r^{'}) e^{- \frac{1}{τ _{k}} (t - t^{'})} cos (ω_{k} (t - t^{'}))] θ (t - t^{'}),

η (t; r) = \mbox R e ⎩ ⎨ ⎧ k = 0 \sum K \frac{2 μ h ( ω _{k} ) k _{B} T}{τ _{k}} \int_{- \infty}^{t} \mbox d t^{'} G_{k} (r) e^{- (\frac{1}{τ _{k}} + i ω_{k}) (t - t^{'})} ξ^{(k)} (t^{'}) ⎭ ⎬ ⎫,

η (t; r) = \mbox R e ⎩ ⎨ ⎧ k = 0 \sum K \frac{2 μ h ( ω _{k} ) k _{B} T}{τ _{k}} \int_{- \infty}^{t} \mbox d t^{'} G_{k} (r) e^{- (\frac{1}{τ _{k}} + i ω_{k}) (t - t^{'})} ξ^{(k)} (t^{'}) ⎭ ⎬ ⎫,

h (ω) = \frac{ℏ ω}{2 k _{B} T} coth (\frac{ℏ ω}{2 k _{B} T}) = {\frac{1}{2} \frac{ℏ ∣ ω ∣}{k _{B} T} 1 for T \to 0 \mbox or ω \to \infty, for T \to \infty \mbox or ω \to 0 .

h (ω) = \frac{ℏ ω}{2 k _{B} T} coth (\frac{ℏ ω}{2 k _{B} T}) = {\frac{1}{2} \frac{ℏ ∣ ω ∣}{k _{B} T} 1 for T \to 0 \mbox or ω \to \infty, for T \to \infty \mbox or ω \to 0 .

⟨ η (t; r (t)) η (t^{'}; r (t^{'})) ⟩ = k_{B} T K^{(q)} (t - t^{'}; r (t), r (t^{'})),

⟨ η (t; r (t)) η (t^{'}; r (t^{'})) ⟩ = k_{B} T K^{(q)} (t - t^{'}; r (t), r (t^{'})),

K^{(q)} (t - t^{'}; r, r^{'}) = μ k = 0 \sum K [h (ω_{k}) G_{k} (r) G_{k} (r^{'}) e^{- \frac{1}{τ _{k}} (t - t^{'})} cos (ω_{k} (t - t^{'}))] θ (t - t^{'}) .

K^{(q)} (t - t^{'}; r, r^{'}) = μ k = 0 \sum K [h (ω_{k}) G_{k} (r) G_{k} (r^{'}) e^{- \frac{1}{τ _{k}} (t - t^{'})} cos (ω_{k} (t - t^{'}))] θ (t - t^{'}) .

K^{(q)} (t - t^{'}; r, r^{'}) = μ k = 0 \sum \infty [h (ω_{k}) G_{k} (r) G_{k} (r^{'}) cos (ω_{k} (t - t^{'}))] θ (t - t^{'}) = \frac{2}{π} \int_{- \infty}^{\infty} \frac{\mbox d ω}{ω} h (ω) J (ω; r, r^{'}) cos (ω (t - t^{'})),

K^{(q)} (t - t^{'}; r, r^{'}) = μ k = 0 \sum \infty [h (ω_{k}) G_{k} (r) G_{k} (r^{'}) cos (ω_{k} (t - t^{'}))] θ (t - t^{'}) = \frac{2}{π} \int_{- \infty}^{\infty} \frac{\mbox d ω}{ω} h (ω) J (ω; r, r^{'}) cos (ω (t - t^{'})),

J (ω; r, r^{'}) = \frac{π μ ω}{2} k = 0 \sum \infty G_{k} (r) G_{k} (r^{'}) δ (ω - ω_{k}) θ (ω) .

J (ω; r, r^{'}) = \frac{π μ ω}{2} k = 0 \sum \infty G_{k} (r) G_{k} (r^{'}) δ (ω - ω_{k}) θ (ω) .

⟨ η (t; r (t)) η (t^{'}; r (t^{'})) ⟩ = \frac{ℏ}{π} \int_{- \infty}^{\infty} \mbox d ω coth (\frac{ℏ ω}{2 k _{B} T}) J (ω; r, r^{'}) cos (ω (t - t^{'})),

⟨ η (t; r (t)) η (t^{'}; r (t^{'})) ⟩ = \frac{ℏ}{π} \int_{- \infty}^{\infty} \mbox d ω coth (\frac{ℏ ω}{2 k _{B} T}) J (ω; r, r^{'}) cos (ω (t - t^{'})),

J (ω; r, r^{'}) = \frac{μ ω}{2} k = 0 \sum K G_{k} (r) G_{k} (r^{'}) [\frac{τ _{k}}{1 + ( ω - ω _{k} ) ^{2} τ _{k}^{2}} + \frac{τ _{k}}{1 + ( ω + ω _{k} ) ^{2} τ _{k}^{2}}] θ (ω)

J (ω; r, r^{'}) = \frac{μ ω}{2} k = 0 \sum K G_{k} (r) G_{k} (r^{'}) [\frac{τ _{k}}{1 + ( ω - ω _{k} ) ^{2} τ _{k}^{2}} + \frac{τ _{k}}{1 + ( ω + ω _{k} ) ^{2} τ _{k}^{2}}] θ (ω)

K (t - t^{'}; r, r^{'}) = \frac{2}{π} \int_{- \infty}^{\infty} \frac{\mbox d ω}{ω} J (ω; r, r^{'}) cos (ω (t - t^{'})) .

K (t - t^{'}; r, r^{'}) = \frac{2}{π} \int_{- \infty}^{\infty} \frac{\mbox d ω}{ω} J (ω; r, r^{'}) cos (ω (t - t^{'})) .

V (r) = V_{b} [1 - (\frac{r}{r _{min}})^{2}]^{2},

V (r) = V_{b} [1 - (\frac{r}{r _{min}})^{2}]^{2},

r_{min} = \frac{8 V _{b}}{m Ω _{0}^{2}},

r_{min} = \frac{8 V _{b}}{m Ω _{0}^{2}},

c_{k} = ⎩ ⎨ ⎧ \frac{1}{ω _{D}} \frac{3}{( 2 K + 1 )} \frac{1}{ω _{D}} \frac{6}{( 2 K + 1 )} \mbox i f k = 0, \mbox i f k > 0 .

c_{k} = ⎩ ⎨ ⎧ \frac{1}{ω _{D}} \frac{3}{( 2 K + 1 )} \frac{1}{ω _{D}} \frac{6}{( 2 K + 1 )} \mbox i f k = 0, \mbox i f k > 0 .

τ_{k} = τ = λ \frac{( 2 K + 1 )}{2 ω _{D}},

τ_{k} = τ = λ \frac{( 2 K + 1 )}{2 ω _{D}},

\frac{λ}{π} = (1 + 2 k = 1 \sum K \frac{1}{1 + k ^{2} λ ^{2} ( 1 + \frac{1}{2 K} ) ^{2}})^{- 1}

\frac{λ}{π} = (1 + 2 k = 1 \sum K \frac{1}{1 + k ^{2} λ ^{2} ( 1 + \frac{1}{2 K} ) ^{2}})^{- 1}

coth (z) = \frac{1}{z} + 2 z k = 1 \sum \infty \frac{1}{z ^{2} + π ^{2} k ^{2}},

coth (z) = \frac{1}{z} + 2 z k = 1 \sum \infty \frac{1}{z ^{2} + π ^{2} k ^{2}},

\frac{λ}{π} = (1 + 2 x^{2} k = 1 \sum \infty \frac{1}{x ^{2} + k ^{2} π ^{2}})^{- 1} = \frac{1}{x coth ( x )} .

\frac{λ}{π} = (1 + 2 x^{2} k = 1 \sum \infty \frac{1}{x ^{2} + k ^{2} π ^{2}})^{- 1} = \frac{1}{x coth ( x )} .

x = \frac{π}{λ} (1 + \frac{1}{2 K})^{- 1} .

x = \frac{π}{λ} (1 + \frac{1}{2 K})^{- 1} .

x = \mbox a r co t h (\frac{π}{λ x})

x = \mbox a r co t h (\frac{π}{λ x})

\frac{π}{λ} = (1 + \frac{1}{2 K}) \mbox a r co t h (1 + \frac{1}{2 K}) = \frac{1}{2} (1 + \frac{1}{2 K}) ln (1 + 4 K),

\frac{π}{λ} = (1 + \frac{1}{2 K}) \mbox a r co t h (1 + \frac{1}{2 K}) = \frac{1}{2} (1 + \frac{1}{2 K}) ln (1 + 4 K),

\mbox a r co t h (x) = \frac{1}{2} ln (\frac{x + 1}{x - 1}) .

\mbox a r co t h (x) = \frac{1}{2} ln (\frac{x + 1}{x - 1}) .

λ \sim \frac{2 π}{ln ( 4 K )},

λ \sim \frac{2 π}{ln ( 4 K )},

τ \sim \frac{2 π K}{ω _{D} ln ( 4 K )}

τ \sim \frac{2 π K}{ω _{D} ln ( 4 K )}

J (ω) = \frac{2Γ ω _{D} ω}{π ( 2 K + 1 )} [\frac{τ}{1 + ω ^{2} τ ^{2}} + k = 1 \sum K (\frac{τ}{1 + ( ω - ω _{k} ) ^{2} τ ^{2}} + \frac{τ}{1 + ( ω + ω _{k} ) ^{2} τ ^{2}})] θ (ω),

J (ω) = \frac{2Γ ω _{D} ω}{π ( 2 K + 1 )} [\frac{τ}{1 + ω ^{2} τ ^{2}} + k = 1 \sum K (\frac{τ}{1 + ( ω - ω _{k} ) ^{2} τ ^{2}} + \frac{τ}{1 + ( ω + ω _{k} ) ^{2} τ ^{2}})] θ (ω),

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

††thanks: This article is dedicated to William G. Hoover on his 80th birthday

Modelling a Bistable System Strongly Coupled to a Debye Bath: A Quasiclassical

Approach Based on the Generalised Langevin Equation

L. Stella

[email protected]

Atomistic Simulation Centre, School of Mathematics and Physics, Queen’s University Belfast, University Road, Belfast BT7 1NN, Northern Ireland, UK

H. Ness

Department of Physics, Faculty of Natural and Mathematical Sciences, King’s College London, Strand, London WC2R 2LS, UK

C.D. Lorenz

Department of Physics, Faculty of Natural and Mathematical Sciences, King’s College London, Strand, London WC2R 2LS, UK

L. Kantorovich

Department of Physics, Faculty of Natural and Mathematical Sciences, King’s College London, Strand, London WC2R 2LS, UK

Abstract

Bistable systems present two degenerate metastable configurations separated by an energy barrier. Thermal or quantum fluctuations can promote the transition between the configurations at a rate which depends on the dynamical properties of the local environment (i.e., a thermal bath). In the case of classical systems, strong system-bath interaction has been successfully modelled by the Generalised Langevin Equation (GLE) formalism. Here we show that the efficient GLE algorithm introduced in Phys. Rev. B 89, 134303 (2014) can be extended to include some crucial aspects of the quantum fluctuations. In particular, the expected isotopic effect is observed along with the convergence of the quantum and classical transition rates in the strong coupling limit. Saturation of the transition rates at low temperature is also retrieved, in qualitative, yet not quantitative, agreement with the analytic predictions. The discrepancies in the tunnelling regime are due to an incorrect sampling close to the barrier top. The domain of applicability of the quasiclassical GLE is also discussed.

Generalised Langevin equation; Quantum fluctuations; Debye bath; Quantum transition rate

I Introduction

The Generalised Langevin Equation (GLE)Zwanzig (1961); Mori (1965a, b); Zwanzig (1973) is a stochastic equation which describes a mechanical system subject to a random force or noise. At variance with the original Langevin equationLemons and Gythiel (1997) used to model Brownian motion, the GLE can deal with a wider range of random forces, e.g., with non-trivial time correlations. The GLE also provides a theoretical framework for the definition of the frequency dependent linear response of a mechanical system through the fluctuation-dissipation theorem.Kubo (1966); Zwanzig (2001)

The GLE has been employed to extend the classical transition rate theory of KramersKramers (1940); Hänggi et al. (1990) to include the coupling to realistic thermal baths with a frequency dependent spectral function.Pollak (1986a); Pollak et al. (1989); Tucker et al. (1991); Carmeli and Nitzan (1982, 1983); Marchesoni and Grigolini (1983) It is even possible to formulate a quantum GLEFord et al. (1965); Benguria and Kac (1981); Ford and Kac (1987); Ford et al. (1988a) to model the deviation from the classical transition rate theory at low temperature due to dissipative tunnelling (i.e., tunnelling without energy conservation).Pollak (1986b, c); Ford et al. (1988b) In fact, the problem of dissipative tunnelling has been solved theoretically by using path-integral techniques.Caldeira and Leggett (1981, 1983a, 1983b); Grabert et al. (1988) Alternative approaches which make use of a $c$ -number 111Following DiracDirac (1926), $c$ -numbers are “classical numbers” which always commute (e.g., complex numbers), while $q$ -numbers are “quantum numbers” which do not commute in general (e.g., linear operators). quantum GLESchmid (1982); Eckern et al. (1990); Stockburger and Grabert (2002); Barik et al. (2003a); Banerjee et al. (2004); Dammak et al. (2009); Kantorovich et al. (2016) are in principle better suited for numerical simulations since they are based on real-time equations of motion. However, the existing approaches are either more computationally demanding than the classical GLE or their applicability to the strong system-bath coupling regime has not been fully demonstrated, yet.

In this article we use a $c$ -number quantum GLE approach similar to the quasiclassical Langevin equation of SchmidSchmid (1982) or the quantum thermal bath of Dammak *et al.*Dammak et al. (2009). On the other hand, the GLE approach considered in this article is able to model a wider class of thermal bath and takes full advantage of the algorithmic development introduced in Ref. Stella et al., 2014. In this way, we are able to investigate the strong coupling regime of a bistable system coupled to a Debye bath, i.e., a thermal bath with a sharp frequency cut-off. Here we employ the adjective “quasiclassical” to distinguish this approximate scheme from the exact $c$ -number quantum GLE introduced in Ref. Kantorovich et al., 2016. In particular, the main topic of this article is the domain of applicability of the quasiclassical GLE in the case of low temperature and strong system-bath coupling, while a detailed discussion of the $c$ -number quantum GLE formalism can be found in Ref. Kantorovich et al., 2016.

As in the case of similar quasiclassical approximations,Eckern et al. (1990) the GLE approach used in this article fails to model tunnelling (with or without energy conservation) at low temperature, while dissipative tunnelling — especially in the weak coupling regime — can be tackled by real-time GLE approaches which include quantum corrections to the force field.Banerjee et al. (2002a, b); Barik et al. (2003a) However, the addition of these quantum corrections comes at a computational price. In this article we demonstrate that the results of the quasiclassical GLE approach — the computational cost of which is essentially equal to that of its classical counterpart — are in surprisingly good agreement with the analytic predictions for the quantum transition ratesPollak (1986c); Grabert et al. (1987) in the strong coupling limit. In particular, the isotope effect and the convergence of the quantum and classical transition rates in the strong coupling limit are correctly modelled.

The article is organised as follows: in Sec. II, the quasiclassical GLE is introduced along with the relevant terminology. In Sec. III, the model bistable potential is defined, the main properties of the Debye bath discussed, and the capabilities of the quasiclassical GLE to model the quantum probability densities demonstrated for a light test particle (hydrogen or deuterium). In Sec. IV the classical and quasiclassical transition rates are investigated as a function of the particle mass and system-bath coupling strength. Finally, in Sec. V and VI the results of the quasiclassical GLE are discussed in detail and the conclusions about the domain of applicability of the quasiclassical GLE approach are drawn.

II Quasiclassical GLE

In this Section, we complete the GLE formalism introduced in Ref. Stella et al., 2014 to include the quantum delocalisation at low temperature. For the sake of simplicity, we consider only the case of one particle in one spatial dimension. The generalisation to many particles in three spatial dimensions is straightforward. This extension is similar to other approaches to the quantum Langevin equation based on the quantum fluctuation-dissipation theorem (QFDT). Ford et al. (1965); Benguria and Kac (1981); Ford and Kac (1987); Ford et al. (1988a); Schmid (1982); Cortes et al. (1985); Banerjee et al. (2004); Dammak et al. (2009); Ceriotti et al. (2009)

The quasiclassical GLE is integrated by means of the following complex Langevin equations:

[TABLE]

where $r$ and $p$ are the physical degrees of freedom (DoFs), $s_{1}^{\left(k\right)}$ and $s_{2}^{\left(k\right)}$ are $K+1$ pairs of auxiliary DoFs, $m$ is the physical mass, $\mu$ is the mass of the auxiliary DoFs, $V\left(r\right)$ is the physical potential,222In fact, $V\left(r\right)$ includes a “polaronic” correction due to the system-bath interaction. For the sake of simplicity, we do not treat this correction explicitly in the rest of this article and we treat $V\left(r\right)$ as if it was independent of the system-bath coupling strength. $G_{k}\left(r\right)$ are the (dimensional) coupling strengths, $\tau_{k}$ are the relaxation times of the pair of auxiliary DoFs, $\omega_{k}\geq 0$ are the frequencies of the auxiliary DoFs, and ${\bf\xi}_{1}^{\left(k\right)}$ and ${\bf\xi}_{2}^{\left(k\right)}$ are pairs of uncorrelated sources of white Gaussian noise, i.e., stochastic processes with zero average,

[TABLE]

and the following 2-point correlation function:333Higher order correlations function are computed by means of the Wick’s theorem.

[TABLE]

Following the derivation used in Ref. Stella et al., 2014, the exact integration of the equations of motion (EoMs) of the complex auxiliary DoFs, $s^{\left(k\right)}=s_{1}^{\left(k\right)}+is_{2}^{\left(k\right)}$ , and its substitution into the second line of Eq. (1) yield the quasiclassical GLE

[TABLE]

with the (classical) memory kernel defined as

[TABLE]

and the (complex) coloured Gaussian noise

[TABLE]

where ${\bf\xi}^{(k)}={\bf\xi}_{1}^{(k)}+i{\bf\xi}_{2}^{(k)}$ . Note that the noise includes the quantum weight444The quantum weight is the ratio between the internal energy of an independent bosonic oscillator, $\frac{\hbar\omega}{2}\coth\left(\frac{\hbar\omega}{2k_{B}T}\right)$ , and the internal energy of an independent classical oscillator, $k_{B}T$ , of equal frequency, $\omega$ . In both cases, “independent” means “in the limit of vanishing small system-bath interaction”.

[TABLE]

The parameters $\tau_{k}$ and $\omega_{k}$ , along with the coupling strengths, $G_{k}\left(r\right)$ , can be either deduced from the first principle system-bath Lagrangian (in the classical case, Kantorovich and Rompotis (2008); Stella et al. (2014)) or fitted to an approximate memory kernel, $\mathcal{K}$ , obtained from benchmark molecular dynamics simulations.Ness et al. (2015, 2016); Gottwald et al. (2015, 2016) While the second case is most useful in practise, the exact mapping between the first principle Lagrangian and the parametrisation of the GLE kernel ensures that both the equilibrium and relaxation of the physical DoFs are correctly modelled, at least in the classical case.

The coloured Gaussian noise defined in Eq. (6) has zero average, $\left\langle{\eta}\left(t;r\left(t\right)\right)\right\rangle=0$ , while the 2-point correlation function is given by

[TABLE]

where the quantum memory kernel is defined as

[TABLE]

A crucial difference between Eq. (5) and Eq. (9) is the presence of the quantum weight in the second equation, although we have that $\mathcal{K}^{\left(q\right)}\to\mathcal{K}$ in the limit of either $T\to\infty$ or $\omega\to 0$ (i.e., in the classical limit, see Eq. (7)).

In order to faithfully reproduce the quantum delocalisation close to a minimum of the physical potential, $V\left(r\right)$ , the QFDT must hold. This is indeed the case in the limit of infinitely many auxiliary DoFs, $K\to\infty$ . In this limit, we have that $\tau_{k}\to\infty$ (see Sec. III) and we can rewrite the quantum kernel asCortes et al. (1985)

[TABLE]

where we have introduced the spectral density:

[TABLE]

By means of Eq. (11) and Eq. (7), the 2-point correlation function of the noise, Eq. (8), can be expressed as

[TABLE]

which is a most familiar form of the QFDT. Note that the noise correlation saturates in the limit of $T\to 0$ , while in the classical case (obtained by fixing $h\mbox{$ \left(\omega\right) $}=1$ ) we have that $\left\langle{\eta}\left(t;r\left(t\right)\right){\eta}\left(t^{\prime};r\left(t^{\prime}\right)\right)\right\rangle\to 0$ .

The QFDT is only approximately satisfied for a finite number of auxiliary DoFs. In this case, one can still define a spectral density

[TABLE]

so that

[TABLE]

However, Eq. (10) does not hold strictly because of the frequency dependence of $h\left(\omega_{k}\right)$ , which cannot be factorised out of the (finite) summation over the index, $k$ , of the auxiliary degrees of freedom. In practise, numerical convergence of the correlation functions and other figures of merit must be verified for each model of the environment. In the case of the Debye bath considered in Sec. III, convergence is quickly achieved (namely, for $K=50$ ) in the weak coupling regime, although extra care must be paid to the strong coupling regime.Stella et al. (2014)

III Model bistable system coupled to a Debye bath

We model the bistable system by means of the quartic double-well potential

[TABLE]

where $V_{b}$ is the barrier height and the two equivalent minima are located at $r=\pm r_{min}$ . To investigate the possible relevance of the isotope effect, the mass of the test particle is taken either as $m=m_{H}=1.0079$ amu (hydrogen) or $m=m_{D}=2.0141$ amu (deuterium). This model is artificial, but simple enough to provide neat results about the transition rates (see Sec. IV). On the other hand, it can also serve as a first step towards the application of the quasiclassical GLE to model hydrogen-bonded solids and liquids.

A natural unit of energy is provided by the Debye energy of the bath, $k_{B}T_{D}=\hbar\omega_{D}$ . The barrier height is then fixed to be $V_{b}=3k_{B}T_{D}$ and the potential minima are defined using

[TABLE]

after the harmonic frequency of the two equivalent minima has been fixed at $\Omega_{0}=0.8\left(\sqrt{m_{H}/m}\right)\omega_{D}$ . The presence of the square root of the mass ratio makes the harmonic constant (i.e., the second derivative of the potential at $r=\pm r_{min}$ ), $m\Omega_{0}^{2}$ , a geometric parameter independent of the particle mass, $m$ , as expected. The selected values of the barrier height and harmonic constant make possible to sample the probability densities (see Fig. 1 and 2) and the transition rates (see Fig. 4 and 4) by direct molecular dynamics simulations. In the case of the hydrogen mass, the choice of the harmonic frequency agrees with the example considered in Ref. Stella et al., 2014.

The Debye bath is defined by means of its Debye temperature, $T_{D}$ , the dimensionless system-bath coupling strength, $\gamma$ , and the auxiliary mass, $\mu$ . In particular, we consider the values $T_{D}=170$ K and $\gamma=0.1$ , $0.2$ , $0.5$ , $1.0$ . Having fixed, $\mu=m$ , the parameters $G_{k}$ ,555For the sake of simplicity, here we assume that the $G_{k}$ are trivial (i.e., constant) function of $r$ . $\tau_{k}$ , and $\omega_{k}$ in Eq. (1) depends on $T_{D}$ and $\gamma$ , only. Following Ref. Stella et al., 2014, we choose a uniform sampling of the frequency interval $\left[-\omega_{D},\omega_{D}\right]$ , i.e., $\omega_{k}=\frac{k}{K}\omega_{D}$ , with $k=0,1,2,\dots,K$ . We then write that $G_{k}=g_{0}c_{k}$ , where $g_{0}=\gamma\left(m/m_{H}\right)\Omega_{0}^{2}$ and

[TABLE]

Once again, the presence of the mass ratio makes the parameters $G_{k}$ independent of the particle mass, $m$ , and, in the case of the hydrogen mass, the choice of the parameters agrees with the example considered in Ref. Stella et al., 2014. For the sake of simplicity, we choose an equal decay time for all the auxiliary DoFs,

[TABLE]

with the auxiliary constant $\lambda$ defined through the self-consistent equation

[TABLE]

in order to retrieve the exact behaviour of the memory kernel in the two limits of $\omega\to 0$ and $\omega\to\infty$ .Stella et al. (2014) This choice of the parameters $\omega_{k}$ , $c_{k}$ and $\tau_{k}$ has been preferred to a least squares fit because it yields a more transparent convergence to the spectral density in the limit of $K\to\infty$ (see below).

By means of the Mittag-Leffler expansion of the hyperbolic cotangent,

[TABLE]

in the limit of $K\to\infty$ , we can express Eq. (19) as

[TABLE]

where

[TABLE]

Hence, the self-consistent equation can be written as

[TABLE]

which yields

[TABLE]

where we have used that

[TABLE]

Solving the last equation for $\lambda$ , we can also estimate the asymptotic behaviour in the limit of $K\to\infty$ ,

[TABLE]

which yields

[TABLE]

and the expected limit of $\tau\to\infty$ if $K\to\infty$ (see Sec. II).

By means of Eq. (13), we write the spectral density of the Debye bath as

[TABLE]

where the effective friction constant, $\Gamma$ , is defined by the equation

[TABLE]

As usual, the presence of the mass ratio makes the parameter $\Gamma/m$ independent of the particle mass, $m$ , and, in the case of the hydrogen mass, the definition of $\Gamma$ agrees with the example considered in Ref. Stella et al., 2014. We also note that the integral of the spectral density

[TABLE]

does not depend on the number of pairs of DoFs, $K+1$ , and that the spectral density has an algebraic asymptotic behaviour

[TABLE]

in the limit of $\omega\gg\omega_{D}$ . It can be also proven that, in the limit of $K\to\infty$ , the spectral density in Eq. (28) converges to the expected

[TABLE]

where $\chi_{[0,\omega_{D}]}$ is the characteristic function of the interval $[0,\omega_{D}]$ .

Despite the apparent simplicity of the Debye model, the limit $K\to\infty$ is not entirely trivial Kemeny et al. (1986). As shown in Ref. Stella et al., 2014, a persistent (i.e., undamped) oscillation with a frequency larger than the Debye frequency, $\omega_{D}$ , is observed in the strong coupling regime. A thorough discussion of this persistent oscillation is neither brief nor pertinent to the main topic of this article and it is then left to a future publication.

In Fig. 1 we show the probability densities obtained by numerical integration of the quasiclassical complex Langevin equations introduced in Eq. (1) with $K=50$ for the case of the hydrogen mass. The numerical integration provides an accurate solution of the equivalent quasiclassical GLE defined in Eq. (4). Details of the integration algorithm can be found in Ref. Stella et al., 2014. 666In fact, in this work we have preferred the “BAOAB” algorithm of Leimkuhler and Matthews which is known to give better estimates of the confrontational averages.Leimkuhler and Matthews (2013) This choice implies a simple rearrangements of the split algorithm employed in Ref. Stella et al., 2014. For each value of the temperature, $T$ , and the dimensionless coupling strength, $\gamma$ , $50$ independent trajectories have been generated to sample the position histograms. Each trajectory is randomly started at rest in either the left or the right minima with equal probability. A time step of $1$ fs and $10^{8}$ steps have been used, while the configurations in the extended phase space $\left(r,p,s_{1}^{\left(k\right)},s_{2}^{\left(k\right)}\right)$ have been recorded every $10^{4}$ steps. In each panel, we have also indicated the classical probability density, $P_{cl}\left(x\right)\propto\exp\left(-V\left(x\right)/k_{B}T\right)$ , and the quantum probability density, $P_{q}\left(x\right)\propto\sum_{n}\left|\phi_{n}\left(x\right)\right|^{2}\exp\left(-E_{n}/k_{B}T\right)$ , where $\phi_{n}$ and $E_{n}$ are the eigenvectors and eigenvalues of the Hamiltonian operator $H=-\left(\hbar^{2}/2m\right)\nabla^{2}+V\left(x\right)$ .

From the results shown in the different panels of Fig. 1, we can conclude that the quasiclassical GLE is rather accurate in modelling the quantum probability density when the temperature is not too low and the coupling is not too strong. This conclusion agrees with previous observations.Ceriotti et al. (2009) The capability of a quantum GLE scheme based on the QFDT to model the quantum delocalisation in a moderately anharmonic potential has been exploited to improve the convergence of path-integral molecular dynamics.Ceriotti et al. (2009, 2011) Discrepancies at low temperature are due to the lack of ergodicity which follows a reduced transition rate, $\kappa_{gle}$ (see Sec. IV). Discrepancies in the strong coupling regime are due to the non-negligible corrections to the quantum probability density caused by the system-bath interaction.Riseborough et al. (1985) A detailed assessment of these corrections depends on the characterisation of the persistent oscillation of a Debye bath (see above) and it is therefore left to a future publication. In fact, the main conclusions about the transition rate in the strong coupling regimes (see Sec. IV) do not depend on the detailed assessment of these corrections.

To investigate the isotope effect, in Fig. 2 we show the probability densities obtained by numerical integration of the quasiclassical complex Langevin equations for a deuterium atom. The same trends with decreasing $T$ and increasing $\gamma$ are observed, even if the quantum probability densities are more localised and the ergodicity breaking is then more severe.

IV Transition rates

The transition rate, $\kappa$ , has been estimated from the decay of the position autocorrelation functionGillan (1987)

[TABLE]

sampled at $T=13$ , $16$ , $20$ , $25$ , $31$ , $40$ , $50$ , $63$ , $79$ , $100$ , $126$ , $159$ , $200$ K and $\gamma=0.1$ , $0.2$ , $0.5$ , $1.0$ . The remaining simulation parameters are the same as for the trajectories used to investigate the probability density (See Sec. III).

In Fig. 3 we show the Arrhenius plots for the different values of the dimensionless coupling strength in the case of the hydrogen mass (left panels) or in the case of the deuterium mass (right panels). The transition rates obtained from the numerical solution of the quasiclassical GLE saturate at very low temperature, while the transition rates from the classical GLE (obtained by fixing $h\left(\omega\right)=1$ in Eq. (1)) display the familiar linear behaviour.

In all the panels of Fig. 3 we also show the analytic estimates of the classical transition rate (Grote-Hynes-PollakGrote and Hynes (1980); Pollak (1986a))

[TABLE]

where $-i\omega_{b}$ is the imaginary frequency of the barrier top (i.e., $-m\omega_{b}^{2}$ is the second derivative of the potential at the barrier top) and

[TABLE]

is the bare estimate of the transition state theory.Eyring (1935) In the case of a Debye bath, the Grote-Hynes-Pollak imaginary frequency, $i\omega_{ghp}$ , is given by the positive solution of the equation

[TABLE]

The numerical values of $\omega_{ghp}$ are reported in Table 1, for the case of the hydrogen mass, or in Table 2 in the case of the deuterium mass.

The quantum transition rate (WolynesWolynes (1981); Pollak (1986b)) is approximated asBell (1959)

[TABLE]

The agreement between the Grote-Hynes-Pollak predictions and the classical GLE is very good in all cases, except the very weak, i.e., $\gamma=0.1$ , coupling regime. A disagreement is expected in this case since the transitions are limited more by energy diffusion than spatial diffusion (the Kramers turnover problem).Hänggi et al. (1990) Corrections to Eq. (34) are known,Mel’nikov and Meshkov (1986); Grabert (1988); Pollak et al. (1989) but are not relevant in the strong coupling regime.

The quasiclassical GLE rates are in clear quantitative disagreement with the Wolynes predictions in all cases, except the very strong, i.e., $\gamma=1.0$ , coupling regime. In particular, the quasiclassical GLE rates saturate at a temperature well above the so-called critical temperature,

[TABLE]

at which Eq. (37) displays an (apparent) divergence.Grabert et al. (1987); Hänggi and Hontscha (1988, 1991)

Following Miller,Miller (1977) one can use the critical temperature, $T_{c}$ , to characterise the tunnelling through the barrier top, approximated as $V\left(r\right)\approx V_{b}-\frac{1}{2}m\omega_{b}^{2}r^{2}$ . In this case, the tunnelling probability is

[TABLE]

where $E$ is the total energy of the system. By means of Eq. (39), the quantum transition rate is estimated as

[TABLE]

Note that Eq. (37) and Eq. (40) differ only in the choice of imaginary frequency of the barrier top. In fact, we have that $\omega_{ghp}\to\omega_{b}$ in the limit of $\gamma\to 0$ (see Eq. (36) and Eq. (29)).

At low temperature, $k_{B}T\ll V_{b}$ , assuming that the zero-point energy is negligible, $\hbar\Omega_{0}/2\ll V_{b}$ , one can substitute $E\approx 0$ into Eq. (39) to approximate $P_{t}\left(0\right)\approx e^{-\frac{V_{b}}{k_{B}T_{c}}}$ . As a consequence, tunnelling is expected to be the dominant transition mechanism for $T<T_{c}$ . The numerical values of $T_{c}$ are reported in Table 1 or Table 2. The critical temperature, $T_{c}$ , is a function of both the dimensionless coupling strength, $\gamma$ , and the mass, $m$ , through its dependence on $\omega_{ghp}$ . The regions corresponding to $T<T_{c}$ have been shaded in the panels of Fig. 3.

To help interpret the quasiclassical GLE results, we model the quasiclassical transition rate by means of the function

[TABLE]

where $A$ , $B$ , and $\omega^{\ddagger}$ are adjustable parameters, the values of which are reported in Table 1,

or Table 2. 777One can set $B=V_{b}$ to reduce the number of adjustable parameters, but the fits generally get worse. The difference between $B$ and $V_{b}$ is due to the quantum fluctuations.Pollak (1986c)

The global accuracy of these fits can be better appreciated from Fig. 4, where we have reported only the quasiclassical GLE rates. By interpreting the exponential in Eq. (41) as a Boltzmann factor, we can define an effective quantum temperature, $T_{q}$ , so that

[TABLE]

The numerical values of $T_{q}$ are reported in Table 1, for the case of the hydrogen mass, or Table 2 in the case of the deuterium mass. The quantum temperature can be used to assess the validity of discrepancy between the quasiclassical GLE and the analytic predictions (see Sec. V).

V Discussion

The quasiclassical GLE formalism considered in this article does not include tunnelling.Eckern et al. (1990) It is then not surprising that it fails to model the transition rates in the deep quantum regime. On the other hand, it is not immediately clear why such a formalism overestimates instead of underestimating — as naively suggested by the absence of tunnelling — the quantum transition rates. In this Section we attempt an answer by discussing in more detail the results of Sec. III and IV.

First of all, from the results shown in Figs. 1 and 2, we know that the quasiclassical GLE reproduces rather accurately the probability density, at least close to the potential minima. In the limit of $T\to 0$ and $\gamma\to 0$ , the quantum fluctuations close to the minima are entirely due to the zero-point motion which can be characterised by the effective temperature $T_{zp}=\hbar\Omega_{0}/2k_{B}$ . In particular, we have that $T_{zp}=136$ K in the case of the hydrogen mass or $T_{zp}=96.2$ K in the case of the deuterium mass. Those temperatures are much larger than $T_{q}$ and $T_{c}$ in both cases. Given the good agreement between the classical GLE rates and the Grote-Hynes-Pollak formula for moderate to strong system-bath coupling, we can also exclude a large contribution from the finite height of the barrier (the condition $k_{B}T_{zp}\ll V_{b}$ is satisfied). It is then plausible that the discrepancies between the quasiclassical GLE and the analytic predictions originate from an incorrect sampling of the region close to the barrier top.

The quasiclassical GLE considered in this article provides an inherently thermal (i.e., classical) description of the random forces, although the physical temperature is weighted by a correcting factor, $h\left(\omega\right)$ , which depends on the frequency of the oscillations, $\omega$ , to mimic the quantum fluctuations. In practise, the quantum temperature, $T_{q}$ , can be used to estimate the effective temperature close to the barrier top.

Interestingly, we observe that the ratio between $T_{q}$ and $T_{c}$ is a relatively constant function of $\gamma$ and $m$ (between $0.4$ and $0.5$ , see Fig. 5). Our results are in agreement with the findings of Eckern et al. (see Ref. Eckern et al., 1990, in particular at the end of Sec. 2.3). This observation suggests that, despite the quantitative discrepancy between the quasiclassical GLE and the analytic predictions, the functional dependence of $\kappa_{gle}\left(T\right)$ on both $\gamma$ and $m$ is qualitatively correct. In particular, the isotope effect and the convergence of the quantum and classical transition rates in the strong coupling limit are in good qualitative agreement with the analytic predictions.Pollak (1986c); Grabert et al. (1987)

VI Conclusions

In this article, we have completed the GLE formalism introduced in Ref. Stella et al., 2014 to include the quantum delocalisation at very low temperature. Our results confirm the applicability of this formalism to model the equilibrium properties (e.g., the probability density) of a bistable system coupled to a Debye bath. In particular, the quasiclassical GLE formalism equally applies to both the weak and strong coupling regimes.

The quantitative discrepancy between the quasiclassical GLE and the analytic predictions for the quantum transition rates has been rationalised as the consequence of an incorrect sampling close to the barrier top. In particular, the quasiclassical GLE predicts a saturation of the transition rate at an effective quantum temperature, $T_{q}$ , which is roughly twice as large as the expected critical temperature, $T_{c}$ . Since the value of $T_{c}$ depends on the imaginary frequency of the barrier top (see Eq. (38)), we can conclude that the quasiclassical GLE effectively samples a different imaginary frequency. On the other hand, the quasiclassical GLE accurately samples the quantum probability distribution (at least for weak system-bath interaction) close to the minima of the potential wells. Since the ratio between $T_{q}$ and $T_{c}$ is roughly constant, the functional dependence of $\kappa_{gle}\left(T\right)$ on both the system-bath coupling strength, $\gamma$ , and the particle mass, $m$ , is also qualitatively correct. This qualitative agreement includes the isotope effect and the convergence of the quantum and classical transition rates in the strong coupling limit.

Our results shed more light on the domain of applicability of a real-time GLE approach to model the relaxation of quantum dissipative system. The simple quasiclassical approach considered in this article ignores both the quantum corrections to the force field Banerjee et al. (2002a); Barik et al. (2003a, b); Banerjee et al. (2004) and the proper treatment of the quantum fluctuations by means of the path integral formalismCeriotti et al. (2009, 2011); Richardson and Althorpe (2009); Hele et al. (2015). However, it is surprisingly accurate in the limit of strong system-bath interaction, at a computational cost essentially equal to that of its classical counterpart.

Acknowledgements

The authors acknowledge the financial support from EPSRC, grant EP/J019259/1. LS is also grateful to Professors Roberto D’Agosta and Ian Ford for many useful conversations.

Bibliography69

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Zwanzig (1961) R. Zwanzig, Phys. Rev. 124 , 983 (1961).
2Mori (1965 a) H. Mori, Prog. Theor. Phys. 33 , 423 (1965 a).
3Mori (1965 b) H. Mori, Prog. Theor. Phys. 34 , 399 (1965 b).
4Zwanzig (1973) R. Zwanzig, J. Stat. Phys. 9 , 215 (1973).
5Lemons and Gythiel (1997) D. S. Lemons and A. Gythiel, Am. J. Phys. 65 , 1079 (1997).
6Kubo (1966) R. Kubo, Rep. Prog. Phys. 29 , 255 (1966).
7Zwanzig (2001) R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, 2001).
8Kramers (1940) H. Kramers, Physica 7 , 284 (1940).