Chapter 1

Computer simulation of time-dependent quantities for quantum systems is limited by the exponential scaling of exact methods with the number of degrees of freedom. Indeed, while algorithms for classical atomic dynamics, such as molecular dynamics based on empirical potentials, allow us to describe systems with thousands of degrees of freedom for times of the order of several nanoseconds, quantum dynamical simulations (in which atoms are usually represented through wave functions) are limited to 5-10 atoms and to short times (of the order of the picoseconds at most). On the other hand, the calculation of quantum time-dependent quantities is key in many interesting problems and in particular to compute, via linear response theory and time correlation functions, the system's response to external perturbations. That is why the development of approximate quantum dynamical methods is a very active and challenging area of research in Theoretical Chemistry and Physics. One of the main directions currently pursued is to adapt efficient techniques for classical simulations, based on molecular dynamics (MD) and/or Monte Carlo, to the quantum case. A reasonable compromise among accuracy and computational cost is realised by the so-called quasi-classical methods [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Hernandez | Quantum time correlation functions and classical coherence[END_REF]3,[START_REF] Sun | Semiclassical theory of electronically nonadiabatic dynamics: Results of a linearized approximation to the initial value representation[END_REF][START_REF] Shi | Semiclassical Theory of Vibrational Energy Relaxation in the Condensed Phase[END_REF][START_REF] Serena Causo | An adiabatic linearized path integral approach for quantum time correlation functions: electronic transport in metal-molten salt solutions[END_REF] for computing time correlation functions. In these methods, sampling of the exact quantum thermal equilibrium density for the system is combined with generalised trajectories that provide an approximation for quantum dynamics, usually valid only for short times. The main goal of my thesis is to develop, test and apply approximate quasi-classical methods for vibrational spectroscopy. The focus is, in particular on the so called Phase Integration Method (PIM) [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF]9]. In the first formulation of PIM the correlation function, calculated via path integrals, was written as the average, with respect to a probability density known analytically, of an observable containing a phase factor. To obtain that result, we calculate the correlation function via path integrals, combine a change of variable and a linearisation approximation for time propagators which allows us to have the exact initial density with classical trajectories and a phase factor coming from the linearisation. This method was introduced recently and it can be shown that it is the lowest order of a controlled approximation that can, in principle, systematically improve, by including higher order quantum corrections to classical dynamics, the calculation of quantum time correlation functions. To tame the phase factor in PIM, which affects convergence for high dimesional system, a cumulant expansion was introduced. The resulting form of the correlation function does not have a phase factor but the probability density becomes a quantity which can be estimated only numerically and which is affected by statistical uncertainty. As a consequence, we use an original algorithm that combines methods previously developed for sampling noisy distributions via Monte Carlo schemes using Penalty [START_REF] Ceperley | The penalty method for random walks with uncertain energies[END_REF] and Kennedy [START_REF] Kennedy | Noise without Noise: A New Monte Carlo Method[END_REF] methods [START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. As other quasi-classical methods [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF], PIM is based on combining MD and MC algorithms to compute averages of appropriate quantities, and has been applied to obtain, for example, the dynamical quantum structure factor of a relatively high dimensional model of a realistic condensed phase system (64 Ne atoms modelled via the Aziz potential) [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. Comparison with experimental data shows that PIM can capture non trivial quantum effects with an efficiency comparable to/or better than most of other approximate numerical schemes. In spite of these encouraging results, progress must be made to further reduce the numerical cost of PIM based calculations and to improve the flexibility and theoretical soundness of the method. In particular, in this thesis, I focused on adapting and extending PIM to vibrational spectroscopy (a direct probe of the atomic dynamics). This is interesting both from a theoretical and an applicative point of view. The availability of a general and efficient algorithm would in fact allow to study and characterize, in particular, systems with hydrogen bonds which are known to be very sensitive to the quantum nature of the dynamic and for which methods based on completely classical calculations fail. The thesis is organized as follows. Chapter 2 is devoted to a general description of the path integral methods. This formalism is key for PIM and all other schemes considered in this work for comparison purposes, so the principle of the path integrals of Feynman [START_REF] Feynman | Space-time approach to non-relativistic quantum mechanics[END_REF] is presented. Then, the current approximate methods based on path integrals to compute quantum time correlation functions are explained. The two most popular methods to calculate these quantities are Centroid Molecular Dynamics (CMD) and Ring Polymer Molecular Dynamics (RPMD). The details of these methods with their ad hoc approximations are presented. Finally, the principle of quasi-classical methods, such as our method PIM, is detailed with the linearisation approximation. Then two other quasi-classical methods are explained: the FK-LPI method introduced by Poulsen et al. [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF] based on the Feynman-Kleinert variational principle [START_REF] Feynman | Effective classical partition functions[END_REF] and the LGA method introduced by Liu and Miller which use a local Gaussian approximation [START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF]. This will complete the review of the state-of-the-art methods that are more closely related to this work and with whom we compare directly PIM's performance. Chapter 3 illustrates how PIM can be adapted to the calculation of the Wigner density. This is a key quantity for quasi-classical methods since it appears naturally in the linearisation approximation. Unfortunetely, it is not possible to interpret it directly as a probability density or to computre it directly from brut force for a general system due to the presence of a phase factor. Approximate methods (FK-LPI and LGA) have only ad hoc schemes to obtain Wigner density such as harmonic approximation (see LGA). What we do with PIM is better because ,as I will explain, we have a consistent way to treat the phase factor in the Wigner density definition and so we can in principle calculate to Wigner dentity for any system. The principle of the PIM's algorithm is also detailed and results are presented and compared to other methods [START_REF] Basire | Computing Wigner distributions and time correlation functions using the quantum thermal bath method: application to proton transfer spectroscopy[END_REF]17]. Chapter 4 focuses on the adaptation of PIM for the infrared spectroscopy. We explained how and why we have to use the dipole-derivative Kubo autocorrelation function to calculate Infrared spectra [9,18]. The results obtained on models (OH and CH 4 ) will be compared to CMD, RPMD. OH spectra with an 3D shifted highly anharmonic Morse potential is calculated and a comparison of this test is performed with LGA [START_REF] Liu | Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initiobased flexible and polarizable force field[END_REF]. This is an interesting model case because exact results are known and classical simulations are far from this exact result. As a consequence it is the perfect model (small size and physically interesting) for investigating the accuracy of PIM for Infrared spectra. Chapter 5 presents the methodology used to calculate rate constants with PIM. This is then applied to the chemical reaction rate for a linear model of the H+H 2 reaction using an Eckart potential [START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF]21,22,[START_REF] Stuart | Derivation of a true quantum transitionstate theory[END_REF]. Results are compared to classical and fully quantum simulations [START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF]22].

Chapter 2

State of the art

It is a curious historical fact that modern quantum mechanics began with two quite different mathematical formulations: the differential equation of Schroedinger, and the matrix algebra of Heisenberg. The two apparently dissimilar approaches, were proved to be mathematically equivalent. [• • • ] a third formulation of non-relativistic quantum theory [• • • ] was suggested by some of Dirac's remarks concerning the relation of classical action to quantum mechanics. A probability amplitude is associated with an entire motion of a particle as a function of time, rather than simply with a position of the particle at a particular time. This formulation is mathematically equivalent to the more usual formulations. There are, therefore, no fundamentally new results. However, there is a pleasure in recognizing old things from a new point of view. R.P. Feynman

Formalism of path integrals

In 1948, Richard Phillip Feynman presented a new formulation of quantum mechanics [START_REF] Feynman | Space-time approach to non-relativistic quantum mechanics[END_REF]. At that time this formulation was conceived more as a pleasant intellectual exercise than as a truly practical new stream in quantum theory. Nevertheless, the circumstances were bound to change in a short time. As computers hove in sight, Feynman's approach became applicable, at least for the calculation of time independent statistical properties as we shall see in the following, and this paved the way to a number of studies in Physics and Chemistry. The "strangeness" of quantum mechanics is still present in Feynman's formulation as it is in the Heisenberg and Schrödinger pictures. To introduce Feynman's representation, we can consider a particle prepared in a state initially localized in x that evolves unobserved to a point x . The wave packet of the initial sate evolves and spreads in time, causing the state to become more and more delocalized spatially until it is measured at the point x . The measurement is responsible for the collapse of the wave function in x and so the state is localized again. Feynman's picture formulate this similar to a classicle particle view in which the particle is evolving unobserved between x and x . However there is a fundamental difference from the classical point of view. Indeed, classically even if we do not observe the particle and therefore we do not know what path it will take, we know that the particle is going to follow an unique definite path. On the other hand in the Feynman point of view, the quantum nature of the particle makes it impossible to specify a path. Instead of following a unique path, the particle goes between x and x using an infinity of paths.

In the following, we formalize this intuitive picture focusing on the path integral representation of the canonical quantum density. This will be one of the key quantities in our calculation of spectra. Also the developments of the next section can be generalized easily (via the so called Wick rotation [START_REF] Mark | Statistical Mechanics: Theory and Molecular Simulation[END_REF]) to the path integral expression of the quantum real time propagator, the starting point of our approximation to the dynamics.

Path integrals and the canonical density matrix

In this section, path integral will be illustrated on the canonical density matrix.

We will give a more precise mathematical formulation and for simplicity of the demonstration we will focus on a single particle moving in one dimension. The generalisation to three dimension and many particles can be done easily.

In the canonical ensemble a state is defined by the fixed quantities N (number of particles), volume V , temperature T . The total energy of the system is denoted E.

The Hamiltonian of the single particle is:

H = p 2 2m + V ( x) = K + V (2.1)
where V ( x) is the potential, p the momentum and m the mass of the particle.

The definition of the quantum density operator is:

ρ(β) = 1 Z(N, V, T ) e -β H (2.2)
where β = 1 k B T and Z(N, V, T ) (denoted Z in the following) is the canonical partition function and has the following definition:

Z = T r e -β H (2.3)
where T r[] represents the trace of the operator inside the squares.

If we consider a coordinate space matrix elements of the density matrix ρ(β):

ρ(x, x ; β) = 1 Z x |e -β H |x (2.4)
We can notice that H is the sum of two operators K( p) and V ( x) which do not commute with each other. As a consequence:

e -β( K+ V ) = e -β K e -β V (2.5)

Chapter 2. State of the art

To avoid this problem we have to find a way to split kinetic and potential operators. We can do that in two steps. The first one is a "time" composition property:

e -β H = e -β P H P (2.6)
The second is a Trotter decomposition and exploits the fact that if P is an integer suffienciently high, we can write: e -β( K+ V )/P ≈ e -β K/P e -β V /P (2.7)

e -β( K+ V )/P ≈ e -β V /2P e -β K/P e -β V /2P

(2.8)

The first expression is an approximation of the second order in 1/P and the second one of the third order in 1/P . Then we can use the Trotter theorem (1958) which allows us to express our density operator as (considering only the second approximation but analogous with the first one): e -β( K+ V ) = lim P →∞ e -β V /2P e -β K/P e -β V /2P P (2.9)

The derivation of the Trotter theorem is technical and will not be shown in this thesis. However, there are several conditions to satisfy. Our operators have to be lower bounded and the operators T, V and T+V must be self-adjoint. In all cases considered, these conditions are satisfied.

Coming back to our density matrix element and inserting the Trotter theorem, we have:

ρ(x, x ; β) = 1 Z lim P →∞
x | e -β V /2P e -β K/P e -β V /2P P |x (2.10)

In order to simplify the notation we can define an operator γ as:

γ = e -β V /2P e -β K/P e -β V /2P (2.11)

ρ(x, x ; β) = 1 Z lim P →∞
x | γ P |x (2.12)

We then insert P -1 times the identity operator between the P factors of γ. This will introduce P -1 integrations over coordinates label:

I = dy |y y| (2.13) ρ(x, x ; β) = 1 Z lim P →∞ dx 2 • • • dx P x | γ|x P x P | γ|x P -1 • • • x 2 | γ|x (2.14)
The advantage of using the Trotter theorem is that we have to evaluate P -1 identical terms (only the index is different) simpler than the general one. V is only a function of the coordinates and as a consequence is diagonal in the coordinate basis:

x k+1 | γ|x k = x k+1 |e -β V /2P e -β K/P e -β V /2P |x k

x k+1 | γ|x k = e -βV (x k+1 )/2P x k+1 |e -β K/P |x k e -βV (x k )/2P (2.15)

Since K is a momentum operator, the matrix element e -β K/P |x k is slightly less trivial to evaluate. We have to insert the identity operator but in the momentum basis this time and use the fact that we known the form of x|p :

I = dp |p p| (2.16) x|p = 1 √ 2π e i px (2.17)
We insert the identity operator to obtain:

x k+1 |e -β K/P |x k = dp x k+1 |e -β K/P |p p|x k (2.18)

The utility of introducing the identity in the momentum basis becomes clear because now we have e -β K/P which is diagonal in the momentum basis and as a consequence we have:

x k+1 |e -β K/P |x k = dp x k+1 |p p|x k e -βp 2 /2mP (2. [START_REF] Liu | Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initiobased flexible and polarizable force field[END_REF])

x k+1 |e -β K/P |x k = 1 2π
dp e ip(x k+1 -x k )/ e -βp 2 /2mP (2.20)

Then it is simply a Gaussian integral because we integrate p between -∞ and +∞. Thus:

x k+1 |e -β K/P |x k = mP 2πβ 2 1/2 exp - mP 2β 2 (x k+1 -x k ) 2 (2.21)
which leads us to:

x k+1 | γ|x k = mP 2πβ 2
P k=1 mP 2β (x k+1 -x k ) 2 + β 2P (V (x k+1 ) + V (x k )) (2.24)
In equation 2.23, the endpoints of the paths are at x 1 and x P +1 and are fixed at the initiation and detection points so x and x respectively. We can see that the kinetic energy is now present in the form of a harmonic coupling between points (x 1 , • • • , x P +1 ) of consecutive indexes. We are here in the limit P → ∞ of a discretized representation for the density matrix and the intermediate integrations over

x 2 • • • x P represent the sum over all the possible paths between x and x . For finite P , which will be our case for simulations, the paths are lines between successive points as it is shown in Figure 2.1. The probability to follow a path depends of the value inside the exponential in the equation 2. [START_REF] Stuart | Derivation of a true quantum transitionstate theory[END_REF].

We can connect the canonical density matrix ρ to the quantum time propagator U . Indeed, we have:

ρ(β) = U (-iβ ) U (t) = ρ(it/ ) (2.25)
This equivalence allows us to have easily a formal expression of the real time propagotor using path integrals. We can also interpret obtain the density matrix as the time propagator at an imaginary time t = -iβ and we can refer to the density matrix as an imaginary time propagator. Using this correspondance, we have a path integral expression for the coordinate-space matrix elements of the real time propagator:

U (x
exp i P k=1 mP 2t (x k+1 -x k ) 2 - t 2P (V (x k+1 ) + V (x k )) (2.27)
Following the steps to obtain the equation 2.23 which gives us the density matrix elements we can derive the expression of the canonical partition function which is necessary to calculate all the equilibrium properties in the canonical ensemble. Indeed, we have the following expression for our partition function:

Z = T r exp(-β H) (2.28) Z = dx x|e -β H |x (2.29)
We need only to evaluate the diagonal elements of the density matrix in order to obtain the canonical partition function. As a consequence we only need to set x 1 = x P +1 = x in the equation 2.23 and then we integrate all these diagonal elements.

Z = lim

P →∞ mP 2πβ 2 P/2 dx 1 • • • dx P exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β 2P (V (x k+1 ) + V (x k )) Z = lim P →∞ mP 2πβ 2 P/2 dx 1 • • • dx P exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β P V (x k ) (2.30) with x P +1 = x 1 .
The integration over cyclic paths in equation 2.30 is illustrated in Figure 2.2. We can notice that in the classical limit (T → ∞), the harmonic spring constant connecting neighbor points becomes infinite and as a consequence all the points collapse in one point which corresponds to the classicle point particle. Thus, the high temperature limit of the path integral formalism is equivalent to the classical limit. Analytical evaluation of path integrals is only possible for free particle (V (x) = 0) and quadratic potentials (V (x) = 1 2 kx 2 ). However, the path integral picture is applicable for quantum statistical mechanical calculations, even for large system where wave function method are impossible to use due to the exponential scaling with the number of degree of freedom of the numerical cost. In order to do that, we use what we call the classical isomorphism of path integral which allows us to use what we call Path Integral Molecular Dynamics (PIMD) or Path Integral Monte Carlo (PIMC). Unfortunately, this simplification to discrete paths is not so easy with the time propagator due to the complex exponentials (see equation 2.26). The latter causes numerical calculations to oscillate widly as different paths are sampled leading to a severe convergence problem known as the dynamical sign problem. Thus, as we are going to see later in this chapter with the presentation of RPMD, CMD and the linearization approximation, the calculation of dynamical properties from path integrals remains one of the most challenging problem in computational science.

Polymer isomorphism and PIMD

As mentioned before, in order to calculate, via computer simulation, time independent equilibrium statistical properties using path integrals, we need to introduce the idea of the "classical isomorphism" [START_REF] Sprik | Staging: A sampling technique for the Monte Carlo evaluation of path integrals[END_REF]. We are still considering for the moment a single particle moving in one dimensional potential V ( x). We recall the expression of the discrete partition function (see equation 2.30), which is a good approximation of the quantum canonical partition function if P is large enough:

Z P = mP 2πβ 2 P/2 dx 1 • • • dx P exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β P V (x k ) (2.31)
with x P +1 = x 1 and Z P denotes that we are at P large but finite. Equation 2.31 is the classical configurational canonical partition function of a cyclic polymer chain moving in the potential V (x) but at an effective temperature of P T in phase space variables (x,p). The classical analogy can be made more complete by introducing a set of momenta via a Gaussian integration. Then we recast the P Gaussian prefactor as a set of Gaussian integrals over the variable p 1 , • • • , p P so that they look like momenta conjugate to x 1 , • • • , x P :

mP 2πβ 2 = 1 2π dp k e -β P p 2 k 2m
(2.32)

Z P = 1 2π P dp 1 • • • dp P dx 1 • • • dx P × exp - P k=1 β 2mP p 2 k + mP 2β 2 (x k+1 -x k ) 2 + β P V (x k ) = 1 2π P dp 1 • • • dp P dx 1 • • • dx P × exp - β P P k=1 1 2m p 2 k + 1 2 mω 2 P (x k+1 -x k ) 2 + V (x k ) (2.33)
with x P +1 = x 1 and ω P = P β which is called the chain frequency. In the exponential, we have replaced the prefactor by Gaussian integrals. In these Gaussian integrals, the physical mass is appearing. However, since the prefactor does not affect any equilibrium averages we are free to choose the mass m as we like, so not necessarly the physical mass of the system. As explained before, once we write the partition function in this way, we make appear the classical Hamiltonian of a cyclic chain polymer at an effective temperature of P T , a chain frequency of ω P in the potential V (x):

H P = P k=1 1 2m p 2 k + 1 2 mω 2 P (x k+1 -x k ) 2 + V (x k ) (2.34) Z P (L, T ) = 1 2π P dp 1 • • • dp P dx 1 • • • dx P e -β P H P (2.35)
The identificaton of the partition function, in equation 2.32 with the one of a cyclic polymer of P points, leads to the so called classical isomorphism. We are going to exploit the isomorphism between the classical and approximate (because P is finite) quantum partition function because it presents many numerical advantages. The classical isomorphism is illustrated in Figure 2.3. Because the cyclic polymer looks like a necklace, we call very often the P points "beads". The coils represent the fact that neighbor beads are connected via spring characterized by ω P . According to the classical isomorphism, we can use classicle techniques such as Molecular Dynamics (MD) or Monte Carlo (MC) to compute approximate quantum properties. This is not the only advantage of the polymer isomorphism. Indeed, this isomorphism is very useful to visualize and deal with quantum systems. Most of the time the quantum reality goes against our natural intuition and path integrals via the classical isomorphism can help us to visualize it. The idea is to replace one quantum atom by a necklace of many classicle particles. If we look closely to the equation 2.32, we notice that ω P = P β . As a consequence, ω P increases linearly (at P fixed) with the temperature which means that the springs of the polymer become stiffer with the temperature and so at high temperature the P beads are almost not moving relative to one another. In this case the polymer evolution can be reduced to one point and we are in the classical limit.

This argument can be made more precise by considering the root mean square radius of gyration ∆ of our ring polymer as:

∆ 2 = 1 P P k=1 x 2 k f (2.36)
Where f will be defined in the equation 2.47 and corresponds to a probability density. Focusing on a free particle (V (x) = 0), we obtain [START_REF] Gillan | Computer Modelling of Fluids Polymers and Solids[END_REF]: Λ is the de Broglie thermal wave lenght. Thus the polymer collapses in one point if we increase the temperature. On the contrary, at low temperature the polymer is spreading more and more when the temperature decreases.

∆ = Λ(T ) √ 8π (2.
In the presence of a potential the spreading is limited by the potential itself. A polymer ring within a real system simulation experiences two opposing forces. From one side, the beads within an atom want to approach each other, since there are springs connecting them but, on the other hand, there is an inter atomic potential attracting beads from different atoms and spreading the necklaces. The thermal equilibrium is reached when the mean value of these two forces balances. However, should the number of beads per atom be not sufficiently high, the polymer would not spread as much as required for the given conditions of mass, temperature and potential. The results obtained in that situation would be biased towards the classical limit. But, if we put a sufficient amount of beads the polymer would not spread out more than necessary for the given conditions, since the length of the necklace is limited. Thus, increasing the number of beads once the sufficient amount has been found would not lead to better results. In equation 2.30, we said that the partition function becomes exact using an infinite number of beads. However, our results will be exact, provided that a sufficient number of beads is taken, and we only have to worry about the statistical error. From a physical point of view, a sufficient number of beads is roughly reached when the root mean squared length of the springs is smaller than the relevant length scale of the external potential. This translates into the following relation [START_REF] Gillan | Computer Modelling of Fluids Polymers and Solids[END_REF][START_REF] Brualla I Barbera | Path integral Monte Carlo Algorithms and applications to quantum fluids[END_REF]:

P β 2 mσ 2 (2.39)
with P the number of beads, and σ the length scale of the potential. To a first order, we can say that the number of beads required is inversely proportional to the temperature. In practice, the value of P is determined by verifying numerical convergence of an appropriate set of observables (the energy for example).

Up to now, we only presented the path integral formalism and the polymer isomorphism for the special case of one particle in one dimension potential. If we consider N particles in 3 dimensions with interaction between particles, we obtain for the polymer Hamiltonian (following the same procedure as the one dimension case) the following expression:

H P = P k=1 N i=1 1 2m i p (k)2 i + 1 2 m i ω 2 P (r (k+1) i -r (k) i ) 2 + V (r (k) 1 , • • • , r (k) N ) (2.40)
with the condition r

(P +1) i = r (1)
i . An important point to notice is that the potential V (r

(k) 1 , • • • , r (k)
N ) only acts between beads of the same index k. The consequence is that all beads with the same index interact with each other but do not interact with the others. This behaviour is illustrated in Figure 2.5 for the case of two quantum particles. This classical isomorphism is the basis of what we call the Path Integral Molecular Dynamics (PIMD) (or Path Integral Monte Carlo: PIMC) method. For example in PIMD, we let the polymer previously described evolving with the following Newtonian equations (in a one dimensional case but the generalization is staigthforward):

ẋk = ∂H P ∂p k ṗk = - ∂H P ∂x k (2.41) ẋk = p k m ṗk = - mP 2 β 2 2 (2x k -x k+1 -x k-1 ) - ∂V ∂x k (2.42)
where m is not necessarly the physical mass.

If the equations 2.41 are coupled to a thermostat (Nosé-Hoover for example), then the dynamics will sample the canonical distribution. Therefore, these classical trajectories are used as a sampling device to explore the configuration space of the ring polymer and calculate the exact thermodynamic and structural properties of the system in the limit of a sufficiently large number of beads.

Time-independent equilibrium properties in the canonical ensemble from path integral formulation

Now that we know how to calculate the density matrix and the partition function, we are able to derive the expressions for equilibrium properties in the canonical ensemble via the path integral formalism. Furthermore, with the classical isomorphism and PIMD (or PIMC) the methods used to compute these quantites are established. The basic relation for the expectation value of an Hermitian operator, A, in the canonical ensemble is:

A = 1 Z T r A exp(-β H) A = 1 Z dx x| Ae -β H |x (2.43)
In the following I will discuss operators with a form that will be interesting in the calculations presented later in the thesis (see Chapters 3 and 4). The first and also the simplest case to consider is the evaluation of an operator A function of x only. In this case, A is diagonal in the coordinate basis, A( x)|x = A(x)|x , and we can express the expectation value of A as:

A = 1 Z dx A(x) x|e -β H |x (2.44)
we are here in a similar case as the partition function because we only need the Chapter 2. State of the art diagonal elements of the density matrix and this leads to the following expresion:

A = 1 Z lim P →∞ mP 2πβ 2 P/2 dx 1 • • • dx P A(x 1 ) × exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β P V (x k ) (2.45) with x P +1 = x 1 = x.
The integral above can be rexpressed in a more convenient form observing that we favor one particular position variable (x 1 ). The fact to evaluate A only in one point makes the convergence slower. Due to the fact that in the case we are considering the paths are cyclic, all the points x 1 • • • x P are equivalent. As a consequence, the previous equation becomes:

A = 1 Z lim P →∞ mP 2πβ 2 P/2 dx 1 • • • dx P 1 P P k=1 A(x k ) × exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β P V (x k ) (2.46) with x P +1 = x 1 .
This expression allows us to treat each point equally and minimizes the statistical error on the evaluation of our expectation value. In order to obtain this quantity numerically (and so with a finite P ), it is interesting to express the equation 2.46 as the average value of A over a probability density. To do so, we define f (x 1 , • • • , x P ) as:

f (x 1 , • • • , x P ) = 1 Z P mP 2πβ 2 P/2 × exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β P V (x k ) (2.47)
(2.48)

with x P +1 = x 1 .
Z P is the partition function for a finite P as mentioned in the equation 2.31. We recover the fact that Z = lim P →∞ Z P . Furthermore, we can show easily that f (x 1 , • • • , x P ) ≥ 0 for all x 1 , • • • , x P and that:

dx 1 • • • dx P f (x 1 , • • • , x P ) = 1 (2.49) so f (x 1 , • • • , x P )
satisfies the conditions of a probability density and we can define a appropriate numerical estimator of A as:

A P (x 1 , • • • , x P ) = 1 P P k=1 A(x k ) (2.50)
Then we can define the expectation value of A for a finite P as the average of the numerical estimator according to the probability density f :

A P = dx 1 • • • dx P f (x 1 , • • • , x P )A P (x 1 , • • • , x P ) = A P (x 1 , • • • , x P ) f (2.51)
and where the limits P → ∞ is satisfied when numerical convergence with respect to P is observed.

A = lim P →∞ A P (2.52)
From a numerical point a view, all we need is to evaluate A P (x 1 , • • • , x P ) if we are able to sample correctly f . This is done using PIMD or PIMC as explained in the previous section.

The second case to consider is A a function of the momentum operator only: A = A( p). In this situation, it is no longer possible to express A in terms of diagonal elements of the density matrix in the coordinate basis. We have to introduce an identity operator between A and exp(-β H) because we cannot evaluate simply x| A:

A = 1 Z dx x| Ae -β H |x = 1 Z dx dx x| A|x x |e -β H |x = dx dx x| A|x ρ(x, x ; β) (2.53)
Then if we substitute the equation 2.23 with x 1 = x and x P +1 = x , in the previous equation, we obtain the following path integral expression for the expectation value of A:

A = 1 Z lim P →∞ mP 2πβ 2 P/2 dx 1 • • • dx P +1 x 1 | A|x P +1 × exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β 2P (V (x k+1 ) + V (x k )) (2.54)
Then we can introduce the identity operator in the momentum basis to look at the 20 Chapter 2. State of the art matrix element x 1 | A|x P +1 :

x 1 | A|x P +1 = dp x 1 | A|p p|x P +1 = 1 2π dp e ip(x 1 -x P +1 )/ A(p) (2.55)
and then we can introduce this result in the previous equation. Calculating expectation value of momentum operators makes appear two problems. The first one is that our paths are no longer cyclic (x 1 = x P +1 ) and in general sums over open paths are more difficult to evaluate than sums over closed paths due more important fluctuations in the endpoints. This leads us to the second problem which is the phase factors in the previous equation. If the endpoints have large fluctuations it could be difficult to make these phase factors converge. However, in certain cases it possible to do the previous integration analytically. We can mention the momentum distribution n(p), which is obtained by taking:

A( p) = δ( p -p I) (2.56)
where p is a number and I the identity operator.

In this particular cases, we have:

x 1 | A|x P +1 = 1 2π e ip (x 1 -x P +1 )/
This distribution can be mesured in neutron scattering experiments and by PIMD algorithms [START_REF] Sebastiani | Proton momentum distribution in water: an open path integral molecular dynamic study[END_REF]. Last but not least, we will focus on thermodynamic functions and more particularly on the expectation value of the energy. These functions are special in the sense that they are functions of both coordinate and momentum operators but can be evaluated via cyclic path integrals. Indeed, we have:

E = H = p 2 2m + V ( x) (2.57)
We can see that E depends of both position and momentum operators but we can avoid this problem via the thermodynamic relation:

E = - ∂ ∂β lnZ = 1 Z ∂Z ∂β (2.58)
We only need cyclic paths to express Z which implies that it is the same to calculate E. Taking the derivative of equation 2.30 with respect to β, we obtain the following expression for the energy:

E = lim P →∞ 1 Z P mP 2πβ 2 P/2 dx 1 • • • dx P E P (x 1 , • • • , x P ) × exp - 1 P k=1 mP 2β (x k+1 -x k ) 2 + β P V (x k ) (2.59) E = lim P →∞ E P (x 1 , • • • , x P ) f (2.60)
where

E P (x 1 , • • • , x P )
is the numerical estimator of the energy and f the probability density previously defined.

E P (x 1 , • • • , x P ) = P 2β - P k=1 mP 2β 2 2 (x k+1 -x k ) 2 + 1 P P k=1 V (x k ) (2.61) E P (x 1 , • • • , x P
) is an estimator (called primitive estimator) of the energy for finite P and so in the limit P → ∞ converges to the correct thermodynamic energy E. However, from a numerical point of view we can notice that our estimator is linearly dependent on P and as a consequence suffers from large fluctuations. So the energy of highly quantum systems which require a large number of discretization are difficult to converge with such an estimator. I will mention here but not demonstrate it the virial theorem for path integrals:

P 2β - P k=1 mP 2β 2 2 (x k+1 -x k ) 2 ρ = 1 2P P k=1 x k ∂V ∂x k ρ (2.62)
which leads to the definition of a new estimator of the energy, called the virial estimator:

E vir (x 1 , • • • , x P ) = 1 P P k=1 V (x k ) + 1 2 x k ∂V ∂x k (2.63)
In this formulation of the energy estimator, the kinetic energy is eliminated and this yields to an estimator with a much smaller variance and better convergence behaviour.

In conclusion, even if there are some subtilities and difficulties, evaluating time-independent equilibrium properties with path integrals formulation is a quite simple task. I did not present in this section the technical numerical aspects of path integrals to calculate these properties. Indeed the main aspects of it will be presented later in this Chapter (during the presentation of RPMD and CMD) and in Chapter 3 when I shall describe our method in detail.

Dynamical properties from quantum correlation functions

Real time quantum dynamics simulation of large molcular systems present a challenge to theoretical physics and chemistry. Indeed, computer simulation of time-dependent quantities for quantum systems is hindered by the exponential scaling of exact methods with the number of degrees of freedom. However, access to these quantities is crucial in many interesting problems and in particular to compute, via linear response theory and time correlation functions, the response of the system to external perturbations. The development of approximate quantum dynamical schemes is therefore an active area of research. One of the main directions currently pursued is to adapt efficient techniques for classical simulations, based on molecular dynamics (MD) and/or Monte Carlo, to the quantum case. A reasonable compromise among accuracy and computational cost is realised by the so-called quasi-classical methods for computing time correlation functions. In these methods, sampling of the thermal equilibrium density for the system is combined with generalised trajectories that provide an approximation for quantum dynamics usually valid for short times.

In this section, we will present different quasi-classical methods, containing different approximations. First, we will introduce the different forms of quantum correlation functions which is possible to compute. These expressions are equivalent in the sense that they can be related to one another analytically and therefore have the same physical content. It is however useful to discuss them all since, depending of the adopted approximation, some may be more advantageous from a computational point of view. Then, we will introduce two of the most popular quasi-classical methods, which, although not formally justified, have the advantage of converging with a small number of trajectories: centroid and ring polymer molecular dynamics (CMD and RPMD, respectively) [START_REF] Hernandez | Quantum time correlation functions and classical coherence[END_REF][START_REF] Ian | Quantum statistics and classical mechanics: real time correlation functions from ring polymer molecular dynamics[END_REF][START_REF] Rossi | How to remove the spurious resonances from ring polymer molecular dynamics[END_REF][START_REF] Markland | Annual Review of Physical Chemistry[END_REF][START_REF] Cao | A new perspective on quantum time correlation functions[END_REF][START_REF] Jang | Path integral centroid variables and the formulation of their exact real time dynamics[END_REF]. Finally, we will describe three methods, which are part of the category of the linearized quasi-classical methods. The Feynman-Kleinert linearized path integral (FK-LPI) method introduced by Poulsen and Rossky [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Aage Poulsen | Static and dynamic quantum effects in molecular liquids: a linearized path integral description of water[END_REF][START_REF] Poulsen | Quantum density fluctuations in liquid neon from linearized path-integral calculations[END_REF]36,[START_REF] Tyler D Hone | Comparison of approximate quantum simulation methods applied to normal liquid helium at 4 K. The journal of physical chemistry[END_REF][START_REF] Kyle | Refinement of the experimental dynamic structure factor for liquid para-hydrogen and ortho-deuterium using semi-classical quantum simulation[END_REF], the linearized semi-classical initial value representation (LSC-IVR) with the local Gaussian approximation (LGA) describe by Liu and Miller [4,[START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF][START_REF] Liu | Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initiobased flexible and polarizable force field[END_REF][START_REF] Liu | Test of the consistency of various linearized semiclassical initial value time correlation functions in application to inelastic neutron scattering from liquid para-hydrogen[END_REF][START_REF] Liu | Linearized semiclassical initial value time correlation functions with maximum entropy analytic continuation[END_REF][START_REF] Liu | Quantum dynamical effects in liquid water: A semiclassical study on the diffusion and the infrared absorption spectrum[END_REF] and finally our method, the phase integration method (PIM) introduced by Bonella and Ciccotti [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF]9,[START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]17,18].

Time correlation functions play a crucial role in relating macroscopic observations accessible in experiments to the microscopic dynamics of physical systems.

Although the statstical mechanics required to establish this link, in particular linear response theory, is valid for classical and quantum systems, our ability to perform calculations differs dramatically in the two situations as explained before. Furthermore, for quantum systems, different definitions of correlation functions are possible due to the fact in quantum statistical we deal with operators which do not commute.

The first and most common possibility is to evaluate the standard quantum expression of the time correlation function:

C AB (t, β) = 1 Z T r e -β H A(0) B(t) (2.64)
where β = 1 k B T , Z is the canonical partition function (defined in the previous section) and A(0) and B(t) are Heisenberg-evolved quantum operators at time 0 and time t:

B(t) = e it H/ Be -it H/
(2.65)

In the next sections, during the description of the different methods, we will explain and show how the dynamical part is approximated.

In the equation 2.64, the correlation function is written in the standard way as the thermal average of the operator product A(0) B(t). However, this is not the only possibility. A more symmetric alternative introduced by Kubo [42,[START_REF] Kubo | Statistical physics II. Non equilibrium statistical mechanics[END_REF][START_REF] Kubo | Statistical physics[END_REF] is the Kubo-transformed correlation function:

K AB (t, β) = 1 βZ β 0 dλ T r e -(β-λ) H A(0)e -λ H B(t) (2.66)
In the equation 2.66, the Boltzmann operator is averaged between A(0) and B(t). This differs from C AB because the operators A and B do not necessarly commute with H. This form of the quantum correlation is the closest to the classical one and as a consequence is well suited to the semi-classical approximations.

A last alternative form of the quantum time correlation function has been proposed by Schofield in 1960 [START_REF] Schofield | Space-Time Correlation Function Formalism for Slow Neutron Scattering[END_REF] originally for the neutron scattering. Schofield s function, also known as the symmetrized correlation function, is defined as:

G AB (t, β) = 1 Z T r Ae it * c H/ Be -itc H/ (2.67)
where t c = ti β/2. In the equation 2.67, the Boltzmann and time operators are equally distributed between the operators A and B. Schofield s form of the correlation function shares some formal properties with classical correlation functions, it is for example a real function by construction. Thus suggests that it could be a good starting point for describing semiclassical systems [21,46,47,3,48,36]. These three different definitions of quantum correlation function are linked together Chapter 2. State of the art via their Fourier transform. Indeed, we have:

CAB (ω) = e β ω/2 GAB (ω) CAB (ω) = β ω 1 -e β ω KAB (ω)
(2.68)

GAB (ω) = β ω/2 sh (β ω/2) KAB (ω)
where the time Fourier transform is:

F (ω) = +∞ -∞ dt e -iωt F (t) (2.69)
To obtain the relations in equation 2.68, we have to evaluate the traces in the different definitions in the basis of energy eigenstates. As an example, we carry out the calculation for the standard and the symmetrized correlation function. We have:

C AB (t, β) = 1 Z T r e -β H Ae it H/ Be -it H/ (2.70) = 1 Z n n|e -β H Ae it H/ Be -it H/ |n = 1 Z n e -βEn n| Ae it H/ B|n e -itEn/
where |n is an eigenvector in the energy basis and E n an eigenvalue in this basis, so H|n = E n |n . Then it follows after inserting the identity operator in the energy basis between the operators A and B (these operators are usually not diagonal in this basis set):

C AB (t, β) = 1 Z n m e -βEn
n| A|m e itEm/ m| B|n e -itEn/ (2.71)

= 1 Z n m
A nm B mn e -βEn e it(Em-En)/

where A nm = n| A|m and B mn = m| B|n . If we proceed in the exact same way for the symmetrized correlation function we obtain:

G AB (t, β) = 1 Z n m A nm B mn e -β 2 (En+Em) e it(Em-En)/ (2.72)
Then we take the Fourier transform of the equation 2.70 which leads us to:

CAB (ω) = 1 Z n m A nm B mn e -βEn δ ω - E m -E n (2.73)
where we used:

+∞ -∞ e -iωt dt e it(Em-En)/ = δ ω - E m -E n (2.74)
and for the Schofield s form of the correlation function we have:

GAB (ω) = 1 Z n m A nm B mn e -β 2 (En+Em) δ ω - E m -E n (2.75) = 1 Z n m A nm B mn e -βEn e -β ω/2 δ ω - E m -E n
It is now straightforward to see that the relation in the equation 2.68 is verified:

CAB (ω) = e β ω/2 GAB (ω) (2.76)

Empirical quasi-classical methods: RPMD and CMD

Now we know what we have to calculate but as explained before we cannot do it fully quantum mechanically. As a consequence we have to perform approximations to calculate these quantum correlation functions. We will present two of the most popular path integral methods to compute quantum time correlation function: RPMD and CMD. As we are going to see, these methods are not formally justified but are very efficient numerically. Then, we will describe the linearization approximation and different methods, including PIM, which are based on this approximation.

Outline of RPMD

In the section 2.1.2, we described the isomorphism between the path integral formalism and a ring polymer which is evolving under a classical Hamiltonian at a temperature P T where P is the number of beads in our polymer. This isomorphism is the basis of what we call the path integral molecular dynamics (PIMD) (or Path Integral Monte Carlo: PIMC), as we have explained, via the following equations of motions:

ẋk = ∂H P ∂p k ṗk = - ∂H P ∂x k (2.77) ẋk = p k m ṗk = - mP 2 β 2 2 (2x k -x k+1 -x k-1 ) - ∂V ∂x k
The idea (and approximation) of the Ring Polymer Molecular Dynamic (RPMD) method, first introduced by Manolopoulos in 2004 [START_REF] Ian | Quantum statistics and classical mechanics: real time correlation functions from ring polymer molecular dynamics[END_REF][START_REF] Rossi | How to remove the spurious resonances from ring polymer molecular dynamics[END_REF][START_REF] Markland | Annual Review of Physical Chemistry[END_REF], is to use these classical trajectories also to approximate Kubo-transformed quantum correlation functions using the physical mass m of the system in the equations 2.77.

For the case in which the correlated operators A and B are local (functions of the coordinate operator x ), this RPMD approximation is simply: First, in the high-temperature limit, the harmonic spring force constant mω 2 P in equation 2.34 becomes so large that the radius of gyration of the ring polymer shrinks to zero as we have already shown. Under these circumstances, it suffices to use just a single bead (P = 1), in which case the equation 2.78 is equivalent to the classical correlation function. Second, the ring-polymer correlation function in equation 2.78 coincides with the exact Kubo-transformed quantum mechanical correlation function in equation 2.66 in the limit as t → 0. To see this in detail, one can expand the correlation functions in the two equations in Taylor series around t = 0. Only the even expansion coefficients survive in this case because the exact K AB (t) and its RPMD approximation are both real and even functions of t. Comparison of the exact and approximate expansion coefficients reveals that the RPMD approximation has a leading error of O(t 8 ) for the position autocorrelation function [START_REF] Mark | Statistical Mechanics: Theory and Molecular Simulation[END_REF][START_REF] Markland | Annual Review of Physical Chemistry[END_REF][START_REF] Perez | A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals[END_REF].

K AB (t, β) = 1 βZ β 0 dλ T r e -(β-λ) H A(0)e -λ
Third, when the simple case of the harmonic oscillator (V (x) = 1 2 mωx 2 ) is considerated, it is possible to derive the equation 2.78 to show that it gives the exact quantum mechanical results (with P → ∞) for correlation functions of the form K Ax (t) and K xB (t). Last but not least when we consider the special case of A as the unit operator, it is possible to show that we are only calculating the equilibrium expectation value of the operator B. The equation 2.78 becomes in this special case: ) and the phase space volume element (d P pd P x = d P p t d P x t ) in the classical ring polymer dynamics. Relabeling p t and x t as p and x, one sees that this is just the P -bead path integral expression for B . The RPMD approximation therefore becomes exact for A = 1 in the limit as P → ∞. The last of these limiting cases is especially important for applications of the RPMD model to condensed-phase systems because it confirms that the classical dynamics of the ring polymers is at least consistent with the quantum mechanical equilibrium distribution, as one would expect from the connection between RPMD and PIMD. This implies that an RPMD simulation at a given NVT thermodynamic state point will sample the correct initial quantum distribution. For the moment only operators of postion have been discussed. However, we can easily show that for operators linear in momentum such as the velocity auto-correlation function, we obtain the following expression within the RPMD approximation:

K 1B (t, β) ≈ 1 (
K vv (t, β) ≈ 1 (2π ) P Z P d P p d P x exp - β P H P (x, p) vv t (2.83)
where v = 1 P m P k=1 p k . In both cases, when we deal with operators linear in position or momentum we only have to evaluate their value for the centroid position at time 0 and time t. The RPMD approximation gives quite good results in many cases but there is no rigorous justification for the adoption of this dynamics (at least not for generic potentials and observables), so results can have uncontrolled features. The first one is encountered in infrared spectroscopy calculations due to resonance with spurious frequencies. This phenomena will be illustrated in Chapter 4 on toy model calculations. The occurence of this phenomena is quite simple to explain after a normal mode transformation of the ring-polymer Hamiltonian. Indeed, if we consider simply a one dimensional harmonic oscillator with a frequency ω (V (x) = 1 2 mω 2 x 2 ), we obtain after a normal mode transformation (which will not be specified but is explain in the following references [START_REF] Mark | Statistical Mechanics: Theory and Molecular Simulation[END_REF][START_REF] Perez | A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals[END_REF]):

Chapter 2. State of the art

H P (x, p) = P k=1 1 2m p2 k + 1 2 mω 2 k x2 k (2.84)
where x, p are the normal mode coordinates and where:

ω k = ω 2 + 2P β 2 sin 2 ((k -1)π/P ) (2.85)
In RPMD, the centroid mode (k = 1) therefore oscillates at the harmonic frequency ω, whereas the other normal modes oscillate at higher frequencies that depend on the number of beads and on the temperature. These internal-mode oscillations are unrelated to the dynamics of the physical potential; they are artificial oscillations that arise from the structure of the extended phase space. In many applications, the high-frequency oscillations of the ring polymer do not matter, because they are well separated from the dynamics of physical interest. However, problems do arise when using RPMD to calculate absorption spectra in systems containing high-frequency physical vibrations, as illustrated in Chapter 4. The lesson to be learned from this is that the spectral information produced by RPMD is not to be trusted at frequencies above the first free ring-polymer excitation frequency, ω 2 = (2P/β)sin(π/P ) ≈ 2π/β , which corresponds to a wave number of approximately 1.300cm -1 at 300 K. This spurious frequency problem has other consequences.Indeed, a second important limitation of RPMD is its poor performance for the correlation functions of nonlinear operators [START_REF] Markland | Annual Review of Physical Chemistry[END_REF]. Consider the calculation of correlation functions for the simple one dimension harmonic oscillator. RPMD gives the exact quantum mechanical result when one of the correlated operators is a linear function of position or momentum. However, in the case of the nonlinear operator x 2 , for example, the exact quantum Kubo-transformed autocorrelation function is:

K x 2 x 2 = 2mω 2 2 β ω coth β ω 2 cos(2ωt) + 2 coth 2 β ω 2 -1 (2.86)
and the RPMD result is:

K x 2 x 2 = 1 βm 2 P k=1 1 ω 4 k (cos(2ω k t) + 1) + P k=1 P l=1 1 ω 2 k ω 2 l (2.87)
where ω k and ω l have the same definition has for the spurious frequency problem. From these equations, we see that, although the RPMD autocorrelation function does have a component (k = 1 which correspond to the centroid of the ring-polymer) that oscillates at the correct frequency (2ω), there are additional, spurious components that arise from the internal modes of the ring polymer and contaminate the time signal. Many drawbacks of RPMD have been presented but do not limit the success of RPMD. Indeed, as we can see in the review of Manolopoulos [START_REF] Markland | Annual Review of Physical Chemistry[END_REF], RPMD has been applied successfully to many applications such as quantum diffusion in liquid parahydrogen, liquid water, proton transfer in polar solvent, gas phase chemical rate ... From a computational point of view, RPMD is quite cheap because it's only a classical dynamic where we multiply the number of degrees of freedom by the number of beads needed. Another advantage is that RPMD is purely Newtonian and the convergences of the dynamical properties are achieved with few trajectories due to the fact that the quantum density is preserved during the dynamic (and so tricks such as Wiener-Khintchine can be used [50,51]). Indeed, Witt et al. [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF] showed that in order to obtain a converged IR spectra you may need as few as 50 trajectories.

Outline of CMD

Another popular semi-classical method to approximate quantum time correlation function is the Centroid Molecular Dynamic (CMD) method [START_REF] Hernandez | Quantum time correlation functions and classical coherence[END_REF][START_REF] Cao | A new perspective on quantum time correlation functions[END_REF][START_REF] Jang | Path integral centroid variables and the formulation of their exact real time dynamics[END_REF]. In order to present in a synthetic way this method, we need to introduce the functional integrals formalism.

If we consider that we are in the P → ∞ limit, we can construct a mathematical picture of path integrals known as a functional integral. For simplification we introduce a parameter ε = β /P and as a consequence P → ∞ implies ε → 0. If we rewrite the density matrix element in terms of ε we obtain:

ρ(x, x ; β) = lim P →∞ ε→0 m 2πε P/2 1 Z P dx 2 • • • dx P × exp - ε P k=1 m 2 x k+1 -x k ε 2 + 1 2 (V (x k+1 ) + V (x k ))
(2.88)

In the limit P → ∞ and ε → 0 the intervals between the points x 1 , x 2 , • • • , x P , x P +1 become infinitively small while the number of points becomes infinite. We need to introduce a continuous function x(s) satisfying x(0) = x and x(β ) = x via the following identification:

x k = x(s = (k -1)ε) (2.89)
In this continuous limit, we can write:

lim ε→0 x k+1 -x k ε = dx ds (2.90)
Chapter 2. State of the art as a consequence we have:

lim ε→0 ε P k=1 m 2 x k+1 -x k ε 2 + 1 2 (V (x k+1 ) + V (x k )) = β 0 ds m 2 dx ds 2 + V (x(s)) (2.91)
and we can define:

S[x] = β 0 ds m 2 dx ds 2 + V (x(s)) (2.92)
S[x] corresponds to the action (the integral of the Lagrangian) in imaginary time.

The density matrix is obtained by integration over all paths x(s) satisfying x(0) = x and x(β ) = x weighted by exp(-S[x]/ ). The points x 1 , • • • , x P +1 comprise all the points of the function x(s) in the limit P → ∞ with x 1 = x and x P +1 = x . Thus, the integration over x 2 , • • • , x P constitues an integration over all possible function x(s) that satisfy the endpoints conditions x(0) = x and x(β ) = x . This type of integration is referred to as functional integration and is mathematically a delicate object. Indeed, the measure is known for imaginary time propagation, not for real time. Symbolically, it is written as follows:

lim P →∞ ε→0 m 2πε P/2 dx 2 • • • dx P = Dx(s) (2.93) ρ(x, x ; β) = 1 Z x(β )=x x(0)=x Dx(s) exp [-S[x]/ ] (2.94)
We can now express the canonical partition function in the functional integral formalism as:

Z(β) = dx x(β )=x x(0)=x Dx(s) exp [-S[x]/ ] = Dx(s) exp [-S[x]/ ] (2.95)
where Dx(s) indicates that the functional integral is taken over all the possible paths with x(0) = x(β ).

After this brief introduction on functional integrals formalism, we can present the Centroid Molecular Dynamics (CMD) introduced by J. Cao in 1993 [START_REF] Cao | A new perspective on quantum time correlation functions[END_REF]. This is a technique for computing real-time quantum correlation functions. The primary object in this approach is the path centroid defined as:

x c = 1 P P k=1 x k = 1 β β 0 x(s) (2.96)
the second equation is true when we are in the continuous limit so when P → ∞.

Feynman and Kleinert introduced in 1986 [START_REF] Feynman | Effective classical partition functions[END_REF] the centroid potential of mean force and the idea of CMD is based on the notion that the time evolution of the centroid on this potential of mean force can be used to calculate approximate quantum dynamical properties of a system. In CMD, coming back to the simple case of one particle in one dimension, the centroid is postulated (approximation without any justification) to evolve in time according to the following equations of motion:

dx c dt = p c m (2.97) dp c dt = - dV 0 (x c ) dx c = F 0 (x c ) (2.98)
where m is the physical mass of the particle, p c is a momentum conjugate to x c and V 0 (x c ) is the centroid potential of mean force given by:

V 0 (x c ) = - 1 β ln 2πβ 2 m 1/2 ρ c (x c ) (2.99) ρ c (x c ) = Dx(s)δ(x 0 [x(s)] -x c )e -S[x(s)]/ (2.100)
where x 0 [x(s)] = 1 β ds x(s). V 0 (x c ) corresponds to the exact excess quantum free energy of the centroid and ρ c (x c ) is the centroid density and we can make the link with the canonical partition function:

Z(β) = dx c ρ c (x c ) (2.101)
In equation 2.100, S[x(s)] is the action in imaginary time defined previously and the δ function restrict the functional integration to cyclic paths whose centroid position is x c . The centroid force at x c is derived by spatial differentiation from the equation 2.100 and given by:

F 0 (x c ) = - Dx(s)δ(x 0 [x(s)] -x c ) 1 β β 0 ds V (x(s )) e -S[x(s)]/ Dx(s)δ(x 0 [x(s)] -x c )e -S[x(s)]/
(2.102)

Although formally exact within the CMD framework equations 2.100 and 2.102 require a full path integral calculation at each centroid configuration, which is computationally very demanding. The technical way to alleviate this problem will not be discussed (is done via a normal mode transformation and an adiabatic approximation [START_REF] Cao | A new perspective on quantum time correlation functions[END_REF][START_REF] Mark | Statistical Mechanics: Theory and Molecular Simulation[END_REF]). In CMD, the centroid potential of mean force is generated "on the fly". However, one important thing to notice at this stage from the equation 2.100 is that the centroid potential is dependent on the temperature which as we will see later (see Chapter 4) is a source of problems in CMD.

The key assumption of CMD is that the Kubo-transformed quantum time correlation function K AB (t, β) can be approximate by the centroid dynamic by (considering A and B operators of position only):

K AB (t, β) = 1 Z(β) dx c dp c A(x c )B(x c (t)) exp -β p 2 c 2m + V 0 (x c ) (2.103)
where B(x c (t)) is evaluated using the time evolved centroid variables generated by the equations 2.97. A similar definition holds for operators function of momentum. CMD is exact in the classical limit and in the purely harmonic case. Away from these limits, position autocorrelation function are accurate for short times up to O(t 6 ) [START_REF] Mark | Statistical Mechanics: Theory and Molecular Simulation[END_REF][START_REF] Markland | Annual Review of Physical Chemistry[END_REF]. Furthermore, CMD is computationally quite cheap because we do the dynamic of the centroid of the polymer only. The quantity expensive to compute is the centroid force but as explained before there are ways to calculate it "on the fly" and so the numerical cost can be controlled. Similar to RPMD, the convergence to calculate dynamical properties is quickly achieved and, for example, for IR spectra of simple molecules, 30 trajectories are enough to obtain a converged spectra [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF]. CMD has been applied successfully to many quantum calculations such as the quantum diffusion liquid para-hydrogen, neutron scattering ... [START_REF] Perez | A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals[END_REF][START_REF] Ramírez | The Schrödinger formulation of the Feynman path centroid density[END_REF]54]. However, the method lacks a rigorous justification in the sense that there is no formal demonstration that this is the correct dynamic for the quantum system. As a consequence, it can produce unexpected and pathological results. For instance, as we are going to explain in more details in Chapter 4, IR spectra computed from CMD are not reliable due to unphysical features originating from the force calculations.

Linearized methods

As we have seen previously, two popular quasiclassical methods exists, which, although not formally justified, have the advantage of converging with a small num-ber of trajectories: centroid and ring polymer molecular dynamics. However, these methods have also the drawback to present intrinsinc limitations which are difficult to predict and mitigate. An alternative family of approaches instead, develops approximate quantum methods starting from the exact formulations of the problem and developing algorithms from controlled formal manipulations. In particular our method (PIM for Phase Integration Method), which I will present in the next Chapter, and the methods described in this section are part of the so called linearized methods. To set the stage for PIM, I will present first the linearization approximation in general and then the Feynman Kleinert Linearized Path Integral (FK-LPI) method introduced by Poulsen and Rossky in 2003 [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF] and the Linearized Semi-Classical Initial Value Representation (LSC-IVR) method using the Local Gaussian Approximation (LGA) introduced by Liu and Miller [START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF][START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF].

The Linearized Path Integral (LPI) representation of quantum correlation functions

I will present here the linearized path integral approximation on the standard quantum correlation function but the idea is stricly the same for the Schofield or the Kubo correlation functions. We first recall the equation 2.64, which corresponds to the definition of the standard quantum correlation function:

C AB (t, β) = 1 Z T r e -β H Ae it H/ Be -it H/ (2.104)
where β = 1 k B T and Z is the canonical partition function. We introduce resolutions of identity in the coordinate basis to isolate matrix elements of the different operators. Thus, we shall make appear in particular a forward (denoted with x + ) real time path and a backward (denoted with x -) real time path.

C AB (t, β) = 1 Z dx + P +1 dx - P +1 dx + 0 dx - 0 x + 0 |e -β H A|x - 0 × x - 0 |e it H/ |x - P +1 x - P +1 | B|x + P +1 x + P +1 |e -it H/ |x + 0 (2.105)
We are going to focus on only the real time evaluation for the moment because the linearized approximation is done on this part. The imaginary time part of the correlation function will be described just after the LPI approximation. Within the different linearized methods, this part is more difficult to calculate and makes the originality of each method previously cited (PIM, FK-LPI and LSC-IVR).

We insert P times (where we consider P suffienctly high to be in the limit P → ∞) the resolution of identity in the coordinate basis to "resolve" the forward path x + P +1 |e -it H/ |x + 0 and the backward path x - 0 |e it H/ |x - P +1 as:

x + P +1 |e -it H/ |x + 0 x - 0 |e it H/ |x - P +1 = dx + 1 • • • dx + P dx - 1 • • • dx - P × x + P +1 |e -iε H/ |x + P • • • x + 1 |e -iε H/ |x + 0 x - 0 |e iε H/ |x - 1 • • • x - P |e iε H/ |x - P +1
(2.106)

where ε = t P +1 . Then we insert P + 1 times the resolution of identity in the momentum basis to the forward and backward paths which allow us to obtain a hybrid coordinatemomentum path integral representation of the real time propagators. If we detail this step of the calculation on one term of the forward path we have:

x + k+1 |e -iε H/ |x + k ≈ x + k+1 |e -iε p 2 2m e -iε V |x + k (2.107) = x + k+1 |e -iε p 2 2m |x + k e -iε V (x + k ) = 1 2π dp + k+1 x + k+1 |e -iε p 2 2m |p + k+1 p + k+1 |x + k e -iε V (x + k ) = 1 2π dp + k+1 e i p + k+1 (x + k+1 -x + k ) e -iε (p + k+1 ) 2 2m e -iε V (x + k )
For the backward path we simply have:

x - k |e iε H/ |x - k+1 = x - k+1 |e -iε H/ |x - k * (2.108) = 1 2π dp - k+1 e -i p - k+1 (x - k+1 -x - k ) e iε (p - k+1 ) 2 2m e iε V (x - k )
If we regroup all the terms, we obtain:

x + P +1 |e -it H/ |x + 0 x - 0 |e it H/ |x - P +1 = 1 (2π ) 2(P +1) dx + 1 • • • dx + P dx - 1 • • • dx - P × dp + 1 • • • dp + P +1 dp - 1 • • • dp - P +1 exp i P k=0 p + k+1 (x + k+1 -x + k ) -p - k+1 (x - k+1 -x - k ) × exp - i ε P k=0 (p + k+1 ) 2 2m - (p - k+1 ) 2 2m exp - i ε P k=0 V (x + k ) -V (x - k ) (2.109)
Then we proceed to a change of variables with a mean and a difference variable:

xk = x + k + x - k 2 (2.110) pk = p + k + p - k 2
(2.111)

∆x k = x + k -x - k (2.112) ∆p k = p + k -p - k (2.113)
For k = 0 to P . In the new variables (the Jacobian of this change is equal to 1), we have the following relations:

(p + k+1 ) 2 -(p - k+1 ) 2 = 2p k+1 ∆p k+1 (2.114) P k=0 p + k+1 (x + k+1 -x + k ) -p - k+1 (x - k+1 -x - k ) (2.115) = pP +1 ∆x P +1 -p1 ∆x 0 - P k=1 ∆x k (p k+1 -pk ) + P k=0 ∆p k+1 (x k+1 -xk )
Substituting in equation 2.109:

x + P +1 |e -it H/ |x + 0 x - 0 |e it H/ |x - P +1 = 1 (2π ) 2(P +1) dx 1 • • • dx P dp 1 • • • dp P +1 × exp i (p P +1 ∆x P +1 -p1 ∆x 0 ) d∆x 1 • • • d∆x P exp - i P k=1 ∆x k (p k+1 -pk ) × exp - i ε P k=0 V (x k + ∆x k 2 ) -V (x k - ∆x k 2 ) d∆p 1 • • • d∆p P +1 × exp - i ε P k=0 pk+1 ∆p k+1 m exp i P k=0 ∆p k+1 (x k+1 -xk ) (2.116)
Then knowing that:

+∞ -∞ dx e -ix(a-b) = 2πδ(a -b) (2.117)
we can notice the following simplification:

d∆p 1 • • • d∆p P +1 exp - i ε P k=0 pk+1 ∆p k+1 m exp i P k=0 ∆p k+1 (x k+1 -xk ) = (2π ) P +1 P k=0 δ ε pk+1 m -xk+1 + xk (2.118)
Then we can explicit the linearization approximation simply by saying that we consider the forward and backward paths very close to each other. The mathematical consequence of this assumption is that we can do a Taylor expansion of the potential around the mean variable in coordinates. This is fully justified if we are at short time (and exact for harmonic potential). Indeed, the paths are closed (with operators diagonal in positions or momenta) at the endpoints, i.e. at time 0 and time t. So, if we are at a time t close to 0 we can assume that the paths remain close to each other between the endpoints which is not the case anymore if t is far from 0. It is also possible to show (see [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Bonella | Linearization approximations and Liouville quantum-classical dynamics[END_REF] for more details) that this is equivalent to an expansion of the potential around the mean path. We have:

V (x k + ∆x k 2 ) = V (x k ) + ∆x k 2 ∇V (x k ) (2.119) V (x k - ∆x k 2 ) = V (x k ) - ∆x k 2 ∇V (x k ) (2.120) V (x k + ∆x k 2 ) -V (x k - ∆x k 2 ) = ∆x k ∇V (x k ) (2.121)
In our path integral formalism, the truncated expansion has the consequence to make appear delta function in the part of the path integrals where the potential is involved:

d∆x 1 • • • d∆x P exp - i P k=1 ∆x k (p k+1 -pk ) × exp - i ε P k=0 V (x k + ∆x k 2 ) -V (x k - ∆x k 2 ) = d∆x 1 • • • d∆x P exp - i P k=1 ∆x k (p k+1 -pk ) × exp - i ε P k=1 ∆x k ∇V (x k ) exp - i ε∆x 0 ∇V (x 0 ) (2.122)
Then it is straightforward to apply the same operation as before:

d∆x 1 • • • d∆x P exp - i P k=1 ∆x k (p k+1 -pk ) exp - i ε P k=1 ∆x k ∇V (x k ) exp - i ε∆x 0 ∇V (x 0 ) = (2π ) P P k=1 δ (p k+1 -pk + ε∇V (x k )) (2.123)
We can now regroup the results above in the correlation function expression:

C AB (t, β) = 1 2π Z dx 0 • • • dx P +1 dp 1 • • • dp P +1 d∆x 0 d∆x P +1 × x0 + ∆x 0 2 e -β H A x0 - ∆x 0 2 xP +1 - ∆x P +1 2 B xP +1 + ∆x P +1 2 (2.124) × exp i (p P +1 ∆x P +1 -p1 ∆x 0 ) P k=1 δ (p k+1 -pk + ε∇V (x k )) P k=0 δ ε pk+1 m -xk+1 + xk
Note that above we have dropped the term expi ε∆x 0 ∇V (x 0 ) . In the limit P → ∞ this amounts to neglecting one term in an infinite sum and does not affect the result appreciably. The variables x1 • • • xP +1 can now, one by one, be integrated out as follows. For some value of x0 and p1 , the delta function δ ε p1 m -x1 + x0 fixes the value of x1 . Then the delta function δ (p 2 -p1 + ε∇V (x 1 )) determines the value of p2 and so on. In this manner all the variables are integrated out. More importantly, they are, in addition, seen to fullfill, to first order in ε, the classical equations of motion where ε is the time step. Thus, the value of xP +1 = x t and pP +1 are simply equal to the classically propagated value of x0 = x 0 and p1 = p 0 . We can rewrite the remaining ∆ variable as ∆x P +1 = ∆x t and we arrive at the linearized path integral representation of correlation functions:

C AB (t, β) = 1 2π Z dx 0 dp 0 d∆x 0 d∆x t exp i (p t ∆x t -p 0 ∆x 0 ) × x 0 + ∆x 0 2 e -β H A x 0 - ∆x 0 2 x t - ∆x t 2 B xt + ∆x t 2 (2.125)
which can be written in a more compact way if we use the Wigner transform [START_REF] Wigner | On the Quantum Correction For Thermodynamic Equilibrium[END_REF] of a generic operator O:

O W (x, p) = d∆e ip∆/ x - ∆ 2 O x + ∆ 2 (2.126) so:
Chapter 2. State of the art

C AB (t, β) = 1 2π Z dx 0 dp 0 e -β H A W (x 0 , p 0 )B W (x t , p t ) (2.127)
This representation of quantum correlation function is exact for harmonic systems, in the high temperature limit and in the short time limit [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Bonella | Linearization approximations and Liouville quantum-classical dynamics[END_REF].

The general implementation of the equation 2.127 can be described as follows. The LPI-CF is evaluated by first sampling x 0 and p 0 from the Wigner transform of e -β H A followed by classical propagation up to a time t. The Wigner transform of B is then evaluated with the phase-space arguments (x t , p t ). This is repeated for a sufficiently large set of sampled initial conditions and the final result obtained as the mean value over all the trajectories. By using the LPI approximation, we have arrived at a simple path integral expression which does not suffer from any sign problem connected to the real time dynamics of the system. However, there is still an oscillatory phase or sign problem connected to the Wigner transform of the Boltzmann operator, which is essentially a many-dimensional Fourier transform. Ways to deal with this problem are the subject of the next sections.

The FK-LPI method

This method was first introduced Poulsen and Rossky in 2003 [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF] and is based on the calculation of the Wigner transform of the Boltzmann operator via the Feynman-Kleinert method [START_REF] Feynman | Effective classical partition functions[END_REF]. Then, this Wigner transform is utilized with the LPI representation of correlation functions described in the previous section. As a consequence, we will present briefly in this section the way used to calculate the Wigner transform of the Boltzmann operator used by Poulsen and Rossky [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Aage Poulsen | Static and dynamic quantum effects in molecular liquids: a linearized path integral description of water[END_REF][START_REF] Poulsen | Quantum density fluctuations in liquid neon from linearized path-integral calculations[END_REF]36].

As we can see in the equation 2.127, we need to evaluate the Wigner transform of e -β H A, which is not, in general, possible. To circumvent this problem Poulsen and Rossky used a method of representing the density operator in a semi-harmonic form which can be Wigner transformed analytically. The approach is built on combining two distinct theories: the first is the variational harmonic frequency representation of path integrals as developed by Feynman and Kleinert [START_REF] Feynman | Effective classical partition functions[END_REF]. The second is the centroid constrained quasi density operator (QDO) formalism of Jang and Voth [START_REF] Jang | Path integral centroid variables and the formulation of their exact real time dynamics[END_REF].

We start by focusing on the Wigner transform of e -β H alone. If we use the hybrid coordinate-momentum path integral representation and the functional integral formalism, we have [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF]:

exp -β H = dxdx |x ρ(x, x ; β) x| (2.128) ρ(x, x ; β) = x(β )=x x(0)=x Dx(s)Dp(s) exp [-S[x, p]/ ] (2.129) S[x, p] = β 0 ds T (p(s)) + V (x(s)) -ip(s) dx ds (2.130)
where T (p(s)) = p(s) 2 2m . Then, we can introduce in the Boltzmann operator a constraint over the centroid of the path and the momentum centroid (x c = 1 β β 0 ds x(s) and a similar definition for p c ):

exp -β H = dx c dp c (2.131) × dxdx x(β )=x x(0)=x Dx(s)Dp(s) δ(x -x c )δ(p -p c ) exp - S[x, p] |x x|
where x and p denote specific values of the path centroids.

The quantity inside the square bracket represent the part of the Boltzmann operator where the centroid paths are constrained to x c and p c . Based on the idea of Feynman and Kleinert [START_REF] Feynman | Effective classical partition functions[END_REF], for these specific value, the action S[x, p] is localized around x c and p c and therefore, a local harmonic approximation is reasonable. We transform transform S[x, p] into a harmonic form using the following approximation of the potential:

1 2 mΩ 2 (x c )(x(s) -x c ) 2 + L 1 (x c ) for V (x(s)).
So the action has now the following form:

S[x, p] = β 0 ds T (p(s)) + 1 2 mΩ 2 (x c )(x(s) -x c ) 2 + L 1 (x c ) -ip(s) dx ds (2.132)
The optimal frequency Ω(x c ) and potential L 1 (x c ) were derived by Feynman and Kleinert [START_REF] Feynman | Effective classical partition functions[END_REF], based on a variational approximation to the partition function. This derivation will not be specified here but all we need to know is that these quantities are optimized iteratively via a variational approach.

We can now obtaine an approximate expression for the thermal propagator:

exp -β H = dx c dp c dxdx x(β )=x x(0)=x Dx(s)Dp(s) |x x|δ(x -x c )δ(p -p c ) × exp - 1 β 0 ds T (p(s)) + 1 2 mΩ 2 (x c )(x(s) -x c ) 2 + L 1 (x c ) -ip(s) dx ds (2.133)
Then, we can define the Feynman-Kleinert approximation to the centroid density as:

ρ F K (x c , p c ) = T r dxdx x(β )=x x(0)=x Dx(s)Dp(s) |x x|δ(x -x c )δ(p -p c ) × exp - 1 β 0 ds T (p(s)) + 1 2 mΩ 2 (x c )(x(s) -x c ) 2 + L 1 (x c ) -ip(s) dx ds (2.134)
The power of this approximation appears now in the sense that due to the fact that the path integral contains a harmonic action, it can be explicitly evaluated, yielding (not detailed here, see [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Feynman | Effective classical partition functions[END_REF] for more details) to:

ρ F K (x c , p c ) = 1 2π exp -β p 2 c 2m + W 1 (x c ) (2.135)
Where the definition of the so-called centroid potential is:

W 1 (x c ) = L 1 (x c ) + 1 β ln sinh 2 βΩ(x c ) 2 βΩ(x c ) (2.136)
The Feynman-Kleinert approximation to canonical partition function is therefore:

Z F K = dx c dp c ρ F K (x c , p c ) (2.137)
We can also define the FK (Feynman-Kleinert) quasi-density operator (QDO) [START_REF] Jang | Path integral centroid variables and the formulation of their exact real time dynamics[END_REF]:

δ F K (x c , p c ) = 1 ρ F K (x c , p c ) dxdx x(β )=x x(0)=x Dx(s)Dp(s) |x x|δ(x -x c )δ(p -p c ) × exp - 1 β 0 ds T (p(s)) + 1 2 mΩ 2 (x c )(x(s) -x c ) 2 + L 1 (x c ) -ip(s) dx ds (2.138) By construction δ F K (x c , p c
) is of unit trace, which allows us to write compactly:

exp -β H = dx c dp c δ F K (x c , p c )ρ F K (x c , p c )
The interest of this approximation stems from the fact that we are not able to perform the Wigner transform of the full density operator but we can do it analytically on the QDO. An explicit formula of this QDO has been provided by Jang and Voth [START_REF] Jang | Path integral centroid variables and the formulation of their exact real time dynamics[END_REF] δ

F K (x c , p c ) = dxdx mΩ(x c ) π α |x x| exp - mΩ(x c )α 4 (x -x) 2 × exp   ip c (x -x) - mΩ(x c ) α x + x 2 -x c 2   (2.139)
where α is defined as:

α = coth Ω(x c ) β 2 - 2 Ω(x c ) β (2.140)
and α takes values between zero and unity. We can now perform analytically the Wigner transform of δ F K (x c , p c ):

[ δ F K (x c , p c )] W (x, p) = 2 α exp - mΩ(x c ) α (x -x c ) 2 - 1 mΩ(x c )α (p -p c ) 2 (2.141)
For large temperature, α → 0 so x and p will be fixed at x c and p c . In that case, all the samplings stay in ρ F K (x c , p c ). As ρ F K (x c , p c ) tends to the classical phase space density in the high temperature limit, one obtains the classical density function. At the other extreme, when T → 0, the x c and p c values sampled from ρ F K (x c , p c ) will lie in the harmonic minimum, and all quantum fluctuations are determined from equation 2.141 with α = 1.

From equation 2.141, the Wigner transform of the Boltzmann operator can be written as:

exp -β H W (x, p) = dx c dp c ρ F K (x c , p c ) 2 α × exp - mΩ(x c ) α (x -x c ) 2 - 1 mΩ(x c )α (p -p c ) 2 (2.142)
This is the key result of the FK-LPI method. From it, the Wigner transforms of operators of the general form exp -β H A are readily analytically evaluated if A is sufficiently simple. For example, if we consider the position operator A = x, we have:

exp -β H x W (x, p) = dx c dp c ρ F K (x c , p c ) 2 α x + i p -p c mΩ(x c )α × exp - mΩ(x c ) α (x -x c ) 2 - 1 mΩ(x c )α (p -p c ) 2 (2.143) Chapter 2.

State of the art

And finally, if we want to calculate the position autocorrelation function, we obtain from the equations 2.127 and 2.143:

C F K-LP I xx (t) = dx c dp c ρ F K (x c , p c ) dx 0 dp 0 2 α x 0 + i p 0 -p c mΩ(x c )α x t × exp - mΩ(x c ) α (x -x c ) 2 - 1 mΩ(x c )α (p -p c ) 2 (2.144)
The practical computational procedure for calculating correlation function using the FK-LPI method then goes as follows. First position and momentum centroids are sampled from ρ F K (x c , p c ) . After that, x 0 and p 0 are sampled from equation 2.142. These are then propagated classically and the Wigner transform B W (x t , p t ) is constructed.

The FK-LPI method has been applied for example on the calculations of linear and nonlinear correlation functions in Helium liquid [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Tyler D Hone | Comparison of approximate quantum simulation methods applied to normal liquid helium at 4 K. The journal of physical chemistry[END_REF], liquid Neon [START_REF] Poulsen | Quantum density fluctuations in liquid neon from linearized path-integral calculations[END_REF] or liquid para-hydrogen [START_REF] Kyle | Refinement of the experimental dynamic structure factor for liquid para-hydrogen and ortho-deuterium using semi-classical quantum simulation[END_REF].

The Linearized Semiclassical Initial Value Representation (LSC-IVR) method with a Local Gaussian Approximation (LGA)

The LSC-IVR [21,57,22] approximation is stricly equivalent to the LPI approximation detailed before. The originality of the method develop by Liu and Miller is to combine this approximation with the Local Gaussian Approximation (LGA) [START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF][START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF] for sampling the Wigner density. As a consequence, I will insist here on the Local Gaussian Approximation and why it is a good method to compute Wigner densities. I remind the here the LSC-IVR approximation (or LPI) result for the standard quantum correlation function in the multidimensional case:

C AB (t, β) = (2π ) -N Z dx 0 dp 0 e -β H A W (x 0 , p 0 )B W (x t , p t ) (2.145)
where the Wigner transform has been defined in the equation 2.126, N is the number of degrees of freedom and x 0 = (x

(1) 0 , • • • , x (N ) 0 ) (idem for p 0 ).
Calculation of the Wigner function for operator B in equation 2.127 is usually straightforward; in fact, B is often a function only of coordinates or only of momenta, in which case its Wigner functions is simply the classical function itself. However, as shown in the previous section, calculating the Wigner transform of e -β H A is far from being trivial. Nonetheless, this is necessary in order to obtain the distribution of initial conditions of coordinates x 0 and momenta p 0 for the real time trajectories. Usually, local harmonic approximations of the potential (such as the Local Harmonic Approximation (LHA) of Shi and Geva [START_REF] Shi | Semiclassical Theory of Vibrational Energy Relaxation in the Condensed Phase[END_REF], the FK-LPI of Poulsen and Rossky [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF], ...) are adapted and give good results to many quantum simulations. These methods, however, fail when the imaginary local frequency and temperature are such that β |ω| ≥ π and so are not reliable for problems dominated by potential barriers at low temperature.

The LGA for the momentum distribution attemps to remedy this problem by modifying the LHA of Shi and Geva [START_REF] Shi | Semiclassical Theory of Vibrational Energy Relaxation in the Condensed Phase[END_REF]. First a factorization of the diagonal matrix element of the Boltzmann operator (which can be evaluated accurately by path integral techniques) is performed in the Wigner function (or density) W (x, p):

W (x, p) = 1 2π Z d∆x x - ∆x 2 e -β H x + ∆x 2 e ip∆x/ (2.146) W (x, p) = 1 2π Z x|e -β H |x d∆x x -∆x 2 e -β H x + ∆x 2 x|e -β H |x e ip∆x/ (2.147)
and then the LHA [START_REF] Shi | Semiclassical Theory of Vibrational Energy Relaxation in the Condensed Phase[END_REF] is done for the ratio of off-diagonal elements to diagonal matrix elements:

x -∆x 2 e -β H x + ∆x 2 x|e -β H |x ≈ exp - mω 4 coth(β ω/2)∆x 2 (2.148)
where ω is the local frequency so:

ω = 1 m d 2 V dx 2 (2.149)
It is quite obvious in this expression that if the curvature of our potential is negative, the local frequency will be imaginary. This the case of problems dominated by barrier potentials such as the 1D problem of Eckart barrier [21,[START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF]. This approximation is exact for harmonic systems. The LHA of equation 2.148 makes the Fourier transform of the equation 2.146 analytical and gives us for the Wigner distribution function within the LHA:

W (x, p) = 1 2π Z x|e -β H |x β 2πmQ(u) 1/2 exp -β p 2 2mQ(u) (2.150)
where u = β ω is a dimensionless parameter and :

Q(u) = u/2 tanh[u/2] (2.151)
is the so-called quantum correction factor.

In the high temperature limit, β → 0, the classical limit, → 0, and the free particle limit, ω → 0, the parameter u → 0 and, as a consequence, the quantum correction factor Q(u) → 1. In these limits, the equation 2.150 gives the classical momentum distribution as expected.

In the region where the local frequency ω is imaginary, the parameter u becomes imaginary and u = iu i where u i = β |ω|. Thus:

Q(u) = Q(iu i ) = u i /2 tan[u i /2]
(2.152) ! Figure 2.6 shows this quantum correction factor Q(u) as the solid line for both real and imaginary frequencies (on the positive and negative u axes, respectively), and one sees the breakdown of the LHA in the imaginary frequency regime when u i ≥ π. One simple ad hoc procedure to deal with the regime u i ≥ π is to set Q(u) = 0 in this domain which also means that the momentum p = 0 since the momentum distribution in equation 2.150 becomes an infinitely sharp Gaussian, i.e., a delta function. However, momentum distribution even at low temperature are not delta functions so this suggests that the quantum correction factor shown in figure 2.6 should decrease smoothly to 0 as u i becomes greater that π as shown by the dashed curve in figure 2.6. The Local Gaussian Approximation defines one way to enforce this smooth approach to zero and constitutes the difference between LHA and LGA. There are obviously many ad hoc choices one can make to modify the quantum correction factor in the imaginary frequency regime . One reasonable requirement is that it should agree with the harmonic result in equation 2.152 for small u so:

Q(u) ≈ 1 + u 2 12 + O(u 4 ) (2.153)
for u real and imaginary. The form that Liu and Miller [START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF] found to be both simple and to obey the previous equation is the following:

Q(u) =        u/2 tanh(u/2) for real u 1 Q(u i ) = tanh(u i /2) u i /2 for imaginary u (u = iu i ) (2.154)
which is plotted as the dashed line in figure 2.6. The utility of this ad hoc post treatment of the LHA arises at low temperature for problems dominated by imaginary frequencies.

Once we perform the LGA on the Wigner transform of the Boltzmann operator (or Wigner density), we can obtain an expression for quantum time correlation functions. Indeed, due to the LGA we can express analytically the Wigner transform of e -β H A as:

e -β H A W = (2π ) x 0 |e -β H |x 0 β 2πmQ(u) 1/2 exp -β p 2 0 2mQ(u) f A (x 0 , p 0 ) (2.155)
where

f A (x 0 , p 0 ) = d∆x x 0 -∆x 2 |e -β H A|x 0 + ∆x 2 x 0 |e -β H |x 0 e ip∆x/ d∆x x 0 -∆x 2 |e -β H |x0+ ∆x 2 x 0 |e -β H |x 0 e ip∆x/
(2.156)

f A (x 0 , p 0 ) can be calculated anatically if A is simple enough. For example, if A = x f A (x 0 , p 0 ) as the following expression:

f A (x 0 , p 0 ) = x 0 + iβ 2mQ(u) p 0 (2.157)
where Q(u) is the quantum correction factor previously defined. Finally, the explicit form of the LSC-IVR correlation function with the LGA in one dimension is given by:

C AB (t) = dx 0 x 0 |e -β H |x 0 Z dp 0 β 2πmQ(u) 1/2 exp -β p 2 0 2mQ(u) f A (x 0 , p 0 )B(x t , p t ) (2.158)
We can summarize the procedure to calculate the LSC-IVR correlation function with the LGA as follows:

1. Use path integral Monte Carlo (PIMC) or path integral molecular dynamics

(PIMD) to calculate x 0 |e -β H |x 0 Z .
2. At specific intervals in the PIMC (or time steps in the PIMD), randomly select one path integral bead as the initial configuration x 0 for the real time dynamics and calculate the local frequency.

3. The LGA gives the Gaussian distribution for the momentum

β 2πmQ(u) 1/2 exp -β p 2 0 2mQ(u)
which is used to sample initial momentum p 0 for real time trajectories.

4. Run real time classical trajectories from phase space points (x0, p0) and estimate the property f A (x 0 , p 0 )B(x t (x 0 , p 0 ), p t (x 0 , p 0 )) for the corresponding time correlation function.

5. Repeat steps 2-4 and sum the property f A (x 0 , p 0 )B(x t (x 0 , p 0 ), p t (x 0 , p 0 )) for all real time classical trajectories until a converged result is obtained.

Conclusion

To sum up this Chapter, we have described the path integral formalism and see how it could be apply to the calculation of time-independent equilibrium properties. Furthermore, the power of this formalism, conceptual and numerical, has been introduced via the classical isomorphism which allows us to describe a quantum particle as a "necklace" of classical particles.

In a second part, we introduced the different approximations used currently to compute quantum time correlation functions. We first presented RPMD and CMD which are two popular methods to calculate quantum correlation functions. The popularity of these methods comes from the fact that they are numerically cheap and can be applied successfully on many problems. However, for both methods, there is no formal justification for the approximate dynamics and so they present limitations which are difficult to predict. An alternative family of approaches instead, develops approximate quantum methods starting from the exact formulations of the problem and developing algorithms from controlled formal manipulations. These are the methods based on the LPI representation of quantum correlation function. We first introduced the FK-LPI [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF], this method is built on combining two distinct theories: the first is the variational harmonic frequency representation of path integrals as developed by Feynman and Kleinert [START_REF] Feynman | Effective classical partition functions[END_REF]. The second is the centroid constrained quasi density operator (QDO) formalism of Jang and Voth [START_REF] Jang | Path integral centroid variables and the formulation of their exact real time dynamics[END_REF]. Then we presented the LSC-IVR with the LGA introduced by Liu and Miller [START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF]. The LGA can be seen as a patch to the LHA [START_REF] Shi | Semiclassical Theory of Vibrational Energy Relaxation in the Condensed Phase[END_REF] to treat correctly local imaginary frequencies. The LSC-IVR or LPI representation of quantum correlation function can treat both linear and nonlinear operators in a consistent way, can be applied to nonequilibrium as well as the above equilibrium correlation functions, and can also be used to describe electronically nonadiabatic dynamics, i.e., processes involving transitions between several potential energy surfaces [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]9,18,[START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Aage Poulsen | Static and dynamic quantum effects in molecular liquids: a linearized path integral description of water[END_REF][START_REF] Poulsen | Quantum density fluctuations in liquid neon from linearized path-integral calculations[END_REF]36,[START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF][START_REF] Liu | Test of the consistency of various linearized semiclassical initial value time correlation functions in application to inelastic neutron scattering from liquid para-hydrogen[END_REF][START_REF] Liu | Linearized semiclassical initial value time correlation functions with maximum entropy analytic continuation[END_REF][START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF][START_REF] Liu | Quantum dynamical effects in liquid water: A semiclassical study on the diffusion and the infrared absorption spectrum[END_REF][START_REF] Liu | Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initiobased flexible and polarizable force field[END_REF]. These merits of LPI representation fo correlation function make it a versatile tool to study a variety of quantum mechanical effects in chemical dynamics of large molecular systems.

Contents

3.1 Wigner densities via the phase integration method (PIM) . 50

Wigner densities via the phase integration method (PIM)

In 1932, Wigner introduced his representation of quantum operators to formulate quantum statistical mechanics in a language with important analogies with the classical case. Indeed, in Wigner's representation, the average of a quantum operator A, defined as we have seen in Chapter 2 as T r[ ρ A], where ρ is the density operator of the system, can be identically rewritten as an integral over momenta and coordinates variables, thus:

A = T r[ ρ A] = dqdp A w (q, p)W (q, p) (3.1)
where we recall the definition of the Wigner transform of the operator A:

A w (q, p) = d∆ q - ∆ 2 A q + ∆ 2 e i p∆ (3.2)
In the equation 3.1, we also introduced the key quantity in this framework: the Wigner transform of the density operator W (q, p), which has been already defined in the previous Chapter (see equation 2.146). We recall here the definition:

W (q, p) = 1 2π Z d∆r ν e i p∆rν q - ∆r ν 2 e -β Ĥ q + ∆r ν 2 (3.3)
where Z is the partition function and Ĥ = P 2 /2m + V (q) the Hamiltonian of the system, which we will assume composed of distinguishable particles of mass m. The subscript on the displacement variable, ∆r ν , has been introduced for future convenience. (One dimensional notation will be used for simplicity here and in the following.) Unfortunately, W (q, p) is a real but not necessarily positive definite function. This prevents its direct interpretation as a probability density and compromises, to a large extent, the conceptual usefulness of the analogy between Wigner's formulation and classical statistical mechanics. Furthermore, while the Wigner transform of common operators (functions of either coordinate or momenta of the system) can be easily computed analytically, there is no analytical expression for W (q, p) for systems subject to generic potentials. As discussed in Chapter 2, the Wigner density is also extremely difficult to obtain via numerical methods mainly due to the phase factor in its definition. In contrast, Feynman's path integral formalism [START_REF] Feynman | Space-time approach to non-relativistic quantum mechanics[END_REF] provides a representation for average values that, thanks to the so called classical isomorphism, is straightforwardly amenable to efficient classical simulation methods [START_REF] Chandler | Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids[END_REF]. This is indeed the reason why path integrals have emerged as the tool of choice for computing quantum time independent averages. In the past two decades, however, the interest in Wigner's formalism has been renewed in the context of developing methods for computing quantum time correlation functions. Exact algorithms to solve this problem scale exponentially with the number of degrees of freedom and cannot be used for realistic applications. Thus, considerable efforts are currently devoted to develop computationally affordable approximations as we have explained in Chapter 2. Many of these are based on the idea to combine sampling of the exact quantum thermal density with classical trajectories. In particular, the so called linearisation schemes [START_REF] Hernandez | Quantum time correlation functions and classical coherence[END_REF][START_REF] Wang | Semiclassical approximations for the calculation of thermal rate constants for chemical reactions in complex molecular systems[END_REF][START_REF] Sun | Semiclassical theory of electronically nonadiabatic dynamics: Results of a linearized approximation to the initial value representation[END_REF][START_REF] Shi | Semiclassical Theory of Vibrational Energy Relaxation in the Condensed Phase[END_REF][START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Serena Causo | An adiabatic linearized path integral approach for quantum time correlation functions: electronic transport in metal-molten salt solutions[END_REF] introduced previously rely on sampling the Wigner transform of the thermal density (or, more precisely, the Wigner transform of the density times an operator) to obtain a set of initial conditions that are then evolved classically. The approximate correlation function is computed by averaging, over the different trajectories, appropriate observables. Although linearisation methods fail in highly quantum conditions (for example when coherence dominates the dynamics), they have been shown to provide surprisingly good results for condensed phase systems [START_REF] Aage Poulsen | Static and dynamic quantum effects in molecular liquids: a linearized path integral description of water[END_REF][START_REF] Poulsen | Quantum density fluctuations in liquid neon from linearized path-integral calculations[END_REF]36,[START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF][START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF]. In fact, these methods derive from an approximation to the full quantum evolution which is valid for short times (the LPI approximation described in the previous chapter). In the condensed phase, where the decay time of correlations can occur on even shorter time scales, linearised methods can indeed prove effective. Unfortunately, as we have seen in the previous Chapter, even when the approximation on the dynamics is reliable, linearised calculations are non trivial, due to the difficulty to compute the Wigner transform of the density to sample initial conditions. Different schemes have been proposed to overcome this problem, as discussed already in Chapter 2. Shi and Geva, for example, introduced a sampling scheme based on a local harmonic approximation (LHA) of the potential that is often used for multidimensional systems [START_REF] Shi | Semiclassical Theory of Vibrational Energy Relaxation in the Condensed Phase[END_REF], which was later modified -via an ansatz to treat imaginary frequencies -by Liu and Miller [START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF] via the LGA. Poulsen et al. [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF] combined the variational harmonic frequency path integral representation of Feynman and Kleinert [START_REF] Feynman | Effective classical partition functions[END_REF] with the centroid quasi density formalism developed by Jang and Voth [START_REF] Jang | Path integral centroid variables and the formulation of their exact real time dynamics[END_REF] to obtain an approximate form of the Wigner density which can be computed via a relatively simple iterative scheme. Alternative approaches compute the Wigner function by solving Bloch's equation in phase space space by propagating in imaginary time Gaussian phase space packets, starting from the classical Boltzmann distribution (i.e. the high temperature limit of the Wigner density). This is the so called Gaussian phase space packet (GPSP) approach recently proposed by Marinica et al. [START_REF] Marinica | Generating approximate Wigner distributions using Gaussian phase packets propagation in imaginary time[END_REF]. More recently, the quantum thermal bath method (QTB), proposed independently by Dammak et al. [START_REF] Dammak | Quantum thermal bath for molecular dynamics simulation[END_REF] and Ceriotti et al. [START_REF] Ceriotti | Nuclear quantum effects in solids using a colored-noise thermostat[END_REF] and in which a generalized Langevin equation with colored noise tuned to satisfy the fluctuation dissipation theorem is used to thermalize the system, was suggested as an alternative method to generate the Wigner density [START_REF] Basire | Computing Wigner distributions and time correlation functions using the quantum thermal bath method: application to proton transfer spectroscopy[END_REF]. Although the methods just mentioned have proved successful on given applications, they all introduce approximations that are often difficult to control both from a theoretical and numerical point of view.

To make progress, we introduce and investigate a new sampling scheme for the Wigner density, based on a non standard Monte Carlo method recently developed within the framework of time dependent quantum correlation function calculations [START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. To that end, we first rewrite the off diagonal density matrix in Eq. 3.3 in a convenient path integral representation, and then employ a cumulant expansion to control the phase factor. The path integral representation adopted and the cumulant expansion are the key, original, features of our approach. We then test the sampling scheme on a set of model problems of increasing complexity, and we will show that it captures highly non classical effects, most notably quantum correlation among positions and momenta (or even momenta of different degrees of freedom) of the system. At variance with alternative schemes, the approach only contains numerical approximations that can be systematically checked and improved upon. In particular, the method relies on the convergence of the cumulant expansion. We will show that, for the most physically meaningful examples considered, truncating this expansion at second order suffices to reproduce benchmark results. However, the convergence of the cumulant expansion (at infinite order) implies that the Wigner density for the system is positive. Thus, since it is known that this property does not always hold, there will be situations in which the cumulant expansion does not converge. We then also consider a variation on the approach which substitutes cumulants with an expansion of the density with respect to a convenient reference probability [START_REF] Cramér | Mathematical methods of statistics[END_REF]. This procedure is the analogous, in Fourier space, to the so-called Edgeworth expansion [START_REF] Cramér | Mathematical methods of statistics[END_REF] of the density, which can account for negative values of the Wigner function and is known to have different, and often better, asymptotic properties with respect to simple cumulants. In the following we shall refer to it as an Edgeworth expansion.

The Chapter is organised as follows. In the Method section we express the Wigner density in a form suitable for direct, even though non standard, Monte Carlo sampling. We then introduce the expansion with respect to the reference probability and show how to write an appropriate estimator for the expression of the Wigner density. Then I will described in a more detailed way how the algorithm of the non-standard Monte Carlo scheme works [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF]. In the section Results, we compare the Wigner density obtained for the set of models we have chosen (ranging from a simple harmonic oscillator to a model system for a bound proton in an A-H-B system) with reference quantum results and with some of the alternative approximate methods mentioned above.

The PIM expression for the thermal Wigner density

Let us begin by inserting a resolution of the identity in the coordinate basis to rewrite equation 3.3 as:

W (q, p) = 1 2π Z d∆r ν dr 0 e i p∆rν q - ∆r ν 2 e -β 2 Ĥ r 0 r 0 e -β 2 Ĥ q + ∆r ν 2 (3.4)
In the exact same way as we described in the previous Chapter, the matrix elements in the expression above can be represented in path integral form, via a symmetric Trotter break up. In fact, repeating the steps that lead to the expression of an off-diagonal element of the density matrix, we obtain:

q - ∆r ν 2 e -β 2 Ĥ r 0 = C ν dy 1 • • • dy ν-1 exp - 1 2σ 2 ν-1 λ=0 (y λ+1 -y λ ) 2 × exp - δ β 2 V (y ν ) exp - δ β 2 V (y 0 ) exp -δ β ν-1 λ=1 V (y λ ) (3.5)
and similarly we have:

r 0 e -β 2 Ĥ q + ∆r ν 2 = C ν dx 1 • • • dx ν-1 exp - 1 2σ 2 ν-1 λ=0 (x λ+1 -x λ ) 2 × exp - δ β 2 V (x ν ) exp - δ β 2 V (x 0 ) exp -δ β ν-1 λ=1 V (x λ ) (3.6)
where

δ β = β/2ν, σ 2 = 2 δ β /m, C ν = mν πβ 2 ν/2
, and with boundary conditions y 0 = x 0 = r 0 , y ν = q + ∆r ν /2 and y ν = q -∆r ν /2. For equations 3.5 and 3.6, the expression above becomes exact in the limit ν → ∞ (ν here is the number of beads, previously denoted as P ). We then introduce the change of variables:

r λ = x λ + y λ 2 (3.7) ∆r λ = x λ -y λ (3.8)
note that λ = 0, • • • , ν and with the definition above, r ν = q and ∆r 0 = 0. In the new variables, equation 3.4 becomes:

W (q, p) = C 2 ν 2π Z dr 0 • • • dr ν-1 d∆r 1 • • • d∆r ν e i p∆rν exp - 1 2σ 2 r ν-1 λ=0 (r λ+1 -r λ ) 2 × exp - σ 2 p 2 ν-1 λ=0 (∆r λ+1 -∆r λ ) 2 exp - δ β 2 [V (r ν + ∆r ν 2 ) + V (r ν - ∆r ν 2 )] × exp -δ β ν-1 λ=1 [V (r λ + ∆r λ 2 ) + V (r λ - ∆r λ 2 )] exp [-δ β V (r 0 )] (3.9)
or

W (q, p) = C 2 ν 2π Z dr d∆r e i p∆rν W(r, ∆r, q) (3.10)
where (in the first equation) we set σ 2 r = δ β /2m, σ 2 p = m/2δ β , and (in the second) we introduced the notation ∆r = {∆r 1 , ..., ∆r ν }, r = {r 0 , ..., r ν-1 }. We also defined W(r, ∆r, q) as the product of exponentials in the curly bracket on lines one,two and three of the equation 3.9. To set the stage for the Monte Carlo scheme we intend to apply, it is convenient to express the Wigner density as an expectation value. To that end, we begin by using the identity:

W (q * , p * ) = dqdp W (q, p)δ(q -q * )δ(p -p * ) (3.11) = C 2 ν 2π Z dqdp
dr d∆r e i p∆rν W(r, ∆r, q)δ(qq * )δ(pp * )

Introducing the δ functions in the expression above is not strictly necessary and most of the considerations that follow could be modified to avoid it. However, writing the Wigner density as in equation 3.11 allows us to proceed in complete analogy with previous work done with PIM [START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF] to first introduce the cumulant expansion mentioned before and detailed below, and then to construct the Monte Carlo algorithm. This simplifies somewhat the discussion and, perhaps more importantly, makes it possible to rely on the performance of the sampling algorithm introduced in those references when discussing the scaling properties of the method. To set the stage for our developments, note that W(r, ∆r, q) is a positive definite quantity which, barring a normalization factor N , can be interpreted as the joint probability density for all coordinate variables {r, ∆r, q}. (At this stage, it is in particular the inclusion of q in the set of variables that simplifies manipulations.) Let us indicate this density as:

ρ(r, ∆r) = W(r, ∆r)/N (3.12)
with a supplementary simplification due to the fact that q = r ν , we have now r = {r 0 , ..., r ν } . or, more explicitely:

ρ(r, ∆r) = 1 N exp - 1 2σ 2 r ν-1 λ=0 (r λ+1 -r λ ) 2 exp - σ 2 p 2 ν-1 λ=0 (∆r λ+1 -∆r λ ) 2 exp [-δ β V (r 0 )] × exp - δ β 2 [V (r ν + ∆r ν 2 ) + V (r ν - ∆r ν 2 )] exp -δ β ν-1 λ=1 [V (r λ + ∆r λ 2 ) + V (r λ - ∆r λ 2 )] (3.13)
We can express the density above in terms of marginal ρ m (r) and conditional ρ c (∆r|r) probabilities defined as::

ρ m (r) = d∆r ρ(r, ∆r) (3.14) = 1 N exp - 1 2σ 2 r ν-1 λ=0 (r λ+1 -r λ ) 2 exp [-δ β V (r 0 )] d∆r exp - σ 2 p 2 ν-1 λ=0 (∆r λ+1 -∆r λ ) 2 × exp - δ β 2 [V (r ν + ∆r ν 2 ) + V (r ν - ∆r ν 2 )] exp -δ β ν-1 λ=1 [V (r λ + ∆r λ 2 ) + V (r λ - ∆r λ 2 )]
and:

ρ c (∆r|r) = ρ(r, ∆r) ρ m (r) (3.15) = e - σ 2 p 2 ν-1 λ=0 (∆r λ+1 -∆r λ ) 2 e - δ β 2 [V (rν + ∆rν 2 )+V (rν -∆rν 2 )] e -δ β ν-1 λ=1 [V (r λ + ∆r λ 2 )+V (r λ - ∆r λ 2 )] d∆r e -σ 2 p 2 ν-1 λ=0 (∆r λ+1 -∆r λ ) 2 e - δ β 2 [V (rν + ∆rν 2 )+V (rν -∆rν 2 )] e -δ β ν-1 λ=1 [V (r λ + ∆r λ 2 )+V (r λ - ∆r λ 2 )]
Using these notations we can rewrite our Wigner density as:

W (q * , p * ) = N C 2 ν 2π Z dp dr ρ m (r) d∆r e i p∆rν ρ c (∆r|r) δ(r ν -q * )δ(p -p * ) (3.16)
Let us now consider more in detail the term in square brackets in the expression above. Its logarithm is, by definition, the generating function of the random variable ∆r ν with respect to ρ c (∆r|r) [START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF]. The integral over the ∆r variables can then be formally expressed as:

d∆r e i p∆rν ρ c (∆r|r) = e -E(p,r) (3.17) with E(p, r) = - ∞ n=1 (-ip/ ) n n! ∆r n ν c ρc(∆r|r) (3.18)
where ∆r n ν c ρc(∆r|r) is the nth cumulant. Since, as it can be seen from the definition of ρ(r, ∆r), ρ c (∆r|r) = ρ c (-∆r|r), only even order terms in equation 3.18 are non zero. Thus, when the cumulant series converges, E(p, r) is a real function and the exponential in equation 3.17 is positive definite. This implies, see equation 3.19 below, that also the Wigner density is positive definite. It is however known that there are systems (e.g. double well potential at low temperature) that have a non positive definite Wigner density. For such systems, the writing above fails because the cumulant series does not converge and our method cannot be applied. Importantly however, this convergence can always be tested numerically (an analytical assessment of the convergence of the cumulant series, on the other hand, is non trivial for generic potentials) ensuring control of the suitability and accuracy of the method. In the following subsection we discuss an alternative scheme that is applicable also for non positive definite Wigner densities. Substituting the cumulant form of the ∆r integral in equation 3.16, we obtain:

W (q * , p * ) = N C 2 ν 2π Z dp dr ρ m (r)e -E(p,r) δ(r ν -q * )δ(p -p * ) (3.19)
By reverse engineering the steps in equations 3.9-3.17 it is not difficult to see that:

N C 2 ν 2π Z dp dr ρ m (r)e -E(p,r) = 1 (3.20)
thus:

W (q * , p * ) = δ(q -q * )δ(p -p * ) P (3.21)
where the average is over the following quantity:

P(r, p) = N C 2 ν 2π Z ρ m (r)e -E(p,r) (3.22) P(r, p) = ρ m (r)e -E(p,r) dp dr ρ m (r)e -E(p,r) (3.23)
Due to the fact that E(p, r) is a real function and the exponential in equation 3.17 is positive definite, P(r, p) fulfills all the conditions to be defined as a probability density.

We will see in the next section how to deal with such a probability density from an algorithmic point of view.

Equation 3.19 can be most effectively read as the definition of an histogram and calculated by generating a set of values of (r ν , p) distributed according to the probability density introduced above. This probability density is, however, unusual since it contains two factors (ρ m (r) and e -E(p,r) ) that are, in general, not known analytically. This poses a problem for the numerical sampling of P for example via the standard Monte Carlo schemes such as the standard Metropolis scheme [START_REF] Teller | Equation of state calculations by fast computing machines[END_REF]. Indeed, in the Metropolis algorithm, the probability density needs to be known analytically in order to sample it. Notably, in general this probability density does not factor in a term depending only on coordinates and one depending only on momenta indicating that, in contrast to the classical distribution but also to other available approximate approaches, our expression of the Wigner density can account for correlation among these degrees of freedom. Perhaps even more strikingly compared to the classical case, since for more than one degree of freedom the cumulant expansion allows for coupling among different momenta, this distribution can also manifest correlation among the momenta. We shall see in the Results section that these features can be physically relevant.

Edgeworth expansion for the Wigner density

The manipulations presented above give, by construction, a positive definite Wigner density. While this is usually not considered a serious problem for most realistic applications (and indeed, essentially all existing approximation schemes have the same limitation), it is interesting to explore alternative representations that can reproduce also negative features of the Wigner density if they exist. To that end, let us consider again the integral over ∆r at the left hand side of equation 3.17.

It is the characteristic function of the conditional probability ρ c (∆r|r). We have previously introduced a cumulant expansion of this quantity, which, for all practical purposes, will be truncated in the following at second order thus approximating this integral by a Gaussian in p, with an r dependent variance (covariance matrix in the multidimensional case). By taking this Gaussian as a reference it is then possible to trivially rewrite the integral as: because ∆r ν ρc(∆r|r) = 0 due to the fact that ρ c (∆r|r) is even in ∆r ν ). Using the expansion [START_REF] Cramér | Mathematical methods of statistics[END_REF]:

d∆r e i p∆rν ρ c (∆r|r) = e -κ 2
e κ 2 2 p 2 2 e i p∆rν = +∞ n=0 H n ∆rν √ κ 2 κ n/2 2 n! ip n , (3.25)
(above H n (x) is the n th Hermite polynomial, and we omitted the dependence of κ 2 on the coordinates to simplify the notation), we obtain:

d∆r e i p∆rν ρ c (∆r|r) = e -κ 2 2 p 2 2 +∞ n=0 C n n! ip n (3.26)
with:

C n = d∆r H n ∆r ν √ κ 2 κ n/2 2 ρ c (∆r|r) (3.27)
Using the symmetry of the conditional probability ρ c , all odd terms in the series at the right hand side of equation 3.26 are equal to zero. By the choice of κ 2 , truncation at n = 2 is identical to the truncation at second order in the cumulant expansion considered in the previous subsection, since C 2 = 0. The next term in the series gives, recalling that H 4 (x) = x 4 -6x 2 + 3:

d∆r e i p∆rν ρ c (∆r|r) = e -κ 2 2 p 2 2 1 + C 4 4! ip 4 (3.28)
with:

C 4 = d∆r H 4 ∆r ν √ κ 2 κ 2 2 ρ c (∆r|r) (3.29)
= d∆r ∆r 4 ν κ 2 2 -6 ∆r 2 ν κ 2 + 3 κ 2 2 ρ c (∆r|r) = d∆r ∆r 4 ν ρ c (∆|r) -3κ 2 2
the fourth order cumulant of ρ c . Thus, truncating the expansion at the fourth power in p, we can approximate the Wigner distribution as:

W (q * , p * ) = N C 2 ν 2π Z dp dr ρ m (r)e -κ 2 (r) p 2 2 2 1 + C 4 4! ip 4 δ(r ν -q * )δ(p -p * ) (3.30)
Since the Wigner distribution is normalized, the following condition should be ensured:

N C 2 ν 2π Z dp dr ρ m (r)e -κ 2 (r) p 2 2 2 1 + C 4 4! ip 4 = 1 (3.31)
This condition is enforced by estimating the Wigner distribution as the ratio of expectation values:

W (q * , p * ) = 1 + C 4 4! ip 4 δ(r ν -q * )δ(p -p * ) P 2 (p,r) 1 + c 4 4! ip 4 P 2 (p,r) (3.32)

Structure of the PIM algorithm 59

where the sampling is performed on the probability distribution:

P 2 (r, p) = ρ m (r)e -κ 2 (r) p 2 2 2
dp dr ρ m (r)e -κ 2 (r) p 2 2 2

(3.33)

The sampling is thus performed on the same distribution as for the original PIM method (see previous subsection) when truncating at the second order in the cumulant expansion. The correction is here included in the estimator, not in the sampling, the more so since the correction term can take negative values and therefore so can our estimator of the Wigner function. The changes in the sign of the estimator may add extra noise in the evaluation of the Wigner distribution but it will be shown in the next section that this noise can be controlled in many cases using the same number of samplings as for the original PIM. When this is not the case, the correction still provides qualitative information on the existence of negative parts in the Wigner density.

It should be noted that the expansion of the characteristic function proposed here is equivalent by Fourier transform to the Gram-Charlier A-series of the conditional probability ρ c (∆r|r) or, after reordering of the terms, to the Edgeworth expansion of this conditional probability, the two being identical for a symmetric distribution when truncating at the fourth order [START_REF] Cramér | Mathematical methods of statistics[END_REF]. Because of better convergence properties [START_REF] Cramér | Mathematical methods of statistics[END_REF], the Edgeworth expansion would probably be more suitable would one want to extend the expansion to higher order terms. Indeed, it has been shown, that the Edgeworth expansion is a true asymptotic series so that the error is, at least in principle, under control after truncation both when it converges and when it does not [START_REF] Cramér | Mathematical methods of statistics[END_REF][START_REF] Juszkiewicz | Weakly non linear Gaussian fluctuations and the Edgeworth expansion[END_REF][START_REF] Blinnikov | Expansions for nearly Gaussian distributions[END_REF]. In the following, we will refer to equation 3.32 as an "Edegworth correction" to PIM.

Structure of the PIM algorithm

Sampling P(r, p) (or P 2 (r, p)) is non trivial since both ρ m (r) and e -E(p,r) can only be estimated numerically. Given their expressions, we will show that however, they can be calculated as an average with an associated variance. Due to this standard Monte Carlo methods, which rely on the analytical knowledge of the density, cannot be directly applied. We are going to describe how to circumvent this problem by combining two schemes for sampling "noisy" probability densities like the one in equation 3.22. These methods are the penalty [START_REF] Ceperley | The penalty method for random walks with uncertain energies[END_REF] and Kennedy [START_REF] Clark | Asymptotics of fixed point distributions for inexact Monte Carlo algorithms[END_REF] Monte Carlo algorithms. Both schemes use two main steps: (1) appropriate, unbiased, numerical estimators of the probability densities are defined to generate trial values for the random variables; (2) the trial values are accepted or rejected based on a generalized acceptance criterion which, on average, corrects for the effect of the noise.

For convenience in the description of the algorithm, we introduce the following notations:

e -Vr(r) = exp - 1 2σ 2 r ν-1 λ=0 (r λ+1 -r λ ) 2 (3.34) e -V ∆ (∆r) = exp - σ 2 p 2 ν-1 λ=0 (∆r λ+1 -∆r λ ) 2 (3.35) e -δ β V (r,∆r) = exp [-δ β V (r 0 )] e - δ β 2 [V (rν + ∆rν 2 )+V (rν -∆rν 2 )] e -δ β ν-1 λ=1 [V (r λ + ∆r λ 2 )+V (r λ - ∆r λ 2 )] (3.36)
which allow us to rewrite the marginal (ρ m (r)), the conditional (ρ c (∆r|r)) and the noisy probability density in a more synthetic way as:

ρ m (r) = d∆r ρ(r, ∆r) (3.37) = 1 N e -Vr(r) d∆r e -δ β V (r,∆r) (3.38) = 1 N e -Vr(r) ρ m (r) (3.39)
ρ c (∆r|r) = ρ(r, ∆r) ρ m (r) (3.40)
= e -V ∆ (∆r) e -δ β V (r,∆r) d∆r e -V ∆ (∆r) e -δ β V (r,∆r) (3.41) P(r, p) = e -Vr(r) ρ m (r)e -E(p,r) dp dr e -Vr(r) ρ m (r)e -E(p,r)

(3.42)

As stated before, our goal is to sample numerically P (note that p and r in this probability density are not independent and have to be treated together in the sampling). We use a Monte Carlo scheme in which the probability to generate a new state of the system (by changing momenta or coordinates of a particle) and accept/reject this new state is generalized to statisfy the detailed balance even in the presence of significant noise in the estimates of ρ m (r) and E(r, p). I will present first the important steps of the algorithm and the two different noisy Monte Carlo schemes (Penalty and Kennedy methods) and then show how we calculate the different estimators that we need for our algorithm.

Structure of the algorithm

The first step of the algorithm is to choose if we move the momentum or coordinates.

Let us consider that we move the momentum. In this case, we choose p = p+δp , δp being an uniform random vector centered at zero (the magnitude of the displacement is chosen so as to optimize the acceptance). With this prescription, keeping into account that the r variables are not being updated, detailed balance takes the form:

P(r, p)T p (p → p )P p acc (p → p ) = P(r, p )T p (p → p)P p acc (p → p) (3.43)
where P p acc is the acceptance probability of the displacement and T p the, uniform, transition probability.

We can manipulate the previous equation by simplifying all the terms which depend exclusively of r in P(r, p) and using the fact that, for uniform displacement, we have This detailed balance is similar to the one considered by Ceperly et al. [START_REF] Ceperley | The penalty method for random walks with uncertain energies[END_REF]when they introduce the penalty method. This is in fact a generalized Monte Carlo for sampling a density given by the exponential of a function known with a statistical error (in our case E(p, r)). In this case the procedure is the following, if we have a probability density P(s) of this form:

T p (p → p ) = T p (p → p) = 1 ∆ if |p -p| < ∆ 2 .
P(s) ∝ e -L(s) with L(s) estimated (3.45)
Then we perform the calculation of the difference of the exponent at the current, s, and proposed,s , state N times (as L(s) is estimated the result will be different each time we calculate this quantity) and take the average and the variance of this quantity. Thus, we compute

δ i (s, s ) = L i (s ) -L i (s) (3.46) D(s, s ) = 1 N N i=1 δ i (s, s ) (3.47) χ 2 (s, s ) = 1 N (N -1) N i=1 (D(s, s ) -δ i (s, s )) 2 (3.48) where i = 1, • • • , N .
According to Ceperly et al. [START_REF] Ceperley | The penalty method for random walks with uncertain energies[END_REF], asymptotic sampling of P(s) is achieved by using an acceptance probability of following form:

a(s → s ) = min 1, exp(-D(s, s ) -u χ 2 ) (3.49) with u χ 2 = χ 2 2 + χ 4 4N + 1 + χ 6 3(N + 1)(N + 3) + • • • (3.50)
This acceptance test differs from the standard Metropolis rule for the presence of u χ 2 and is valid when χ 2 4 ≤ 1: in the limit of an infinitively precise estimate of the difference u χ 2 → 0 we are in the case of the Metropolis algorithm, when it is non 0, u χ 2 corrects, on average, for the effect of the noise.

In our situation, we have:

P(r, p) ∝ e -E(r,p) with E(r, p) estimated (3.51)
In order to implement the penalty method, we need a numerical estimator, ∆E p (p, p ; r), of the difference E(r, p ) -E(r, p) to be used in equation 3.51. This quantity will be discussed in the section 3.2.2. We now move to see what is required instead to sample a coordinate change in our Monte Carlo. In this case, indicating with T r (r → r ) and A r (r → r ) the probability to generate and accept a new configuration for the system (respectively), detailed balance becomes:

P(r, p)T r (r → r )A r (r → r ) = P(r , p)T r (r → r)A r (r → r) (3.52)
Using the explicit form of P and simplifying we have:

e -Vr(r) e -E(r,p) ρ m (r)T r (r → r )A r (r → r ) = e -Vr(r ) e -E(r ,p) ρ m (r )T r (r → r)A r (r → r) (3.53)
As for the p displacement, we have a probability density which is known within some numerical estimates. Nevertheless, the numerical complexity of the calculation of P is increasing in this case. Indeed, we have: r,p) with E(r, p) and ρ m (r) estimated (3.54)

P(r, p) ∝ ρ m (r)e -E(
Given the non exponential form of the probability in detailed balance, we cannot simply apply the Penalty method previously described but we have to combine this method with another method introduced by Kennedy et al. [START_REF] Kennedy | Noise without Noise: A New Monte Carlo Method[END_REF] who adapted Monte Carlo procedure for the following case:

P(s) ∝ f (s)e -L(
T (s → s ) ∝ e -L(s ) (3.56)
Detailed balance is simplified to:

f (s)A(s → s ) = f (s )A(s → s) (3.57)
This can be satisfied by defining the acceptance probability as:

A(s → s ) = c U(s → s ) if "s > s " c if "s ≤ s " (3.58)
where U(s → s ) is a numerical estimate of f (s ) f (s) . Note that equation 3.58 depends on the ratio f (s ) f (s) and a specific ordering criterion s > s (its particular form will be specified for our case) and guarantees that detailed balance is satisfied even in the presence of noise with a parameter to adjust the acceptance probability between 0 and 1. Indeed, c < 1 is a constant that ensures A(s → s ) ∈ [0, 1]. In our case the Kennedy method is implemented defining:

T r (r → r ) ∝ e -Vr(r ) e -E(r ,p) (3.59)
And the detailed balance is satisfied if the acceptance probability is of the following form and indicating the ordering criteria as "r > r " and "r ≤ r ":

A r (r → r ) = c U(r → r ) if "r > r " c if "r ≤ r " (3.60)
In the equation 3.60, U(r → r ) is an unbiased estimator of the ratio ρ m (r ) ρ m (r) . Numerical tests have shown that, for our calculations, we can set c = 0.9 (c is choosen in a way that we do not have often, around 1% maximum, A r (r → r ) > 1). The optimal choice of the ordering criteria, "r > r " and "r ≤ r " ", depends on the problem, here we adopted "r > r " equivalent to e -δ β V(r,∆r=0) > e -δ β V(r ,∆r=0) . When introducing the Kennedy procedure, I specified that e -L(s) is assumed to be known analytically. This is postulated in the original method to allow sampling of the transition probability either analitycally or via standard Monte Carlo. In our case, the condition is not met because e -E(r,p) is only known numerically. We can however solve this problem by using first the Penalty method to generate configurations according to T r (r → r ) ∝ e -Vr(r ) e -E(r ,p) . These configurations are thus generated using a Monte Carlo with transition probability t(r → r ) and acceptance probability a(r → r ): e -Vr(r) e -E(r,p) t(r → r )a(r → r ) = e -Vr(r ) e -E(r ,p) t(r → r)a(r → r)

(3.61)

t(r → r ) ∝ e -Vr(r ) (3.62)
e -E(r,p) a(r → r ) = e -E(r ,p) a(r → r)

(3.63) a(r → r ) = min 1, exp(-∆E r (r , r; p) -u χ 2 r ) (3.64)
where ∆E r (r , r; p) is a unbiased numerical estimate of E(r, p) -E(r , p) and u χ 2 r is defined in analogy with equations 3.49 and 3.50.

Here we have t(r → r ) ∝ e -Vr(r ) which means that we have to sample our transition probability along a spring chain in the r variables and so we have to use a Gaussian sampling and the staging variables [START_REF] Sprik | Staging: A sampling technique for the Monte Carlo evaluation of path integrals[END_REF]. I present the technical aspects of that in the Appendix B. In figure 3.1, the scheme of our Monte Carlo moves is summarized. In the next subsection we define the different numerical estimators ∆E r (r , r; p), U(r → r ) and ∆E p (p, p ; r) required by the algorithm.

Definition and evaluation of the numerical estimators

As it will become clear from equations 3.69, 3.72 and 3.75, in order to compute the unbiased numerical estimators ∆E r (r , r; p), U(r → r ) and ∆E p (p, p ; r), we need an auxiliary Monte Carlo procedure on the ∆r variables over the conditional ρ c (∆r|r).

In this procedure we satisfy the detailed balance when:

ρ c (∆r|r)T (∆r → ∆r )A(∆r → ∆r ) = ρ c (∆r |r)T (∆r → ∆r)A(∆r → ∆r) (3.65)
e -V ∆ (∆r) e -δ β V (r,∆r) T (∆r → ∆r )A(∆r → ∆r ) = e -V ∆ (∆r ) e -δ β V (r,∆r ) T (∆r → ∆r)A(∆r → ∆r) (3.66)

We generate the new configuration according to T (∆r → ∆r ) ∝ e -V ∆ (∆r ) . This can be done with the staging algorithm [START_REF] Sprik | Staging: A sampling technique for the Monte Carlo evaluation of path integrals[END_REF] (see Appendix B). The detailed balance can be now simplified to:

e -δ β V (r,∆r) A(∆r → ∆r ) = e -δ β V (r,∆r ) A(∆r → ∆r) (3.67)

In this case, we can use the standard Metropolis algorithm and we have an acceptance probability with the following expression:

A(∆r → ∆r ) = min 1, e -δ β V (r,∆r )-V (r,∆r) (3.68)
The need of this auxiliary Monte Carlo becomes clear if we look at the definition of U(r → r ). Indeed, U(r → r ) is defined as the ratio of the marginal probabilities:

ρ m (r ) ρ m (r) = d∆re -V ∆ (∆r) e -δ β V(r ,∆r) d∆re -V ∆ (∆r) e -δ β V(r,∆r) ρ m (r ) ρ m (r) = d∆re -V ∆ (∆r) e -δ β V(r,∆r) e -δ β V(r ,∆r)-V(r,∆r) d∆re -V ∆ (∆r) e -δ β V(r,∆r) ρ m (r ) ρ m (r) = e -δ β V(r ,∆r)-V(r,∆r) ρc(∆r|r) (3.69)
This quantity is sampled via the auxiliary Monte Carlo which allows us to calculate our unbiased estimator U(r → r ) as:

U(r → r ) = 1 N N k=1 e -δ β V(r ,∆r k )-V(r,∆r k ) (3.70)
where k = 1, • • • , N and N is the number of Monte Carlo step used to sample ∆r according to ρ c (∆r|r).

Let us now consider ∆E p (p, p ; r) = E(r, p ) -E(r, p). Since in all calculations discuss in the following it was sufficient to consider a second order cumulant expansion, we shall limit the discussion of the estimator to this case. Up to second order we have:

E(r, p) ≈ 1 2 C 2 (r)p 2 (3.71) with C 2 (r) = (∆r ν ) 2 ρc(∆r|r)
(3.72)

Structure of the PIM algorithm 67

The auxiliary Monte Carlo is required to compute C 2 (r) because it is a mean value over ρ c (∆r|r). Indeed:

C 2 (r) = 1 N N k=1 (∆r ν,k ) 2 (3.73)
and

∆E p (p , p; r) = 1 2 C 2 (r)(p 2 -p 2 ) (3.74)
where N is still the number of Monte Carlo steps used to sample ∆r.

It remains to specify the last estimator that we need for our algorithm, ∆E r (r , r; p) = E(r , p) -E(r, p). From its definition, this is given by:

∆E r (r , r; p) = 1 2 p 2 (C 2 (r ) -C 2 (r)) (3.75)
where C 2 (r ) and C 2 (r) are sampled from two different independant averages over the probability densities ρ c (∆r|r ) and ρ c (∆r|r) respectively. Nevertheless, if we proceed in this way the final variance is the sum of the two independent variances. This may severely reduce the acceptance probability in other cases. However, to avoid this issue, we can define an alternative estimator which will be presented in the Appendix A for the multidimensional calculations. Indeed, this is not an issue for all the calculations presented in this Chapter due to the fact that models for which we compute Wigner densities in the following are 1D and 2D model cases.

Nonetheless, rigorously for 2D model case, we have to refer to the multidimensional notations introduce in the Appendix A. This conclude our remarks on the algorithm section.

In the following, we will quantify the numerical cost of our tests by considering the total number of samplings computed as N T = N M × N , where N M is the number of Monte Carlo steps in the main cycle, i.e. the number of generated samples for (r, p). Convergence was obtained using N M of the order of 5 × 10 6 steps for 1D systems and 10 × 10 6 steps for 2D systems. N , the number of steps in the auxiliary Monte Carlo, was set to 500 steps except when noted otherwise. In all cases we have performed convergence tests on the number of beads, ν in equations 3.34. Depending on the potential and the temperature, convergence was achieved with ν equal to 16 beads (or slices) at 800 K, 33 beads at 300 K, and ν in the range 64 to 133 beads at 100 K. Sampling of the r and ∆r was performed using staging variables, setting the Monte Carlo parameters such that the acceptance ratio is about 40%.

Results

In this section we compare the Wigner density computed with PIM and with PIM+Edgeworth correction to benchmark results and to the results of calculations performed with alternative approximate schemes. With PIM, assuming that the cumulant expansion converges, the only approximation we introduce is the truncation of the series in equation 3.18 at second order. Similarly, in calculations with Edgeworth, a truncation of the series in equation 3.26 to order 4 is enforced, so we use equation 3.32 as estimator of the Wigner density. To obtain benchmarks, for one dimensional problems, we employed the Numerov method [START_REF] Nørsett | Solving ordinary differential equations I: Nonstiff problems[END_REF] to solve the one-dimensional time-independent Schrödinger equation. The calculated eigenfunctions and eigenenergies were then used to compute the density matrix and its Wigner transform on a two dimensional regular grid. For systems with two degrees of freedom, the eigenfunctions and eigenenergies were obtained using a numerically complete (256 basis set elements) expansion on the 2D harmonic oscillator eigenfunctions and a Gauss-Hermite quadrature to evaluate the necessary integrals. Also in this case the density matrix was then computed on a regular grid and its Wigner transform obtained using a fast-Fourier transform procedure.

Harmonic oscillator

We begin by testing the PIM algorithm on a simple one dimensional harmonic oscillator of the form V (q) = 1 2 mω 2 q 2 , with m = 1 amu and ω = 3000 cm -1 . For this system, the method is exact (including truncation of the cumulant expansion) so the calculations we present are intended as a test of the algorithm. The exact expression for the Wigner density

W ex (q, p) = tanh(f ) π e -tanh(f )[ mω q 2 + 1 mω p 2 ] (3.76)
is also available for this system (f = 1 2 β ω) and will be used for comparisons. We first try the method at the relatively high temperature T=1000 K. This is still not the classical limit for this oscillator, but at this temperature convergence of our estimators with respect to the number of beads is not hard simplifying equilibration and sampling. Indeed, we find that ν = 16 beads is enough to converge the result. In Figure 3.2, we show on the left/right panel the results for the Wigner density integrated over momenta/coordinate, so, for example,

W p (p) = dq W (q, p) (3.77)
(in the following, we shall refer to these integrated Wigner transforms as "marginal" densities). In the figure, the solid black curve is the exact result, while the solid green curve is the result of our calculation with N T ≈ 2×10 9 Monte Carlo steps. We

Results

69

also report, as the dashed black curve, the Boltzmann distributions of coordinates and momenta to show that we are far from the classical regime. As it should be, we are in excellent agreement with the exact results. The same quality of agreement is obtained, with the same number of Monte Carlo steps, at the considerably lower temperature T=100 K, as shown in Figure 3.3. In the figure, where the dashed curve is the classical result as before and the solid black curve is the exact result, we show convergence our estimator with the number of beads. Interestingly, W q (q) seems to converge faster than W p (p), requiring about half the number of beads. More in detail, for W q (q) we obtain results very close to the exact result indifferently with 64, 133 or 256 beads (number of beads for polymer in r so the total number of beads in the r and ∆r variables is 2×64, 2×133 and 2×256). Nevertheless, results for W p (p) show that we need to use 133 beads to have a result that matches the exact one, while 64 beads give a W p (p) that is slightly too narrow.

Morse potential

The next one dimensional system that we consider is the (anharmonic) Morse potential

V (q) = D × e -2αq -2e -αq + Dζ(q) (3.78)
with D = 20 Kcal mol -1 , α = 2.5 Å -1 , ζ(q) = 1 (q ≤ 2.5 Å), ζ(q) = e 10×(q-2.5) (q > 2.5 Å). This potential is a simple model for a bound proton with an additional term introduced to prevent dissociation. With the chosen parameters, the anharmonicity of the potential is quite strong and our method is no longer exact. We begin by considering a high temperature case in which the system is essentially in the classical limit. Figure 3.4 shows a comparison among our results (blue curve) and reference calculations (black curve). The two curves are superimposed for both "marginal" Wigner densities. We also show in red the results obtained with the QTB method. The W p (p) and W q (q) computed with QTB show a small, but visible, shift with respect to the exact results suggesting possible inaccuracies of the method also close to the classical limit. The problem with QTB becomes even more apparent if we consider the results shown in Figure 3.5, which refer to the same potential but at the lower temperature T=100 K. In this highly non classical case (see black dashed curve, Boltzmann distribution for the system), we find that the marginals computed via QTB, red curves in both panels, are close but appreciably different from the exact result, black curve. In particular, the shift in the position of the maximum of the coordinate distribution is now clearly visible. On the other hand, results obtained with PIM are closer to the exact ones both when the standard cumulant approximation is used (blue curve) and when the Edgeworth correction is included. Indeed, this correction is essentially irrelevant for the coordinate Wigner density, while it improves slightly the shape of the momentum distribution on the tails. The most striking improvement introduced by our method can, however, be appreciated considering the results shown in Figure 3.6 where we report the Wigner density computed as a function of the momenta for three different, fixed, values of the coordinates: q = 0.1 Å (data in green), q = 0.0 Å (data in red), q = -0.1 Å (data in black). In this figure, for each value of the coordinate, dashed lines indicate reference results while solid lines (of the same color) are the corresponding PIM results. Results for QTB are shown as filled squares. The reference results show a clear dependence between the p and q variables in the density. This dependence is captured by PIM, which is quite close to the exact result in all cases, but completely missed by QTB which displays a factorized form for the probability in coordinate and momenta. As mentioned in the previous section, we do not make this assumption and correlation, if it exists, is maintained at least approximately via the use of the cumulant expansion. This is the first example in which such correlation is significant.

Proton transfer model

We now test PIM on the following generic potential for symmetric/asymmetric A-H-B hydrogen bonds introduced by Marinica et al. [START_REF] Marinica | Generating approximate Wigner distributions using Gaussian phase packets propagation in imaginary time[END_REF]. Its form is a mixture of the H-bond potential introduced by Lippicott and Schroeder [START_REF] Lippincott | One dimensional model of the hydrogen bond[END_REF][START_REF] Reid | Semiempirical treatment of the hydrogen bond[END_REF] and the double-Morse model of Matsushita and Matsubara [START_REF] Matsushita | Note on isotope effect in hydrogen bonded crystals[END_REF] :

V (q, Q) = D e -2α( Q 2 +q-d) -2e -α( Q 2 +q-d) + 1 + Dξ 2 e -2 α ξ ( Q 2 -q-d) -2e -α ξ ( Q 2 -q-d) + Ae -BQ - C Q 6 (3.79)
The parameters in the potential are taken from Ref. [START_REF] Marinica | Generating approximate Wigner distributions using Gaussian phase packets propagation in imaginary time[END_REF]: D = 60 kcal mol -1 , α = 2.52 Å -1 , d = 0.95 Å, A = 2.32 × 10 5 kcal mol -1 , B = 3.15 Å -1 and C = 2.31 × 10 4 kcal mol -1 Å 6 , irrespective of the symmetric/asymmetric nature of the molecule. Within this model, q represents the distance of the proton from the center of the A-B bond and Q is the A-B distance. For ξ = 1, the proton potential is symmetric and the parameters chosen for the Morse potentials result in a potential without barrier for Q less than 2.45 Å and a proton frequency of 3750 cm -1 at infinite O-O distances. The other parameters are adjusted to have an equilibrium O-O distance of 2.4 Å, an energy of formation of 35 kcal.mol -1 , and a O-O force constant of 320 kcal.mol -1 .Å -2 , as for a H 5 O + 2 dimer. Such parameters are typical of O-H + -O proton transfer system. In Figure 3.8, left panel, we show the 2d contour plot of the potential (top part) and the shape of the potential energy for the transferring proton at a set of relevant values of the A-B distance (bottom part). In the bottom part of the figure, we see that at the equilibrium distance Q = 2.45 Å (red curve) the potential resembles a quartic interaction with no barrier to proton transfer, while at larger distances, e.g. Q = 2.75 Å (green curve), a barrier to proton transfer appears. At shorter distances the symmetric configuration, q = 0 Å, is a narrow minimum with quadratic character. A different situation can be represented by setting ξ = 0.707, all other parameters unchanged. In this case, the potential mimics an asymmetric O-H-N hydrogen bond with an equilibrium O-N distance of Q = 2.75 Å. The right panel of Figure 3.8 displays the contour plot (upper part of the figure) and the potential energy for the proton for fixed values of the A-B distance (lower part of the figure). It can be seen that in this model localisation of the proton on the oxygen site is favoured by about 15 kcal/mol. This simple two dimensional model of strong hydrogen bonds is highly anharmonic and manifests non trivial quantum effects as will be shown below. Let us consider first the symmetric case (ξ = 1) at fixed value of Q. Simulations were performed at 300 K with ν = 33 beads for two different values of Q. In the following we compare results obtained with PIM with reference quantum calculations, and results obtained with the QTB and GPSP methods [START_REF] Basire | Computing Wigner distributions and time correlation functions using the quantum thermal bath method: application to proton transfer spectroscopy[END_REF]. Classical Boltzmann distributions will also be shown for comparison. Figure 3.9 shows results for Q=2.45 Å. For the coordinate marginal distribution, QTB and GPSP are quite far from the reference. The QTB curve, in particular, has the same shape as the classical one with a rather flat distribution near the maximum which overemphasises the contribution of the flat region around the minimum in the almost quartic potential. In contrast, the curves obtained via PIM with the standard cumulant expansion (blue curve) and the Edgeworth correction (orange curve) are very close to the reference result and capture correctly the influence of the confining walls. The differences are less striking for the momentum distribution. However, in this case too, PIM is more accurate than the alternatives, and it can be observed that including the Edgeworth correction slightly improves agreement with exact results, especially on the tails.

For Q = 2.75 Å, tunnelling effects make an important contribution to proton transfer as can be seen from the much higher probability to find the proton at the top of the potential energy barrier in the reference calculation (solid black curve) as compared to the classical Boltzmann distribution in Figure 3.10. PIM is able to recover this tunnelling effect and performs better than QTB, which overestimates tunnelling, and GPSP, which underestimates it. It is interesting to observe that truncating the cumulant expansion at second order in PIM is a reasonably good approximation even when tunnelling plays a significant role and in regions of negative curvature such as those at the top of the barrier. Indeed, including the Edgeworth correction has a small effect on the PIM results, slightly improving agreement in the barrier region but essentially confirming convergence of the cumulant expansion. The effect of the correction is again more relevant on the momentum distribution function for the proton, Figure 3.10 (right), where it slightly improves agreement with the maximum and the tails of the distribution. In this case too, however, the main features of the overall distribution are captured quite well already by the second order cumulant expansion. This quality of agreement between the two calculations is a reliable numerical convergence test on the cumulant expansion. At this value of Q, we also performed calculations at the higher temperature T=800 K. As expected, Figure 3.11 shows even better agreement between PIM and Edgeworth results and the reference quantum calculations, with the disappearance of the small discrepancy in the tunnelling region that could be observed in Figure 3.10.

The impact of temperature on the performance of our method can be better appreciated from Figure 3.12 where we show the contour lines of the two dimensional Wigner density for this potential at T=300 K (left panel) and T= 800 K (right panel). At both temperatures, the classical (blue curve) and reference quantum results (black curves) are very different, indicating that, even at the higher temperature, the system is quite far from the classical limit. The exact densities show two main traits. The first is the "bow tie" feature extending from q ≈ -0.45 Å to q ≈ 0.45 Å , typical of the double well structure of the potential at this value of Q (see also the lower left panel of Figure 3.8). At T=300 K, our method fails to completely capture the deep tunnelling connecting the lobes of the distributions in the left and right potential wells, but the penetration in the barrier region is quite accurately captured. At the higher temperature, the agreement becomes very good and for these lobes the reference, PIM and PIM+Edgeworth results are essentially indistinguishable. The second characteristic feature of this double well potential are the contour lines that develop along the p axis. These are in fact negative regions of the Wigner density and, as expected, they are missed by the PIM calculations. On the other hand, the Edgeworth corrected results do reveal the presence of these regions. At the lowest temperature, this is reflected in the presence of some density in the qualitatively correct regions of the plot. At the highest temperature, however, the accuracy of the Edgeworth results improves considerably and both the shape and the location of the negative signal are better resolved in the region considered here. This can be further checked on Figure 3.13 where a cut at q = 0 along the Wigner distribution is shown. It can be seen that, even at this temperature, an oscillation is present due to tunnelling. The fourth order Edgeworth correction is able to reproduce both the position and depth of this oscillation, within the noise level. Around p = 0.09 g.mol -1 .Å.fs -1 , the exact Wigner distribution is again slightly positive; this can not be catched at fourth order correction and a higher degree of polynomial would be needed. This is however minimal at this temperature. It can be seen also, that the Edgeworth corection leads to a improvement of the Wigner distribution along this cut on the main peak height and width, and not just in the tails of the Wigner distribution.

As a final test, we investigated the full two dimensional problem both for the symmetric and asymmetric model of the hydrogen bond of Eq. 3.79 at the rather low temperature, T = 100 K. In Figure 3 significant role on the momenta distribution where it provides a result essentially indistinguishable from the exact. The effect of this correction can be appreciated more in detail by looking at the proton's marginal Wigner densities in Figure 3.15 where the improvement on the maximum and the tails of the momentum distribution is visible, while the coordinate distribution is already in perfect agreement with exact results when the second order cumulant expansion is used. The same can be seen on the heavier degree of freedom, Figure 3.16, where, however, a small shift of the coordinate distribution with respect to the exact result appears. It is also interesting to consider the marginal joint distributions for the coordinates and momenta of the light and heavy degree of freedom separately, see Figure 3.17 left and right panel, respectively. By construction, the Boltzmann densities show no correlation, while the Wigner densities do for both degrees of freedom. In particular, the momentum distribution for the light degree of freedom broadens when the coordinate deviates from the q = 0 Å equilibrium position. In this case too, the Edgeworth correction slightly improves the accuracy in regions of lower probability. More strikingly, for larger values of the momentum, the exact result shows negative parts of the marginal Wigner density for the light degree of freedom visible (these are the symmetric features for |p| > 0.1 g/mol. Å/fs in Figure 3.18). As in the previous example, the Edgeworth correction captures at least in part their presence and correctly identifies the onset of the negative regions of the Wigner density. However, this correction is not sufficient to fully account for the depth and the extent of the negative region so the overall shape of the contour level differs from the exact result. Also, in this low probability region, the numerical noise associated to the estimator of the Wigner density becomes quite large. This behaviour is more apparent in Figure 3.19 where we show a cut of the surface for q = 0 Å. The localisation and scale of the negative contributions (which are about two order of magnitudes smaller than the intensity of the main peak) can be appreciated from the inset, which also shows that the Edgeworth correction does produce a (very noisy with our statistics) negative signal in the correct location.

Let us now consider the asymmetric model (ξ = 0.707. In Figure 3.20 we show on the left the two dimensional marginal coordinate Wigner density. While the main quantum effect is broadening of the distribution as in the symmetric case, here the most striking difference with respect to the classical calculation is the larger cross correlation between the heavy and light degree of freedom reflected in the different orientation of the axis of the classical and quantum "ellipsoids". As the O-N distance decreases, the potential energy surface for the proton becomes more anharmonic and the quantum distribution is shifted to the right with respect to the minimum of the potential where the classical density remains more localized. This is similar to what we observed for the Morse potential calculations discussed above. PIM is once again able to account for this effect, giving a very reasonable agreement with the exact calculation. The difference between classical and quantum correlations among the two degrees of freedom in this system is even more striking if we consider the momenta. The right panel of Figure 3.20 shows the two dimensional momentum Wigner density. Both quantum and PIM curves exhibit a clear correlation between Red: exact, black: PIM, orange: PIM with Edgeworth correction and blue: classical Boltzmann distribution. For W q,Q (q, Q) the contour levels are separated by 5 Å -2 for the quantum distribution functions and 25 Å -2 for the Boltzmann distributions. For W p,P (p, P ), the levels are separated by 10 mol.g -1 .Å -2 .fs for the quantum distributions and 50 mol.g -1 .Å -2 .fs for the Boltzmann distribution. Exact results are shown in red, PIM with Edgeworth correction in orange. In the inset, the region with negative contributions is enlarged.

these momenta, with correlation coefficient equal to 0.34 in the PIM calculation and 0.27 in the reference calculation. Given the almost quantitative agreement of the contour lines shown in the plot, this difference in the calculated correlation coefficients must arise from regions of low probability where noise is more important. Aside from the existence of correlation between the momenta, the plot has some rather counterintuitive features. In particular, the sign of the correlation coefficient: classically the shape of the potential, see Figure 3.8, would suggest that as the O-N distance increases (positive velocity), the proton is displaced to the left moving as if it were bound to the oxygen (negative velocity) resulting in a negative correlation. This classical reasoning would then suggest that the correlation between momenta follows the shape of the two dimensional potential well (red part in Figure 3.8). However, a positive correlation is observed. This can be explained as follows: along the longest axis of the two dimensional potential well, the confinment is less than in the orthogonal direction such that Heisenberg inequality leads to a smaller momentum along the longest axis and a large momentum along the shortest axis of the potential well, hence the shape observed for the momentum correlation. The quantum nature of this effect is underlined also by the comparison with the Boltzmann distribution which, by construction, shows no correlation. PIM results, although noisier than in the symmetric case, are in remarkable agreement with the quantum curves. The second order cumulant matrix thus faithfully preserves the correlation between the momenta for this system. As observed for the symmetric case, the Edgeworth correction does not change significantly the results but (as in the two dimensional momentum Wigner density) noise on the tails can become significant. To complete the comparison we show also the individual degree of freedom Wigner densities in Figure 3.21. The overall agreement with exact results remains very good for both degrees of freedom, even if there is a slight decrease in accuracy compared to the symmetric case.

Results

85

-0.6 -0.5 -0.4 -0.

.20: Two dimensional probability distributions functions W q,Q (q, Q) (left) and W p,P (p, P ) (right) of the asymmetric model of proton transfer at T = 100 K. Red: exact, black: PIM, orange: PIM with Edgeworth correction and blue: classical Boltzmann distribution. For W q,Q (q, Q) the contour levels are separated by 5 Å -2 for the quantum distribution functions and 25 Å -2 for the Boltzmann distributions. For W p,P (p, P ), the levels are separated by 10 mol.g -1 .Å -2 .fs for the quantum distributions and 100 mol.g -1 .Å -2 .fs for the Boltzmann distribution.

-0.5 -0. The levels are separated by 10 mol.g -1 .Å -2 .fs for the quantum distributions and 100 mol.g -1 .Å -2 .fs for the Boltzmann distribution. (right) heavy degree of freedom, Q. The levels are separated by 10 mol.g -1 .Å -2 .fs for the quantum distributions and 20 mol.g -1 .Å -2 .fs for the Boltzmann distribution.

Conclusions

In this Chapter, we presented the Phase Integration Method (PIM) for the calculation of the Wigner density. As we have seen in Chapter 2, the Wigner density is a key quantity for methods based on the LPI representation of quantum correlation function. We showed that, using a second order cumulant expansion we can mitigate the effects of the phase factor present in the definition of the Wigner density. We also explored the possibility to obtain negative contributions to the Wigner density by modifying our sampling scheme so as to include the first non trivial term of the Edgeworth series.

We then introduced the algorithm used in PIM to compute the Wigner density. Indeed, the probability density, which has to be sampled from PIM in order to calculate the Wigner function, is not known analytically. Due to this standard Monte Carlo methods (such as the Metropolis scheme), which rely on the analytical knowledge of the density, cannot be directly applied. We described how to deal with this problem by combining two schemes for sampling "noisy" probability densities. These methods are the penalty [START_REF] Ceperley | The penalty method for random walks with uncertain energies[END_REF] and Kennedy [START_REF] Clark | Asymptotics of fixed point distributions for inexact Monte Carlo algorithms[END_REF] Monte Carlo algorithms. Finally, we applied PIM to calculations of the thermal Wigner density for a set of increasingly challenging low dimensional systems. PIM is able to capture quantum correlation effects among the different degrees of freedom (coordinate-momenta and momenta-momenta). Essentially all other available approaches (presented in this Chapter so QTB and GPSP) assume a factorized form of the momenta and coordinates reminiscent of the Boltzmann density. We have shown on a model for asymmetric A-H-B hydrogen bonds that this assumption can be violated in physically significant situations. The Edgeworth corrected approach is interesting since it allows to capture (at least in part) negative features of the Wigner density, a challenging if not impossible task for all alternative schemes. In this case, the effort necessary to capture the main (positive) features of the Wigner density is similar to that of a standard PIM calculation but the situation is more delicate for the negative parts of the density. These are usually much smaller than the main features, so the convergence of their Monte Carlo estimate is more problematic simply due to a worse signal to noise ratio. As a final comment, the applicability of the method discussed in this work depends on the convergence of the cumulant or Edgeworth series expansions employed. The convergence of these expansions can always be checked numerically thus ensuring control first of the validity and then of the accuracy of a given calculation. For low dimensional systems, this can be done by comparing the probability density (or an exhaustive set of marginal probabilities) obtained with truncation of the series at successive orders. For high dimensional systems, the convergence of average values of functions of the coordinate and momenta can be used as a test (as we will see in next Chapter). In Chapter 3, we have shown that PIM was a suitable method to compute Wigner densities. As we have explained in Chapter 2, this quantity plays a key role for methods based on a LPI representation of the quantum correlation function. It allows us to know if we have the correct set up of initial conditions before propagating them by classical dynamics. In this Chapter, we are interested in adapting PIM to the calculation of quantum time correlation function for computing infrared spectra. The Chapter is organised as follows. We first express the symmetrised correlation function for operators linear in postions. We then show, by direct comparison on a simple harmonic oscillator calculation with the Kubo correlation function, why the symmetrized one is not the most adapted correlation function for the calculation of infrared spectra. The next step will be then to describe the derivation to calculate the Kubo dipole-derivative correlation function for a point charge model of the dipole. Finally, in the Results section, we will compare our results, on relevant model calculations, to methods such as RPMD, CMD or LGA.

Symmetrised correlation function for Infrared spectroscopy

Let us begin by recalling the key aspects of the PIM method for calculating linearized quantum time-correlation functions. In its original formulation [START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF], PIM starts from the symmetrised (or Schofield) [START_REF] Schofield | Space-Time Correlation Function Formalism for Slow Neutron Scattering[END_REF] form of the quantum time correlation function of operators  and B, which is given by (see equation 2.67):

G AB (t; β) = 1 Z Tr Âe i Ĥt * c Be -i Ĥtc (4.1)
where Ĥ is the Hamiltonian, t c = tiβ 2 and Z = Tr e -β Ĥ is the partition function with β = 1 k B T , k B being the Bolzmann constant and T the temperature of the system. In Chapter 2, we showed that this function can be related in Fourier space to alternative quantum expressions of time correlation functions (see equations 2.68). In particular:

CAB (ω; β) = e β ω 2 GAB (ω; β) (4.2)
where CAB (ω; β) is the Fourier transform of the standard quantum correlation function (see equation 2.64). Furthermore,

KAB (ω; β) = sinh β ω 2 β ω GAB (ω; β) (4.3)
where KAB (ω; β) is the Fourier transform of the Kubo form of the time correlation function (see equation 2.66).

The relationships above ensure that physical observables, usually written in terms of the Fourier transform of the standard or Kubo correlation, can be obtained by computing the appropriate symmetrised correlation, taking a Fourier transform, and then multiplying by the appropriate factor. This is indeed the route employed before my work in calculations based on the PIM quasi-classical approximation [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. The linearized PIM form of equation 4.1 is obtained as follows. We introduce resolutions of identity in the coordinate basis to isolate matrix elements of the different operators in the equation 4.1 in the same manner as we did in the equations 2.105.

G AB (t, β) = 1 Z dr + 0 dr - 0 dr + 0 dr - 0 r + 0 | A|r - 0 × r - 0 |e it * c H/ |r - 0 r- 0 | B|r + 0 r+ 0 |e -itc H/ |r + 0 (4.4)
where, if N a is the number of atoms, r + 0 = (r

+(1) 0 , • • • , r +(3Na) 0
). With this notation, the upper index in the variables indicates a particular degree of freedom. Usually the operators A and B are straightforward to evaluate so the non trivial part in the previous equation is the product of complex time propagators:

K(r - 0 , r + 0 , r- 0 , r+ 0 ) = r - 0 |e it * c H/ |r - 0 r+ 0 |e -itc H/ |r + 0 (4.5)
Inserting two more resolutions of the identity a thermal and a dynamical part can be isolated as follows:

K(r - 0 , r + 0 , r- 0 , r+ 0 ) = dr + ν dr - ν r - 0 |e -β 2 H |r - ν r - ν |e -it H/ |r - 0 r+ 0 |e it H/ |r + ν r + ν |e -β 2 H |r + 0 (4.6)
The thermal part is calculated in the same way as we did for the Wigner density in Chapter 3 and the dynamical part derivation follows the steps described in Chapter 2 for the LPI approximation. Using the notation defined in the previous equation and writing the forward and backward path integrals in terms of mean and difference variables ( r 0 = (r - 0 + r + 0 )/2 and ∆r 0 = r + 0 -r - 0 for example), we end up with the following linearized expression: where:

r ν - ∆r ν 2 |e -it
U (r ν , p, r t , p t ) = n l=1 δ p l+1 -p l δ t -∇V (r ν+l-1 ) n+1 l=1 δ r ν+l -r ν+l-1 δ t - p l m (4.8)
where δ t = t/(n + 1), l = 1, • • • , n + 1 and r ν+n+1 = r0 = r t . We have also simplified the notation with p 1 = p and p n+1 = p t . The δ functions above collapse the mean path to a classical trajectory from phase space points (r ν , p) to (r t , p t ). Combining the thermal and dynamical calculations, we can show [START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF] that for operators A and B linear in positions the PIM symmetrised correlation function is simply:

G AB (t; β) = e -i p∆rν A(r 0 )B(r ν (t)) ρ e -i p∆rν ρ (4.9)
where ρ is identical (with multidimensional notations) to the probability density defined in equation 3.13. In the expression 4.9 a phase factor is still present. As we did for the Wigner density in the previous Chapter, we use a cumulant expansion to tame this phase factor. The steps of the derivation are identical to those in Chapter 2 and we finally obtain:

G AB (t; β) = A(r 0 )B(r ν (t)) P (4. 10
)
where P is the same probability density as in the equation 3.22 in a multidimensional point of view. Using the notations r = (r

(1) i , • • • , r (3Na) i
) (for i between 0 and ν) and p = (p (1) , • • • , p (3Na) ), we have:

P(r, p) = ρ m (r)e -E(p,r) dp dr ρ m (r)e -E(p,r) (4.11)
We present, in the Appendix A, the modification of the algorithm for the multidimensional case using the notations introduced above. As we can see in that Appendix, the structure of the algorithm to determine the initial conditions remains the same. The modifications occur on the evaluation of the estimators. We can see on Figure 4.1, the schematic sketch to compute correlation functions from PIM simulations. Indeed, compare to Chapter 3, we have to incorporate the dynamics in our algorithm. We perform a total of N M C Monte Carlo steps to sample initial conditions and N traj classical dynamics runs so we have N M C /N traj Monte Carlo steps between two trajectories. I present more in detail, in the Figure 4.2, the scheme of a single trajectory in the linearized path integral representation in PIM. For operators linear in positions, the symmetrized correlation function is obtained by taking the value over the first beads in the polymer in r for the time 0 quantity and then propagated classically the system with (r ν , p) as initial conditions. The sampling of the polymers in r and ∆r is done by the Monte Carlo algorithm presented in Chapter 3. It is also important to notice that the system is evolved classically only from the last bead of the polymer, the other beads have no interest from the dynamical point of view (but essential for the sampling).

Let us now consider more in detail the case of IR spectra. The experimental IR spectrum is given by the product of the absorption coefficient α(ω) and the refractive index n(ω) of the medium [72]. Using Fermi's Golden Rule, this quantity can be directly related to the Fourier transform of the symmetrised dipole autocorrelation function, see Appendix C:

n(ω)α(ω) = 4πω 3c sinh β ω 2 Gµµ (ω; β) (4.12)
For simple applications, we consider that we have a point charge model for our Chapter 4. Application of PIM to the Infrared spectroscopy dipole and so:

µ = Na i=1 q i r (i) (4.13)
where q i is the charge of the atom i, this leads to:

G µµ (t; β) = µ(r 0 ).µ(r ν (t)) P (4.14)
From equations 4.12 and 4.14, we can see two problems for the calculation of IR spectrum using the symmetrized PIM correlation function. The first one, discussed in detail in the next section, is the fact that we calculate a correlation between the first beads (A(r 0 )) for t = 0 and the last bead (B(r ν (t))) at time t. As we are going to see later, it requires a much higher statistic to converge the calculation in this case than in a case where we have a correlation between the last bead at time 0 and time t.

The second problem is more serious. Indeed, for IR calculations the prefactor sinh in the absorption coefficient in the symmetric case (see equation 4.12) is really hard to handle numerically. To see this, imagine a situation where we have to obtain two vibrational peaks with the same intensities at 1000 cm -1 and 3000 cm -1 (similar to water). These frequencies are the typical frequencies in IR calculations. The ratio for the prefactor for these two frequencies is equal to:

r = ω 2 sinh β ω 2 2 ω 1 sinh β ω 1 2 (4.15)
where ω 1 = 1000 cm -1 and ω 2 = 3000 cm -1 . At high temperature the problem is not very important but at 300 K, we obtain r ≈ 370 and this get even worse when we decrease the temperature with r ≈ 5.3 × 10 6 at 100 K. The statistic needed to converge the calculation of G µµ (t; β) in this case is unreachable due to the fact that we have to multiply the Fourier transform of G µµ by the sinh prefactor. So at low temperature, we can only obtain the positions of the peaks but not their intensities [9].

From the equations 2.68, the easiest version to compute seems to be the Kubo correlation function. Indeed, in this case, the IR spectra can be obtained from the dipole autocorrelation function and the dipole-derivative autocorrelation function, via:

n(ω)α(ω) = ω 2 βπ 3cV ε 0 Kµµ (ω; β) (4.16) n(ω)α(ω) = βπ 3cV ε 0 K μ μ(ω; β) (4.17)
In the dipole autocorrelation function K µµ (t; β) we only have to handle an ω 2 prefactor and, even better, there is no frequency dependent prefactor if we are able to calculate K μ μ(t; β).

In the next section, we will present how we can calculate PIM Kubo correlation functions which depend linearly of positions or momentum and how we can apply these to Infrared spectroscopy.

Kubo correlation functions for Infrared spectroscopy

We will derive in this section the PIM version of the Kubo correlation functions.

The first case will be for operators linear in positions (it is possible to do it for more general position operators but we will not use it). Indeed, this particular case will show us how the computational cost associated to the Kubo expression is reduced compared to the symmetrised case on the most simple test case: the position autocorrelation function of a 1D harmonic oscillator. Linear operators in the position, as we can see on the equations 4.16, are relevant for IR spectra and are appropriate for the isolated molecules considered here. Application of the methods to bulk systems in periodic boundary conditions is, however, simpler starting from the dipole-derivation expression in equation 4.17. I will then discuss how we can calculate Kubo correlation functions with operators linear in momentum.

Operators linear in positions

Derivation of position autocorrelation function

To simplify the notation and the demonstration, I will present the calculation of the position autocorrelation function. Then the generalisation to any operators linear in position is straightforward.

Let us consider the Kubo position autocorrelation function:

K rr (t) = 1 βZ β 0 dλ T r e -(β-λ) H re -λ H r(t) (4.18) (4.19)
with

r(t) = e it H/ re -it H/ (4.20)
To set the stage for the derivation, we begin by identically rewriting the expression above in a more convenient form. This is done by first using the change of variable λ = λβ/2 to obtain:

K rr (t; β) = 1 βZ β/2 -β/2 d λ T r e -(β/2-λ) H r e -(β/2+ λ) H r(t) (4.21)
and then reorganising the integral above as follows

K rr (t; β) = 1 βZ β/2 0 d λ T r e -(β/2-λ) H re -(β/2+ λ) H r(t) + 1 βZ 0 -β/2 d λ T r e -(β/2-λ) H re -(β/2+ λ) H r(t) = 1 βZ β/2 0 dλ T r r(t)e -(β/2-λ) H re -(β/2+λ) H + 1 βZ β/2 0 dλ T r r(t)e -(β/2+λ) H re -(β/2-λ) H ≡ 1 βZ β/2 0 dλ A λ + 1 βZ β/2 0 dλ B λ (4.22)
To go from the first to the second equality, we used the invariance under cyclic permutation of the trace to move the time evolved coordinate to the left of the product under trace and made the change of variable λ = -λ in the second integral. The third equality defines A λ and B λ so, for example,

A λ = T r r(t) e -(β/2-λ) H re -λ H e -β/2 H (4.23)
We are now ready to start the manipulations that will lead to the PIM expression of the correlation function. The strategy that we shall follow is to obtain an exact, computable, expression for the thermal part of A λ and B λ and then approximate the time evolution of the momentum operator on the left of the expression via classical trajectories. This second step corresponds to the linearisation approximation described in the Chapter 2 [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Sun | Semiclassical theory of electronically nonadiabatic dynamics: Results of a linearized approximation to the initial value representation[END_REF] and it is completed following the same steps described in equations 2.105 to 2.127 (and the same as we did for the symmetrized correlation function in this Chapter). We shall instead focus on the treatment of the thermal contributions since this is where the difference between the Kubo and the standard linearisation for the symmetrised form of the correlation function plays a significant role. Furthermore, since the operations to be performed on A λ and B λ are completely analogous, we shall limit ourselves to discuss them in detail for A λ . To proceed, let us introduce a coordinate representation for eq.(4.23) by inserting resolutions of the identity to conveniently isolate matrix elements of operators. Thus,

A λ = dr + ν dr - 0 dr - λ dr - ν r + ν | r(t) |r - ν r - ν | e -(β/2-λ) H |r - λ × r - λ | r |r - λ r - λ | e -λ H |r - 0 r - 0 | e -β/2 H |r + ν (4.24)
For further convenience, we can identically rewrite this as:

A λ = dr + ν dr + 0 dr - 0 dr - λ dr - ν r + ν | r(t) |r - ν r - ν | e -(β/2-λ) H |r - λ × r - λ | r |r - λ r - λ | e -λ H |r - 0 r - 0 | I |r + 0 r + 0 | e -β/2 H |r + ν (4. 25
)
where I is the identity operator. Since we consider an operator which depends only on the position, we can simplify the previous equation:

A λ = dr + ν dr + 0 dr - 0 dr - λ dr - ν r + ν | r(t) |r - ν r - ν | e -(β/2) H |r - 0 r - λ δ(r + 0 -r - 0 ) r + 0 | e -β/2 H |r + ν (4.26) C λ = C 2 ν dr ν • • • dr 0 d∆r 1 d∆r ν r ν + ∆r ν 2 r(t) r ν - ∆r ν 2 2r λ × exp - 1 2σ 2 r ν α=1 (r α -r α-1 ) 2 exp - σ 2 p 2 ν α=1 (∆r α -∆r α-1 ) 2 × exp - δ β 2 r ν + ∆r ν 2 ) + V (r ν - ∆r ν 2 ) exp - δ β 2 V (r 0 + ∆r 0 2 ) + V (r 0 - ∆r 0 2 ) × exp -δ β ν-1 α=1 V (r α + ∆r α 2 ) + V (r α - ∆r α 2 ) (4.32)
with σ 2 r = δ β /2m and σ 2 p = m/2 δ β , and ∆r 0 = 0 by definition. The relatively simple expression for C λ depends crucially on the fact that we are considering the position operator and is preserved for operators linear in the coordinates. This results, in particular, in the absence of terms in ∆r λ outside of the exponential which will prove convenient in the following (see equation 4.34). The complications arising from operators non linear in the coordinates are illustrated in the next Chapter. Substituting the equation above in the expression for the Kubo transform of the time correlation function, we then obtain (using the equations 3.34 to simplify the notation):

K rr (t; β) = C 2 ν βZ β 2 0 dλ C λ = C 2 ν βZ β 2 0 dλ d∆r dr r ν + ∆r ν 2 r(t) r ν - ∆r ν 2 2r λ e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (4.33)
This expression can be further simplified by evaluating the integral over λ. This can be done by discretising the integral consistently with the discretization of the path (i.e. using a discretization step ∆λ = δ β = β 2ν ), as follows :

K rr (t; β) = 1 βZ β 2 0 dλ C λ ≈ ∆λ βZ ν j=0 C j (4.34)
This leads to:

K rr (t; β) ≈ C 2 ν ∆λ βZ ν j=0 d∆r dr r ν + ∆r ν 2 r(t) r ν - ∆r ν 2
× 2r j e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (4.35)

≈ C 2 ν νZ ν j=0 d∆r dr r ν + ∆r ν 2 r(t) r ν - ∆r ν 2 × r j e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r)
where we used the expression of ∆λ to simplify the prefactor multiplying the integral.

The expression that we have obtained for the Kubo correlation function has the nice property that the thermal propagators are now explicit functions of a finite set of variables. However, the real time propagation in the matrix element

r ν + ∆rν 2 r(t) r ν -∆rν 2
is still fully quantum and that makes the result incalculable. If we use as described in the Chapter 2 the linearisation approximation, however, we can write this matrix element as:

r ν + ∆r ν 2 r(t) r ν - ∆r ν 2 ≈ dp r ν (t)e -i p.∆rν (4.36)
where r ν (t) is the coordinate computed at time t along a classical trajectory with initial conditions (r ν , p). Using this result the quasi classical expression for the Kubo transform of the correlation function is obtained as:

K rr (t; β) ≈ C 2 ν νZ ν j=0
d∆r dr dp r ν (t)e -i p.∆rν r j e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (

The next steps of the calculation are similar in a multidimensional case to the calculation of the Wigner density. We did not specified them for the symmetrized correlation function in order to simplify the discussion in this section. To be clear and have a complete demonstration in the multidimensional case, we specify these steps in the following. All factors in the integrand above are now explicit and calculable. However, the presence of the phase factor e -i p.∆rν , which can cause very rapid oscillations in the integrand, hinders a brute force evaluation of this quantity. To obtain an expression more suitable for numerical purposes, it is convenient to manipulate further the quasi classical expression of the Kubo correlation function. We are going to proceed in the same way as we did for the Wigner density or the symmetrized correlation function, so we now introduce the joint probability density for the semisum and difference variables:

ρ(r, ∆r) = C 2 ν Z e -V ∆ (∆r) e -δ β V (r,∆r) e -Vr(r) (4.38)
and observe that this quantity can be expressed as ρ(r, ∆r) = ρ m (r)ρ c (∆r|r) where the marginal and conditional probabilities ρ m (r) and ρ c (∆r|r) are given by:

ρ m (r) = d∆r ρ(r, ∆r) = C 2 ν Z e -Vr(r) d∆r e -V ∆ (∆r) e -δ β V (r,∆r) (4.39)
and

ρ c (∆r|r) = e -V ∆ (∆r) e -δ β V (r,∆r)
d∆r e -V ∆ (∆r) e -δ β V (r,∆r) (

It is quite easy to show that we have also:

Z = C 2 ν d∆r dr dp ρ(r, ∆r)e -i p.∆rν (4.41)
With this notation, the quasi classical Kubo correlation can be written as:

K rr (t; β) = C 2 ν νZ ν j=0 d∆r dr dp ρ m (r)ρ c (∆r|r)e -i p.∆rν r j r ν (t) (4.42)
K rr (t; β) = C 2 ν νZ ν j=0
dr dp ρ m (r) d∆rρ c (∆r|r)e -i p.∆rν r j r ν (t) (4.43)

In the second line of the equation above, we have reorganised the expression so as to isolate, on the far right of the expression, the integrals over the ∆r variables. These integrals can be recognised as the generating function for the cumulant moments for ∆r ν on the conditional probability density ρ c (∆r|r) (similar for the Wigner density and symmetrized correlation function) [START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. Thus, d∆rρ c (∆r|r)e -i p.∆rν = e -E(r,p) (

where (in the multidimensional situation):

E(r, p) = +∞ |ne|≥2 (-1) |ne| 2 -1 n e ! |ne| p ne C ne (r) (4.45) with |n| = 3Na i=1 n i , n! = n 1 !...n 3Na !, and n e = {n if |n| even}, p n = p n 1 1 × ... × p n 3Na
3Na . C ne (r) are the cumulant moments of the random variables ∆r ν so:

C ne (r) = C n 1 ...n 3Na (r) = ∆r (1) ν n 1 ... ∆r (3Na) ν n 3Na c ρc(∆r|r) (4.46)
where the superscript c indicates a cumulant moment (see Appendix A). Substituting eq.(4.44) in eq.(4.43) we obtain:

K rr (t; β) = 1 ν ν j=0
dr dp ρ m (r)e -E(r,p) r j r ν (t) dr dp ρ m (r)e -E(r,p) (

As we have seen in Chapter 3, E(r, p) is a real quantity. We can then define a new probability density which the multidimensional identical to the one defined in equation 4.11 for the symmetrized correlation function:

P(r, p) = ρ m (r)e -E(p,r) dp dr ρ m (r)e -E(p,r) (4.48)
and interpret the expression for the Kubo transform of the correlation function as the following expectation value:

K rr (t; β) = 1 ν ν j=0 r j r ν (t) P (4.49)
We present, in the Figure 4.2, the scheme of the linearized path integral representation in PIM. For operators linear in positions, the Kubo correlation function is obtained by taking the mean value over all the beads in the polymer in r for the time 0 quantity and then propagating classically the system with (r ν , p) as initial conditions.

Application to the 1D harmonic oscillator

We are going to test the PIM algorithm for the position autocorrelation functions (symmetrised and Kubo) on a 1D harmonic oscillator (V (x) = 1 2 mω 2 0 x 2 ) and compare it to the exact result. For this system, PIM results must coincide with the exact one due to the fact that in this particular case we do not do any approximation (linearization being exact for quadratic potentials). Comparison of the symmetrised and Kubo correlation functions obtained via PIM will allow us to understand why the use of the Kubo version is recommanded for IR calculations. The exact symmetric position autocorrelation function for a 1D harmonic oscillator is:

G xx (t) = 2mω 0 sinh(β ω 0 /2) cos(ω 0 t) (4.50) (4.51)
The dynamics, in this case, being exact, we are going to focus in particular at value at time 0.

G xx (0) = 4mω 0 sinh(β ω 0 /2) (4.52)
with ω 0 the frequency of the harmonic oscillator. The PIM version (the one that we compute) of the symmetrized position autocorrelation function is simply (see equation 4.10):

G P IM xx (t) = x 0 x ν (t) P (4.53)
The exact Kubo position autocorrelation function is:

K xx (t) = k B T mω 2 0 cos(ω 0 t) (4.54)
K xx (0) = k B T mω 2 0 (4.55)
The PIM version of the Kubo correlation is as detailed before:

K P IM xx (t) = 1 ν ν i=0 x i x ν (t) P (4.56)
It would be convenient to express the exact Kubo correlation function as a sum over all the beads in order to have a direct comparison of all the x i x ν (0) to an exact result. We will not demonstrate in detail the following result but sketch the proof. First, we recall some basic results of the 1D harmonic oscillator:

H = ω 0 ( a † a + 1 2 ) (4.57)
x = 2mω 0 1/2 ( a † + a) (4.58)
H|n = ω 0 (n + 1 2 ) (4.59) x|n = 2mω 0 1/2 √ n|n -1 + √ n + 1|n + 1 (4.60)
with a † the creation operator, a the annihilation operator and |n the basis of eigenstates of the harmonic oscillator. Then calculating the trace in the Kubo position autocorrelation function in the basis of the harmonic oscillator, we can show that:

T r e -(β-λ) H xe -λ H x(0) = cosh ( β 2 -λ) ω 0 mω 0 e -β ω 0 (1 -e -β ω 0 ) 2 (4.61)
and knowing that for the harmonic oscillator the partition function as the following analytical result:

Z = 1 2 sinh β ω 0 2 (4.62)
We end up with the following expression (integration over λ):

K xx (0) = βmω 0 sinh β ω 0 2 β/2 0 dλ cosh (λ ω 0 ) (4.63) = 2mω 0 sinh β ω 0 2 1 ν ν i=0 cosh β ω 0 2ν i
As a consequence, PIM has to satisfy the following relation for each bead:

x i x ν (0) = 2mω 0 sinh β ω 0 2 cosh β ω 0 2ν i (4.64)
We can notice that we recover x 0 x ν (0) =

2mω 0 sinh β ω 0 2
. This is interesting because by taking only the value on the first bead we find the PIM result of the symmetrized correlation function.

To test this property, we performed calculations at different temperatures and fixing the value of ω 0 to be in the IR range for typical molecules. In particular, we chose ω 0 = 0.01 ua (ω 0 ≈ 2000 cm -1 ) and performed the calculations at the temperatures 1000 and 300 K. We did these simulations with a high statistic using 10 5 main Monte Carlo steps and 500 auxiliary Monte Carlo steps so a total of 5 × 10 7 Monte Carlo moves (very high for such a simple model). The convergence with the number of beads has been systematically checked. In Figure 4.3, we can see that at 1000 K only 8 beads are enough to obtain a converged result. Furthermore, we can notice that the PIM calculations are very close to the exact result for all the beads. However, we also observe that x 0 x ν (0) (so the symmetrized correlation function) and x ν x ν (0) show small discrepencies. This is highlighted in the table 4.1. The error on these two values are almost equal but the value of x 0 x ν (0) is 2.5 times smaller than x ν x ν (0) . So the relative error on this value is in proportion more important. Finally, if we look at the integration over all the beads in order to have K xx (0), we notice that the PIM result is very close to the exact one and so the error on the exact value is mitigated by the integration. This justifies the fact that it is more convenient to calculate the Kubo correlation function instead of the symmetrised one. The impact of averaging over the beads is even more striking if we look at the results at 300 K. Indeed, as we can see in the

PIM exact

x 0 x ν (0) = G xx (0) 1.13 × 10 -2 ± 5.10 -4 1.18 × 10 -2 x ν x ν (0) 2.96 × 10 -2 ± 5.10 -4 2.99 × 10 -2 K xx (0) 1.73 × 10 -2 ± 5.10 -4 1.74 × 10 -2
Table 4.1: Comparison of x 0 x ν (0) which correspond to the symmetrized correlation function at time 0, x ν x ν (0) and K xx (0) between the PIM results and the exact ones for the 1D harmonic oscillator at 1000 K.

PIM 32 beads PIM 64 beads exact x 0 x ν (0) = G xx (0) 5.51 × 10 -4 ± 3.10 -4 4.58 × 10 -4 ± 2.10 -4 2.84 × 10 -4

x ν x ν (0) 2.60 × 10 -3 ± 3.10 -4 2.78 × 10 -3 ± 2.10 -4 2.74 × 10 -3 K xx (0) 5.17 × 10 -3 ± 3.10 -4 5.35 × 10 -2 ± 2.10 -4 5.22 × 10 -3 Table 4.2: Comparison of x 0 x ν (0) which correspond to the symmetrized correlation function at time 0, x ν x ν (0) and K xx (0) between the PIM results (32 and 64 beads) and the exact ones for the 1D harmonic oscillator at 300 K.

small values of the bead index i is very hard to obtain due to the fact that these values are close to 0. Furthermore, we can see that if we want to have all beads (for high i) converged we need to use 64 beads. However, as we can see in table 4.2, if we want only to have a good value for K xx (0) 32 beads are enough even though the error is more important than with 64 beads. The most important thing is that if we look at the value of x 0 x ν (0) (which is the simmetrised correlation function), the relative error between the PIM results and the exact one are very important. Indeed, we have a factor 2 if we use 32 beads and a factor 1.5 with 64 beads. For the Kubo correlation function (integration over all the beads), the relative error is much smaller (less that 5%). The analysis performed in this section then demonstrates and rationalizes the numerical advantages of the Kubo over the symmetrized PIM expression of cor-relation functions depending linearly on the coordinates. The same benefits hold for operators depending on the momentum, as we shall illustrate in the following subsection.

Operators linear in momentum

In this section, we focus on obtaining the PIM approximation for the Kubo transform of the momentum autocorrelation function. As discussed at the end of the previous subsection, this is motivated by the specific application that we intend to pursue, i.e. computing infra red spectra. The relevant observable for these calculations, in fact, is the dipole derivative time correlation function which, for simple expressions of the dipole, can be related to the momenta time correlation function. The performance of the Kubo PIM algorithm will be demonstrated by applying it to obtain the spectra for two interesting model systems: the OH and CH 4 molecules. The anharmonic potentials adopted in these calculations have been used as benchmarks for other approximate schemes and have proved surprisingly challenging for some of them [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF]9,18,[START_REF] Liu | Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initiobased flexible and polarizable force field[END_REF]. Thus, in spite of the relatively low number of degrees of freedom, these systems represent a non trivial test set for our method.

As mentioned in Section 4.1, the experimental IR spectrum is given by the product of the absorption coefficient α(ω) and the refractive index n(ω) of the medium. This quantity can be directly related to the Fourier transform of the Kubo dipole-derivative autocorrelation function [START_REF] Liu | Quantum dynamical effects in liquid water: A semiclassical study on the diffusion and the infrared absorption spectrum[END_REF], K μ μ(ω; β), via the relation (see equation 4.17):

n(ω)α(ω) = 2βπ 3c K μ μ(ω; β) (4.65)
For systems in the gas phase, n(ω) ≈ 1. In the following, we are going to consider simple models for which a point charge description is adequate, so the dipole of the system can be written as µ = Na i=1 q i r i , where N a is the number of atoms and q i the charge of atom i. The dipole-derivative operator will then be a linear combination of the momenta of the particles:

μ = Na i=1 q i m i p i (4.66)
To simplify the notation, in this section we shall focus on obtaining the PIM expression for the momentum autocorrelation function of a system of many particles of equal mass. The generalization of this result to the dipole derivative case in then straightforward.

Let us then consider:

K pp (t; β) = 1 βZ β 0 dλ T r e -(β-λ) H pe -λ H p(t) (4.67)
To set the stage for the derivation, we proceed as we did for the position autocorrelation function. We begin by retracing the steps in equations 4.21 and 4.22 to obtain:

K pp (t; β) ≡ 1 βZ β/2 0 dλ A λ + 1 βZ β/2 0 dλ B λ (4.68)
Where now A λ , for example, is defined as:

A λ = T r p(t) e -(β/2-λ) H pe -λ H e -β/2 H (4.69)
The presence of the momentum operator between the thermal propagators introduces important differences with respect to the previous case (position operator).

To discuss them, let us begin by identically rewriting A λ as:

A λ = dr + ν dr - 0 dr - λ dr - λ dr - ν r + ν | p(t) |r - ν r - ν | e -(β/2-λ) H |r - λ × r - λ | p |r - λ r - λ | e -λ H |r - 0 r - 0 | e -β/2 H |r + ν (4.70)
where, as usual, the boldface notation indicates 3N a dimensional vectors, so, for example r + ν = r

+(1) ν , • • • , r +(3Na) ν
. The structure of the integrand above is shown in Figure 4.5 and can be read as a sequence of propagations in imaginary time (i.e. evaluation of matrix elements of operators of the form e -α H , with α real and positive) and evaluation of matrix elements of p. Starting from the bottom right corner of the triangle in Figure 4.5 and moving counterclockwise, we see in fact that computing the integrand requires to first obtain the matrix element of e -β/2 H , which amounts to computing the probability amplitude for the system to evolve from |r + ν to |r - 0 in an imaginary time equal to β/2. This is followed by two new legs of propagation in imaginary time, from |r - 0 to | r - λ in a time λ, and from | r - λ to |r - ν in a time β/2 -λ. These legs are interrupted by the evaluation of the matrix element of p between | r - λ and | r - λ . The momentum operator being local in the coordinate representation:

r - λ | p |r - λ = i δ(r - λ -r - λ ) ∂ ∂r - λ , (4.71)
however, the presence of this matrix element does not introduce a discontinuity in the propagation on the left side of the triangle in Figure 4.5 so that the total length of propagation in imaginary time along this side is also β/2. This property will prove advantageous in the following (see in particular the change of variables in eq.(4.74) that requires equal number of beads along the two thermal paths). The last operation to be performed in the integrand is the evaluation of the matrix element The imaginary time evolution on the right is direct, while the one on the left is obtained, starting from the top vertex, as the composition of segments of length λ and β/2 -λ. The evaluation of the matrix element of the momentum operator between r λ -and r λ -, the red segment, interrupts the evolution on the left, but does not change the total time. The real time evolution necessary to evaluate the matrix element of p(t) is indicated as the green horizontal line. of the (real) time evolved momentum operator between positions | r - ν and | r + ν , represented in the figure as the green segment. Let us focus first on obtaining a more convenient representation for the propagations and the evaluation of the matrix element shown in the blue and red segments of Figure 4.5. The first step is to introduce a path integral representation of the matrix elements of the imaginary time propagators to rewrite [START_REF] Feynman | Space-time approach to non-relativistic quantum mechanics[END_REF][START_REF] Kleinert | Path Integrals in Quantum Mechanics, Statistics, and Polymer Physics, and Financial Markets[END_REF]:

A λ = dr + 1 • • • dr - ν dr - 0 • • • dr - λ dr - λ • • • dr - ν r + ν | p(t) |r - ν e -m 2 2 δ β ν α=λ+1 (r - α -r - α-1 ) 2 × e -δ β ν α=λ+1 V (r - α-1 ) r - λ | p |r - λ e -m 2 2 δ β (r - λ -r - λ-1 ) 2 e -m 2 2 δ β λ-1 α=1 (r - α -r - α-1 ) 2 e -δ β λ α=+1 V (r - α-1 ) × e -m 2 2 δ β ν α=1 (r + α -r + α-1 ) 2 e -δ β ν α=1 V (r + α-1 ) (4.72)
with the boundary condition r + 0 = r - 0 , and where as always δ β = β 2ν with ν the number of beads. In the second line, we have isolated the dependence of the integrand on r - λ (the Gaussian term e

-m 2 2 δ β (r - λ -r - (λ-1) ) 2
) to prepare the next step, which is to evaluate the matrix element of the momentum operator in parenthesis. This is done by substituting eq.(4.71) in the integrand above, taking the derivative of e

-m 2 2 δ β (r - λ -r - (λ-1) ) 2
with respect to r - λ , and then using the δ(r - λ -r - λ ) to perform the integral over r - λ . Reorganising the terms in the exponents in parenthesis, these operations result in the following form for A λ :

A λ = dr + 1 • • • dr - ν dr - 0 • • • dr - λ • • • dr - ν r + ν | p(t) |r - ν e -m 2 2 δ β ν α=1 (r - α -r - α-1 ) 2 × e -δ β ν-1 α=1 V (r - α ) e - δ β 2 (V (r - 0 )+V (r - ν )) - im δ β (r - λ -r - λ-1 ) × e -m 2 2 δ β ν α=1 (r + α -r + α-1 ) 2 e -δ β ν-1 α=1 V (r + α ) e - δ β 2 V (r + ν ) (4.73)
Proceeding in analogy with the developments in the Phase Integration Method for the position autocorrelation function, we now introduce semisum and difference paths by performing the following change of variables according to equation 4.30:

r 0 = r - 0 (4.74) r j = r + j + r - j 2 ∆r j = r + j -r - j
for j = 1, ..., ν. In these new variables, A λ can be identically rewritten as (using the notation in equation 3.34):

A λ = - im δ β d∆r 1 • • • d∆r ν dr 0 • • • dr ν r ν + ∆r ν 2 | p(t) |r ν - ∆r ν 2
× e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (r λ -r λ-1 ) -1 2 (∆r λ -∆r λ-1 ) (4.75)

with r = (r 0 • • • r ν ) and ∆r = (∆r

1 • • • ∆r ν ).
The same set of operations leading to eq.( 4.75) can be performed on the term B λ = T r p(t)e -(β/2+λ) H pe -(β/2-λ) H in eq.( 4.68) to obtain:

B λ = im δ β d∆r 1 • • • d∆r ν dr 0 • • • dr ν r ν + ∆r ν 2 | p(t) |r ν - ∆r ν 2
× e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (r λ -r λ-1 ) + 1 2 (∆r λ -∆r λ-1 ) (4.76)

So that

A λ + B λ = im δ β d∆r dr r ν + ∆r ν 2 | p(t) |r ν - ∆r ν 2 × e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (∆r λ -∆r λ-1 ) (4.77)
Substituting the equation above in the expression for the Kubo transform of the time correlation function, last line of equation (4.68), we then obtain:

K pp (t; β) = 1 βZ β 2 0 dλ (A λ + B λ ) = im δ β 1 βZ β 2 0 dλ d∆r dr r ν + ∆r ν 2 | p(t) |r ν - ∆r ν 2
× e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (∆r λ -∆r λ-1 ) (4.78)

It is interesting to compare the equations 4.33 and 4.78, to note that the change of observable (coordinate vs momentum) has a big influence in the expression. Indeed, in 4.33, we do not have any more ∆ variables in the result contrary to 4.78. As we shall see later this case is less straightforward than the coordinate case. However, as with the coordinate PIM Kubo, this expression can be further simplified by evaluating the integral over λ. Using again the discretization step ∆λ = δ β = β 2ν , we obtain:

K pp (t; β) ≈ im δ β ∆λ βZ ν j=1 d∆r dr r ν + ∆r ν 2 | p(t) |r ν - ∆r ν 2
× e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) (∆r j -∆r j-1 ) . (4.79)

Observing that the exponents in the integrand are invariant with respect to the value of λ and that ν j=1 (∆r j -∆r j-1 ) = ∆r ν (we used ∆r 0 = 0), the Kubo correlation function can be written as:

K pp (t; β) ≈ 1 2νZ im δ β d∆r dr r ν + ∆r ν 2 | p(t) |r ν - ∆r ν 2
× e -δ β V (r,∆r) e -Vr(r) e -V ∆ (∆r) ∆r ν (4.80)

As we did previously, we are going to use the linearisation approximation (see [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Monteferrante | Linearized symmetrized quantum time correlation functions calculation via phase pre-averaging[END_REF] and Chapter 2 for detailed derivation) to simplify the time evolution in the previous equation to get:

r ν + ∆r ν 2 | p(t) |r ν - ∆r ν 2 ≈ dp p(t)e -i p.∆rν (4.81)
where p(t) is the momentum computed at time t along a classical trajectory with initial conditions (r ν , p). Using this result the quasi classical expression for the Kubo transform of the correlation function is obtained as:

K pp (t; β) = 1 2νZ im δ β
d∆r dr dp p(t)e -i p.∆rν × e -V ∆ (∆r) e -δ β V (r,∆r) e -Vr(r) ∆r ν (4.82)

All factors in the integrand above are now explicit and calculable. However, as discussed previously, the presence of the phase factor e -i p.∆rν , which can cause very rapid oscillations in the integrand, hinders a brute force evaluation of this quantity.

To obtain an expression more suitable for numerical purposes and an expression amenable, with relatively minor modifications, to the algorithm employed in our previous work, it is convenient to manipulate further the quasi classical expression of the Kubo correlation function. Let us begin by observing that the factors ∆r ν and e -i p.∆rν in the integrand can be recombined as:

∆r ν e -i p.∆rν = i ∂ ∂p e -i p.∆rν (4.83) so that:

K pp (t; β) = 1 2νZ - m δ β d∆r dr dp p(t) ∂ ∂p e -i p
.∆rν e -V ∆ (∆r) e -δ β V (r,∆r) e -Vr(r) (4.84)

The same marginal and conditional probabilities defined in equations 4. [START_REF] Liu | Test of the consistency of various linearized semiclassical initial value time correlation functions in application to inelastic neutron scattering from liquid para-hydrogen[END_REF]

The expectation value above can be computed via a the algorithm that has been described in detail in Chapter 3 and Appendix A. 1I present in detail in the Appendix D an alternative demonstration to obtain the result in the equation above. Furthermore, as we can see in the Appendix D, we can write formally:

K pp (0; β) = Na i=1 m i β p i (0) • ∂E(r, p) ∂p i P (4.90) and Na i=1 1 β p i (0) • ∂E(r, p) ∂p i P = 3 2 N a k B T (4.91)
where N a is the number of atom and m i the mass of the atom i. This relationship provides a formal connection between the PIM Kubo momentum autocorrelation function at time 0 and the temperature of our system. For all our simulations, we calculate K pp (0; β) and compare it to 3 2 N a k B T to check the convergence of our results with respect to the number of beads ν and truncation of the cumulant series.

To conclude this section, note that although we have derived our result for the position and momentum autocorrelation function, the method can be applied to any operator that can be expressed (or approximated) as a linear combination of powers of the coordinates or momenta. The path integral representation of the thermal propagators can in fact be repeated as long as the matrix element measured at imaginary time λ (the red segment in Figure 4.5) is local in the coordinate representation. This implies that the derivation of equation (4.88) can be repeated and the specific form of the operator affects only the form of the observable to be averaged. While this procedure is not difficult, the new form of the observable, which in the general case may depend on all r λ and ∆r λ , can complicate the numerical scheme and future work will focus on developing appropriate algorithms.

Application to the Infrared spectroscopy

In this section, we test the PIM approximation of the Kubo transform by computing three interesting benchmark infrared spectra. We will compare our results with exact ones (when available) but also with those of alternative quasi-classical schemes to asses relative accuracy and computational efficiency. To that end we have chosen two non harmonic model systems for which numerically exact (and experimental) results exist. These systems have also been adopted as non trivial test cases for alternative quasi classical methods such as RPMD and CMD. In all the calculations reported in this work, we have used the Eckart procedure [START_REF] Eckart | Some Studies Concerning Rotating Axes and Polyatomic Molecules[END_REF][START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF]9] (see Appendix E) to remove the effects of rigid translations and rotations on the system. For the first and second systems, OH -and a CH 4 molecules with quadratic type potentials for the bonds and angles (due to the presence of a non zero equilibrium distance, these potentials are not harmonic in Cartesian coordinates), results will be compare to the CMD and RPMD ones. We will see how and why these two methods fail even for simple systems. Then we will consider an OH molecule but with a strongly anharmonic Morse potential. The results for this molecule are going to be compared with spectra obtained via the LGA method.

OH -molecule

In the calculation considered here, the OH -molecule is described via the potential

V 1 (r) = 1 2 k b (r -r eq ) 2 (4.92)
where r = |r O -r H | and the values of the parameters, equal to those in [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF], are listed in Table 4.3. Formal point charges equal to q H = +1 and q O = -2 were assigned to the respective nuclear positions to compute the dipole of the molecule. With a non neutral molecule, the definition of the dipole depends on the origin of the reference system. We then employed the usual convention of placing the origin in the center of mass of the molecule.

• CMD and RMPD results

Due to the presence of a non zero equilibrium distance, the potential in 4.92 is not harmonic in Cartesian coordinates (in which all calculations discussed in the following will be performed). As a consequence, the prescriptions for CMD, RPMD and PIM are no longer exact and the potential is a non trivial test for these, and our, methods. Let us begin by summarizing the results obtained in [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF] which are shown, for several temperatures, in left (CMD results) and right panels (RPMD results) of Figure 4.7. The CMD spectra were computed at T=300 K (black curve), T=200 K (red curve), and T=100 K (blue curve). In the figure, the classical result at 300 K is also reported in green for comparison. The shape and location of the peaks show some striking features: as the temperature decreases, the peak shifts to lower frequencies, it broadens and the intensity diminishes. A discrepancy with respect to the classical limit is noticeable even at the highest temperature where the two calculations should converge to the same result. The exact spectrum, on the other hand, consists of a single peak at the physical frequency k b /µ ≈ 3730 cm -1 (µ is the reduced mass of the molecule) independent of the temperature. The difference between the CMD and exact results has been explained as follows [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF]. As mentioned in Chapter 2, in this method classical time evolution is imposed on the centroid to compute dynamical properties. The evolution is governed by an effective potential defined as

U ef f (x c ) = - 1 β ln ρ c (x c ) (4.93)
where ρ c (x c ) ≡ dx 1 ...dx P δ 1 P P i=1 x i -x c e -βU (x 1 ,...,x P ) is the centroid density. This effective potential depends on the temperature and for OH -it becomes strongly anharmonic when the temperature is lowered, with a shift of the minimum towards shorter distances. Therefore, the dynamics of the centroids is effectively a classical (Newtonian) dynamics in the effective potential U ef f (x c ). The consequences of this behavior have been analyzed by Witt et al. [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF] and they show that the quality of CMD results depends on the centroid's displacement, which becomes larger the more curved in space the ring polymer is (which is the case when we decrease the temperature). This phenomenon is referred to as the "curvature problem". In terms of dynamics, the curvature problem means that the U ef f (x c ) becomes skewed towards the interaction center, leading to different (and artificial) dynamics. Indeed, the shape of the effective potential in Figure 4.6 clearly explains the unphysical features of the CMD spectra. Let us now consider the spectra calculated via RPMD recalling first that in this scheme each bead in the path integral representation of the system becomes a classical dynamical variable moving on the potential U (x 1 , ..., x P ) = 1 2 mω 2 P P i=1 (x i+1 -x i ) 2 + 1 P P i=1 V (x i ) as we described in the section on RPMD in the Chapter 2. The first term in this potential represents a harmonic chain connecting the beads in the ring polymer. A normal mode analysis (see Chapter 2) shows [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF] that the chain frequencies depend on the temperature and vary in the range ω m = 2π/β and ω M = 2P/β . The value of ω m for temperatures of typical experimental interest places it in the low to intermediate IR frequency range which means that spectra calculated from evolution on the effective potential will include peaks due to these (unphysical) modes in addition to those due to the physical potential. The presence of these spurious peaks is evident in the right panel of Figure 4.7 where we reproduce from [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF] the results at three different temperatures and in different resonance conditions among the chain mode frequencies and the physical frequency. In the top panel of the figure, we show results for the relatively high temperature T=350 K in a case where the physical and chain frequencies are off-resonance (in all panels of the figure, the dashed lines indicate the values of the chain frequencies). As expected, in addition to the physical peak at approximately 3730 cm -1 , another (lower intensity) peak is visible at about 2700 cm -1 , which is one of the chain frequencies with this set up. The agreement among RPMD and exact results worsens if one of the chain frequencies happens to coincide with the physical frequency. This is shown in the middle and bottom panels of the figure where the temperatures (T=109.2 K and T=436.5 K, respectively) were tuned so as to produce resonance. In both cases (one more quantum, the other in what should be the classical limit for the system) the spectra show several artifacts. At low temperature, the spurious peaks of low intensity originating from the chain mode are visible (this is the same phenomenon as in the upper panel of the figure), but in addition the physical peak at 3700 cm -1 is split in two. This splitting is amplified at the highest temperature (bottom panel), where there is also a broadeninig of the two features. On the other hand, at high temperature the intensity of the spurious peaks becomes essentially negligible on the scale of the intensity of the physical peak.

• PIM results

We performed calculations of the spectra with our method at three of the tem- peratures where either CMD or RPMD showed problems. These temperatures were chosen as T=300 K, T=100 K (to explore the range in which the red shift of the CMD peaks was evident) and T = 436 K, one of the resonant temperatures for RPMD. In all calculations the protocol followed was the same. First, convergence of the initial condition sampling with respect to ν, the number of beads in our path integral representation of the thermal quantum density, was established. In order to do this, we computed relevant observables (the Kubo kinetic energy and the zero time value of our correlation function) for increasing number of beads and verified convergence of the mean value (see below) by comparing with standard path integral estimates of the same quantity2 . The number of beads necessary to converge was ν(100) = 64, ν(300) = ν(436) = 32 (the number in parenthesis indicates the temperature of the calculation). The Eckart transformation was applied to this set of variables, and classical trajectories propagated from the transformed (r ν , p) i where i is the number of the trajectory. In all cases 5000 trajectories were enough to converge our spectra. The classical propagation was 10 ps long, as in the CMD and RPMD simulations, and was performed with a standard velocity Verlet algorithm, with forces obtained from the physical potential, with a time step dt=0.1 fs. The Fourier transform of the correlation function along each trajectory was computed and the final result was obtained by averaging these Fourier transforms. The spectra computed with this procedure are shown in Figure 4.8. In accordance with theoretical expectations, our spectra are all at the same -physical -frequency and, for this feature, there are no noticeable differences with the varying temperature. Indeed, the spectra are completely superimposed in the Figure 4.8. This is to be expected since none of the pathologies listed above for CMD or RPMD affects our calculations. In particular, since the positions and momenta of only one of the beads are propagated and the dynamics happens on the physical potential, none of the spurious features induced by the effective centroid potential or the harmonic part of the effective potential in RPMD can arise. In our calculations there is a computational overhead which arises from the auxiliary Montecarlo necessary to sample the noisy probability. In this case, we required a number of auxiliary Montecarlo steps equal to N aux = 250 to generate each initial condition. This means, in particular, that we required 250 calls to the potential function to generate each initial condition. Some computational overhead with respect to the declared number of trajectories is present also in the CMD and RPMD simulations, since in both cases initial conditions were obtained from uncorrelated phase space points collected along a thermostatted NVT trajectory for the system. It should also be pointed out that, for this model, also a purely classical calculation can determine the location of the peak in this system which is considered purely to test the algorithms.

CH 4 harmonic molecule

The next system that we consider is a point charge model for methane, with

V (r, θ) = bonds 1 2 k b (|r| -r eq ) 2 + angles 1 2 k a (θ -θ eq ) 2 (4.94)
where k b =0.30345 hartree/bohr 2 and r eq =1.09 for the C-H bond, k a =3.1068×10 -5 hartree/degree 2 and θ eq =107.8 degree for the H-C-H angle. These parameters have been chosen to reproduce the main features of the experimental IR spectrum, when the charges of q H =1.0e for the hydrogen and q C =-4.0e for the carbon are used.

• CMD and RMPD results

Results for CMD and RPMD are shown in Figure 4.9, left panel and right panel, respectively. The CMD spectra show again a shifting and broadening with decreasing temperature (the origin of this effect is the same as in the previous case). This effect is more pronounced for the high frequency peak due to the stretching modes, which red shifts by more than 200 cm -1 as the temperature goes from T=300 K (top panel) to T=100 K (bottom panel) [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF]. Also RPMD manifests several problems. In the top part of right panel , Witt et al. in [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF] show the power spectrum computed from the velocity autocorrelation function to identify both the IR active and inactive modes. Of course, the latter should not appear in the IR spectra computed from the dipole autocorrelation function and shown -for two different temperatures with different resonance conditions with the chain modes -in the lower panels. Contrary to this (correct) expectation, however both the middle and bottom panels show a resonance around 2700 cm -1 in correspondence with the highest IR inactive peak of the top panel. This is a new pathology with respect to the ones found for OH -. Furthermore, the stretching peak visible at 3000 cm -1 in the top panel splits and there is a spurious signal to its right (according to the analysis in [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF], this signal originates from a combination band involving one of the chain frequencies and the slowest stretching mode).

• PIM results

In the following, we show that the main features of the CH 4 spectrum can be reproduced using PIM adapted to the Kubo expression for the dipole derivative which is, in fact, again free from the problems apparent in CMD and RPMD. To that end, we performed calculations at two different temperatures, T=300 K (with ν = 64 in the path integral representation of the thermal propagators) and T=100 K (with ν = 128). At these temperatures, the system is quite far from the classical limit. 25000 classical trajectories of duration 10 ps were used to converge the spectrum at both temperatures. The results are shown in Figure 4.10 as the red and blue curves respectively. In the figure we also report the classical spectrum, black dashed line, to show the discrepancies with respect to the quasi classical calculations. We find that Figure 4.9: In left panel, IR spectra for CH 4 molecule obtained with CMD (with only the streching vibration); in right panel the result calculated using RPMD. Results are shown for several different temperatures. For centroid MD, as indicated in the insets, T goes from 300 K to T=100 K (lower panel). As for the RPMD we show the spectrum for T=400 K (as explained in ref. [START_REF] Witt | On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy[END_REF] from where the figure is reproduced, this spectrum is calculated via the velocity autocorrelation function to identify IR inactive modes, see that reference for discussion), in the middle T=136.5 K, at the bottom T=273 K. shows two main features: a peak resulting from the bending mode of the molecule at about 1300 cm -1 and a second peak due to stretching at about 3020 cm -1 . For the main features of the spectrum, the PIM method, along with the parametrization of the empirical force field, gives results in good agreement with the experiment, see first column in the table. The high frequency stretching peak calculated via PIM is shifted by only 15 cm -1 with respect to the experiment. The peak corresponding to the bending is less accurate, with a shift of about 100 cm -1 , but the agreement is still satisfactory. The experimental spectrum [75] also shows a smaller feature centred at 1550 cm -1 . This has been assigned as a combination band between an IR inactive deformation mode of the molecule and the bending centred around 1200 cm -1 . PIM is able to qualitatively capture also this rather subtle effect (the experiment shows one broad band, while our spectrum has a double band structure). As shown in the inset in Figure 4.10, the calculated spectrum at T=300 K also shows two low intensity features at about 2600 cm -1 and 2900 cm -1 (red curves in the inset). Similar peaks are present also in the experimental spectrum and have been assigned to an overtone of the bending mode centred around 1300 cm -1 and to a combination band due to the modes at 1500 cm -1 and 1300 cm -1 , respectively. In contrast, the classical calculation at the same temperature (black curve) barely shows signal at 1500 cm -1 , the peak around 2600 cm -1 is shifted and less intense compared to the PIM result, and there is essentially no signal at 2900 cm -1 . It seems then that the PIM method is able to capture these non-harmonic effects.

OH Morse molecule

The last system, that we consider for our comparison with alternative schemes, is a point charge model for the OH molecule. The potential is a shifted Morse given by and the masses are m H =1837 a.u (atomic units) and m O = 16 × 1823 a.u. The charges assigned to the atoms are q H =1.1128e for the hydrogen and q O =-1.1128e for the oxygen. This model has been proposed by Liu and Miller [START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF] as a standard benchmark for testing the accuracy of quasi classical methods. These authors used it, in particular, to test the so called LSC-IVR with LGA method. As discussed in Chapter 2, this is a scheme in which propagation of classical trajectories is combined with an approximate, ad hoc, method for approximate sampling of initial conditions from the Wigner transform of the thermal density [START_REF] Wigner | On the Quantum Correction For Thermodynamic Equilibrium[END_REF] (see Chapter 2).

V (r) = D e 1 -e -α(
In their paper [START_REF] Liu | Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initiobased flexible and polarizable force field[END_REF], they provide numerically exact results for the IR spectrum. We shall also use that reference for comparison with LSC-IVR, a state-of-the-art quasi classical approach. We performed calculations of the PIM Kubo IR spectrum significantly varying the temperature to verify the efficiency and accuracy of our result in different conditions. We begin by showing, in Figure 4.11, the result obtained at T=4800K. This calculation was performed as a consistency check: at this high temperature, the system is classical and indeed the results of a classical and PIM evaluation of the spectrum, the black and red curves respectively, are in excellent agreement. We now move away from the classical regime by lowering the temperature. In Figure 4.12, we show PIM spectra at T=800K (red curve), T=300K (green curve), and T=100K (blue curve). The number of beads used in the path integral representation of the thermal propagators, see eq. The position of the maximum of these spectra was estimated by direct inspection with an error approximately equal to 5 cm -1 , as inferred from the spectrum error bars.

The classical calculation is also reported as the black dashed curve. The results show that at these temperatures, the system is well in the non classical regime. The PIM spectra are essentially identical within noise and reproduce the quantum stability of the position of the maximum at the different temperatures. A comparison between PIM, classical, exact and LSC-IVR results is given in Table 4.5 where we present data for the position of the maximum of the IR peak at the two lowest temperatures. The first observation is that the classical calculation results in a considerable blue shift (about 160 cm -1 ) of the peak. The classical frequency is in fact quite close to that corresponding to an harmonic approximation of the minimum of the Morse, which gives ω e =3886.56 cm -1 . This result can be explained, by observing that the classical sampling of the initial conditions does not account for thermal quantum effects such as zero point energy. Consequently, the system has a very low energy and is confined to explore only the bottom of the Morse well

Kubo correlation functions for Infrared spectroscopy 123

where the harmonic approximation applies. The table also shows that PIM and LSC-IVR give very similar results, with a peak position of about ω qc = 3782 cm -1 in closer agreement with the exact result, ω ex ≈ 3700 cm -1 . This improvement is a direct consequence of the more accurate sampling of the quantum thermal density which allows to explore, at least in part, the anharmonic region of the potential. The residual blue shift of the peak compared to the exact result (about 80 cm -1 ) can be explained via the failure of the classical approximation to reproduce the long time behaviour of the quantum time correlation function. To verify this, Liu and Miller suggested to analyse the shift of the position of the peak as a function of the maximum time in which the correlation function is computed. They found that better agreement with the quantum result could be found for LSC-IVR if only the first few periods of the time correlation function were used to fit an oscillating function to use in the evaluation of the mean frequency. They argued that this was due to the fact that the classical approximation of the propagation becomes exact for short times. To compare with their result, we then computed the position of the spectrum at T=300K as a function of the number of periods included in the fit of the time correlation function. The results are summarised in table 4.6 and, in agreement with the observation of Liu and Miller [START_REF] Liu | Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initiobased flexible and polarizable force field[END_REF], we also find that a very good peak position can be obtained by evaluating the mean frequency via the fitting of an oscillating function to the short time dynamics. With only two periods we obtain essentially the correct quantum peak position, while an increasing blue shift is observed for longer calculations. This is in contrast with the classical simulation in which the frequency is completely stable (and incorrect) as a function of periods included in the fitting procedure. Thus, this section shows that there is essentially perfect agreement between our results and those of the LSC-IVR method, confirming that, within the limits of the quasi classical treatment, PIM is a reliable tool for computing IR spectra.

Let us conclude with some comments on the numerical cost of the method. The calculations discussed above require essentially the same effort necessary to converge LSC-IVR. Indeed, the results presented here were obtained by propagating N=50,000 trajectories, which is the same ensemble size used in ref. [START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF]. In our runs, we have imposed rather strict convergence criteria for the (numerically exact) sampling of the initial condition, leading to a relatively large number of steps, N aux =200, in the auxiliary Monte Carlo calculations related to the sampling of noisy probabilities. This overhead is however of the same order of, if not smaller than, the one necessary to decorrelate successive initial conditions using the Liu Miller sampling. Furthermore, the standard PIM approach (i.e. PIM calculations based on the symmetrised form of the quantum time correlation function) was recently applied to compute the dynamic structure factor for a realistic model of a condensed phase system in the quasi classical regime [START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. This calculation, in which the semiempirical Aziz potential was employed to simulate the dynamics of 64 Neon atoms, showed that PIM results match, both in accuracy and efficiency, those obtained with the Feynman-Kleinert linearised method proposed by Poulsen and Rossky [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF] and were in very good agreement with available experimental data. This comparison is relevant for two reasons. Firstly, the method by Poulsen and Rossky is the closest to PIM in terms of rigour and formal properties (for example, it is in principle able to capture quantum correlations among the momenta of the particles). Secondly, the initial condition sampling used in the standard and Kubo version of the PIM method is identical. Therefore, the good scaling of standard PIM with number of degrees of freedom observed for the Neon calculation will be preserved by the Kubo variation. As a final consideration, note that the computational cost of PIM is essentially dominated by the initial condition sampling. The time evolution of the system is in fact completely classical and involves only one of the beads introduced in the path integral representation of the thermal density. This implies, in particular, that the effort necessary to increase the length of the (real) time for which the correlation is computed is identical to that of a classical molecular dynamics calculation. The considerations above suggest that the scaling of the cost of our algorithm with number of degrees of freedom will be similar to that of other linearised approaches.

Conclusions

In this Chapter, we started by the derivation of the PIM symmetrized position autocorrelation function and showed why this function is not numerically adapted to the calculation of IR spectra. Then we extended the Phase Integration Method to compute the quasi classical approximation of the Kubo expression for the position autocorrelation function. This extension improves on the original formulation of the approach because it allows to access directly experimental quantities, often expressed in terms of the Fourier transform of the Kubo representation. In contrast, in the original approach, the method approximated the symmetrised form of the correlation function. To connect with experiments, it was then necessary to multiply the result by exponential factors often difficult to control numerically. In addition, on the calculation of the position autocorrelation function of a 1D harmonic oscillator, we have seen that at identical numerical cost the Kubo correlation function was converged contrary to the symmetrized one. Finally, we specialised the derivation to the case of operator linear in the momenta and coordinates and investigated the numerical performance of the approach in calculating the infra red spectrum of interesting benchmark systems. Comparison of our results with LSC-IVR show that, as for the original approach, "Kubo PIM" is as accurate as state-of-the-art quasi classical methods for low dimensional models. The calculation of the OH -and CH 4 spectra confirms not only that our approach does not suffer from the pathologies of popular schemes such as centroid and ring polymer molecular dynamics, but also that it can treat systems with a larger number of degrees of freedom reliably. The numerical cost of our calculations is very similar to that of alternative methods. Furthermore,the derivation presented here can be extended to the rather general case of operators that can be expressed (or approximated) as linear combinations of powers of coordinates or momenta. Future investigations will focus on using this observation to compute other observables of interest as we are going to see in the next Chapter. As a final technical point, to perform more efficiently the simulations just discussed and pave the way for higher dimensional calculations, during my PhD I developed and tested a parallel (using MPI) version of the PIM code. The first step is to communicate the input to all the processors. Then if we have to perform 50000 trajectories, we will run on the 96 processors (for example), 50000/96 trajectories on each one. Finally, they will communicate their results to do the final average. using basic quantum mechanics as a transmission probability T (E). On the reactant side, the rate constant corresponds to the states, initially present on the reactant side of the energetic surface, pondered by their probability of crossing the barrier. The exact quantum mechanical expression for the thermal rate coefficient of this "reaction" is [77]:

k(T ) = 1 2π Q r (T ) ∞ 0 dE e -βE N (E) (5.1)
Here Q r (T ) is the reactant partition function per unit length (it would become the partition function per unit volume for a reaction in three-dimensional space), which is given by statistical mechanics as:

Q r (T ) = 1 Λ(T ) = m 2πβ 2 1/2
(5.2) and N (E) is the "cumulative reaction probability" [21] which in the present case of a one dimensional barrier transmission problem is simply the barrier transmission probability, N (E) = |T (E)| 2 , so we obtain:

k(T ) = 1 2π Q r (T ) ∞ 0 dE e -βE |T (E)| 2
(5.3)

One can show that from this expression, we can derive a rate constant formula which involves quantum correlation function. Indeed, Miller et al. [21,77] provide an extremely elegant (and entirely rigorous) way to calculate the rate coefficient in equation 5.1 using flux-side correlation function. They show that the cumulative reaction probability N (E) can be calculated from the steady-state flux through any point (or "dividing surface") x = s as:

N (E) = 2π m p Ψ p | F |Ψ p (5.4)
where F is the flux operator for a dividing surface s defined as:

F (s) = 1 2m
[ p δ( xs) + δ( xs) p] (5.5)

After some algebra [21,77], we end up with the following flux-side correlation function expression for the rate constant: (5.20)

C f s (t) = i 2δ β 1 Q(s) dx - 0 • • • dx - ν dx + 0 • • • dx + ν x - ν | h(t)|x + ν δ(x + 0 -s)δ(x + 0 -x - 0 ) × x + 1 -x - 1 e -m 2
where we used the fact that x + 0 = x - 0 to simplify the equation. The next step has been presented many times and corresponds to the change of variables according to the equations 4.30. Using the simplified notations introduced in the equations 3.34, we can write:

C f s (t) = i 2δ β 1 Q(s) dxd∆x x ν - ∆x ν 2 h(t) x ν + ∆x ν 2 δ(x 0 -s)
× ∆x 1 e -Vx(x) e -V ∆ (∆x) e -δ β V (x,∆x) (5.21) where x = (x 0 , • • • , x ν ) and ∆x = (∆x 1 , • • • , ∆x ν ).

Then we proceed as we did for the Kubo correlation functions. We use the linearisation approximation which allows us to express our correlation function as:

C f s (t) = i 2δ β 1 Q(s)
dxd∆xdp h(t) δ(x 0 -s)∆x 1 e -Vx(x) e -V ∆ (∆x) e -δ β V (x,∆x) e -i p∆xν

(5.22)

Then using the equations 3.13, we can express our correlation function as:

C f s (t) = i 2δ β
dxdp h(t) δ(x 0 -s)ρ m (x; x 0 = s) d∆x ∆x 1 ρ c (∆x|x; x 0 = s)e -i p∆xν (5.23)

In this case, the cumulant expansion, previously introduced in Chapters 3 and 4, cannot be done in a straightforward way due to the fact that we have now ∆x 1 in the observable. We have to use a mathematical trick to express our correlation function in a more computable way. First, we define N as: (5.25)

N = d∆x
Then we do a joint cumulant expansion of the two independant variables ∆x 1 and ∆x ν , which leads us to:

N = d dε exp ∞ n 1 =0 ∞ n 2 =0 ε n 1 n 1 ! (-ip) n 2 n 2 n 2 ! (∆x 1 ) n 1 (∆x ν ) n 2 c
ρc(x 0 =s) ε=0

(5.26) with (n 1 , n 2 ) = (0, 0).

In the following we consider a cumulant expansion of order 2. The first thing to notice is that all the terms in ε with n 1 > 1 will be equal to zero because we take the derivative in ε = 0 (the derivative of N in this case will be proportional to ε n 1 -1 ). Given that ρ c is even in ∆x, we can observe that (n 1 , n 2 ) = (0, 1) and (n 1 , n 2 ) = (1, 0) give a 0 contribution. So the only couples of (n ≡ -ip ∆x 1 ∆x ν ρc e -E(x,p) (5.27)

where ≈ indicates the truncation to the second order of the cumulant expansion. After this treatment of the phase, we can now express our flux-side correlation function as:

C f s (t) = 1 2δ β 2
dxdp h(t) δ(x 0 -s)ρ m (x; x 0 = s)p ∆x 1 ∆x ν ρc e -E(x,p)

C f s (t) = ν β 2 p ∆x 1 ∆x ν ρc δ(x 0 -s) h(t) P(x,p;x 0 =s) (5.28) At this point it is necessary to discuss this expression. The first thing to notice is that in our simulation we have to constrain the first bead of our polymer. Then our result includes a correlation on the ∆ variables between the second beads (denoted with the label 1) and the last one (label ν). As discussed for the Kubo position autocorrelation function in Chapter 4, the convergence for the average of x 0 x ν is numerically difficult. We then expect difficulties with the observable in the equation 5.28. The relevence of this problem will be discussed in the following. where E ‡ 0 represents the height of the energetic barrier. The quantum (exact calculations) values were taken from [57].

At first, we point out the excellent agreement between our results and the quantum reference for temperatures above 300 K. Our objective of computing rate constants using the PIM method is therefore successful in this temperature range, and with very good precision up to 400 K. Nonetheless, at 300 K and at 200 K, we overestimate the exact results as we shown in more detail in table 5.2.

A possible explanation of the failure of PIM is based on the observations of the distributions of ∆x 1 as we can see in Figures 5.3,5.4. At low temperature, two maxima for the distribution of ∆x 1 when the first bead is fixed are visible which is not the case at high temperature where this distribution seems to be Gaussian. Furthermore, if we look at the histogram of the value of ∆x 1 ∆x ν (see Figure 5.5), we have a maximum for a value close to 0 (0.005) but we have also a local maximum around 0.3 (so ∆x 1 ∆x ν of the second maximum approximatively 70 times higher than the first one) at 300 K and only one maximum at 1000 K. The distribution at low temperature is more complicated to converge due to the presence of the two maxima. Failure to obtain this convergence leads to overestimating the rate 138 Chapter 5. Application of PIM to the calculation of rate constants The numerical cost for the calculations of our correlation function is the following. At high temperature (above 300 K in our case), we need 100000 trajectories of at least 10 fs to obtain the plateau of the correlation function which correspond for us to the limit t → ∞, as we can see on figure 5.6. Nevertheless, the difficult task in this calculation is the evaluation of ∆x 1 ∆x ν . As a consequence we have to use for each main Monte Carlo step, 5000 auxiliary Monte Carlo steps. Due to the problem of convergence at low temperature, even if we use 60000 auxiliary Monte Carlo steps (which is huge for such a simple system), the result seems to be not fully converged. These observations lead us to two possible situations. The first one is that the configuration where we have two maxima is very hard to handle numerically and we need to increase even more the statistics. The second one is more fundamental and corresponds to the fact that the ∆x distributions are far from being Gaussian. Indeed, in this case, the truncation of the cumulant expansion at the order 2 can be inadequate for such calculations.

Conclusion

In this Chapter, we have seen a rigourous derivation of the PIM symmetrised fluxside correlation function in order to calculate rate constants. We applied PIM to the symmetric Eckart barrier with parameters which allow us to do a comparison with classical and fully quantum simmulation (and which is a model to the linear H+H 2 reaction). The results are in very good agreement with the exact reference until a temperature of 300 K. Below this temperature we overestimate the rate constant. Our first impression is that it comes from numerical problems at low temperatures. Indeed, as we can see on Figures 5.3 and 5.5, we have to converge a quantity which is very sensitive numerically and far from being Gaussian (binodal distribution). In this case, the cumulant expansion at the order 2 might be not enough. We first tried to increase the statistic but this didn't give us more reliable results and the numerical cost was too high for such a system. The first possibility is to go the next order in the cumulant expansion, so the order 4. However, an Edgeworth expansion (see Chapter 3) is numerically harder in this case because we cannot include it so easily as we did for the Wigner distribution and in the case of the rate constant the effect of Edgeworth is mitigate between the difference of two terms. The other possibility is to include directly in the sampling the order 4 of the cumulant expansion but this is very hard to handle from a numerical point of view. Given the difficulties, future work will focus on using the Kubo form of the flux side correlation function to compute the rate constants to see if the numerical advantages demonstrated by it for IR spectra transfer also to this case. The starting point of the work discussed in this thesis is the path integral formalism. This was introduced in 1948 by Richard Feynman and it gives us a different formal way to look at quantum mechanics. The increasing power of computer made this formal idea applicable for the calculation of quantum properties. From a numerical point of view, the path integral formalism allows immediately to compute time-independent equilibrium properties. The calculation is done using the path integral molecular dynamics or Monte Carlo (PIMD or PIMC) methods based on the so-called classical isomorphism discussed in Section 2.1.2. Path integrals are not directly effective for real time propagation due to the dynamical sign problem. Due to their intuitive appeal, however, several path integral based methods have been developed to provide approximate information at a reasonable numerical cost. In Section 2.3, we first presented the Ring Polymer Molecular Dynamics (RPMD) and the Centroid Molecular Dynamics (CMD), two popular methods to calculate quantum correlation functions. However there is no formal justification for the approximate dynamics in these methods and so they present limitations which are difficult to predict. An alternative family of approaches, the so-called quasi-classical methods, are based on the Linearized Path Integral (LPI) representation of quantum correlation function. These schemes, discussed in Section 2.4, start from a path integral representation of the correlation function and approximate the quantum time evolution via classical trajectories, while preserving a fully quantum treatment for the thermal density. The methods introduced in Section 2.4, FK-LPI and LSC-IVR with a Local Gaussian Approximation (LGA), differ in the way in which the thermal density is represented. These state of the art methods are more closely related to the one developed in this thesis (Phase Integration Method) and introduced to compare them directly with PIM's performance. As we explained in Chapter 2, the Wigner density is a key quantity for the quasiclassical methods to obtain a set of initial conditions that are then evolved classically. In Chapter 3, we presented how this quantity can be calculated via PIM. The key step in PIM is the cumulant expansion to mitigate the effects of the phase factor present in the definition of the Wigner density. Sampling the probability density resulting from this step is non trivial since it can only be estimated numerically. We have described how to circumvent this problem by combining two schemes for sampling "noisy" probability densities. These methods are the penalty [START_REF] Ceperley | The penalty method for random walks with uncertain energies[END_REF] and Kennedy [START_REF] Clark | Asymptotics of fixed point distributions for inexact Monte Carlo algorithms[END_REF] Monte Carlo algorithms. Using a second order cumulant expansion, PIM is consistently able to reproduce the results of numerically converged reference quantum calculations even very far from the classical regime. We also explored the possibility to obtain negative contributions to the Wigner density by modifying our sampling scheme so as to include the first non-trivial term of the Edgeworth series. Our tests show that this method does provide improvements over current alternatives, most notably the presence of detectable signal in regions of negative Wigner density. Moreover, PIM is able to capture quantum correlation effects among the different degrees of freedom (coordinate-momenta and momenta-momenta). This is an improvement compare to the other approaches presented (QTB for example) which assume a factorized form of the momenta and coordinates reminiscent of the Boltzmann density. In particular, we have shown on a Morse potential and on a model for asymmetric A-H-B hydrogen bonds that this assumption can be violated in physically significant situations. In Chapter 4, we developed, tested and applied PIM for vibrational spectroscopy. As other quasi-classical methods [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF], PIM is based on combining MD and MC algorithms to compute quantum correlation functions. We first showed, on the very simple example of the position autocorrelation function of an harmonic oscillator, why it is more suitable for infrared spectra calculations to compute the Kubo correlation function instead of the symmetrized one. Then, we specialised the derivation to the case of operators linear in the momenta and investigated the numerical performance of the approach in calculating the infrared spectrum of interesting benchmark systems. Comparison of our results with LSC-IVR showed that, as for the original approach, "Kubo PIM" is as accurate and numerically efficient as state-of-the-art quasi classical methods for low dimensional models. Furthermore, the calculation of the OH and CH 4 spectra confirms not only that our approach does not suffers from the pathologies of popular schemes such as centroid and ring polymer molecular dynamics, but also that it can treat systems with a larger number of degrees of freedom reliably. The derivation presented here can be extended to the rather general case of operators such as the dipole-derivative correlation function but with a non linear dipole operator. Future investigations will focus on using this observation to compute other spectra of interest such as those for the water dimer [START_REF] Bowman | Vibrational analysis of an ice Ih model from 0 to 4000 cm(-1) using the ab initio WHBB potential energy surface[END_REF] or CH + 5 [START_REF] Joel | Quantum Deconstruction of the Infrared Spectrum of CH5+[END_REF][START_REF] Asvany | Understanding the infrared spectrum of bare CH5+[END_REF] with accurate potentials fitted on ab initio surfaces. In Chapter 5, we presented the methodology used to calculate rate constants with PIM. In this case, we have to calculate the symmetrized flux-side correlation function in which the operators are non linear in both coordinates and momenta. The characteristics of these operators have an influence on the way the cumulant expansion must be implemented. In order to compute rate constants, we have to separate the calculation in two parts. One part is a free energy calculation using an Umbrella sampling and the second is the calculation of the symmetrized flux-side correlation function using a constraint on a bead. We applied the algorithm to the chemical reaction rate for a linear model of the H+H 2 reaction using an Eckart potential, a known benchmark for such calculations. Results have been compared to classical and fully quantum simulations [START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF]. We successfully obtain the correct results on a broad range of temperature with accuracy similar to the LGA results which is the reference for quasi-classical methods on this system. We are, however, not able to reproduce the correct rate constants at low temperature (below 300K) even if PIM represents a definite improvement compared to purely classical calculations. In order to circumvent this limitation, future work will explore computing rate constants from the Kubo flux-side correlation function rather than the symmetrized one.

1 Figure 2 . 1 :

 121 Figure 2.1: Possible paths in the path sums of equation 2.23 using P = 5 slices.

1 Figure 2 . 2 :

 122 Figure 2.2: Possible paths in the discrete path sums for the canonical function in equation 2.30 using P = 5 slices.

Chapter 2 .Figure 2 . 3 :

 223 Figure 2.3: Intuitive picture of the classical isomorphism with P = 7.

1 Figure 2 . 4 :

 124 Figure 2.4: Evolution of the ring polymer with the temperature.

Figure 2 . 5 :

 25 Figure 2.5: Interaction pattern between two quantum particles represented by a ring polymer with P = 7.

Figure 2 . 6 :

 26 Figure 2.6: Quantum correction factor Q(u). Imaginary frequencies i|ω| are shown as -|ω| on the negative axis. Solid line: LHA and dashed line: LGA. The imaginary frequency for u i = β |ω| ≥ π is where the LHA's Q(u) becomes negative.

κ 2 =

 2 κ 2 (r) = ∆r 2 ν c ρc(∆r|r) , the second order cumulant (κ 2 (r) = ∆r 2 ν ρc(∆r|r)

  Thus equation 3.43 becomes e -E(p,r) P p acc (p → p ) = e -E(p ,r) P p acc (p → p) (3.44)

Figure 3 . 1 :

 31 Figure 3.1: Schematic representation of the Monte Carlo moves in the algorithm (for the multidimensional algorithm but it is identical for the 1D case): the penalty scheme for momenta moves and the combined penalty+Kennedy scheme for coordinate moves are represented on the left and right of the diagram, respectively. The numerical estimators which require auxiliary Monte Carlo on the ∆r variables are doubly underlined.

Figure 3 . 2 :Figure 3 . 3 :

 3233 Figure 3.2: (left) Integrated Wigner function W q (q) for a harmonic oscillator at T = 1000 K, computed exactly (plain line in black) and with PIM for different beads values. (right) Same for the integrated Wigner function W p (p).

Figure 3 . 4 :

 34 Figure 3.4: (left) Integrated Wigner function W q (q) for a Morse potential at T = 3000 K, computed exactly (plain line in black), with PIM (blue curve), and QTB (red curve). The insets zooms in the region near the maximum. (right) Same for the integrated Wigner function W p (p).

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: (left) Integrated Wigner function W q (q) for a Morse potential at T = 100 K, computed exactly (plain line in black), with PIM (blue curve), PIM+Edgeworth (orange curve), and QTB (red curve). (right) Same for the integrated Wigner function W p (p).

Figure 3 . 7 :

 37 Figure 3.7: Schematic representation of double Morse potential on a generic A-H-B configuration. Q represents the A-B distance and q the distance of the hydrogen from the center of the A-B bond.

Figure 3 . 8 :

 38 Figure 3.8: (left) Symmetric double Morse potential (ξ = 1) for different values of Q. (right) Asymmetric double Morse potential (ξ = 0.707) for different values of Q. The top panels show contour plots of the potentials as a function of the two coordinates. The horizontal lines on these plots correspond to the values of Q that were fixed in the calculations discussed in the text. The shape of the potential's cuts corresponding to those values of Q are shown in the bottom part of the Figure.

Figure 3 . 9 :

 39 Figure 3.9: (left) Integrated Wigner function W q (q) for a double Morse potential with Q=2.45 Å at T = 300 K (equilibrium distance), computed exactly (plain line in black) and with the different approximate schemes discussed in the text. (right) Same for the integrated Wigner function W p (p).

Figure 3 . 10 :

 310 Figure 3.10: (left) Integrated Wigner function W q (q) for a double Morse potential with Q=2.75 Å at T = 300 K. (right) Same for the integrated Wigner function W p (p).

Figure 3 .Figure 3 . 12 :Figure 3 . 13 :

 3312313 Figure 3.11: (left) Integrated Wigner function W q (q) for a double Morse potential with Q=2.75 Å at T = 800 K. (right) Same for the integrated Wigner function W p (p).

Figure 3 . 14 :

 314 Figure3.14: Two dimensional probability distributions functions W q,Q (q, Q) (left) and W p,P (p, P ) (right) of the symmetric model of proton transfer at T = 100 K. Red: exact, black: PIM, orange: PIM with Edgeworth correction and blue: classical Boltzmann distribution. For W q,Q (q, Q) the contour levels are separated by 5 Å -2 for the quantum distribution functions and 25 Å -2 for the Boltzmann distributions. For W p,P (p, P ), the levels are separated by 10 mol.g -1 .Å -2 .fs for the quantum distributions and 50 mol.g -1 .Å -2 .fs for the Boltzmann distribution.

Figure 3 . 15 :

 315 Figure 3.15: (left) Integrated Wigner function W q (q) for a symmetric 2D model of proton transfer at T = 100 K, computed exactly (black), with PIM (red) and PIM with Edgeworth correction (orange). The Boltzmann distribution is shown in dashed line. (right) Same for the integrated Wigner function W p (p).

Figure 3 .Figure 3 . 17 :Figure 3 . 18 :Figure 3 . 19 :

 3317318319 Figure 3.16: (left) Integrated Wigner function W Q (Q) for a symmetric 2D model of proton transfer at T = 100 K, computed exactly (black), with PIM (red) and PIM with Edgeworth correction (orange). The Boltzmann distribution is shown in dashed line. (right) Same for the integrated Wigner function W P (P).

Figure 3 . 21 :

 321 Figure 3.21: Partial Wigner distribution functions for the light and heavy degrees of freedom for the 2D asymmetric model of proton transfer at T = 100 K. Red: exact, black: PIM, orange: PIM with Edgeworth correction and blue: classical Boltzmann distribution. (left) light degree of freedom, q.The levels are separated by 10 mol.g -1 .Å -2 .fs for the quantum distributions and 100 mol.g -1 .Å -2 .fs for the Boltzmann distribution. (right) heavy degree of freedom, Q. The levels are separated by 10 mol.g -1 .Å -2 .fs for the quantum distributions and 20 mol.g -1 .Å -2 .fs for the Boltzmann distribution.

1 Figure 4 . 1 :

 141 Figure 4.1: Schematic representation of our algorithm for the calculation of correlation functions. It requires N M C Monte Carlo steps to sample initial conditions and N traj classical dynamics runs. Then an average is done to obtain our correlation function.

Figure 4 . 2 :

 42 Figure 4.2: Schematic representation of the polymers in r and ∆r used to represent the thermal propagators in our method and of the time propagation. The interactions between the thermal polymers are indicated in dashed lines and the red line represents the classical dynamic evolved from the initial conditions (r ν , p).

Figure 4 . 3 :

 43 Figure 4.3: Comparison of the Kubo PIM position autocorrelation function at 1000 K using 8 beads (red curve) and the exact result (black curve) at the same temperature.

Figure 4 . 4 :

 44 Figure 4.4: Comparison of the Kubo PIM position autocorrelation function at 300 K using 32 beads (green curve) or 64 beads (red curve) and the exact result (black curve) at the same temperature.

Figure 4 . 5 :

 45 Figure 4.5: Schematic representation of the sequence of propagations in imaginary time and evaluation of the matrix elements for the Kubo transformed momentum time correlation function. In blue segments of imaginary time propagation; note that the two sides of the triangle have equal total length, β/2 in imaginary time.The imaginary time evolution on the right is direct, while the one on the left is obtained, starting from the top vertex, as the composition of segments of length λ and β/2 -λ. The evaluation of the matrix element of the momentum operator between r λ -and r λ -, the red segment, interrupts the evolution on the left, but does not change the total time. The real time evolution necessary to evaluate the matrix element of p(t) is indicated as the green horizontal line.

Figure 4 . 6 :

 46 Figure 4.6: Effective potential for the OH -molecule computed via CMD. T=300 K (black curves), T=200 K (red curve, left panel only), T=100 K (blue curve) and we can see the anharmonicity of the effective potential when the temperature is decreasing.

Figure 4 . 7 :

 47 Figure 4.7: IR spectra for the OH -molecule computed via CMD, left panel , and RPMD, right panel. Results are shown for several different temperatures. In left panel, T=300 K (black curves), T=200 K (red curve, left panel only), T=100 K (blue curve); the classical result is shown in green for comparison. The number of beads, P , employed in the different temperature ranges is indicated in the inset. In right panel, T=350 K (top panel), T=109.2 K (middle panel), and T=436.5 K (bottom panel). The dashed lines indicate the chain frequencies.

4. 2 .Figure 4 . 8 :

 248 Figure 4.8: IR spectra of the 3D harmonic potential for the OH -molecule. Comparison of the Kubo-transformed PIM spectra at 100 (green curve),300 (red curve) and 436 K (black curve) spectrum.

Figure 4 . 10 :

 410 Figure 4.10: IR spectra of the 3D harmonic potential for the CH 4 molecule. Comparison of the Kubo-transformed PIM spectra at 100 and 300 K and the classical spectrum at 300 K. The inset shows a zoom of the region between 1500 and 3000 cm-1.

  (4.72), is ν = 32 at T=800K, ν = 64 at T=300K, and ν = 128 at T=100K. All calculations were converged by propagating 50000 classical trajectories of duration 10 ps. The number of auxiliary Monte Carlo steps was N aux =200. The convergence of the simulations has been checked with the number of beads and according to the equipartition of the energy (see equation 4.91). The inset of Figure4.11 displays the error bars on the PIM spectrum determined from a block average with 5 blocks of 10000 trajectories each.

Figure 4 . 11 :

 411 Figure 4.11: IR spectra of the 3D Morse potential for the OH molecule. Comparison of the Kubo-transformed PIM spectra at 4800 K and the classical spectrum at the same temperature.

2 e

 2 -δ β ν α=1 V (x - α ) e -δ β ν α=1 V (x + α )

  ∆x 1 ρ c (∆x|x; x 0 = s)e -i p∆xν ≡ ∆x 1 e -i p∆xν ρc(∆x|x;x 0 =s)(5.24) We can then introduce a variable ε and write N formally as: 134 Chapter 5. Application of PIM to the calculation of rate constants N = d dε e ε∆x 1 -i p∆xν ρc(∆x|x;x 0 =s) ε=0

Figure 5 . 1 :

 51 Figure 5.1: Free energy profile for an Eckart potential for different temperatures : 1000 K (blue curve), 500 K (red curve) and 300 K (orange curve). These free energy profiles are compare to the Eckart potential (black curve).

Figure 5 . 2 :

 52 Figure 5.2: Arrhenius plot : rate constants calculated using classical simulation (red curve), using a fully quantum code (black curve), and using the PIM flux-side correlation function (gree curve), as a function of 1000/T with T in Kelvins.

Figure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.3: Distributions of ∆x 1 with the first bead fixed at 300 K.

Figure 5 . 5 :Figure 5 . 6 :

 5556 Figure 5.5: Histograms of ∆x 1 ∆x ν at two different temperatures: 300 K on the left and 1000 K on the right.

  , Rev. Mod. Phys., 367, 20, 1948 Contents 2.1 Formalism of path integrals . . . . . . . . . . . . . . . . . . . 6 2.1.1 Path integrals and the canonical density matrix . . . . . . . . 7 2.1.2 Polymer isomorphism and PIMD . . . . . . . . . . . . . . . . 12 2.1.3 Time-independent equilibrium properties in the canonical ensemble from path integral formulation . . . . . . . . . . . . . 17 2.2 Dynamical properties from quantum correlation functions . 22 2.3 Empirical quasi-classical methods: RPMD and CMD . . . . 25 2.3.1 Outline of RPMD . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3.2 Outline of CMD . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.4 Linearized methods . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4.1 The Linearized Path Integral (LPI) representation of quantum correlation functions . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.2 The FK-LPI method . . . . . . . . . . . . . . . . . . . . . . . 38 2.4.3 The Linearized Semiclassical Initial Value Representation (LSC-IVR) method with a Local Gaussian Approximation (LGA) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

  Z P , and the ring polymer Hamiltonian H P (x, p). We use also the simplified notation d P p = dp 1 • • • dp P and the same for x. Finally, the functions A P (x) and B P (x t ) are averaged over the beads of the ring-polymer necklace at times 0 and t as it has already been explained for the determination of equilibrium properties:

							H B(t)	(2.78)
	K AB (t, β) ≈	1 (2π ) P Z P	d P p	d P x exp -	β P	H P (x, p) A P (x)B P (x t ) (2.79)
	where x t is the time evolved of x under the equations 2.77 and we keep the notations
	previously introduced for the canonical partition function with a finite number of
	beads P , A P (x) =	1 P	P k=1	A(x k )	(2.80)
			B P (x t ) =	1 P	P k=1	B(x k (t))	(2.81)
	The correlation function in equation 2.78 is just a classical correlation function in
	the extended phase space of the P -bead imaginary-time path integral, and this is
	the central feature of RPMD. A systematic derivation of the approximation in Equa-
	tion 2.78 does not exist. However, the RPMD approximation can be motivated or
	'partially' justified by showing that it gives the correct quantum mechanical results
	in various limiting cases.

  Structure of the PIM algorithm . . . . . . . . . . . . . . . . . 59 3.2.1 Structure of the algorithm . . . . . . . . . . . . . . . . . . . . 60 3.2.2 Definition and evaluation of the numerical estimators . . . . 64 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.3.1 Harmonic oscillator . . . . . . . . . . . . . . . . . . . . . . . . 68 3.3.2 Morse potential . . . . . . . . . . . . . . . . . . . . . . . . . . 70

	3.1.1 The PIM expression for the thermal Wigner density . . . . .	53
	3.1.2 Edgeworth expansion for the Wigner density . . . . . . . . .	57
	3.2

3.3.3 Proton transfer model . . . . . . . . . . . . . . . . . . . . . . 70 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

  s) with f (s) estimated and L(s) known analytically (3.55) Kennedy et al. showed that choosing the transition probability as:

  Operators linear in momentum . . . . . . . . . . . . . . . . . 105 4.2.3 Application to the Infrared spectroscopy . . . . . . . . . . . . 112 4.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

	Contents
	4.1 Symmetrised correlation function for Infrared spectroscopy 88
	4.2 Kubo correlation functions for Infrared spectroscopy . . . . 93
	4.2.1 Operators linear in positions . . . . . . . . . . . . . . . . . .	94
	4.2.2

  p t ∆r 0 e -i p∆rν U (r ν , p, r t , p t )

	H/ |r 0 -	∆r 0 2	r0 +	∆r 0 2	|e it H/ |r ν +	∆r ν 2	≈ dpdp t e
							(4.7)

i

Table 4 .

 4 3: Equilibrium structures and force constants for the O-H stretches and for the C-H stretches and the H-C-H angles in the OH -and CH4 model potentials.

	Type	r eq (Å)	k b (Ha/bohr 2 )
	O-H	1.00	0.49536
	C-H	1.09	0.30345
	Type θ eq (deg) k a (Ha/deg 2 )
	H-C-H	107.8	3.1068×10 -5

Table 4 .

 4 4: Peak positions of the IR active modes in the spectrum of CH 4 as obtained with PIM and exact calculations, compared to experimental results from ref.[75]. Units are cm -1 . the position and the shape of the main peaks is stable with respect to variations of the temperature. As reported in Table4.4, experimentally the IR spectrum of CH 4

		PIM Classical Exp.
	bending	1204	1195	1304
	deformation	1600		1550
	overtone of deformation 2600	2585	2601
	combination band	2900		2826
	stretch	3005	2995	3020

  1 , n 2 ) giving non zero term are (n 1 , n 2 ) = (0, 2) and (n 1 , n 2 ) = (1, 1) leading us to:

	N ≈	d dε	exp -	iεp ∆x 1 ∆x ν ρc -	p 2 2 2 (∆x ν ) 2	ρc	ε=0
	≡	-	ip ∆x 1 ∆x ν ρc exp -	iεp ∆x 1 ∆x ν ρc -	p 2 2 2 (∆x ν ) 2	ρc	ε=0

Table 5 .

 5 T (K) k quant (s -1 ) k class (s -1 ) k P IM (s -1 ) 3000 9.921 × 10 12 9.539 × 10 12 9.820 × 10 12 2000 3.697 × 10 12 3.423 × 10 12 3.610 × 10 12 1000 2.672 × 10 11 2.056 × 10 11 2.430 × 10 11 500 2.830 × 10 9 1.048 × 10 9 2.818 × 10 9 400 3.869 × 10 8 7.962 × 10 7 3.429 × 10 8 300 2.278 × 10 7 1.131 × 10 6 5.886 × 10 7 200 5.413 × 10 5 2.506 × 10 2 5.778 × 10 6 2: Value of the rate constant in s -1 for the symmetric Eckart barrier at different temperatures.

	5.4. Results

In actual calculations, the cumulant expansion in eq.(3.17) is truncated at the order

Numerical evidence shows that this ensures convergence of the series.

Testing against an independent calculation is necessary since, our probability density is obtained in part via a cumulant expansion truncated at second order. This truncation can only be validated via comparison to an independent calculation except when analytical results such as 4.91 are available.

Chapter 3

Wigner densities with the Phase Integration Method

Similarly, we have the following expression for B λ :

Given the structure of A λ and B λ , we can write

Following exactly the same steps described in Chapter 3, the matrix elements r + 0 | e -β/2 H |r + ν and r - ν | e -(β/2) H |r - 0 can be represented in path integral form, via a symmetric Trotter break up. In this way, we obtain an identical result as for the Wigner density:

where I kept the notation from the previous Chapter so

and the expression above becomes exact in the limit ν → ∞. Note that, to simplify the notation, we consider only particles withe same mass m. In the case where the mass are different, we have to associate the mass and the coordinates of each particle. Proceeding in analogy with the developments in the Phase Integration Method for the symmetrised form of the correlation function and, in general, linearization schemes (see [START_REF] Aage Poulsen | Practical evaluation of condensed phase quantum correlation functions: A Feynman-Kleinert variational linearized path integral method[END_REF][START_REF] Liu | Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions[END_REF] for example), we now introduce semisum and difference paths by performing the following change of variables [START_REF] Bonella | Path integral based calculations of symmetrized time correlation functions[END_REF][START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]17]: Comparison of the Kubo-transformed PIM spectra at 100, 300 and 800 K and the classical spectrum at 300 K. The structure towards the maximum of the position for the spectrum at T=100 K is most likely due to some residual noise. The black solid line is the exact numerical result reported in Table 1. The inset shows a zoom in the region from 3400 cm -1 to 4000 cm -1 of the PIM spectrum at 300 K, green curve. The curves in black are the plus and minus one σ standard deviation of this spectrum obtained from a block average (see text). Studying chemical reactions requires knowledge of the various rate constants that are involved in a reaction mechanism. A classical description can be adequate providing that the mechanism involves heavy atoms exchanges (or transfer) and/or that the reaction temperatures are high enough. In such cases, simple classical simulations can be reliable for computing the rate constants, using the transitionstate theory (TST) [START_REF] King | Development of transition-state theory[END_REF] for example. However, for systems expected to show a more quantum behaviour -encountered at lower temperatures and for light particles exchanges, such as electrons or protons -, classical simulation, which does not take into account the quantum nature of the phenomena (such as tunnelling effect or zero point energy), is inadequate. The derivation of a transition-state theory is in this case much more complicated (see [START_REF] Stuart | Derivation of a true quantum transitionstate theory[END_REF] for the derivation a quantum TST), and the numerical experiments required to investigate the rate constants values are more difficult. In this Chapter, we describe a scheme to calculate rate constants with PIM. The theoretical frame for this development -the model of the thermal rate constant and the expression obtained via PIM will be presented, and some preliminary test results will be discussed.

T (K) PIM LSC-IVR

General expression for chemical rate constant

The calculation of thermal rate constants, k(T ), for chemical reactions remains one of the central tasks of theoretical chemistry. The simplest thermodynamic model for a bimolecular chemical reaction would be crossing an energetic barrier (with the potential V (x) tending to zero as x ± ∞), separating the configuration space in two limit domains : the reactant side, for x → -∞ (on the left side of the barrier), and the product side, for x → +∞ (on the other side). The framework for our study requires a quantum treatment of the particles -the crossing of the barrier can hence also occur via tunnelling effect, for example. Let us call Ψ p the scattering wavefunction -composed of an incident, a reflected and a transmitted component -defined such that Ψ p (x → +∞) ∼ φ p (x)T (E), where φ p is a freeparticle wavefunction with momentum p and E its associated energy. The general probability for a scattering wave function to cross the barrier is then calculated 5.1. General expression for chemical rate constant

T r e -β H/2 F (s)e -β H/2 e -i Ht/ h( xs)e i Ht/ (5.6)

where h(s) is the Heaviside operator and h( x-s, t) ≡ e -i Ht/ h( x-s)e i Ht/ . The flux operator characterises the probability to go from the reactant side to the product side at t = 0 and the side operator allows to know if we are in the reactant or the product side. Formally, we have:

To simplify the notation, we denote h( xs, t) as h(t). It is also possible to derive expressions with the side-side correlation function or the flux-flux correlation function of the chemical rate constant [21,77]:

These two expressions are more difficult to handle numerically, due to the integration and the derivation, and thus, the flux-side is the most commonly used to calculate rate constants. Two important things have to be noticed. The first one is that formally the rate constant has to be evaluated from a symmetrised correlation function (it has been shown [START_REF] Ian | Quantum statistics and classical mechanics: real time correlation functions from ring polymer molecular dynamics[END_REF] that the standard and the Kubo correlation functions could also be used).

The second point is that our observables will be quite hard to compute. Indeed, F is an observable of momenta and coordinates and in addition non linear in coordinates. Furthermore, as we can see in the definition of the flux (see equation 5.5), we have to deal with δ-functions centered on the dividing surface (which for our 1D potential will be the top of the barrier). Then the probability for our system to be localized at this energetic maximum is extremely low. To extract a non-zero value for our correlation function, the simulations performed should therefore be very long to hope to crossings of the barrier, and inefficient from a numerical point of view. Moreover, for general multidimensional case defining the dividing surface is non-trivial. In the next Section, we will see how we can compute rate constants via the PIM procedure.

PIM Flux-side correlation function

For more convenience, we express the rate constant as:

where Q(s) is the partition function with the particle constrained at the dividing surface and thus, is defined as:

As we are going to see in the following, the numerical evaluation of

is simply the flux-side symmetrised correlation function with a constraint. The term Q(s) Qr(T ) is , as we will see in the next Section, a free energy term. We will present here how we can calculate this flux-side correlation function using PIM. We pose:

where we used the cyclic permutation property of the trace.

In analogy to what we did for the Kubo correlation functions, we introduce a coordinate representation by inserting resolutions of the identity to conveniently isolate matrix elements of operators. Thus:

The first step is to introduce a path integral representation of the matrix elements of the imaginary time propagators to rewrite (as we did for the Kubo correlation function):

(5.13)

Chapter 5. Application of PIM to the calculation of rate constants

Then we decompose the flux operator in two terms (see equation 5.5).

Considering, to isolate the relevant steps of the calculation, only the terms of the imaginary time propagators which depend of x - 0 and x + 0 (note that F is evaluated between thiese two coordinates), we have:

(5.15)

We recall that the momentum operator is local in the coordinate representation:

Given 5.17 we can express A in the coordinate basis:

Similarly we obtain the following result for B:

Then if we include the rest of the thermal paths and the prefactor 1/2m coming from the flux operator (see equation 5.5), we obtain for our symmetrised flux-side correlation function:

Free energy calculations

The second ratio involved in our calculations of rate constants is Q(s)/Q r (T ). We shall give this expression a more physical meaning. Firstly, and following the definition already used, we can rewrite the two partition functions as configuration integrals :

where x r designs the "reactant position", that is, when x → ∞ and we constrain ourselves in this reactant position. This expression corresponds, by definition, to the exponential of a free-energy difference. To say it otherwise, if we note A the free-energy, then:

where P designs a probability. We can hence rewrite as

With expression (5.31), we have transformed our problem to a free-energy profile calculation. Our goal is thus to express and compute the free-energy for our restrained (with one bead fixed) system.

I will describe the procedure used to calculate the free energy profile without enter in the detail since it is fairly standard. Indeed, for the simple 1D model that we are looking at, we simply use an Umbrella sampling method [78] combined with a standard Weighted Histogram Analysis Method (WHAM) [79].

In our path integral calculation method, the way to compute free energy using Umbrella sampling is the following: we add to the physical potential an harmonic potential (V umb = 1 2 k(x i -x min ) 2 ) on a chosen bead i (i is not relevant because in this case paths are cyclic) with a value for the harmonic spring k sufficiently high to keep the bead around the minimum of the harmonic potential. We perform this harmonic restraint for different value of x min in order to scan the physical potential. For every value of x min , we calculate the histogram of x i in this biased potential. Then we unbiase the histograms using the WHAM procedure which gives us our free energy profile. As we are going to see in the next section, if we have a potential with an energetic barrier, around the barrier we have to increase the value of k in the biased potential in order to be able to restrain x i around x min .

Results

In our model 1D system, the chemical reaction is modelled by a single particle crossing an energetic barrier (the particle was given the mass of a hydrogen atom).

Chapter 5. Application of PIM to the calculation of rate constants

The barrier chosen is the symmetric Eckart barrier, whose analytical form is

(5.32)

with parameters V 0 = 0.425 eV, a = 0.734 bohr and m = 1061 a.u chosen to mimic the H+H 2 reaction. This model was quite extensively studied by Miller and co-workers, which gives us numerous references to compare our results with (see [57,22,[START_REF] Liu | A simple model for the treatment of imaginary frequencies in chemical reaction rates and molecular liquids[END_REF]). We choose the diving surface to be on the top of the barrier so s = 0 a.u.

Free energy profiles

Free-energy profiles were computed using a 51-points grid between q = -4 and q = 4 atomic units for x min in the biased potential. Then we unbiased these 51 histograms using the WHAM code of Grossfield (see http://membrane.urmc.rochester.edu/content/wham). As output from this code, we obtain our free energy profile. In figure 5.1, we can see that at high temperature the free energy profile is close to the Eckart potential. However, already at 1000 K, the maximum of the free energy profile is below the top of the Eckart barrier. This behaviour becomes more and more important when we decrease the temperature and the free energy profiles become more different from the Eckart potential. This is due to the fact that at low temperature the tunneling effect becomes more important and lower the free energy. We can see this effect quantitatively in the table 5.1, where the maximum of the energy profile for different temperature is compared to the value of the Eckart potential at the top of the barrier.

Temperature (K)

Rate constants

Combining the free energy profiles with the PIM calculation of the flux-side correlation function (see equation 5.28), we have access to the rate constant. Figure 5.2 compares the performance of the PIM method to the quantum and classical results. The classical values were calculated through the usual transition-state

Modification of the algorithm for the multidimensional case

In the Chapter 3, we presented results for very low dimensional systems (1D or 2D). In this Appendix, we are going to look at how we have to adapt the algorithm presented in Chapter 3 for the multidimensional case.

The first difference is that before choosing if we move the coordinates or the momenta, we decide which atom we are going to move as we can see on the figure 3.1.

In fact the main difference is not in the structure of the algorithm but in the definition of the numerical estimators needed during the Monte Carlo moves to sample the density in the equation 4.11 (which is simply the mutidimensional version of the probability density written in the equation 3.22). Indeed, these changes are introduced to limit the numerical cost [START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. In multidimensional notation, the estimators are indicated as ∆E r (r , r; p) (see equation 3.75), U(r → r ) (see equation 3.69) and ∆E p (p, p ; r) (see equation 3.74). Here boldface indicates N a dimensional vectors, where N a is the number of atoms in the system, of vector components (r i = r i x , r i y , r i z ) for example (i = 1, • • • , N a ). The definition of the estimator U(r → r ) remains the same in multidimensional case. The first change with the method describes in Chapter 3 is in the evaluation of the cumulant C 2 (r). The definition becomes (truncating the cumulant expansion at the order 2):

where

If we move the momentum of the atom l (p = (p (1) , • • • , p (l) , • • • , p (Na) ) and p = (p (1)

We separate in the previous equation the terms involving the atom that is moved from the others. This allows us to not take into account many terms which simplify analytically to zero. C ij (which corresponds to an estimator of the cumulant) is equal to:

with ∆r

ν,α sampled from the conditional probability density ρ c (∆r|r) with the auxiliary Monte Carlo described in the Chapter 3. N is the number of auxiliary Monte Carlo moves. If we move the coordinates of the atom l (r = (r (1)

ν )), we formally have:

As explained in Chapter 3, the simplest way to proceed is to compute directly C ij (r ) and C ij (r) via two different independent averages (over the conditional probabilities ρ c (∆r|r ) and ρ c (∆r|r)). Nevertheless, if we proceed in this way the final variance is the sum of the two independent variances. This quantity is in general quite large and this severely reduces the acceptance probability. This was not a problem for the calculations described in Chapter 3 due to the low dimensionality of the systems studied. However, for the calculations of infrared spectra performed in Chapter 4, it becomes an issue due to the larger dimensionality of the systems. To avoid this issue, we can define an alternative estimator [START_REF] Monteferrante | Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method[END_REF]. To see how, let us start from the definition:

The trick in our development is to express ∆r

as a mean quantity over ρ c (∆r|r ) instead of ρ c (∆r|r):

The ratio of the marginals can turn be written as an expectation value:

= d∆r e -V ∆ (∆r) e -δ β V(r ,∆r) e -δ β V(r,∆r)-V(r ,∆r) d∆r e -V ∆ (∆r) e -δ β V(r ,∆r)

We can use the result above to express the equation A.10 as:

Furthermore, ∆x and ∆r are interchangable integration variables so:

d∆r e -V ∆ (∆r) e -δ β V(r ,∆r) d∆x e -V ∆ (∆x) e -δ β V(r,∆x)

d∆r e -V ∆ (∆r) e -δ β V(r ,∆r) e -δ β V(r,∆r)-V(r ,∆r) ∆r

So our estimator can be written as:

Now if we go back to the expression for the numerical estimator of equation A.7 (where we move only the coordinates of the particle l), we can write it as:

with

Although the expression of the new estimator is more complicated, it requires averaging over one distribution instead of two and this leads to a smaller variance than the direct one. Even though the algorithm works almost in the same way as for the 1D case, we can see that the multidimensional case is more technical.

Appendix B

Monte Carlo sampling of the polymer chains

B.1 Staging variables

In the Chapter 3, we needed to sample movements according to the equations 3.61 for the r and 3.65 for the ∆r displacements respectively. I recall here the relevantl quantities in the equations 3.61 and 3.65:

To improve efficiency in the sampling of these Gaussian densities, we can use staging variables [START_REF] Sprik | Staging: A sampling technique for the Monte Carlo evaluation of path integrals[END_REF][START_REF] Mark | Statistical Mechanics: Theory and Molecular Simulation[END_REF][START_REF] Brualla I Barbera | Path integral Monte Carlo Algorithms and applications to quantum fluids[END_REF]. In staging a collective move of ν -1 beads belonging to a fixed ends segment of the polymer. The anchors (first and last beads that are not moved) are choosen randomly. I will show an example of the staging variables on the r variables but the procedure remains the same for the ∆r variables. Define the end point for a segment of lenght ν as:

Then we can introduce the staging variables, denoted by u = (u 0 , • • • , u ν ), based on these two endpoints:

which gives us in the inverse definition:

Appendix B. Monte Carlo sampling of the polymer chains

Then we can apply this transformation to our transition probability and:

In this way we have decoupled the beads of the polymer chain. As a consequence, we can sample r k by using a simple Gaussian sampling on u k with a variance of σ r k k+1 and then use the inverse relation to obtain the new r k .

The Gaussian variables can be sampled via the Box-Muller method [83].

Appendix C

Infrared spectroscopy

C.1 Absorption coefficient from the Fermi Golden rule

H is the Hamiltonian of the system without any perturbation :

A time dependent perturbation of the following form is imposed:

We suppose that the perturbation

Where f and i represent the final and initial states of the transition. The transition rate is given by the Fermi golden rule:

where E f and E i are the energy of the final and the initial state of the transition. Then, we can express the time derivative adsorption energy of this transition according to the following expression:

where ρ i is the probability to be in the state i. This probability follows a Boltzmann distribution and we have:

Using the definition of delta functions, we obtain that:

We can use the fact that the matrix elements are symmetric in (i, f ) and so it is possible to switch the indices in the second term:

Appendix C. Infrared spectroscopy

We use again the delta function and the fact that we have a Boltzmann distribution:

Now if we replace the delta function by its integral representation, the absorption energy becomes:

In the Heisenberg representation : V (t) = e itE/ V e -itE/ and a| V (t)|b = e itEa/ a| V |b e -itE b / which allows us to write the absorption energy as:

where i ρ i i| V V (t)|i corresponds to the standard quantum correlation function (see equation 2.64):

Now if we consider that our perturbation is an electric field of the following form:

where ε represents the direction of the electric field and M the electric dipole is defined as follows in a point charge model:

This leads us using the mechanics describe before to the following expression:

If the absorbing system is isotropic, the polarization of the electric field is irrelevant and so:

C.1. Absorption coefficient from the Fermi Golden rule 155 and:

We can then introduce the spectral density as:

I stand (ω) means that it is the spectral density calculated from the Fourier transform of the "standard" quantum correlation function.

The definition of the Poynting vector S is:

Finally we obtain α the absorption coefficient via the relations:

We can also make the connection with the different definitions of quantum time correlation function with I K (ω) and I sym (ω) the Fourier transform of the Kubo correlation function and the symmetrized one respectively for the electric dipole.

Appendix D

Alternative demonstration for the Kubo momentum autocorrelation function

D.1 Notations of Kubo

In the next section, we will present an alternative demonstration for the momentum autocorrelation function obtained via PIM. However, in order to be clear and rigourous, we will introduce few notations. They come from the formalism of Kubo which has been introduced in his paper on the fluctuation-dissipation theorem [START_REF] Kubo | The fluctuation dissipation theorem[END_REF] and more particularly when he presents the linear response theory. The first notation is the Kubo correlation function:

where ρ = e -β H Z , β = 1 k B T and Z is the canonical partition function.

Then the notation of the time derivative expression of an quantum operator X is defined as follows:

Finally from the linear response developped by Kubo using these notations [START_REF] Kubo | The fluctuation dissipation theorem[END_REF], we have the following relation:

D.2 Momentum autocorrelation function

Using the previous notations, the relation D.

where N dof is the number of degrees of freedom of our system. Using the equation D.3, we obtain:

Then we can express [r i (0) ; p i (t)] in a more convenient way:

We used the permutation properties of the trace to go from the first equality to the second one. Then we will skip the details but using the linearisation approximation introduced in Chapter 2, we end up with the following expression:

where the W means Wigner transform and we have r

to have a correspondance with the notation used to present PIM. We can explicit [ρ; r (i) ν (0)] W . Indeed, we have: Here we only need to go to the order one because we have an operator linear in position only. So:

And we finally obtain:

where the PIM ρ W is nothing else that the Wigner density that we have derived in Chapter 3 [17]. So according to the equation 3.21, we have: The delta functions can be integrated out and in this case the PIM Wigner density can be expressed as:

If we combine equations D.6 and D.16, we obtain the following relation for momentum autocorrlation function:

Then we only need to take the derivation of ρ W along p (i) which gives us:

The same expression as the one in the equation 4.88 is recovered and a direct connexion with the equipartition of the energy is possible via this derivation.

Appendix D. Alternative demonstration for the Kubo momentum autocorrelation function

D.3 Equipartition of the energy

From equation D.21, we can obtain very simply K pp m :

To make the connection with the equipartition of the energy, we only need to have the time t = 0. Then:

If we come back to the equation D.6, we have:

which making the connexion with the PIM expression [17,18] gives us:

∂E(r, p) ∂p (i) p (i) (0)

The fact that from a time-independent quantity, we have access with PIM to the equipartition energy is very important because this will allow us to check the convergence of our simulation.

Appendix E

Eckart transformation

We have used the Eckart procedure [START_REF] Eckart | Some Studies Concerning Rotating Axes and Polyatomic Molecules[END_REF] to remove the effects of rigid translations and rotations on the system for our IR sepctra calculations. Power spectra obtained from MD simulations of finite molecular systems are dominated by the very broad zero frequency peak arising from the global translation and rotation of the molecule.

To isolate the relevant vibrational modes, the contribution of these modes must be eliminated from the signal. This is done via the so-called Eckhart transformation.

In this scheme, commonly applied in single trajectories calculations, the coordinate system is changed so as to employ a reference frame in which the total linear and angular momentum are zero. In our method, in which more than one trajectory is needed to compute the approximate correlation function, this change in coordinates is applied -after sampling of the initial conditions -to each member of the ensemble. The overall transformation is performed in two steps: first a change in variables to the centre of mass system is considered, and then the momenta are transformed to obtain zero total angular momentum. To describe these steps, let us introduce the following notation. The initial coordinates and momentum of atom i in the molecule are indicated as r (i) ν and p (i) . The vectors after the first and second transformation are denoted with a single and double prime, respectively. Let us begin by going to the reference frame of the centre of mass:

M tot (E.1)

where N a is the number of atoms in the molecule, m j the mass of the atom j and M tot = Na j=1 m j the total mass of the molecule. In this reference system, the total angular momentum (in general non-zero) is:

To obtain zero angular momentum, we perform the so-called rotation correction to the momenta:

Appendix E. Eckart transformation

In the above expression, the angular velocity ω is unknown and is determined by imposing the condition:

L is the new total angular momentum. Substituting equation E.4 in equation E.5 and using the properties of the vector product, the condition of zero angular momentum can be recast as:

where

ν ) T is the tensor of inertia of the molecule. The equation E.6 can now be solved for ω to obtain:

giving us the last ingredient to perform the Eckhart transformation.