Historiquement, le mouvement d'un réseau cristallin était traité du point de vue de la mécanique classique. Cependant, les effets quantiques nucléaires (NQE), à savoir l'énergie du point zéro et l'effet tunnel à travers la barrière de potentiel, peuvent modifier radicalement le comportement d'un cristal qui peut être analysé via diverses propriétés thermodynamiques. Ces propriétés peuvent être calculées exactement en utilisant des techniques Monte-Carlo d'intégrale de chemin (PIMC). Cependant, en utilisant cette approche pour déterminer des grandeurs dynamiques, on se heurte à des difficultés intrinsèques à la méthode. Dans le formalisme de Green et Kubo, qui est un outil commun pour l'analyse des propriétés de transport, les fonctions de réponse linéaires peuvent en principe être calculées par continuation analytique de la fonction de corrélation temporelle imaginaire obtenue avec PIMC. En pratique, cela correspond à une transformée de Laplace inverse, qui devient mal définie pour des données numériques souffrant d'une précision finie. Dans le présent travail, nous abordons ces questions et indiquons la voie possible pour le problème d'inversion qui convient à divers schémas de calcul. Nous démontrons l'utilité de cette approche sur plusieurs modèles simples de type oscillateur. A l'aide de cette machinerie, nous tournons notre attention vers l'analyse du NQE dans un cristal. Une description, qui tient correctement compte des phénomènes quantiques, revêt une importance pratique particulière pour les propriétés de transportȃ: contrairement aux métaux, le transport dans les semi-conducteurs et les isolants est régi par les vibrations de son réseau cristallin, qui sont décrites par les modes normaux. Dans le cristal réel, les modes interagissent et se dispersent les uns avec les autres, ce qui donne une durée de vie finie qui, en fin de compte, détermine la conductance thermique finie. Nous démontrons que la présence de NQE modifie le comportement en température de chaque mode ainsi que son taux de décroissance. Nous étudions également les changements correspondants de la conductance thermique par rapport aux prédictions classiques.

List of Figures

4.1

The dependence of the mean squared displacement x 2 on the imaginary time discretization ∆τ = β/M for a single quantum harmonic oscillator. Here we compare analytical results obtained with Eq.(4.5)(orange line) and Eq.(4.18)(green line) against the results of a PIMC calculation (blue dots). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Position correlation x(τ)x(0) as a function of imaginary time for a single quantum harmonic oscillator. Solid line corresponds to the diagonalization of the Hamiltonian (eq.(4.5)) and dashed line is given by eq.(4.18). Dots correspond to the MC simulations for different imaginary time discretizations β/M. Value of β is fixed and equal to 1 for all the data on the figure. Errorbars are smaller than the size of the dots. 4.3 MC results for the normalized correlation function(dots) as a function of τ plotted for several values of β. We compare them with the exact result(dashed line) which is obtained from Eq.(4.18) . . . . . . . . . . . 4.4 Spectral reconstruction of position correlation function (eq.4.18). Calculation was performed on a range of ω = [0, 2] with N ω = 50 δpeaks. We used 2 types of constraint for spectral density: "hard" cutoff A i ≥ 0 and "soft" cutoff A i ≥ -0.01 . . . . . . . . . . . . . . . . . . . . . 4.5 Comparison of different estimators for momentum correlation function for the same MC simulation of harmonic oscillator. Here we can see results for "naive" estimator eq.(4.29) (orange curve) and the calculation made with the estimator from eq.(3.32) written in the form (3.37) adopted to the momenta correlations (plotted in green). . . . . . 4.6 Comparison of relative errors between different estimators (naive and virial ones) of momentum correlation function. Figure (4.6a) shows the result for the whole range of imaginary time [0, β] and figure (4.6b) demonstrates the behaviour of the error for small times. With correlation function C(τ)

τ→β/2 ----→ 0 accuracy of the measurements becomes worse, most notably for the naive estimator. . . . . . . . . . . . . . . . . x 4.7 Reconstruction of spectra of momentum correlation function (eq. 4.21).

Calculation was performed on a range of ω (0,5) with discretization N ω = 50. We compare the results obtained with 2 different estimators: "naive" and "virial" ones. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Comparison of the reconstructed imaginary time correlation functions with the MC data used for the analysis. We plot the result as the deviation from the exact momenta correlation function (4.21). The errorbars are estimated from the analysis of the MC data. The reconstructed correlation functions are obtained from the spectra depicted on the Fig.(4.7) with the same color notations. . . . . . . . . . . . . . . 4.9 Relative discretization error, 1 -C PI pV (τ)/C exact pV (τ), between the path integral, C PI pV (τ), and the exact results, C exact pV (τ), for the energy current correlation function, as a function of hω 0 β/M. We show the data corresponding to the imaginary times τ = 0 and τ = β/2, and indicate with symbols and solid lines the results for hβ = 3 and 10, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Difference between the exact correlation function C exact pV (τ) and the values obtained by PIMC sampling, C MC pV (τ), of a path with M = 100 time slices, for hω 0 β = 1, illustrating the variance reduction obtained by the improved estimator discussed in Sect. 3.3.2. We show with line-points the primitive estimator and with the continuous line the improved estimator, both using the same path-integral Monte Carlo data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Difference between the exact correlation function C exact pV (τ) and the values obtained by PIMC sampling, C MC pV (τ), of a path with M = 30 time slices, for hω 0 β = 3, illustrating the variance reduction obtained by the improved estimator discussed in Sect. 3.3.2. We show in blue the calculations the 'primitive' estimator and in red -the improved estimator, both using the same path-integral Monte Carlo data. . . . . 4.12 Reconstruction of the spectral function associated to C pV (τ) at β = 10 corresponding to the indicated values for the number of delta functions in the model, N ω , and effective temperature Θ = 1. The area of the filled rectangles indicate the weight of the two delta-functions of the exact spectrum centered at ω 1 = 1 and ω 2 = 3, corresponding to the ∆ω = 1 discretization. As indicated in the text, we set ω 0 = 1. . . 4.13 Reconstructed spectra for the energy current correlation function C pV (τ) at β = 10, with N ω = 25 and at the indicated values of Θ. The filled rectangles are centered at the positions of the two delta-functions of the exact spectrum, with an area corresponding to their respective weights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 N ω -dependence of the χ 2 val extracted from the validation step of the reconstructed spectral functions for C pV (τ) at β = 10 and with Θ = 1. Squares and triangles correspond to shifted grids: for N ω = 5 red square shows the shift δω = 0.25 and the green one δω = 0.5, for N ω = 10 plot shows δω = 0.1 as a red triangle and δω = 0.25 as a green one. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . pV (τ) for the continuous distribution of oscillator frequencies, at the indicated values of the discretization N ω , at fixed Θ = 1. The shaded area indicates the exact spectral function. . 5.8 Spectral reconstruction of C cont pV (τ) for the continuous distribution of oscillators, for N ω = 10 and Θ = 1 and 6, respectively. Here we compare the results pertaining to a grid shifted by δω = 0.25 to those with δω = 0, the usual (not shifted) case. The shaded area indicates the exact spectral function. . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Main panel: χ 2 val from the validation procedure of spectral function at the corresponding values Θ and N ω for C cont pV (τ). The area of the circles is proportional to the value of χ 2 val . Purple circles correspond to the results for an ω-grid shifted by δω = ∆ω/2. Inset: χ 2 val as a function of the effective temperature Θ, at the indicated values of N ω .

6.1 Comparison [START_REF] Maitland | The intermolecular pair potential of argon[END_REF] of experimentally measured interaction between atoms of argon (solid line) and LJ potential(dashed line) (reprint from [START_REF] Allen | Computer simulation of liquids[END_REF]) . . . . . . . . . . . . . . 6.2 In nearest image convention one considers several copies of the original system and accounts only the interaction that corresponds to the smallest separation between particles. In this example only the interaction between green particle and blue one from the left-upper image of the system is accounted. . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Kinetic energy per particle e K of neon for several imaginary time discretization numbers M. The system is studied at T = 5K and ρσ 3 = 0.965. We compare the results of our simulation (red) with the data obtained from QMC simulation of analogous system in the work [START_REF] Cuccoli | Monte Carlo computations of the quantum kinetic energy of rare-gas solids[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Kinetic energy per particle e K of neon as a function of temperature T.

The values of density of crystal at each temperature are as indicated in Table 6.2. The red data indicate the present simulation, whilst the blue ones -reference data from paper [START_REF] Cuccoli | Monte Carlo computations of the quantum kinetic energy of rare-gas solids[END_REF]. The reference data were approximated in the limit ∆τ → 0. . . . . . . . . . . . . . . . 6.5 Pressure as a function of temperature for neon crystal described by quantum LJ potential. The figure shows the data for systems at indicated constant densities ρσ 3 . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Pressure as a function of temperature for neon crystal described by LJ potential. The figure shows the data for systems at constant densities ρσ 3 = 0.965 and ρσ 3 = 0.944. For comparison we plot results for the classical (squares) and quantum (triangles) systems. . . . . . . . . . . . 6.7 Example of a mode correlation function C AA (t) obtained from the MD calculation for the classical system with Lennard-Jones pair interaction. Various colors represent 12 different modes, all having the same degerate frequency. The black line shows the correlation that is obtained by averaging over all these degenarate modes. . . . . . . . . . . xiii 6.8 Example of a fit performed on the average mode correlation functions C AA (t). Figure shows the spectral density for the same mode calculated for the system at various temperatures. Different dots correspond to the spectra obtained by transforming C AA (t) whereas dashed line denotes the fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.9 Examples of various fit of a mode correlation function which corresponds to the lowest frequency in the LJ spectrum. The fit is performed with the stochastic reconstruction on a discretized range of frequencies, the different colours show few distinct effective inverse temperatures, whilst the number of δ-peaks is fixed N ω = 50. For the purpose of comparison we show the approximation of the spectrum with a gaussian distribution S(ω), which is obtained by a minimization routine. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.10 Comparison between different calculation of spectra. The dots show the stochastic spectral reconstruction with different effective inverse temperature Θ with ∆ω = 0.06. The blue line demonstrates the fit result obtained from minimizing the functional C AA iα (τ) -dωK(ω)S(ω, ω 0 , γ) with a Gaussian spectrum S(ω, ω 0 , γ). . . . . . . . . . . . . . . . . . . . 95 6.11 The comparison of mean frequency for two sets of modes corresponding to the {1, 0, 0} direction in k-space. On the left we show the transverse branch and on the right -longitudinal. Results of the Molecular dynamics calculation are plotted in blue with the linear approximation of temperature behaviour depicted as a dashed line. The linear regression is made with a fixed point at T = 0 corresponding to the mode frequency of a perfectly harmonic crystal. The output of the QMC simulation is showed in red. As the indication of insensitivity of the QMC data with respect to the change of the system size, we depict in green the result for the crystal of N = 864 particles at T * = 0.54.

In order to avoid overlapping with other data we shifted the point to the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.12 The total energy of the system per particle for the QMC simulation (data are shown in blue). The orange dashed line shows the polynomial fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.13 The comparison between various calculations of heat capacity for quantum neon crystal. The dashed line is calculated via the energy calculation [Fig.(6.12)], the blue and red dots are computed with the mode heat capacity given by expression (6.38). For the calculation of red points the frequency were used frequencies obtained from the QMC simulation whilst for the blue ones -the results of the MD computation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.14 The comparison between various calculations of heat capacity for quantum neon crystal. The dashed line shows the heat capacity calculated with the Debye law at Θ D = 65.8K, or in LJ units T * D = 1.87. The blue and red dots are computed with the mode heat capacity given by expression (6.38). As for the Fig. (6.13), for the calculation of red points the frequency were used frequencies obtained from the QMC simulation whilst for the blue ones -the results of the MD computation. For comparison we show the experimental data [START_REF] Fenichel | Low-temperature specific heats of solid neon and solid xenon[END_REF]] 100 xiv 6. 15 We demonstrate here the shift in mode frequencies ω k for the quantum Lennard-Jones crystal from the harmonic ones ω 0 k plotted against the ω 0 k . The different colors represent longitudinal and transverse branches. The data as well as the uncertainty are shown divided by the harmonic value of the same mode. . . . . . . . . . . . . . . . . . . . 6. 16 The figure shows the shift in mode frequencies ω k for the classical Lennard-Jones crystal from the harmonic ones ω 0 k plotted against the ω 0 k . The different colors again represent longitudinal and transverse branches. The data along with the uncertainty are shown divided by the corresponding temperature and harmonic frequency of the same mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.17 The phonon dispersion relationship for several directions of wavevector k. Here we plot the results of QMC (red triangles) and MD (blue dots) simulations and demonstrate the prediction of equations (6.39) and (6.40) for frequencies of quantum and classical Lennard-Jones crystal. The shown results correspond to the system of N = 108 particles and T * = 0.41. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.18 The phonon dispersion relationship for several directions of wavevector k. The results of QMC (in red ) and MD (in blue) simulations along with the prediction of equations (6.39) and (6.40) are compared against the experiments [START_REF] Endoh | Lattice dynamics of solid neon at 6[END_REF]. . . . . . . 6.19 The comparison of peak widths for two sets of modes corresponding to the {1, 0, 0} direction in k-space. On the left we show the transverse branch and on the right -longitudinal. Results of the Molecular dynamics calculation are plotted in blue whilst the output of the QMC simulation is showed in red. As the indication of insensitivity of the QMC data with respect to the change of the system size, we depict in green the result for the crystal of N = 864 particles at T * = 0.54. In order to avoid overlapping with other data we shifted the point to the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.20 The calculation of heat conductivity from the data of Sec.(6.6). In blue we show the results for the Molecular dynamics calculation whilst in red we plot the quantum Monte-Carlo results. The mode heat capacity is calculated with the same formula given by Eq.(6.42) with the replacement of ω by the mode peak frequency in both cases. . . . . . . xv

List of Tables

Introduction

Since its inception the physics aims to describe the nature around us in all its beauty.

In order to do it in a quantitatively accurate way physicists make use of various mathematical tools and methods. However, in order to apply these tools one must recast the physical system in a mathematically expressible model. Thus, the study of any model or process that resembles the actual one is performed with a certain set of approximations. These approximations are fundamentally connected to our ignorance about the underlying governing laws and arise due to our limitations in calculating the properties of interest for any realistic system. In materials study, i.e. the field of quantum chemistry or molecular dynamics, one usually adopts several such conventions [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF]. One of them is the Born-Oppenheimer approximation [M. [START_REF] Born | Zur Quantentheorie der Molekeln[END_REF], which is commonly applied in particular for the purpose of numerical analysis. This approach assumes that the wave function describing the molecule can be separated into nuclear and electronic parts that can then be studied independently. The assumption is based on the fact that the nuclei are much heavier than the electrons which implies the different time scale for the characteristic processes of these two subsystems [START_REF] Scherrer | On the Mass of Atoms in Molecules: Beyond the Born-Oppenheimer Approximation[END_REF]. Unlike Born-Oppenheimer approximation whose limitations were known for long time [START_REF] Beratan | Failure of the Born-Oppenheimer and Franck-Condon approximations for long distance electron transfer rate calculations[END_REF], another commonly used approximation which typically gets much less attention -that the behaviour of the nuclei can be understood in terms of classical mechanics [START_REF] Kushwaha | Acoustic band structure of periodic elastic composites[END_REF][START_REF] Datta | Photonic band gaps in periodic dielectric structures: The scalar-wave approximation[END_REF]. In the pursuit of more easily performed calculations one discards the machinery to account for the quantum effects, for example zero-point energy or tunneling, hoping that the contribution of said effects to the measurement is sufficiently small. This, however, is not the case for many systems of interest, which can be illustrated on several examples. The thermal wavelength of a proton at 300 K is about 0.1 nm, which is comparable to typical intermolecular distances [R. D. [START_REF] Shannon | Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides[END_REF]. By consequence, proton delocalization certainly affects vibrational properties in many molecular systems. Another consideration is that even at room temperature the zero-point energy (∼ hω/2) for a typical chemical bond with frequency ω is by an order of magnitude larger than the thermal energy attributed to the corresponding degree of freedom (∼ k B T) [START_REF] Brock | 1.25 -Noble-Gas Chemistry[END_REF]. These indications hint that these effects, commonly referred to as nuclear quantum effects (NQE) in the literature, may alter considerably the structure and dynamics of the system observables.

In the rest of this section we shall briefly overview several (semi-)analytical methods that are used for the study of nuclear quantum effects in application to the physics of crystals before giving a brief outline of the methods used here and the organization of the manuscript.

Harmonic crystal

The problem of calculating heat capacity in solids was one of the questions that classical mechanics was unable to address and which in the end ushered quantum mechanics. Let us briefly review the problem.

At temperatures below melting temperature the atoms are arranged in a certain pattern (called crystal lattice) with the position of atoms being at the minimum of the interaction energy. At low enough temperatures the motion of these atoms is fluctuations around their minima and it can be described with a small displacement u.

The resulting model (schematically depicted on the right) can be studied already at the level of classical mechanics: for a system of N particles the Hamiltonian reads as

H = N ∑ i=1 p 2 i 2m + 1 2 N ∑ i,j=1 ∂ 2 V ∂x i ∂x j u i u j
where p i -moment of the i-th particle, u i -its displacement from the minimum, and m -mass of a particle. The potential in the quadratic form can be diagonalised with N modes = 3N normal modes which renders the Hamiltonian to the sum of independent oscillators with frequencies

ω i H = N modes ∑ k=1 p 2 k 2m 3 2 k B T + N modes ∑ k=1 mω 2 k y 2 k 2 3 2 k B T
Upon averaging each mode has a total energy of k B T, thus, the total energy of a crystal is E int = 3Nk B T. This leads to the renown law of DulongPetit [START_REF] Dulong | Recherches sur la mesure des températures et sur les lois de la communication de la chaleur[END_REF] for heat conductivity c class v = 3k B N which is independent of temperature and is uniform for any material. The problem with this simple picture arises from the experimental observations (see Fig .(1.1) for example of solid neon). The classical theory is unable to capture even qualitatively the behaviour at low temperaturesthe limit where it is justified.

The solution to this conundrum was proposed by Debye [START_REF] Debye | Vorträge über die kinetische Theorie der Materie und Elektrizität, Gottinger Wolfskehlvortrage[END_REF] who was able to solve it by using quantized energy excitations. The resulting formula for the heat capacity is

c v = 9Nk B T T D 3 T D T 0 x 4 e x dx [e x -1] 2 (1.1)
. This equation contains the important parameter -the Debye temperature

T D = hv s 2k B 3 6ρ
π which is expressed in terms of Planck's constant h, crystal sound velocity v s and the crystal number density ρ. It controls the number of phonons participating in the calculation. Comparing this result with experimental observations (see Fig .(1.1)) one concludes, that the harmonic model, despite its simplicity, is able to explain behaviour of internal energy (and heat capacity) for many materials at temperatures below the melting one.

Quasiharmonic approach

There several limitations to the famous Debye law. The formula (1.1) was derived in the assumption of harmonic approximation for potential energy. This approach is valid at small temperatures, however, with the increasing temperatures the anharmonicities of the potential become more important thus deviating from the Debye law. The prominent example are the semiconductors, and we plot the specific heat at various temperatures for the case of silicon (Fig.(1.2)). This figure shows that the deviations are significant for a relatively low temperatures of around 200K which is smaller than the Debye temperature of 600K and considerably below melting temperature ( 1400K).

In order to make use of the harmonic model one typically treats the anharmonic interactions as a small perturbation. Following the ideas of Einstein [START_REF] Einstein | Die Plancksche Theorie der Strahlung und die Theorie der spezifischen Wärme[END_REF] the vibrations inside a crystal are described by the normal modes which as we have seen above in the case of pure quadratic interactions is equivalent to the analysis of ensemble of harmonic oscillators. Energy of each normal mode is quantized and the excitations from the ground energy is called phonons [START_REF] Tamm | Über die Quantentheorie der molekularen Lichtzerstreuung in festen Körpern[END_REF]. This formalism is not quite suited for the proper description because it involves only second order terms in the expansion of the particles interaction energy. The popular recipe is to use the quasiharmonic approach, which adopts the harmonic oscillator viewpoint and modifies the frequencies in order to include thermal expansion corrections [START_REF] Leibfried | Theory of anharmonic effects in crystals[END_REF][START_REF] Brüesch | Phonons: Theory and experiments I: Lattice dynamics and Models of interatomic forces[END_REF]. It is usually assumed that the anharmonicities of the potential affect the observables of the system at high temperatures so that the quasiharmonic approach is valid at low temperatures when the phonon occupation number is low, which is strictly speaking not the case, since the zero-point energy allows particles to fluctuate from the equilibrium positions thus 200 400 600 800 1000 1200 1400 increasing high-frequency occupation numbers. Following the quasiharmonic approach, anharmonic effects are considered weak so that atomic forces and phonon frequencies are renormalized by accounting only for thermal expansion. It is motivated by the fact that the anharmonic effects make an important contribution to the thermal expansion of the crystal which can be analyzed in experiment. One typically describes the frequency shift by a mode Grüneisen parameter, γ k = -(V∂ω k )/(ω k ∂V), where V is the volume and ω k being the frequency of the kth mode. The parameter appears naturally when writing the equation of state in the quasiharmonic approximation.

T(K) 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 c V / k B N
Quasiharmonic approach was reasonably successful in describing certain properties of the crystals. In particular it was able to provide with reasonable accuracy the thermal expansion of silicon at varying temperatures [START_REF] Xu | Theory of the thermal expansion of Si and diamond[END_REF][START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF]] so it became the common tool for the analysis [START_REF] Toher | High-throughput computational screening of thermal conductivity, Debye temperature, and Grüneisen parameter using a quasiharmonic Debye model[END_REF]. However, observations of frequency shift in various semiconductors already from early papers demonstrated that this shift cannot be fully accounted by the thermal expansion alone [Nilsson and Nelin, 1972] and various numerical studies supported the role of anharmonicities in the temperature shift of phonon frequencies [START_REF] Lang | Anharmonic line shift and linewidth of the Raman mode in covalent semiconductors[END_REF][START_REF] Narasimhan | Anharmonic self-energies of phonons in silicon[END_REF].

Computational approaches

Computational approaches

Due to the complexity of the problem for the exact analytical solution one is pushed to the numerical computations.

A promising approach that allows treatment of the NQE is the use of the quantum thermal bath [START_REF] Dammak | Quantum Thermal Bath for Molecular Dynamics Simulation[END_REF][START_REF] Ceriotti | Nuclear quantum effects in solids using a colored-noise thermostat[END_REF]. In this approach the delocalization of a particle is represented by the stochastic motion caused by the random external force. In this regard the method is similar to the classical Langevin thermostat with the force constructed as a specific power spectrum obtained from the quantum dissipation-fluctuation theorem [START_REF] Callen | Irreversibility and Generalized Noise[END_REF]. The method provides reasonable results for certain range of systems and can be used to refine approximate models for other numerical techniques, yet as many semiclassical approaches is prone to the zero-point energy leakage -the issue when part of the energy is unphysically transferred from the high frequency modes to the low frequency modes [START_REF] Brieuc | Zero-Point Energy Leakage in Quantum Thermal Bath Molecular Dynamics Simulations[END_REF].

In our work we shall follow the path-integral Monte-Carlo (PIMC) method for the calculation of imaginary time correlation functions. The approach necessitates the knowledge of the interaction energy between atoms. In insulator systems the interaction between atoms is governed by the behavior of particles electron clouds -at short atom separation non-bonded overlap between electron clouds lead to the strong repulsion whilst at long distances the correlation between these clouds result in attraction of atoms. Thus, the use of the exact potential in the PIMC simulation would require solution of the system of Shrödinger equations at each time of the crystal evolution. Instead, typically one adopts another approach and uses the empirically obtained potential for description of the particles interaction [START_REF] Allen | Computer simulation of liquids[END_REF] and in this thesis we restrict our attention to such potentials.

Thesis overview

In the current thesis we are interested in computing transport coefficients, taking into account NQE, without resorting to the approximate schemes. As a main tool for the analysis we shall use path-integral Monte-Carlo approach which allows to obtain time correlation functions for the quantities of interest in an at least in principle exact way. In practice, however, the calculations are affected by statistical noise, and the accuracy of the computation is determined by the simulation length which results in the finite precision. In general, path-integral Monte-Carlo calculations of real-time evolution is numerically infeasible due to a large statistical error. However, as we shall see in the following chapters, PIMC can access directly correlations in imaginary time. The problem of getting real time response function is then migrated to analytically continue the imaginary correlation functions to real times. Alternatively, one can study the spectral density that describes the correlations. In the case of imaginary time functions, however, one is faced with the ill-defined Laplace transform due to the finite errorbars of the PIMC calculation. We will employ a stochastic reconstruction approach for this inversion problem and introduce validation prescription in order to uniquely fix the parameters of the reconstruction. To our knowledge, combining the Bayesian stochastic reconstruction with the validation procedure has never been reported before. We will show that the combination actually provides very encouraging results.

The manuscript is organized as follows. In the first two chapters we introduce general formalism and machinery for the calculation of imaginary time correlation functions and spectral inversion. In chapter 4 and 5 we then apply these techniques to the study of simple harmonic oscillator and several oscillator-like models in order to test their validity. Then, in chapter 6 we analyze the solid neon with the interparticle interaction described by pairwise Lennard-Jones potential.

Methods

Heat transport properties are usually obtained by measuring heat transfer caused by a small temperature gradient. In this chapter we lay basis for the linear response theory which allows to compute changes of thermodynamic properties when a small external perturbation is applied. The Green-Kubo formalism relates the linear response of the system to correlation functions obtained at thermodynamic equilibrium. Therefore, we start with the calculation of the time correlation function for several operators. In the Green-Kubo approach [START_REF] Green | Markoff Random Processes and the Statistical Mechanics of TimeDependent Phenomena. II. Irreversible Processes in Fluids[END_REF] the correlation functions are composed of nothing else but the heat current operators. In the framework of single-mode relaxation time approach [START_REF] Srivastava | The Physics of Phonons[END_REF] one can analyze the properties of each phonon from the mode correlation functions. Whilst there exist various approximate methods to perform such calculations, no analytical solution can be expected for a general interacting system. As a consequence, one is forced to resort to numerical approaches in order to obtain the result and to this end we will use Monte-Carlo techniques [START_REF] Metropolis | The Monte Carlo Method[END_REF]Ulam, 1949, Bernu and[START_REF] Bernu | Quantum simulations of complex many-body systems: From theory to algorithms 10[END_REF].

In order to obtain the heat conductivity from Green-Kubo we need to have access to real-time correlation functions, the problem too hard, in general, even for numerical approaches. However, path-integral Monte-Carlo methods allow calculation of imaginary times correlation function which are the analytical continuation of their real-time analogs. Obtaining the spectral function of the corresponding correlations necessitates solving the inverse problem for the underlying linear transformation. However, the inversion of the corresponding kernel is an ill defined problem [START_REF] O'sullivan | A statistical perspective on ill-posed inverse problems[END_REF], so that the output of the calculation is extremely sensitive to the accuracy of the input data. In the present thesis, we employ a stochastic reconstruction approach [START_REF] Sandvik | Stochastic method for analytic continuation of quantum Monte Carlo data[END_REF][START_REF] Habershon | Quantum mechanical correlation functions, maximum entropy analytic continuation, and ring polymer molecular dynamics[END_REF] in order to produce the desired spectrum. In order to improve the reliability of the result, we will introduce a recipe to increase the accuracy of the time correlation function calculation. Namely, we will introduce improved numerical estimators for correlation functions that involve momentum operators that are in general very noisy. Inspired by the so-called virial estimator [START_REF] Barker | A quantum-statistical Monte Carlo method; path integrals with boundary conditions[END_REF][START_REF] Parrinello | Study of an F center in molten KCl[END_REF] for kinetic energy, we describe how a low-variance estimator can be obtained for any imaginary time correlation function that involves operator of momentum. We split this part of the thesis into two segments: in the first one we set the ground for the calculation by describing the Green-Kubo approach and the related inverse problem. In the next chapter we lay the basics of the path integral method for the calculation of time correlation function and overview the numerical tools that we will use to evaluate said correlations. We will also show the prescription on reducing the errors of the computations.

Linear response theory

The first building block of our approach is the theory of linear response introduced by Green [START_REF] Green | Markoff Random Processes and the Statistical Mechanics of TimeDependent Phenomena. II. Irreversible Processes in Fluids[END_REF] and Kubo [START_REF] Kubo | Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems[END_REF]. This formalism is used to calculate Chapter 2. Methods the response of an observable to a small perturbation of the system's Hamiltonian of a system in the linear approximation. In this chapter we give a brief review of general linear response theory [START_REF] Baroni | Heat transport in insulators from ab initio Green-Kubo theory[END_REF].

Heat transport in condensed matter system is an example of a macroscopic process. For this type of processes one can typically attribute an extensive variable A which can be expressed, from it's definition, as a sum of its values for each subsystem. In the limit of a macroscopic system one can introduce a density a such that:

A V = V a(r, t)dr, (2.1)
where V denotes the volume of the system and bold indexes denote vectors, that is r ≡ r = {r x , r y , r z }. For a locally conserved extensive quantity one can write an equation of continuity for the density a:

∂a(r, t) ∂t + div j A (r, t) = 0, (2.2)
in the absence of a total current through the system boundary. In the last equation j A is the density current through the system surface. Performing a Fourier transform over space coordinates in the form ã(k, t) = a(r, t)e ikr dr, one writes:

∂ ã(k, t) ∂t + ik jA (k, t) = 0 (2.3)
This equation shows that the quantities with longer wavelength k change slower in comparison to the short wavelength components and one may expect the adiabatic decoupling of these Fourier densities, which are typically called hydrodynamic variables. Without loss of generality one can shift the minimum (equilibrium) value of conserved quantities A to zero. Then, at equilibrium all the corresponding densities and currents vanish. Treating a system away from equilibrium we shall assume that the the wavelength along with the time scale of characteristic processes are sufficiently long for the local thermal equilibrium to settle in. A local equilibrium is characterized by certain value of the local temperature, pressure, chemical potential etc satisfying local equation of state. For small enough deviation from the equilibrium one can express the time derivative of ã(k, t) as the linear combination of the density itself:

∂ ãi (k, t) ∂t = ∑ j B ij ãj (k, t) (2.4)
with i and j running over all the conserved densities. Combining the equation (2.4) with the previous one (2.3) we obtain for the current:

ji (k, t) = i k k 2 ∑ j B ij ãj (k, t) (2.5)
One expects that in isotropic media coefficients B ij are spherically symmetric functions of k and the current vanishes at k = 0, since the a non-vanishing value would indicate the dependence of current at long distances on the density fluctuations which contradicts the idea of local equilibrium. Thus, one can suggest that in the limit of k → 0 the coefficients behave as B ij ∝ k 2 b ij and one can study a macroscopic (k = 0) components of stationary current, J i = 1

V j i (r)dr. These currents are 2.1. Linear response theory 9 related to the corresponding density gradients W i = 1 V ∇a i (r)dr as:

J i = ∑ j b ij W j (2.6)
This expression is reminiscent of celebrated Onsager's relation [START_REF] Onsager | Reciprocal Relations in Irreversible Processes. I[END_REF]:

J i = ∑ j L ij F j (2.7)
where F j denotes thermodynamic forces, or affinities, whilst linear coefficients L ik satisfy the celebrated Onsager's reciprocity relations L ik = L ki . Indeed, equation (2.6) can be brought to the form of (2.7) by defining susceptibilities κ such that density gradients can be related to thermodynamic forces as:

W j = ∑ k κ jk F k (2.8)
thus leading to the Onsager's form with

L ij = ∑ k b ik κ kj .
Let us consider a system that consists of N classical particles interacting via potential U 0 . In 3D the system is described with 3N degrees of freedom. We write the initial (unperturbed) Hamiltonian for this system as

H 0 = N ∑ i p 2 i 2m i + U 0 (r 1 , r 2 , r 3 . . . ).
(2.9)

One does not need to restrain the form of the potential and thus it can include all various combinations of particles. Suppose now that the potential interaction between atoms is modified with a time-dependent part:

U(R, t) = U 0 (R) + U (R, t) (2.10)
We assume that the perturbation U can be expressed in terms of densities a i (r, t):

U (R, t) = ∑ i V u i (r, t)a i (r, t)dr (2.11)
Then, according to [START_REF] Kubo | Statistical-Mechanical Theory of Irreversible Processes. II. Response to Thermal Disturbance[END_REF] one can write the linear response to the perturbation for the i-component of the current:

j γ i (r, t) = 1 k B T ∑ k t -∞ dt V dr j γ i (r, t) ȧk (r , t ) 0 u k (r , t ) (2.12)
Here γ denotes the component of a vector. We used . . . 0 to show the averaging over energy density states of the unperturbed Hamiltonian. Using the continuity equation (2.2) one can rewrite this expression

j γ i (r, t) = - 1 k B T ∑ k,σ t -∞ dt V dr j γ i (r, t) ∂ ∂r σ j σ k (r , t ) 0 u k (r , t ) = = 1 k B T ∑ kσ t -∞ dt V dr j γ i (r, t)j σ k (r , t ) 0 ∂ ∂r σ u k (r , t ) (2.13)
where we integrated by parts and neglected any surface contributions, assuming Chapter 2. Methods that the fields in the expression decay sufficiently fast. In order to make the comparison with Onsager's relation (2.7) we integrate both sides of the equation (2.12) over space and take into an account the relation L ik = L ki . The resulting expressions reads as

J γ i = 1 V drj γ i (r) (2.14) L γσ ik = V k B ∞ 0 dt J γ i (t)J σ k (0) 0 (2.15) F σ k = 1 VT dr ∂ ∂r σ u k (r)
(2.16)

Time correlation function of two operators

From the previous section it is clear that quantities that will be most interesting to us in the following are equilibrium time correlation functions, which will play an essential role in application of the linear response theory. In general correlation function for the system described by the Hamiltonian Ĥ that involves observables A at time t 1 and B at time t 2 can be written as the thermodynamic average of a product of operators Â(t 1 ) = e it 1 Ĥ/h Â(0)e -it 1 Ĥ/h and operator B(t 2 ) = e it 2 Ĥ/h B(0)e -it 2 Ĥ/h :

C AB (t 1 -t 2 ) = Â(t 1 ) B(t 2 ) = 1 Z β tr[ Â(t 1 ) B(t 2 )e -β Ĥ ] (2.17)
where

Z β = tr[e -β Ĥ ].
From its definition we observe several properties of the correlation function:

-Invariance to time translation Correlations between operators does not depend on the absolute value of time, rather it depends on the time difference

Â(t + t ) B(t ) = 1 Z tr[e -β Ĥ Â(t + t ) B(t )] = = 1 Z tr[e -β Ĥ e i Ĥ(t+t )/h Â(0)e -i Ĥt/h B(0)e -i Ĥt /h ] = 1 Z tr[e -β Ĥ Â(t) B(0)] = Â(t) B(0) -Inversion of time The correlation function is symmetric around time t = ihβ/2 C AB (-t) = Â(-t) B(0) = 1 Z tr[e -β Ĥ e -i Ĥt/h Â(0)e i Ĥt/h B(0)] = = 1 Z tr[e i Ĥt/h B(0)e -β Ĥ e -i Ĥt/h Â(0)] = 1 Z tr[e -β Ĥ e β Ĥ e i Ĥt/h B(0)e -i Ĥ(t/h-iβ) Â(0)] = = B(t -ihβ) Â(0) = C BA (t -ihβ),
where in deriving both equations we used the cyclicity of a trace operation tr(ABC) = tr(CAB).

Spectral function

We now introduce another important quantity that we will use in the following.

Let us define a spectral density function S AB (ω) that corresponds to the Fourier transform of the correlation function C AB between operators  and B:

S AB (ω) = 1 2π ∞ -∞ C AB (t)e -iωt dt (2.18)
We notice from the definition and using the properties of the correlation function from above that

S AB (-ω) = 1 2π ∞ -∞ e iωt C AB (t)dt = e -βω S BA (ω) (2.19)

Imaginary time correlation functions

For the purpose of numerical analysis with PIMC one deals with the correlations defined for imaginary times, t = ihτ, where we also incorporated inside the definition the factor of h. The operators in the Heisenberg representation will then take a form Â(τ 1 ) = e -τ 1 Ĥ Â(0)e τ 1 Ĥ , with the correlation function defined in imaginary times as

C AB (τ) = 1 Z tr[A(τ)B(0)e -β Ĥ ]
as the analytical continuation t = ihτ. For such correlation functions one can modify the definition for C(τ) :

C AB (τ) = ∞ 0 dω S AB (ω)e -hωτ + S BA (ω)e -hω(β-τ) ,
(2.20)

where we used property (2.19). Using this spectral function that one can reexpress the coefficients L ik from equation (2.15) as spectral density of the current-current correlation function at ω → 0:

L γσ = V 2k B S γσ (ω = 0) (2.21)

Inverse Problem

As we mentioned above, the PIMC provide us with the correlations defined for imaginary times. Omitting for the moment the question of actual calculation of the correlation functions (we shall return to this issue in the next chapter) we now consider the following problem related to the spectral functions. When trying to obtain the spectral function from imaginary time correlations as indicated in the previous section from inverting Eq.(2.28), one is quickly faced with a so-called illdefined problem [START_REF] O'sullivan | A statistical perspective on ill-posed inverse problems[END_REF] meaning that the reconstructed spectrum is very sensitive to the accuracy of the input data. There exist many different approaches to tackle this problem [START_REF] Alifanov | Inverse heat transfer problems[END_REF]. We start this section by briefly reviewing the classical method which make use of a kernel parameter regularization [START_REF] Beck | Inverse heat conduction: Ill-posed problems[END_REF][START_REF] Kaltenbacher | Iterative regularization methods for nonlinear ill-posed problems[END_REF]. However, regularized solutions are somewhat unsatisfactory for our purposes. First, it is often difficult to impose physical constraints on the solution, e.g. a non-negative spectral function. Second, the solutions can depend very sensitively on the regularization parameter, such that their uncertainty is very difficult to control. Therefore, in the following however we will adopt a different route, namely we will use a so-called stochastic, or maximum entropy, approach. As we shall see, its formulation suggests that it is a parameter-free method of spectral reconstruction which makes it very attractive. However, as we will point out, this is, for most of the cases, an illusion, and the method nevertheless depends implicitly on several external parameters, some of which are not immediately apparent. In order to uniquely fix the parameters, both explicit as well as implicit ones, we propose the procedure for verification of the reconstruction based on statistically independent data. As we will show, this step is essential to obtain reliable and accurate results.

General inverse problem

Let us start with a general definition of the problem. The equation of our interest defined in the previous section -see equation (2.20)) -is an example of the inverse problem [START_REF] Kabanikhin | Theory and Applications[END_REF]. In general this problem can be postulated as a linear map between two Hilbert spaces X and Y y = Kx (2.22)

where x and y belong to X and Y respectively, whilst K is a some compact linear operator, typically referred to as a kernel. The name of the problem comes from the fact that for such problems one knows the set of measurements y along with a given operator K and the goal is to reconstruct set x. Moreover, in our case samples y that we obtained are not exact data points but they are calculated within the uncertainty , which we we will discuss later. Thus, the problem we aim to solve is

y = Kx + (2.23)
Before proceeding further, we recall the conditions for a stable inversion:

-Existence Existence of the solution is guaranteed if the map K is surjective [Bourbaki, 1970]:

∀y ∈ Y ∃x ∈ X : y = Kx (2.24)
that is, each element from the space of values x is mapped to at least one element of Y

-Uniqueness Solution is unique if the map K is injective [Bourbaki, 1970]: ∀x, x ∈ X : y = y =⇒ x = x (2.25)
i.e. for each element x ∈ X there exists only one element y ∈ Y such that y = Kx. The condition for this is that ker(K) = 0 where ker(K) = {x ∈ X|Kx = 0} With this assumption we are able to define an inverse map: K -1 which maps the elements of Y back to the space X. If accuracy of the data does allows for more than one solution a usual prescription is to impose additional constraints.

-Stability In order to guarantee stability of the solution with respect to errors we need to ensure continuity of inverse map K -1 . Violation of stability condition is the most difficult task to address

An inverse problem that satisfy all these criteria is called well-posed [START_REF] Hadamard | Sur les problèmes aux dérivés partielles et leur signification physique[END_REF]. If some of them are not met then the problem is ill-posed. In the example that we are interested, namely in solving the equation (2.20), in the data that we are using for correlation function are determined with certain accuracy, and the illdefinition of the inversion means that there exist a huge number of spectral functions that would produce the result within this accuracy.

Calculation of a unique solution can be done in a rather straightforward way by regularization of the kernel [START_REF] Tikhonov | Numerical methods for the solution of illposed problems[END_REF]. From the mathematical point of view, a regularization for the inverse problem y = Kx is defined as a family of continuous mappings R λ : Y → X that converges to the inverse of K in the limit λ → 0. That is, the original problem is replaced with well-posed problem y = R -1 λ x. The regularization parameter λ determines how well R λ approximates the inverse operator. Let us demonstrate how one can implement regularization in practice on a simple example.

Regularization of the kernel

In order to preform analysis of the equation (2.20) we use numerical calculations. Thus, we are obliged to introduce some parametrization of the spectral density S(ω).

For our example we consider a spectral function to be described by the collection of δ-peaks:

Ŝ(ω) = P ∑ p=1 A p δ(ω -ω p ) (2.26)
where P is some arbitrary number of fitting parameters A p . The frequencies ω p are assumed to be fixed, for example, they form a linear progression on a range [ω min , ω max ]. In our problem vector y consists of correlation functions at M discrete imaginary time points. Thus, the kernel K from the general equation y = Kx can be considered as a map between P-dimensional space of ω and M-dimensional space of y. In a mathematically rigorous sense one cannot demonstrate in general the existence of the solution for this problem. However, one can prove [START_REF] Horn | Matrix Analysis[END_REF] that the problem admits a solution in a "least squares" sense -one can minimize the difference |y -Kx| on a subset of frequency space. One can devise different approaches for this minimization [START_REF] Golub | Matrix Computations[END_REF], here we show one such example. We use the form for spectral function (2.26) to solve eq.(2.20) by minimizing the functional χ 2

χ 2 = M-1 ∑ α=0 [C(τ α ) -C(τ α )] 2 σ 2 (τ α )
Here we used

C(τ α ) = P ∑ p=1 A p e -ω p τ α + e -(β-τ α )ω p
Minimization of χ 2 means that we need a first derivative of χ 2 to be equal to 0:

0 = dχ 2 dA p = M-1 ∑ α=0 1 σ 2 (τ α ) C(τ α ) - P ∑ p=1 A p e -ω p τ α + e -(β-τ α )ω p e -ω p τ α + e -(β-τ α )ω p
(2.27) One can satisfy this set of equations by demanding that each in [. . . ] should vanish:

C(τ α ) = P ∑ p=1 A p e -ω p τ α + e -(β-τ α )ω p (2.28)
Chapter 2. Methods

We can write this expression in a matrix form as C = KA, with

C = (C(τ 1 ), C(τ 2 ), . . . , C(τ M )) † ≡ (C 1 , C 2 . . . C M ) † , (2.29) A = (A(ω 1 ), A(ω 2 ), . . . A(ω P )) † ≡ (A 1 , A 2 , . . . , A P ) † .
(2.30)

and

K =     K 1 (τ 1 ) K 2 (τ 1 ) . . . K P (τ 1 ) K 1 (τ 1 ) K 2 (τ 1 ) . . . K P (τ 1 ) . . . . . . . . . . . . K 1 (τ M-1 ) K 2 (τ M-1 ) . . . K P (τ M-1 )     where K i (τ α ) = e -ω i τ α + e -(β-τ α )ω i
and we have i, j

= (K T K) -1 K T C (2.31)
or if we write all the indexes:

A f = (K T K) -1 f j K T jα C α (2.32)
In the last line we assumed the summation over the repeating indexes. In general the error can be estimated from the calculation of variance, which in our case reads as

Var = A 2 -A 2 (2.33)
Due to the nonvanishing errors on values of C the first term entering equation (2.33) could be written as follows:

A m (A T ) m = (K T K) -1 mj K T jα C α C T β K βn (K T K) -1T nm (2.34) Since C α C T β = σ 2 δ αβ and (K T K) -1T = (K T K) -1 we can simplify this expression into A m (A T ) m = (K T K) -1 mj K T jα δ αβ σ 2 α K βn (K T K) -1 nm (2.35)
For constant σ 2 we can simplify this formula to:

A m (A T ) m = (K T K) -1 mj K T jα δ αβ K βn (K T K) -1 nm σ 2 = = (K T K) -1 mj (K T K) jn (K T K) -1 nm σ 2 = δ mn (K T K) -1 nm σ 2 = (K T K) -1 mm σ 2 (2.36)
The reason for the instability of the inversion is the linear dependence between elements of the kernel K which renders the determinant det(K T K) in the matrix inversion equal to 0. A usual regularization scheme that produces a well-defined inversion of K T K involves some non-zero parameter. For example, one may consider adding a constant term λ to the diagonal elements of K T K, thus making the determinant positive defined. Then the expression for A(λ) can be written as a function of λ:

A f (λ) = (K T K + λI) -1 f j K T jα C α (2.37)
The solution for A is then formally obtained in the limit λ → 0: A = lim λ-→0 A(λ).

In practice however this limit often cannot be reached as the results can get very sensitive to the numerical precision involved in the operations, and one is forced to use a finite parameter which one chooses sufficiently small such that the solution is insensitive to the changes of λ. On Fig.

(2.1) we show example of fitting for different regularization values λ. Here we perform fit for a correlation function data obtained in the Sec.(4.1). The fit is done for a fixed values of ω p with 3 adjustable (but not constrained) parameters A p . We see that for a large range of parameters λ the sum rule for the coefficients ∑ p A p = 1 is quite closely fulfilled automatically, however, for sufficiently small λ's some of the coefficients A p become negative which is clearly not acceptable for a positively defined functional S(ω).

10 9 10 7 10 5 10 3 0.0

0.2 0.4 0.6 0.8 A A 1 A 2 A 3 FIGURE 2
.1: Example of the regularization procedure for a kernel of the inverse problem. Here we perform a 3 parameter fit and plot values of χ 2 for the fixed set of ω p and corresponding A p as defined by eq.(2.37) as a function of regularization parameter λ. Set of correlation function data points and corresponding errors are taken from the first test case from Sec.(4.1)

Stochastic method for spectral reconstruction

Over the years the problem of obtaining spectral information from imaginary time correlation functions was tackled from various perspectives. One of the most recent emerging topic of research is the use of neural network for this task [START_REF] Kades | Spectral reconstruction with deep neural networks[END_REF]. An obvious drawback of this approach is that it requires construction of a training data set -data for a model with a known output. The model should be taken such that it captures the important features of the wanted spectral function [START_REF] Fournier | Artificial neural network approach to the analytic continuation problem[END_REF], yet it should be sufficiently general in order to not introduce artifacts to the measurements.

For our calculation we adopt the approach that tackles the inversion problem from a statistical point of view [START_REF] Sandvik | Stochastic method for analytic continuation of quantum Monte Carlo data[END_REF][START_REF] Habershon | Quantum mechanical correlation functions, maximum entropy analytic continuation, and ring polymer molecular dynamics[END_REF][START_REF] Levy | Implementation of the maximum entropy method for analytic continuation[END_REF]. In this method the spectral function S ME is obtained as an average over possible spectra S(ω) that describe the correlation function data within their accuracy. In practice one uses some parametrization in order to Chapter 2. Methods define spectral configurations. These spectra are summed with a weight determined by the probability that they are the exact model given the data set (C, σ 2 ), S(ω) ME = DS p(S|C, σ 2 )S(ω).

(2.38)

Here DS the phase space element associated with the parametrization of S(ω), whilst the vector C = (C(τ 1 ), C(τ 2 ), . . . , C(τ M )) † ≡ (C 1 , C 2 , . . . C M ) † denotes the totality of the correlation function data points, and σ 2 describes the statistical uncertainty of these data in the form of a covariance matrix. Making use of the Bayes formula, one can write

p(S|C, σ 2 ) = p(C, σ 2 |S) p(C, σ 2 ) p(S), (2.39)
where p(Cσ 2 |S) is the likelihood function to obtain C given S, p(C, σ 2 )) = dSp(C|S) a normalization factor, and p(S) is a prior probability density that incorporates any information on the spectrum that has before performing the reconstruction. Assuming that the error of the data C is given by the central limit theorem the likelihood function obeys Gaussian statistics so that we can write

p(C|S, σ) ∝ e -1 2 (C-C[S])(σ 2 ) -1 (C-C[S]) = e -1 2 χ 2 [S]
(2.40)

where we introduced χ 2 (S) = (C -C[S])(σ 2 ) -1 (C -C[S]
) and σ 2 as the covariance of the data. Here vector C[S] is computed from the definition of correlation function

C(τ) = ∞ 0 dωS(ω) e -hωτ + e -hω(β-τ) (2.41)
and it is obtained by inserting a trial spectrum S into the r.h.s. of equation (2.20) and computing the resulting M correlation values. In the case of a trial spectrum consisting of the amplitudes A(ω p ) which are defined for a set of N ω discrete frequencies on a regular grid, namely

Ŝ(ω) = N ω ∑ p=1 A p δ(ω -ω p )
trial correlation function C[S](τ α ) can be written using equation (2.41) in the following form

C[S](τ α ) = N ω ∑ p=1 A(ω p ) e -hω p τ α + e -h(β-τ α )ω p .
(2.42)

In traditional maximum entropy methods, equation (2.38) is solved with saddle point approximation by minimizing the functional

F = 1 2 χ 2 [S] -H[S].
Here, H[S] is an entropic functional, which punishes irregular solutions which presumably would lead to an overfitting of the resulting spectrum -feature of the reconstruction which appears when statistical errors contained in the data C(τ α ) influence to a significant extent the result and is well-known in the machine learning [START_REF] Cawley | Preventing Over-Fitting during Model Selection via Bayesian Regularisation of the Hyper-Parameters[END_REF]. A usual prescription for H[S] is to use the associated Shannon entropy [C. E. [START_REF] Shannon | A mathematical theory of communication[END_REF]:

H[S] = γ N ω ∑ p=1 A p ln A p ,
with a coefficient γ controlling the regularisation of the solution [START_REF] Linden | Maximum-entropy data analysis[END_REF].

In this work, however, we follow a slightly different route and employ the stochastic analytical inference, or stochastic maximum entropy [START_REF] Fuchs | Analytic continuation of quantum Monte Carlo data by stochastic analytical inference[END_REF] method. In this approach expression(2.38) is evaluated by Monte-Carlo sampling over space of possible spectral configurations, which are constrained to positive values of S since the spectral density is strictly non-negative. The term χ 2 [S] can hence be considered as an effective energy functional. In order to make this analogy more apparent, one can introduce an additional parameter in the form of an effective inverse temperature Θ as,

S(ω, Θ) ME = Z(Θ) -1 DS S(ω)e -1 2 Θχ 2 [S] .
(2.43)

Here the normalization factor Z(Θ) = dS exp[-Θχ 2 /2] is an effective partition function. Such approach allows to integrate traditional maximum entropy method as a mean field version of Eq.(2.43), where one uses as an estimate of the spectrum the minimum of the mean field free energy

F MF (θ) = 1 2 χ 2 [S] -Θ -1 H[S].
In the following we make the simplifying assumption that data points are uncorrelated (which is strictly speaking is not the case, but the end conclusion is not affected), so that the covariance matrix is diagonal. As a result, we can write the energy functional χ 2 [S] in the form

χ 2 = M-1 ∑ α=0 [C(τ α ) -C[S](τ α )] 2 σ 2 (τ α ) , (2.44)
with σ 2 (τ α ) the statistical uncertainty of the data point α. This approach, however, is not a parameter free calculation -in fact, various simulation parameters, such as grid density in spectrum space, effective inverse temperature Θ etc play a role analogous to the explicit regularization parameters [START_REF] Ghanem | Average spectrum method for analytic continuation: Efficient blocked-mode sampling and dependence on the discretization grid[END_REF].

Validation procedure

Despite the vast body of literature on the subject, there does not exist a prescription that allows one to fix unambiguously the value of Θ. Some arguments [START_REF] Fuchs | Analytic continuation of quantum Monte Carlo data by stochastic analytical inference[END_REF] have been proposed for choosing Θ = 1. In contrast, other arguments were put forward that suggested picking for Θ the value Θ * that maximizes Z(Θ), which is argued to also maximize the posterior probability P(θ|C) [START_REF] Calvetti | Inverse problems: From regularization to Bayesian inference[END_REF]. The latter proposal, which corresponds to a balance between energy and entropy dominated solutions, requires however a full free energy calculation. One can also note that increasing the value of Θ is effectively equivalent to rescaling the uncertainties on the data points, a procedure that may lead to overfitting. The corresponding effect on the validation procedure is discussed in following parts where we consider test models. Contrary to this proposals, we do not try to justify an a priori choice of effective temperature. Instead, we propose an unambiguous and efficient method to determine the optimal value of Θ. We start with an initial data set, C(τ α ), which is generated with known statistical uncertainty σ 2 (τ α ) by means of, for instance, a path integral simulation of the considered model. For another source of data we can consider cases where C(τ) is known analytically, so that "synthetic" data could be generated from the exact solution, and a controlled uncertainty can be introduced. Starting from any of these initial data, we reconstruct the spectrum S ME (Θ), according to the scheme introduced above. The spectrum is described by P degrees of freedom A(ω p ), and we obtain it through a Monte-Carlo sampling of Eq.(2.43) Chapter 2. Methods for a given value of Θ. We denote CΘ (τ α ) the correlation function associated with this average spectrum. Incidentally, this argument can be extended to include any other possible parameter, such as number of discretization points N ω , different discretizations on an ω-grid, etc. In order to determine the optimal choice of parameters and discriminate between different models for S(ω), we combine the maximum entropy approach with a validation procedure borrowed from the statistical learning theory [START_REF] Mehta | A high-bias, low-variance introduction to machine learning for physicists[END_REF]. One generates P new sets of validation data, C val,i (τ α ) (i = 1, . . . , P ), by using the same technique (even not necessarily with the same accuracy, analogously to neural network studies, see for example [START_REF] Barry-Straume | An evaluation of training size impact on validation accuracy for optimized convolutional neural networks[END_REF]) that we used to produce the initial data set, and determine the associated χ 2 val as

χ 2 val = 1 P P ∑ i=1 M-1 ∑ α=0 [ CΘ (τ α ) -C val,i (τ α )] 2 .
(2.45)

This quantity can be interpreted as a measure of the difference between the estimate CΘ (τ α ) and the exact correlation function, denoted by C exact (τ α ). We can see explicitly by writing

χ 2 val = 1 P P ∑ i=1 M-1 ∑ α=0 [ CΘ (τ α ) -C exact (τ α ) + C exact (τ α ) -C val,i (τ α )] 2 (2.46)
In the limit of large P and assuming that the average over the validation data returns the exact correlation function, we obtain

χ 2 val = M-1 ∑ α=0 [ CΘ (τ α ) -C exact (τ α )] 2 + M-1 ∑ α=0 σ 2 val (τ α ). (2.47)
The first term is the distance of the estimate to the exact data, while the second is the variance of the validation data, which is independent of Θ. The choice of parameters will therefore be eventually influenced only by the behavior of the first term.

Calculation of imaginary time correlation functions

Path Integral formalism

In this chapter we return to the problem of calculating imaginary time correlation functions, the essential block for the analysis of the dynamical properties. Pathintegral Monte-Carlo (PIMC) techniques, a method of choice for the computation of the general observables at thermal equilibrium [START_REF] Feynman | Statistical Mechanics: A Set Of Lectures[END_REF], allows to obtain them in exact, that is systematically unbiased, way within statistical error. However, a straightforward application of PIMC to the calculation of real time operator evolution is unpractical, since the statistical error is not controlled. Nevertheless, PIMC yields exact correlation functions in imaginary time which contains, at least in principle, exactly the same analytical information, as we have seen in the previous chapter. Linear response in imaginary time by means of the analytical continuation provides the real time response, which constitutes a recipe for calculating dynamical observables. But, since the reconstruction of the spectral function is an ill-posed problem, controlling (and reducing) statistical error of PIMC computations is essential, and we address this problem by constructing estimators with reduced variance. We start this chapter by introducing path-integral calculation procedure for several properties of the system, most notably for a general imaginary time correlation function for two operators defined at distinct imaginary times. We shall see that the formalism of path integral lends itself to the numerical calculation and we then follow this chapter by presenting Monte-Carlo techniques that will implement the ideas of the path integral formalism and allow us to finally obtain the correlations of interest. Next, we pass to a more technical problem of representing the operator in the form that is most suitable for the calculation and we compare several such representation for a couple of quantities of interest and discuss the performance of the computation.

General approach

Let us begin by obtaining an expression for a partition function which is arguably the most important quantity that describes the system [START_REF] Landau | Course of Theoretical Physics[END_REF]. For simplicity we consider a model which contains only one single particle placed in the external potential and which is described by one spatial degree of freedom, which we denote r. The approach can be easily generalized to systems consisting of many particles and characterized by numerous degrees of freedom. The Hamiltonian of the model in question is written as

Ĥ = K + Û(r) = p2 2m + U(r) (3.1)
Chapter 3. Calculation of imaginary time correlation functions

Here U(r) characterizes external potential, K is a kinetic energy operator, m is a particle mass, r and p denote the position and momentum operators respectively. We now write the expression for a partition function -for a quantum system one obtains the partition function by calculating the trace of the density operator ρβ = exp(-β Ĥ):

Z β = tr( ρβ ) ≡ tr(e -β Ĥ ) = dr 0 r 0 |e -β Ĥ |r 0 , (3.2)
where β ≡ 1/(k B T), with k B denoting the Boltzmann's constant and T is the temperature. The trace is taken over any complete basis state, and in our particular case we select a basis of a position operator |r .

Let us look at the expression for ρβ in more details. Operators p and r that enter the Hamiltonian do not commute and neither do any function of these operators, hence we cannot write e β K+β Û = e β K e β Û straightforwardly as it would be the case in the classical limit. Instead, one uses the Trotter's [START_REF] Trotter | On the product of semi-groups of operators[END_REF][START_REF] Suzuki | Generalized Trotter's formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems[END_REF] theorem in the form

e ( Â+ B) = lim M-→∞ 1 + Â + B M M = = lim M-→∞ 1 + Â M 1 + B M M = lim M-→∞ e ( Â/M) e ( B/M) M (3.3)
In numerical calculations one chooses M to be a large positive integer number. The decomposition (3.3) for a finite M is not exact, but the expansion (3.3) is absolutely convergent, so that the extrapolation to infinite M can be done numerically.

In practice one performs a computation by splitting the interval [0, β] into M slices of equal length ∆τ = β/M. Naturally, one can generalize this separation to the intervals of varying size, provided the constraint ∑ M p ∆ p τ = β is satisfied, which is, however, not done in the following and we employ the former separation. Up to the lowest order in ∆τ one can then write the density operator [START_REF] Chandler | Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids[END_REF] as

e -β Ĥ = (e -∆τ Ĥ ) M = (e -∆τ p2 2m e -∆τ Û(r) ) M + O(∆τ) (3.4)
Depending on the interaction potential, the calculation of error made by this approximation might be quite involved [START_REF] Lloyd | Universal Quantum Simulators[END_REF][START_REF] Childs | Theory of trotter error with commutator scaling[END_REF]. We then make use of eigenvectors for operators of potential and kinetic energy |r and |p respectively:

U(r)|r = U(r i )|r p2 2m |p = p 2 2m |p
Each of sets |r and |p forms a complete basis, so that one can write identities

1 = dr|r r| , 1 = 1 2π dp|p p| (3.5)
Let us introduce the first set of identities into the expression of Z β . Acting with the operators that enter the partition function on these states, we recast the expression

3.1. Path Integral formalism 21 (3.2) into Z β = dr 0 r 0 |e -β Ĥ |r 0 = lim M-→∞ dr 0 dr 1 dr 2 dr 3 . . . dr M-1 r 0 |e -∆τ K e -∆τ Û |r 1 r 1 |e -∆τ K e -∆τ Û |r 2 r 2 |e -∆τ K e -∆τ Û |r 3 . . . r M-1 |e -∆τ K e -∆τ Û |r 0
In the following in order to improve readability and to simplify formulas we will denote the terms that enter the expression

ρ(r 0 , r 1 , ∆τ) = r 0 |e -∆τ Ĥ |r 1
Using these notations we write the expression for Z(β)

Z β = dr 0 dr 1 dr 2 . . . dr M-1 ρ(r 0 , r 1 , ∆τ)ρ(r 1 , r 2 , ∆τ) . . . ρ(r M-1 , r 0 , ∆τ) (3.6)
We see that the partition function decomposes into a product of density matrix elements ρ(r i , r i+1 , ∆τ), and each of these elements can be calculated independently.

In order to compute matrix element ρ(r 1 , r 2 , ∆τ) we use the second set of identities from Eq.(3.5). After inserting it into expression for ρ(r 1 , r 2 , ∆τ) we get

ρ(r 1 , r 2 , ∆τ) = r 1 |e -∆τ K e -∆τ Û |r 2 = 1 (2π) dp 1 r 1 |e -∆τ K |p 1 p 1 |e -∆τ Û |r 2
Next we remember the relationship between vectors of different basis |r and |p , which has the form r|p = e -ipr and write:

ρ(r 1 , r 2 , ∆τ) = 1 (2π)
dp 1 e -ip 1 r 1 e -∆τ p 2 1 /2m e ip 1 r 2 e -∆τU(r 2 )

The potential energy is independent of momentum, which means that we can perform the integration over p 1 and get the expression for a matrix element:

ρ(r 1 , r 2 , ∆τ) = m 2πh 2 ∆τ 1/2 e -m (r 2 -r 1 ) 2 2h 2 ∆τ e -∆τU(r 2 ) (3.7)
Let us write the resulting formula for a partition function explicitly:

Z β = m 2πh 2 ∆τ M/2 dr 0 dr 1 . . . dr M-1 e -m (r 1 -r 0 ) 2 2h 2 ∆τ e -∆τU(r 1 ) e -m (r 2 -r 1 ) 2 2h 2 ∆τ e -τU(r 2 ) e -m (r 3 -r 2 ) 2 2h 2 ∆τ e -∆τU(r 3 ) . . . e -m (r 0 -r M-1 ) 2 2h 2 ∆τ e -∆τU(r 0 ) (3.8)
In the schematic depiction of this formula on the right the interaction part of each matrix element is calculated within a slice, and the kinetic part is represented as springs that connect same particles on different slices. Slice M + 1 is associated with slice 0 (shown in green). The utility of the expression (3.8) lies in the fact that we managed to translate calculation of purely quantum quantity (3.2) into "classical" computation of expression (3.8) which is suitable for the numerical calculations. However, the price to pay is the necessity to perform multivariable integrals.

Thermodynamic average of an observable

In the similar spirit we can calculate thermodynamic average of any observable. In quantum mechanics the average is calculated as

 = 1 Z β tr Âe -β Ĥ = 1 Z dr 0 r 0 | Âe -β Ĥ |r 0 (3.9)
In order to apply path integral formalism we need to introduce Trotter separation for density matrix. For simplicity we consider here the observable  that depends on the coordinates, however the calculations can be applied to momentadependent quantities in the analogous way. The average is calculated by considering action with the operator  on a matrix element ρ(r i , r i+1 , ∆τ) in the decomposition of partition function. Due to cyclicity of the trace operator we can act with  on any matrix element.

(k + 1)∆τ k∆τ (k -1)∆τ 0 β ∆τ = β/M  The resulting path integral expression is  = 1 Z β dr 0 r 0 | Âe -β Ĥ |r 0 = 1 Z β dr 0 • • • dr M-1 ρ(r 0 , r 1 , ∆τ)A(r 1 )ρ(r 1 , r 2 , ∆τ) . . . ρ(r M-1 , r 0 , ∆τ) (3.10)
Similar to the calculation of partition function in the previous section, we managed to replace computation of trace operator for quantum operators with "classical" computation of integrals. The common way to calculate multidimensional integrals is to use Monte-Carlo approach. We will see in the following sections how it can be performed. Before proceeding with the numerical calculations, however, we need to address important question of computing correlation function between operators at different times.

Correlation function between pair of operators at different imaginary times

We have seen in Section 2.1.1 that in order to calculate heat conductivity we require the correlation function between heat current operators. In the sections above we saw that applying path integral approach we can calculate equilibrium properties.

In order to apply PI formalism we remember the calculation of real time correlation function for two quantum operators. First, we write the evolution of operator  using Heisenberg representation:

Â(t) = e i Ĥt/h Â(0)e -i Ĥt/h
Then, the correlation function between operators  and B separated by the time interval t is defined as

C AB (t) = Â(t) B(0) = 1 Z β tr Âe -it Ĥ/h Be -(β+it/h) Ĥ (3.11)
What one can calculate in the path integral formalism is the correlation between operators at different slices. Again, here we consider operators that depend on position of particles, however one can do similar analysis for a general operator. First, we act with the operators on separate slices. Then, associating quantity τ with the "distance" between them, i.e. τ = k∆τ ≡ kβ/M in the example on the right, we can write the PI expression for the correlation as

1 Z β dr 0 dr 1 • • • dr M-1 ρ(r 0 , r 1 , ∆τ) . . . ρ(r k-1 , r k , ∆τ)A(r k )ρ(r k , r k+1 , ∆τ) . . . ρ(r M-1 , r 0 , ∆τ)B(r 0 ) = 1 Z β tr[ Âe -τ Ĥ Be -(β-τ) Ĥ ] ≡ Â(τ) B(0) (3.12) τ k+1 τ k τ k-1 0 β ÂB τ = k∆τ
Comparing this expression with the previous one (3.11) we observe that we can now translate one equation into another by replacing real time t with imaginary time τ by performing a Wick rotation [START_REF] Wick | Properties of Bethe-Salpeter wave functions[END_REF] and replacing the real time t with an imaginary time t = iτh.

Alternatively, one can observe that the reason to work with imaginary time lies in the fact that imaginary part for the statistical weight of the matrix element, which makes it not suitable for usual numerical sampling methods [START_REF] Makri | Monte carlo integration with oscillatory integrands: implications for feynman path integration in real time[END_REF].

Numerical methods

In the previous section we have established the method that allows to represent the density operator as a path integral. The later has the advantage, that it can be explicitly written down in the position representation in form of an integration over paths, without having to diagonalize the Hamiltonian. In the following we show how, at least in principle, one can calculate thermodynamic averages and time correlation functions for any given system and observables. In practice, however, the highly dimensional integral that one ends up with after writing the expression for the correlation functions (even with momenta degrees of freedom integrated out) are often quite elaborated to tackle down exactly. There are several possibilities to get around this issue -either to use various approximations (as for example [START_REF] Erkoç | Path-integral solution for pseudoharmonic potential[END_REF]) or to resort to numerical methods. For the purpose of this work -namely, to obtain imaginary time correlations exactly and control the accuracy of the computation -we adopt the numerical approach in the form of Monte-Carlo path integral [START_REF] Binder | Monte Carlo Simulation in Statistical Physics: An Introduction[END_REF]. The rest of this chapter is constructed in the following way. First, we provide a short introduction to a general Monte-Carlo calculation of the average quantities and corresponding errors. Then, we discuss several expressions for the numerical representations of correlation functions that involve momentum operator and compare their effect on the accuracy of the result.

Monte Carlo integration

Monte-Carlo methods (which we will occasionally refer to as MC for simplicity) have a great variety of applications, however, originally MC calculations were developed as a set of techniques for evaluating integrals (namely in application to nuclear physics). The general idea can be illustrated by the following simple example [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]. Let us consider an integral of a single variable defined on the interval [a, b]: I = b a f (x)dx where f (x) is some function (not necessarily regular). In order to evaluate such integral one can employ different strategies, and Monte-Carlo approaches tackle this problem from the stochastic point of view. Let us take at random N points x i uniformly distributed on [a, b) and then calculate the quantity

I est = b -a N N ∑ i=0 f (x i )
It is clear from the definition, that in the limit of N → ∞ the quantity I est , which we will call estimator, is equal to I. However, depending on the function and the chosen range of variables, the convergence rate of the estimator to the actual value as a function of sample points N might be quite slow (for example, see illustrations by Figs.

3.1). For a function in our illustration it is then reasonable to divide interval [a, b) into smaller intervals and evaluate integral I on each of these intervals separately:

I est = c -a N 1 N 1 ∑ i f (x i ) + b -c N 2 N 2 ∑ j f (y i )
where points {x i }, {y j } were plucked from intervals [a, c), [c, b) respectively. This is an example of importance sampling [START_REF] Kloek | Bayesian Estimates of Equation System Parameters: An Application of Integration by Monte Carlo[END_REF]. In general, the importance sampling is the strategy to reduce the variance of the Monte-Carlo measurement. The idea behind the method is quite straightforward. One assumes that FIGURE 3.1: Calculation of integral certain ranges of simulated quantity have more importance for the output than the others and finds such distribution that would bias the measurements in order to emphasize these ranges in the sampling. Naturally, this would bias the result if it is applied directly to the calculation of quantities. However, by weighting the simulation output to account for the bias, one ensures that the importance sampling estimator is unbiased. The weight is given by the likelihood ratio between true underlying distribution with respect to the introduced biased distribution.

Let us start with calculation of the average value. Consider N statistically independent random samples X = x 1 , x 2 , . . . , x N obtained from the distribution P X . In order to calculate the average of function g(X) defined for values of x one writes

g(X) ≡ ∑ x g(x)p X (x) (3.13)
In the computer simulations one has only limited access to the "pure" random numbers -for example, in the form of thermal noise [START_REF] Gabriel | A generator for unique quantum random numbers based on vacuum states[END_REF][START_REF] Bierhorst | Experimentally generated randomness certified by the impossibility of superluminal signals[END_REF]. However, at the present the collection of such data is quite slow, so instead one typically employs for any practical purpose a set of pseudo-random numbers constructed via numerical algorithm from initial ("seed") number which can be taken to be "true" random number [START_REF] Niederreiter | Study of the Homology between Silicon and Germanium by Thermal-Neutron Spectrometry[END_REF].

The most common such algorithm currently used is the Mersenne twister [START_REF] Matsumoto | Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator[END_REF]].

Let us now define another random variables Y with a probability p Y (y) such that p Y (x) > 0 when p X (x) > 0.

Then, the average g(X) can be cast as

g(X) = ∑ x g(x)p X (x) = ∑ x g(x) p X (x) p Y (x) p Y (x) = = ∑ y g(y) p X (y) p Y (y) p Y (y) = p X (Y) p Y (Y) g(Y) (3.14)
In deriving last equation we took into account that p X (y) = 0 for elements y that do not belong to the space of x.

One can immediately apply this expression to the computation of average values by Monte-Carlo integration. In the numerical calculation one generates N individual random variables x1, . . . , x N . Then, in order to evaluate the mean value of a certain function g(x), one uses the sample mean:

g N = 1 N N ∑ i=1 g(x i ) (3.15)
Taking the expression above, one can find alternative approximation for g N by sampling N values y 1 , . . . y N from another distribution. The mean value then can be Chapter 3. Calculation of imaginary time correlation functions estimated by

g N = 1 N N ∑ i=1 p X (Y i ) p Y (Y i ) g(Y i ) (3.16)

Metropolis algorithm

In the practical computation of correlation functions the normalization factor (which is the partition function Z) is usually unknown. The calculation of Z poses at least the same difficulties as the the calculation of correlations themselves. Furthermore, the partition function is not required for the study of transport properties. In order to sample the configuration space via Monte-Carlo methods, one typically adopts dynamic process: one prescribes the procedure that the system follows to evolve from one configuration to another. In general this procedure does not correspond to the actual time evolution of a system: the purpose of this procedure is merely to efficiently generate states of the system according to a weight proportional to p i . The advantage of this approach is that the transition from state i to state j can be generated knowing only with the ratio p j /p i instead of p i , thus eliminating the need to know the partition function Z. The resulting chain of states is called Markov chain [START_REF] Feller | An introduction to probability theory and its applications[END_REF].

Let us introduce the transition probability T ij ≡ T(x i → x j ) -probability of moving system from x i to x j . This probability depends only on the current state and the previous dynamics does not influence it. In order that the resulting states are indeed distributed according to p i we must impose several conditions on T ij [START_REF] Diaconis | The markov chain monte carlo revolution[END_REF]:

-Detailed balance
The transition probability should satisfy the relation

p i T ij = p j T ji ∀i, j (3.17)
-Ergodicity It must be possible to reach any configuration starting from a given one in a finite number of steps.

Together with normalization of T ij , ∑ j T ij = 1, these conditions ensure that the system in state i already distributed with a probability p i actually evolves into state j distributed with probability p j

∑ i p i T ij = p j ∑ i T ji = p j
There are however many ways to construct a Markov chain. One of the vastly used approaches is the Metropolis method. The algorithm works in a following way. First, we need to make initial configuration: we assign to parameters of the system (position of the particles, for example) some values and calculate the total energy of the system E 0 . Then we propose a simple move: change value of one or several parameters and recalculate energy -E j . If E j < E i then we accept the move. If E j > E i then we either accept the move with probability P = e -β(E j -E i ) or reject it. In general, if the transition matrix is not symmetric, i.e. T ij = T ji , the acceptance probability is written as P = min 1,

p j T ji p i T ij .
After that we propose another step and so the algorithm proceeds. One can check explicitly the detailed balance for this algorithm.

Supposing that p j T ji ≤ p i T ij , one can write probability of moving the system from state i to j as

Π ij = T ij min 1, p j T ji p i T ij = p j p i T ji
The probability of the inverse process is given by

Π ji = T ji min 1, p i T ij p j T ji = T ji
Combining these expressions one gets

p i Π ij = p i p j p i T ji = p j T ji = p j Π ji (3.18)
which demonstrates that the Metropolis' rule satisfies detailed balance.

Estimating of error

Monte Carlo calculation is quite often referred to as a "numerical experiment". Similar to a real-life experiment we need to perform a statistical analysis on the results of our simulation in order to define errors that correspond to our observables [START_REF] Krauth | Statistical Mechanis: Algorithms and Computations[END_REF]. For a set of N uncorrelated measurements of O we can write the error as:

δ O est = σ 2 N = ∑ N i=1 (O i -O est ) 2 N(N -1) , (3.19)
provided that O i is gaussian distributed. However, data that we get from Monte Carlo calculations are highly correlated due to the nature of the algorithm for Markov chain. One can use autocorrelation function to quantify the degree of correlation:

φ(t) = O t+i O i -O i 2 O 2 i -O i 2 , (3.20)
where we denoted O i and O t+i -measurements of O taken after i and t + i MC steps respectively. For large t, φ(t) is proportional to exp(-t/T). We start a simulation from a given configuration and then the system relaxes to a global energy minimum and the observables fluctuate around their average values. Initial state and correlation time T play an important role at how long we have to run the simulation before taking samples.

Binning analysis of error

Another issue associated with large correlation time is that we can no longer use simple expression (3.19) while evaluating the errors. There exist several recipes to take the correlation between data points into account [START_REF] Geyer | Practical markov chain monte carlo[END_REF][START_REF] Tibshirani | An introduction to the bootstrap[END_REF][START_REF] Koehler | On the Assessment of Monte Carlo Error in Simulation-Based Statistical Analyses[END_REF]]. We will focus on binning analysis. The set of samples N that we get from MC simulation is divided into N b batches of size L without overlap. Then we calculate block averages:

O b k = 1 L L ∑ i=1 O i+(k-1)L , k = 1, .., N b (3.21)
Chapter 3. Calculation of imaginary time correlation functions

Providing that these block averages are uncorrelated we can then estimate the average

O Ō = 1 N b N b ∑ k=1 O b k (3.22)
If the distribution of O i has a finite variance, then the distribution of the block averaged values is guaranteed to approach a Gaussian distribution for large L by the central limit theorem. Therefore the error can be estimated by the variance of the block averaged data using the formula (3.19):

δ O est = σ 2 b N b = 1 N b (N b -1) N b ∑ k=1 ( O b k -Ō) 2 (3.23)
The practical question -how to determine the size L of the block? For small values of L block averages O b k are still strongly correlated. Increasing size of a block we can get more uncorrelated points, however we can construct fewer blocks N b . To avoid the effect of separation into finite bins, on practice one typically performs binning analysis repeatedly: after calculating block averages on smaller block bins, we treat these values as a new data set and perform this procedure several times. After that we reach true error and repeating this analysis further does not change it.

3.3. Calculation of PIMC estimators for a correlation functions that involve momenta 29

Calculation of PIMC estimators for a correlation functions that involve momenta

Computation of transport coefficients typically necessitates use of correlation functions that involve momentum operator, the simplest example of such function is just a correlation between momenta C pp (τ) = p(τ) p(0) . In the path integral representation one can use different forms of the momentum operator with varying statistical accuracy. In the derivation of Eq.(??) in the Section 3.1.1 we used a so called primitive approximation for momentum operator. Starting from this equation, we can then apply the momentum operator to the density matrix. The resulting estimator is well known classical form, "mass times velocity", i.e. the difference of coordinates of consecutive imaginary times divided by the imaginary time difference.

C pp (τ k ) = - 1 ∆τ 2 (r k+1 -r k )(r 1 -r 0 ) ,
where we denoted by r k ≡ r(τ k ) value of r on the k-th time slice with

τ k = k∆τ = k β M
being discretized imaginary time. The formula is valid for τ = 0, whilst for the zerotime correlation function we write

C pp (0) = 1/(2β) - 1 ∆τ 2 (r 1 -r 0 ) 2
In the next part we will consider a more general form of the correlation function, namely correlations that include position operators.

"Naive" estimator for a current correlation function

Let us study a more complicated quantity, such as p f (r), and calculate its autocorrelation function. We consider f (r) to be dependent only on coordinate operator. In the derivation we shall follow the same ideas as in the Section 3.1.3 and make use of the primitive approximation for the momentum operator.

We start with a general formulation of the correlation function p(t)F(r(t)) p(0)F(r(0)) as a path integral:

p(t)F(r(t)) p(0)F(r(0)) = 1 Z dr r|e -(β-τ) Ĥ p(0)F(r(0))e -τ Ĥ p(0)F(r(0))|r
(3.24) Inserting identities with position and momenta basis states from Eq.(3.5) we get

p(t)F(r(t)) p(0)F(r(0)) = 1 Z dRdp 1 dp 2 . . . r 0 |e -(β-τ) Ĥ |r i p 1 F(r j ) r j |e -τ Ĥ |r k e -ip 2 r k /h p 2 e ip 2 r 0 /h F(r 0 )
where dR = dr 0 dr i dr j dr k . . . By using the shorthand notation

r 1 |e -τ Ĥ |r 2 ≡ ρ(r 1 , r 2 ; τ)
Chapter 3. Calculation of imaginary time correlation functions we can rewrite the expression in the form:

p(t)F(r(t)) p(0)F(r(0)) = 1 Z dRdp 1 dp 2 ρ(r 0 , r i ; β -τ)p 1 e -ip 1 (r i -r j )/h F(r j ) ρ(r j , r k ; τ)p 2 e -ip 2 (r k -r 0 )/h F(r 0 )
Treating p 1 as -i ∂ ∂r j e ip 1 (r j -r i )/h we can write:

p(t)F(r(t)) p(0)F(r(0)) = 1 Z dRdp 1 dp 2 ρ(r 0 , r i ; β -τ)(-i) ∂ ∂r j e ip 1 (r j -r i )/h F(r j )ρ(r j , r k ; τ)(-i) ∂ ∂r 0 e ip 2 (r 0 -r k )/h F(r 0 ) = = 1 Z dRdp 2 (ρ(r 0 , r j ; β -τ)F (r j )ρ(r j , r k ; τ) ∂ ∂r 0 e ip 2 (r 0 -r k )/h F(r 0 )+ + ρ(r 0 , r j ; β -τ)F(r j ) ∂ ∂r j ρ(r j , r k ; τ) ∂ ∂r 0 e ip 2 (r 0 -r k )/h F(r 0 )) (3.25)
Then taking the derivative we finally obtain

p(t)F(r(t)) p(0)F(r(0)) = 1 Z dRdp 2 (ρ(r 0 , r j ; β -τ)F (r j )ρ(r j , r k ; τ) ∂ ∂r 0 e ip 2 (r 0 -r k )/h F(r 0 )+ + ρ(r 0 , r j ; β -τ)F(r j )ρ(r j , r k ; τ) m h r j+1 -r j ∆τ ∂ ∂r 0 e ip 2 (r 0 -r k )/h F(r 0 )) = = - 1 Z dR(ρ(r 0 , r j ; β -τ) m h r 1 -r 0 ∆τ F (r j )ρ(r j , r 0 ; τ)F(r 0 )+ + ρ(r 0 , r j ; β -τ)F (r j )ρ(r j , r 0 ; τ)F (r 0 )+ + ρ(r 0 , r j ; β -τ)F(r j ) m h r 1 -r 0 ∆τ ρ(r j , r 0 ; τ) m h r j+1 -r j ∆τ F(r 0 )+ + ρ(r 0 , r j ; β -τ)F(r j )ρ(r j , r 0 ; τ) m h r j+1 -r j ∆τ F (0)) (3.26)
We note that p f (r) is not hermitian, so that it cannot represent an observable. One of the possible observables is the symmetric combination C {pF} (τ) of operators p and f (r),

C {pF} (τ) = 1 4 p(τ)F(r(τ)) p(0)F(r(0)) + p(τ)F(r(τ))F(r(0)) p(0)+ + F(r(τ)) p(τ)F(r(0)) p(0) + F(r(τ)) p(τ) p(0)F(r(0)) (3.27)
Computation of other terms from C {pF} (τ k ) goes in a similar way as the one above, the difference being the order of the function involved in the calculation.

Calculation of PIMC estimators for a correlation functions that involve momenta 31

In the end we get

C {pF} (τ k ) = 1 Z dr 0 dr j (ρ(r 0 , r j ; β -τ)F(r j ) - r 1 -r 0 ∆τ ρ(r j , r 0 ; τ) m 2 h2 r j+1 -r j ∆τ F(r 0 )+ + 1 2 ρ(r 0 , r j ; β -τ)F(r j )ρ(r j , r 0 ; τ) m h r j+1 -r j ∆τ F (r 0 )+ + 1 2 ρ(r 0 , r j ; β -τ) m h r 1 -r 0 ∆τ F (r j )ρ(r j , r 0 ; τ)F(r 0 )- - 1 4 ρ(r 0 , r j ; β -τ)F (r j )ρ(r j , r 0 ; τ)F (r 0 )) (3.28)
For the case when τ = 0 formula for the estimator C {pF} (τ k ) is slightly modified:

C {pF} (τ = 0) = 1 Z dr 0 ρ(r 0 , r 0 ; β) -1 4 r 1 -r 0 ∆τ r 0 -r M-1 ∆τ m 2 h2 F 2 (r 0 ) - 1 2 r 0 -r M-1 ∆τ m h F(r 0 )F (r 0 )+ + 1 2 F(r 0 )F (r 0 ) m h r 1 -r 0 ∆τ - 1 4 F (r 0 )F (r 0 ) (3.29)

Improved estimators for Monte Carlo calculations

In order to calculate C {pF} (τ k ) numerically we will employ MC approach. However, the MC evaluation of C pF (τ k ) is hampered by the fact that, when ∆τ is getting small, relative fluctuations in (x i+1x i ) become large and the variance of the measured observable grows rapidly (in fact it diverges for ∆τ → 0). As the uncertainty δ MC of the MC estimate of an observable A is related to its variance σ 2 A by δ MC ∝ σ A / √ τ sim , one is therefore forced to increase the simulation time, τ sim , in order to achieve a given precision.

This problem was identified early in the development of PIMC, when trying to estimate the atoms kinetic energy, which is ∝ C pp (τ = 0). A solution was proposed in [START_REF] Herman | On path integral Monte-Carlo simulations[END_REF]: instead of directly using the above expression for C pp (τ k ), the integrals entering the correlation function can be rearranged in order to producing a new estimator for C pp (τ k ), with identical average but smaller variance. The new expression, known in the case of the kinetic energy as the "virial estimator", does not depend explicitly on terms of ∆τ divergent in the limit ∆τ → 0, and therefore does not suffer from the diverging variance associated with the "naive" estimator.

Kinetic energy estimator

Let us first recall the derivation of "virial" form for the kinetic energy [D. M. [START_REF] Ceperley | Path integrals in the theory of condensed helium[END_REF]. Later we will use the same approach in order to obtain the estimator for the desired correlation function. We start by considering the following identity

0 = 1 Z M-1 ∑ j=0 dr 0 dr 1 . . . dr M-1 ∇ j r j M-1 ∏ i=0 e -Mm 2βh 2 (r i+1 -r i ) 2 e -β M U(r i ) (3.30)
Chapter 3. Calculation of imaginary time correlation functions assuming that surface terms vanish. Here Z is the usual partition function

Z = dr 0 r 0 |e -βH |r 0 = dr 0 dr 1 . . . dr M-1 M-1 ∏ i=0 e -Mm 2βh 2 (r i+1 -r i ) 2 e -β M U(r i )
Taking the derivative of both sides of the last equation we obtain

1 Z M-1 ∑ j=0 dr 0 dr 1 . . . dr M ∇ j r j M-1 ∏ i=0 e -Mm 2βh 2 (r i+1 -r i ) 2 e -β M U(r i ) = = 3M + 1 Z M-1 ∑ j=0 dr 0 dr 1 . . . dr M x j ∇ j M-1 ∏ i=0 e -Mm 2βh 2 (x i+1 -x i ) 2 e -β M U(r i ) = = 3M - 1 Z M-1 ∑ j=0 dr 0 dr 1 . . . dr M β M r j U (r j ) M-1 ∏ i=0 e -Mm 2βh 2 (r i+1 -r i ) 2 e -β M U(r i ) + 1 Z M-1 ∑ j=0 dr 0 dr 1 . . . dr M Mm βh 2 r j (r j+1 -r j ) - Mm βh 2 r j (r j -r j-1 ) M-1 ∏ i=0 e -Mm 2βh 2 (r i+1 -r i ) 2 e -β M U(r i ) = 0
If we impose periodicity of r: r j+M = r j , we can rearrange terms in brackets into following expression:

M-1 ∑ j=0 r j (r j+1 -r j ) -r j (r j -r j-1 ) = - M-1 ∑ j=0 (r j+1 -r j ) 2
Our integral then can be rewritten as

0 = 3M - 1 Z M-1 ∑ j=0 dr 0 dr 1 . . . dr M Mm βh 2 (r j+1 -r j ) 2 M-1 ∏ i=0 e -Mm 2βh 2 (r i+1 -r i ) 2 e -β M U(r i ) - - β M r j U (r j ) M-1 ∏ i=0 e -Mm 2βh 2 (r i+1 -r i ) 2 e -β M U(r i ) = = 3M - Mm βh 2 M-1 ∑ j=0 r j+1 -r j 2 - β M M-1 ∑ j=0 r j U (r j ) = = 3M - M 2 m βh 2 (r 1 -r 0 ) 2 -β r 0 U (r 0 ) = 0
Thus, we can reexpress the kinetic term that has large fluctuations:

3M 2β - m 2∆τ 2 h2 (r 1 -r 0 ) 2 = 1 2 r 0 U (r 0 ) (3.31)
The main characteristic of new estimator is its independence from the imaginary time discretization, which ensures lower fluctuations at large M. In the next part we will consider a more general form of a correlation function and we will derive similar estimator with reduced variance.

Current correlation function

In this part we show that the strategy used to obtain the virial estimator can be generalized to any correlation function involving the momentum operator. We consider 3.3. Calculation of PIMC estimators for a correlation functions that involve momenta 33 the same correlation function as the one from the previous part which is a building block in the calculation of transport coefficients:

C {pF} (τ) = ( p(τ) F(τ)) s ( p(0) F(0)) s .
Here F(τ) is a shorthand notation for a generic local function F( x(τ)), which in the case of heat transport would be related to the potential energy. The subscript s indicates that the operator product, which represents an observable quantity, is by convention made Hermitian by symmetrizing the operator, as ( p F) s = 1 2 ( p F + F p). Within the primitive approximation we obtained in the previous part the expression for C {pF} (τ k ) (Eq.(3.28)):

C {pF} (τ k ) = - m 2 h2 ∆τ 2 (r k+1 -r k )F(r k )(r 1 -r 0 )F(r 0 ) + m 2h∆τ (r k+1 -r k )F(r k )F (r 0 ) + + m 2h∆τ (r 1 -r 0 )F(r 0 )F (r k ) - 1 4 F (r k )F (r 0 ) , (3.32)
This expression is valid for k ≥ 1, while the case k = 0 requires slightly different expression as was treated by the Eq.(3.29).

The MC calculation of Eq. (3.32) suffers from the same numerical problem as the momentum correlations, the variance of the leading term in 1/∆τ diverging as ∆τ approaches zero. In order to improve the estimator, we generalize the procedure originally used for the kinetic energy calculations (C pp (0)), and obtain a new estimator with reduced variance for general correlation functions. Let us first consider non-zero imaginary times. We start from the term quadratic in 1/∆τ in Eq. (3.32), since it has the strongest dependence on ∆τ. It can be expressed as

m 2 h2 ∆τ 2 F(r k )(r k+1 -r k )F(r 0 )(r 1 -r 0 ) = = m 2 h2 ∆τ 2 Z dr 0 dr 1 • • • dr M F(r k )(r k+1 -r k )F(r 0 )(r 1 -r 0 ) ρ 0 (r 1 -r 0 ; ∆τ) . . . ρ 0 (r M -r M-1 ; ∆τ) exp -∆τ M ∑ j=0 V(x i ) , (3.33)
where ρ 0 (xy; ∆τ) = x|e -∆τ K |y ∼ exp{-m (x-y) 2 2h 2 ∆τ }. We now transform the set of coordinates {r 0 , r i } to {r 0 , y i }, such that y i = r i+1r i . The constraint r M ≡ r 0 is accounted for by introducing a term δ ∑ M-1 i=0 y i , leading to

m 2 h2 ∆τ 2 F(r k )(r k+1 -r k )F(r 0 )(r 1 -r 0 ) = = m 2 h2 ∆τ 2 Z dr 0 dy 0 • • • dy M-1 δ M-1 ∑ i=0 y i F k-1 ∑ i=0 y i + r 0 y k F(r 0 )y 0 ρ 0 (y 0 ; ∆τ) . . . ρ 0 (y M-1 ; ∆τ) exp[-∆τW], (3.34) with W = M-1 ∑ j=0 V j ∑ i=0 y i + r 0 (3.35)
By using the identity:

m h∆τ y k ρ 0 (y k ; ∆τ) = -∂ y k ρ 0 (y k , ∆τ), (3.36)
we can integrate by parts for the integration over y k . Our next step is based on the observation that the derivative of the δ function w. r. t. to y 0 can be distributed over all coordinates, i.e., ∂ y k δ ∑ y j = 1 M ∑ i ∂ y i δ ∑ y j . A second integration by parts over each of the y i variables eventually leads to

m 2 h2 ∆τ 2 F(r k )(r k+1 -r k )F(r 0 )(r 1 -r 0 ) = m h F(r k )(r 1 -r 0 )F(r 0 ) 1 M M-1 ∑ j=1 jV (r j ) - M-1 ∑ j=k+1 V (r j ) - - mk (h∆τM) F (r k )(r 1 -r 0 )F(r 0 ) - m (h∆τM) F(r k )F(r 0 ) . (3.37)
Repeating the procedure for the terms linear in 1 ∆τ , such as the second term in Eq.(3.37), we can write the correlation in a form that does not depend on ∆τ explicitly:

m h∆τ F(r k )(r k+1 -r k )F (r 0 ) = = ∆τ - M-1 ∑ i=k+1 F(r k )V (r i )F (r 0 ) - k ∆τM F (r k )F (r 0 ) + + 1 M F(r k )F (r 0 ) M-1 ∑ i=1 iV (r i ) , (3.38) and m h∆τ F (r k )(r 1 -r 0 )F(r 0 ) = = ∆τ 1 ∆τ F (r k )F(r 0 ) - M-1 ∑ i=1 F (r k )V (r i )F(r 0 ) - - k ∆τM F (r k )F(r 0 ) + 1 M F (r k )F(r 0 ) M-1 ∑ i=1 iV (r i ) (3.39)
When computing with Eq.(3.37), we need to keep in mind that the expression is valid only for τ k = 0. For the case k = 0, we can apply a similar trick finding the "virial" form.

Again, starting with the general expression:

m 2 h2 ∆τ 2 (r 1 -r 0 )F(r 0 ) 2 (r 1 -r 0 ) = = m h2 ∆τ 2 Z dr 0 dy 0 • • • dy M-1 δ( M-1 ∑ i=0 y i )F(r 0 ) 2 y 0 y 0 ρ 0 (y 0 ; ∆τ) . . . ρ 0 (y M-1 ; ∆τ) exp(-∆τW), (3.40)
where as before we denoted

W = M-1 ∑ j=0 V( j ∑ i=0 y i + r 0 ). ( 3
m 2 h2 ∆τ 2 (r 1 -r 0 )F(r 0 ) 2 (r 1 -r 0 ) = = - m h∆τZ dr 0 dy 0 • • • dy M-1 δ( M-1 ∑ i=0 y i )F(r 0 ) 2 y 0 ∂ y 0 ρ 0 (y 0 ; ∆τ) . . . ρ 0 (y M-1 ; ∆τ) exp(-∆τW) = = - m h M-1 ∑ k=1 F(r 0 ) 2 (r 1 -r 0 )V (r k )F(r 0 ) + 1 ∆τ F(r 0 )F(r 0 ) + m h∆τZ dr 0 dy 0 • • • dy M-1 ∂ y k δ( M-1 ∑ i=0 y i )F(r 0 )y 0 F(r 0 )ρ 0 (y 0 ; ∆τ) . . . ρ 0 (y M-1 ; ∆τ) exp(-∆τW) (3.43)
By substituting,

∂ y 0 δ( M-1 ∑ i=0 y i ) = 1 M M-1 ∑ j=0 ∂ y j δ( M-1 ∑ i=0 y i ), (3.44)
we finally obtain,

m 2 h2 ∆τ 2 (r 1 -r 0 )F(r 0 ) 2 (r 1 -r 0 ) = = ( m h∆τ - m h∆τM ) F(r 0 )F(r 0 ) - m h M-1 ∑ k=1 F(r 0 ) 2 (r 1 -r 0 )V (r k )F(r 0 ) + m hM F(r 0 )(r 1 -r 0 )F(r 0 ) M-1 ∑ j=1 jV (r j ) (3.45)
In contrast with the initial expression Eq. (3.32), all terms are now well-defined as ∆τ → 0. We note, however, that the number of terms involved in the first part of Eq. (3.37) increases linearly with M = β/∆τ, so that the gain following our manipulation is not immediately obvious. The argument that Eq. (3.37) indeed leads to a variance reduction is the following: If all the M contributions to the first term were independent, its variance would scale as ∆τ × M, where ∆τ comes from the term |r 1r 0 | , and the factor M accounts for the M contributions in the sum. As the segments in the path are correlated, even if this estimate is only approximate it still indicates that the variance remains finite even for ∆τ → 0. We will explicitly verify the variance reduction numerically for the test case of harmonic oscillator in the following sections.

We note, that the above derivation of the improved estimator can be adapted beyond the use of the primitive approximation. Similar refinements can be expected to work when employing improved actions, as well as within improved sampling schemes, e.g., PIMD methods based on staging or normal modes M. Tuckerman, 2010, as the variance of the estimator is entirely determined by the analytical form of the kinetic energy part of the action.

Path integral Monte Carlo and stochastic reconstruction study of a single quantum harmonic oscillator

In the last chapters, we have briefly introduced the two main methods that we employed in this thesis. PIMC is used to provide unbiased calculations of imaginary time correlation functions from which one expects to obtain the underlying spectral function via a stochastic reconstruction. We have further addressed the main bottlenecks of this approach: the need of very precise data of the imaginary time correlation functions for the spectral reconstruction, and, the remaining sensitivity of the stochastic reconstruction to unavoidable intrinsic parameters, such as the precise representation of the spectral function and the effective temperature due to correlation of the noise. For that, we have shown how to construct variance-reduced estimators involving momentum operators to be used in the PIMC calculation for our ultimative purpose of calculating current-current correlations. Further, we proposed to combine the stochastic reconstruction with a validation procedure which selects the best reconstruction based on its unbiased predictive power. In this chapter we apply both methods, together with our improvements, to study the harmonic oscillator. Although elementary, the quantum oscillator is the standard test in this context as our numerical methods, PIMC and stochastic reconstruction, are generic for many degres of freedom, e.g. spatial dimension and number of particles, and do not assume any specific interaction potential

In this part we focus in details on the classical example of a single quantum harmonic oscillator and analyze computation of various correlation functions as well as their respective spectra. In the next chapter we consider more elaborate models, that are nonetheless based on the single oscillator and thus readily solvable analytically. When studying these models we shall concentrate primarily on the spectral reconstruction, as it poses in these examples more difficult problem.

We start this chapter by analyzing the correlation function for coordinate operator as the simplest example of correlations. We then pass to the correlations of momenta and correlation function for current operator. For these quantities we then compute the respective spectral densities. The simplicity of the model allows a straightforward analytical computation of properties of interest both in real and imaginary times. Thus, the model constitutes a good benchmark test for our numerical tools, particularly for the spectral reconstruction. Furthermore, the example will serve as a preliminary step to more complicated systems such as solids and disordered systems, although we have not addressed the latter one in the thesis.

Chapter 4. Path integral Monte Carlo and stochastic reconstruction study of a single quantum harmonic oscillator

Position correlation function in imaginary time

Analytical expressions

Hamiltonian of the quantum harmonic oscillator has the from:

Ĥ = h2 2m p2 + mω 2 0 2 x2 (4.1)
First of all, we consider a correlation function between position operators x which we write as C xx (τ) = x(τ) x(0) . We notice that position correlation function at t = 0 is just a mean squared displacement x2 . It is related to the average potential energy E pot which can be computed in several different ways, as we will demonstrate. In this part we treat position correlation functions (and related quantities) in an analytical way and then we compare the results with numerical calculations in the next part. We start by writing the expression for the partition function Z h.o. of quantum harmonic oscillator

Z h.o. = ∞ ∑ i=0 e -hβω 0 (n+1/2) = 1 2 sinh hβω 0 2
For the harmonic oscillator one can relate the value of the potential energy to the total energy E pot = 1 2 E tot which in turn can be calculated from the partition function Z:

E pot = - ω 0 2 ∂ log Z h.o. ∂ω 0 = - 1 2 ∂ ∂β log Z h.o. = - 1 2Z h.o. ∂ ∂β Z h.o.
After differentiating Z h.o. with respect to β we obtain the result

E pot = hω 0 2 coth hβω 0 2
However, this expression is valid only in the limit of continuous imaginary time, that is in the limit M → ∞ of the corresponding path integral. In the following we will derive the same quantities calculated within a discretized path integral for finite M to provide a more direct comparison with numerical PIMC evaluations.

Computation of a correlation function within a discretized path integral

We will now derive the analog expressions within path integral formalism outlined in the Section 3.1.1. We start with the general formula for a partition function (Eq.(3.8)) written for the primitive approximation in imaginary time

Z = dx 1 dx 2 . . . dx M M-1 ∏ k=0 e -m 2h 2 ∆τ (x k+1 -x k ) 2 e - mω 2 0 ∆τ 2 x 2 k (4.2)
and reexpress it as a gaussian integral in a matrix form . The formal expression for the partition function in these notations is

Z = dXe -1 2 X T AX 4
Z = (2π) M 2 [det A] -1/2
For a gaussian integral the pair correlation function x(τ k )x(0) ≡ x i x 0 written for a discretized time τ k = k∆τ can be readily calculated as

x k x 0 = (A -1 ) k0 (4.4)
where A -1 is the inverse matrix of A. In general, one performs the inversion of matrix A by numerical methods. Diagonal elements of A -1 , however, can be calculated in a straightforward manner fully analytically. First of all we determine the set of eigenvalues {λ} of the matrix A. This amounts to solving the equation

λx j = ax j + bx j-1 + bx j+1
We look for solutions of the form

x (k) j = e 2iπ kj M
where k ∈ [0, . . . M -1].

The corresponding eigenvalues are

λ (k) = a + 2b cos 2πk M
Since the determinant of a matrix is invariant with respect to orthonormal rotations the determinant can be expressed as a product of the eigenvalues det

A = ∏ k λ (k) : det A = M-1 ∏ k=0 a + 2b cos 2πk M
For the calculation of thermodynamic properties we will also require an expression for log(det A):

log(det A) = M-1 ∑ k=0 log a + 2b cos 2πk M
With this preparations we can finally obtain ∑ j x 2 j = 1 Z tr ∑ j x 2 j e -β Ĥ . By differentiating the partition function Z with respect to β we will get

x 2 j ≡ x 2 0 = 1 M ∑ k x 2 k = 1 M ∂ ∂a log(det A) = = M-1 ∑ k=0 1 mω 2 0 β + 4m M 2 h2 β sin 2 πk M = M/2 ∑ k=-M/2 1 mω 2 0 β + 4m M 2 h2 β sin 2 πk M , (4.5)
Chapter 4. Path integral Monte Carlo and stochastic reconstruction study of a single quantum harmonic oscillator where in the last line we used the invariance of the expression with respect to the translation k → k + M. In the following part we shall obtain the same correlation function by means of the PIMC and by comparison with this formula we will estimate the dependence of the result on the imaginary time discretization. One can already get an intuition for this formula by considering the limit of large number of time slices M (we remind ourselves that the value of β is fixed). In this limit one can approximate the expression by: 1

mω 2 0 β + 4m M 2 h2 β sin 2 πk M ---→ M→∞ const + O(M -2 ) (4.6)

Analytical calculation of the real-time position correlation function

Yet another way to find the position correlation function x(t) x( 0) is to use canonical representation for the position and momenta operators. We go back to the Hamiltonian (Eq.(4.1)) and write it in terms of operators x and p. Then we apply the general formula for a thermodynamic average of a quantum observable Ô = 1 Z tr[ Ôe -β Ĥ ] which we can write in the basis of Hamiltonian operator as

x(t) x(0) = 1 Z ∑ n n| x(t) x(0)|n e -βE n , (4.7)
where E n and |n are the energy eigenstates and eigenvectors respectively. In the Heisenberg representation time evolution for operators is embodied by the transformation of the form x(t) = e i Ĥt/h xe -i Ĥt/h . To simplify slightly the notations we denote as usual x(0) ≡ x. For a quantum oscillator it is possible to introduce creation and annihilation operators â † and â that will facilitate the computation of correlation functions. These new operators are chosen such that momentum and position operators can be expressed as:

x = h 2mω 0 â + â † (4.8) p = i hmω 0 2 â † -â (4.9)
Operators â and â † satisfy commutation relation

[ â, â † ] = 1 (4.10)
Using the operators â and â † one can simplify the expression for the Hamiltonian and write it as :

Ĥ = hω 0 ( â † â + 1 2
).

(4.11)

One can also define N ≡ â â † as the occupation number operator. It will be useful to calculate commutators between Ĥ with â and â † :

[ Ĥ, â] = hω 0 [ â † â, â] = hω 0 ( â † â â -â â † â) = hω 0 [ â † , â] â = -hω 0 â (4.12) [ Ĥ, â † ] = hω 0 â † (4.13)
Now let us consider eigenvectors of the energy operator:

Ĥ|n = E n |n
and derive some of their properties. From commutators between Hamiltonian and creation-annihilation operators we can see a following properties:

Ĥ â|n = (E n -hω 0 ) â|n Ĥ â † |n = (E n + hω 0 ) â † |n
Hence vectors â|n and â † |n belong to the linear space of eigenvectors of Hamilton operator with eigenvalues (E nhω 0 ) and (E n + hω 0 ) respectively. We denote the corresponding eigenvectors as |n -1 and |n -1 . From the commutation relation one can also find the explicit action of operators â and â † on any given state |n . In the end we get:

â|n = √ n|n -1 (4.14) â † |n = √ n + 1|n + 1 (4.15)
Using this formalism one can write an expression for the correlation function x(t) x(0) :

x(t) x(0) = 1 Z ∑ n n|e i Ĥt/h xe -i Ĥt/h x|n e -βE n = = 1 Z h 2mω 0 ∑ n n|e i Ĥt/h xe -i Ĥt/h √ n|n -1 + √ n + 1|n + 1 e -βE n = = 1 Z h 2mω 0 ∑
n e iE n t/h ne -iE n-1 t/h + (n + 1)e -iE n+1 t/h e -βE n (4.16)

For the energy operator of the quantum oscillator (4.11) one can find the eigenvalues explicitly E n = hω 0 (n + 1 2 ) along with partition function Z h.o. = 1 2 sinh(hβω 0 /2) , which we can use to simplify the expression for the correlation function

x(t) x(0) = h 2mω 0 e iω 0 t + h mω 0 1 Z h.o. ∑ n n cos(ω 0 t)e -hω 0 (n+1/2)β
Using the formula

∑ n ne -xn = d dx 1 1 -e -x = -e -x (1 -e -x ) 2 = e x/2 4 sinh(x/2 -cosh(x/2)) sinh 2 (x/2)
we can rewrite the second term in the previous equation and obtain

x(t) x(0) = h 2mω 0 e iω 0 t + h mω 0 cos(ω 0 t)e -hω 0 β/2 2 sinh hω 0 β 2 e hω 0 β/2 4 1 sinh(hω 0 β/2) - cosh(hω 0 β/2) sinh 2 (hω 0 β/2) = = h 2mω 0 (cos(ω 0 t) + i sin(ω 0 t)) - h 2mω 0
1 -cosh(hω 0 β/2) sinh(hω 0 β/2) cos(ω 0 t) (4.17) single quantum harmonic oscillator Since our computer simulations are performed in the imaginary time we can easily update our equation accordingly and make the substitution t = ihτ which results in a following result:

x(τ) x(0) = = h 2mω 0 (cosh(hω 0 τ) -sinh(hω 0 τ)) - h 2mω 0 1 - cosh(hω 0 β/2) sinh(hω 0 β/2) cosh(hω 0 τ) = = h 2mω 0 cosh(hω 0 (τ -β/2)) sinh(hω 0 β/2) (4.18)
The correlation function is explicitly symmetric around τ = β/2. Using the definition of spectral function 0) and comparing Eq.(4.18) with the inversion for imaginary times correlations (2.28) written in the form β-hτ) one observes that the corresponding spectral function has only one non-zero frequency. Thus, using a discrete representation for the spectrum S(ω) = ∑ i A(ω i )δ(ωω i ) performing a spectral reconstruction we should end up with only one coefficient which corresponds to the frequency ω = ω 0 . The value of this coefficient is also defined by the Eq.(4.18) up to corrections of order M -2 .

S xx (ω) = 1 2π ∞ -∞ e -iωt C xx (t)dt with C xx = x(t) x(
C xx (τ) = ∞ 0 dωS xx (ω) e -hωτ + e -ω(

PIMC calculations of the imaginary time correlation function

In this part we present the result for the Monte-Carlo calculation of the position correlation function and compare it with an exact analytical result provided by Eq.(4.18). The dependence of x 2 on a discretization scheme (Fig.(4.1)) demonstrates that, as we have already observed from analyzing Eq.(4.6), the mean squared displacement x 2 converges to the continuous time value (orange line on the plot) as M -2 for large M. The formula (4.5) provides us with an exact expression for this convergence. Plotting numerical results on the same figure we see that PIMC data follow this behavior. These results prove that discretization of imaginary time in Eq.(3.2) and our numerical representation of this formula play an important role in computation of thermodynamical properties of the system. This figure also illustrates how to choose the imaginary time discretization in practice. Namely, ratio ∆τ = β/M needs to be sufficiently small so that the error in the approximation (3.4) becomes negligible in comparison with the statistical error, associated with the Monte-Carlo algorithm itself. This plot also shows us that in order to obtain quite accurate data, that is, such that aiming for an accuracy of order 10 -3 for our PIMC data, it is sufficient to choose β/M ≈ 0.1.

In Figs. 2) in this section we observe that the difference between numerical outputs and analytical results decreases rapidly with ∆τ getting smaller and the accuracy of these data reaches relative error of 10 -3 for comparatively small values of ∆τ ∝ 0.1. From these simulations we also notice that the part of the correlation function that is most susceptible to the stochastic and discretization error is the region around τ = β/2, whilst the region around τ = 0 and τ = β is affected the least. This is a common pattern and we will see the same behavior for all the following computations.

Stochastic reconstruction of the spectral function

We now proceed to the calculation of a spectral density that corresponds to the correlation function C xx (τ) from the previous part (Fig 4 .2). To this end we shall follow the average spectrum approach outlined in chapter (2.2.3). We sample configuration space of spectral densities and then average them with weight e -Θχ 2 /2 , where we define functional χ 2 with inaccurate data coming from the MC simulation of a correlation function.

Since we approximate spectral density with S(ω) = ∑ P i A(ω)δ(ωω i ), in practice sampling of the configuration space involves probing space of coefficients A i ≡ A(ω i ). Due to positivity of a spectral function we know that these coefficients should also be positive A i ≥ 0. Another constraint can be obtained by considering the correlation function at a given fixed time τ k , for simplicity of the argument let us put k = 0. Then it is clear from the definition of the spectral function C(τ) = ∞ 0 dωS(ω) e -hωτ + e -hω(β-τ) that the coefficients A i are not independent and in the case τ = 0 their sum is equal to the value of C(0) within the error σ 2 (0). These restrictions may however result in not sampling the configuration space efficiently. It might then be desirable to slightly weaken these conditions. From our computations we find that a sum rule is reflected in the output data, so that it is not necessary to impose it explicitly. The positivity constraint is more severe, however, and one need to incorporate it explicitly. In order to test the effect of this requirement one can introduce a "soft" cutoff: we include in the calculation of average spectrum such configurations that allow small negative value of coefficients A i . We expect that after averaging we recover positive coefficients and the condition S(ω) ≥ 0 will be satisfied. The reconstructed spectral density for the position correlation function Chapter 4. Path integral Monte Carlo and stochastic reconstruction study of a single quantum harmonic oscillator

Momentum correlation function

Analytical results

Now we would like to obtain the correlation between momenta operators C pp = p(t) p(0) of the quantum harmonic oscillator. Following the canonical formalism that we used in the previous section to derive position correlation function, we again use canonical creation and annihilation operators a † and a for expressing the momentum operator as follows:

p = i hmω 0 2 â † -â
In this formalism the correlation function can be written in the Heisenberg representation as:

p(t) p(0) = 1 Z ∑ n n|e i Ĥt/h pe -i Ĥt/h p|n e -βE n (4.19)
where as before E n and |n are the eigenvalues and eigenvectors of the energy operator Ĥ. We denoted as usual p(0) ≡ p.

Writing the action of the operators p on this states and simplifying the result in the same way as we did in derivation in the previous section, we get:

p(t) p(0) = hmω 0 2Z ∑ n e iE n t/h (n + 1)e -iE n+1 t/h + ne -E n-1 t/h e βE n (4.20)
This formula is similar to the expression for position correlation function (eq.(4.16)).

Putting in the values of energy for the harmonic oscillator and replacing real time with the imaginary one t = ihτ we end up with the resulting formula for momenta correlations :

p(τ) p(0) = hmω 0 2 cosh(hω 0 (τhβ/2)) sinh(hω 0 β/2) = m 2 ω 2 0 x(τ) x(0) (4.21)

One observes the same behaviour with imaginary time as for the position correlation function (Eq.(4.18)). Furthermore, setting m = ω 0 = 1 one gets exactly the same expressions for two quantities. Despite the fact that the analytical expression is the same, this correlation is an interesting benchmark for our numerical calculation. As we shall see in the next section, the path-integral Monte-Carlo estimator that one straightforwardly obtains from the path integral calculation is not the same as for the position correlation function. Whilst both these estimators produce the same average value, their variance differs and the PIMC algorithm produces data with different accuracy, which in turn affects the spectral reconstruction.

Path-integral approximation of the momentum operator

Let us look in detail at the numerical expression that one uses in the computational algorithm of path-integral calculations. In the position representation, the momentum operator simply acts as a derivative operator, such that we can obtain an estimator for the momentum correlation function when applied inside the discretized path-integral expression (see chapter 3):

p(τ) p(0) = 1 Z dx 0 x 0 |e -(β-τ) Ĥ p(0)e -τ Ĥ p|x 0 (4.22)

Momentum correlation function 47

Then we insert into this equation identities from the Eq.(3.5) and obtain p(τ) p(0) = 1 Z dxdp 1 dp 2 x 0 |e -(β-τ) Ĥ |x i e -ip 1 x i /h p 1 e ip 1 x j /h x j |e -τ Ĥ |x k e -ip 2 x k /h p 2 e ip 2 x 0 /h (4.23) To simplify notations we use the notation:

x a |e -τ Ĥ |x b ≡ ρ(x a , x b ; τ) (4.24) Then p(τ) p(0) = 1 Z dx 0 dx i dx j dx k dp 1 dp 2 ρ(x 0 , x i ; β -τ)e ip 1 (x j -x i )/h p 1 ρ(x j , x k ; τ)e ip 2 (x 0 -x k )/h p 2 (4.25)
From this formula we see that momentum can be written as a derivative. For example:

p 1 e ip 1 (x j -x i ) = -ih ∂ ∂x j e ip 1 (x j -x i )/h (4.26) Thus we have p(τ) p(0) = 1 Z dx 0 dx i dx j dx k dp 1 dp 2 ρ(x 0 , x i ; β -τ) (-i)h ∂ ∂x j e ip 1 (x j -x i ) ρ(x j , x k ; τ)(-i)h ∂ ∂x 0 e ip 2 (x 0 -x k ) (4.27)
In the simple approximation of the momentum in the path integral calculation we represented it as a difference of coordinates. Using this, we can take this integral by parts and swap derivative with respect to x 0 p(τ) p(0) = 1 Z dx 0 dx 1 dx 2 . . . dp 1 dp 2 x 1x 0 β mM ρ(x 0 , x 1 ; ∆τ)ρ(x 1 , x 2 ; ∆τ) . . . ρ(x i-1 , x i ; ∆τ) ∂ ∂x j e ip 1 (x j -x i ) ρ(x i , x i+1 ; ∆τ)e ip 1 (x 0 -x k ) . . . (4.28)

In derivinig this equation we took into account that the surface terms vanish. Now we can do the integral over p 1 which results in delta function δ(x 0x k ).

Repeating this procedure with respect to x j we end up with

p(τ) p(0) = 1 Z dx 0 dx 1 dx 2 . . . x 1 -x 0 β mM ρ(x 0 , x 1 ; ∆τ)ρ(x 1 , x 2 ; ∆τ) . . . x j+1 -x j β mM ρ(x i , x i+1 ; ∆τ) . . . ρ(x M-1 , x 0 ; ∆τ) = = x 1 -x 0 ∆τ m x j+1 -x j ∆τ m (4.29)
The case τ = 0 requires some additional calculations.

p(0) p(0) = 1 Z dx 0 x 0 |e -β Ĥ p(0) p|x 0 (4.30)

Chapter 4. Path integral Monte Carlo and stochastic reconstruction study of a single quantum harmonic oscillator

Following the same ideas we just showed we write: p(0) p(0) = 1 Z dx 0 dx i dx j dx k dp 0 dp 1 ρ(x 0 , x j ; β)p 2 1 e ip 1 (x 0 -x j ) (4.31)

To avoid differentiating twice with respect to x 0 we will employ instead mixed derivative

p 2 1 e ip 1 (x 0 -x j ) = ∂ ∂x 0 ∂ ∂x j e ip 1 (x 0 -x j ) (4.32) p(0) p(0) = 1 Z dx 0 dx j dx k dp 0 ρ(x 0 , x j ; β) ∂ ∂x 0 ∂ ∂x j e ip 1 (x 0 -x j ) = = - 1 Z dx 0 dx 1 dx 2 . . . dp 0 x 1 -x 0 β mM ρ(x 0 , x 1 ; ∆τ)ρ(x 1 , x 2 ; ∆τ) . . . ∂ ∂x j e ip 1 (x 0 -x j ) = = - 1 Z dx 0 dx 1 . . . dp 0 x 1 -x 0 β mM x j -x j-1 β mM ρ(x 0 , x 1 ; ∆τ)ρ(x 1 , x 2 ; ∆τ) . . . e ip 1 (x 0 -x j ) - 1 Z dx 0 dx 1 . . . dp 0 x 1 -x 0 β mM β M V (x j )ρ(x 0 , x 1 ; ∆τ)ρ(x 0 , x 1 ; ∆τ) . . . e ip 1 (x 0 -x j ) (4.33)
As we mentioned already in the previous section, these estimators are different, even after setting m = ω 0 = 1, from the ones we used to calculate position correlation function, namely simple x(τ) x(0) = x j x 0 . Yet, it could be brought to the same expression.

Let us go back to Eq.(4.27): p(τ) p(0) = 1 Z dx 0 dx i dx j dx k dp 0 dp 1 ρ(x 0 , x i ; βτ)

(-i) ∂ ∂x j e ip 1 (x j -x i ) ρ(x j , x k ; τ)(-i) ∂ ∂x 0 e ip 2 (x 0 -x k ) = = 1 Z dx d dx 0 e -1 2∆τ (x 1 -x 0 ) 2 e -1 2∆τ (x 2 -x 1 ) 2 . . . d dx j e -1 2∆τ (x j+1 -x j ) 2 .
. . e -∆τV(x 0 )... , (4.34)

where we have swapped the derivative from the δ-function by taking the integral by parts. But instead of differentiating kinetic parts of the density element ρ we take one of this integrals by parts again to swap derivative on the potential part:

p(τ)p(0) = - 1 Z dx (x 1 -x 0 ) ∆τ 1 ∆τ (x j -x j-1 ) -∆τV (x j ) e -S = = 1 ∆τ 2 (x 1 -x 0 )(x j -x j-1 ) -(x 1 -x 0 )V (x j ) (4.35)
where e -S includes all the density matrices factors. Let us consider term x 1 V (x j ) in more detail:

x 1 V (x j ) = - 1 Z dXx 1 1 ∆τ d dx j (e -∆τV(x j
) )e -S j (4.36)

where S j = S -∆τV(x j ). After integration we get

x 1 V (x j ) = 1 ∆τ 2 x 1 (x j-1 -x j + x j+1 -x j ) (4.37)
In order to help ourselves to see a connection with the equation (4.35) more easily we denote x i x 0 := C(i). Then we can write

x 1 V (x j ) = C(j -2) -2C(j -1) + C(j)
If reexpress first term in Eq.(4.35) we get

1 ∆τ 2 (x 1 -x 0 )(x j -x j-1 ) = 2C(j -1) -C(j) -C(j -2)
which cancels exactly against x 1 V (x j ) . We end up with expression for an estimator

p(τ)p(0) = x 0 V (x j ) ≡ x(0)V (x(τ)) (4.38)
which in the case of harmonic oscillator simplifies to

p(τ)p(0) = x 0 x j ≡ x(τ)x(0)
which demonstrates the correspondence between momentum correlation and position correlation at the level of the estimators, which we already observed from canonical computation of these functions. Note that in deriving the last equation we have used for simplification m = ω 0 = 1.

Numerical comparison of the estimators

The example of p(τ) p(0) allows us to compare results for different types of estimators: a simple, or "naive", expression that one obtains directly from path integral calculation such as the one provided by Eq.(4.29) and a "virial" one that is obtained after some recalculations in the starting expression. The example of the second type is the Eq.(4.38). We have already inadvertently made the simulation with a "virial" estimator of the momentum correlation of a harmonic oscillator when we computed x(τ) x(0) . Another example of the virial estimator can be found in the introduction where we obtained a general expression for the estimators of the correlation functions that involve momentum operator (see Eq.(3.28)). Here we can see results for "naive" estimator eq.(4.29) (orange curve) and the calculation made with the estimator from eq.(3.32) written in the form (3.37) adopted to the momenta correlations (plotted in green).

The actual "virial" estimator ( p(τ)p(0) = x(0)V (τ)) produces both the smallest error and the relative one. Naturally, we would be interested in using similar estimators for a more complicated correlation functions. The problem, however, comes from the fact that the cancellation that occurred in the derivation of Eq.(4.38) in case of a more complicated operators, for example of the form p(τ) f (x(τ))p(0) f (x(0))

in which we will be interested in the following, would involve computation of some additional terms with large variance. For this reason we employed slightly different approach in deriving the more general expression of a correlation function in the introduction. The formula contains terms that are proportional to the number of time slices M which potentially can lead to a large variance. From the results presented here, however, we observe that the variance, while larger than the one from the x(τ)x(0) estimator, is considerably lower than the one produced by the naive estimator p(τ)p(0) . Thus, we expect that such estimators will perform significantly better than the straightforward ones, and we shall check this explicitly while calculating current correlation function in the next section. Let us note, that the large variance that Eq.(4.29) produces is due to the M 2 term. Just by replacing it with terms linear in M we already drastically improved the output data, particularly for small imaginary times. ----→ 0 accuracy of the measurements becomes worse, most notably for the naive estimator.

Chapter 4. Path integral Monte Carlo and stochastic reconstruction study of a single quantum harmonic oscillator

Spectral reconstruction

Using the correlation functions C pp (τ) from the previous part we can now compute the corresponding spectral density. As before we follow the scheme of the Section 2.2.3 and sample the space of the spectral configurations. This time we restrict the configuration space to the strictly non-negative values S(ω) ≥ 0. Since we already considered the reconstruction from the virial estimator R i V (R i ) while considering the position correlation function, here we restrict our attention only to the second virial estimator, namely the one from the Eq.(3.28) with F = 0, and the naive one, obtained from the expression (4.29).

We compare two results on the Fig. ). Despite that both data produce a spectra peaked around ω = 1 (here we put ω 0 = 1), these figures demonstrates that for the model under consideration the larger errorbars on the data results in a less defined peak. As we will in the subsequent studies for the models with several peaks it might become difficult to distinguish different frequencies. In these future cases we will resort to a more complicated reconstruction schemes which can described by several degrees of freedom (one can find examples for such reconstruction in literature, for example, see [START_REF] Tsujii | Evaluation and extension of maximum entropy models with inequality constraints[END_REF]).

Current-current correlation function

Analytical calculation of the current-current

Now we would like to employ methods we introduced in previous parts to calculate a more complex correlation function. Since we are still working with a quantum harmonic oscillator with single degree of freedom, we will be able to obtain an exact result which we will use later for the purpose of comparison.

In the following we will be interested in a current response function. In general, such quantities require computation of correlation functions composed by the operators of a form Ĵ ∝ f (x(t)) p, where f (x) is some operator that depends only on the particle position. Due to the non-commutativity of momentum and position operators functions that depend on these operators, such as f (x) also do not commute with one another. In fact, one can write a commutation relationships, for example between f (x) and p: [ f (x), p] = i f (x). Hence, we have that operator Ĵ is not self-adjoint Ĵ † = Ĵ, which means that the operator Ĵ is not suitable to serve as an observable, so that we use the following symmetrized one

ĴS = 1 2 ( p f (x) + f (x) p) (4.39)
We can write the expression for a correlation function for this operator:

C J S J S (t) = ĴS (t) ĴS (0) = 1 Z tr e -(β-it/h) Ĥ ĴS e -it Ĥ/h
ĴS , (4.40) single quantum harmonic oscillator

In order to calculate C J † J (t) for the quantum harmonic oscillator system we shall use canonical formalism. To simplify the calculation we will consider the case f (x) = V(x), where V(x) is a harmonic potential energy operator. Similar to the previous calculation we will make use of canonical ladder operators â and â † introduced earlier, that allows to rewrite operators of position and momentum

x = h 2mω 0 â + â † (4.41) p = i hmω 0 2 â † -â (4.42)
With these operators potential energy will take the form V(0):

V|n = mω 2 0 2 x2 |n = mω 2 0 2 h 2mω 0 x √ n|n -1 + √ n + 1|n + 1 = = hω 0 4 n|n + n(n -1)|n -2 + (n + 1)|n + (n + 1)(n + 2)|n + 2 (4.43)
where we used simplified notation Â(0) ≡ Â.

The computation of the correlations proceed in a straightforward but a tedious way, and here we just mention few steps in its derivation. Let us act with the operator ĴS on a state |n

ĴS |n = 1 2 ( p V + V p)|n = i hmω 0 2 hω 0 4 -n(n -1)(n -2)|n -3 -n √ n|n -1 + (n + 1) √ n + 1|n + 1 + (n + 1)(n + 2)(n + 3)|n + 3 (4.44)
Then we substitute this expression into the correlation function (4.40):

ĴS (t) ĴS = 1 Z h.o. ∑ n n|e -β Ĥ e -i Ĥt ĴS e i Ĥt ĴS |n (4.45)
Writing down explicitly the action of the operators ĴS we get

ĴS (t) ĴS = 1 Z ∑ n n|e -β Ĥ e -i Ĥt ĴS e i Ĥt ĴS |n = = hω 0 4 i hmω 0 2 1 Z ∑ n e -βE n e -iE n t n| ĴS e i Ĥt -n(n -1)(n -2)|n -3 -n √ n|n -1 + + (n + 1) √ n + 1|n + 1 + (n + 1)(n + 2)(n + 3)|n + 3 (4.46)
4.3. Current-current correlation function 55 Let us look separately at terms that enter this expression. We compute as an example the first one.

ĴS (t) n(n -1)(n -2)|n -3 = hω 0 4 i hmω 0 2 n(n -1)(n -2) -(n -3)(n -4)(n -5)|n -6 -(n -3) √ n -3|n -4 + + (n -2) √ n -2|n -2 + (n -2)(n -1)n|n
The only non-zero term in this expression is the |n term

ĴS e i Ĥt n(n -1)(n -2)|n -3 → hω 0 4 i hmω 0 2 (n -2)(n -1)n|n e iE n-3 t
In analogy one can compute all other terms in Eq.(4.46). Combining the results and doing some cancellation we end up with

ĴS (t) ĴS (0) = mh 3 ω 3 0 32 1 Z ∑ n e -βhω 0 (n+1/2) [(n -2)(n -1)ne -3iω 0 t + n 3 e -iω 0 t + (n + 1) 3 e iω 0 t + (n + 3)(n + 2)(n + 1)e 3iω 0 t ] = = mh 3 ω 3 0 32 1 Z ∑ n e -βhω 0 (n+1/2)
(n 3 -3n 2 + 2n)e -3iω 0 t + n 3 e -iω 0 t + (n 3 + 3n 2 + 3n + 1)e iω 0 t + (n 3 + 6n 2 + 11n + 6)e 3iω 0 t

= mh 3 ω 3 0 32 1 Z ∑ n e -βhω 0 (n+1/2)
n 3 (e -3iω 0 t + e -iω 0 t + e iω 0 t + e 3iω 0 t ) + n 2 (-3e -3iω 0 t + 3e iω 0 t + 6e 3iω 0 t )+

+ n(2e -3iω 0 t + 3e iω 0 t + 11e 3iω 0 t ) + e iω 0 t + 6e 3iω 0 t

ĴS (t) ĴS (0) = mh 3 ω 3 0 32 1 Z ∑ n e -βhω 0 (n+1/2)
n 3 (e -3iω 0 t + e -iω 0 t + e iω 0 t + e 3iω 0 t ) + n 2 (-3e -3iω 0 t + 3e iω 0 t + 6e 3iω 0 t )+ + n(2e -3iω 0 t + 3e iω 0 t + 11e 3iω 0 t ) + e iω 0 t + 6e 3iω 0 t (4.47)

We can perform the summation if we remember that n 3 (e -3iω 0 t + e -iω 0 t + e iω 0 t + e 3iω 0 t ) + n 2 (-3e -3iω 0 t + 3e iω 0 t + 6e 3iω 0 t )+ + n(2e -3iω 0 t + 3e iω 0 t + 11e 3iω 0 t ) + e iω 0 t + 6e 3iω 0 t = = mh 3 ω 3 0 32 1 Z e -βhω 0 /2 e -βhω 0 + 4e -2βhω 0 + e -3βhω 0 (1e -βhω 0 ) 4 (e -3iω 0 t + e iω 0 t + e -iω 0 t + e 3iω 0 t )+

+ e -βhω 0 + e -2βhω 0 (1 -e -βhω 0 ) 3 (-3e -3iω 0 t + 3e iω 0 t + 6e 3iω 0 t )+ + e -βhω 0 (1 -e -βhω 0 ) 2 (2e -3iω 0 t + 3e iω 0 t + 11e 3iω 0 t ) + 1 1 -e -βhω 0 (e iω 0 t + 6e 3iω 0 t ) = = mh 3 ω 3 0 32
1 Z e -βhω 0 /2 (1e -βhω 0 ) 4 (e -βhω 0 + 4e -2βhω 0 + e -3βhω 0 )(e -3iω 0 t + e iω 0 t + e -iω 0 t + e 3iω 0 t )+ + (1e -βhω 0 )(e -βhω 0 + e -2βhω 0 )(-3e -3iω 0 t + 3e iω 0 t + 6e 3iω 0 t )+ + (1e -βhω 0 ) 2 e -βhω 0 (2e -3iω 0 t + 3e iω 0 t + 11e 3iω 0 t )+ + (1e -βhω 0 ) 3 (e iω 0 t + 6e 3iω 0 t ) Simplifying the result we get ĴS (t) ĴS (0) = mh 3 ω 3 0 256 1 sinh 3 (βhω 0 /2) 6e 3iω 0 t + (4e -βhω 0 + e -2βhω 0 + 1)e iω 0 t + (e -βhω 0 + 4e -2βhω 0 + e -3βhω 0 )e -iω 0 t + 6e -3βhω 0 e -3iω 0 t (4.48)

For comparison with MC calculation we need to substitute t = ihτ

C pV (τ) = ĴS (τ) ĴS (0) = mh 3 ω 3 0 256 1 sinh 3 (βhω 0 /2) 6e -3ω 0 τ + (4e -βhω 0 + e -2βhω 0 + 1)e -ω 0 τ + (e -βhω 0 + 4e -2βhω 0 + e -3βhω 0 )e ω 0 τ + 6e -3βhω 0 e 3ω 0 τ = = mh 3 ω 3 0 256 1 sinh 3 (βhω 0 /2) 12 cosh( 3βhω 0 2 -3ω 0 τ) + 2(4e -βhω 0 + e -2βhω 0 + 1)e βhω 0 cosh( βhω 0 2 -ω 0 τ) (4.49)
From the expression (4.49) it is clear that our spectral density will have 2 non-zero frequencies: ω = ω 0 and ω = 3ω 0 . We study numerically the effect of different variables that enter this formula in more details in section on the spectral reconstruction.

PIMC calculation of the correlation function

The result of any PIMC calculation is given with uncertainty. The error can be attributed to two general categories: systematic bias and stochastic fluctuations. For a harmonic oscillator, the systematic deviation due to the discretization of the imaginary time ∆τ = β/M can be assessed directly, by comparing the result expected from the path integral approach (which in this case can be obtained exactly) with the analytical expression for the correlation function C pV (τ) above, which corresponds to the continuous limit M → ∞.

In order to calculate the exact expression of the correlation function within the primitive approximation of the discretized path integral, we note again that all the integrals involved in the calculation are Gaussian. By using the discretized representation for the momentum operator, one writes C pV (τ) as a thermodynamic average of products of the variables x. Leaving out for the moments all the prefactors, one can write:

C pV ∝ (x k+1 -x k )x 2 k (x 1 -x 0 )x 2 0
For a gaussian integral, Wick's theorem allows to recast such correlations into products of pair correlation functions x i x j . The pair correlations in turn can be easily accessed from the following consideration. Let us write the pair correlation function in the form x i x j = dXx i x j e -X T AX , where A is a symmetric M × M matrix which incorporates the Hamiltonian. Then the required pair correlation can be calculated as the element of the inverse matrix x i x j = A -1 ij . We can therefore use numerical methods to calculate the matrix elements, as we discussed in Section 4.1.1. The relative difference between the two calculations is single quantum harmonic oscillator illustrated in Fig. 4.9. We observe that, for a sufficiently small value of β/M, the deviation is virtually not affected by a change of β. Also, our case studies below are performed by employing an imaginary time discretization ∆τ = β/M = 0.1ω -1 0 , a choice primarily dictated by the need to control the error associated to the primitive approximation. It also limits, however, the resolution of the imaginary time correlation function and, consequently, that of the reconstructed spectral function, especially at high frequencies. We will comment below how this potential bias can be addressed within the verification process. In general, since the high frequency asymptotics is governed by sum rules, it is often most conveniently dealt with by computing leading terms of the short time Taylor expansion. We next focus on the second source of error affecting the PIMC calculation: limited sampling. Indeed, error bars corresponding to average values are obtained by estimating the variance of the observable, which decreases as τ -1/2 sim , with τ sim the simulation time. For a given τ sim , the quality of the result therefore crucially depends on the variance of the estimator. We illustrate this point in Fig. 4.10 and Fig. 4.11, by comparing calculations for the energy current correlation function, C pV , using the naive estimator, Eq. (3.32), and the improved version of Eq. (3.37). The data of Fig. 4.10 clearly show that the virial estimator leads to a spectacular improvement compared to the naive one, with a statistical error that is now comparable to the systematic one resulting from the discretization. Comparing this figure with Fig. 4.11 we also observe that the variance of the improved estimator less sensitive to the decrease of ∆τ as we predicted from the derivation of Eq. (3.37). We show in blue the calculations the 'primitive' estimator and in red -the improved estimator, both using the same path-integral Monte Carlo data.

Stochastic reconstruction combined with validation

We now use the reconstruction procedure outlined in Sect. 2.2.3 to extract the frequency spectrum for the correlation functions obtained in the previous section. In order to perform a reconstruction one needs both to define the set of parameters that expresses the spectral density in Eq. (2.42) and in the integration measure of Eq. (2.43), and to chose the effective inverse temperature Θ. (To simplify the notation in the considered examples, we set ω 0 = m = h = 1.) In the following, we use a discretized model of the spectral density, which is described as a sum of N ω delta-functions in the ω-space, see Eq. (2.38). Specifically, we consider a regular grid of ω-values defined on the interval [0, 5], with a fixed spacing between points, ∆ω = 5/N ω . In addition, we will consider the possibility of a global shift of the grid by δω < ∆ω. Unless specified otherwise, δω = 0, and we fix the origin of the grid in ω = 0. The exact expression for the time correlation function, Eq. (4.49), implies that C pV (τ) decays exponentially with τ in the interval [0, β/2], with a decay rate O(1). Larger values of β therefore lead to a larger amplitude in the decay, with the consequence that the contribution of different frequencies can be more easily resolved for larger β's. In short, a correlation function of the form [exp(-τ) + exp(-3τ)] will be hard to distinguish from 2 exp(-2τ) if data are only available in the interval [0, 1]. Resolving the two frequencies ω 1 = 1 and ω 2 = 3 is therefore essentially impossible if β/2 < 1.

In order to illustrate this point, we calculate and analyze the spectral function for the energy current correlation functions at the two inverse temperatures β = 3 and 10, with an imaginary time discretization ∆τ = 0.1. With this value of ∆τ, the does not further improve the result. We remark that in this Section the spectral reconstruction has been based on discretized imaginary time correlation functions. For convenience, the discretization usually coincides with the imaginary time step controlling the Trotter error of the path integral. In order to estimate the influence and potential bias of the discretization on the spectral reconstruction, one can perform the verification step involving subsets of C (τ k ) at no additional cost. For the test case of the double well potential discussed later, we have investigated this point explicitly. In particular, we have observed that the verification is not qualitatively affected by the change of data discretization for a reasonable range of ∆τ, apart from a shift of the minimum in Θ. The filled rectangles are centered at the positions of the two deltafunctions of the exact spectrum for ∆ω = 1, with an area corresponding to their respective weights.

Case studies for stochastic reconstruction

In our previous study of the harmonic oscillator in chapter 4, we have addressed two conceptually different difficulties. First, we have used the example of the oscillator to study and improve the convergence of various correlation functions calculated by Path-Integral Monte Carlo methods. Then, we have benchmarked the stochastic reconstruction of the spectral function from the actual imaginary time correlations.

In the following, we focus on the second aspect and test our reconstruction method combined with validation on models having spectral functions qualitatively different from the single oscillator: a double well potential and a continuous density of oscillators. For that, we do not require to perform a full Path-Integral Monte Carlo calculation. These models are still sufficiently simple to be solved exactly by other means. In order to analyze realistic cases, we then simply create an artificial noise on top of the exactly calculated correlation functions. This allows us to concentrate only on the stochastic reconstruction scheme we developed. We stress that to our knowledge, the validation procedure on top of the stochastic reconstruction has never been proposed and studied in the present context. Our results are encouraging that the validation allows for eliminating the remaining influence of explicit and implicit parameters entering into the reconstruction.

Double well potential

Exact analytical spectrum

In order to illustrate that our discussion is not intrinsically limited to harmonic oscillators, we examine an additional classical benchmark example [Perez, M. E. Tuckerman, and Muser, 2009] with strongly anharmonic features. Let us consider the Hamiltonian describing a particle trapped in a double well potential,

Ĥ = p2 2m - α 2 x2 + γ 4 x4 , (5.1)
and we choose α = 1 and γ = 1. We can obtain the energy spectrum and the correlation functions of interest by numerical diagonalisation of Ĥ. We show the resulting potential along with first few energy levels on the Fig.

(5.1). In the following we consider the position correlation function,

C β xx (τ) = x(τ) x(0) . C xx (τ) = x(τ) x(0) = 1 Z ∑ n,m e -E m τ e -(β-τ)E n | n| x|m | 2 (5.2)
The corresponding spectral function can be expressed analytically as,

S(ω) = 1 Z ∑ n,m e -βE n | n| x|m | 2 δ(ω -E m + E n ). (5.3) 2 1 0 1 2 x 0.0 0.5 1.0 V V min E 0 = 0.147 E 1 = 0.872 E 2 =

Spectral reconstruction and validation

Starting from the exact correlation function, we generate an extended data set C(τ k ) with an artificial Gaussian noise of variance σ k = 10 -3 × Cβ xx (τ k ), which is comparable to the typical stochastic error of the MC calculation that we performed for a single harmonic oscillator. We next employ the calculation scheme detailed in Sect. 2.2.3, and apply the χ 2 val validation criteria in order to determine the optimal reconstruction. As discussed in the previous section, we are interested in the effect of the modification of the simulation hyper-parameters on the quality of the spectral reconstruction.

As above, we focus in particular on the effective inverse temperature, Θ, the number, N ω , of δ-functions considered in the ω-interval [0, 5], and the uniform shift, δω. Note that, contrary to the case of the harmonic oscillator, the spectral function is now expressed in terms of a set of non-equally-spaced δ-peaks (see Table 5.2). As a consequence, none of the nodes of the uniform ω-grid we consider in our calculations coincides with the ω i . This is an intentionally non-optimal choice which, however, allows us to further illustrate important features of the reconstruction and associated validation, and provides hints towards possible improvements.

We first consider the inverse temperature β = 8, and choose ∆τ = 0.1. At this low temperature, all coefficients at frequencies larger than ω 1 = E 1 -E 0 in Eq.( 5.3) are suppressed exponentially, and the spectrum practically consists of a single δpeak. As for any discrete spectrum, the reconstruction becomes quite sensitive to the positions of the frequencies ω i if the grid is sparse. One can therefore study in detail the accuracy of the reconstruction when modifying the shift, δω, and the distance between the grid points. Here, however, we do not perform such a detailed analysis, and for the sake of illustration we consider only a few different models with a regular grid shifted by a constant amount.

In Fig. 5.2 we show the resulting spectral functions at the indicated values of N ω and δω, while Fig. 5.3 displays the corresponding χ 2 val . The data for N ω = 10 and δω = 0 correspond to values of χ 2 val approximately ten times larger than those associated with the shifted lattice case, δω = 0.25, and have been thus omitted. If we consider the inaccuracy in the alignment of the ω-grid with the exact peak positions, it comes as no surprise that the model with N ω = 10 and δ = 0.25 provides a worse reconstruction than what we have obtained with the non-shifted N ω = 25 grid, despite visually resembling more closely the exact spectrum. Reconstruction of a single peak spectrum is performed essentially with only a few coefficients A(ω i ) that are closest to ω 1 , and because with denser grids we have more fitting parameters at our disposal, it is natural that one obtains a better validation with N ω = 25.

We now consider the more challenging case of the correlation function C β=1 xx at β = 1 (calculated with ∆τ = 0.02), whose spectrum displays several peaks (Fig. 5.4). We perform the reconstruction with the same parameters as in the case above, the resulting values of χ 2 of the reconstruction are plotted in Fig. 5.5. We observe that the use of sparse grids (even "inaccurately" placed) proves advantageous in this case. Both models with N ω =10 perform considerably better than in the previous case, whilst the shifted lattice shows the best fit among considered models at Θ = 10 (Fig. 5.6). This figure also clearly demonstrates the effect of the overfitting, which manifests itself in the plot as the increase of the validation χ 2 val after some value of Θ for each configuration of parameters.

Another practical issue that is interesting to address is the influence of the imaginary time discretization of the test set of correlation function. Using sparser data set for calculation of χ 2 val we do not observe qualitative effect on the result. Value of Θ, at which the best fit (minimum of χ 2 val ) is achieved, however, becomes smaller with less accurate test data.

Continuum distribution of oscillators

Calculation of the spectrum

We now move to our next test model, and study the potential energy current correlation function of a system containing a large number of independent, non interacting harmonic oscillators. Considering the C pV of Eq. (4.49) as a function of ω 0 , the correlation function for an ensemble of oscillators with a continuum of frequencies can be written as,

C cont pV (τ) = ω cut 0 dω 0 C exact pV (τ; ω 0 )g(ω 0 ). (5.4)
The form of the density of states, g(ω 0 ), and the value of the frequency cutoff, ω cut , are arbitrary. In the following we consider a Debye-like g(ω 0 ) ∝ ω 2 0 , with ω cut = 1, and fix β = 10. With this choice, the exact spectrum for the energy current correlation is a superposition of two functions with a compact support, assuming non zero values in the range [0, ω cut ] and [0, 3ω cut ], respectively. As a result, it will display two sharp discontinuities, at ω cut and 3ω cut , respectively. As we already done in our study of the double well potential, here we do not generate the data by Monte Carlo simulation, but we rather employ the exact analytical expression, subsequently adding a Gaussian random noise with a variance proportional to the data themselves, σ k = 10 -2 × C cont pV (τ k ). This variance is also used as the uncertainty to compute the χ 2 of Eq. (2.44). By following the same workflow discussed above for the single oscillator, we reconstruct the spectral densities for different values of Θ and number of delta functions in the model, N ω . In Fig. 5.7, we show the influence of the discretization N ω by fixing the canonical value Θ = 1. Following the same procedure as above, we calculate again χ 2 val for the validation set by generating test correlation function from the exact result of Eq. (5.4), with the same variance σ k . The values of χ 2 val , shown in Fig. 5.9, indicate again a more statistically sound reconstruction corresponding to sparse grids. Unfortunately, none of the curves of Fig. 5.7, convincingly captures the sharp edges of the exact spectral density, which rather resemble two symmetrically broadened peaks. Considering shifted grids (Fig. 5.8), however, as also quantitatively supported by the validation procedure, results in contrast in more asymmetric features, clearly improving the reconstruction towards the exact spectrum. Note, however, that employing sparse ω-grids considerably limits frequency resolution, so that the reconstruction in the case of the continuous spectrum with its sharp discontinuities remains quite difficult.

Extension of the reconstruction approach

We tested the stochastic reconstruction method coupled with the validation step on several simple models. We note that the approach can be developed further and tested on the spectra of real solids, for example, by including spectra with van Hove singularities. It is also worth testing the reconstruction approach on other more complex models such as the models with non-harmonic Hamiltonians. The prominent example of such models is the model with Fermi resonances or the system of two coupled harmonic oscillators (see, [START_REF] Plé | Anharmonic spectral features via trajectory-based quantum dynamics: A perturbative analysis of the interplay between dynamics and sampling[END_REF]. One should study carefully the dependence of the reconstruction results for various parameters and adapt the discretization of spectral density. We remark that the validation procedure introduces an unambiguous criteria for selecting the most optimal spectrum. Here we compare the results pertaining to a grid shifted by δω = 0.25 to those with δω = 0, the usual (not shifted) case. The shaded area indicates the exact spectral function.

0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 S( ) N =10 =0, =4 =0.25, =10 N =25 =0, =15
N = 10 N = 25 N = N = N = 10 N = 25 N = N = FIGURE 5
.9: Main panel: χ 2 val from the validation procedure of spectral function at the corresponding values Θ and N ω for C cont pV (τ). The area of the circles is proportional to the value of χ 2 val . Purple circles correspond to the results for an ω-grid shifted by δω = ∆ω/2. Inset: χ 2 val as a function of the effective temperature Θ, at the indicated values of N ω .

6 Phonon lifetimes in rare gas solids, implication for heat conductivity

After having introduced and tested our methods on simple models with known, analytical solutions, we will apply them in the present chapter to the Lennard-Jones solid towards first, preliminary results on the heat conductivity. As we have seen in the second chapter, linear response theory relates the heat conductivity to the heat (or energy) current-current correlation function. Its behavior in imaginary time can be directly obtained from PIMC, and the corresponding spectrum from the stochastic reconstruction. Although this program is rather straigthforward, we will not pursue it in the following. Since we have leave the playground of exactly solved problems, comparison with approximate theories or experiment are the only possibility to judge the quality and validity of our results. However, the finite size of our simulation cell used in PIMC usually prevents direct comparisons and results must be properly extrapolated to the thermodynamic limit. Therefore, we have started to focus first on more "elementary" quantities, like phonon frequencies and life times. These quantities are rather well understood, and provide the basic quantities underlying many thermal and transport properties of the solid. Still, a quantitative description of those quantities, in particular of lifetimes of the quantum crystal obtained from PIMC is lacking in literature. Here we provide a first study, comparing our results with those of obtained from classical molecular dynamics. We use our results of the phonon frequencies and lifetimes to estimate the heat conductivity directly via elementary transport theory without invoking Green-Kubo The question of heat transport is a long standing one and it has been addressed for various systems. The heat transport in solids can be understood in terms of several separate contributions. One of the component comes from the motion of free electrons inside the bulk and their interactions with the lattice structure. In the case of metal systems these effect plays the crucial role where the electrons define the properties of material. Another contribution to the heat transport comes from the vibrations of the underlying lattice itself.

For metals it is smaller than the electron part and can be safely neglected for most calculations. However, starting with the semiconductors systems and continuing with insulating solids composed of rare gas atoms the heat transport is conducted mainly by means of lattice vibrations. The pure crystal vibrations is an idealized concept as seen from the fact that for a system with quadratic potential normal modes cannot interact with one another and hence there is no thermal equilibrium [START_REF] Debye | Vorträge über die kinetische Theorie der Materie und Elektrizität, Gottinger Wolfskehlvortrage[END_REF]. In this system the heat would travel with the speed of sound yielding infinite (lattice) conductivity. On the contrary, in the real system the anharmonicity in the potential energy and the defects of the lattice force modes to interact and decay which results in a finite heat conductivity. We will choose as a model that captures the most prominent features of interparticle interactions inside the insulator the Lennard-Jones crystal [Jones, 1924] (Fig.(6.1)). It is a common practice for calculations to consider the motion of lattice nodes as classical one. As we shall see in this chapter for the case of insulating systems it is not always justified, because the quantum effects can modify drastically the parameters of these vibrations thus affecting the predictions for the heat transport coefficients.

The outline of this chapter is following. We shall provide a brief description of the model, highlighting several prominent issues that will be important in the following analysis. Then, we shift our focus to the classical theory of harmonic solids, putting more attention to the calculation of mode correlation function and its spectral density. We then generalize this approach to systems with pair-wise anharmonic interactions. We will briefly explain the role of anharmonicities on the quantities of interest and rederive the expression for the heat conductivity using this method, which is our ultimate objective for the analysis. We then proceed with the actual calculation of the mode correlations, first in the classical approximation and then in the quantum regime. We shall mention the techniques that we use in both of these cases and provide the results of the calculation. We will concentrate on the spectral function which we will study individually for each mode, particularly focusing on a peak frequency of the spectrum and the broadening of the peak which indicates the phonon lifetime. We then will discuss the role of quantum effects for these quantities and demonstrate the different predictions for the heat conductivity arising from classical and quantum computations.

Lennard-Jones model of rare gas

The expansion of the interaction energy

Powered by TCPDF (www.tcpdf.org) FIGURE 6.1: Comparison [START_REF] Maitland | The intermolecular pair potential of argon[END_REF] of experimentally measured interaction between atoms of argon (solid line) and LJ potential(dashed line) (reprint from [START_REF] Allen | Computer simulation of liquids[END_REF])

As we outlined above we are interested in studying the purely insulating material, represented by the rare gases.

The distinctive feature of the noble gases is the filled electron shells which leads to the absence of free electrons inside the crystal bulk. Thus, we can model the system by ensemble of simple atomic particles placed on a lattice nodes and interacting via some potential. In a real system the constituents interact via various different processes. In general the potential energy can be written as a sum of interactions between various number of particles -single atoms, pairs, triplets, quadruples etc. The potential energy in this case depends on the positions of atoms {r n } and have a form:

V = N ∑ i=1 v 1 (r i ) + N ∑ i,j>i v 2 (r i , r j ) + N ∑ i,j>i,k>j>i v 3 (r i , r j , r k ) + . . . , (6.1
) where N denotes the total number of particles in the system. The sums over pairs, triplets ... are taken in such way in order to avoid counting the interaction several times. The first term in the expression (6.1) encompasses the effects that the external field has on the system. Among the interaction terms the leading contribution is provided by the pair-wise potential. As a side remark, we note that from the symmetry one can see that the pair-wise interaction can only depend on the distance between particles, that is v 2 (r i , r j ) ≡ v(r ij ), where we introduced the notation r ij ≡ |r ir j |.

6.1. Lennard-Jones model of rare gas 77 For the purpose of numerical analysis the terms in the expansion (6.1) that involve more than two particles are generally excluded from the analysis due to the increasing difficulty of the simulation.

The interaction between two particles of the rare gas is well studied experimentally [START_REF] Tanaka | Absorption Spectrum of the Argon Molecule in the Vacuum-uv Region[END_REF][START_REF] Buck | Determination of the Interatomic Potential of Krypton[END_REF][START_REF] Aziz | The Ne-Ne interatomic potential revisited[END_REF]]. An estimation of this potential is shown on the Fig.( 6.1). One can see several distinct properties of the potential. First, one observes the quickly rising repulsion at short separation due to the overlap of the particle electron clouds which do not form a bond configuration. Then, there is the weakly attractive tail at sufficiently large distances between particles which arises essentially due to the correlation between the electron clouds. There is also a potential well which allows to form a bound state.

In numerical simulation one typically employs a rather simple and idealized potential that approximates the empirical one. A common choice is the Lennard-Jones model with appropriately defined parameters:

V = 4ε ∑ i<j σ r ij 12 - σ r ij 6 (6.2)
Here the value of σ sets the coordinate scale, whilst sets the energy scale of the model. In solid phase the atoms form a lattice whose structure is defined exactly by the interactions between atoms and certain external parameters. At sufficiently low temperatures and pressures the particles position themselves in a face-centered cubic (fcc) lattice. Certain freedom, however, remains in defining the exact value of crystal density and the related lattice spacing. We shall treat the crystal density as an external parameter. This model captures the prominent features of the interaction inside the material, including the energy well and the repulsion barrier.

Conventions

For the purpose of numerical analysis, we introduce several conventions that are commonly used in the field. The Lennard-Jones potential approaches zero value quite fast with increasing particle separation, so that one can neglect the interaction between particles far away. Furthermore, the particle-particle interaction becomes screened by the other particles in the system which means that the Lennard-Jones potential is in fact an effective potential. In our simulations we select the value of cutoff of the potential at R cut = 2.5σ which allows to account for the effect of nearest and next nearest neighbours on a particle. That is, we replace the (at least in principle) infinite range interaction with

v(r ij ) =    4ε σ r ij 12 -σ r ij 6 , if r ij ≤ R cut 0, if r ij > R cut , (6.3)
In order to make this potential continuous at the point r ij = R cut we will shift all the values of interaction energy upwards by some constant, so that v(R cut ) ≡ 0.

Clearly, the value of the shift is defined by the cutoff separation:

e cut ≡ 4ε σ R cut 12 -σ R cut 6 .
The system that one can study numerically is necessarily of the finite size and particle number. In view of minimizing the effect of boundaries we employ the periodic boundary conditions and simulate the ensemble of particles on a torus: r i = r i + L with L being the linear box size. It may however occur that, even within 78 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity the cutoff range, the interaction between certain particles is taken more than once in the calculation of Eq.(6.2) due to the various possible images of the particles. To avoid the overcounting we introduce another important simplification -the nearest image convention (Fig.( 6.2)). The idea is as simple as the name suggests -one chooses R cut ≤ L/2, such that one needs to account each interaction only once keeping track of the separations between various copies of the system, only those which correspond to the smallest distance r ij modulo L.

PIMC of the Lennard-Jones solid

FIGURE 6.2: In nearest image convention one considers several copies of the original system and accounts only the interaction that corresponds to the smallest separation between particles. In this example only the interaction between green particle and blue one from the left-upper image of the system is accounted.

In the following, we consider N particles interacting via a Lennard-Jones potential. In the path-integral representation, our configuration of particles will have a classical Boltzmann weight ∝ exp[-S] where the action S is given by

S = M ∑ n=1 N ∑ i=1 mMσ 2 2h 2 β r n+1 i -r n i σ 2 + + 4ε β M ∑ j>i   σ |r n i -r n j | 12 - σ |r n i -r n j | 6   , (6.4)
The indices i and j run over the total number of particles whereas the index m goes over the M time slices. Looking at the expression 6.4 we can now introduce several dimensionless variables.

Working with the Lennard-Jones potential, v(r) = 4ε[(r/σ) -12 -(r/σ) -6 ], it is natural to measure length in units of σ and energies in units of ε. Namely, one can consider quantity V(r)/ε as a dimensionless function of the dimensionless argument x = r/σ. Next, we define a characteristic time t 0 for our system by setting the kinetic energy ≈ mσ 2 /t 2 0 equal to ε. We get t 2 0 = mσ 2 /ε which can be used as unit of time, e.g. ωt 0 = ω √ mσ 2 /ε will be dimensionless for frequencies ω.

For quantum system, we have additional dimensionless combination Q = h2 /εmσ 2 which allows to rewrite characteristic time as t 0 = h/εQ, so that we have hω/ε = ωt 0 h/(εt 0 ) = ωt 0 h2 /εmσ 2 = Qωt 0 . Therefore, when referring to "Lennard-Jones units", one writes all formulas in terms of T * = k B T/ε, ω * = ωt 0 , x = r/σ and hω/k B T = Qω * /T We denote the combination β as β * and the product h2 /( mσ 2 ) as Q 2 . Then the formula 6.4 will take the form

S = M ∑ m=1 N ∑ i=1 M 2Q 2 β * x m+1 i -x m i 2 + + 4β * M ∑ j>i   1 |x m i -x m j | 12 - 1 |x m i -x m j | 6   . (6.5)
where we have also introduced the dimensionless variable for coordinates x = r/σ.

The ratio Q (which essentially characterizes how "quantum" the system is) the effective inverse temperature β * , along with the length of the box L/σ and the cutoff distance R cut /σ for the Lennard-Jones potential constitute the total set of parameters that one needs for the simulation. The size of the system is related to the density of the crystal ρ which we select instead as the parameter of the simulation. A reasonable choice for the density is dictated by the consideration that one should have approximately 1 particle per unit of space, that is, the density needs to satisfy the relation ρσ 3 ≈ 1. In the following we fix the crystal density for all the calculation at value close to 1 in the parameters of our simulation (more exactly, we choose density ρσ 3 = 0.965). There is more freedom in choosing values of β * , we select such temperatures for which the crystal is well below the melting point. In order to choose Q we have following guidance. One expects that the lighter a particle the more pronounced NQE, which means that among the rare gases neon suits the best for the analysis. The choice of neon translates to the value Q ≈ 0.092.

Harmonic analysis

Expansion of the potential

It is not possible to obtain exact analytical results for a many-particle system interacting with a Lennard-Jones potential, so one is forced to resort to approximations.

A particularly useful approach is the harmonic approximation [START_REF] Ashcroft | Solid State Physics[END_REF]. In this part we shall briefly overview the main concepts and results for Lennard-Jones system in the harmonic approximation. Let us define the set of coordinates for particles in the crystal (at zero temperature, i.e. when particles form the fcc lattice) as {X 0 }. For clarity of notations we omit the time indexes. One can then expand the potential up to the second order as

V = V 0 + N ∑ i (x i -x 0 i )(∇ i V X 0 ) + N ∑ ij (x i -x 0 i )(∇V i ∇V j X 0 )(x j -x 0 j ) (6.6)
Here the V 0 is the minimum value of energy and does not play a role in the following analysis so we will consider the energies as a shift from the minimum. Each term in the first sum is nothing else but the force exerted on particle in the position x 0 i at the equilibrium. Since at the equilibrium there is no total force acting on any particle, this sum must be equal to 0. The second sum can be rewritten as

V harm = 1 2 ∑ ij (x iα -x 0 iα )D αβ ij (x jβ -x 0 jβ ), (6.7)
where we put the labels α, β . . . for the vector components. The matrix

D ij ≡ ∂ 2 V ∂x i ∂x j X 0
is the hessian matrix with values calculated at the lattice coordinates {X 0 }:

     D αβ ij = V ij (x 0 ij ) 2 - V ij (x 0 ij ) 3 x 0 ij,α x 0 ij,β - V ij x 0 ij δ αβ , for i = j D αβ ii = ∑ j =i D αβ ij , otherwise (6.8)
A useful approach to describe excitations in the system is to represent them as quasiparticles. These excitation can be characterized by the energy hω and the quasimomentum hk defined for each particle inside the Brillouin zone. The simplest way to calculate the frequencies is to diagonalize the matrix D and to obtain the eigenvalues. In the following, we will use translational symmetry to blockdiagonalize the matrix, such that ω(k) can be obtained much more efficiently.

Normal modes

We remind ourselves that the primitive vectors of the fcc lattice can be written

       a 1 = {0, a 2 , a 2 } a 2 = { a 2 , 0, a 2 } a 3 = { a 2 , a 2 , 0}
, where a is the separation between nodes of the lattice. These vectors determine the reciprocal lattice (bcc) with the corresponding basis which satisfies the relation

a i b j = 2πδ ij :        b 1 = {-2π a , 2π a , 2π a } b 2 = { 2π a , -2π a , 2π a } b 3 = { 2π a , 2π a , -2π a }
The idea of employing normal modes for the analysis of crystal was introduced in the works [START_REF] Born | Theoretical investigations on the relation between crystal dynamics and x-ray scattering[END_REF] and [START_REF] Begbie | Thermal scattering of X-rays by crystals -I. Dynamical foundation[END_REF]. We have a system of 3N equations describing the behaviour of each of the particle component:

m üiγ = - ∂V harm ∂u iγ = -∑ j D γβ ij u jβ (6.9)
where we put for brevity u iα ≡ x iαx 0 iα . One can recast these equations in a vector form combining all components of a particle together

m üi = - ∂V harm ∂u i = -∑ j D ij u j (6.10)
We are looking for a regular solution in the form of plane waves:

u l (t) = κe -i(kR 0 l -ωt) (6.11)
Vector κ is the polarization of the normal mode which describes the direction of wave propagation. Periodic boundaries of the crystal translate into the Born-von Karman condition on the solution u: u(R 0 k + N γ a γ ) = u(R 0 k ) for each of the 3 primitive vectors a γ . Here N γ are quite large integers taken such that N 1 N 2 N 3 = N. The choice of boundary condition puts 6.2. Harmonic analysis 81 a restriction on the allowed values of wave vectors k: (6.12) with n c being the integrals whilst b c denote the basis of reciprocal lattice. The translational invariance of the crystal imposes periodicity on the wavevector which means that only values of k corresponding to the primitive cell of the reciprocal lattice will produce a distinctive result. Hence there exists only N solutions It is typically convenient to assign this cell as a first Brillouin zone. We can now substitute the solution (6.11) into the equation of motion: (6.13) with D(k) = ∑ j D ij e -ikR 0 j , which is independent of value of i due to the translational symmetry. The matrix D(k) is typically referred to as a dynamical matrix.

k = n 1 N 1 b 1 + n 2 N 2 b 2 + n 3 N 3 b 3 ,
mω 2 κ = D(k)κ,
Let us look a little closer at the dynamical matrix. The matrix D αβ ij has several properties of interest. First, it is symmetric under the exchange of pair of indexes:

D αβ ij = D αβ ji and D αβ ij = D βα ij .
Second, its sum over all particles is equal to zero: ∑ i D αβ ij = 0 which is simply the embodiment of the fact that if all the particles in the crystal are shifted by the same u α i ≡ u then the energy of the system remains constant. With this observations one write the matrix D(k) as

D(k) = ∑ j D ij e -ikR 0 j = 1 2 ∑ j D ij e -ikR 0 j + e ikR 0 j -2 = = ∑ j D ij [cos(kR 0 j ) -1] = ∑ j D ij sin 2 (kR 0 j /2) (6.14)
We see that the sum is strictly positive and symmetric. From the linear algebra one knows that the symmetric positive-defined matrices has the same number of eigenvalues as the rank of the matrix. They satisfy the equation:

D(k)κ s = λ s (k)κ s (6.15)
The relation with normal mode frequency is then straightforward:

ω γ (k) = λ γ (k) m (6.16)
Thus, putting in different values of k from (6.12) one recovers all the normal frequencies of the Lennard-Jones potential.

Effect of higher order terms in the expansion of potential

In the previous analysis we considered only quadratic terms in the expansion of Lennard-Jones potential. For many practical purposes this is not sufficient as we shall see in the example of heat conductivity in the next part. The full analysis is a challenging task and we refrain ourselves to a qualitative description. The formalism in terms of normal modes is still applicable, yet there will be some noticeable changes.

In the following analysis we will compare two distinct setups: a classical system governed by the Hamiltonian H = ∑ i p 2 i 2m + ∑ ij v LJ (r ij ) and its quantum equivalent 82 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity with the usual replacement by the operators. In the classical case the behavior of the ω is dictated predominantly by the thermal fluctuations from the equilibrium (lattice) configuration. Thus, at the limit of zero temperature the fluctuations are negligible and the mode frequency approaches the one of the harmonic crystal. With higher temperature the system is able to probe configurations further away from the minimum and the difference between full potential and the quadratic approximation becomes more significant leading to a larger shift in frequencies. Since the thermal fluctuations grow linearly with T one can expect that the normal modes will follow the same pattern.

The picture is different for the "quantum" scenario where already at zero temperature the system has some nonvanishing ground energy. In this case the behaviour is more difficult to analyse and we will discuss it in more details later in this chapter. One can anticipate however that the effect of temperature on a normal mode will be less pronounced than in the classical case and for high enough ω even irrelevant due to the fact that in this limit the temperature acts mostly on the density of state operator 1 e βhω -1 which is small in the same limit. Another important consequence of treating full potential is that there will be non-zero interaction of any given mode with large number of other modes. This will show itself in the broadening of the peaks δ-peaks of the spectrum. This behaviour is similar in both classical and quantum systems, however in the classical regime one expects that the broadening vanishes with T → 0 and grows with T increasing. In the quantum case the zero-point motion should play a role apparently leading to a non width of the peak at the same limit of T → 0.

Heat conductivity

The scattering of phonons and their finite lifetime give rise to thermal resistivity which manifests itself in the (finite) heat conductivity. For the analysis of lattice heat transport one typically employs one of the general approaches: Green-Kubo approximation, variational method and the relaxation time approach. Here we show derivation using the latter method. The relaxation time approach (as well as variational method) is based on the Boltzmann transport equation [START_REF] Peierls | Zur Kinetischen Theorieder Wärmeleitung in Kristallen[END_REF] and it was successfully used to describe heat transfer in semiconductors [START_REF] Broido | Intrinsic lattice thermal conductivity of semiconductors from first principles[END_REF][START_REF] Zhou | First-principles calculations of thermal, electrical, and thermoelectric transport properties of semiconductors[END_REF]] and whilst the Green-Kubo formalism tackles the problem from a more fundamental angle using the correlations of quantum operators.

Kinetic theory approach

In order to apply the Boltzmann equation we require the existence of the distribution n k (r, t) which shows the mode occupation number in the vicinity of r at time t. We assume that in the system there exists a temperature gradient across the sample. Then, one can study two main processes that cause the change in the phonon occupation number -diffusion and scattering (decay) into other phonons.

One can describe the effect of diffusion of phonons with the use of temperature gradient:

∂n k ∂t diff = -v k ∇T ∂n k ∂T (6.17)
where v k is the phonon group velocity.

Harmonic analysis 83

Another process that contributes to evolution of n k (r, t) is a phonon scattering. The rate of variation due to decay is given by the relaxation time τ k as (6.18) which means that the rate is proportional to the deviation of n k from the equilibrium distribution. Since in the steady state of heat flow through the crystal the total time derivative of the occupation number vanishes, one can write:

∂n k ∂t scat = - n k -nk τ k ,
∂n k ∂t scat + ∂n k ∂t diff = 0
Substituting into this expression formulas from above one get:

-

n k -nk τ k = v k ∇T ∂n k ∂T (6.19)
In the steady state with a temperature gradient ∇T the heat flow per unit area in the direction normal to the said gradient is given by the macroscopic expression

J = -λ∇T (6.20)
where λ denotes the heat conductivity. For dielectric materials one can express the microscopic heat current as the summation of energy per phonon mode over all possible modes of the crystal:

J = 1 V ∑ k hω k (n k -nk )v k (6.21)
We notice that only modes that contribute to the heat flux are those that have the nonequlibrium phonon occupation number. Using the Eq.( 6.19) we bring this expression to the form of Eq.(6.20) with

λ = 1 V ∑ k v 2 k τ k ∂ k ∂T = 1 V ∑ k l k ∂ k ∂T ∂ω k ∂k (6.22)
Here ∂ k ∂T is the modal heat capacity, ∂ω k ∂k is the definition of the mode group velocity and l k is mean free path for phonon. The mean free path is the characteristic distance that phonon travels between different scattering events. In the absence of anharmonic interactions normal modes do not interact with one another and they propagate through the whole crystal unperturbed thus producing an infinite mean free path l k which in turn results in an infinite heat conductivity. This is clearly unphysical which means that any real system ought to have some sort of processes that forces modes to decay.

Let us rewrite the formula (6.22) in units of Lennard-Jones as discussed in Section (6.1.3):

λ = 1 V ∑ k v 2 k τ k ∂ ¯ k ∂T = 1 σ 3 (V/σ 3 ) ∑ k (σ/t 0 ) 2 (v k t 0 /σ) 2 t 0 (τ k /t 0 )k B ∂(h ωk ) ∂(k B T) = = k B σt 0 V * ∑ k v * 2 k τ * k ∂ ¯ * k ∂T * (6.23)
84 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity where ¯ ≡ hω k [(exp[hω k /k B T] -1) -1 + 1/2] is the mean energy per mode. From the last line we see that quantity λσt 0 /k B is dimensionless.

Green-Kubo linear response

In the first chapter we showed the general formulation of the linear response theory.

In that treatment we used a system at the thermal equilibrium being subject to the external (mechanical) perturbation. Despite the fact that the heat transport cannot be associated with any mechanical perturbations, one can still formally derive similar expression for heat current as we did before [START_REF] Kadanoff | Hydrodynamic equations and correlation functions[END_REF]. The approach is based on the assumption that the system is locally at thermal equilibrium, i.e. the temperature changes sufficiently slow on the length and time scale in comparison with atomic processes. For a smooth function T(r) that can be written as a small perturbation from the constant temperature T(r) = T + ∆T(r) with |∆T| T. Then the effect of the perturbation can be described by This potential resembles the general mechanical perturbation that we considered earlier and thus we can adopt the same attitude and write for the system whose only conserved quantity is energy, the heat current via the the temperature gradient: J = -λ∇T.

P ∝ e
(6.26)

The heat conductivity that enters this equation can be expressed from the Green-Kubo relation as the current-current correlation function:

λ αβ = V k B T 2 ∞ 0 J α (t)J β (0) 0 dt (6.27)
The macroscopic current J has the form (6.28) where we used the continuity equation in going to the second expression. Expanding this formula one can write

J(Γ) = 1 V j(r)dr = 1 V ė(r)rdr,
J = 1 V ė(r)rdr = 1 V ∑ i ∂e ∂R i dR i dt + ∂e ∂P i dP i dt rdr (6.29)
The formula for the heat current is not uniquely defined -one can, for example, consider antisymmetric terms in the particle interaction which does not change the total Hamiltonian of the system whilst modifying the current-current correlations. These various expression for the heat current, however, result in a singly defined 6.3. Analysis of equilibrium properties 85 heat capacity [START_REF] Ercole | Gauge invariance of thermal transport coefficients[END_REF] which stresses the fact that the heat capacity is a proper physical observable unlike the heat flux.

Analysis of equilibrium properties

Before moving to the calculation of the dynamical properties, let us first analyze the time independent thermodynamic quantities and compare the results with already existing studies, numerical and experimental. In this section we look into the behavior of kinetic energy and pressure for a range of parameters of the simulation.

Kinetic energy

We have seen in the Chapter 3 that calculation of kinetic energy can pose certain difficulty, thus it serves as a good benchmark for our approach. We compare the values of kinetic energy with reference data from the work [START_REF] Cuccoli | Monte Carlo computations of the quantum kinetic energy of rare-gas solids[END_REF]. For our analysis we use the "virial" estimator for the kinetic energy that we introduced in Chapter 3. In this section we simulate the crystal of N = 108 particles. We start by studying behavior of kinetic energy at a fixed temperature (T = 5K) and fixed density of crystal (ρσ 3 = 0.965) for various time discretization M. Fig.( 6.3) shows the results of the calculation. To make comparison with reference analysis we use the same separations M = 4, 8, 16, 24, and 32 as indicated in the paper. The result is then plotted as a function of 1/M. The data for sufficiently large M usually can be approximated by a smooth curve (in the current example -straight line) which tends to a well defined limit as 1/M → 0. our simulation reference data FIGURE 6.3: Kinetic energy per particle e K of neon for several imaginary time discretization numbers M. The system is studied at T = 5K and ρσ 3 = 0.965. We compare the results of our simulation (red) with the data obtained from QMC simulation of analogous system in the work [START_REF] Cuccoli | Monte Carlo computations of the quantum kinetic energy of rare-gas solids[END_REF].

On Fig.( 6.4) we also compare the calculation of kinetic energy at different temperatures. Again, we follow the work [START_REF] Cuccoli | Monte Carlo computations of the quantum kinetic energy of rare-gas solids[END_REF] and compare the quantity at the indicated values of parameters, namely we adopted the values of crystal density following the Table 6.2. In both these examples we observe good correspondence between current simulation and the previously performed analysis.

We remark that in the reference study the crystal density (as indicated in the Table 6.2) was adapted with temperature such that the system was at zero pressure at all values of T. On the contrary, in our simulation we observed non-zero pressure at indicated values of parameters. The discrepancy is caused by the fact that the pressure calculations, unlike computations of kinetic energy, are quite sensitive to the form of the potential energy and in particular to the smoothness of the interaction at the cut-off distance. Since in our work we did not consider smoothed potential it is natural to expect the deviations from reference data, for example obtained here [START_REF] Müser | Path-integral Monte Carlo study of crystalline Lennard-Jones systems[END_REF]. Nevertheless, as we shall see in the following test, even with such shape of the potential we are able to obtain results for the thermal expansion of crystal lattice, which are reasonably close to the experimentally measured values. 6.2. The red data indicate the present simulation, whilst the blue ones -reference data from paper [START_REF] Cuccoli | Monte Carlo computations of the quantum kinetic energy of rare-gas solids[END_REF]. The reference data were approximated in the limit ∆τ → 0.

Thermal expansion of lattice constant

Before moving to the study of phonon time correlation function let us look in more details into the equilibrium properties of the system. An interesting quantity to compare classical and quantum calculations is the thermal expansion of lattice. In the experiment the quantity is typically measured at constant pressure with varying temperature. The QMC algorithm that we use allows us to obtain properties at constant volume, and in order to study the properties at constant pressure we need to simulate the system at different densities. On ρ for N = 108 particles. Increasing the temperature from 5K to 20K we observe the lattice parameter changes by approximately +1.2% which is reasonably close to the experimentally observed value of +0.9% [START_REF] Batchelder | Measurements of Lattice Constant, Thermal Expansion, and Isothermal Compressibility of Neon Single Crystals[END_REF] for the same temperature change. We note that the classical calculation of lattice parameter produce higher shift of around 7%. 6.4. Molecular Dynamics simulation of mode correlation function 89

Molecular Dynamics simulation of mode correlation function

As we mentioned above, we perform the analysis of classical and quantum systems of Lennard-Jones crystal. In order to separate the "pure" quantum effects from those that can be understood without evoking quantum mechanics, in this section we perform a calculation of a crystal that is governed by the classical Hamiltonian.

We shall focus primarily on the temperature dependence for quantities of interest, i.e. we will look at the spectral function for each mode and study the frequency of the peak and the broadening of the peak caused by the interaction of a normal mode with other modes. We take the temperatures well below the melting point and study the range T = [5K, 20K], which translates to the range [0.136, 0.545] in the dimensionless LJ units. We will also consider systems consisting of different number of particles (N = 108 and N = 864) in view of estimating size effects. In the simulation we will focus on neon due to the fact that among the rare gases neon possess the lightest mass. Keeping in mind the Hamiltonian of the system, one can expect that the Nuclear Quantum effects will be the strongest in comparison to other rare gases. Before discussing of results let us say a few words about the numerical procedures and algorithms that we used in order to obtain the desired correlation functions and to study the respective spectral densities. Without going into much details we shall focus on several questions of the simulations, both MD and QMC, and make some notes considering the error estimation in both of the cases. In this part we treat Molecular Dynamics calculation and then in the next part we move to the Quantum Monte-Carlo. Due to the fact that for QMC methods it is notoriously difficult to treat real time correlations and that for a classical system there exists many reliable well-developed tools, the machinery that we introduced in the first chapter of this thesis is not particularly convenient to use here for the simulation. In the present study we will employ Molecular Dynamics (MD) for the calculation of modes time correlation functions. We do not develop any particularly method, instead we opt for the use of LAMMPS algorithm [START_REF] Thompson | LAMMPS -a flexible simulation tool for particlebased materials modeling at the atomic, meso, and continuum scales[END_REF] for this task.

At its base the MD approaches are build around the idea of solving numerically equations of motion. In practice it is typically done by calculating the forces acting on particles at each point of discretized time and then using these forces to displace the particles according to the second Newton's law. By adjusting a time step suitably one can access a wide range of time dependent quantities on different time scales. In our case the time step is determined by the frequency of oscillations t * ∝ 1/ω which we estimate from the quadratic expansion of the potential. 90 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity The mode correlation function that we are interested in is calculated through the position correlation function which is simply the product

C xx iα (t) = x iα (t 1 + t)x iα (t 1 ) , (6.30)
where t 1 is some arbitrary time. For the calculation of modes we use diagonalization of Hessian matrix D in the form (6.31) where the matrix Λ is diagonal and the O is the orthonormal matrix composed of eigenvectors of D. In this case the modes can then be computed and treated separately via application of the transformation O as A p = ∑ j O lj x j which leads to a correlation function of the form:

D = OΛO -1 ,
C AA pγ (t) = A pγ (t 1 + t)A pγ (t 1 ) (6.32)
In the simulation the average (6.32) is obtained by summing over different values of t 1 in such fashion that preserves t constant. On Fig.(6.7) we show the sample of correlation function for 12 degenerate modes that correspond to the second lowest frequency. In black is plotted the average of all these degeneracies. One clearly observes the oscilations with certain frequency modulated by the decaying rate. The longer time oscilations in the correlation function that are visible on the figure is simply the data noise that can be reduced by taking more values of t 1 , and one can see the reduction of these oscilations already from the averaging over degeneracies. In the end, the correlations at long times are not relevant for the present analysis as the most information comes from the small times. This figure also shows that whilst shown mode correlations should have precisely the same frequencies along with decay rate (at least in theory), in practice the result is affected by some noise. The main source of the noise is the stochastic fluctuations which can be reduced in the same way we just mentioned -by averaging over more data. This, however, would require increasingly more data storage, which means that the compromise is necessary.

Thus, one cannot avoid having uncertainty on the simulation output. However, the proper error analysis is tricky since the error distribution does not obey the Gaussian law so we will not go into this elaborated calculation. Instead, we estimate the errorbars from the following consideration. We will not study the results for each component separately, but we will group them by degeneracy and analyze the configuration averaged over such block. The results for each individual mode of the block will be then plotted as the uncertainty. Unlike the simulation of the quantum system where one is forced to work with imaginary times, the correlations for the classical system can be readily calculated at real times. This makes the correlation function easily accessible for the spectral analysis which is a straightforward Fourier transform. One can then approximate the spectrum with a Gaussian or Lorentzian function in order to evaluate the frequency of the peak and its width. The example of a spectrum and the corresponding fits are shown on Fig. (6.8). Making a fit for a given mode possess an uncertainty in defining the position of the peak and, mainly, in its width and amplitude (Fig.(6.8)). This uncertainty, however, is smaller than the uncertainty across the degenerate modes and thus already accounted by the error estimation.

As we expected earlier there is a broadening of the peak due to the interaction of a mode with other modes which is in a stark contrast with a single δ-peak of a 92 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity harmonic potential. Whilst the position of the peak is well defined and corresponds to the maximum overlap of the eigenvectors (which is clearly κ iα κ iα interaction with itself), the spectral density does not decay exponentially fast and one observes heavy tails indicating that there is a nonzero overlap in eigenvectors even for modes far apart in frequency space.

Quantum Monte-Carlo calculation

After discussing the simulation of a classical system we now move on to the quantum case. As a main tool for the analysis of the system governed by the quantum Hamiltonian in this section we will again use the QMC techniques developed in the first chapter. This section is dedicated to a brief description of the simulation with emphasis on several prominent issues. In the next section we will collect all the results from this and the previous section and discuss them together.

For the purpose of comparison with the classical results we perform the simulations on the systems described by the same parameters introduced in the previous section. Similar to the above studied problem of a single oscillator the output of the MC simulation is affected by different biases, systematic and stochastic. One deals with the stochastic uncertainty in the usual way, i.e. by controlling the length of the simulation. Since the observable does not include momentum dependence there is no reason that one can reduce errors by rewriting the MC estimators. The systematic deviations come from the discretization of time in the Trotter decomposition and one can reliably control these discrepancy as we showed in the Chapter 3. For the analysis there we used the exact expression for the correlation function of interest and tested it against the numerical calculation of the path-integral formula. One can do similar analysis for the purely harmonic approximation of Lennard-Jones interaction. For the case of a single oscillator the position correlation function have the form:

x(τ)x(0) = h 2mω cosh(hω(τ -β/2)) sinh(hωβ/2)
The mode correlation function in the harmonic approximation have exactly the same form. Translating the expression into Lennard-Jones units one writes

(C harm ) AA iα = A iα (τ)A iα (0) harm = Qσ 2 2ω * iα cosh(Q ωiα (τ * -β * /2)) sinh(Q ωiα β * /2) , (6.33)
we remind that Q = h2 m σ 2 and β * = β . In writing this formula we used the expression for frequency in LJ units which comes from the dimensional analysis

ω * = mσ 2 ω (6.34)
The analysis presented in the Chapter 5 for a single oscillator let us estimate that the systematic deviations are small when time discretization satisfies the relation βhω M ≈ 0.1. In the present simulation one can use for the estimate the highest frequency (which in our case is approximately 20 in the Lennard-Jones units) and write this relation as

β * Q ωmax M ≈ 0.1 (6.35)
which sets the number of time slices M ≈ 200β * Q.

Quantum Monte-Carlo calculation 93

In order to make use of the formula for the spectral reconstruction one also needs first to translate the expression (4.18) into the proper units. Since we are not interested in spectral function unless it can be compared with Molecular dynamics result, we will write it immediately in the LJ units. We also denote explicitly that the spectral density depends on the broadening of the peak:

C(τ) = dωS( ω, ω0 , γ) e -τhω + e -(β-τ)hω = = mσ 2 d ωS( ω, ω0 , γ) e -τ * Q ω + e -(β * -τ * )Q ω (6.36)
For the analysis of the spectral density we can again use the method that we developed in Chapter 3. However, the study of the harmonic oscillator in the Chapter 4 showed us that the fitting of spectral function with the discrete peaks (with a trial spectral density S(ω) = ∑ j A(ω j )δ(ωω j ))is not particularly useful for the description of continuous spectra. In order to get an accurate fit it is preferable to introduce several auxiliary parameters and then study the behaviour for different parameter values. Since we are interested in the peak frequency and the width of the distribution, in the end we need to approximate the resulting spectra with a Gaussian (or Lorentzian) in analogy to what we did for MD calculation in the previous part. Instead of doing this, in this section we search for a spectral density in the form

S( ω, ω0 , γ) = A exp - ( ω -ω0 ) 2 2 γ2 , (6.37)
which depends only on few parameters. We then make a simple minimization procedure on a functional C AA iα (τ) -CAA iα (τ) with mode correlations CAA iα (τ) = dωK(ω)S( ω, ω0 , γ) being constructed from the trial spectrum.

Let us compare this minimization to the stochastic reconstruction we used before. We remind ourselves that following the recipe from the Chapter 3 we need to discretize the frequency range with N ω delta-functions each contributing some weight A i ≡ A(ω i ) to the spectral function. The resulting frequency discretization ∆ω along with several others parameters constitute the parameter space for the model. We do not provide here an exhaustive analysis and restrict our attention to the discussion of the reconstruction for a given ∆ω with varying effective inverse temperature Θ which defines the weight of the configuration. On Fig.( 6.9) we demonstrate several such spectra and overlay them with the result of the fit with a Gaussian distribution. For concreteness, we used in this example N ω = 50 on the range of frequencies [0, 30], the relative error of the input correlation function is of order 10 -2 . 94 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity = 10 0 = 10 1 = 10 2 = 10 3 FIGURE 6.9: Examples of various fit of a mode correlation function which corresponds to the lowest frequency in the LJ spectrum. The fit is performed with the stochastic reconstruction on a discretized range of frequencies, the different colours show few distinct effective inverse temperatures, whilst the number of δ-peaks is fixed N ω = 50. For the purpose of comparison we show the approximation of the spectrum with a gaussian distribution S(ω), which is obtained by a minimization routine.

Similar to the single harmonic oscillator model, increase in Θ considerably narrows the peak. Provided the peak lies within the chosen range of frequencies, the minimization procedure is equivalent to taking Θ → ∞ in the limit N ω → ∞, and thus the algorithm selects the most probable configuration. For N ω finite this is not the case, because of the nonzero ∆ω terms contribution and misalignment of deltapeaks with the maximum of spectral function, unless artificially constructed.

On the Figure (6.10) we show the efficiency of different fits demonstrated by the calculation of validation parameter:

χ 2 val = ∑ m C val iα (τ m ) -Ciα (τ m )
2 , where we put back explicitly the time slice indexes. The validation data set C val iα (τ) is obtained through the separate calculation of the mode correlation function with the same QMC algorithm which we used to obtain original data. The efficiency of the reconstruction, as usual, increases with Θ as demonstrated by the decrease in the value of χ 2 val which we plot here corrected by some constant χ 2 0 that we discussed in Chapter 3. As we expected at large Θ the performance of the reconstruction is quite close to the Gaussian fit of the spectrum. At smaller Θ, however, one can find a spectra which agree better with test data. This could also be expected simply from the fact that the Gaussian function does not reflect properly the features of the underlying spectrum which is asymmetric and have heavy tails. On the other hand, with sufficiently many δ-peaks one can redistribute the weights A(ω i ) associated with them and gain a better description. As we mentioned above, the analysis of the discrete spectrum might be quite cumbersome and in the end one would still need to approximate the resulting spectrum with a smooth function in order to get an estimate for the frequency of the maximum of the spectral density as well as the width of the peak. Keeping also in mind that the Gaussian fit provides a reasonable description for the spectral density, as illustrated by this discussion, we will use it for the analysis of other correlation functions. 96 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity

Discussion of the results

In this section we collect together the results of the classical and quantum simulations announced in the previous parts. We show the calculations for the ensemble of N = 108 particles. The choice of the system size is dictated by the constraints of the fcc lattice (the number of atoms needs to satisfy the relation N = 4p 3 with p integer). Another limitation on quantity N comes from the error estimation analysis, in particular for Monte-Carlo simulation. The QMC calculation of ensemble at low T * requires increasingly more time slices M which for the case of large systems drastically slows down the convergence of observables to the mean values. Thus, we want to keep the value of N sufficiently small. In this section we are interested in obtaining the dispersion relation for ω(k, T * , ω 0 ) as a function of wave-vector k, temperature T * and harmonic frequencies ω 0 . With this objective in mind, the reasonable requirement for a crystal size is to allow at least few discretizations of k vector inside the Brillouin zone in several directions. The ensemble of 108 particles allows to do it while enabling to well control the error for various temperatures. The Monte-Carlo calculations are prone to system size effects which might affect the observables. In order to control it we make a few test calculations of a larger system (namely, with N = 864 in order to be able to compare the frequencies at the same wavevectors) at few temperatures. As we shall see immediately, within the errorbars of calculation the result does not depend on the crystal size. We show it explicitly on the following figures on the example of QMC simulation of systems of different sizes for T * = 0.54.

Mean frequency

We start by the calculation of a mode peak frequency. For each set of degenerate modes we obtain a temperature dependence along with uncertainty estimation. On Fig.( 6.11) we provide the examples for several of these modes. In this example we choose modes that correspond to the wavevectors along the {1, 0, 0} direction and we plot the data for two different wave polarizations, longitudinal and transverse. We put also for the reference the value of frequency for the same mode for strictly harmonic Hamiltonian. As can be seen from the plot, the frequencies of the classical system (depicted in blue) grow linearly with temperature. At the same time the frequency of the quantum system (shown in red) does not exhibit strong temperature dependence and for high enough frequency does not change with temperature.

We found that the system size effect for both simulations is irrelevant and does not affect the results. We demonstrate it on the same set of figures on the example of QMC calculation for the ensemble of N = 864 particles (plotted in green).

We can also gather the data for each mode on a single plot. We do it separately for MD and QMC to avoid cluttering the figure. From Fig. (6.11) it is clear that all the frequencies should be calculated with respect to the mode frequency of a harmonic crystal ω 0 k which in turn depends on few other quantities. On

F) k x = 2π a FIGURE 6 (
.11: The comparison of mean frequency for two sets of modes corresponding to the {1, 0, 0} direction in k-space. On the left we show the transverse branch and on the right -longitudinal. Results of the Molecular dynamics calculation are plotted in blue with the linear approximation of temperature behaviour depicted as a dashed line. The linear regression is made with a fixed point at T = 0 corresponding to the mode frequency of a perfectly harmonic crystal. The output of the QMC simulation is showed in red. As the indication of insensitivity of the QMC data with respect to the change of the system size, we depict in green the result for the crystal of N = 864 particles at T * = 0.54. In order to avoid overlapping with other data we shifted the point to the right.

98 Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity

Specific heat

Using the data from the previous part we are now in a position to analyze the validity of the single-phonon approach. To this end, we study the behaviour of the specific heat. As the baseline for our comparison we measure this quantity in a rather straightforward fashion from the calculation of the total energy of the system [Fig.(6.12)]. Since we are dealing here with the quantum system, energy estimation, in particular the calculation of kinetic energy, is more involved as we discussed in preceding sections, see Sec.(3.3.2). We then approximate the data points with the 4th order polynomial in temperature and calculate its derivative in order to obtain the result for specific heat. Next, we compute the heat capacity in the phonon approximation making use of the formula: 6.38) where in the last line we adopted the expression to reduced units of LJ formalism. With this expression, we compare both sets of data, from QMC as well as from classical MD simulations, the latter having the meaning of the "quasi"-quantum approximation. We show the comparison on Fig. (6.13). For neon, which corresponds to our choice of parameters, we find that the Debye law for heat capacity with the temperature Θ D = 68.5K (similar to [START_REF] Moyano | Lattice dynamics for fcc rare gas solids Ne, Ar, and Kr from ab initio potentials[END_REF]), or T * D = Θ D / = 1.87 describes well our simulation data. This value is reasonably close to the experimentally measured Θ D = 74.6K or Θ * D = Θ D / = 2.03 [START_REF] Fenichel | Low-temperature specific heats of solid neon and solid xenon[END_REF]. We note that the temperature range we have considered is far from the Debye temperature, more precisely 0.05 < T/Θ D < 0.3. On fig.(6.13) we illustrate different calculations of c V and compare them with the experimentally obtained data from [START_REF] Fenichel | Low-temperature specific heats of solid neon and solid xenon[END_REF]. FIGURE 6.13: The comparison between various calculations of heat capacity for quantum neon crystal. The dashed line is calculated via the energy calculation [Fig.(6.12)], the blue and red dots are computed with the mode heat capacity given by expression (6.38). For the calculation of red points the frequency were used frequencies obtained from the QMC simulation whilst for the blue ones -the results of the MD computation. 6.13), for the calculation of red points the frequency were used frequencies obtained from the QMC simulation whilst for the blue ones -the results of the MD computation. For comparison we show the experimental data [START_REF] Fenichel | Low-temperature specific heats of solid neon and solid xenon[END_REF]

c V = 1 N ∂ ∂T ∑ n hω n e hω n /k B T -1 = k B N ∑ n hω n k B T 2 e hω n /k B T [e hω n /k B T -1] 2 = = k B N ∑ n (Qω * n β * ) 2 e Qω * n β * [e Qω * n β * -1] 2 (

Dispersion relation

One can analyze the behavior of frequency shifts using the data on Fig. (6.15) and Fig.(6.16) and deduce an empirical expression for the mode frequency for classical and quantum systems. For simplicity we distinguish frequencies only by their polarization and look for the description which is otherwise uniform for all modes.

In the case of quantum system one observes that there is a clear dependence on ω 0 which we estimate as:

∆ω QMC k ω 0 k = a (L,T) ω 0 k + b (L,T) (6.39)
with a and b being some external parameters that we approximate for a given polarization, transverse (T) or longitudinal(L), denoted on the plots. From our analysis we find that the parameters are 102Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity One can study in the same way the classical system (Fig.(6.16)). In this figure we accounted for the temperature dependence of the quantity by separating the variables. Similarly to the quantum system, the shift in frequency mode also grows with ω 0 . We can now write the expression for frequency shift as (1,0,0) (0,0,0) (1,1,0)

(1/2,1/2,1/2) k harmonic QMC MD FIGURE 6.17: The phonon dispersion relationship for several directions of wave-vector k. Here we plot the results of QMC (red triangles) and MD (blue dots) simulations and demonstrate the prediction of equations (6.39) and (6.40) for frequencies of quantum and classical Lennard-Jones crystal. The shown results correspond to the system of N = 108 particles and T * = 0.41.

These estimations can also be used in order to compare dispersion relationship with the one from the experiment. In order to make the comparison we translate the measured quantities from LJ units to the units used in the experiment. For neon, we have ε = 36.68K/k B = 3.16meV, σ = 2.79Å = 2.79 • 10 -10 m and Q = 0.0918. The experimental measurements in the Ref. [START_REF] Endoh | Lattice dynamics of solid neon at 6[END_REF]] are taken at T = 6.5K, however, as we have seen above, there is no temperature dependence for the phonon frequency for the quantum system, thus we are allowed to juxtapose them. The MD data are, however, more sensitive to temperature, yet they systematically considerably underestimate the frequency values and thus they are shown for illustration purpose. (1,0,0) (0,0,0) (1,1,0)

(1/2,1/2,1/2) QMC MD experiment FIGURE 6.18: The phonon dispersion relationship for several directions of wave-vector k. The results of QMC (in red ) and MD (in blue) simulations along with the prediction of equations (6.39) and (6.40) are compared against the experiments [START_REF] Endoh | Lattice dynamics of solid neon at 6[END_REF].

With this estimations we can now calculate the mode frequencies for every value of wavevector k inside the Brillouin zone using the computation of harmonic mode frequency from Sec.(6.2.2). We show this approximation on Fig.( 6.17) where we show the result for values of k for several commonly used directions, namely {1, 0, 0}, {1, 1, 0} and {1, 1, 1}.

Estimation of the peak width

In the same spirit one can analyze the width of the peak which we define as the parameter γ in the spectral function S(ω) ∝ e -(ω-ω) 2 2γ 2 . We study the same systems as in the previous part, namely, classical and quantum Lennard-Jones crystal of 108 particles at fixed density and for various temperatures. On Fig.( 6.19) we again show the example of modes corresponding to values of wavevector in the direction {1, 0, 0}. For the classical system the width is going to zero with decreasing temperature since the anharmonic corrections become small. For the quantum system it is generally not the case due to the zero-point motion. As we observe, however, for small frequencies this effect is quite small and the temperature behaviour is similar to the classical crystal. At larger frequencies, on the other hand, the quantum effects are more pronounced which leads to the finite width even at zero temperature. As before we also study the effect of the size of the crystal and do not observe significant changes. 104Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity The comparison of peak widths for two sets of modes corresponding to the {1, 0, 0} direction in k-space. On the left we show the transverse branch and on the right -longitudinal. Results of the Molecular dynamics calculation are plotted in blue whilst the output of the QMC simulation is showed in red. As the indication of insensitivity of the QMC data with respect to the change of the system size, we depict in green the result for the crystal of N = 864 particles at T * = 0.54. In order to avoid overlapping with other data we shifted the point to the right.

6.7. Heat conductivity calculation 105

Heat conductivity calculation

We can combine the results for the peak frequencies and the width of the peak in order to compute heat conductivity λ. Let us remind ourselves of the equation for λ (Eq.(6.23)):

λ = k B σt 0 V * ∑ k v * 2 k τ * k ∂ ¯ * k ∂T * (6.41)
As we can see it combines several ingredients. First, one needs the group velocity ∂ω ∂k for each mode in the system. We can easily obtain this quantity using the approximations from Sec.(6.6.3) which reduce the derivative of ω MD,QMC to the derivative of harmonic frequencies ω 0 which one can straightforwardly calculate for any k inside the Brillouin zone using the algorithm from Section (6.2.2). The second constituent is the decay rate of mode τ k which can be calculated via the width of the peak as τ k = 1/2γ k for the case of Lorentzian peak or as τ k = 1/2 √ 2γ k for the case of Gaussian spectral function. The last ingredient to the calculation of heat conductivity is the modal heat capacity c k given by the expression (6.38)

c k = k B (Qω * k β * ) 2 e Qω * k β * [e Qω * k β * -1] 2 .
(6.42)

Again, the last equation is written in Lennard-Jones units. The mode heat capacity is calculated with the same formula given by Eq.(6.42) with the replacement of ω by the mode peak frequency in both cases.

106Chapter 6. Phonon lifetimes in rare gas solids, implication for heat conductivity

Combining together the data we obtain the result for λ for various temperatures and plot them on Fig.(6.20) along with the experimental data from [START_REF] Weston | Temperature and volume dependence of the thermal conductivity of solid neon[END_REF] (taken for molar volume v M = 13.35cm 3 /mol which translates to density ρσ 3 ≈ 0.98 -reasonably close to the value ρσ 3 = 0.965 used in our simulation).

In order to compare our results with experimental values given in Ref. [START_REF] Weston | Temperature and volume dependence of the thermal conductivity of solid neon[END_REF], we have to go back from Lennard-Jones units to conventional ones for experiments. Our highest temperature, T * = 0.55 corresponds to T = 20K. There, the heat conductivity λ * = λσt 0 /k B is around 15 in LJ units translates to λ = 0.027λ * [Wm -1 K -1 ] which gives ≈ 0.4Wm -1 K -1 . This value compares approximately to the experimental one obtained for v M = 13.36 in Ref. above. However, for lower temperatures, there is a qualitative disagreement: whereas our calculated heat conductivity drops with decreasing temperature, the experimental one increases! Indeed, the drop of the heat conductivity lowering temperature is unphysical, as the mean-free path should rather rise at low temperature.

In order to avoid introducing artifacts in the calculation from the small size of the crystal we study the heat conductivity for systems of larger sizes. We observe an increase of around 20% in values of κ which is clearly not sufficient to account for the deviations of the calculated conductivity from the experiment. The expression (6.23) consists of three types of ingredients: phonon heat capacity, phonon group velocity and phonon lifetimes. We found that the dispersion relationship for frequencies that we obtained in Sec.6.6.1 is in good correspondence with experimentally observed values. Furthermore, we checked that the heat capacity calculated as the sum over individual modes (eq.(6.38)) produces the result which is close to specific heat obtained from the total energy of the crystal. This also indicates that the study of system properties in terms of phonons is a reasonably good approximate approach. All these considerations indicate that the discrepancy in the calculation of heat conductivity comes from our estimation of phonon life time. Going back to the calculation of τ k which we performed in the Sec.(6.6.4) we note that we estimated the phonon life time from the decay of mode correlation function. Thus, we included implicitly in this calculation all possible scattering processes of mode A k . However, the only processes that contribute to the energy transport inside the crystal are the phonon scatterings to the modes outside of the Brillouin zone. This is the Umklapp process, and clearly, at high temperature this process becomes dominant, and phonon mean free path approaches the mean free path of phonons scattering outside of the Brillouin zone. At low temperatures, however, it is not the case and the phonon life time is smaller than the life time of "transport phonons". In our calculation we do not distinguish between various scattering processes which results in incorrect values of τ k . The proper way to study heat conductivity is by following Green-Kubo formalism and analyze heat current correlation functions.

Conclusion and outlook

The Nuclear Quantum effects are important for interatomic interactions inside the insulating crystal and have to be treated carefully. In order to analyze quantitatively their effect we introduce the computational scheme based on the calculation of imaginary time correlation functions obtained by path integral Monte Carlo simulations. One can tackle the problem from several angles and in general it may involve computation of correlation functions containing current or momentum operators. Thus, we have described a general strategy for wisely expressing improved estimators with reduced statistical variance for such correlations. Next, we have introduced an inversion procedure in order to obtain corresponding spectral functions. The algorithm is based on a stochastic maximum entropy method, a Bayesian approach commonly used for such problems. The outcome of these procedures is, in general, strongly dependent on the parameters of the simulation, as we have illustrated on several examples of the oscillator spectra employing different values for the effective inverse temperature, Θ, as well as different choices for the grid discretization, N ω , or offset, δω. These models proved a challenging benchmarks for the spectral reconstruction due to the sharp undamped delta-functions they contain and the necessity of high accuracy of the imaginary time correlation function.

Pure Bayesian approaches suggest to eliminate the parameters dependence by using the most general and flexible model for the spectral density, e.g., a large value for N ω , together with Θ = 1 to encompass all possible solutions consistent with the data. On the contrary, in our case studies of oscillator models we have shown that the spectra corresponding to these standard choices exceedingly suffer from the usual problems of all maximum entropy reconstructions: broadening or merging of peaks, smoothing out any sharp features in the underlying exact spectrum.

In fact, path integral Monte Carlo data are strongly correlated in imaginary time, which undermines a true justification of the Bayesian choice Θ = 1. Different values of Θ may therefore be considered to efficiently approximate the true, unknown likelihood function. On the other hand, the use of flexible models for the spectral function, containing a large number of parameters, possibly introduces a large amount of entropy into the Bayesian inversion, such that different parametrizations in general strongly modify the results.

In order to address these difficulties we developed a validation procedure to quantitatively control any parameter dependence of the Bayesian inversion. Our proposal is based on the quantity χ val which measures the accuracy of the fitted correlations with respect to independent data (i.e. not involved in the maximum entropy inversion), which provides an efficient and readily applicable method to select the optimal choice of parameters, corresponding to the lowest value of χ val .

We have shown that the validation step identifies a discrete set of two delta functions in the case study of the single harmonic oscillator, and provides indications towards the correct asymmetric sharp edges in the case of an underlying continuous frequency spectrum. In the case of the double well potential we also demonstrated the utility of sparse frequency grids with varying (non uniformly) the spacing between nodes for optimal description of discrete spectral functions.

Chapter 7. Conclusion and outlook

We then put this machinery to use for the analysis of insulating crystal modeled by the nearest-neighbor Lennard-Jones potential. Using the normal mode description we analyze mode correlation functions along with the corresponding spectra. For the conclusive study we perform the quantum Monte-Carlo simulation and compare the results with classical molecular dynamics simulation. Non-quadratic terms in the interaction potential led in both cases to the shift and broadening of delta-peaks of the harmonic crystal spectrum which indicated the interaction and decay of modes. For both types of system we obtained well-defined characteristic Gaussian peaks and their respective frequencies of the spectral function maxima and widths. Using the validation procedure we verified that the usual minimization of the spectral function provides good results both for classical and quantum cases thus allowing for more accurate evaluation of peaks and broadening.

For the mentioned systems we observed two distinct behaviours. The atoms of neon (at the observed parameters) used in the simulation constitute a system that displays strong quantum effects. For the mean mode frequency this leads to the weak dependence on temperature, whilst its classical counterpart grows linearly with T. For the broadening of the peak this results in a different temperature behaviour and non-vanishing of the decay at low T. Taking into account both of these effects along with the change of density of state, we calculate the corrections to the classical heat conductivity.

A possible extension of this work includes the application of the Green-Kubo formalism outlined in the work for computation of heat current correlation function. This would then allow for more accurate analysis of the utility of the single mode relaxation time approach. Following this route, another possible application of this method would be in the domain of glassy materials: recent studies [START_REF] Simoncelli | Unified theory of thermal transport in crystals and glasses[END_REF][START_REF] Isaeva | Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach[END_REF] showed that the Green-Kubo method coupled with the perturbative treatment of anharmonic effects can be remarkably efficient for the determination of heat conductivity at low temperature in systems such as amorphous silicon. One can envision that the current work can be extended to arbitrary region of temperatures and stronger anharmonic effects, using the path integral approach to go beyond perturbative treatment of anharmonicity, and employing the spectral reconstruction techniques developed here.

1. 1

 1 Comparison of heat conductivity for solid Ne obtained from experiment with the analytical calculations -in red classical approach and in black -Debye model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Heat capacity of silicon: comparison of experiment with Debye model 2.1 Example of the regularization procedure for a kernel of the inverse problem. Here we perform a 3 parameter fit and plot values of χ 2 for the fixed set of ω p and corresponding A p as defined by eq.(2.37) as a function of regularization parameter λ. Set of correlation function data points and corresponding errors are taken from the first test case from Sec.(4.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Calculation of integral . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi 4 .

 4 15 Main panel: Comparison of the χ 2 val obtained from our validation of the reconstructed spectral function for various values of Θ and N ω for C pV (τ) at β = 10. The area of the circles is proportional to the corresponding value of χ 2 val . Inset: χ 2 val as a function of Θ, at the indicated values of N ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Spectral reconstruction for C pV (τ) at β = 3, obtained at the indicated values of the discretization, N ω , for a fixed Θ = 1. The filled rectangles are centered at the positions of the two delta-functions of the exact spectrum for ∆ω = 1, with an area corresponding to their respective weights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Spectral reconstructions from C pV (τ) at β = 3 for N ω = 5 using different values of Θ. The filled rectangles are centered at the positions of the two delta-functions of the exact spectrum, with an area corresponding to their respective weights. . . . . . . . . . . . . . . . . . . . 4.18 Main panel: χ 2 val from the validation procedure of the reconstructed spectral function at the corresponding values Θ and N ω for C pV (τ) at β = 3. The area of the circles is proportional to the value of χ 2 val . Inset: χ 2 val as a function of the effective temperature Θ, at the indicated values of N ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 The double well potential described by the Hamiltonian Eq.(5.1). The values of parameters m, α and γ are fixed to 1. We also indicate the values of first energy levels of the system. . . . . . . . . . . . . . . . . 5.2 Spectral reconstruction of C β=8 xx (τ) for the double well potential, at the indicated values of the discretization N ω and Θ reconstructed on the range of ω = [0, 5]. Value of Θ corresponds to the minimum of χ 2 val for N ω = 25 model as shown in the Fig. 5.3. Red spectral function corresponds to a N ω = 10 lattice shifted by δω = 0.25. The exact spectrum of C β=8 xx (τ) at β = 8 is plotted as gray hatched area. For the purpose of illustration all delta-functions are shown with the finite width ∆ω = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Θ-dependence of the χ 2 val calculated from the validation step of the reconstructed spectral functions for C β=8 xx (τ) for frequency discretization N ω = 10 and N ω = 25. Spectral functions for non shifted N ω = 10 grid (shown at Θ = 15 as blue spectrum on Fig. 5.2) gives considerably worse result and the corresponding χ 2 val are not plotted here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Spectral reconstruction of C β=1 xx (τ) for the double well potential, at the indicated Θ. Values of Θ correspond to the minima of χ 2 val for each model respectively as shown in the Fig. 5.6. We compare models with frequency discretization N ω = 10 and N ω = 25, red spectral function corresponds to a N ω = 10 lattice shifted by δω = 0.25. The gray hatched area indicates the exact spectral function at β = 1. All delta-functions are shown with the finite width ∆ω = 0.01. . . . . . . 5.5 Θ-dependence of the χ 2 calculated from the reconstruction procedure for C β=1 xx (τ). Here we compare shifted (δω = 0.25) and non shifted lattices which have frequency discretization N ω = 10 and N ω = 25. . . xii 5.6 Θ-dependence of the validation ∆χ 2 val ≡ chi 2 valchi 2 0 calculated from the validation step of the reconstructed spectral functions for C β=1 xx (τ). Value of chi 2 0 corresponds to a reconstruction with N ω = 10 shifted by δω = 0.25 at Θ = 10. Here we compare shifted (δω = 0.25) and non shifted lattices which have frequency discretization N ω = 10 and N ω = 25. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Spectral reconstruction of C cont

5. 1

 1 Ground energy and lowest excited energy states of the particle inside the double well potential. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Lowest excited states frequency values, ω n = E n -E n-1 , of the spectrum of C β xx (τ) for the double well potential. . . . . . . . . . . . . . . . 6.1 Parameters of the simulation of the neon crystal. . . . . . . . . . . . . . 6.2 Parameters of the test simulation of the neon crystal for the comparison with [Cuccoli et al., 1993]. . . . . . . . . . . . . . . . . . . . . . . . .

  FIGURE 1.1: Comparison of heat conductivity for solid Ne obtained from experiment with the analytical calculations -in red classical approach and in black -Debye model.

  notations a = 2m h2 ∆τ + mω 2 0 ∆τ and b = -m h2 ∆τ

FIGURE 4 . 1 :

 41 FIGURE 4.1: The dependence of the mean squared displacement x 2 on the imaginary time discretization ∆τ = β/M for a single quantum harmonic oscillator. Here we compare analytical results obtained with Eq.(4.5)(orange line) and Eq.(4.18)(green line) against the results of a PIMC calculation (blue dots).

FIGURE 4 . 2 :

 42 FIGURE 4.2: Position correlation x(τ)x(0) as a function of imaginary time for a single quantum harmonic oscillator. Solid line corresponds to the diagonalization of the Hamiltonian (eq.(4.5)) and dashed line is given by eq.(4.18). Dots correspond to the MC simulations for different imaginary time discretizations β/M. Value of β is fixed and equal to 1 for all the data on the figure. Errorbars are smaller than the size of the dots.

FIGURE 4

 4 FIGURE 4.3: MC results for the normalized correlation function(dots) as a function of τ plotted for several values of β. We compare them with the exact result(dashed line) which is obtained from Eq.(4.18)

FIGURE 4 . 4 :

 44 FIGURE 4.4: Spectral reconstruction of position correlation function(eq.4.18). Calculation was performed on a range of ω = [0, 2] with N ω = 50 δ-peaks. We used 2 types of constraint for spectral density: "hard" cutoff A i ≥ 0 and "soft" cutoff A i ≥ -0.01 is plotted on the Fig.(4.4). As we discussed earlier the exact spectrum of correlation function consists of only one frequency at ω = 1 in the units of the plot. With the accuracy of the data for x(τ) x(0) we were able to identify a single peaked function centered around ω = 1. The width of the curve is determined by the accuracy of our sample data. For the considered model both constraint ("hard" and "soft" cutoff) result in a practically identical spectral density and the convergence of the MC algorithm used for the stochastic reconstruction was not affected by the change of this constraint.

  Fig.(4.5) demonstrates the comparison between these results. As can be seen, mean values that are obtained with all three of them converge to the result predicted with analytical calculation. The difference between them appears when we consider variance of the output data and as consequence accuracy of these data. The simple equation (4.29) produces the most noisy sample points, with a sizeable portion of samples lying outside of the errorbars from the theoretical curve(Fig.(4.6)). The relative error σ(τ i )/F(τ i ) is huge in comparison to the other estimators.

FIGURE 4 . 5 :

 45 FIGURE 4.5: Comparison of different estimators for momentum correlation function for the same MC simulation of harmonic oscillator.Here we can see results for "naive" estimator eq.(4.29) (orange curve) and the calculation made with the estimator from eq.(3.32) written in the form (3.37) adopted to the momenta correlations (plotted in green).

  FIGURE 4.6: Comparison of relative errors between different estimators (naive and virial ones) of momentum correlation function. Figure (4.6a) shows the result for the whole range of imaginary time [0, β] and figure (4.6b) demonstrates the behaviour of the error for small times. With correlation function C(τ) τ→β/2

  (4.7) and Fig.(4.8). Plotted are the reconstructions for the correlation functions from the previous section, on the first figure we demonstrate the spectral densities whilst the second one shows the corresponding imaginary time correlation functions. For the latter figure we plot the deviation of momenta correlation functions (both for MC data and reconstruction) from the exact function given by the equation (4.21

FIGURE 4 . 7 :

 47 FIGURE 4.7: Reconstruction of spectra of momentum correlation function (eq. 4.21). Calculation was performed on a range of ω (0,5) with discretization N ω = 50. We compare the results obtained with 2 different estimators: "naive" and "virial" ones.

FIGURE 4 . 8 :

 48 FIGURE 4.8: Comparison of the reconstructed imaginary time correlation functions with the MC data used for the analysis. We plot the result as the deviation from the exact momenta correlation function (4.21). The errorbars are estimated from the analysis of the MC data. The reconstructed correlation functions are obtained from the spectra depicted on the Fig.(4.7) with the same color notations.

FIGURE 4 . 9 :

 49 FIGURE 4.9: Relative discretization error, 1 -C PI pV (τ)/C exact pV (τ), between the path integral, C PI pV (τ), and the exact results, C exact pV (τ), for the energy current correlation function, as a function of hω 0 β/M. We show the data corresponding to the imaginary times τ = 0 and τ = β/2, and indicate with symbols and solid lines the results for hβ = 3 and 10, respectively.

  FIGURE 4.10: Difference between the exact correlation function C exactpV (τ) and the values obtained by PIMC sampling, C MC pV (τ), of a path with M = 100 time slices, for hω 0 β = 1, illustrating the variance reduction obtained by the improved estimator discussed in Sect. 3.3.2. We show with line-points the primitive estimator and with the continuous line the improved estimator, both using the same path-integral Monte Carlo data.

  FIGURE 4.11: Difference between the exact correlation function C exact pV (τ) and the values obtained by PIMC sampling, C MC pV (τ), of a path with M = 30 time slices, for hω 0 β = 3, illustrating the variance reduction obtained by the improved estimator discussed in Sect. 3.3.2.We show in blue the calculations the 'primitive' estimator and in red -the improved estimator, both using the same path-integral Monte Carlo data.

Chapter 4 .FIGURE 4

 44 FIGURE 4.14: N ω -dependence of the χ 2val extracted from the validation step of the reconstructed spectral functions for C pV (τ) at β = 10 and with Θ = 1. Squares and triangles correspond to shifted grids: for N ω = 5 red square shows the shift δω = 0.25 and the green one δω = 0.5, for N ω = 10 plot shows δω = 0.1 as a red triangle and δω = 0.25 as a green one.

  FIGURE 4.15: Main panel: Comparison of the χ 2 val obtained from our validation of the reconstructed spectral function for various values of Θ and N ω for C pV (τ) at β = 10. The area of the circles is proportional to the corresponding value of χ 2 val . Inset: χ 2 val as a function of Θ, at the indicated values of N ω .

FIGURE 4 .

 4 FIGURE 4.16: Spectral reconstruction for C pV (τ) at β = 3, obtained at the indicated values of the discretization, N ω , for a fixed Θ = 1. The filled rectangles are centered at the positions of the two deltafunctions of the exact spectrum for ∆ω = 1, with an area corresponding to their respective weights.

  FIGURE 5.2: Spectral reconstruction of C β=8 xx (τ) for the double well potential, at the indicated values of the discretization N ω and Θ reconstructed on the range of ω = [0, 5]. Value of Θ corresponds to the minimum of χ 2 val for N ω = 25 model as shown in the Fig. 5.3. Red spectral function corresponds to a N ω = 10 lattice shifted by δω = 0.25. The exact spectrum of C β=8 xx (τ) at β = 8 is plotted as gray hatched area. For the purpose of illustration all delta-functions are shown with the finite width ∆ω = 0.01.

  FIGURE 5.3: Θ-dependence of the χ 2 val calculated from the validation step of the reconstructed spectral functions for C β=8 xx (τ) for frequency discretization N ω = 10 and N ω = 25. Spectral functions for non shifted N ω = 10 grid (shown at Θ = 15 as blue spectrum on Fig. 5.2) gives considerably worse result and the corresponding χ 2 val are not plotted here.

FIGURE 5 . 4 :

 54 FIGURE 5.4: Spectral reconstruction of C β=1 xx (τ) for the double well potential, at the indicated Θ. Values of Θ correspond to the minima of χ 2 val for each model respectively as shown in the Fig. 5.6. We compare models with frequency discretization N ω = 10 and N ω = 25, red spectral function corresponds to a N ω = 10 lattice shifted by δω = 0.25. The gray hatched area indicates the exact spectral function at β = 1. All delta-functions are shown with the finite width ∆ω = 0.01.

FIGURE 5 FIGURE 5

 55 FIGURE 5.5: Θ-dependence of the χ 2 calculated from the reconstruction procedure for C β=1 xx (τ). Here we compare shifted (δω = 0.25) and non shifted lattices which have frequency discretization N ω = 10 and N ω = 25.

-

  (r) dr where e(r) is the density of the Hamiltonian such that H 0 = e(r)dr. Expanding T(r) one can write P ∝ e -

  FIGURE 6.4: Kinetic energy per particle e K of neon as a function of temperature T. The values of density of crystal at each temperature are as indicated in Table6.2. The red data indicate the present simulation, whilst the blue ones -reference data from paper[START_REF] Cuccoli | Monte Carlo computations of the quantum kinetic energy of rare-gas solids[END_REF]. The reference data were approximated in the limit ∆τ → 0.

  Fig.(6.6) we illustrate the dependence of pressure p on temperature on the example of systems with densities ρσ 3 = 0.965 and ρσ 3 = 0.944. For the purpose of comparison we indicate on the same plot the values of pressure at the same parameters for the system described by the classical Hamiltonian. Using these data one can now analyze the thermal expansion of lat-

88FIGURE 6 . 5 :

 65 FIGURE 6.5: Pressure as a function of temperature for neon crystal described by quantum LJ potential. The figure shows the data for systems at indicated constant densities ρσ 3 .

  FIGURE 6.8: Example of a fit performed on the average mode correlation functions C AA (t).Figure shows the spectral density for the same mode calculated for the system at various temperatures. Different dots correspond to the spectra obtained by transforming C AA (t) whereas dashed line denotes the fit.

  FIGURE 6.8: Example of a fit performed on the average mode correlation functions C AA (t).Figure shows the spectral density for the same mode calculated for the system at various temperatures. Different dots correspond to the spectra obtained by transforming C AA (t) whereas dashed line denotes the fit.

  FIGURE 6.10: Comparison between different calculation of spectra. The dots show the stochastic spectral reconstruction with different effective inverse temperature Θ with ∆ω = 0.06. The blue line demonstrates the fit result obtained from minimizing the functional C AA iα (τ) -dωK(ω)S(ω, ω 0 , γ) with a Gaussian spectrum S(ω, ω 0 , γ).

  Fig.(6.16) and Fig.(6.15) we plot the shifts of mode frequencies with respect to ω 0 as a function of the same harmonic frequency ω 0 for Monte-Carlo and Molecular dynamics simula-

100Chapter 6 .

 6 Phonon lifetimes in rare gas solids, implication for heat conductivity 0.05 0.10 0.15 0.20 0.25 0.30 0.35 T / T D

FIGURE 6 .

 6 FIGURE 6.14: The comparison between various calculations of heat capacity for quantum neon crystal. The dashed line shows the heat capacity calculated with the Debye law at Θ D = 65.8K, or in LJ units T * D = 1.87. The blue and red dots are computed with the mode heat capacity given by expression (6.38). As for the Fig.(6.13), for the calculation of red points the frequency were used frequencies obtained from the QMC simulation whilst for the blue ones -the results of the MD computation. For comparison we show the experimental data[START_REF] Fenichel | Low-temperature specific heats of solid neon and solid xenon[END_REF]

  FIGURE 6.15: We demonstrate here the shift in mode frequencies ω k for the quantum Lennard-Jones crystal from the harmonic ones ω 0 k plotted against the ω 0 k . The different colors represent longitudinal and transverse branches. The data as well as the uncertainty are shown divided by the harmonic value of the same mode.

  FIGURE 6.16: The figure shows the shift in mode frequencies ω k for the classical Lennard-Jones crystal from the harmonic ones ω 0 k plotted against the ω 0 k . The different colors again represent longitudinal and transverse branches. The data along with the uncertainty are shown divided by the corresponding temperature and harmonic frequency of the same mode.

  data coefficients c (L,T) we obtain c (T) = 0.2 and c (L) = 0.3.

  FIGURE 6.19: The comparison of peak widths for two sets of modes corresponding to the {1, 0, 0} direction in k-space. On the left we show the transverse branch and on the right -longitudinal. Results of the Molecular dynamics calculation are plotted in blue whilst the output of the QMC simulation is showed in red. As the indication of insensitivity of the QMC data with respect to the change of the system size, we depict in green the result for the crystal of N = 864 particles at T * = 0.54. In order to avoid overlapping with other data we shifted the point to the right.

  FIGURE 6.20: The calculation of heat conductivity from the data of Sec.(6.6). In blue we show the results for the Molecular dynamics calculation whilst in red we plot the quantum Monte-Carlo results. The mode heat capacity is calculated with the same formula given by Eq.(6.42) with the replacement of ω by the mode peak frequency in both cases.

  • • • = 1 . . . P and α, β, • • • = 1 . . . M -1 In this form we can write a solution for A : A

  2.127

	E 0	0.14719
	E 1	0.87197
	E 2	2.12692
	E 3	3.58846
	E 4	5.2399

FIGURE 5.1: The double well potential described by the Hamiltonian Eq.(5.1). The values of parameters m, α and γ are fixed to 1. We also indicate the values of first energy levels of the system.

TABLE 5

 5

.1: Ground energy and lowest excited energy states of the particle inside the double well potential.

TABLE 5 .

 5 2: Lowest excited states frequency values, ω n = E n -E n-1 , of the spectrum of C

β xx (τ) for the double well potential.

  Example of a mode correlation function C AA (t) obtained from the MD calculation for the classical system with Lennard-Jones pair interaction. Various colors represent 12 different modes, all having the same degerate frequency. The black line shows the correlation that is obtained by averaging over all these degenarate modes.

		0.0010
		0.0005
	C AA (t)	0.0000
		0.0005
		0.0010
		0	5	10	15 t	20	25	30
		FIGURE 6.7:

Powered by TCPDF (www.tcpdf.org)Powered by TCPDF (www.tcpdf.org)

systematic discretization error is smaller than the statistical error for our simulation time, so it can be safely neglected. The main constraint for the reconstruction comes from the imaginary time interval [0, 1]. The relative error of the PIMC data corresponding to these values of τ is of O(10 -2 ). For larger τ the relative error becomes comparable with the data due to the fact that C pV (τ) approaches 0 with τ → β/2. We start by considering the case β = 10. First, we evaluate the effect of the grid size, N ω , on the reconstruction. In Fig. 4.12 we show the spectra obtained for various values of N ω , keeping a fixed Θ = 1. As mentioned above, there is no apriori argument guiding the most appropriate parametrization of the spectrum. In the following we analyze the accuracy of the spectral reconstruction by comparing the values of χ 2 val defined in Eq. (2.45), using an independent test data set. This is obtained within an additional MC simulation of the correlation function, with the same parameters as the original one. We also consider a data set of the same size, P , as the one that was used to produce C pV (τ k ).

In Fig. 4.14 we show χ 2 val as a function of the number of grid points. Clearly, increasing the number of coefficients A(ω i ) of Eq. (2.38) does not lead to a better spectral reconstruction. In contrast, by introducing more degrees of freedom, one increases the entropy, and the spectral weight is smeared out excessively. In Fig. 4.14 we also show the effect on χ 2 val of a shift δω. As expected, shifting the nodes away from ω 1 = 1 and ω 2 = 3, which are the only frequencies present in the exact spectrum determined by Eq. (4.49), deteriorates the accuracy of the spectrum obtained through the validation step. The second parameter determining the quality of the statistical maximum entropy reconstruction is the effective temperature, Θ. In Fig. 4.13 we show the behaviour of the spectral function for a chosen ω-grid at the indicated values of Θ. As expected from Eq. (2.38), by increasing Θ the result approaches the most probable configuration that describes the correlation function C pV (τ), reducing entropic effects. In Fig. 4.15 we combine the above results for different pairs of parameters (Θ, N ω ), and plot the corresponding χ 2 val . Our validation procedure therefore strongly points to using models with a smaller number of delta functions combined with large values of Θ 1 for the spectral reconstruction. Based on the comparison with the exact spectrum, this choice is also clearly the one that leads to the description of the spectrum in closest agreement with the exact prediction. We conclude that the use of χ 2 val indeed seems to provide an unbiased estimate of the quality of the reconstruction.

We now consider the spectral reconstruction for C pV (τ) at β = 3, again clarifying the influence of Θ and of the lattice discretization N ω . In Figs. 4.16 and 4.17 we show selected examples of the resulting spectra. In contrast to the case β = 10, we now observe in general a much stronger broadening of the peaks, which prevents us from resolving the two peak structure for Θ = 1, even for sparse ω-grids. However, when combining sparse grids with sufficiently large Θ in the inversion, one improves towards the correct two peaks structure, as can be seen in Fig. 4.17. The data shown in Fig. 4.18 also indicate that this choice indeed corresponds to the lowest values of χ 2 val , confirming the validity of this indicator. We also note that, for large Θ, the values of χ 2 val tend to exhibit a minimum or weak oscillations, that are probably indicative of overfitting. As a consequence, considering larger values of Θ