# Local Harmonic Approximation to Quantum Mean Force Gibbs State

Prem Kumar<EMAIL_ADDRESS><EMAIL_ADDRESS>Optics and Quantum
Information Group, The Institute of Mathematical Sciences, C.I.T. Campus,
Taramani, Chennai 600113, India.  Homi Bhabha National Institute, Training
School Complex, Anushakti Nagar, Mumbai 400094, India.

###### Abstract

When the strength of interaction between a quantum system and bath is non-
negligible, the equilibrium state can deviate from the Gibbs state. Here, we
obtain an approximate expression for such a mean force Gibbs state for a
particle in an arbitrary one dimensional potential, interacting with a bosonic
bath. This approximate state is accurate when either the system-bath coupling
or the temperature is large, or when the third and higher derivatives of the
potential are small compared to certain system-bath specific parameters. We
show that our result recovers the ultra strong coupling and high temperature
results recently derived in literature. We then apply this method to study
some systems like a quartic oscillator and a particle in a quartic double-well
potential. We also use our method to analyze the proton tunneling problem in a
DNA recently studied in literature [Slocombe et al., Comm. Phys., vol. 5, no.
1, p. 109, 2022], where our results suggest the equilibrium value of the
probability of mutation to be orders of magnitude lower than the steady state
value obtained there ($10^{-8}$ vs $10^{-4}$).

##  I Introduction

Continuous variable (CV) quantum systems have been immensely studied in the
context of fields like quantum information processing, communication and
computation braunstein2005quantum ; Ralph2006Quantum ; Djordjevic2022Hybrid ;
Inoue2022Toward ; Yang2015Electromagnetically ; Gross2011Atomic ;
Huang2019Generation , quantum resource theories and thermodynamics
Takagi2018Convex ; Ferrari2020Asymptotic ; Narasimhachar2019Thermodynamic ,
quantum technology Hoskinson2008Quantum ; Ferreira1997Tunnelling , quantum
chaos Monteoliva2000Decoherence ; B.C.Bag1998A ; Jt1992Semiclassical ;
Chakrabarty2019Out  and open quantum system Liu2007Non-markovian ;
Vasile2011Quantifying ; Hoerhammer2007Environment-induced ;
ruostekoski1998bose ; fisher1985dissipative ; grabert1985quantum ;
topaler1994quantum ; Oliveira2006Quantum-classical . Thermodynamic equilibrium
properties of these CV systems are hence of great interest in many of these
studies and applications. But when the strength of interaction between a
quantum system and bath is non-negligible, the equilibrium state can deviate
from the Gibbs state and is given by the mean force Gibbs state
(MFGS)cresser2021weak , defined as,

$\displaystyle\rho=Tr_{R}\left[\frac{e^{-\beta H_{tot}}}{Z}\right]$ (1)

where $H_{tot}$ is the full system-bath Hamiltonian and the partial trace is
taken over the bath degrees of freedom. The MFGS has found several
applications in fields like strong coupling quantum thermodynamics
seifert2016first ; philbin2016thermal ; jarzynski2004nonequilibrium ;
jarzynski2004nonequilibrium ; campisi2009fluctuation  and quantum chemistry
allen2006molecular ; maksimiak2003molecular ; roux1999implicit ;
roux1995calculation  and is relevant to many other problems in quantum
chemistry and biology.

For example, consider a particle in a double well (DW) potential
ruostekoski1998bose ; gillan1987quantum ; liu2021understanding ;
fisher1985dissipative ; grabert1985quantum ; topaler1994quantum  interacting
with a bosonic bath which finds application, for example, in quantum chemistry
and biology to study processes like chemical reactions occurring in the
presence of a solvent slocombe2022open ; slocombe2022quantum ; zhang2020proton
; bothma2010role . Finding an analytical expression MFGS for such a system is
hence of interest. But most of these physical processes happen in the
intermediate coupling regime, where an analytical expression for MFGS is not
available. By intermediate coupling regime, we mean the regime in which the
bath induced timescales on the system are of the same order as the timescales
associated with the bath itself (e.g., $\hbar/k_{B}T$) BRE02 .

Analytical expressions for MFGS in weak and ultra strong coupling (USC) limits
were derived recently cresser2021weak . Further, corrections to the USC regime
MFGS have been evaluated using dynamics trushechkin2022quantum  and
perturbative expansion methods latune2021steady . Additionally, high
temperature expansion for the MFGS was also derived recently
timofeev2022hamiltonian ; gelzinis2020analytical .

If the CV system is a harmonic oscillator (HO), then the expression for the
MFGS is known analytically grabert1988quantum ; weiss2012quantum . But for a
particle in an arbitrary 1D potential $V(q)$, analytical expression for MFGS
in intermediate coupling regime remains a challenge. Recently, two numerically
exact approaches, “Time Evolving matrix Product Operator” (TEMPO)
chiu2022numerical ; strathearn2018efficient ; strathearn2020modelling ;
feynman2000theory ; makri1995tensor1 ; makri1995tensor2  and ”The hierarchical
equations of motion” (HEOM) tanimura2020numerically ; liu2021understanding ;
zhang2020proton , have made these MFGS and steady state calculations much more
computationally efficient. But for a particle in an arbitrary 1D potential
$V(q)$, the system is infinite dimensional and a typical numerical technique
proceeds by truncating the system in energy eigenbasis topaler1993system . In
the intermediate coupling regime, this is an issue because the effective
dimension of the system increases with the coupling, causing the calculation
to become computationally very expensive.

Here, we provide an approximate analytical expression for the MFGS, which we
call local harmonic approximation (LHA) to MFGS, for a particle in a 1D
potential $V(q)$ interacting with a non interacting bosonic bath. This method
is accurate when either the system-bath coupling or the temperature is large,
or when the third and higher derivatives of the potential are small compared
to certain system-bath specific parameters. We find that LHA reduces to the
exact USC and high temperature results in the appropriate limits
cresser2021weak ; timofeev2022hamiltonian .

In Sec. II, we provide our main results. In Sec. III, we use LHA to study some
CV open quantum systems, like a quartic oscillator and a quartic double well
potential. Sec. III.3 provides an application of LHA to the proton tunneling
problem in DNA that has been recently investigated by Slocombe et al.
slocombe2022open ; slocombe2022quantum . Sec. IV provides conclusion and
discussions.

##  II Main results

We consider a particle of mass $m$ in an arbitrary potential $V(q)$, linearly
coupled to a bath of harmonic oscillators. This system is described by the
Hamiltonian,

$H=\frac{p^{2}}{2m}+V(q)\\\
+\sum_{k}\left[\frac{p_{k}^{2}}{2m_{k}}+\frac{1}{2}m_{k}\omega_{k}^{2}\left(q_{k}-\frac{c_{k}q}{m_{k}\omega_{k}^{2}}\right)^{2}\right],$
(2)

where $\omega_{k}$, $m_{k}$, and $c_{k}$ are the frequency, mass, and coupling
strength of the $k$th oscillator, respectively. For simplicity, we assume a
one dimensional system, though our results are easily generalized to higher
dimensions. Here and throughout, we assume $\hbar=1$. Throughout this article
we write operator matrix elements in position basis as
$O(q,\eta)=\bra{q+\eta}\hat{O}\ket{q-\eta}$. Diagonal elements are obtained
when $\eta=0$ and will be denoted by $O(q)=O(q,0)$.

Upon integrating out the bath degrees of freedom, the MFGS for the particle
can be expressed as a Euclidean path integral weiss2012quantum ,

$\rho(q,\eta)=Z_{R}^{-1}\int_{q+\eta}^{q-\eta}\mathcal{D}q(\tau)e^{-S_{0}[q(\tau)]-\Phi[q(\tau)]}$
(3)

where, $\int_{x}^{y}\mathcal{D}q(\tau)$ signifies $q(0)=x$ and $q(\beta)=y$ in
the functional integral. The reduced partition function is the ratio of bath
and system+bath partition functions, $Z_{R}=Z/Z_{B}$, and

$\displaystyle S_{0}[q(\tau)]=$
$\displaystyle\int_{0}^{\beta}d\tau\left(\frac{1}{2}m\dot{q}(\tau)^{2}+V\left(q(\tau)\right)\right)$
(4) $\displaystyle\Phi[q(t)]=$
$\displaystyle\frac{1}{4}\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}K\left(\tau-\tau^{\prime}\right)\left(q(\tau)-q(\tau^{\prime})\right)^{2}$
(5)

are the action for the free particle and the influence action, respectively.
The bath correlation function is defined as,

$\displaystyle K(\tau)$ $\displaystyle=\int_{0}^{\infty}d\omega
J(\omega)\frac{\cosh{[\omega(\beta/2-\tau)]}}{\sinh{[\omega\beta/2]}}$ (6)

in terms of the spectral density

$J(\omega)=\alpha\sum_{k}\frac{c_{k}^{2}}{2m_{k}\omega_{k}}\delta(\omega-\omega_{k}),$
(7)

where we have introduced the dimensionless number $\alpha$ to characterize the
strength of the coupling.

At large system-bath coupling the influence action [Eq. 5] dominates over the
free action, such that amplitudes for paths with $q(\tau)\neq
q(\tau^{\prime})$ are exponentially small in the path integral Eq. 3. In the
limit $\alpha\to\infty$, the ultra-strong coupling (USC) limit, the amplitudes
for these path vanish and only the constant path contributes,
$q(\tau)=q(\tau^{\prime})=q$ with $\dot{q}(\tau)=0$. Similarly, at high
temperature (but arbitrary coupling) the path length $\beta=1/k_{B}T$ becomes
small, so the kinetic term in the free action, $\propto\dot{q}(\tau)^{2}$,
increasingly suppresses amplitudes for paths that deviate significantly from
$q(\tau)=q$. In the limit $T\to\infty$ (zero path length) we once again find
that the only relevant path is the constant path with $\dot{q}(\tau)=0$. Thus,
in either the USC limit or infinite temperature limit we have,

$\displaystyle\rho(q,\eta)=\delta(\eta)Z^{-1}e^{-\beta V(q)},$ (8)

where now $Z=\int dqe^{-\beta V(q)}$. This is the same USC result obtained in
cresser2021weak  and is CV equivalent of the high temperature result given in
timofeev2022hamiltonian . Note that the argument used above to obtain the USC-
MFGS is equally applicable for discrete systems. See Appendix. D for more
details.

Given the discussion above, when either the coupling or temperature is large,
but not infinite, we expect that only paths with small deviation around the
constant path will contribute significantly to the path integral. In other
words, reparameterising the path variable as $q(\tau)=q+\delta(\tau)$, we
expect that only small values of $\delta(\tau)$ will be relevant to the final
result. Therefore, we propose a local harmonic approximation (LHA) for
$\rho(q,\eta)$, where for large coupling or temperature, the effect of the
potential locally around $q$ is effectively harmonic. Then, to second order in
$\delta(\tau)$, we can write

$V(q+\delta(\tau))\\\ \approx
V(q)+\delta(\tau)V_{1}(q)+\frac{1}{2}\delta(\tau)^{2}V_{2}(q)\\\
=V(q)-\frac{1}{2}m\omega_{q}^{2}Q_{q}^{2}+\frac{1}{2}m\omega_{q}^{2}\left(\delta(\tau)+Q_{q}\right)^{2}$
(9)

where

$\displaystyle\omega_{q}$ $\displaystyle=\sqrt{V_{2}(q)/m}$ (10)
$\displaystyle Q_{q}$ $\displaystyle=\frac{V_{1}(q)}{V_{2}(q)}$ (11)

are the frequency and the displacement of the effective harmonic potential,
respectively, and $V_{n}(q)$ is the $n$th derivative of $V(q)$. In Fig. 1 we
show the effective harmonic potentials found for different values of $q$ using
a double well potential for $V(q)$ as an example.

With Eq. 9, the path integral in Eq. 3 becomes Gaussian and can be computed
exactly to obtain

$\rho(q,\eta)=Z_{R}^{-1}\frac{\exp\left(-\beta
V(q)\right)}{\exp\left(-\beta\frac{1}{2}m\omega_{q}^{2}Q_{q}^{2}\right)}\tilde{Z}_{q}\tilde{\rho}_{q}(Q_{q},\eta)$
(12)

where $\tilde{\rho}_{q}(Q_{q},\eta)$ is the MFGS for a HO with frequency
$\omega_{q}$, coupled to the same bath as the original system, and
$\tilde{Z}_{q}$ is the corresponding partition function for the HO-bath system
which depends on $V(q)$ only through the effective frequency, $\omega_{q}$. We
highlight that both $\tilde{\rho}_{q}$ and $\tilde{Z}_{q}$ are straightforward
to compute since the HO-bath system is Gaussian and thus exactly solvable
grabert1988quantum ; weiss2012quantum . Although we have expressed
$\rho(q,\eta)$ in terms of the Gaussian state of a HO, we stress that Eq. 12
is generally not Gaussian with respect to $q$. This is because the potential,
$V(q)$, remains completely arbitrary and all quantities on the right hand side
of Eq. 12 are explicit or implicit functions of this potential.

While we have derived the LHA [Eq. 12] under the assumption of either strong
coupling or large temperature, we note that it is exact for arbitrary coupling
and temperature when the potential is actually harmonic, i.e. when
$V(q)=aq+bq^{2}$ for some $a$ and $b$, no matter how large $\delta(\tau)$
becomes. We therefore expect Eq. 12 to also be accurate even at weak coupling
and low temperature, where $\delta(\tau)$ may not be small, as long as the
third and higher derivatives of the potential are relatively small. This gives
our result great flexibility in its range of validity.

Figure 1: Fitting a harmonic potential at each point $q$ of a general
potential $V(q)$.

One possible problem with the LHA is that when $V_{2}(q)<0$, e.g., near the
maxima of the potential $V(q)$, the effective frequency of the HO defined in
Eq. 9 becomes imaginary, so the HO is inverted barton1986quantum ;
subramanyan2021physics ; baskoutas1993quantum  (see the orange curve in Fig.
1). For an inverted HO, the potential is unbounded from below as
$q\to\pm\infty$, so the MFGS for the HO, $\tilde{\rho}_{q}$ in Eq. 12, is
unphysical, unnormalisable, and gives divergent observables, e.g.
$\braket{\hat{q}^{2}}=\infty$. However, we show in Appendix. 2 that as long as
$V_{2}(q)>-V_{P}$, for some $V_{P}>0$, the combination
$\tilde{Z}_{q}\tilde{\rho}_{q}$ can still be finite, where the quantity
$\tilde{Z}_{q}\tilde{\rho}_{q}$ amounts to an unnormalised path integral which
has its important contribution still coming from a finite deviation from the
constant path $q(\tau)=q$, as required for LHA to work for general potentials.

This number, $V_{P}$, is given by the smallest solution of the following
equation in $x$,

$\sum_{n=1}^{\infty}\frac{x}{y_{n}-x}=\frac{1}{2},$ (13)

where

$y_{n}=m\nu_{n}^{2}+\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega}\frac{2\nu_{n}^{2}}{\omega^{2}+\nu_{n}^{2}},$
(14)

and $\nu_{n}=2\pi nk_{B}T$ are the Matsubara frequencies. Thus, provided that
$V_{2}(q)>-V_{P}$, the LHA remains finite even for $V_{2}(q)<0$ and then can
be applied to nontrivial potentials, like the one depicted in Fig. 1. We note
that in the large coupling or temperature limits we have assumed in deriving
the LHA, for $x<y_{1}$, the $n=1$ term in Eq. 13 dominates and we obtain the
approximate solution $V_{P}\approx y_{1}/3$. Since $y_{n}$ diverges with both
coupling and temperature, we find the condition $V_{2}(q)>-V_{P}$ is always
satisfied in this limit, consistent with LHA becoming exact there.

In order to estimate the error associated with LHA, let us first define two
length scales. For the effective quadratic potential approximated at point
$q$, let $q_{cl}(\tau)$ denote the classical path starting and ending at the
point $q+\eta$ and $q-\eta$. Then, let $\Delta(q,\eta)$ denote the maximum
deviation of the classical path from the constant path $q^{\prime}(\tau)=q$
and let $\mathcal{E}(q)$ denote the typical quantum deviation of all the paths
about $q_{cl}(\tau)$ that have significant contribution to the action. An
estimate for $\Delta(q,\eta)$ and $\mathcal{E}(q)$ is derived in Appendix. 1.

Then, the leading order correction to LHA is given by
$\mathcal{T}(q,\eta)\rho(q,\eta)$, where we have,

$\displaystyle\mathcal{T}(q,\eta)$
$\displaystyle\equiv-\beta\frac{V_{3}(q)}{3!}\left(\Delta(q,\eta)^{3}+3\mathcal{E}(q)^{2}\Delta(q,\eta)\right)$
(15)

The relative error associated with the expectation value of an observable
$\hat{O}$ is given as [Appendix. C],

$\displaystyle\epsilon_{\hat{O}}=\frac{1}{\braket{\hat{O}}}\int dq\int
d\eta\;\hat{O}(q,\eta)\epsilon(q,\eta)|\rho(q,\eta)|$ (16)

where, we have,

$\displaystyle\epsilon(q,\eta)$
$\displaystyle=|\mathcal{T}(q,\eta)|+\epsilon_{T}$ (17)
$\displaystyle\epsilon_{T}$ $\displaystyle\equiv\int
dq\rho(q)|\mathcal{T}(q)|$ (18)

For LHA to be valid for the matrix element $\rho(q,\eta)$, we need
$\mathcal{T}(q,\eta)\ll 1$. We also need $\epsilon_{T}\ll 1$, which translates
into,

$|V_{3}(q)|<<\frac{3!}{\beta\rho(q)\left|\Delta(q)^{3}+3\mathcal{E}(q)^{2}\Delta(q)\right|}$
(19)

If, for a system, Eq. 19 is invalid in some region, this does not make LHA
totally inapplicable. As an example, if it is valid in two regions of space, A
and B, then LHA can still predict the ratio of populations
${\rho_{A}}/{\rho_{B}}$ in these two regions.

Finally, for an off-diagonal density matrix element $\rho(q,\eta)$, the value
of the classical deviation, $\Delta(q,\eta)$, is at least as large as $\eta$,
and hence as the value of $\eta$ is increased, a stronger constraint on
$V_{3}(q)$ is imposed in Eq. 19. This sets a limit on the length scale on
which LHA can accurately capture spatial coherences.

We summarize by commenting that each term in the action in Eq. 3 gives rise to
a regime in which LHA is accurate, i.e. the kinetic, dissipative and potential
terms each give the high temperature, large coupling and the ‘slowly varying
potential’ limit, respectively.

##  III Application

###  III.1 Quartic Oscillator

Figure 2: The 2nd and 4th cumulant as a function of the system-bath coupling
parameter, $\gamma$, for anharmonicity (a) $a=2.5\times
10^{-3}(E_{h}/a_{0}^{4})$ and (b) $a=5\times 10^{-2}(E_{h}/a_{0}^{4})$. Orange
and blue curves indicate results obtained using TEMPO and LHA, respectively,
while the horizontal black lines indicate the exact results for $\gamma=0$ and
$\gamma=\infty$. The vertical dashed line corresponds to the point to the
right of which $\epsilon_{\kappa_{2/4}}<0.1$ [Eq. 22, Eq. 23]. Convergence for
TEMPO calculations for both (a) and (b) at $\gamma=0.5E_{h}$ was observed for
the number of path integral timesteps equal to 60, the system truncation
dimension of 9, and the TEMPO SVD precision parameter of $10^{-9}$.

We are now going to use LHA to study various CV systems. Let us first study a
quartic oscillator whose potential is given by,

$V(q)=\frac{1}{2}m\omega^{2}q^{2}+aq^{4}$ (20)

A quartic potential is the simplest bounded potential that provides a
deviation from the quadratic HO potential, latter being the case for which LHA
works exactly. A quartic oscillator has been studied in the context of quantum
to classical transition Oliveira2006Quantum-classical , quantum chaos
B.C.Bag1998A ; Jt1992Semiclassical ; Chakrabarty2019Out  and as a model of
bath oscillators interacting with a system Montenegro2014Nonlinearity . Here,
we study how the non-Gaussianity that naturally emerges in the MFGS of a
quartic oscillator changes as the system bath coupling is increased.

We use an Ohmic spectral density of exponential cutoff form (as some integrals
required to do the TEMPO calculations done in this section are analytically
known in this cutoff),

$\displaystyle
J(\omega)=\frac{2m\gamma}{\pi}\omega\exp\left(-\frac{\omega}{\Omega}\right)$
(21)

Using Hartree atomic units, we set $\omega=1E_{h}$, $\Omega=5E_{h}$,
$m=1m_{e}$ and $k_{B}T=0.5E_{h}$.

To characterize the MFGS, we calculate the second and fourth order cumulants
of position, i.e. $\kappa_{2}=\langle q^{2}\rangle$ and $\kappa_{4}=\langle
q^{4}\rangle-3\langle q^{2}\rangle^{2}$, respectively. The former of these
quantifies the width of the spatial probability density while the latter
quantifies how non-Gaussian the state is, since all cumulants of order three
and higher vanish for Gaussian states. The system is symmetric around $q=0$ so
the odd-order cumulants of $q$ vanish identically, regardless of whether the
state is Gaussian.

To benchmark our results we use the analytically known weak and USC results,
as well as the numerically exact method, TEMPO strathearn2018efficient ;
strathearn2020modelling ; chiu2022numerical , applicable in weak and
intermediate coupling regime.

The TEMPO is a tensor network technique to efficiently evaluate the path
integral given in Eq. 3. It has been previously used to calculate the MFGS for
discrete systems chiu2022numerical . Here we use it to study a CV system like
a particle in a generic 1D potential $V(q)$. Although the system in question
is infinite dimensional, at any given temperature and coupling, it can be
approximated as an effectively finite dimensional system by truncating the
system Hilbert space in the energy eigenbasis topaler1993system . The value of
this effective dimension can be determined through a convergence test. In this
case for example, for $\gamma=0.5\;E_{h}$, we needed to treat the system as an
effectively 9 dimensional discrete system to get sufficient convergence.

Figure 2 (a) and (b) plot both the cumulants, $\kappa_{2}$ and $\kappa_{4}$,
calculated using LHA, as a function of coupling strength for $a=2.5\times
10^{-3}$ and $5\times 10^{-2}$ [Eq. 20], respectively. The vertical red line
marks the value of $\gamma$ such that for $\gamma$ larger than this, the
associated relative error $\epsilon_{\kappa_{2}},\epsilon_{\kappa_{4}}<0.1$,
where, from Eq. 16 we have,

$\displaystyle\epsilon_{\kappa_{2}}$ $\displaystyle=\epsilon_{\hat{q}^{2}}$
(22) $\displaystyle\epsilon_{\kappa_{4}}$
$\displaystyle=\frac{\braket{\hat{q}^{4}}\epsilon_{\hat{q}^{4}}+6\braket{\hat{q}^{2}}\epsilon_{\hat{q}^{2}}}{\kappa_{4}}$
(23)

As expected, the accuracy of LHA falls as the value of $\gamma$ is decreased
from infinity while for TEMPO, the computational resources required increase
as the value of $\gamma$ is increased from the weak coupling regime. Because
of computational constraints, here we only show the TEMPO results up to
$\gamma\leq 0.5$. In Figure 2 (b), we observe a nontrivial rise in the TEMPO
curve for $\kappa_{4}$ as $\gamma$ is increased, followed by an eventual decay
to the USC value as suggested by the LHA curve, suggesting a peak at about
$\gamma=10$.

The results for MFGS in weak and USC limit in Figure 2 (a) and (b) suggests
that the non-gaussianity, that naturally emerges in a quartic oscillator, may
either get enhanced or decay as $\gamma$ is increased, depending on the value
of the parameter $a$ used in the potential $V(q)$ [Eq. 20]. TEMPO and LHA
complement each other to calculate the $\gamma$ dependence of $\kappa_{4}$ in
the low and high $\gamma$ limit, respectively. Note that non-Gaussianity is
considered a resource for quantum information processing in CV systems
Takagi2018Convex ; weedbrook2012gaussian , and LHA and TEMPO can help quantify
how the non-Gaussianity depends upon the system-bath coupling strength for a
CV system.

###  III.2 Asymmetric double well potential

Let us now discuss an asymmetric quartic DW potential of the form,

$V(q)=a_{4}q^{4}-a_{2}q^{2}+a_{1}q$ (24)

where we set $a_{4}=1/2=a_{2}$. In this and following sections, we use an
ohmic spectral density of Drude-Lorentz form,

$\displaystyle
J(\omega)=\frac{2m\gamma}{\pi}\omega\frac{\Omega^{2}}{\omega^{2}+\Omega^{2}}$
(25)

as this form of spectral density admits a closed form expression for the LHA
weiss2012quantum . We again set $\omega=1E_{h}$, $\Omega=5E_{h}$, $m=1m_{e}$
and $k_{B}T=0.5E_{h}$.

As the accuracy of LHA increases with increasing system bath coupling, we will
instead test its applicability in weak coupling limit. At $\gamma=0$, we are
interested in the LHA value of the probability, $\mathcal{P}$, of finding the
particle in the right well, given as,

$\mathcal{P}=\int_{q_{b}}^{\infty}\rho(q)dq$ (26)

where, $q_{b}$ represents the position of the peak of the barrier. We also
want to verify the associated estimated relative error for LHA [Eq. 16]
against the analytically derivable value for the same at $\gamma=0$.

Since LHA will work at lower coupling only if $V(q)$ varies slowly, let us
define,

$\bar{V}_{b}(q)\equiv b\times V(q/b)$ (27)

where for a fixed value of $b>1$, $\bar{V}_{b}(q)$ is just a new rescaled
potential with smaller values for the higher derivatives with respect to $q$.
Note that $\bar{V}_{b}(q)$ has been mathematically constructed just to verify
the validity of the error estimate [Eq. 16].

Figure 3: Plot for relative error (solid lines) [Eq. 28] for LHA value of the
Gibbs state probability $\mathcal{P}$ [Eq. 26] for the particle to be found in
the right well of an asymmetric double well potential $\bar{V}_{b}(q)$ [Eq.
27], as a function of $b$, a dimensionless parameter of $\bar{V}_{b}(q)$. The
dotted lines represent the estimated relative error for the same [Eq. 16].
Plots are for $a_{1}=0.1$ (red lines) and $0.3$ (green lines), where $a_{1}$
quantifies the degree of asymmetry in $\bar{V}_{b}(q)$. Note that both the
dotted lines are almost overlapping.

At a given value of $b$, let $\rho_{b}$ represent the LHA state and
$\sigma_{b}=Z_{b}^{-1}\exp\\{-\beta H_{0}\\}$ represent the Gibbs state with
$\mathcal{P}_{l}(b)$ and $\mathcal{P}_{g}(b)$ being the corresponding values
of the probability of finding the particle in the right well, respectively.
Here, $Z_{b}$ is an appropriate normalization factor. Then, let the estimated
and analytical value of the relative error be denoted as $\epsilon_{l}(b)$
[Eq. 16] and $\epsilon_{g}(b)$, respectively, where $\epsilon_{g}(b)$ is
defined as,

$\displaystyle\epsilon_{g}(b)$
$\displaystyle=\frac{|\mathcal{P}_{l}(b)-\mathcal{P}_{g}(b)|}{\mathcal{P}_{g}(b)}$
(28)

In Figure 3, we plot $\epsilon_{l}(b)$ and $\epsilon_{g}(b)$ as a function of
$b$ for $a_{1}=0.1$ and $0.3$ [Eq. 24]. As expected, as the value of $b$
increases and the potential becomes smoother, both $\epsilon_{l}(b)$ and
$\epsilon_{g}(b)$ decrease in value. We find that $\epsilon_{l}(b)$ provides
an upper bound on $\epsilon_{g}(b)$, roughly overestimating it by about two
orders of magnitude in this case. One possible explanation for this could be
that the actual error associated with the expectation value of an observable
$\hat{O}$ can get canceled out during the act of taking the the trace,
$\operatorname{Tr}{\left(\hat{O}\rho\right)}$.

Figure 4: The left y-axis plots the LHA value for the probability
$\mathcal{P}$ (solid blue line) for a particle in an asymmetric double well
potential $\bar{V}_{b}(q)$ [Eq. 27] to be found in the right well [Eq. 26], as
a function of coupling parameter $\gamma$, for $b=4$. The dashed blue and red
lines mark the value of $\mathcal{P}$ corresponding to $\gamma=0$ and
$\infty$, respectively. The right y-axis plots the estimated relative error
for the corresponding value of $\mathcal{P}$ [Eq. 16].

In Figure 4, we fix $b=4$ and plot $\mathcal{P}_{l}(b=4)$ as a function of
$\gamma$. We note that throughout the range of $\gamma$, the relative error,
$\epsilon_{l}$, is always less than or equal to $0.07$, indicating that LHA
can accurately calculate the value of $\mathcal{P}$ for this DW system
throughout the range of the coupling parameter. Note that for this DW system,
the values for the distance between the wells, the barrier height from the
potential minima and the energy difference between the well minimas are given
as $5.6a_{0}$, $0.81E_{h}$ and $0.56E_{h}$, respectively. Hence, although
these parameters are not extreme, yet we note that this DW potential was
constructed to be slowly varying with respect to $q$ so that LHA would give
accurate results for it, and hence at best, provides an evidence for the
accuracy of LHA when $V(q)$ varies slowly.

On the right y-axis of Figure 4, we plot the relative error estimate
$\epsilon_{l}(\gamma)$. We note that, as expected, $\epsilon_{l}(\gamma)$
monotonically decreases with increasing value of $\gamma$. Note that for
$\gamma\lesssim 0.1E_{h}$, the value of $\epsilon_{l}(\gamma)$ plateaus,
suggesting that in low $\gamma$ limit, the kinetic term in the action [Eq. 3]
dominates over the dissipative term, and is the sole ($\gamma$ independent)
factor that is still responsible for the accuracy of LHA.

###  III.3 Proton tunneling in DNA

Slocombe et al. slocombe2022open ; slocombe2022quantum  performed an open
quantum system analysis of a proton in Guanine-Cytosine (G-C) pair in DNA to
study the probability for it to tunnel to a mutated state. They treated the
problem as a particle in an asymmetric double well potential $V(q)$ (see
Figure 1), interacting with a non-interacting bosonic bath through an Ohmic
spectral density [Eq. 25]. Here, $V(q)$ is taken to be a back-to-back Morse
potential,

$V(q)=v_{1}\left(e^{-2a_{1}(q-r_{1})}-2e^{-a_{1}(q-r_{1})}\right)\\\
+v_{2}\left(e^{-2a_{2}(r_{2}-q)}-2e^{-a_{2}(r_{2}-q)}\right)$ (29)

where the lower and higher well represent the canonical and mutated state for
the proton to be in, respectively. Here, in Hartree atomic units,
$v_{1}=0.1617E_{h}$, $v_{2}=0.082E_{h}$, $a_{1}=0.305a_{0}^{-1}$,
$a_{2}=0.755a_{0}^{-1}$, $r_{1}=-2.7a_{0}$ and $r_{2}=2.1a_{0}$. The
parameters of the spectral density [Eq. 25] was taken to be,
$\gamma=0.018E_{h}\equiv\gamma_{p}$ and $k_{B}T=0.00095E_{h}$ while an
infinite bath cutoff was assumed as a Markovian approximation (i.e.,
$\Omega=\infty E_{h}$) slocombe2022quantum . Here, $\gamma_{p}$ denotes the
biologically relevant value of $\gamma$ for this system. They calculated the
steady state for this system using Caldeira Leggett master equation (CLME) and
found the probability for the proton to tunnel to the mutated state, defined
as

$\mathcal{M}=\int_{q_{b}}^{\infty}\rho(q)dq$ (30)

to be of the order of $10^{-4}$, which, as a point of comparison, is several
orders of magnitude higher than the zero coupling value of about $2.7\times
10^{-8}$. Here $q_{b}$ denotes the position of the barrier maxima.

We remark that CLME is a weak coupling master equation, which requires
$\gamma_{p}<<K_{B}T,\Omega$. But, for this system, $\gamma_{p}>K_{B}T$ and
hence the system is well outside the weak coupling regime. This motivates us
to study the MFGS of this system using LHA.

In Figure 5, we plot LHA value for $\mathcal{M}$ as a function of $\gamma$
along with the estimated relative error, $\epsilon_{\mathcal{M}}$ [Eq. 16],

$\displaystyle\epsilon_{\mathcal{M}}=\frac{1}{\mathcal{M}}\int_{q_{b}}^{\infty}\epsilon(q)\rho(q)dq$
(31)

We note that at $\gamma=\gamma_{p}$, $\epsilon_{\mathcal{M}}\approx 0.19$.
while for $\gamma\geq 10\gamma_{p}$, $\epsilon_{\mathcal{M}}\ll 1$. Hence, in
the region where LHA is applicable ($\gamma\geq 10\gamma_{p}$), as the value
of $\gamma$ is lowered from infinity, we see the value of $\mathcal{M}$
reducing from the USC value, suggesting that $\mathcal{M}=10^{-4}$ is
unlikely. But at $\gamma=\gamma_{p}$, we can’t speculate on the true value of
$\mathcal{M}$ since the error is high.

This disagreement could have a couple of possible reasons. Note that in the
intermediate coupling regime, as in this case, it is still unclear whether
MFGS is actually the steady state of the system or not cresser2021weak ;
trushechkin2022open . On the other hand, CLME, being a weak coupling master
equation, could be responsible for this disagreement.

Figure 6 presents a plot of $\mathcal{M}$ as a function of temperature $T$,
across a temperature range of approximately $300K\pm 20\%$. As expected, we
find that $\mathcal{M}$ rises roughly exponentially with the value of $T$. We
also plot the value of $\mathcal{M}$ corresponding to $\gamma=0$ and
$\gamma=\infty$, finding it to closely follow the same trend. At the
temperature of about $87^{\circ}\mathrm{C}$, we find an approximately 10 times
rise in the value of $\mathcal{M}$ as compared with corresponding room
temperature value. The accuracy of LHA is observed to increase with a rise in
temperature, as expected [Appendix. 2], as observed from the decay in the plot
for the relative error $\epsilon_{\mathcal{M}}$ with increasing temperature.

Figure 5: Left y-axis plots the LHA value for mutation probability of the
proton, $\mathcal{M}$ (solid blue line), as a function of the coupling
parameter $\gamma$. The dashed blue and red lines mark the value of
$\mathcal{M}$ corresponding to $\gamma=0$ and $\infty$, respectively. The
vertical orange dotted line marks $\gamma_{p}=0.018E_{h}$ (the biologically
relevant value of $\gamma$). The green dotted curve corresponding to the right
hand side y-axis quantifies the estimated relative error
($\epsilon_{\mathcal{M}}$ [Eq. 31]) in the LHA value for $\mathcal{M}$. Figure
6: Left y-axis plots the LHA value for the mutation probability,
$\mathcal{M}$, as a function of temperature $T$ (in Kelvin) for biologically
relevant value of the coupling parameter ($\gamma_{p}=0.018\;E_{h}$) (blue
scatter points). The dashed blue and red lines mark the value of $\mathcal{M}$
corresponding to $\gamma=0$ and $\infty$, respectively. The vertical orange
dotted line marks $T=300K$ (i.e., room temperature). The green dotted curve
corresponding to the right hand side y-axis quantifies the estimated relative
error ($\epsilon_{\mathcal{M}}$ [Eq. 31]) in the LHA value for $\mathcal{M}$.

##  IV Conclusion

For a particle in a 1D potential $V(q)$, we derived an approximate expression
for the MFGS [Eq. 1]. This method is accurate when either the system-bath
coupling or the temperature is large, or when the third and higher derivatives
of the potential are small compared to certain system-bath specific
parameters. This result thus has applications in a wide variety of problems in
quantum thermodynamics, technology, chemistry and quantum biology.

We have found that for a CV system, the USC cresser2021weak  and high
temperature timofeev2022hamiltonian  results obtained recently are a special
case of LHA [Eq. 12]. As a future direction, it would be interesting to apply
a similar approach to use the analytically known exact expression for the
master equation for a harmonic oscillator hu1992quantum  to derive an
approximate master equation for a particle in generic potential $V(q)$.
Moreover, determining the origin of mismatch between MFGS and the steady state
expression for the mutated population $\mathcal{M}$ for the proton tunneling
problem considered in Sec. III.3 will be of interest slocombe2022open .

Finally, different ways to improve on LHA and get better error bounds for it
will enhance its applicability to physical problems.

##  V Acknowledgement

I would like to thank my PhD supervisor, Dr. Sibasish Ghosh, for his comments,
suggestions and several extended discussions. I also thank Aidan Strathearn,
who has contributed significantly to the preparation of this manuscript.

## References

  * (1) S. L. Braunstein and P. Van Loock, “Quantum information with continuous variables,” _Reviews of modern physics_ , vol. 77, no. 2, p. 513, 2005.
  * (2) T. Ralph, “Quantum optical systems for the implementation of quantum information processing,” _Reports on Progress in Physics_ , vol. 69, pp. 853 – 898, 2006.
  * (3) I. Djordjevic, “Hybrid cv-dv quantum communications and quantum networks,” _IEEE Access_ , vol. 10, pp. 23 284–23 292, 2022.
  * (4) A. Inoue, T. Kashiwazaki, T. Yamashima, N. Takanashi, T. Kazama, K. Enbutsu, K. Watanabe, T. Umeki, M. Endo, and A. Furusawa, “Toward a multi-core ultra-fast optical quantum processor: 43-ghz bandwidth real-time amplitude measurement of 5-db squeezed light using modularized optical parametric amplifier with 5g technology,” _Applied Physics Letters_ , 2022.
  * (5) X. Yang and M. Xiao, “Electromagnetically induced entanglement,” _Scientific Reports_ , vol. 5, 2015.
  * (6) C. Gross, H. Strobel, E. Nicklas, T. Zibold, N. Bar-Gill, G. Kurizki, and M. Oberthaler, “Atomic homodyne detection of continuous-variable entangled twin-atom states,” _Nature_ , vol. 480, pp. 219–223, 2011.
  * (7) J. Huang, Y. Liu, J.-F. Huang, and J. Liao, “Generation of macroscopic entangled cat states in a longitudinally coupled cavity-qed model,” _Physical Review A_ , 2019.
  * (8) R. Takagi and Q. Zhuang, “Convex resource theory of non-gaussianity,” _Physical Review A_ , 2018.
  * (9) G. Ferrari, L. Lami, T. Theurer, and M. Plenio, “Asymptotic state transformations of continuous variable resources,” _Communications in Mathematical Physics_ , vol. 398, pp. 291 – 351, 2020.
  * (10) V. Narasimhachar, S. Assad, F. Binder, J. Thompson, B. Yadin, and M. Gu, “Thermodynamic resources in continuous-variable quantum systems,” _npj Quantum Information_ , vol. 7, pp. 1–7, 2019.
  * (11) E. Hoskinson, F. Lecocq, N. Didier, A. Fay, F. Hekking, W. Guichard, O. Buisson, R. Dolata, B. Mackrodt, and A. Zorin, “Quantum dynamics in a camelback potential of a dc squid.” _Physical review letters_ , vol. 102 9, p. 097004, 2008.
  * (12) R. Ferreira and G. Bastard, “Tunnelling and relaxation in semiconductor double quantum wells,” _Reports on Progress in Physics_ , vol. 60, pp. 345–387, 1997.
  * (13) Monteoliva and Paz, “Decoherence and the rate of entropy production in chaotic quantum systems,” _Physical review letters_ , vol. 85 16, pp. 3373–6, 2000.
  * (14) B.C.Bag, S.Chaudhuri, J. Chaudhuri, and D.S.Ray, “A semiclassical theory of quantum noise in open chaotic systems,” _Physica D: Nonlinear Phenomena_ , vol. 125, pp. 47–64, 1998.
  * (15) A. Jt, “Semiclassical chaos in quartic anharmonic oscillators.” _Physical Review A_ , vol. 45, pp. 5373–5377, 1992.
  * (16) B. Chakrabarty and S. Chaudhuri, “Out of time ordered effective dynamics of a quartic oscillator,” _SciPost Physics_ , 2019.
  * (17) K.-L. Liu and H. Goan, “Non-markovian entanglement dynamics of quantum continuous variable systems in thermal environments,” _Physical Review A_ , vol. 76, p. 022312, 2007.
  * (18) R. Vasile, S. Maniscalco, M. Paris, H. Breuer, and J. Piilo, “Quantifying non-markovianity of continuous-variable gaussian dynamical maps,” _Physical Review A_ , vol. 84, p. 052118, 2011.
  * (19) C. Hörhammer and H. Buettner, “Environment-induced two-mode entanglement in quantum brownian motion,” _Physical Review A_ , vol. 77, p. 042305, 2007.
  * (20) J. Ruostekoski and D. F. Walls, “Bose-einstein condensate in a double-well potential as an open quantum system,” _Physical Review A_ , vol. 58, no. 1, p. R50, 1998.
  * (21) M. P. Fisher and A. T. Dorsey, “Dissipative quantum tunneling in a biased double-well system at finite temperatures,” _Physical review letters_ , vol. 54, no. 15, p. 1609, 1985.
  * (22) H. Grabert and U. Weiss, “Quantum tunneling rates for asymmetric double-well systems with ohmic dissipation,” _Physical review letters_ , vol. 54, no. 15, p. 1605, 1985.
  * (23) M. Topaler and N. Makri, “Quantum rates for a double well coupled to a dissipative bath: Accurate path integral results and comparison with approximate theories,” _The Journal of chemical physics_ , vol. 101, no. 9, pp. 7500–7519, 1994.
  * (24) A. C. Oliveira, J. G. P. de Faria, and M. Nemes, “Quantum-classical transition of the open quartic oscillator: The role of the environment.” _Physical review. E, Statistical, nonlinear, and soft matter physics_ , vol. 73 4 Pt 2, p. 046207, 2006.
  * (25) J. Cresser and J. Anders, “Weak and ultrastrong coupling limits of the quantum mean force gibbs state,” _Physical Review Letters_ , vol. 127, no. 25, p. 250601, 2021.
  * (26) U. Seifert, “First and second law of thermodynamics at strong coupling,” _Physical review letters_ , vol. 116, no. 2, p. 020601, 2016.
  * (27) T. G. Philbin and J. Anders, “Thermal energies of classical and quantum damped oscillators coupled to reservoirs,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 49, no. 21, p. 215303, 2016.
  * (28) C. Jarzynski, “Nonequilibrium work theorem for a system strongly coupled to a thermal environment,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2004, no. 09, p. P09005, 2004.
  * (29) M. Campisi, P. Talkner, and P. Hänggi, “Fluctuation theorem for arbitrary open quantum systems,” _Physical review letters_ , vol. 102, no. 21, p. 210401, 2009.
  * (30) T. W. Allen, O. S. Andersen, and B. Roux, “Molecular dynamics—potential of mean force calculations as a tool for understanding ion permeation and selectivity in narrow channels,” _Biophysical chemistry_ , vol. 124, no. 3, pp. 251–267, 2006.
  * (31) K. Maksimiak, S. Rodziewicz-Motowidło, C. Czaplewski, A. Liwo, and H. A. Scheraga, “Molecular simulation study of the potentials of mean force for the interactions between models of like-charged and between charged and nonpolar amino acid side chains in water,” _The Journal of Physical Chemistry B_ , vol. 107, no. 48, pp. 13 496–13 504, 2003.
  * (32) B. Roux and T. Simonson, “Implicit solvent models,” _Biophysical chemistry_ , vol. 78, no. 1-2, pp. 1–20, 1999.
  * (33) B. Roux, “The calculation of the potential of mean force using computer simulations,” _Computer physics communications_ , vol. 91, no. 1-3, pp. 275–282, 1995.
  * (34) M. Gillan, “Quantum-classical crossover of the transition rate in the damped double well,” _Journal of Physics C: Solid State Physics_ , vol. 20, no. 24, p. 3621, 1987.
  * (35) Y. Liu, Y. Yan, T. Xing, and Q. Shi, “Understanding the large kinetic isotope effect of hydrogen tunneling in condensed phases by using double-well model systems,” _The Journal of Physical Chemistry B_ , vol. 125, no. 22, pp. 5959–5970, 2021.
  * (36) L. Slocombe, M. Sacchi, and J. Al-Khalili, “An open quantum systems approach to proton tunnelling in dna,” _Communications Physics_ , vol. 5, no. 1, p. 109, 2022.
  * (37) L. Slocombe, “Quantum effects in dna replication fidelity,” Ph.D. dissertation, University of Surrey, 2022.
  * (38) J. Zhang, R. Borrelli, and Y. Tanimura, “Proton tunneling in a two-dimensional potential energy surface with a non-linear system–bath interaction: Thermal suppression of reaction rate,” _The Journal of Chemical Physics_ , vol. 152, no. 21, 2020.
  * (39) J. P. Bothma, J. B. Gilmore, and R. H. McKenzie, “The role of quantum effects in proton transfer reactions in enzymes: quantum tunneling in a noisy environment?” _New Journal of Physics_ , vol. 12, no. 5, p. 055002, 2010.
  * (40) H. P. Breuer and F. Petruccione, _The theory of open quantum systems_. Great Clarendon Street: Oxford University Press, 2002.
  * (41) A. Trushechkin, “Quantum master equations and steady states for the ultrastrong-coupling limit and the strong-decoherence limit,” _Physical Review A_ , vol. 106, no. 4, p. 042209, 2022.
  * (42) C. L. Latune, “Steady state in ultrastrong coupling regime: perturbative expansion and first orders,” _arXiv preprint arXiv:2110.02186_ , 2021.
  * (43) G. Timofeev and A. Trushechkin, “Hamiltonian of mean force in the weak-coupling and high-temperature approximations and refined quantum master equations,” _International Journal of Modern Physics A_ , vol. 37, no. 20n21, p. 2243021, 2022.
  * (44) A. Gelzinis and L. Valkunas, “Analytical derivation of equilibrium state for open quantum system,” _The Journal of Chemical Physics_ , vol. 152, no. 5, 2020.
  * (45) H. Grabert, P. Schramm, and G.-L. Ingold, “Quantum brownian motion: The functional integral approach,” _Physics reports_ , vol. 168, no. 3, pp. 115–207, 1988.
  * (46) U. Weiss, _Quantum dissipative systems_. World Scientific, 2012.
  * (47) Y.-F. Chiu, A. Strathearn, and J. Keeling, “Numerical evaluation and robustness of the quantum mean-force gibbs state,” _Physical Review A_ , vol. 106, no. 1, p. 012204, 2022.
  * (48) A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, “Efficient non-markovian quantum dynamics using time-evolving matrix product operators,” _Nature communications_ , vol. 9, no. 1, p. 3322, 2018.
  * (49) A. Strathearn, _Modelling non-Markovian quantum systems using tensor networks_. Springer Nature, 2020.
  * (50) R. P. Feynman and F. Vernon Jr, “The theory of a general quantum system interacting with a linear dissipative system,” _Annals of physics_ , vol. 281, no. 1-2, pp. 547–607, 2000.
  * (51) N. Makri and D. E. Makarov, “Tensor propagator for iterative quantum time evolution of reduced density matrices. i. theory,” _The Journal of chemical physics_ , vol. 102, no. 11, pp. 4600–4610, 1995.
  * (52) ——, “Tensor propagator for iterative quantum time evolution of reduced density matrices. ii. numerical methodology,” _The Journal of chemical physics_ , vol. 102, no. 11, pp. 4611–4618, 1995.
  * (53) Y. Tanimura, “Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (heom),” _The Journal of chemical physics_ , vol. 153, no. 2, 2020.
  * (54) M. Topaler and N. Makri, “System-specific discrete variable representations for path integral calculations with quasi-adiabatic propagators,” _Chemical physics letters_ , vol. 210, no. 4-6, pp. 448–457, 1993.
  * (55) G. Barton, “Quantum mechanics of the inverted oscillator potential,” _Annals of Physics_ , vol. 166, no. 2, pp. 322–363, 1986.
  * (56) V. Subramanyan, S. S. Hegde, S. Vishveshwara, and B. Bradlyn, “Physics of the inverted harmonic oscillator: From the lowest landau level to event horizons,” _Annals of Physics_ , vol. 435, p. 168470, 2021.
  * (57) S. Baskoutas, A. Jannussis, and R. Mignani, “Quantum tunnelling of a damped and driven, inverted harmonic oscillator,” _Journal of Physics A: Mathematical and General_ , vol. 26, no. 23, p. 7137, 1993.
  * (58) V. Montenegro, A. Ferraro, and S. Bose, “Nonlinearity induced entanglement stability in a qubit-oscillator system,” _Physical Review A_ , vol. 90, p. 013829, 2014.
  * (59) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” _Reviews of Modern Physics_ , vol. 84, no. 2, p. 621, 2012.
  * (60) A. Trushechkin, M. Merkli, J. Cresser, and J. Anders, “Open quantum system dynamics and the mean force gibbs state,” _AVS Quantum Science_ , vol. 4, no. 1, 2022.
  * (61) B. L. Hu, J. P. Paz, and Y. Zhang, “Quantum brownian motion in a general environment: Exact master equation with nonlocal dissipation and colored noise,” _Physical Review D_ , vol. 45, no. 8, p. 2843, 1992.
  * (62) R. Feynman and A. Hibbs, “The path integral formulation of quantum mechanics,” _McGraw-Hiill, New York_ , 1965.
  * (63) W. Sherman, “Summation of harmonics with random phase angles,” in _Proceedings of the Institution of Electrical Engineers_ , vol. 119, no. 11. IET, 1972, pp. 1643–1648.

##  Appendix A Summary of main results

In this section, we will summarize the basic results of this paper. For
details, see the sections mentioned here.

The path integral representation for the MFGS,
$\rho(q,\eta)=\braket{q+\eta}{\hat{\rho}}{q-\eta}$, is given as,

$\rho(q,\eta)=Z_{R}^{-1}\int_{q+\eta}^{q-\eta}\mathcal{D}q(\tau)e^{-S[q(\tau)]}$
(32)

Parameterising the path variable in terms of its deviation from the constant
path, $q(\tau)=q+\delta(\tau)$, we expand the potential to second order in
$\delta$ [Eq. 9], resulting in Eq. 32 becoming,

$\displaystyle\rho(q,\eta)$ $\displaystyle=Z_{R}^{-1}\exp\left(-\beta
V(q)+\beta\frac{1}{2}m\omega_{q}^{2}Q_{q}^{2}\right)\int_{\eta}^{-\eta}\mathcal{D}\delta(\tau)e^{-S_{q}[\delta(\tau)+Q_{q}]}$
(33)

where $S_{q}[\delta(\tau)+Q_{q}]$ is the action for a HO with frequency
$\omega_{q}$ [Eq. 10]. Reparametrize $\delta(\tau)+Q_{q}\to\delta(\tau)$ to
get,

$\displaystyle\int_{\eta}^{-\eta}\mathcal{D}\delta(\tau)e^{-S_{q}[\delta(\tau)+Q_{q}]}$
$\displaystyle=\int_{Q_{q}+\eta}^{Q_{q}-\eta}\mathcal{D}\delta(\tau)e^{-S_{q}[\delta(\tau)]}$
(34) $\displaystyle=\tilde{Z}_{q}\tilde{\rho}_{q}(Q_{q},\eta)$ (35)

where, $\tilde{\rho}_{q}(q,\eta)$ is the MFGS for a HO with frequency
$\omega_{q}$. That is, the functional integral in this equation is identical
to that of a HO but for the reduced HO partition function $\tilde{Z}_{q}$.
This leads us to [Eq. 12],

$\rho(q,\eta)=Z_{R}^{-1}\frac{\exp\left(-\beta
V(q)\right)}{\exp\left(-\beta\frac{1}{2}m\omega_{q}^{2}Q_{q}^{2}\right)}\tilde{Z}_{q}\tilde{\rho}_{q}(Q_{q},\eta)$
(36)

We can write the expression for LHA more explicitly as,

$\displaystyle\rho(q,\eta)$
$\displaystyle=Z^{-1}\mathcal{G}(q)\exp\left(-\frac{1}{2\braket{\hat{X}^{2}(q)}}\left(\frac{V_{1}(q)}{V_{2}(q)}\right)^{2}-2\braket{\hat{P}^{2}(q)}\eta^{2}-\beta\left(V(q)-\frac{V_{1}(q)^{2}}{2V_{2}(q)}\right)\right)$
(37)

Here, $\braket{\hat{X}^{2}(q)}$ and $\braket{\hat{P}^{2}(q)}$ are the
corresponding effective $\braket{\hat{q}^{2}}$ and $\braket{\hat{p}^{2}}$
values of the harmonic potential approximated at point $q$, whose value is
given as weiss2012quantum ; grabert1988quantum ,

$\displaystyle\braket{\hat{X}^{2}(q)}$
$\displaystyle=\frac{1}{m\beta}\sum_{n=-\infty}^{+\infty}\frac{1}{\mathcal{A}_{n}(q)}$
(38) $\displaystyle\braket{\hat{P}^{2}(q)}$
$\displaystyle=\frac{m}{\beta}\sum_{n=-\infty}^{+\infty}\frac{\omega_{q}^{2}+\zeta_{n}}{\mathcal{A}_{n}(q)}$
(39) $\displaystyle\mathcal{A}_{n}(q)$
$\displaystyle=\omega_{q}^{2}+\nu_{n}^{2}+\zeta_{n}$ (40)
$\displaystyle\nu_{n}$ $\displaystyle\equiv 2\pi nk_{B}T$ (41)
$\displaystyle\zeta_{n}$
$\displaystyle\equiv\frac{1}{m}\int_{0}^{\infty}d\omega^{\prime}\frac{J(\omega^{\prime})}{\omega^{\prime}}\frac{2\nu_{n}^{2}}{\omega^{\prime
2}+\nu_{n}^{2}}$ (42)

To arrive at Eq. 37, we did the following replacement in Eq. 36. The
expression for the reduced partition function, $\tilde{Z}_{q}$, is given as
[Appendix. 1],

$\displaystyle\tilde{Z}_{q}$
$\displaystyle=\mathcal{I}(q)\sqrt{2\pi\braket{\hat{X}^{2}(q)}}$ (43)

where $\mathcal{I}(q)$ is defined as,

$\displaystyle\mathcal{I}(q)$
$\displaystyle\equiv\int_{0}^{0}\mathcal{D}\delta(\tau)e^{-S_{q}[\delta(\tau)]}$
(44) $\displaystyle=\frac{Z_{R}}{Z}\mathcal{G}(q)$ (45)

Here, $Z$ [Eq. 37] is an appropriate constant independent of $q$, and hence is
a normalization factor common for all the effective quadratic approximations
done at all points $q$. This would imply, for example, that $Z$ can depend
upon $k_{B}T$, $\Omega$ and $\gamma$ but not on $\omega_{q}$ because the
latter explicitly depends upon $V(q)$ by definition [Eq. 10]. The non-trivial
part is to then show that $\mathcal{G}(q)$ is given as [Appendix. 2],

$\displaystyle\mathcal{G}(q)$
$\displaystyle=\frac{1}{\sqrt{m\omega_{q}^{2}\beta\braket{\hat{X}^{2}(q)}}}\prod_{n=1}^{\infty}\frac{\nu_{n}^{2}+\zeta_{n}}{\mathcal{A}_{n}(q)}$
(46)

Furthermore, the expression for $\tilde{\rho}_{q}(q,\eta)$ is given as
weiss2012quantum ; grabert1988quantum ,

$\displaystyle\tilde{\rho}_{q}(q,\eta)$
$\displaystyle=\frac{1}{\sqrt{2\pi\braket{\hat{X}^{2}(q)}}}\exp{\left(-\frac{q^{2}}{2\braket{\hat{X}^{2}(q)}}-2\braket{\hat{P}^{2}(q)}\eta^{2}\right)}$
(47)

With these replacements, we can move from Eq. 36 to Eq. 37. We note that, for
a specific spectral density, like Ohmic spectral density with Drude-Lorentz
cutoff [Eq. 25], the expressions for Eq. 38, Eq. 39 and Eq. 46 can be
evaluated in closed form weiss2012quantum .

Note that when $\omega_{q}^{2}<0$, $\braket{\hat{X}^{2}(q)}$ and
$\braket{\hat{P}^{2}(q)}$, as defined in Eq. 38 and Eq. 39, can no more be
interpreted as expectation values $\braket{\hat{q}^{2}}$ and
$\braket{\hat{p}^{2}}$ for the corresponding inverted HO, and must instead be
taken as mathematical objects as defined above. In Appendix. 2, we show that
as long as

$m\omega_{q}^{2}>-V_{P}$ (48)

where $V_{P}$ is defined in Sec. II, the unnormalized path integral for an
inverted HO,

$\displaystyle\tilde{Z}_{q}\tilde{\rho}_{q}(Q_{q},\eta)=\int_{Q_{q}+\eta}^{Q_{q}-\eta}\mathcal{D}\delta(\tau)e^{-S_{q}[\delta(\tau)]}$
(49)

is still convergent, which is what we require in the expression of LHA [Eq.
36]. We note here that $\tilde{\rho}_{q}(Q_{q},\eta)$ and $\tilde{Z}_{q}$ are
both individually undefined for any inverted HO, regardless of Eq. 48. Given
Eq. 48, to show that Eq. 49 converges, it is sufficient to show the
convergence of $\mathcal{G}(q)$, $\braket{\hat{X}^{2}(q)}^{-1}$ and
$\braket{\hat{P}^{2}(q)}$ that occur in the expression for LHA [Eq. 37].

The expression for the error estimate associated with LHA [Eq. 16] is derived
in Appendix. C. The expression for $\Delta(q,\eta)$ and $\mathcal{E}(q)$,
required during the error estimation, are given as [Appendix. 1],

$\displaystyle\Delta(q,\eta)=\max\left\\{\left|\frac{V_{1}(q)}{V_{2}(q)}\left(\frac{W(q)}{\braket{\hat{X}^{2}(q)}}-1\right)\right|,|\eta|\right\\}$
(50)
$\displaystyle\mathcal{E}(q)=\sqrt{\braket{\hat{X}^{2}(q)}-\frac{1}{V_{2}(q)\beta}}\left(\sqrt{\frac{1}{V_{2}(q)\beta\braket{\hat{X}^{2}(q)}}}+1\right)$
(51)

where,

$\displaystyle W(q)$
$\displaystyle=\frac{1}{m\beta}\sum_{n=-\infty}^{+\infty}\frac{(-1)^{n}}{\mathcal{A}_{n}(q)}$
(52)

Finally, LHA in USC, high temperature and infinite cutoff limit is studied in
Appendix. D.

##  Appendix B Expression for LHA

###  1 Harmonic oscillator partition function

In this section, we are going to derive the expression for the HO partition
function, $\tilde{Z}$ [Eq. 43] for a generic HO potential given as,

$V_{HO}(q)=\frac{1}{2}m\omega^{2}q^{2}$ (53)

Let us write any path $q(\tau)$ starting and ending at $q+\eta$ and $q-\eta$,
respectively, in terms of the classical path for HO and deviation over it,

$\displaystyle q(\tau)$ $\displaystyle=q_{\textbf{cl}}(\tau)+\delta(\tau)$
(54)

where $\delta(0)=0=\delta(\beta)$. Then, since the classical path is a local
extremum, the action given in Eq. 3 becomes feynman1965path ,

$\displaystyle S_{HO}[q(\tau)]$
$\displaystyle=S_{HO}[q_{\textbf{cl}}(\tau)]+S_{HO}[\delta(\tau)]$ (55)
$\displaystyle=S_{HO}[q_{\textbf{cl}}(\tau)]+\int_{0}^{\beta}d\tau\left(\frac{1}{2}m\dot{\delta}(\tau)^{2}+V_{HO}\left(\delta(\tau)\right)\right)+\frac{1}{4}\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}K\left(\tau-\tau^{\prime}\right)(\delta(\tau)-\delta(\tau^{\prime}))^{2}$
(56)

where $S_{HO}[q_{\textbf{cl}}(\tau)]$ is the action for the corresponding
classical path $q_{\textbf{cl}}(\tau)$ for the HO, given by weiss2012quantum ,

$S_{HO}[q_{\textbf{cl}}(\tau)]=\frac{q^{2}}{2\braket{\hat{q}^{2}}}+2\braket{\hat{p}^{2}}\eta^{2}$
(57)

Then, in Eq. 3, we have,

$\displaystyle\rho_{HO}(q,\eta)=\tilde{Z}^{-1}e^{-S_{HO}[q_{\textbf{cl}}(\tau)]}\mathcal{I}$
(58)

where, analogous to Eq. 44, we have,

$\displaystyle\mathcal{I}=\int_{0}^{0}\mathcal{D}\delta(\tau)e^{-S_{HO}[\delta(\tau)]}$
(59)

Then, using the expression for $\rho_{HO}(q,\eta)$ from Eq. 47 and using Eq.
58 and Eq. 57, we can obtain an expression for $\tilde{Z}$ as,

$\tilde{Z}=\mathcal{I}\sqrt{2\pi\braket{\hat{X}^{2}}}$ (60)

###  2 The path integral for the quantum deviation about the classical path

In this section, we are going to derive the expression for $\mathcal{I}(q)$
[Eq. 44] and consequently $\mathcal{G}(q)$ [Eq. 46] and also derive the
criteria for the convergence of these quantities. Let us do a Fourier series
decomposition of $\delta(\tau)$ [Eq. 54] as,

$\displaystyle\delta(\tau)$
$\displaystyle\equiv\frac{1}{\beta}\sum_{n=-\infty}^{\infty}b_{n}e^{i\nu_{n}\tau}$
(61)

####  1 Action in terms of fourier components

The main result of this subsection, given in Eq. 74, is to express
$S_{HO}[\delta(\tau)]$ [Eq. 56] in terms of the Fourier components $b_{n}$ of
$\delta(\tau)$ [Eq. 61]. What follows is some routine but tedious math and one
can skip directly to Eq. 74.

Given Eq. 61, the kinetic part of the action [Eq. 56] can be evaluated using,

$\int_{0}^{\beta}d\tau\dot{\delta}(\tau)^{2}=\frac{1}{\beta^{2}}\int_{0}^{\beta}d\tau\left(\frac{d}{d\tau}\sum_{n=-\infty}^{\infty}b_{n}e^{i\nu_{n}\tau}\right)^{2}=\frac{1}{\beta}\sum_{n}\nu_{n}^{2}|b_{n}|^{2}$
(62)

Here, we have used the fact that $b_{-n}=b_{n}^{*}$. The potential part of Eq.
56 can be evaluated as,

$\displaystyle\int_{0}^{\beta}\delta(\tau)^{2}=\frac{1}{\beta^{2}}\left(\sum_{n=-\infty}^{\infty}b_{n}e^{i\nu_{n}\tau}\right)^{2}=\frac{1}{\beta}\sum_{n}|b_{n}|^{2}$
(63)

In order to evaluate the dissipative term, let us first expand the following
expression,

$\displaystyle(\delta(\tau)-\delta(\tau^{\prime}))^{2}$
$\displaystyle=\frac{1}{\beta^{2}}\left(\sum_{n}b_{n}e^{i\nu_{n}\tau}-\sum_{m}b_{m}e^{i\nu_{m}\tau^{\prime}}\right)^{2}$
(64)
$\displaystyle=\frac{1}{\beta^{2}}\sum_{p}\sum_{q}b_{p}b_{q}e^{i(\nu_{p}+\nu_{q})\tau}+\sum_{p}\sum_{q}b_{p}b_{q}e^{i(\nu_{p}+\nu_{q})\tau^{\prime}}-2\sum_{p}\sum_{q}b_{p}b_{q}e^{i(\nu_{p}\tau+\nu_{q}\tau^{\prime})}$

Furthermore, let us first Fourier expand the influence functional as,

$K(\tau-\tau^{\prime})=\frac{m}{\beta}\sum_{k}\xi_{k}e^{i\nu_{k}(\tau-\tau^{\prime})}$
(65)

where we have weiss2012quantum ,

$\displaystyle\xi_{n}$ $\displaystyle=\frac{1}{m}\int_{0}^{\infty}d\omega
J(\omega)\frac{2\omega}{\nu_{n}^{2}+\omega^{2}}$ (66)

The dissipative term can then be evaluated as,

$\displaystyle\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}K\left(\tau-\tau^{\prime}\right)(\delta(\tau)-\delta(\tau^{\prime}))^{2}=\frac{1}{\beta^{2}}\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}\frac{m}{\beta}\sum\xi_{k}e^{i\nu_{k}(\tau-\tau^{\prime})}\left(\sum
b_{n}e^{i\nu_{n}\tau}-\sum b_{m}e^{i\nu_{m}\tau^{\prime}}\right)^{2}$ (67)

$\displaystyle=\frac{m}{\beta^{3}}\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}\sum\xi_{k}e^{i\nu_{k}(\tau-\tau^{\prime})}\left(\sum_{p}\sum_{q}\left(b_{p}b_{q}e^{i(\nu_{p}+\nu_{q})\tau}+b_{p}b_{q}e^{i(\nu_{p}+\nu_{q})\tau^{\prime}}-2b_{p}b_{q}e^{i(\nu_{p}\tau+\nu_{q}\tau^{\prime})}\right)\right)$
(68)
$\displaystyle=\frac{m}{\beta^{3}}\left(\beta\int_{0}^{\beta}d\tau\sum_{p,q}2\xi_{0}b_{p}b_{q}e^{i(\nu_{p}+\nu_{q})\tau}-2\int_{0}^{\beta}d\tau\int_{0}^{\beta}d\tau^{\prime}\sum_{p,q,k}\xi_{k}b_{p}b_{q}e^{i(\nu_{p}+\nu_{k})\tau+i(\nu_{q}-\nu_{k})\tau^{\prime}}\right)$
(69)
$\displaystyle=\frac{m}{\beta}\left(\sum_{p}2\xi_{0}b_{p}b_{-p}-2\sum_{p,q,k}\delta(p+k)\delta(q-k)\xi_{k}b_{p}b_{q}\right)$
(70)
$\displaystyle=\frac{2m}{\beta}\sum_{p}(\xi_{0}-\xi_{p})|b_{p}|^{2}\quad\text{sum
over $p$ and $q$ then rename $k\to p$}$ (71)
$\displaystyle=\frac{2m}{\beta}\sum_{p}(\zeta_{p})|b_{p}|^{2}$ (72)

where we used $\zeta_{p}=\xi_{0}-\xi_{p}$ [Eq. 42], which can be derived from
Eq. 66 as,

$\displaystyle\xi_{0}-\xi_{n}$
$\displaystyle=\frac{1}{m}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega}\left(\frac{2\nu_{n}^{2}}{(\nu_{n}^{2}+\omega^{2})}\right)=\zeta_{n}\qquad\text{from
\autoref{eqn:HO_zeta}}$ (73)

Finally, using the simplified expressions for the kinetic [Eq. 62], potential
[Eq. 63] and dissipative [Eq. 72] part of the action Eq. 56, we have,

$\displaystyle S_{HO}[\delta(\tau)]$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{1}{2}\frac{m}{\beta}|b_{n}|^{2}\left(\nu_{n}^{2}+\omega^{2}+\zeta_{n}\right)$
(74)

####  2 Calculating the deviation integral

Let us define,

$\displaystyle\alpha_{n}$
$\displaystyle\equiv\frac{m}{\beta}(\omega^{2}+\nu_{n}^{2}+\zeta_{n})$ (75)

The factor $\frac{1}{\alpha_{n}}$ gives an effective cutoff for the value of
$|b_{n}^{2}|$ for each value of $n$. For any $n$, if $|b_{n}|^{2}$ crosses
this cutoff, then the influence of the corresponding path will decay off
rapidly. Let us move to the Cartesian coordinates, defining $x_{n}$ and
$y_{n}$ as the real and imaginary parts of the complex number $b_{n}$ (i.e.,
$x_{n}\equiv\Re(b_{n})$ and $y_{n}\equiv\Im(b_{n})$). Since $\delta(0)=0$, the
Fourier component of any function of this type must satisfy,

$\displaystyle\sum_{n=-\infty}^{\infty}x_{n}$ $\displaystyle=0$ (76)
$\displaystyle\implies x_{0}$ $\displaystyle=-2\sum_{n=1}^{\infty}x_{n}$ (77)

Hence, $x_{0}$ is not an independent parameter in the path integral.
Similarly, since $x_{-n}=x_{n}$ and $y_{-n}=-y_{n}$, hence the only
independent variables are those corresponding to $n\geq 1$. The path integral
$\mathcal{I}$ [Eq. 44] then becomes,

$\displaystyle\mathcal{I}$
$\displaystyle=J\left(\prod_{n=1}^{\infty}\int_{-\infty}^{\infty}dy_{n}\int_{-\infty}^{\infty}dx_{n}\right)\exp\left\\{-\frac{1}{2}\alpha_{0}\left(2\sum_{n=1}^{\infty}x_{n}\right)^{2}-2\sum_{n=1}^{\infty}\frac{1}{2}\alpha_{n}\left(x_{n}^{2}+y_{n}^{2}\right)\right\\}$
(78)

Here, $J$ is the Jacobian coming from the linear transformation given in Eq.
61. We will now analytically evaluate this integral. But we should first note
that this integral will diverge if for $n\geq 1$, $\alpha_{n}\leq 0$, which
would break down LHA.

Before proceeding, note that we have
$\mathcal{\alpha}_{|n+1|}>\mathcal{\alpha}_{|n|}$. This is true because
$\nu_{|n+1|}>\nu_{|n|}$ follows from definition and immediately gives
$\zeta_{|n+1|}>\zeta_{|n|}$ since $x/(a+x)$ monotonically increases for
positive $x$ and $a$. Hence, as a necessary condition for Eq. 78 to be
convergent, we just need to impose that $\alpha_{1}>0$. We will hence restrict
ourselves to this constraint. (As a side remark, note that this assumption
automatically ensures the convergence of $\braket{\hat{P}^{2}(q)}$ [Eq. 39].)
But it is still unclear whether $\alpha_{0}\leq 0$ will lead to divergence in
Eq. 78 or not. In this section, we will see that $\mathcal{I}$ can be
convergent even if $\alpha_{0}\leq 0$.

After doing the $y_{n}$ integral in Eq. 78, we get,

$\displaystyle\mathcal{I}=J\left(\prod_{m=1}^{\infty}\sqrt{\frac{\pi}{\alpha_{m}}}\right)\left(\prod_{n=1}\int_{-\infty}^{\infty}dx_{n}\right)\exp\left\\{-\left(2\alpha_{0}\left[\sum_{n=1}^{\infty}x_{n}\right]^{2}+\sum_{n=1}^{\infty}\alpha_{n}x_{n}^{2}\right)\right\\}$
(79)

Let us define the negative of the exponent above as,

$\displaystyle\mathcal{F}\equiv\left(2\alpha_{0}\left[\sum_{n=1}^{\infty}x_{n}\right]^{2}+\sum_{n=1}^{\infty}\alpha_{n}x_{n}^{2}\right)$
(80)

If $\min_{x_{n}}\mathcal{F}>0$, then $\mathcal{I}$ will converge. We can think
of $\mathcal{F}$ as a diagonal component of a matrix $M$, i.e.
$\mathcal{F}\equiv\langle v|M|v\rangle$ where, $\langle v|=\langle
x_{1},x_{2},x_{3},...|$ and $M=D+2\alpha_{0}J$ where, for $i,j\geq 1$,
$J_{ij}=1$ and $D_{ij}=\delta_{ij}\alpha_{i}$. Let $\lambda_{n}$ be the
eigenvalues of matrix $M$. Then we have,

$\displaystyle\mathcal{I}=J\left(\prod_{n=1}^{\infty}\frac{\pi}{\sqrt{\alpha_{n}\lambda_{n}}}\right)$
(81)

Before proceeding, let us first simplify an expression we will need
immediately. From the definition of $\braket{\hat{q}^{2}}$ [Eq. 38] and
$\alpha_{n}$ [Eq. 75], we have,

$\displaystyle\braket{\hat{q}^{2}}$
$\displaystyle=\frac{1}{\beta^{2}}\sum_{n=-\infty}^{\infty}\frac{1}{\alpha_{n}}$
(82) $\displaystyle\implies 2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}$
$\displaystyle=\beta^{2}\braket{\hat{q}^{2}}-\frac{1}{\alpha_{0}}$ (83)

Proceeding with the derivation of $\mathcal{I}$ again, we next need to
calculate the value of the product $\prod_{n=1}^{\infty}\lambda_{n}$. We note
that this is just given by the determinant of the matrix $M$, i.e.,

$\displaystyle\prod_{n=1}^{\infty}\lambda_{n}$ $\displaystyle=|M|$ (84)
$\displaystyle=\left(\prod_{n=1}^{\infty}\alpha_{n}\right)+2\alpha_{0}\left(\prod_{n=1}^{\infty}\alpha_{n}\right)\left(\sum_{m=1}^{\infty}\frac{1}{\alpha_{m}}\right)$
(85)
$\displaystyle=\left(\prod_{n=1}^{\infty}\alpha_{n}\right)\left(1+\alpha_{0}\left(\beta^{2}\braket{\hat{q}^{2}}-\frac{1}{\alpha_{0}}\right)\right)\qquad\text{from
\autoref{eqn:sum_as_q2}}$ (86)
$\displaystyle=\alpha_{0}\beta^{2}\braket{\hat{q}^{2}}\prod_{n=1}^{\infty}\alpha_{n}$
(87)

We hence have,

$\displaystyle\mathcal{I}$
$\displaystyle=J\left(\frac{1}{\sqrt{\alpha_{0}\beta^{2}\braket{\hat{q}^{2}}}}\prod_{n=1}^{\infty}\frac{\pi}{\alpha_{n}}\right)$
(88)

This infinite product will give an infinitesimally small result because
$\alpha_{n}$ scales as $n^{2}$. But to get a reasonable result for LHA, we can
proceed as follows. In the expression for LHA [Eq. 36], we effectively have a
factor of, $\frac{\mathcal{I}(q)}{Z_{R}}$. Let us transform it as,

$\displaystyle\frac{\mathcal{I}(q)}{Z_{R}}\equiv\frac{\mathcal{G}(q)}{Z}$ (89)

where we define $Z$ as,

$\displaystyle
Z=Z_{R}\left(\frac{1}{J}\prod_{n=1}^{\infty}\frac{\frac{m}{\beta}(\nu_{n}^{2}+\zeta_{n})}{\pi}\right)$
(90)

We hence have, for $\mathcal{G}(q)$,

$\displaystyle\mathcal{G}$
$\displaystyle=\frac{Z}{Z_{R}}\mathcal{I}=\left(\frac{1}{J}\prod_{n=1}^{\infty}\frac{\frac{m}{\beta}(\nu_{n}^{2}+\zeta_{n})}{\pi}\right)\mathcal{I}$
(91)
$\displaystyle=\frac{1}{\sqrt{\alpha_{0}\beta^{2}\braket{\hat{q}^{2}}}}\prod_{n=1}^{\infty}\frac{\alpha_{n}-\frac{m\omega^{2}}{\beta}}{\alpha_{n}}$
(92)

This quantity can be shown to converge. On the other hand, $Z$ serves as an
overall normalization factor for LHA since, as we remarked in Appendix. A, it
depends only upon $k_{B}T$, $\Omega$ and $\gamma$, but not on $\omega$, and is
hence independent of $q$. Note that the Jacobian $J$, coming from the linear
transformation given in Eq. 61, is also $\omega$ independent.

Now, we are going to investigate the condition for Eq. 92 to converge. Since
we have assumed $\alpha_{n\geq 1}>0$, we note that for Eq. 92 to be real, we
need,

$\displaystyle\alpha_{0}\braket{\hat{q}^{2}}$ $\displaystyle>0$ (93)
$\displaystyle\implies\alpha_{0}\left(\frac{1}{\alpha_{0}}+2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}\right)$
$\displaystyle>0\qquad\text{from \autoref{eqn:sum_as_q2}}$ (94)

If $\alpha_{0}$ is also positive, then Eq. 94 is automatically satisfied. If
$\alpha_{0}<0$, then we are interested in the highest value of $\alpha_{0}$
such that we have,

$\displaystyle\left(\frac{1}{\alpha_{0}}+2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}\right)$
$\displaystyle=0$ (95)
$\displaystyle\implies\sum_{n=1}^{\infty}\frac{\alpha_{0}}{\alpha_{n}}$
$\displaystyle=-\frac{1}{2}$ (96)
$\displaystyle\implies\sum_{n=1}^{\infty}\frac{m\omega^{2}}{m\omega^{2}+m\nu_{n}+m\zeta_{n}}$
$\displaystyle=-\frac{1}{2}$ (97)

Here, we used the definition of $\alpha_{n}$ [Eq. 75]. Redefining
$m\omega^{2}=-x$ and $m\nu_{n}^{2}+m\zeta_{n}=y_{n}$, we get Eq. 13, where now
we are interested in the smallest solution for $x$ in the following equation,

$\displaystyle\implies\sum_{n=1}^{\infty}\frac{x}{y_{n}-x}$
$\displaystyle=\frac{1}{2}$ (98)

The smallest solution for $x$ gives us the value of $V_{P}$, defined in Sec.
II, such that LHA will give convergent results for $V_{2}(q)>-V_{P}$.

Note that, from Eq. 93, for $\alpha_{0}<0$, we are equivalently interested in
the highest negative value of $\omega^{2}$ such that $\braket{\hat{q}^{2}}=0$.
Hence, for the value of $\omega^{2}$ higher than this, convergence of
$\braket{\hat{X}^{2}(q)}^{-1}$ [Eq. 38] is also guaranteed, as required by LHA
[Eq. 37] to be well defined. Also, since some spectral densities admit closed
form expression for $\braket{\hat{q}^{2}}$, this condition can easily be
checked for in these cases.

In summary, we can say that LHA gives convergent results if,

  * •

$\alpha_{0}>0$,

  * •

$\alpha_{0}<0$ and $\alpha_{1}>0$, and Eq. 94 is true. The value of $V_{P}$ is
given as the smallest solution for $x$ in Eq. 98.

##  Appendix C Error estimation

In this section, we are going to derive the expression for the estimation of
the error associated with LHA [Eq. 16]. We start with Eq. 3. Parameterising
the path variable as $q(\tau)=q+\delta(\tau)$, where $q(\tau)$ starts and ends
at points $q+\eta$ and $q-\eta$, respectively, we can write,

$V(q+\delta)=V(q)+V_{1}(q)\delta+\frac{1}{2!}V_{2}(q)\delta^{2}+\frac{1}{3!}V_{3}(q)\delta^{3}+O(\delta^{4})$
(99)

Let us write the full action in terms of the quadratic part and higher order
terms in the potential, as,

$\displaystyle
S[q+\delta(\tau)]=S_{q}[q+\delta(\tau)]+\int_{0}^{\beta}d\tau\sum_{n=3}\frac{1}{n!}V_{n}(q)\delta(\tau)^{n}$
(100)

Then, writing the exact state $\sigma(q,\eta)$ in terms of harmonic action
$S_{q}[q+\delta(\tau)]$ and higher order corrections over it, we get,

$\displaystyle\sigma(q,\eta)=Z^{\prime-1}\left\\{\int_{-\eta}^{\eta}\mathcal{D}\delta(\tau)\exp\left(-\int_{0}^{\beta}d\tau\sum_{n=3}\frac{1}{n!}V_{n}(q)\delta(\tau)^{n}\right)e^{-S_{q}[q+\delta(\tau)]}\right\\}$
(101)

Here, $Z^{\prime}$ is the normalization factor defined as,

$\displaystyle
Z^{\prime}=\int_{-\infty}^{\infty}dq\left\\{\int_{0}^{0}\mathcal{D}\delta(\tau)\exp\left(-\int_{0}^{\beta}d\tau\sum_{n=3}\frac{1}{n!}V_{n}(q)\delta(\tau)^{n}\right)e^{-S_{q}[q+\delta(\tau)]}\right\\}$
(102)

Let us approximate Eq. 101 as,

$\displaystyle\sigma(q,\eta)$
$\displaystyle\approx\rho(q,\eta)-\sum_{n=3}\frac{V_{n}(q)}{n!}\left\\{\int_{-\eta}^{\eta}\mathcal{D}\delta(\tau)\left(\int_{0}^{\beta}d\tau\delta(\tau)^{n}\right)e^{-S_{q}[q+\delta(\tau)]}\right\\}$
(103)
$\displaystyle\approx\rho(q,\eta)-\frac{V_{3}(q)}{3!}\left\\{\int_{-\eta}^{\eta}\mathcal{D}\delta(\tau)\left(\int_{0}^{\beta}d\tau\delta(\tau)^{3}\right)e^{-S_{q}[q+\delta(\tau)]}\right\\}$
(104)

Here, $\rho(q,\eta)$ is the LHA state and in the last step, we have made a
first order approximation. We now need to calculate the quantity,

$\displaystyle
I_{3}=\int_{-\eta}^{\eta}\mathcal{D}\delta(\tau)\left(\int_{0}^{\beta}d\tau\delta(\tau)^{3}\right)e^{-S_{q}[q+\delta(\tau)]}$
(105)

Here, we are going to estimate this path integral using the following
argument. Reparametrize $q+\delta(\tau)\to
q_{\textbf{cl}}(\tau)+\delta(\tau)$. Then using $S_{q}[q+\delta(\tau)]\to
S_{q}[q_{\textbf{cl}}(\tau)+\delta(\tau)]=S_{q}[q_{\textbf{cl}}(\tau)]+S_{q}[\delta(\tau)]$,
we get,

$\displaystyle I_{3}$
$\displaystyle=e^{-S_{q}[q_{\textbf{cl}}(\tau)]}\int_{0}^{0}\mathcal{D}\delta(\tau)\left(\int_{0}^{\beta}d\tau\left(q_{\textbf{cl}}(\tau)-q+\delta(\tau)\right)^{3}\right)e^{-S_{q}[\delta(\tau)]}$
(106)

Now, if $n$ is an odd number, then we have
$(-\delta(\tau))^{n}=-\delta(\tau)^{n}$. But since
$S_{q}[-\delta(\tau)]=S_{q}[\delta(\tau)]$, when we take the full path
integral, we get,

$\displaystyle\int_{0}^{0}\mathcal{D}\delta(\tau)\left(\int_{0}^{\beta}d\tau(q_{\textbf{cl}}(\tau)-q)^{m}\delta(\tau)^{n}\right)e^{-S_{q}[\delta(\tau)]}$
$\displaystyle=0\qquad\text{for odd}\;n$ (107) $\displaystyle\implies I_{3}$
$\displaystyle=e^{-S_{q}[q_{\textbf{cl}}(\tau)]}\int_{0}^{0}\mathcal{D}\delta(\tau)\left(\int_{0}^{\beta}d\tau\left((q_{\textbf{cl}}(\tau)-q)^{3}+3(q_{\textbf{cl}}(\tau)-q)\delta(\tau)^{2}\right)\right)e^{-S_{q}[\delta(\tau)]}$
(108)

Now, we are going to replace $q_{\textbf{cl}}(\tau)-q$ with its maximum value,
given by $\Delta(q,\eta)$, and $\delta(\tau)$ with its typical value, given by
$\mathcal{E}(q)$, to obtain,

$\displaystyle
I_{3}\approx\beta\left(\Delta(q,\eta)^{3}+3\mathcal{E}(q)^{2}\Delta(q,\eta)\right)\int_{-\eta}^{\eta}\mathcal{D}\delta(\tau)e^{-S_{q}[q+\delta(\tau)]}$
(109)

We will estimate $\Delta(q,\eta)$ and $\mathcal{E}(q)$ in Sec. 1.

Let’s define a quantity, $\mathcal{T}(q,\eta)$, as a measure of the relative
error corresponding to the matrix element $\rho(q,\eta)$ as,

$\displaystyle\mathcal{T}(q,\eta)$
$\displaystyle\equiv-\beta\frac{V_{3}}{3!}\left(\Delta(q,\eta)^{3}+3\mathcal{E}(q)^{2}\Delta(q,\eta)\right)$
(110)

Also, for convenience, let us define,

$\displaystyle\epsilon_{T}$ $\displaystyle\equiv\int
dq\rho(q)|\mathcal{T}(q)|$ (111) $\displaystyle=\braket{}{\mathcal{T}(q)}{}$
(112)

Then, Eq. 101 becomes,

$\displaystyle\sigma(q,\eta)$
$\displaystyle\approx\frac{\rho(q,\eta)\pm\rho(q,\eta)|\mathcal{T}(q,\eta)|}{Z^{\prime}}$
(113)
$\displaystyle=\frac{\rho(q,\eta)\pm\rho(q,\eta)|\mathcal{T}(q,\eta)|}{1\pm\epsilon_{T}}$
(114) $\displaystyle\approx\rho(q,\eta)\pm\epsilon(q,\eta)\rho(q,\eta)$ (115)

where,

$\displaystyle\epsilon(q,\eta)=|\mathcal{T}(q,\eta)|+\epsilon_{T}$ (116)

gives a measure of the total relative error for the matrix element
$\rho(q,\eta)$. To get an expression for the relative error associated with an
observable $\hat{O}$, we define,

$\displaystyle\epsilon_{\hat{O}}=\frac{1}{\braket{\hat{O}}}\int dq\int
d\eta\;\hat{O}(q,\eta)\epsilon(q,\eta)|\rho(q,\eta)|$ (117)

###  1 Length Scale estimation

In the following subsections, we will estimate $\Delta(q,\eta)$ and
$\mathcal{E}(q)$. Let $q_{\textbf{cl}}(\tau)$ be the classical path for an HO,
starting and ending at points $q+\eta$ and $q-\eta$, respectively, given by,

$q_{\textbf{cl}}(\tau)=\frac{1}{\beta}\sum_{n=-\infty}^{\infty}a_{n}e^{i\nu_{n}\tau}$
(118)

Then, $a_{n}$ is given as weiss2012quantum ,

$\displaystyle a_{n}$
$\displaystyle=\left\\{\left(\frac{i\nu_{n}}{\omega^{2}+\nu_{n}^{2}+\zeta_{n}}\right)2\eta+\frac{1}{\omega^{2}+\nu_{n}^{2}+\zeta_{n}}\frac{q}{m\braket{\hat{q}^{2}}}\right\\}$
(119)

We conjecture that for $\eta=0$, $|q_{\textbf{cl}}(\tau)-q|$ will have its
global maxima at $\tau=\beta/2$. This intuition rests on the fact that in the
general path integral equation [Eq. 3], although the potential term will push
the extremum classical path towards the potential minima, the kinetic and
dissipative terms will help remove any oscillations from it. Hence, the curve
will have a single extrema, and since for $\eta=0$, the classical path can be
seen to be made up of only cosine Fourier components and hence is symmetric
about $\tau=\beta/2$, therefore this extrema will be at $\tau=\beta/2$. We
have also verified this numerically. When $\eta\neq 0$, we note that although
the classical path still does not have any oscillations, $\tau=\beta/2$ does
not give the extremum point anymore, but we assume that the maximum among
$|q_{\textbf{cl}}(\beta/2)-q|$ and $|\eta|$ still gives a rough estimate for
the maximum deviation of the classical path from the constant path
$q^{\prime}(\tau)=q$.

We will hence use the quantity $q_{\textbf{cl}}(\beta/2)-q$ to estimate
$\Delta(q,\eta)$. We simplify this expression as,

$\displaystyle q_{\textbf{cl}}(\beta/2)-q$
$\displaystyle=-q+\frac{1}{\beta}\sum_{n=-\infty}^{\infty}(-1)^{n}a_{n}$ (120)
$\displaystyle=q\left(\frac{1}{m\beta\braket{\hat{q}^{2}}}\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{\nu_{n}^{2}+\omega^{2}+\zeta_{n}}-1\right)$
(121)

Hence, the estimate for the maximum deviation of the classical path is then
given by,

$\displaystyle\Delta_{HO}(q,\eta)=\mathop{\mathrm{max}}\left\\{\left|q\left(-1+\frac{1}{m\beta\braket{\hat{q}^{2}}}\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{\nu_{n}^{2}+\omega^{2}+\zeta_{n}}\right)\right|,|\eta|\right\\}$
(122)

Relating $q\equiv V_{1}(q)/V_{2}(q)$ [Eq. 11], we directly recover the
expression for the maximum deviation of the classical path claimed earlier
[Eq. 50].

Next, in order to estimate $\mathcal{E}(q)$, rename $r_{n}\equiv|b_{n}|$ [Eq.
74] and write the path integral for $\mathcal{I}$ [Eq. 78] in polar coordinate
as,

$\displaystyle\mathcal{I}=J\left(\prod_{n=1}\int_{0}^{\infty}dr_{n}\int_{0}^{2\pi}r_{n}d\theta_{n}\right)e^{-\left(\frac{1}{2}\alpha_{0}b_{0}^{2}+2\sum_{n=1}^{\infty}\frac{1}{2}\alpha_{n}r_{n}^{2}\right)}$
(123)

Note that $\theta_{n}$ dependence in the integrand is in $b_{0}$ and comes
through the following constraint [Eq. 76],

$\displaystyle b_{0}=-2\sum_{n=1}^{\infty}\Re(b_{n})$ (124)

We are now in the position to immediately do all the $d\theta_{n}$ integrals
in the Eq. 123. For this, we note that for a fixed value of $b_{n}$’s, where
$n=0,1,2,...$, the $d\theta_{n}$ integrals in Eq. 123 give us a measure of all
the configurations that satisfy the constraint in Eq. 124.

Now consider a 2D random walk starting from origin, where the nth step size
for the walk is given by the function $\mathcal{S}(n)=2|b_{n}|=2r_{n}$. Then
after $N\to\infty$ steps, the probability distribution $P(x)$ for this walker
to have the x-axis position value of $-b_{0}$ gives a measure that is
equivalent to doing all the $d\theta_{n}$ integrals in Eq. 123. Note that
there is a factor of 2 in the definition for $\mathcal{S}(n)$ because of the
corresponding factor of 2 in Eq. 124.

This probability distribution is given as sherman1972summation ,

$P(x)=\frac{1}{\sigma}\sqrt{\frac{1}{2\pi}}e^{-\frac{x^{2}}{2\sigma^{2}}}$
(125)

where, the variance $\sigma^{2}$ is given by,

$\displaystyle\sigma^{2}$
$\displaystyle=\frac{1}{2}\sum_{n=1}^{\infty}(\mathcal{S}(n))^{2}$ (126)
$\displaystyle=2\sum_{n=1}^{\infty}r_{n}^{2}$ (127)

Note that the Gaussian probability distribution for the random walk is true
only if a few sets of steps, $\mathcal{S}(n)$, do not dominate over the others
in size. If this is not the case, the actual distribution is more complicated
(although it is still a decaying function with some width), and is given in
terms of Bessel functions sherman1972summation . Here, for simplicity, we
assume Gaussian probability distribution.

We hence replace all the $d_{\theta_{n}}$ integrals with a single integral
over the variable $b_{0}$, and Eq. 123 becomes,

$\displaystyle\mathcal{I}$
$\displaystyle=J\left(\prod_{n=1}^{\infty}\int_{0}^{\infty}r_{n}dr_{n}\right)\int_{-R}^{R}db_{0}\frac{1}{\sigma}\sqrt{\frac{1}{2\pi}}e^{-\left(\frac{1}{2}\left(\alpha_{0}+\frac{1}{\sigma^{2}}\right)b_{0}^{2}+2\sum_{n=1}^{\infty}\frac{1}{2}\alpha_{n}r_{n}^{2}\right)}$
(128)

where $R=2\sum_{n=1}^{\infty}r_{n}$ is put here to ensure that Eq. 124 is not
violated.

Note that this equation gives a typical length-scale for the value for $r_{n}$
and $b_{0}$ as following,

$\displaystyle r_{n}^{2}=|b_{n}|^{2}$
$\displaystyle\approx\frac{1}{\alpha_{n}}$ (129)
$\displaystyle\implies\sigma^{2}$ $\displaystyle\approx
2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}\qquad\text{from \autoref{eqn:sigma
squared}}$ (130) $\displaystyle\implies b_{0}^{2}$
$\displaystyle\approx\frac{1}{\alpha_{0}+\frac{1}{2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}}}$
(131)

Using the typical value of $|b_{n}|$ from here, we can go back in the position
space representation of $\delta(\tau)$ [Eq. 61] that corresponds to the paths
that have important contribution to $\mathcal{I}$, as,

$\displaystyle\delta(\tau)=\frac{1}{\beta}\left(\frac{1}{\sqrt{\alpha_{0}+\frac{1}{2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}}}}+\sum_{n\neq
0}\frac{1}{\sqrt{\alpha_{n}}}e^{i\theta_{n}}e^{i\nu_{n}\tau}\right)$ (132)

where, $\theta_{-n}=-\theta_{n}$ are random phases with the condition, coming
from Eq. 124, that $\delta(\tau=0)=0$.

We are interested in the typical deviation of $\delta(\tau)$. To simplify the
problem, we now relax the condition $\delta(\tau=0)=0$ so that $\theta_{n}$
become truly random and this again reduces to the random walk problem with
fixed step size sherman1972summation . To be precise, the mathematical problem
we are now dealing with is- given a Fourier series with fixed coefficients but
random phases $\phi_{n}$,

$\displaystyle h(x)=c_{0}+\sum_{n=1}^{\infty}c_{n}\sin{\left(2\pi
nx+\phi_{n}\right)}$ (133)

What is the probability density, $P(y)$, that $h(x)=y$ for any given value of
$x$?

Using similar random walk arguments as used before, $P(y)$ can be shown to be
approximately given by,

$\displaystyle
P(y)=\frac{1}{\sigma}\sqrt{\frac{1}{2\pi}}e^{\frac{-(y-c_{0})^{2}}{2\sigma^{2}}}$
(134)

where $\sigma$ is given by,

$\displaystyle\sigma$
$\displaystyle=\sqrt{\frac{1}{2}\sum_{n=1}^{\infty}c_{n}^{2}}$ (135)

An estimate for $\mathcal{E}_{HO}(q)$ can now be given by,

$\displaystyle\mathcal{E}_{HO}(q)=c_{0}+\sigma$ (136)

For our case, from Eq. 132, we have,

$\displaystyle c_{0}$
$\displaystyle=\frac{1}{\beta}\frac{1}{\sqrt{\alpha_{0}+\frac{1}{2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}}}}$
(137) $\displaystyle c_{n\geq 1}$
$\displaystyle=\frac{1}{\beta}\frac{2}{\sqrt{\alpha_{n}}}$ (138)

Let us first simplify the following expression,

$\displaystyle\frac{1}{\alpha_{0}+\frac{1}{2\sum_{n=1}^{\infty}\frac{1}{\alpha_{n}}}}$
$\displaystyle=\frac{1}{\alpha_{0}+\frac{1}{\beta^{2}\braket{\hat{q}^{2}}-\frac{1}{\alpha_{0}}}}\qquad\text{from
\autoref{eqn:sum_as_q2}}$ (139)
$\displaystyle=\frac{\beta^{2}\braket{\hat{q}^{2}}-\frac{1}{\alpha_{0}}}{\alpha_{0}\beta^{2}\braket{\hat{q}^{2}}}$
(140)

Hence, from Eq. 135, Eq. 136, Eq. 137 and Eq. 138, we have,

$\displaystyle\mathcal{E}_{HO}(q)$
$\displaystyle=\frac{1}{\beta}\left(\sqrt{\frac{\beta^{2}\braket{\hat{q}^{2}}-\frac{1}{\alpha_{0}}}{\alpha_{0}\beta^{2}\braket{\hat{q}^{2}}}}+\sqrt{\beta^{2}\braket{\hat{q}^{2}}-\frac{1}{\alpha_{0}}}\right)\qquad\text{from
\autoref{eqn:sum_as_q2}}$ (141)
$\displaystyle=\sqrt{\braket{\hat{q}^{2}}-\frac{1}{\beta^{2}\alpha_{0}}}\left(\sqrt{\frac{1}{\alpha_{0}\beta^{2}\braket{\hat{q}^{2}}}}+1\right)$
(142)
$\displaystyle=\sqrt{\braket{\hat{q}^{2}}-\frac{1}{m\omega^{2}\beta}}\left(\sqrt{\frac{1}{m\omega^{2}\beta\braket{\hat{q}^{2}}}}+1\right)$
(143)

where we have used $\alpha_{0}=\frac{m\omega^{2}}{\beta}$ from Eq. 75. This
recovers the expression for $\mathcal{E}(q)$ [Eq. 51], the typical quantum
deviation of all the paths about $q_{cl}(\tau)$ that have significant
contribution to the action.

##  Appendix D Special cases

###  1 Ultra strong coupling limit

In the ultra strong coupling limit, i.e., $\lim_{J(\omega)\to\infty}$, first
note that, from the definition of $\zeta_{n}$ [Eq. 42], we have,

$\displaystyle\lim_{J(\omega)\to\infty}\zeta_{n\neq 0}=\infty$ (144)

From the definition of $\braket{\hat{X}^{2}(q)}$ [Eq. 38], we have,

$\displaystyle\lim_{J(\omega)\to\infty}\braket{\hat{X}^{2}(q)}$
$\displaystyle=\lim_{J(\omega)\to\infty}\frac{1}{m\beta}\sum_{n=-\infty}^{+\infty}\frac{1}{\frac{V_{2}(q)}{m}+\nu_{n}^{2}+\zeta_{n}}$
(145) $\displaystyle=\frac{1}{\beta V_{2}(q)}$ (146)

Similarly, from definition of $\braket{\hat{P}^{2}(q)}$ [Eq. 39], we have,

$\displaystyle\lim_{J(\omega)\to\infty}\braket{\hat{P}^{2}(q)}$
$\displaystyle=\frac{m}{\beta}+2\sum_{n=1}^{+\infty}\lim_{\zeta_{n}\to\infty}\frac{\frac{V_{2}(q)}{m}+\zeta_{n}}{\frac{V_{2}(q)}{m}+\nu_{n}^{2}+\zeta_{n}}$
(147) $\displaystyle=\frac{m}{\beta}+2\sum_{n=1}^{+\infty}1$ (148)
$\displaystyle=\infty$ (149)

Finally, from definition of $\mathcal{G}(q)$ [Eq. 46], we have,

$\displaystyle\lim_{J(\omega)\to\infty}\mathcal{G}(q)$
$\displaystyle=\lim_{J(\omega)\to\infty}\frac{1}{\sqrt{V_{2}(q)\beta\braket{\hat{X}^{2}(q)}}}\prod_{n=1}^{\infty}\frac{\nu_{n}^{2}+\zeta_{n}}{\mathcal{A}_{n}(q)}$
(150)
$\displaystyle=\lim_{J(\omega)\to\infty}\frac{1}{\sqrt{V_{2}(q)\beta\frac{1}{\beta
V_{2}(q)}}}\prod_{n=1}^{\infty}1$ (151) $\displaystyle=1$ (152)

Hence expression for USC-LHA [Eq. 37] becomes,

$\displaystyle\lim_{J(\omega)\to\infty}\rho(q,\eta)$
$\displaystyle=\delta(\eta)Z^{-1}\exp\left(-\frac{(V_{1}(q)/V_{2}(q))^{2}}{2\braket{\hat{X}^{2}(q)}}-\beta\left(V(q)-\frac{V_{1}(q)^{2}}{2V_{2}(q)}\right)\right)$
(153) $\displaystyle\implies\lim_{J(\omega)\to\infty}\rho(q)$
$\displaystyle=Z^{-1}\exp\left(-\beta V(q)\right)\exp\left(-\frac{\beta
V_{1}(q)^{2}}{2V_{2}(q)}+\frac{\beta V_{1}(q)^{2}}{2V_{2}(q)}\right)$ (154)
$\displaystyle=Z^{-1}\exp\left(-\beta V(q)\right)$ (155)

which is the CV version of the USC result cresser2021weak . To see this
equivalence, note that the USC result is usually stated as,

$\rho=Z_{USC}^{-1}\exp{\left(-\beta\sum_{n}P_{n}H_{s}P_{n}\right)}$ (156)

where $P_{n}=\ket{x_{n}}\bra{x_{n}}$ are projection operators on the non-
degenerate eigenstates $\ket{x_{n}}$ of the system coupling operators $X$ and
$Z_{USC}$ is a normalization factor.

We note that in the CV limit,

$\displaystyle\sum_{n}P_{n}H_{s}P_{n}$ $\displaystyle=\int
dq\left(\frac{\bra{q}\hat{P}^{2}\ket{q}}{2m}+V(q)\right)\ket{q}\bra{q}$ (157)

Now, $\sum_{q}\bra{q}\hat{P}^{2}\ket{q}\ket{q}\bra{q}$ is proportional to the
identity operator, and hence will commute with the potential operator.
Therefore, in the expression for the USC-MFGS, the contribution of the
momentum operator can be absorbed in the trace and we recover Eq. 155 as the
USC result.

Now, let us calculate $\Delta(q,\eta)$ in this limit. Let us first evaluate
$W(q)$ [Eq. 52] that comes up in the definition of $\Delta(q,\eta)$ as,

$\displaystyle\lim_{J(\omega)\to\infty}W(q)$ $\displaystyle=\frac{1}{\beta
V_{2}(q)}+\frac{1}{m\beta}\sum_{n=1}^{+\infty}\lim_{\zeta_{n}\to\infty}\frac{(-1)^{n}}{\frac{V_{2}(q)}{m}+\nu_{n}^{2}+\zeta_{n}}$
(158) $\displaystyle=\frac{1}{\beta V_{2}(q)}$ (159)

Hence, from the definition of $\Delta(q,\eta)$ [Eq. 50], we have,

$\displaystyle\lim_{J(\omega)\to\infty}\Delta(q,\eta)$
$\displaystyle=\max\left\\{\left|\lim_{J(\omega)\to\infty}\frac{V_{1}(q_{0})}{V_{2}(q_{0})}\left(1-\frac{W(q_{0})}{\braket{\hat{X}^{2}(q_{0})}}\right)\right|,\left|\eta\right|\right\\}$
(160)
$\displaystyle=\max\left\\{\left|\frac{V_{1}(q_{0})}{V_{2}(q_{0})}\left(1-\frac{\frac{1}{\beta
V_{2}(q)}}{\frac{1}{\beta V_{2}(q)}}\right)\right|,\left|\eta\right|\right\\}$
(161) $\displaystyle=\max\left\\{0,\left|\eta\right|\right\\}$ (162)

Similarly, let us calculate $\mathcal{E}(q)$ in this limit. From the
definition of $\mathcal{E}(q)$ [Eq. 51], we have,

$\displaystyle\lim_{J(\omega)\to\infty}\mathcal{E}(q)$
$\displaystyle=\lim_{J(\omega)\to\infty}\sqrt{\braket{\hat{X}^{2}(q)}-\frac{1}{V_{2}(q)\beta}}\left(\sqrt{\frac{1}{V_{2}(q)\beta\braket{\hat{X}^{2}(q)}}}+1\right)$
(164)
$\displaystyle=\sqrt{\frac{1}{V_{2}(q)\beta}-\frac{1}{V_{2}(q)\beta}}\left(\sqrt{\frac{1}{V_{2}(q)\frac{1}{V_{2}(q)\beta}}}+1\right)$
(165) $\displaystyle=0$ (166)

Hence, we note that both $\Delta(q,\eta)$ and $\mathcal{E}(q)$ have expected
behavior in the USC limit.

###  2 High temperature limit

To understand why LHA becomes accurate in the high temperature limit, consider
an arbitrary path $q(\tau)$ starting and ending at the points $q+\eta$ and
$q-\eta$, respectively. Now, for $a>1$, let us rescale $\beta\to\beta/a$ such
that we have $q(\tau)\to q_{a}(\tau)\equiv q(a\tau)$, where we have identified
the paths $q(\tau)$ and $q_{a}(\tau)$ with each other through a one-to-one and
onto map. Note that we have $q(0)=q+\eta=q_{a}(0)$ and
$q(\beta)=q-\eta=q_{a}(\beta/a)$. Now, define the new action,
$S_{a}[q_{a}(\tau)]$ as,

$\displaystyle S_{a}[q_{a}(\tau)]$
$\displaystyle=\int_{0}^{\beta/a}d\tau\left(\frac{1}{2}m\dot{q}_{a}(\tau)^{2}+V\left(q_{a}(\tau)\right)\right)+\frac{1}{4}\int_{0}^{\beta/a}d\tau\int_{0}^{\beta/a}d\tau^{\prime}K_{a}\left(\tau-\tau^{\prime}\right)\left(q_{a}(\tau)-q_{a}(\tau^{\prime})\right)^{2}$
(167)

where,

$K_{a}(\tau)=\int_{0}^{\infty}d\omega
J(\omega)\frac{\cosh{[\omega(\beta/2a-\tau)]}}{\sinh{[\omega\beta/2a]}}$ (168)

Then, we have,

$\displaystyle\int_{0}^{\beta/a}\frac{d^{2}}{d\tau^{2}}q_{a}(\tau)d\tau$
$\displaystyle=a^{2}\int_{0}^{\beta/a}\frac{d^{2}}{d(a\tau)^{2}}q(a\tau)d\tau$
(169) $\displaystyle=a\int_{0}^{\beta}\frac{d^{2}}{d\tau^{2}}q(\tau)d\tau$
(170)

Similarly, we have,

$\displaystyle\int_{0}^{\beta/a}V(q_{a}(\tau))d\tau=\frac{1}{a}\int_{0}^{\beta}V(q(\tau))d\tau$
(171)

Before moving to the dissipative term, we write the hyperbolic factor coming
in the correlation function [Eq. 168] in large $a$ limit for
$0\leq\tau\leq\beta/a$ as,

$\displaystyle\lim_{a>>1}\frac{\cosh{[\omega(\beta/2a-\tau)]}}{\sinh{[\omega\beta/2a]}}=\frac{2a}{\omega\beta}$
(172)

For the dissipative term, $T_{D}$, in Eq. 167, we then have,

$\displaystyle T_{D}$
$\displaystyle=\lim_{a>>1}\int_{0}^{\beta/a}d\tau^{\prime}\int_{0}^{\beta/a}d\tau
K_{a}(\tau-\tau^{\prime})(q_{a}(\tau)-q_{a}(\tau^{\prime}))^{2}$ (173)
$\displaystyle=\lim_{a>>1}\int_{0}^{\infty}d\omega
J(\omega)\int_{0}^{\beta/a}d\tau^{\prime}\int_{0}^{\beta/a}d\tau\frac{\cosh{[\omega(\beta/2a-(\tau-\tau^{\prime}))]}}{\sinh{[\omega\beta/2a]}}(q(a\tau)-q(a\tau^{\prime}))^{2}$
(174) $\displaystyle\approx\lim_{a>>1}\int_{0}^{\infty}d\omega
J(\omega)\int_{0}^{\beta/a}d\tau^{\prime}\int_{0}^{\beta/a}d\tau\frac{2a}{\omega\beta}(q(a\tau)-q(a\tau^{\prime}))^{2}$
(175) $\displaystyle\approx\lim_{a>>1}\frac{1}{a}\int_{0}^{\infty}d\omega
J(\omega)\int_{0}^{\beta}d\tau^{\prime}\int_{0}^{\beta}d\tau\frac{2}{\omega\beta}(q(\tau)-q(\tau^{\prime}))^{2}$
(176)

Hence, for large $a$, while the potential and dissipative terms scale as
$1/a$, the kinetic term scales as $a$, dominating in high temperature limit
and causing the paths with large variation to contribute less.

Analogous to the USC limit [Appendix. 1], the high temperature limit gives us
the following values for the relevant quantities,

$\displaystyle\lim_{\beta\to 0}\braket{\hat{X}^{2}(q)}$
$\displaystyle=\frac{1}{\beta V_{2}(q)}$ (177) $\displaystyle\lim_{\beta\to
0}\braket{\hat{P}^{2}(q)}$ $\displaystyle=\frac{m}{\beta}$ (178)
$\displaystyle\lim_{\beta\to 0}\mathcal{G}(q)$ $\displaystyle=1$ (179)

Hence the high temperature limit for the expression of LHA [Eq. 37] is given
as,

$\displaystyle\lim_{\beta\to 0}\rho(q,\eta)$ $\displaystyle=\lim_{\beta\to
0}Z^{-1}\exp\left(-\frac{\beta
V_{2}(q)}{2}\left(\frac{V_{1}(q)}{V_{2}(q)}\right)^{2}-\frac{m}{2\beta}\eta^{2}-\beta\left(V(q)-\frac{V_{1}(q)^{2}}{2V_{2}(q)}\right)\right)$
(180) $\displaystyle\implies\lim_{\beta\to 0}\rho(q)$
$\displaystyle=Z^{-1}\exp\left(-\beta V(q)\right)$ (181)

We will now show that this high temperature MFGS is equivalent to the ones
recently derived in literature cresser2021weak ; gelzinis2020analytical ;
timofeev2022hamiltonian . Timofeev et al.timofeev2022hamiltonian  have shown
that in the high temperature limit, at arbitrary coupling, the general
Hamiltonian of mean force, $H_{MF}$, is given as,

$\displaystyle H_{MF}=\sum_{n}\langle n|H_{S}|n\rangle+\sum_{n\neq
m}\bra{n}H_{S}\ket{m}\exp\\{-\beta\Lambda(q_{n}-q_{m})^{2}/6\\}\ket{n}\bra{m}$
(182)

where $\Lambda\equiv\int_{0}^{\infty}d\omega J(\omega)/\omega$ and $q_{n}$ are
the eigenvalues of the system coupling operator $\hat{q}$. Note that Eq. 182
differs from the one obtained by the Timofeev et al. by a term of
$-\Lambda\hat{q}^{2}$ because their definition of the free Hamiltonian
includes the counterterm. For a CV system, Eq. 182 translates into,

$\displaystyle H_{MF}=\left(\int dxV(x)\ket{x}\bra{x}\right)+\left(\int\int
dxdy\frac{\bra{x}\hat{P}^{2}\ket{y}}{2m}\exp\left\\{-\beta\Lambda\frac{(x-y)^{2}}{6}\right\\}\ket{y}\bra{x}\right)$
(183)

Let us first set $\Lambda=0$ for simplicity. Then, we have,

$\displaystyle\implies\exp\\{-\beta\hat{H}_{MF}\\}$
$\displaystyle\approx\exp\\{-\beta
V(x)\\}\times\exp\left\\{-\beta\frac{\hat{P}^{2}}{2m}\right\\}+O(\beta^{2})$
(184)

Let us focus on the second term,

$\displaystyle\bra{x}\exp(-\beta\hat{P}^{2}/2m)\ket{y}$
$\displaystyle=\sqrt{\frac{m}{2\pi\beta}}\exp\left\\{-\frac{m}{2\beta}(x-y)^{2}\right\\}$
(185) $\displaystyle\approx\delta(x-y)\qquad\text{in high temperature limit}$
(186)

We hence recover the high temperature MFGS result [Eq. 181]. Note that when
$\Lambda\neq 0$, then the off-diagonal elements of $H_{MF}$ [Eq. 183] get
further suppressed and it can be formally shown that the exponential of the
second term in Eq. 183 gets closer to $\delta(x-y)$ as compared with the
$\Lambda=0$ case, and hence Eq. 181 remains the high temperature MFGS.

We will now calculate $\Delta(q,\eta)$ [Eq. 50] and $\mathcal{E}(q)$ [Eq. 51]
in the high temperature limit. For $\Delta(q,\eta)$, we first evaluate,

$\displaystyle\lim_{\beta\to
0}\frac{V_{1}(q)}{V_{2}(q)}\left(1-\frac{W(q)}{\braket{\hat{X}^{2}(q)}}\right)$
$\displaystyle=\lim_{\beta\to
0}\frac{V_{1}(q)}{V_{2}(q)}\left(1-\frac{\sum_{n=-\infty}^{+\infty}\frac{(-1)^{n}}{\frac{V_{2}(q)}{M}+\nu_{n}^{2}+\zeta_{n}}}{\sum_{n=-\infty}^{+\infty}\frac{1}{\frac{V_{2}(q)}{M}+\nu_{n}^{2}+\zeta_{n}}}\right)$
(187) $\displaystyle=\lim_{\beta\to
0}\frac{V_{1}(q)}{V_{2}(q)}\left(1-\frac{\frac{M}{V_{2}(q)}}{\frac{M}{V_{2}(q)}}\right)$
(188) $\displaystyle=0$ (189)

We then have,

$\displaystyle\lim_{\beta\to 0}\Delta(q,\eta)=\left|\eta\right|$ (190)

Similarly, for $\mathcal{E}(q)$, we have,

$\displaystyle\lim_{\beta\to 0}\mathcal{E}$ $\displaystyle=\lim_{\beta\to
0}\sqrt{\braket{\hat{q}^{2}}-\frac{1}{m\omega^{2}\beta}}\left(\sqrt{\frac{1}{m\omega^{2}\beta\braket{\hat{q}^{2}}}}+1\right)$
(191) $\displaystyle=2\lim_{\beta\to
0}\sqrt{\frac{2}{m\beta}\sum_{n=1}^{\infty}\frac{1}{\omega^{2}+\nu_{n}^{2}+\zeta_{n}}}$
(192) $\displaystyle=0$ (193)

Since these length scales become smaller as the temperature is increased, we
infer that LHA becomes more accurate in the process.

###  3 High bath cutoff limit

Let $\Omega$ denote the bath cutoff frequency. Then in the limit
$\Omega\to\infty$, we have,

$\displaystyle\lim_{\Omega\to\infty}\int_{0}^{\infty}\frac{J(\omega)}{\omega}d\omega=\infty$
(194)

From the definition of $\zeta_{n}$ [Eq. 42], we have,

$\displaystyle\lim_{n\to\infty}\zeta_{n}$
$\displaystyle=\frac{1}{m}\int_{0}^{\infty}\frac{J(\omega^{\prime})}{\omega^{\prime}}d\omega^{\prime}$
(195) $\displaystyle=\infty$ (196)
$\displaystyle\implies\lim_{n\to\infty}\frac{\omega^{2}+\zeta_{n}}{\omega^{2}+\nu_{n}^{2}+\zeta_{n}}$
$\displaystyle=1$ (197)
$\displaystyle\implies\lim_{n\to\infty}\sum_{n=1}^{\infty}\frac{\omega^{2}+\zeta_{n}}{\omega^{2}+\nu_{n}^{2}+\zeta_{n}}$
$\displaystyle=\infty$ (198) $\displaystyle\implies\braket{\hat{P}^{2}(q)}$
$\displaystyle=\infty$ (199)

Here we have used the definition of $\braket{\hat{P}^{2}(q)}$ [Eq. 39]. Since
$\braket{\hat{P}^{2}(q)}$ diverges, in this limit, LHA [Eq. 37] will not have
coherence in the position basis, and therefore will be given by,

$\rho(q,\eta)=\delta(\eta)Z^{-1}\mathcal{G}(q)\exp\left(-\frac{1}{2\braket{\hat{X}^{2}(q)}}\left(\frac{V_{1}(q)}{V_{2}(q)}\right)^{2}-\beta\left(V{(q)}-\frac{V_{1}{(q)}^{2}}{2V_{2}{(q)}}\right)\right)$
(200)