Thomas Plé

Simon Huppert
email: simon.huppert@sorbonne-universite.fr

Fabio Finocchi

Philippe Depondt

Sara Bonella

Anharmonic spectral features via trajectory-based quantum dynamics: a perturbative analysis of the interplay between dynamics and sampling

The performance of different approximate algorithms for computing anharmonic features in vibrational spectra is analysed and compared on model and more realistic systems that present relevant nuclear quantum effects. The methods considered combine approximate sampling of the quantum thermal distribution with classical time propagation and include Matsubara dynamics, path integral dynamics approaches, linearized initial value representation and the recently introduced adaptive quantum thermal bath. A perturbative analysis of these different methods enables to account for the observed numerical performance on prototypes for overtones and combination bands and to draw qualitatively correct trends for the numerical results obtained for Fermi resonances. Our results prove that the unequal performances of these approaches often derive from the method employed to sample initial conditions and not, as usually assumed, from the lack of coherence in the time propagation. Furthermore, as confirmed by the analysis reported in J. Chem. Phys. 130, 194510 (2021), we demonstrate, both via the perturbative approach and numerically, that path integral dynamics methods fail to reproduce the intensities of these anharmonic features and follow purely classical trends with respect to their temperature behaviour. Finally, the remarkably accurate performance of the adaptive quantum thermal bath approach is documented and motivated.

I. INTRODUCTION

This paper reports a systematic comparison, based on analytical results and numerical investigations, of the performance of a set of approximate methods for simulating quantum nuclear dynamics, with specific focus on the calculation of anharmonic spectral features in infrared spectra. As experiments and technology increasingly investigate phenomena where nuclear quantum effects (NQEs) are at play, the need for accurate and effective algorithms to solve the quantum time-evolution of condensed matter systems becomes apparent. The prototypical example is Hydrogen, the lightest and most abundant element in the Universe, and a basic constituent of inorganic matter and living organisms. It is well known that, for this key element and for molecules containing it, NQEs such as tunnelling, coherence or zero-point energy, are often relevant to describe processes of experimental significance occurring at low temperature 1,2 or high pressure [3][4][5] . Moreover, since at room temperature the De Broglie wavelength of Hydrogen is of the order of the Angström and therefore often comparable to interatomic distances, NQEs are observed also close to ambient conditions. This occurs, for example, in reactions involving proton transfer [6][7][8] or showing significant isotope effects [9][10][11][12][13][14] . The spectroscopy of ice and water is also affected by NQEs [15][16][17] . On a different scene, the rapid developments in the area of quantum computers motivate accurate simulations, for example, of coherence effects in qbits and of the influence of the environment on quantum components.

The key challenge for computing nuclear quantum timedependent properties is that numerically exact solution methods scale exponentially with the number of degrees of freedom. On the other hand, dynamical properties, captured in particular by time correlation functions (TCFs), are relevant as they can be probed through different experimental techniques such as infrared/Raman spectroscopy or measurements of rate constants and transport coefficients. Several approximate schemes have then been designed to mitigate the exponential scaling. In this paper, we will focus on trajectorybased methods as they provide computationally efficient approximations for treating mildly quantum particles. In these approaches, sampling of initial conditions from the quantum thermal distribution is combined with the propagation of classical-like trajectories and observables are computed as appropriate averages over this ensemble of trajectories. These methods have been proven quite effective since, in common conditions of temperature and pressure, the most relevant NQEs derive from zero-point energy and are typically captured, at least approximately, via the initial sampling. Furthermore, in systems containing large numbers of degrees of freedom, the interactions within a complex environment usually reduce the time scales in which other quantum effects, such as coherence, can be observed. This enables a reliable description with methods that use some form of classical propagation to derive short-time approximations of the quantum propagator. For the purposes of this work, we classify trajectorybased methods in two broad families. The first exploits the framework of imaginary-time PI in order to approximate the nuclear quantum dynamics via a ring-polymer representation. The most common representatives of this family are Centroid Molecular Dynamics [18][19][20][21] (CMD) and (Thermostatted) Ring-Polymer Molecular Dynamics 22,23 ((T)RPMD). While these methods were introduced as somewhat ad hoc, they can be put on firmer grounds as approximate forms of the so-called Matsubara dynamics 24,25 , which provides a unified framework to assess their merits. The second family is composed by the semi-classical methods that are derived from stationary phase approximations of the quantum time propagators appearing in the expression of time-correlation functions. One of the simplest and most widely used approaches in this family is the Linearized Semi-Classical Initial Value Representation (LSC-IVR, sometimes also referred to as classical Wigner) 26 . This method relies on the Wigner thermal distribution 27 combined with classical trajectories to rewrite quantum TCFs in a (generalized) phase space, with strong analogies to the classical one. The explicit form of the Wigner function for multidimensional systems is, however, not known and brute-force sampling of its path integral representation is hindered by the presence of phase factors. Consequently, several variants of LSC-IVR have been proposed based on different schemes to approximate and sample the Wigner density as will be discussed in more detail in section III.

Trajectory-based methods have been successfully applied to simulate the behaviour of condensed matter systems, but it is commonly assumed that they are intrinsically unable to capture phenomena such as overtones, combination bands or Fermi resonances because these are believed to derive predominantly from quantum coherence, and thus to be impossible to capture via classical propagation. Only very recently, numerical tests have been performed to compare the performance of different PI or semi-classical methods in reproducing combination bands in the important case of water and ice spectroscopy 16 . These exploratory calculations indicate that an interesting development of CMD, the so-called quasi-CMD 28 (QCMD) method, and LSC-IVR accurately describe the fundamental bands in the spectrum of gas and liquid-phase water and in ice. Furthermore, LSC-IVR can capture the intensity of the high-frequency combination bands relatively well. As mentioned by the authors of Ref. 16, these results raise the question of why ring-polymer methods (QCMD, CMD, and RPMD) are unable to reproduce overtones and combination bands. This failing cannot be attributed to a lack of real-time coherence, since this approximation is also present in LSC-IVR. More recently 29 a perturbation approach was employed to further analyze this effect and show that the inaccuracy of the most common PI methods originates from their improper treatment of the dynamical coupling between the centroid and the Matsubara fluctuation modes.

In this paper, we present a systematic study of the ability of different trajectory-based methods to describe anharmonic spectral features such as overtones, combination bands and Fermi resonances. To that end, we study simple prototype models obtained by adding anharmonic perturbations to a harmonic reference potential. We present a general perturbation approach to derive analytical expressions for the integrated intensity of the relevant spectral features obtained via Matsubara dynamics, CMD, RPMD, and LSC-IVR with different initial conditions sampling schemes. Our results enable to rationalize the different characteristics (and to some extent the reliability) of these approaches. Moreover, we clearly demonstrate that the lack of coherence in the dynamics cannot be invoked to explain inaccuracies in the description of anharmonic spectral features. In particular, we show that the observed trends in the performance of different LSC-IVR variants can be traced back to the treatment of position-momentum correlations, rather than to coherence effects. The general trends established from the perturbative approach are then completed by a numerical exploration of Fermi resonances and of anharmonic spectral features in a more realistic system: gas phase water. As discussed throughout this paper, the results of Ref. 29 -obtained independently from our work -support our analysis for Matsubara and PI based dynamics.

We also present, for each system, numerical results obtained via an alternative, more recent, approach: the adaptive Quantum Thermal Bath (adaptive QTB -adQTB). (ad)QTB combines Newtonian evolution for the nuclei with a generalized Langevin thermostat to include zero-point energy effects 30,31 . The adaptive version of the algorithm was recently proposed 32 to eliminate the so-called zero-point energy leakage problem (i.e. a nonphysical flow of energy driven by the classical evolution). The idea of the QTB was originally introduced as a single particle alternative to PI for computing static averages in the quantum canonical ensemble. Similarly to RPMD and CMD, it has then been proposed also as an approximate tool to describe dynamical properties, and in particular vibrational spectra 4,5,[START_REF] Mauger | Nuclear quantum effects in liquid water at near classical computational cost using the adaptive quantum thermal bath[END_REF] . Providing an analytical justification for this type of calculations is problematic, but we demonstrate numerically that adQTB spectra are remarkably accurate and can therefore offer an interesting and cost effective alternative to PI-based simulations.

The paper is organized as follows. Section II, presents the perturbation strategy that is used to derive analytical expressions for classical and trajectory-based quantum approximations of time correlation functions. In section III, we focus on the problem of initial conditions sampling for the LSC-IVR approach and present three different approximations to the Wigner thermal density. Section IV and V analyze the results obtained from perturbation theory on simple prototype models for overtones and combination bands, respectively. In section VI, we explore numerically the performance of the different methods to describe the Fermi resonance phenomenon, while section VII illustrates how the trends derived from perturbation theory can provide qualitative and quantitative insight for the description of a more realistic system: gas-phase water. Finally, section VIII summarizes the main conclusions of this study.

II. TRAJECTORY-BASED QUANTUM TIME CORRELATION FUNCTIONS

In this section we present the perturbative approach employed to derive analytical expressions for the momentum autocorrelation functions (ATCF) from which the vibrational spectra can be computed, and summarize the approximate quantum approaches that we shall examine. To better illustrate our perturbative analysis, we first present it in the relatively simple case of the purely classical ATCF and then introduce the main adjustments necessary to treat the trajectorybased approximations of the quantum dynamics. Unidimensional notations are used for simplicity but the procedure can equally be applied to multidimensional systems.

A. Perturbative expansion to the classical dynamics For a classical system at thermal equilibrium, the momentum ATCF is defined as:

c pp (t) = dqdp ρ(q, p) pe iL t p (1)
with ρ(q, p) = e -β H(q,p) /Z c the Boltzmann density, and iL the classical Liouville operator. This operator characterizes the Newtonian dynamics so that e iL t p is the momentum obtained after propagating this dynamics for a duration t, starting from the phase-space point (q, p). We consider systems described by a potential of the following form:

V (q) = V 0 (q) + λ v(q) (2)
composed of a main harmonic component V 0 (q) plus an anharmonic time-independent perturbation proportional to the control parameter λ . We then expand equation ( 1) to second order in λ (the first-order contribution is zero for the anharmonic spectral features on which we focus). The expansion of the probability density is:

ρ(q, p) = ρ 0 (q, p) 1 + λ β v(q) -v(q) 0 + λ 2 β 2 2 v(q) 2 + v(q) 2 0 -2v(q) v(q) 0 + O λ 3 (3)
where the symbols . . . 0 denote the averaging over the Gaussian distribution ρ 0 (q, p) associated to the harmonic potential V 0 . The Liouvillian is split as:

iL = iL 0 -λ dv dq ∂ ∂ p (4)
with iL 0 the harmonic Liouvillian, for which time propagation can be performed analytically. The perturbed time propagator can be expanded using the Dyson series:

e iL t = e iL 0 t -λ t 0 ds e iL 0 (t-s) dv dq ∂ ∂ p e iL 0 s + λ 2 t 0 ds 1 s 1 0
ds 2 e iL 0 (t-s 1 ) dv dq

∂ ∂ p e iL 0 (s 1 -s 2 ) dv dq ∂ ∂ p e iL 0 s 2 + O λ 3 (5)
We then inject both expansions (3) and ( 5) in the expression (1) of the ATCF and truncate terms of order higher than 2 in λ (more details can be found in Appendix A). In sections IV and V, this approach is applied to polynomial perturbations v(q). We use symmetries and the properties of the harmonic Liouvillian iL 0 and Gaussian distribution ρ 0 in order to perform the time and phase-space integrals and we derive the secondorder perturbative expansion of c pp (t). Notably, since the unperturbed Liouvillian iL 0 is purely harmonic, its time propagation does not produce anharmonic resonances. As a result, the first order in the expansion (3) of the density ρ(q, p) is sufficient to describe overtones and combination bands up to second order in λ .

B. Linearized approximations

The extension of the perturbative strategy to the linearized IVR approaches is relatively straightforward. We start by expressing the real part of the standard quantum ATCF as 24,[START_REF] Hillery | [END_REF] :

C pp (t) = Re[C pp (t)] = dqdp ρ w (q, p) pe iL Q t p (6)
where C pp is the (complex) standard TCF and ρ w (q, p) is the thermal Wigner distribution:

ρ w (q, p) = 1 2π hZ d∆ e ip∆ h q - ∆ 2 e -β Ĥ q + ∆ 2 (7)
with Z the quantum partition function of the system. The quantum Liouville operator iL Q reads:

iL Q = p m ∂ ∂ q - ∞ ∑ ν=1,odd 1 ν! ih 2 ν-1 d ν V dq ν ∂ ν ∂ p ν (8)
The real part of the quantum ATCF is related to the Kubotransformed ATCF by the following equation in the Fourier domain:

K pp (ω) = C pp (ω) β Θ(ω, β ) (9)
with Θ(ω, β ) = hω/2 tanh(β hω/2) the average thermal energy of a quantum harmonic oscillator of frequency ω. The Kubotransformed momentum ATCF characterizes the vibrational density of state (VDOS) of a quantum system, which is given by β m K pp (ω). This quantity can be directly compared with results from exact and path integral dynamics methods. For linear observables such as the ones used in this work, equation ( 6) is the exact expression for the real part of the Wigner ATCF and does not imply any approximation at this point [START_REF]non-linear, observables, the exact Wigner representation of the TCF of operators  and B involves the Moyal product 61,62 between the Wigner density and the Wigner transform of Â, defined as [ ρ Â] w (q, p) = ρ w (q, p)A w (q, p) + (ih/2)(∂ q ρ w (q, p)∂ p A w (q, p) -∂ p ρ w (q, p)∂ q A w (q, p)) + ... It has been shown[END_REF] . This expression, however, requires to perform the time-evolution using iL Q , an impossible task for realistic systems. The LSC-IVR approximation is then obtained by replacing iL Q by its classical counterpart iL (i.e. by truncating the sum in equation ( 8) at order ν = 1). This classical evolution can be integrated numerically using standard molecular dynamics methods. The approximation eliminates any effect associated with quantum coherence and interference.

In sections IV and V, we consider 3rd-order polynomial potentials V (q), so that the sum in equation ( 8) is actually finite and contains only one additional term with respect to iL , corresponding to ν = 3. It is then possible to apply the perturbative approach to both the exact momentum ATCF of ( 6) and its LSC-IVR approximation. By comparison, we show that, up to second order in λ , the additional term in iL Q does not contribute to the the momentum ATCF and thus to the intensity of the anharmonic resonances in the vibration spectrum, computed as its Fourier transform. Importantly, this implies that, for these resonances, the accuracy of LSC-IVR results does not depend on the time propagation but only on the sampling of the Wigner thermal density. In particular, as discussed more in detail in the following, if the first order term in the λ -expansion of ρ W is correct, the observables will be exact within the perturbative analysis.

C. Path integral dynamics methods

We also consider three other approaches that are based on the imaginary-time path integral framework: Matsubara dynamics, RPMD and CMD. These methods directly provide approximations to the Kubo-transformed ATCF. In each of the three PI-based methods, the expression for K pp (t) is similar to equation (1), except that the probability density and the Liouville operator are expressed in the extended PI phase space 24,25 . The perturbative strategy presented above, based on the expansion of the density in powers of λ and on the Dyson series for the Liouville operator, can therefore also be applied. The derivation of the perturbative expansion of the ATCF then makes use of the analytical expressions that can be derived for harmonic systems in the PI framework, as detailed in appendix B. This procedure is technically more involved, but conceptually very similar to the classical calculation. The results derived in this work for PI-based methods are further corroborated by the calculations by Benson and Althorpe in Ref. 29, via a different approach based on classical perturbation theory expressed in terms of action-angle variables.

D. Quantum thermal bath

We also include in our comparisons the quantum thermal bath approach 30 . This method (as the closely related quantum thermostat 31 ) relies on a generalized Langevin equation to approximate statistical NQEs:

d p dt = - ∂V ∂ q -γ p + F(t) (10)
where γ is a friction coefficient and F(t) is a colored random force with the following correlation spectrum:

C FF (ω) = 2mγ r (ω)Θ(ω, β ) (11)
with γ r (ω) the random force amplitude that is simply given by γ r (ω) = γ in the original QTB method. The QTB acts as a frequency-dependent thermostat that tends to thermalize the vibration modes of the system with an average energy Θ(ω, β ) -therefore accounting for zero-point contributionsinstead of the classical thermal energy β -1 . Though simple and efficient, the QTB suffers from zero-point energy (ZPE) leakage: the ZPE injected at high frequencies tends to unphysically flow towards low frequencies. In Ref. 32, a solution to this issue was proposed by noting that in a Langevin simulation, the linear susceptibility χ qp (ω) is related to the momentum-random force correlation function by:

χ qp (ω) = Re [C pF (ω)] C FF (ω) (12)
In generalized Langevin simulations, there are therefore two different estimators for the VDOS, that can either be computed from C pp (ω) via eq. ( 9), or from the linear susceptibility estimated via [C pF (ω)]. These two estimators are equivalent if the quantum fluctuation-dissipation theorem [START_REF] Kubo | [END_REF] holds:

Re [C pF (ω)] = γ r (ω)C pp (ω) (13)
However, as shown in Ref. 32, this relation is not generally respected in standard QTB simulations. Indeed, deviations from (13) can be used as a criterion for a quantitative assessment of ZPE leakage. The adaptive QTB (adQTB) method exploits this criterion to adjust the random force amplitude γ r (ω) on the fly in order to compensate for ZPE leakage in a systematic way 32 and restore the fluctuation-dissipation theorem, eq. ( 13). Beyond the assessment of ZPE leakage, we show in the following that the fluctuation-dissipation theorem also proves useful in calculating anharmonic resonances.

III. SAMPLING THE WIGNER DISTRIBUTION

As we noted in section II B, the accuracy of the LSC-IVR approximation for the prediction of anharmonic resonance intensities depends entirely on the quality of the Wigner distribution sampling. This sampling represents a difficult problem for realistic high-dimensional systems as the rapidly oscillating phase factor e ip∆ h in equation ( 7) prevents a direct calculation of the integral. Furthermore, ρ W can take negative values and it is therefore not strictly a probability distribution. Different approximations have been designed to allow for a practical sampling of ρ W , usually relying on (generalized) Gaussian functional expressions 37,38 . In this paper, we analyze three of these approximations: (i) an equilibrium harmonic approximation (that we will denote IVR0), (ii) the Local Harmonic Approximation 39,40 (LHA) and (iii) the Edgeworth Conditional Momentum Approximation 41 (ECMA). In these methods, the Wigner distribution is split as the product of the marginal position probability density and the conditional momentum (pseudo)distribution:

ρ w (q, p) = q| e -β Ĥ Z |q × d∆ e i p∆ h q -∆ 2 e -β Ĥ q + ∆ 2 2π h q|e -β Ĥ |q = ρ m w (q) × ρ c w (p|q) (14)
The marginal position density ρ m w (q) = q| e -β Ĥ Z |q can be sampled exactly using standard PI molecular dynamics or Monte Carlo methods, and the three LSC-IVR approximations that we consider differ only in their expression for ρ c w (p|q).

A. Equilibrium harmonic approximation IVR0

The equilibrium harmonic approximation is obtained by expanding the potential operator V ( q) to second order around the configuration q 0 of minimal potential energy:

V ( q) ≈ V (q 0 ) + 1 2 mΩ 2 0 ( q -q 0 ) 2 (15)
with m the mass of the particle and where the squared frequency Ω 2 0 is defined as:

Ω 2 0 = 1 m d 2 V dq 2 q 0 ( 16
)
As we will apply this method to potentials that are primarily harmonic, there is no ambiguity on the definition of q 0 . We note that the expansion of the potential is only performed in the conditional momentum part of the Wigner density (the marginal position distribution is still exact). The approximate potential (15) is quadratic so that matrix elements of the density operator e -β Ĥ can be written explicitly (see for example ref. 42). The IVR0 momentum distribution is then independent of position and given by:

ρ c IVR0 (p) = exp - p 2 2mΘ(Ω 0 ,β ) dp exp - p 2 2mΘ(Ω 0 ,β ) (17)
The IVR0 expression is the simplest Gaussian approximation to ρ c w and, contrary to LHA and ECMA (see below), the corresponding momentum distribution is not position-dependent. We will show in sections IV and V that this poses a problem, since a precise treatment of the intrinsically quantum and anharmonic position-momentum correlation is actually critical to correctly capture anharmonic resonances.

B. Local Harmonic Approximation

The local harmonic approximation 39,40 is similar to IVR0 except that the quadratic expansion of the potential is performed locally around each position q. This defines the local squared frequency:

Ω 2 (q) = 1 m d 2 V dq 2 (18)
The LHA conditional momentum distribution is then given by:

ρ c LHA (p|q) = exp - p 2 2mΘ(Ω(q),β ) dp exp - p 2 2mΘ(Ω(q),β ) (19)
In this conditional probability, some correlation between momenta and position is present via the dependence on q of Θ(Ω(q), β ). This correlation, however, is purely local in the particle's position and, as we shall see in the following, this can be a major drawback of the LHA. In sections VI and VII, to deal with potentials beyond perturbed harmonic oscillators, where imaginary frequencies may arise, we use the Local Gaussian Approximation (LGA) described in ref. 40 which generalizes the LHA to cases where Ω 2 (q) can be locally negative.

C. Edgeworth Conditional Momentum Approximation

Building on the formalism proposed in ref. 43, the ECMA uses the Edgeworth expansion of ref. 41 to approximate the conditional momentum distribution (details in appendix C).

The resulting distribution is expressed as a position-dependent Gaussian (similarly to LHA) multiplied by a correction factor:

ρ c ECMAn (p|q) = exp -κ 2 (q)p 2 2h 2 C EWn (q, p) dp exp -κ 2 (q)p 2 2h 2 C EWn (q, p) (20)
where κ 2 (q) (and more generally κ m (q)) denotes the second order (m th order) cumulant of the distribution f (∆|q) defined as:

f (∆|q) = q -∆ 2 e -β Ĥ q + ∆ 2 d∆ q -∆ 2 e -β Ĥ q + ∆ 2 ( 21
)
For each position q, κ 2 (q) is computed via an auxiliary openchain path integral calculation, following the procedure of ref. 41. Therefore, in the ECMA, the position-dependent width of the Gaussian momentum distribution is the result of a fully non-local calculation and it can differ significantly from that obtained in LHA from the local curvature of the potential. The Edgeworth correction factor takes the form of a truncated expansion in powers of p:

C EWn (q, p) = 1 + n ∑ m=4,even κ m (q) m! ip h m (22)
It can be positive or negative, which allows recovering the negative values of ρ w , a feature that is totally absent in both IVR0 and LHA. Furthermore, the ECMA converges asymptotically to the exact Wigner distribution in the limit of large truncation orders n. In most cases, ECMA0 already provides a very good approximation of the Wigner distribution, that can be further improved by including the first few terms of the correction factor 41 . Practical sampling of the ECMA goes as follows: (1) sample the exact marginal position distribution using standard PIMD. (2) For each position sample q, compute κ 2 (q) and the successive cumulants of f (∆|q) via an auxiliary (open chain) PI calculation. (3) Sample the Gaussian part of the conditional momentum distribution ρ κ 2 (p|q) ∝ exp -κ 2 (q)p 2 /2h 2 . (4) For n ≥ 4, reweight the samples with the normalized Edgeworth correction C EWn (q, p)/ C EWn (q, p) ρ κ 2

. Appendix C provides more details on each of these sampling steps.

The three sampling schemes presented here were chosen as representatives of increasingly accurate methods but other procedures exist such as the Feynman-Kleinert variational linearized path integral method 44 that also introduces approximate position-momentum correlations.

IV. VIBRATIONAL OVERTONES IN THE PERTURBATIVE REGIME

The simplest model for investigating overtones is given by the following potential:

V (q) = 1 2 mω 2 0 q 2 + λ 3 q 3 ( 23
)
which consists of a harmonic oscillator at frequency ω 0 and a cubic perturbation term.

The dynamics of this model system is mainly harmonic with a small anharmonic contribution -the overtone -at angular frequency 2ω 0 . In this section, we study the temperature dependence of the intensity of this overtone, as obtained from different approximate quantum dynamics methods. To quantify this intensity, we introduce the parameter η, defined from the integrated overtone VDOS:

η(β ) = 2ω 0 +∆ω 2ω 0 -∆ω β m K pp (ω; β )dω (24)
with ∆ω a frequency interval chosen to encompass the whole overtone spectral region. The parameter η(β ) is dimensionless and measures the fraction of the system vibration energy that is contained in the overtone.

For small enough λ , the perturbative strategy outlined in section II can be applied to derive analytical expressions for K pp (ω; β ) and η(β ). In the classical case, the VDOS is related to the classical ATCF c pp (instead of the Kubotransformed K pp ), and, following the approach outlined in section II, the overtone contribution η cl is obtained, at second order in λ , as:

η cl (β ) = 4λ 2 9m 3 ω 6 0 β + O λ 3 (25)
This quantity varies linearly with the temperature T and in particular, at low T , the classical overtone disappears as the system freezes in the vicinity of the potential minimum where the anharmonicity is negligible. For the model system (23), it is also possible to derive a perturbative expression for the quantum VDOS and for the corresponding overtone intensity η Q (β ):

η Q (β ) = 4λ 2 9m 3 ω 6 0 Θ(ω 0 , β ) + O λ 3 ( 26
)
which can be derived in two different ways, either using standard Rayleigh-Schrödinger perturbation theory based on the harmonic oscillator eigenstates 29,45 , or from the perturbative expansion of the exact Wigner distribution provided in appendix E. Equation ( 26) is similar to (25), but the classical thermal energy 1/β is replaced by the quantum energy Θ(ω 0 , β ). As shown in Fig. 1, in the high temperature limit, the classical and quantum results coincide. However, at low temperature, Θ(ω 0 , β ) → hω 0 /2 and thus, contrary to its classical counterpart, η Q tends to a non-zero constant, due to zeropoint motion contributions.

In the following subsections, we use the perturbative analysis introduced in section II to investigate the performance of some commonly adopted trajectory-based methods for approximating the quantum η(β ). We first consider the linearized semi-classical initial value representation (LSC-IVR) and emphasize the importance of the initial conditions by considering the different approximations to the initial Wigner distribution introduced in section III. We then turn to ringpolymer-based approximations, before finally examining the results of the QTB method. Our findings are summarized in Fig. 1. 23) computed using different methods. Lines are perturbative expressions while bullet points are simulation results. In all the calculations, m is the proton mass, ω 0 = 1500 cm -1 and λ = 21 kcal/mol/ 3 , which is small enough for the perturbative analysis to be highly accurate.

A. Linearized dynamics

As we pointed out in section II, for the simple cubic perturbation of equation ( 23), the classical and quantum propagation of p coincide up to second order in λ . A Dyson expansion of the quantum Liouville time propagator gives the following propagation of p: e iL Q t p = e iL t p + λ h2 12 t 0 ds e iL 0 (t-s) ∂ 3 ∂ p 3 e iL 0 s p

+ λ 2 h4 144 t 0 ds 1 s 1 0 ds 2 e iL 0 (t-s 1 ) ∂ 3 ∂ p 3 e iL 0 (s 1 -s 2 ) ∂ 3 ∂ p 3 e iL 0 s 2 p - λ 2 h2 12 t 0 ds 1 s 1 0 ds 2 e iL 0 (t-s 1 ) q 2 ∂ ∂ p e iL 0 (s 1 -s 2 ) ∂ 3 ∂ p 3 e iL 0 s 2 p - λ 2 h2 12 t 0 ds 1 s 1 0 ds 2 e iL 0 (t-s 1 ) ∂ 3 ∂ p 3 e iL 0 (s 1 -s 2 ) q 2 ∂ ∂ p e iL 0 s 2 p + O λ 3 ( 27
)
Using the analytic expressions available for the propagation according to the harmonic operator e iL 0 t , it is easy to show that the third derivative ∂ 3 ∂ p 3 , that appears in the first and second order terms of eq. ( 27) is always null, so that e iL Q t p = e iL t p + O λ 3 (28) Therefore, the LSC-IVR approximation to the momentum ATCF is exact at second order in perturbation. Quantum effects thus only arise from the sampling of initial conditions, and the accuracy of the LSC-IVR estimate for η(β ) is entirely determined by the approximation used for the Wigner distribution ρ w (q, p). More specifically, we show in the following that accounting for momentum-position correlations is essential to describe the overtone accurately.

As we noted above, the first-order expansion of the density ρ w (q, p) is sufficient to describe overtones up to second order in λ . For the model of equation ( 23), the perturbative expansion of the IVR0, LHA, and ECMA0 approximations for ρ c w (p|q) can all be expressed as:

ρ c IVR (p|q) = ρ 0 (p) 1 + α(ω 0 , β ) λ q 3mω 2 0 p 2 mΘ(ω 0 , β ) -1 + O λ 2 (29)
where ρ 0 (p) is the Gaussian momentum distribution of the harmonic oscillator, and α(ω 0 , β ) is a coefficient that differs between the three methods and modulates positionmomentum correlations.

In the IVR0 approximation, α(ω 0 , β ) = 0, which yields the integrated overtone in the VDOS:

η IVR0 (β ) = η Q (β ) Θ(ω 0 , β ) Θ(2ω 0 , β ) + O λ 3 ( 30
)
The approximation is therefore asymptotically correct in the classical high-temperature limit, where the factor that multiplies η Q tends to one, but η IVR0 → 1 2 η Q at low T . The neglect of position-momentum correlations in the initial distribution thus leads to a global underestimation of the overtone intensity.

For non-zero α(ω 0 , β ), the general form of the overtone intensity can be shown to be:

η IVR (β ) = η Q (β ) Θ(ω 0 , β ) Θ(2ω 0 , β ) 1 + α(ω 0 , β ) (31)
In LHA,

α LHA (ω 0 , β ) = 3 2 1 + β Θ(2ω 0 , β ) -2β Θ(ω 0 , β ) (32)
This approximation tends to overestimate α(ω 0 , β ) due to the local nature of the approximation for the Gaussian conditional momentum distribution. As a consequence, the overtone intensity is higher than the exact quantum result, in particular in the zero-temperature limit, where η LHA → 5 4 η Q . Finally, in the ECMA0 approximation,

α ECMA0 (ω 0 , β ) = Θ(2ω 0 , β ) Θ(ω 0 , β ) -1 (33)
which coincides with the first order expansion of the exact Wigner density for the potential (23). This is sufficient to ensure that the LSC-ECMA0 approximation for C pp and for the overtone intensity η(β ) is exact at second order in λ . Moreover, higher order terms in the Edgeworth expansion cancel at this order of perturbation (see appendix E). Therefore the fully non-local description of the position-momentum correlations provided by ECMA allows for an accurate description of the overtone intensity η(β ), even though the linearization of the dynamics suppresses any quantum coherence effects.

B. Path integral based approximations

The perturbative analysis detailed in appendix B reveals that, at second order in λ , Matsubara dynamics yields the exact overtone intensity η Q (β ) of equation (26). In contrast, the results of its most common approximations -RPMD and CMD -equal that of classical simulations η cl (β ). We can analyze this difference as follows: in Matsubara dynamics, the overtone is produced from constructive anharmonic interactions between the centroid of the ring-polymer and the fluctuation Matsubara modes. In CMD, this interaction is averaged out so that only the classical centroid contribution remains in the overtone. On the other hand, the spring terms in the RPMD potential shift the frequencies of the fluctuation Matsubara modes so that they do not interact constructively with the centroid at the overtone frequency but produce spurious peaks instead. Therefore, in RPMD as well, only the centroid contributes to the overtone. The dynamical interaction between the centroid and the fluctuation modes is thus essential to correctly capture anharmonic resonance effects. This analysis is confirmed by the calculations of Benson et al. that, using perturbation theory in the action angle variables representation of the Matsubara dynamics, reached the same conclusions 29 .

C. Numerical experimentation with the Quantum Thermal Bath

Fig. 1 presents also numerical experiments performed with the quantum thermal bath. Remarkably, in a standard QTB dynamics, the overtone intensity depends quite strongly on the estimator chosen for the VDOS. In Fig. 1, we show results for the VDOS obtained from the C pp (ω) (results indicated as QTB-a) and from C pF (ω) (QTB-b). The QTB-a results are in close agreement with the LSC-IVR0 curve as can be explained by the following qualitative argument. Similarly to LSC, QTB combines an approximate sampling of the Wigner distribution (it is exact for purely harmonic systems) with an essentially classical dynamics (at least on timescales shorter than γ -1 ). Furthermore, similarly to the IVR0, the QTB was shown not to capture position-momentum correlations 46 . The QTB-b results, on the other hand, closely follow the exact η Q (β ) curve. This agreement is remarkable and indicates that, contrary to the C pp (ω) ATCF, the estimation of C pF (2ω 0 ), being proportional to χ qp (2ω 0 ), directly probes the system's response at the overtone frequency and provides an accurate estimate of the anharmonic intensity.

For this prototype model, we performed standard QTB simulations, in which the two estimators for the VDOS can differ. In the more realistic cases presented in sections VI and VII, we used the adQTB algorithm, in which case the two estimators coincide (see appendix D).

V. COMBINATION BANDS

In section IV, we treated a simple unidimensional overtone model, but in molecular and condensed matter systems, anharmonic resonances also arise from the coupling between different vibration modes. In particular, when two modes are coupled, combination bands (CB), corresponding to the sum and the difference of the harmonic frequencies, appear in the spectrum. This behavior can be modeled by using the following two-dimensional potential:

V (q 1 , q 2 ) = 1 2 m 1 ω 1 q 2 1 + 1 2 m 2 ω 2 q 2 2 + λ q 1 q 2 2 ( 34
)
consisting of two harmonic oscillators at frequencies ω 1 and ω 2 coupled by a cubic interaction term with amplitude con-trolled by the parameter λ . Three anharmonic resonances appear in the momentum autocorrelation spectrum of potential (34): an overtone at 2ω 2 and combinations bands at ω 1 + ω 2 and ω 1ω 2 (the overtone contribution at 2ω 1 is zero at second order in λ ). We characterize the different resonances by the intensities η(β ; 2ω 2 ), η(β ; ω 1 + ω 2 ) and η(β ; ω 1ω 2 ) of their respective contributions to the vibrational density of states.

The perturbative analysis shows that, in the classical case, these intensities are given by:

η cl (β ; ω 1 + ω 2 ) = λ 2 (ω 1 + ω 2 ) 2 β -1 m 1 m 2 2 (ω 1 + 2 ω 2 ) 2 ω 1 4 ω 2 2 + O λ 3 ( 35
)
η cl (β ; ω 1 -ω 2 ) = λ 2 (ω 1 -ω 2 ) 2 β -1 m 1 m 2 2 (ω 1 -2 ω 2 ) 2 ω 1 4 ω 2 2 + O λ 3 ( 36
)
η cl (β ; 2ω 2 ) = 2 λ 2 β -1 m 1 m 2 2 (ω 1 + 2 ω 2 ) 2 (ω 1 -2 ω 2 ) 2 ω 2 2 + O λ 3 (37)
As in section IV, the combination bands are second order in λ and vary linearly with temperature. Moreover, the configuration ω 1 = 2ω 2 appears as a particular case (the Fermi resonance) where the intensities diverge, indicating that the perturbative expansion of section II does not apply. The specific case of the Fermi resonance is examined numerically in section VI.

In the quantum case, the exact contributions of the resonances to the VDOS are obtained perturbatively as:

η Q (β ; ω 1 + ω 2 ) = (ω 1 Θ (ω 2 , β ) + ω 2 Θ (ω 1 , β ))λ 2 (ω 1 + ω 2 ) m 1 m 2 2 (ω 1 + 2 ω 2 ) 2 ω 1 4 ω 2 2 + O λ 3 (38) η Q (β ; ω 1 -ω 2 ) = (ω 1 Θ (ω 2 , β ) -ω 2 Θ (ω 1 , β ))λ 2 (ω 1 -ω 2 ) m 1 m 2 2 (ω 1 -2 ω 2 ) 2 ω 1 4 ω 2 2 + O λ 3 (39) η Q (β ; 2ω 2 ) = 2 λ 2 Θ (ω 2 , β ) m 1 m 2 2 (ω 1 + 2 ω 2 ) 2 (ω 1 -2 ω 2 ) 2 ω 2 2 + O λ 3 (40)
These expressions are very similar to the classical ones (eqs. 35-37), except that the factor β -1 at the numerator is replaced with weighted combinations of Θ(ω 1 , β ) and Θ(ω 2 , β ). Indeed, contrary to the classical case where equipartition of energy holds, in the quantum case, the different amounts of ZPE in each mode affects the combination bands amplitude.

The comparison, as a function of temperature, between the quantum (solid black curve) and classical results (dashed black curve) is shown in Figure 2. We clearly observe two different behaviors: η Q (ω 1 + ω 2 ) and η Q (2ω 2 ) saturate at a nonzero value at low temperature due to zero-point energy contri-butions. On the other hand, η Q (ω 1ω 2 ) goes to zero at low temperature faster than the classical system. This can be explained as follows. In the Rayleigh-Schrödinger perturbation picture, the difference-frequency component of the quantum ATCF arises from transitions that start from an excited vibrational state. At low temperature, when only the ground state is populated, these transitions are efficiently suppressed and the corresponding resonance vanishes. The figure also presents the perturbative expressions derived below for the different LSC-IVR approximations as well as the results of numerical experiments with the QTB and RPMD methods. In this case, the combination band intensities are obtained by inte-gration of the VDOS over a frequency window chosen to encompass the appropriate peak while excluding contributions from other peaks. As in section IV, we choose λ so that the second-order perturbation expansion is accurate. For the comparisons reported here, we use ω 1 = 2.5ω 2 = 1500 cm -1 and λ = 8.5 kcal mol -1 Å -3 .

A. Linearized dynamics

Using the Dyson expansion and similar arguments as in section IV.A (eq. 27), it is possible to show that, for the cubic perturbation of eq. ( 34), quantum propagation of the momenta is equivalent to classical propagation up to second order in λ . Therefore a relation similar to (28) holds for each degree of freedom and quantum effects only stem from the initial Wigner distribution in the LSC-IVR framework.

As in the one-dimensional case of section IV, ECMA0 yields the exact Wigner distribution at first order in λ so that the anharmonic resonance intensities are correct up to second order. On the other hand, in the IVR0 approximation, momentum-momentum and position-momentum correlations are neglected leading to incorrect intensities of the combination bands. These are given by:

η IVR0 (β ; ω 1 + ω 2 ) = η Q (β ; ω 1 + ω 2 ) Θ(ω 2 , β ) Θ(ω 1 + ω 2 , β ) (41)
η IVR0 (β ; ω 1 -ω 2 ) = η Q (β ; ω 1 -ω 2 ) Θ(ω 2 , β ) Θ(ω 1 -ω 2 , β ) (42)
η IVR0 (β ; 2ω 2 ) = η Q (β ; 2ω 2 ) Θ(ω 1 , β ) Θ(2ω 2 , β ) (43)
In LSC-IVR0, the quantum result η Q is thus multiplied by the ratio between the harmonic thermal energy Θ at frequency ω 2 (or ω 1 for eq. ( 43)) and at the frequency of the combination band considered. This can lead to non-negligible errors, particularly in the evaluation of the sum-frequency component, as shown in Fig. 2. The perturbative form of the LSC-LHA anharmonic intensities is more involved and reported in appendix F. Numerical results for this approach are shown in Figure 2 and that the local approximation for the position-momentum and momentum-momentum correlations improves the estimate of the intensity of the overtone at 2ω 2 with respect to IVR0, but it leads to an overestimation of the sum-frequency combination band at ω 1 + ω 2 (similarly to the observations in section IV). More intriguing is the behavior of η LHA (ω 1ω 2 ), that becomes negative at low temperature. Indeed, in numerical experiments with the potential of eq. ( 34), the momentum autocorrelation spectrum computed in LHA displays a negative dip at frequency ω 1ω 2 . This result is unphysical since, according to the Wiener-Khinchin theorem, the momentum autocorrelation spectrum is always positive for stationary dynamics. Therefore, the negative intensity of the differencefrequency feature is a clear manifestation that the classical dynamics does not conserve the initial distribution.

B. Other approximations

Using the perturbative analysis, we were able to show that Matsubara dynamics yields the exact anharmonic resonance intensities up to second order in λ . On the other hand, both RPMD and CMD yield classical intensities. Numerical results obtained in RPMD simulations are shown in figure 2 as an illustration. These findings are in agreement with those of Benson and Althorpe in ref. 29 and can be interpreted similarly as in section IV B.

The QTB numerical experiments are also consistent with the trends reported in IV: the anharmonic intensities closely follow the exact curve when the VDOS is computed from the linear susceptibility χ(ω) via the C pF (ω) correlation spectrum (QTB-b), whereas the naive evaluation from the C pp (ω) autocorrelation spectrum (QTB-a) leads to significant inaccuracies.

Finally we mention that eqs. 35-40 can also be used to derive a closed-form quantum correction factor (QCF) to be applied to classical vibration spectra in order to recover the quantum intensity of overtones and combination bands. QCFs have been explored before to approximate NQEs in the calculations of different dynamical observables [47][48][49] . Use of these factors, however, can be problematic when dealing with vibrational features since different choices of QCF can lead to very different estimates 39,50,51 . Furthermore, these factors may amplify the effect of numerical noise, especially on the signal at high frequencies. Nonetheless, Benson and Althorpe show in Ref. 29 that QCFs derived from perturbation theory can provide satisfactory results for simple molecules, among which the water model presented in Section VII.

VI. FERMI RESONANCES

In section V, we noted that the configuration ω 1 ≈ 2ω 2 appears as a particular case, the Fermi resonance, in which the perturbative approach developed above for combination bands breaks down. In this configuration, the interaction between the 2ω 2 -overtone and the ω 1 -mode causes the high-frequency peak in the vibration spectrum to split. In ref. 52, Basire et al. developed a different perturbative approach and showed that the temperature-dependence of this splitting is qualitatively different if the system is treated classically or quantummechanically, accounting for zero-point motion. Indeed, the classical Fermi splitting goes to zero as T 1/2 at low temperature while the quantum one saturates at a non-zero value. This behavior can be related to that of the overtone intensity in section IV: because of zero-point fluctuations, the 2ω 2overtone does not disappear at low temperatures and its interaction with the high-frequency mode and the resulting splitting persist. To study the Fermi resonance, we use the simple 3-dimensional model for gas-phase CO 2 -the historic example for which Fermi developed the theoretical explanation of Gaussian function so that the height of the peaks approximately match the other methods. Center panel: Fermi splitting as a function of temperature approximated with TRPMD, LSC-ECMA and adQTB. This splitting is measured as the distance between the center of mass of each of the two peaks that constitute the Fermi dyad. Arbitrary units are used for the intensities of the spectra.

this phenomenon 53 -also adopted in ref. 52:

V (x, y, z) = 1 2 m 1 ω 2 1 x 2 + 1 2 m 2 ω 2 2 y 2 + 1 2 m 2 ω 2 2 z 2 + 1 2 χ 12 x(y 2 + z 2 ) √ m 1 m 2 ( 44
)
where x is the symmetric stretching coordinate and y, z correspond to the two degenerate bending modes. The antisymmetric streching mode is on a different frequency range and does not interact with the relevant modes; it is thus not included in the model. The masses m 1 and m 2 are chosen as the reduced masses of the corresponding modes and the harmonic frequencies are ω 1 = 1261 cm -1 and ω 2 = 634 cm -1 . The anharmonic coupling parameter is taken equal to χ 12 = 1.479 × 10 -7 a.u.. Figure 3 shows the Fermi splitting obtained in this model using TRPMD (blue), LSC-ECMA (red) and adQTB (green) compared to exact results (grey) obtained by solving the Schrödinger equation on a basis of harmonic wavefunctions. For this low-dimensional system, the exact spectrum consists in a sum of Dirac δ -functions, centered at the frequencydifferences between pairs of eigenstates of the system. For visualization purposes, these peaks are displayed with an artificial (and temperature-independent) broadening in the left and in the right panels of Fig. 3. Note that for this system, use of the adaptive version of the QTB was critical due to the massive zero-point energy leakage occurring in the Fermi resonant configuration 32 . Moreover, the thermostatted RPMD algorithm was preferred to prevent the appearance of spurious peaks in the low temperature spectra. At high temperature T = 800 K (right panel), all simulation methods agree in showing a doublet with comparable splitting and peak shapes. At 150 K (left panel), TRPMD clearly underestimates the splitting, while both LSC-ECMA and adQTB yield slightly broader peaks that are centered at the correct quantum frequencies. Interestingly, all LSC-IVR spectra are similar (i.e. IVR0, LGA or ECMA), with splitting values very close to the exact result (thus only ECMA is presented). This is relatively surprising and indicates that position-momentum correlations, although decisive in the overtone and in the combination bands intensities, play only a minor role in determining the Fermi splitting. Indeed, combining classical propagation with a sampling that correctly captures the zero-point energy in each mode is enough for an accurate description of the Fermi resonance. For the same reason, the adQTB also yields accurate results, and the adQTB spectra, shown as green lines in Figure 3, are very similar to those obtained via LSC-IVR. In these calculations, the QTB spectra are deconvoluted using the iterative procedure of Ref. 54 in order to suppress the spectral broadening caused by the friction force and make the splitting more visible (though the peaks are not fully separated in this case and a residual amplitude persists between them).

A more quantitative illustration of the performance of the different methods in predicting the Fermi splitting is provided in the central panel of Figure 3. Our simulations are compared to the perturbative expressions derived by Basire et al. in Ref. 52 for the amplitude of the Fermi splitting in the classical (dashed line) and in the quantum (continuous line) frameworks. As in ref. 52, to compare the analytical trends and the simulations, in the latter, the splitting is measured as the distance between the center of mass of each of the two peaks that constitute the Fermi dyad. The perturbative result is in close agreement with the splitting obtained from the numerical solution of Schrödinger equation (black diamond shapes in Figure 3). The comparison clearly shows that both LSC-ECMA (and more generally LSC-IVR for all three Wigner approximations considered) and adQTB provide an almost exact description of the Fermi resonance, while TRPMD and CMD (not shown in the figure) follow the classical curve and yield a vanishing splitting at low temperature. This inaccuracy of the PI-based methods can be related to the displacement of the spectral intensity towards spurious peaks in (T)RPMD and to the averaging out of the fluctuation modes dynamics in CMD. Interestingly, Benson and Althorpe show in Ref. 29 that Matsubara dynamics allows recovering the correct spectral intensities and therefore yields almost exact predictions for the Fermi splitting.

VII. INFRARED SPECTRA FOR GAS-PHASE WATER

To conclude this paper, we compare the performance of the considered methods for the simulation of the IR spectrum of a water molecule in gas phase. Interactions are modeled with the accurate potential energy surface due to Partridge and Schwenke 55 . This system is of particular interest as coupling between the bending and stretching modes of the molecule naturally gives rise to overtones and combination bands in the vibration spectrum. Furthermore, at room temperature, the ZPE of these high-frequency modes strongly dominates over their thermal energy k B T . This corresponds to a quantum regime in which the various trajectory-based methods perform very differently, according to the analysis in sections IV-VI. This system has also been considered for a similar analysis by Benson et al. in ref. 16 and 29, that provide a comparison of the IR absorption spectra obtained from PI-based methods (TRPMD, CMD and the more recent QCMD), from LSC-IVR (using LGA for the Wigner density sampling) and from a quasi-exact wavefunction-based resolution (using the DVR method 56 ). The authors show that PI approaches do not capture the intensity of overtones and combination bands correctly while LSC-LGA is in much better agreement with the quasi-exact reference (although with broadened spectral lineshapes).

In this section, we extend the set of methods considered by Benson et al by exploring the performance of LSC-ECMA and adQTB. However, contrary to ref. 16 we compute the spectrum using the qTIP4p-f linear dipole moment 57 instead of the non-linear Partridge-Schwenke dipole moment. In the linear case, the real part of the standard dipole-derivative ATCF can be expressed exactly as:

C μ μ = dqdp ρ w (q, p) μ e iL Q t μ ( 45
)
while a similar expression would introduce an approximation for a non-linear dipole moment which would directly impact the LSC-IVR results. Figure 4 (top panel) shows the three main peaks of the IR absorption spectrum simulated using classical MD (dashed grey), TRPMD (solid blue), adQTB (dashed yellow) and LSC-ECMA (dotted dashed green) at 300 K. The spectra are also compared to the numerically exact DVR results (filled curve). Note that the lineshape of the DVR spectrum has been artificially broadened for visualization purposes (by damping the TCF with a Hann window of length 750 fs, consistently with ref. 16). All approximate quantum methods yield similar results in this spectral region and agree reasonably well with the DVR spectra, although none of the methods is able to capture the fine structure of the stretching band (≈ 3700cm -1 ). The bottom panel of Figure 4 shows the anharmonic resonances that appear in the high-frequency region of the spectrum (from 4500cm -1 to 8500cm -1 ). They consist in two main features around 5100-5500cm -1 (sum of the bending and stretching frequencies, which we will label as combination band) and 7000-7400cm -1 (twice the stretching frequency, which we will label as overtone). Consistently with the perturbative analysis, TRPMD strongly underestimates these resonances, with intensities comparable to the classical results (although an accurate estimation of the combination band intensity is difficult due to the strong broadening of the stretching peak). LSC-IVR and adQTB methods, on the other hand, yield more intense peaks, closer to the exact result (though with broader lineshapes) and satisfactory agreement with the perturbative results for combination bands. Indeed, for LSC-LGA, both resonances intensities are overestimated by approximately 40%. The LSC-ECMA0 sligtly underestimates the combination band (∼ -7%) and underestimates more significantly the overtone (∼ -25%). Contrary to the perturbative case, the Edgeworth correction plays an important role in this more anharmonic case: weighting the samples with the first term in the Edgeworth correction (ECMA4) increases the resonance intensities which are both slightly overestimated (∼ 10%). The inclusion of higher orders of the Edgeworth correction, although more computationally demanding, might further improve the results and their effect should be investigated in more details in future work. This impact of the Edgeworth correction is a strong signal that the perturbative approach does not hold quantitatively for this very anharmonic potential.

Despite this strong anharmonicity, the adaptive QTB only slightly underestimates the combination band (∼ -5%, although difficult to estimate accurately due to the broadening of the stretching peak) and is able to accurately capture the overtone intensity (less than 1% discrepancy). For this model, the random force intensity is adapted separately for O and H so that the quantum fluctuation-dissipation relation of eq. ( 13) is enforced for both species (details in appendix D). The accuracy of the results seems to indicate that the quantum fluctuation-dissipation theorem provides a robust criterion to correctly thermalize anharmonic resonances.

VIII. CONCLUSIONS

Several trajectory based methods for approximate calculation of quantum time-correlation functions were compared via analytical and numerical work. These approaches combine approximate sampling of the quantum thermal density with propagation of (generalised) classical trajectories. In our tests, we considered simple low dimensional models endowed with spectral features such as overtones, combination bands and Fermi resonances, that are problematic for trajectory-based methods. The observed failure to correctly reproduce these features, often interpreted as a manifestation of quantum in- terference effects, is typically attributed to lack of coherence in the classical time propagation. Our results contradict this interpretation, pointing to an alternative explanation that highlights the role of accurate initial condition sampling and the subtle interplay of different terms in the classical(like) propagators. In particular:

1. For the weakly anharmonic models of overtones and combination bands considered in this work, the classical evolution of IVR-LSC is exact. The relative accuracy of IVR-LSC approaches then depends critically on the method employed to sample the Wigner thermal density. In particular, approximations of this function that do not take into proper account quantum correlations between positions and momenta (such as IVR0 and LHA) fail to reproduce the correct trends with tem-perature of the anaharmonic features intensities. The performance of the different methods is confirmed numerically also for more complex systems.

2. Matsubara dynamics is also exact for the model overtone.

3. Path integral dynamics methods such as (T)RPMD and CMD follow fully classical behaviour for the temperature dependence of the anharmonic features' intensities. Perturbative analysis based on the Matsubara formalism demonstrates (for the overtone) that this is due to inaccurate description of the coupling between the centroid and fluctuation modes.

The results listed above for path integral based dynamics are confirmed by a recent study 29 employing classical perturbation theory on the action angle representation of the Matsubara dynamics. Furthermore, our work puts on firmer ground previous empirical observations 16 , indicating good accuracy of LSC-IVR for anharmonic spectral features in liquid water. The relative merits of LSC-IVR and path integral dynamics approaches have been examined also in other relevant cases. Problematic features in the vibrational spectra of simple models computed via RPMD (spurious peaks) and CMD (temperature dependence of the peaks' positions) where reported and explained in Ref. 58, and prompted developments, such as TRPMD and the Quasi-CMD 28 to mitigate them. In contrast, LSC-IVR type of methods, while more expensive, were directly able to produce accurate results for these test cases 59 . In the context of vibrational energy relaxation, typically computed from the Fourier transform of the force autocorrelation function (a strongly non-linear observable), it was shown that quantum effects in the high-frequency region of the spectrum can be captured by LSC-IVR, but not by CMD 49 . Interestingly in view of our analysis of the role of the coupling of the centroid mode with the fluctuation Matsubara modes, in ref. 49, it was observed that the centroid correlation function could be obtained from the LSC-IVR one by decoupling the zero-frequency normal mode from the other modes of the ring polymer. The importance of accounting for these higher normal modes in vibrational energy relaxation was also highlighted. This relationship between LSC-IVR and CMD was later formally derived, for time-correlation functions in which one of the two operators is the position 60 . The conclusions of this paper are consistent with these previous results and reinforce them via the adopted perturbative approach. Of course, the pathologies of LSC-IVR stemming from the lack of conservation of the quantum probability density can also lead to dramatic failures of the method, such as ice melting within 1ps of simulation 16 , when the resulting unphysical transfer of energy from high to low frequency modes is faster than the typical time scale of the phenomenon under investigation.

We have also shown that, in contrast with path integral based methods, the (ad)QTB approach gives surprisingly accurate results for the non trivial observables considered in this work, provided that the momentum force correlation function is considered for the calculation of the spectra. The stochastic nature of this approach, that lacks a clear formal justifica-tion both as a sampler of the thermal density and as a generator of time dependent quantities, complicates an analytical investigation of its properties. A qualitative explanation of the remarkable performance of the QTB can be given observing some analogies with IVR based methods. The relevant frequencies for overtones and resonances are typically more than an order of magnitude greater than the damping frequency in the generalized Langevin QTB evolution. In these conditions, the friction and stochastic forces play a negligible role in the short-term dynamics which is therefore essentially classical as for the IVR. The generalized Langevin thermostat only ensures that the configurations that are explored follow (at least approximately) the correct quantum distribution. The results of the exploratory calculations presented here, and complementary work on more realistic systems [START_REF] Mauger | Nuclear quantum effects in liquid water at near classical computational cost using the adaptive quantum thermal bath[END_REF] , indicate that that this method deserves careful consideration as an affordable simulation tool for highly dimensional systems with relevant NQEs.

The analysis described in this paper is based on relatively simple models and further tests are necessary to determine if our results indicate general trends for highly anharmonic condensed phase systems or for more complex (non linear) correlation functions. In spite of this, they present clear evidence that interference effects cannot always "be blamed" for failures of trajectory-based methods. The interplay between sampling and propagation is quite subtle and quantum correlations, rather than coherence, can play a surprisingly important role.

Appendix A: More details on the classical perturbation theory

In this appendix, we provide more details on the classical perturbative calculation of the overtone intensity for the unidimensional potential of equation (23).

The methodology that we use to compute η in the classical case goes as follows: (1) we expand the Boltzmann density ρ in powers of the anharmonicity parameter λ . (2) We write the Dyson expansion of the time propagator e iL t at second order in λ (the order at which the overtone contribution appears, as we will see in the following) and compute the propagated momentum e iL t p. (3) we combine the two previous results by averaging the second-order expansion of pe iL t p on the distribution ρ, which yields c pp (t) at second order in λ . ( 4) Finally, we evaluate the Fourier transform of c pp (t) to obtain the overtone contribution to the VDOS η defined in eq. (24).

The perturbed Boltzmann density, expanded to second order in λ , is given by: ρ(q, p) = ρ 0 (q, p) 1 -λ β q 3 3 + λ 2 β 2 18 q 6q 6 0 +O λ 3 (A1) where ρ 0 (q, p) is the harmonic Boltzmann density and • 0 denotes the expectation value over the Gaussian density ρ 0 .

We use the Dyson expansion of the classical time propagator to rewrite the perturbed propagator e iL t as a series of harmonic propagations (that are represented by the operator e iL 0 t and can be computed analytically) with intermediate applications of the perturbation Liouvillan -λ q 2 ∂ ∂ p at different times. The Dyson series of the classical time propagator truncated at second order is given by: e iL cl t =e iL 0 tλ t 0 ds e iL 0 (t-s) q 2 ∂ ∂ p e iL 0 s + λ 2 t 0 ds 1 e iL 0 (t-s 1 ) q 2 ∂ ∂ p

s 1 0 ds 2 e iL 0 (s 1 -s 2 ) q 2 ∂ ∂ p e iL 0 s 2 + O λ 3 (A2)
The zeroth order of the expansion corresponds to a purely harmonic propagation. The first order involves a single application of the perturbation at a time s combined with harmonic propagations from time 0 to s and s to t. The second order corresponds to two applications of the perturbation at different times. The effect of the harmonic time propagator e iL 0 s on q and p is: e iL 0 s q = q cos(ω 0 s) + p mω 0 sin(ω 0 s) (A3) e iL 0 s p = p cos(ω 0 s)mω 0 q sin(ω 0 s) (A4)

Using (A1) and (A2), the TCF of the perturbed system can be computed entirely from the properties of the harmonic density and time propagator. In the following, we denote as c (mn) pp (t) the contribution obtained from the expansion at order n of the density and order m for the Dyson series. At second order in λ , we thus obtain:

c pp (t) = m+n≤2 ∑ m,n c (mn) pp (t) + O λ 3 (A5)
Note that, since the unperturbed harmonic Liouvillian does not produce overtones, the c (m0) pp terms do not contribute to the intensity of the overtone and they do not appear in the calculation of η. In particular, since c (20) pp is not relevant, the first-order expansion of the density is sufficient to compute η(β ) up to second order in λ . In this section, we nonetheless compute all c (mn) pp terms in order to derive the full second-order expression of the TCF.

The unperturbed contribution for the classical harmonic oscillator reads:

c (00) pp (t) = m β cos(ω 0 t) (A6)
Using the parity and separability of the harmonic Boltzmann density, one can easily show that:

c (01) pp = c (10) pp = c (20) pp = 0 (A7)
The first anharmonic contribution to the correlation function is thus of second order in λ and comes from the two contributions c

pp and c

pp . Using first orders of (A1) and (A2), we obtain:

c (11) pp (t) = λ 2 β
3 dqdp ρ 0 (q, p) q 3 p t 0 ds e iL 0 (t-s) q 2 ∂ ∂ p e iL 0 s p (A8)

After applying the time propagators and exploiting the symmetries of the density, this term reduces to:

c (11) pp (t) = 2λ 2 β 3mω 0 q 4 0 p 2 0 × t 0 ds cos(ω 0 s) cos(ω 0 (t -s)) sin(ω 0 (t -s)) = 2λ 2 3m 2 ω 6 0 β 2 [cos(ω 0 t) -cos(2ω 0 t)] (A9)
Finally, the second order of the Dyson series is:

c (02) pp (t) =λ 2 dqdp ρ 0 (q, p) p t 0 ds 1 e iL 0 (t-s 1 ) q 2 ∂ ∂ p s 1 0
ds 2 e iL 0 (s 1 -s 2 ) q 2 ∂ ∂ p e iL 0 s 2 p (A10)

Applying the propagators and averaging over the harmonic density, we finally obtain:

c (02) pp (t) = 10λ 2 9m 2 ω 6 0 β 2 [cos(2ω 0 t) -cos(ω 0 t) + 3 2 ω 0 t sin(ω 0 t)] (A11)
The full expression of c pp (t) at second order in λ is then:

c pp (t) = m β cos(ω 0 t) + λ 2 m 2 ω 6 0 β 2 4 9
cos(2ω 0 t)cos(ω 0 t)

+ 5 3 ω 0 t sin(ω 0 t) + O λ 3 (A12)
By Fourier transform, we then obtain the overtone contribution to the classical VDOS, as defined in equation ( 24):

η cl (β ) = 4λ 2 9m 3 ω 6 0 β -1 + O λ 3 (A13)
(A11) contains a term proportional to ω 0 t sin(ω 0 t) which diverges when t → ∞. This divergent term, called a secular term, is a well known artefact of this type of Dyson series (see refs. 64 and 65) and can be removed using multiplescale perturbation [START_REF] Bender | [END_REF]67 , resummation techniques 68,69 or by using an appropriate change of variables (for example actionangle variables as in ref. 29). As this term does not affect the overtone intensity (as confirmed by numerical tests), we simply chose to discard it and to discard similar terms appearing in the perturbative calculation for the approximate quantum methods.

Appendix B: Perturbation theory for ring-polymer based methods

In this appendix, we extend the perturbation approach presented in appendix A to ring-polymer methods. We show that Matsubara dynamics yields the exact result, at second order in λ , for the overtone intensity in the unidimensional potential of equation (23). We also show that RPMD and CMD yield the classical result.

Matsubara dynamics

We recall the Matsubara dynamics approximation to the Kubo-transformed momentum ATCF 24 :

K Mats pp (t) = lim M→∞ dQdP ρ M (Q, P) P 0 e iL M t P 0 (B1)
where Q = (Q n ) n∈M and P = (P n ) n∈M are the M first Matsubara modes (with the notation n ∈ M denoting n = -(M -1)/2, . . . , (M -1)/2 ) among which (Q 0 , P 0 ) is the centroid, and ρ M (Q, P) is the Matsubara quasi-density defined as:

ρ M (Q, P) = 1 Z M e -β ∑ n∈M P 2 n 2m +U M (Q)-iθ M (Q,P) (B2)
with the phase factor:

θ M (Q, P) = ∑ n∈M P n ωn Q -n (B3)
and the Matsubara potential:

U M (Q) = 1 2 mω 2 0 ∑ n∈M Q 2 n + λ 3 ∑ i, j,k∈M A i jk Q i Q j Q k (B4)
The A i jk are constants that can be derived for the cubic potential following the procedure described in the supplementary material of ref. 24. The Matsubara partition function Z M normalizes the distribution. The operator iL M is the Matsubara Liouvillian defined, for the potential (23), as:

iL M (Q, P) = ∑ n∈M iL 0 (Q n , P n ) -λ ∑ i, j,k∈M A i jk Q i Q j ∂ ∂ P k (B5)
where iL 0 is the classical harmonic Liouville operator. Note that, for a purely harmonic system, the Matsubara modes evolve independently, each following a classical harmonic motion. As a consequence, coupling between the centroid and the fluctuation modes only arises from the cubic perturbation term. We denote the unperturbed Matsubara distribution as:

ρ M,0 (Q, P) = e iβ θ M (Q,P) ∏ n∈M e -β H 0 (Q n ,P n ) Z n,0 (B6)
with H 0 the unperturbed Matsubara Hamiltonian,

Z n,0 = dQ n dP n e -β (H 0 (Q n ,P n )+ 1 2 m ω2 n Q 2 n
) and ωn = 2nπ/β h. Once the ATCF is written under the form (B1), the same perturbation approach can be applied as for the classical dynamics. First, the Matsubara distribution is expanded as:

ρ M (Q, P) = ρ M,0 (Q, P) 1 - β λ 3 ∑ i, j,k∈M A i jk Q i Q j Q k + O λ 2 (B7)
Then, the time propagator is expanded to second order using the Dyson series of the perturbed Matsubara Liouvillian (B5).

As in the classical case, we decompose the Matsubara time correlation function (B1) as:

K M pp (t) = i+ j≤2 ∑ i, j K M(i j) pp (t) + O λ 3 (B8)
Since the fluctuation Matsubara modes are decoupled from the centroid in the unperturbed dynamics, it can easily be shown that:

K M(00) pp (t) = m β cos(ω 0 t) (B9)
which is the standard result for the harmonic oscillator. Furthermore, using the parity of the harmonic distribution it follows that:

K M(10) pp = K M(01) pp = 0 (B10)
For second order terms, one can show that:

K M(11) pp (t) = 2λ 2 β 3mω 2 0 [cos(ω 0 t) -cos(2ω 0 t)] × dQdP ρ M,0 (Q, P)P 2 0 Q 2 0 ∑ i =0 Q 2 i + Q 4 0 3 (B11)
Since the integrand does not involve P i (i =0) , one can integrate out the momenta and obtain:

K M(11) pp (t) = 2λ 2 3ω 2 0 cos(ω 0 t) -cos(2ω 0 t) × Q 2 0 0 ∑ i =0 Q 2 i 0 + Q 4 0 0 3 (B12)
Since the unperturbed density for Q 0 is Gaussian, we have

Q 4 0 0 = 3 Q 2 0 2 0 therefore K M(11) pp simplifies to: K M(11) pp (t) = 2λ 2 3mω 4 0 β cos(ω 0 t) -cos(2ω 0 t) ∑ n∈M Q 2 n 0
(B13) When the number M of Matsubara modes tends to infinity, the sum appearing in equation (B13) can be related to the average harmonic potential energy:

lim M→∞ ∑ n∈M 1 2 mω 2 0 Q 2 n 0 = 1 2 Θ(ω 0 , β ) (B14) Finally, lim M→∞ K M(11) pp (t) = 2λ 2 Θ(ω 0 , β ) 3m 2 ω 6 0 β [cos(ω 0 t) -cos(2ω 0 t)]
(B15) Using a similar method (but with slightly more involved calculations), it can be shown that the remaining second order term in the expansion of the TCF takes the following form:

lim M→∞ K M(02) pp (t) = λ 2 Θ(ω 0 , β ) m 2 ω 6 0 β 10 9 cos(2ω 0 t)-cos(ω 0 t) + 5 3 ω 0 t sin(ω 0 t) (B16)
Furthermore, it can be shown that K M(02) pp = 0 for the cubic perturbation under study, so that the Matsubara TCF is given by: lim

M→∞ K M pp (t) = m β cos(ω 0 t) + 4λ 2 Θ(ω 0 , β ) 9m 2 ω 6 0 β cos(2ω 0 t) -cos(ω 0 t) + O λ 3 (B17)
where we discarded the secular term that originates from eq. (B16). The contribution of the overtone to the VDOS, approximated using Matsubara dynamics, is then:

η Mats (β ) = 4λ 2 9m 3 ω 6 0 Θ(ω 0 , β ) + O λ 3 = η Q (β ) + O λ 3 (B18)

Ring-polymer molecular dynamics

The procedure for RPMD is almost identical to that for Matsubara dynamics, except that we replace the Matsubara distribution and Liouvillian by the RPMD ones given respectively by:

ρ RPMD (Q, P) = 1 Z RPMD e -β ∑ n∈M P 2 n 2m + 1 2 m ω2 n Q 2 n +U M (Q) (B19) iL RPMD (Q, P) = ∑ n∈M iL 0 (Q n , P n ) -∑ n∈M m ωn Q n ∂ ∂ P n -λ ∑ i, j,k∈M A i jk Q i Q j ∂ ∂ P k (B20)
In the limit of an infinite number of modes, the final result is given by:

K RPMD pp (t) = m β cos(ω 0 t) + 4λ 2 9m 2 ω 6 0 β 2 [cos(2ω 0 t) -cos(ω 0 t)] + 4λ 2 m 2 ω 6 0 β 2 ∞ ∑ n=1 ω 6 0 (ω 2 0 + ω2 n )(3ω 2 0 + 4 ω2 n ) 2 × cos 2 ω 2 0 + ω2 n t -cos(ω 0 t) + O λ 3 (B21)
The overtone intensity at frequency 2ω 0 is therefore identical to the classical intensity. The last line corresponds to the so-called spurious resonances of the RPMD method, appearing at frequencies whose expression combines the physical frequency ω 0 and the ring-polymer internal frequencies ωn .

Centroid molecular dynamics

Finally, for CMD, we first obtain a perturbative expression of the centroid potential as:

V C (Q 0 ) = V (Q 0 ) + λ Q 0 β Θ(ω 0 , β ) -1 β mω 2 0 + λ 2 Q 2 0 α C (ω 0 , β ) β m 2 ω 4 0 + O λ 3 (B22) with α C (ω 0 , β ) = 1 + β Θ(ω 0 , β ) 2 [β Θ(2ω 0 , β ) -2β Θ(ω 0 , β ) -1] ( B23
) and then apply the classical procedure by replacing the standard potential by V C which yields the following momentum ATCF:

K CMD pp (t) = m β cos(ω 0 t) + 4λ 2 9m 2 ω 6 0 β 2 [cos(2ω 0 t) -cos(ω 0 t)] + O λ 3 (B24)
As in RPMD, the only overtone contribution at frequency 2ω 0 is the identical to the classical expression (A12) (Note that in both CMD and RPMD expressions of the ATCF, we discarded the secular terms).

Appendix C: The Edgeworth Conditional Momentum Approximation

In this appendix, we summarize the derivation of the Edgeworth conditional momentum approximation (ECMA). This method makes use of the Edgeworth expansion developed in Ref. 41 to provide an optimal Gaussian approximation for the conditional momentum distribution, together with a correction term that takes the form of an infinite series converging asymptotically to the exact Wigner distribution (including, in particular, its negative parts).

Let us begin by rewriting the Wigner conditional momentum distribution of eq.( 14) as

ρ c w (p|q) = d∆ e i p∆ h q -∆ 2 e -β Ĥ q + ∆ 2 d∆ 2π h δ (∆) q -∆ 2 e -β Ĥ q + ∆ 2 (C1)
with δ the Dirac function. We then insert the relation 2π hδ (∆) = dp e i p∆ h to obtain:

ρ c w (p|q) = d∆ e i p∆ h q -∆ 2 e -β Ĥ q + ∆ 2 dp d∆ e i p∆ h q -∆ 2 e -β Ĥ q + ∆ 2 (C2)
thus making clear the normalization of the conditional momentum distribution. Finally, we trivially rewrite the pseudodensity as:

ρ c w (p|q) = e -κ 2 (q) p 2 2h 2 d∆ e i p∆ h +κ 2 (q) p 2 2h 2 f (∆|q) dp e -κ 2 (q) p 2 2h 2 d∆ e i p∆ h +κ 2 (q) p 2 2h 2 f (∆|q) (C3)
with the probability density

f (∆|q) = q -∆ 2 e -β Ĥ q + ∆ 2 d∆ q -∆ 2 e -β Ĥ q + ∆ 2 (C4)
and κ 2 (q) an arbitrary function of position. Expanding e ip h ∆+ κ 2 p 2 2h 2 as a power series with respect to p, and choosing κ 2 as the second cumulant of the probability density f ,

κ 2 (q) = d∆ ∆ 2 f (∆|q) = ∆ 2 f |q (C5)
yields the so-called Edgeworth Conditional Momentum Approximation:

ρ c ECMAn (p|q) = e -κ 2 (q) p 2 2h 2 C EWn (q, p) dp e -κ 2 (q) p 2 2h 2 C EWn (q, p) (C6) with C EWn (q, p) = 1 + n ∑ m=4,even κ m (q) m! ip h m (C7)
where κ m (q) is the m-th order cumulant of the density f (∆|q), for example, κ 4 (q) = ∆ 4 f |q -3 ∆ 2 2 f |q . In equation (C6), ρ c ECMAn (p|q) is approximated as the (normalized) product of a positive function times the Edgeworth correction factor C EWn (q, p), that takes the form of an expansion in powers of p truncated at order n. The coefficients of this expansion, the cumulants of f , can be computed numerically without facing a sign problem 41 . Note that the odd orders of C EWn (q, p) vanish by parity of f with respect to ∆.

The choice of κ 2 (q) as the second cumulant of f (∆|q) cancels the second order term in the Edgeworth correction series so that the distribution ρ c ECMA0 (p|q) is an optimal Gaussian approximation to the momentum distribution (in the sense of the Edgeworth expansion). In particular, it ensures that ECMA0 is exact when the momentum distribution is actually Gaussian, as for the harmonic potential or in the classical limit 41 .

In order to sample momenta from the distribution (C6), one has to numerically estimate κ 2 (q) (and higher-order cumulants if necessary) as it does not have an analytical expression. This is done by sampling configurations from the pathintegral representation of the density f (∆|q) in an auxiliary PI calculation (see ref. 41 for details). This procedure must be repeated at each configuration q sampled from the marginal Wigner distribution and constitutes most of the computational cost of the method.

The numerical estimation also implies that κ 2 (q) is subject to statistical noise. This noise can induce a bias on the conditional momentum distribution, which in turn can affect the estimation of momentum-dependent observables. In practice, the noise tends to broaden the momentum distribution and leads to an overestimation of the kinetic energy. This effect is corrected by reweighting the distribution with a factor depending on the variance σ 2 of the numerical estimate for κ 2 . For the ECMA0 momentum distribution, it can be shown that the reweighting factor is given, at first order in σ , by:

C EW0 noise (p ; κ 2 , σ ) = 1 -σ 2 p 4 8h 4 - p 2 4h 2 κ 2 - 1 8κ 2 2 (C8)
In practice, σ 2 is estimated by computing the empirical variance of κ 2 (q) from multiple independent calculations (i.e. multiple auxiliary PI simulations of (C4)). A similar correction can be derived for multidimensional system.

In this appendix, we derive a perturbative expression for the matrix elements of the density operator w(q, ∆; β ) = q -∆ 2 e -β Ĥ q + ∆ 2 for the potential (23). Such matrix elements appear both in the definition of the Wigner density and in its ECMA approximation. We follow ref. 42 (section 12.4) to express w(q, ∆; β ) in its continuous path integral form. We then decompose the fluctuations around the classical harmonic path in Fourier components and we obtain analytical expressions for w(q, ∆; β ) at first order in λ .

For the potential (23), the continous path integral form of the off-diagonal density matrix elements are given by the following functional integral: w(q, ∆; β ) = with S 0 [x] the harmonic imaginary-time action. With this expansion in λ , the average over anharmonic paths is translated into the much simpler problem of averaging the cubic potential term on harmonic paths. To that end, we expand the path x(τ) about the classical harmonic path 42 :

x(τ) = x 0 (τ) + y(τ) (E4) where x 0 (τ) is the classical harmonic path, solution of the differential equation m ẍ0 = mω 2 0 x 0 with boundary conditions x 0 (0) = q -∆/2 and x 0 (β h) = q + ∆/2. It can be written as:

x 0 (τ) = (q - ∆ 2
)(e -ω 0 (τ-β h)e ω 0 (τ-β h) )

+ (q + ∆ 2
)(e ω 0 τe -ω 0 τ ) × 1 e β hω 0e -β hω 0 (E5)

The fluctuation path y(τ) (defined as the difference between the full path and the classical harmonic path) has the property y(0) = y(β h) = 0. One can show (see ref. Moreover, expanding y(τ) in a Fourier series and rewriting the functional integral in terms of its Fourier coefficients as in ref. [START_REF] Tuckerman | Statistical Mechanics: Theory and Molecular Simulation[END_REF] , it can be shown that the term proportional to x 0 (τ) does not depend on ∆, so that we can write:

w(q, ∆; β ) = w 0 (q, ∆; β ) 1 -λ 3h [a(q, ∆) + b(q)] + O λ 2 (E8) where w 0 is the harmonic off-diagonal element given by w 0 (q, ∆; β ) = I 0 e Analytical expressions of the functions a(q, ∆) and b(q) can be obtained using the Fourier decomposition of the fluctuation path:

a(q, ∆) = qh 2Θ(ω 0 , β ) The marginal position density defined as ρ m w (q) = w(q, 0; β )/ dq w(q, 0; β ) is then given by:

ρ m w (q) = e - mω 2 0 q 2 2Θ(ω 0 ,β ) 2πΘ(ω 0 , β )/mω 2 0 1- 2λ q 3mω 2 0 Θ(2ω 0 , β ) Θ(ω 0 , β ) -1 - 2λ q 3 9Θ(ω 0 , β ) 5 2 - Θ(2ω 0 , β ) Θ(ω 0 , β ) + O λ 2 (E14)
We also obtain the exact conditional momentum distribution (via the Fourier transform of w(q, ∆; β ) with respect to ∆) :

ρ c w (p|q) = ρ 0 (p) 1 + λ q 3mω 2 0 Θ(2ω 0 , β ) Θ(ω 0 , β ) -1 × p 2 mΘ(ω 0 , β ) -1 + O λ 2 (E15)
To obtain the ECMA0 expression, we compute the variance of the density f (∆|q) = d∆ ∆ 2 w(q,∆;β ) d∆ w(q,∆;β ) at first order in λ :

κ 2 (q) = h2 mΘ(ω 0 , β ) 1 - 2λ q 3mω 2 0 Θ(2ω 0 , β ) Θ(ω 0 , β ) -1 + O λ 2 (E16)
which yields (through eq. (C6)) the exact conditional momentum distribution (E15) at this order in perturbation. Finally, the successive cumulants of ρ c can be computed to show that the Edgeworth correction does not modify the ECMA distribution at first order. The fourth order moment of f reads:

∆ 4 f |q = 3 h2 mΘ(ω 0 , β ) 2 × 1 - 4λ q 3mω 2 0 Θ(2ω 0 , β ) Θ(ω 0 , β ) -1 + O λ 2 (E17)
so that we obtain for the fourth order cumulant κ 4 (q) = ∆ 4 f |q -3κ 2 2 (q) = O λ 2 . Furthermore, assuming that, for any even integer m, all cumulants κ j for 4 ≤ j ≤ m -2 are null at first order in λ , we can use the recursion formula for the order-m cumulant:

κ m (q) = ∆ m f |q - m-2 ∑ j=2, even m -1 j -1 κ j (q) ∆ m-j f |q = ∆ m f |q -(m -1)κ 2 (q) ∆ m-2 f |q (E18)
From this expression, using similar reasoning as above for κ 4 , one can show that κ m (q) = O λ 2 . Therefore, by the recursion principle, all cumulants of ρ c of order higher than two are null at first order in λ and thus the Edgeworth correction does not contribute to the first-order expansion of the ECMA distribution.

Appendix F: Full perturbative expression of the combination bands intensities with the Local Harmonic Approximation

The contributions of the combination bands to the LSC-LHA vibrational density of states are:

η LHA (β ; ω 1 + ω 2 ) = η Q (β ; ω 1 + ω 2 ) γ LHA ω 1 +ω 2 η LHA (β ; ω 1 -ω 2 ) = η Q (β ; ω 1 -ω 2 ) γ LHA ω 1 -ω 2 η LHA (β ; 2ω 2 ) = η Q (β ; 2ω 2 ) γ LHA

FIG. 1 :

 1 FIG.1: Contribution η of the overtone to the VDOS as function of temperature for the potential (23) computed using different methods. Lines are perturbative expressions while bullet points are simulation results. In all the calculations, m is the proton mass, ω 0 = 1500 cm -1 and λ = 21 kcal/mol/ 3 , which is small enough for the perturbative analysis to be highly accurate.

FIG. 2 :

 2 FIG. 2: Contributions of the three combination bands to the VDOS for the potential (34) as a function of temperature for the different methods. Lines are analytical expressions and crosses are simulation results.

FIG. 3 :

 3 FIG.3: Left and right panels: Kubo-transformed momentum autocorrelation spectrum at 150K (left) and 800K (right) in the region of the Fermi dyad for TRPMD, LSC-ECMA and adQTB, compared to the exact results. The exact spectra are obtained by solving the Schrödinger equation on a basis of harmonic wavefunctions and are artificially broadened by convolution with a Gaussian function so that the height of the peaks approximately match the other methods. Center panel: Fermi splitting as a function of temperature approximated with TRPMD, LSC-ECMA and adQTB. This splitting is measured as the distance between the center of mass of each of the two peaks that constitute the Fermi dyad. Arbitrary units are used for the intensities of the spectra.

FIG. 4 :

 4 FIG.4: Infrared spectrum of gas-phase water (Partridge-Schwenke potential energy and qTIP4P-f linear dipole moment) at 300K. Arbitrary units are used for the peak intensities.

  x(β h)=q+∆/2x(0)=q-∆/2 Dx(τ) exp -S λ [x] h (E1)with the imaginary-time actionS λ [x] =) can be expanded to first order in λ as:w(q, ∆; β ) = x(β h)=q+∆/2 x(0)=q-∆/2Dx(τ)e -S 0 [x]

  42) that the harmonic action is separable between the classical and fluctu-ation paths i.e. e -S 0 [x] h = e -S 0 [x 0 ] h e -S 0 [y] h .We thus obtain:w(q, ∆; β ) = e -S 0 [x 0 ] h Dy(τ)e -S 0 [y] h × 1 -λ 3h β h 0 dτ (x 0 (τ) + y(τ)) 3 + O λ 2 (E6)where the closed integral indicates a sum over paths y(τ) starting and ending at 0. Due to the symmetry of the fluctuation paths the terms containing odd powers of y(τ) cancel when integrated over all paths and we obtain:w(q, ∆; β ) = e -S 0 [x 0 ]h Dy(τ)e -S 0 [y] ) + 3x 0 (τ)y 2 (τ) + O λ 2 (E7)

  mω 4h tanh(β hω/2) ∆ 2 -mω h tanh(β hω/2)q 2 (E9)with the constant I 0 = Dy(τ)e -S 0 [y] h and the functions a and b are defined as: a(q, ∆) =

  0 , β ) Θ(ω 0 , β ) -1 (E13) with Θ(ω 0 , β ) = hω 0 /2 tanh(β hω 0 /2) .

ACKNOWLEDGEMENTS

We thank Raz Benson and Stuart Althorpe for interesting discussions and for sharing with us the data for the DVR correlation functions of gas phase water and the contents of ref. 29 before publication.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In this Appendix, we illustrate the effect of the adaptation process of the adQTB method on the intensity of anharmonic resonances with the particular case of the H 2 O molecule examined in section VII. The top panel of figure 5 compares the VDOS obtained from the two different estimators defined in section II D and IV C. In the standard QTB method, the two estimators strongly differ with the QTB-a spectrum (obtained from C pp (ω)) being less intense, in agreement with the results obtained for the prototype perturbative models. In the adQTB simulations, on the other hand, the spectra obtained from the two estimators are essentially identical and coincide with the standard QTB-b result. The converged γ r (ω) coefficients obtained at the end of the adaptation process are shown in the bottom panel of figure 5). They display sharp peaks at the anharmonic resonance frequencies. The increase of γ r (ω) in this region causes an increase of the adQTB-a spectrum, that eventually reaches coincidence with adQTB-b and therefore a much better agreement with the numerically exact results obtained from the DVR calculations (see section VII). Note that the anharmonic resonances in the VDOS appear on top of a relatively large background that corresponds to the high-frequency tail of the O-H stretching peak. The deconvolution procedure of Ref.54 allows to essentially suppress this background and recover almost unaffected spectral features as those presented in Fig. 4.

Appendix G: List of abbreviations