Markus Holzmann

Vladislav Efremkin

Jean-Louis Barrat

Stefano Mossa

Time correlation functions for quantum systems: validating Bayesian approaches for harmonic oscillators and beyond

come

Time correlation functions for quantum systems: validating Bayesian approaches for harmonic oscillators and beyond Vladislav Efremkin, 1, a) Jean-Louis Barrat, 2, b) Stefano Mossa, [START_REF] Mizuno | Relation of vibrational excitations and thermal conductivity to elastic heterogeneities in disordered solids[END_REF] and Markus Holzmann 4 1) Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France 2) Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France 3) CEA, IRIG-MEM, Univ. Grenoble Alpes, 38000 Grenoble, France 4) Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France (Dated: 8 September 2021)

The quantum harmonic oscillator is the fundamental building block to compute thermal properties of virtually any dielectric crystal at low temperatures in terms of phonons, extended further to cases with anharmonic couplings, or even disordered solids. In general, Path Integral Monte Carlo (PIMC) or Molecular Dynamics (PIMD) methods are powerful tools to determine stochastically thermodynamic quantities without systematic bias, not relying on perturbative schemes. Addressing transport properties, for instance calculating thermal conductivity from PIMC, however, is substantially more difficult. Although correlation functions of current operators can be determined by PIMC from analytic continuation on the imaginary-time axis, Bayesian methods are usually employed for the numerical inversion back to real-time response functions. This task not only strongly relies on the accuracy of the PIMC data, but also introduces noticeable dependence on the model used for the inversion. Here, we address both difficulties with care. In particular, we first devise improved estimators for current correlations which substantially reduce the variance of the PIMC data. Next, we provide a neat statistical approach to the inversion problem, blending into a fresh workflow the classical stochastic maximum entropy method together with recent notions borrowed from statistical learning theory. We test our ideas on a single harmonic oscillator and a collection of oscillators with a continuous distribution of frequencies, and provide indications of the performance of our method in the case of a particle in a double well potential. This work establishes solid grounds for an unbiased, fully quantum mechanical calculation of transport properties in solids.

I. INTRODUCTION

Atomic nuclei in condensed phases behave, in many cases, as quantum objects. For instance, Nuclear Quantum Effects are responsible for the heat capacity problem, i.e., the deviation from the classical Dulong and Petit law for the heat capacity of solids at low temperatures. The solution of this issue eventually led to the development of the harmonic theory of solids, an accurate quantum theory that lets us to compute their thermal properties at temperatures lower than the Debye temperature, and can be corrected to account for anharmonic effects [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF] . By reducing the description of an insulating solid to a set of independent harmonic oscillators, the phonons, weakly interacting through anharmonic couplings, this theory also provides a framework for the computation of transport properties, in particular of heat conductivity. In contrast to the very high accuracy that can be achieved for thermal properties, however, the computation of transport features is sensibly more delicate and often requires ad hoc approximations for the lifetime of phonons, which is limited by phononphonon scattering processes and the presence of defects. The general framework of the harmonic theory of solids, originally developed for crystals, can be adapted to disordered solids. This is at the expenses of employing a numerical approach to characterize the harmonic eigenmodes, that replace phonons and are no longer determined by symmetries. Again, this procedure can be efficiently employed to determine thermal properties, while its application to transport is much more limited. a) Electronic mail: vladislav.efremkin@univ-grenoble-alpes.fr b) Electronic mail: jean-louis.barrat@univ-grenoble-alpes.fr Very often these properties are indeed calculated via classical statistical mechanics approaches (based on classical Molecular Dynamics simulations), whose results are next empirically corrected to account for quantum effects (see, among others, [START_REF] Mizuno | Relation of vibrational excitations and thermal conductivity to elastic heterogeneities in disordered solids[END_REF] ). We also note that, in systems (ordered or disordered) involving light nuclei (e.g., hydrogen in solid ice), the large wavelength associated with light atoms makes the harmonic approximation itself inappropriate. Therefore, an exact calculation should in general be considered even for thermal properties, or for the determination of phase boundaries [START_REF] Bronstein | Quantum versus classical protons in pure and salty ice under pressure[END_REF] .

The harmonic theory of crystalline solids undoubtedly constitutes a remarkable achievement, as many results can be obtained based on an almost fully analytical approach. However, the above limitations in computing transport properties or in applying the theory to disordered structures, point to the necessity of numerical approaches. It would therefore be highly desirable to develop a numerical methodology that could fully take into account the quantum nature of atomic nuclei, allowing us to determine without approximations both thermal and transport properties of any insulating solid.

When interested in thermal properties, an exact numerical method that encompasses all quantum aspects and is valid at any temperature, independently of the strength of anharmonic effects, involves the path integral representation of the partition function [START_REF] Barker | Quantum-statistical Monte-Carlo method -path integrals with boundary conditions[END_REF][START_REF] Chandler | Exploiting the isomorphism between quantum-theory and classical statistical-mechanics of polyatomic fluids[END_REF][START_REF] Herman | On path integral Monte-Carlo simulations[END_REF][START_REF] Pollock | Simulation of quantum many-body systems by path-integral methods[END_REF] . In the absence of exchange effects (a reasonable hypothesis in most common solids), the determination of thermodynamic properties at the inverse temperature β = (k B T ) -1 involves the sampling of an equivalent system where each quantum particle is replaced by a discretized "path" consisting of M "imaginary time slices". The method becomes exact in the limit of large M, and the sampling of N quantum degrees of freedom at temperature T turns out to arXiv:2105.09575v2 [cond-mat.stat-mech] 7 Sep 2021 be equivalent to that of N × M classical degrees of freedom at temperature M × T . This sampling can be achieved efficiently using Monte Carlo or Molecular Dynamics methods, leading to the PIMC and PIMD methods, respectively [START_REF] Tuckerman | Statistical mechanics: theory and molecular simulation[END_REF] .

Computation of transport properties is more problematic. The standard Green and Kubo statistical mechanics approach to transport coefficients [START_REF] Green | Markoff random processes and the statistical mechanics of time-dependent phenomena[END_REF][START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF][START_REF] Luttinger | Theory of thermal transport coefficients[END_REF] , obtains the heat conductivity tensor κ in a system of volume V at temperature T from a time correlation function of the energy current operator J as,

κ αβ = 1 V k B T 2 ∞ 0 dt J α (t)J β (0) . (1)
Unfortunately, the path integral method provides directly static (time-independent) quantities only. A possible solution to this problem has been identified long ago [START_REF] Thirumalai | On the calculation of time correlationfunctions in quantum-systems -path integral techniques[END_REF] , by noting that the PIMC approach can rather supply the analytical continuation of the correlation functions on the imaginary time axis, simply by computing the correlation between two imaginary time slices along the path. The power spectrum, S AB (ω), of a real time correlator, C AB , between two operators A and B can then be obtained in an apparently straightforward manner by using the identity,

C AB (iτ) = ∞ 0 dω S AB (ω)e -hωτ + S BA (ω)e -hω(β -τ) . ( 2
)
While Eq. ( 2) in principle allows one to obtain S based on the data for C(iτ), with τ in [0, β ], it is well known that the inversion problem is ill-posed, in the sense that determining S with high precision is an extremely difficult task, even if C is known with excellent accuracy. For this reason, the approach pioneered by a few groups in the eighties within the framework of path integral calculations did not spread widely. As an alternative, real-time correlation functions of the centroid, or ring-polymer dynamics computed within PIMD, often provide accurate spectra [START_REF] Perez | A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals[END_REF] . They lack, nonetheless, a rigorous justification of the underlying semi-classical approximation for obtaining low temperature transport properties, and we will not go this route in what follows.

Many recent studies obtained in various fields [START_REF] Sandvik | Stochastic method for analytic continuation of quantum monte carlo data[END_REF][START_REF] Bertaina | Statistical and computational intelligence approach to analytic continuation in quantum monte carlo[END_REF][START_REF] Levy | Implementation of the maximum entropy method for analytic continuation[END_REF][START_REF] Goulko | Numerical analytic continuation: Answers to well-posed questions[END_REF][START_REF] Kora | Dynamic structure factor of superfluid 4 He from quantum monte carlo: Maximum entropy revisited[END_REF] , however, indicate that the present computing capabilities should by now allow us to carry out the original program of analytical continuation from imaginary to real time spectrasatisfactorily, by addressing the two major (and related) difficulties: i) to obtain with high accuracy the imaginary time correlation, in particular for current operators which suffer from the well known issue of diverging variance 7 in the limit of large M; and ii) to solve the ill-posed problem of extracting the frequency spectrum from the imaginary time correlation functions.

Here we address these two issues based on numerical and analytical calculations of very simple examples, namely a single harmonic oscillator or an ensemble of oscillators with a continuum distribution of frequencies. The interest of this choice is twofold. First, due to its simplicity, we can obtain exact analytical expressions for most quantities of interest, including all time dependent correlations and exact expressions for the discretized path integrals. The availability of these expressions enables a precise control of the different sources of error, which can be both of statistical origin or associated with the discretization itself. Second, the harmonic oscillator is at the heart of the harmonic theory of solids, the natural starting point for any calculation of transport in insulating solids. Completely controlling this case is, therefore, crucial for any serious step forward in this direction. We conclude the study by providing promising indications of the performance of our method when applied to the strongly anharmonic case of a particle evolving in a double well potential.

The manuscript is organized as follows: in Sect. II we introduce the general formalism of the path integral and imaginary time correlations, while in Sect. III we present the procedure that we have developed to cope with the inversion problem. In Sect. IV we next describe a new approach that circumvents the issue of the diverging variance for current-current correlators. Finally, in Sects V and VI we illustrate the application of these methods to a single harmonic oscillator, followed by the case of a collection of oscillators with a continuum distribution of frequencies, mimicking the density of states of a crystalline solid. In Sect. VII we draw our conclusions.

II. THE PATH INTEGRAL FORMALISM FOR TIME CORRELATIONS

The path integral Monte Carlo method provides a numerically exact route to the evaluation of thermodynamic properties of quantum systems at finite temperature, T . If we consider, for simplicity, a system described by a single degree of freedom X of mass m, with Hamiltonian Ĥ = P2 /2m +U( X), the average value of an observable  is

 = 1 Z(β ) Tr[  e -β Ĥ ], (3)
where Z(β ) = Tr[e -β Ĥ ]. In the PIMC approach, the trace is evaluated by expressing the density operator as e -β Ĥ = (e -β Ĥ/M ) M . In the position representation |X , and using the notation ρ(X,Y, τ) = X|e -τ Ĥ |Y , we can write

 = 1 Z(β ) dX 0 . . . dX M X 0 | Â|X 1 ρ(X 1 , X 2 , β /M) . . . ρ(X M , X 0 , β /M). (4)
If an expression for ρ(X,Y, β /M) is known, the observable can be evaluated by sampling the "path" {X 0 . . . X M } with a statistical weight proportional to ρ(X 1 , X 2 , β /M) . . . ρ(X M , X 0 , β /M).

As the matrix element X 0 | Â|X 1 of a local operator  involves in general a term δ (X 0 -X 1 ), the sampling is actually performed over a closed path of M points. In the following we will repeatedly consider the "primitive" approximation, based on the factorization of the kinetic and potential parts of the density operator and valid in the limit of small τ 6 ,

ρ p (X,Y, τ) m 2π h2 τ exp -m (X -Y ) 2 2h 2 τ - τ 2 [U(X) +U(Y )] . (5)
This simplified expression can be replaced by a more accurate one if needed, and if the exact value of ρ is known, as it is the case for the harmonic oscillator, the latter can be used to sample the path more efficiently [START_REF] Feynman | Statistical Mechanics: A Set Of Lectures[END_REF] .

Here, we are interested in equilibrium time correlation functions that determine the linear response properties of the system. A time correlation involving the observables A at time t and B at time t = 0 is the equilibrium average of the product of the operators Â(t) = e itH/h Âe -itH/h , and B(0) = B, which we can write as,

C AB (t/h) = Â(t) B(0) = 1 Z(β ) Tr[ Â(t) B(0)e -β Ĥ ]. (6)
Obviously, the splitting method could be applied to the operators exp(it Ĥ/h). Unfortunately, the statistical weight associated with the resulting path is imaginary, and therefore it is not suitable for usual sampling methods. If, however, the real time t is replaced by an imaginary time t = iτ h, we can write,

C AB (iτ) = 1 Z(β ) Tr[ Âe -τ Ĥ Be -(β -τ) Ĥ ] = 1 Z(β ) dXdX dY dY X| Â|X ρ(X ,Y, τ) Y | B|Y ρ(Y , X, β -τ), (7)
which is defined for 0 ≤ τ ≤ β , and verifies

C AB (iτ) = C BA (i(β -τ)).
Partitioning again the interval [0, β ] into M slices of width ∆τ = β /M, the correlation function can be sampled for discrete values of τ of the form τ k = k∆τ, with k = 0 . . . M -1, at a computational cost that is similar to that needed to calculate the thermodynamic observables of Eq. ( 4), obtaining

C AB (iτ k ) = 1 Z(β ) dXdY dX 1 ...dX M X| Â|X 1 ρ(X 1 , X 2 , ∆τ)... ρ(X k-1 , X k , ∆τ) X k | B|Y ρ(Y, X k+1 , ∆τ)...ρ(X M , X, ∆τ). (8)
As in Eq. ( 4), here the sampling must be performed over the {X 1 . . . X M } coordinates of the path, the X and Y variables being eliminated by the δ -functions contained in the matrix elements of  and B.

III. A STATISTICAL APPROACH TO THE INVERSION PROBLEM

Once the imaginary time correlations, denoted by C(τ) from now on, have been obtained for a set of M discrete values {τ 0 ...τ M-1 } in the interval [0, β ], the real time correlation functions relevant to describe the system physical response can, in principle, be obtained by inverting Eq. (2). This is common to many studies of quantum systems, and generally described as the "analytical continuation" procedure. It is, however, ill-posed, in the sense that if the spectrum S(ω) 21 is described by a set of parameters (such as the values of S on a discrete ω-grid, or the coefficients of an expansion in terms of some basis set), and the C(τ k ) are affected by statistical errors, a very large number of solutions for S compatible with the original data will be found.

This topic is the subject of a vast literature, and it is fair to conclude that no single method emerges as a preferred solution. Generally speaking, most current solutions employ some particular version of a "maximum entropy" approach [START_REF] Jarrell | Bayesian inference and the analytic continuation of imaginary-time quantum monte carlo data[END_REF][START_REF] Boninsegni | Density fluctuations in liquid 4 he. path integrals and maximum entropy[END_REF] . In the context of PIMC and to obtain real time data in combination with real time approximate methods, this procedure was used for instance in Refs. [START_REF] Krilov | Real time quantum correlation functions. II. maximum entropy numerical analytic continuation of path integral monte carlo and centroid molecular dynamics data[END_REF][START_REF] Krilov | Quantum time correlation functions from complex time monte carlo simulations: A maximum entropy approach[END_REF][START_REF] Habershon | Quantum mechanical correlation functions, maximum entropy analytic continuation, and ring polymer molecular dynamics[END_REF] . The spectral function, S ME , is therefore obtained as an average over the possible S(ω)'s (defined by some finite set of parameters), weighted by the probability that they are the exact model given the data set (C, σ 2 ),

S(ω) ME = DS p(S|C, σ 2 )S(ω). (9)
Here DS indicates the phase space element associated with the parametrization of S(ω), C = (C(τ 1 ),C(τ 2 ), . . . ,C(τ M )) † ≡ (C 1 ,C 2 , . . .C M ) † is a line vector that contains the data points, and σ 2 describes the statistical uncertainty of these data in the form of a covariance matrix. By using the Bayes formula, one has,

p(S|C, σ 2 ) = p(C, σ 2 |S) p(C, σ 2 ) p(S), (10)
where p(S) encompasses any prior information on the spectrum.Making the assumption of Gaussian statistics for the likelihood we can write,

p(C|S, σ ) ∝ e -1 2 (C-C[S])(σ 2 ) -1 (C-C[S]) = e -1 2 χ 2 [S] , (11)
which we can interpret as the definition of χ 2 [S]. Here C[S] is the expression of the vector C, obtained by inserting a known spectrum S into the r.h.s. of Eq. ( 2) and computing the resulting M correlation values. In the case of a spectrum defined by the amplitudes A(ω p ) for a set of N ω discrete frequencies on a regular grid, using Eq. ( 2) we obtain,

C[S](τ α ) = N ω ∑ p=1 A(ω p ) e -hω p τ α + e -h(β -τ α )ω p . (12)
In traditional maximum entropy methods, Eq. ( 9) is solved at the saddle point level, by minimizing the functional F =

1 2 χ 2 [S] -H[S].
Here, H[S] is an entropic functional, which assigns a penalty to irregular solutions that would lead to an overfitting of the statistical errors contained in the data. For a positive spectrum, H[S] is usually chosen as the associated Shannon entropy, with a coefficient controlling the strength of the regularisation. In this work we employ the so-called "stochastic analytical inference" or "stochastic maximum entropy" [START_REF] Fuchs | Analytic continuation of quantum monte carlo data by stochastic analytical inference[END_REF] method, where Eq. ( 9) is sampled by Monte-Carlo methods over DS, which can be constrained to positive values of S through the prior probability p(S). (In the case studies described below we have employed a flat prior.)The term [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF] 2 χ 2 [S] can hence be considered as an effective energy functional, and the method can be refined by introducing an additional parameter in the form of an effective inverse temperature Θ as,

S(ω, Θ) ME = Z(Θ) -1 DS S(ω)e -1 2 Θχ 2 [S] . (13)
Here the normalization Z(Θ) = 1/ exp {ΘF(Θ)} is an effective partition function. Note that the traditional maximum entropy approach corresponds to a mean field version of Eq. ( 13), 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 view of the following analysis, we make the simplifying assumption of uncorrelated data points, 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 (τ α ) , (14)
with σ 2 (τ α ) the statistical uncertainty on the data point α. Several arguments [START_REF] Fuchs | Analytic continuation of quantum monte carlo data by stochastic analytical inference[END_REF] have been evoked for fixing Θ = 1. In contrast, in [START_REF] Fuchs | Analytic continuation of quantum monte carlo data by stochastic analytical inference[END_REF] it has been proposed to pick for Θ the value Θ * that maximizes Z(Θ), which is argued to also maximize the posterior probability P(θ |C). This possibility, which corresponds to a balance between energy and entropy dominated solutions, requires however a full free energy calculation. We further 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 our validation procedure is discussed in Sect. V B. At variance with these proposals, we optimize the value of Θ employing the following strategy. An initial data set, C(τ α ), is generated with known statistical uncertainty σ 2 (τ α ) by using, for instance, a path integral simulation of the considered model. In cases were C(τ) is known analytically, synthetic data could also be generated from the exact solution, and introducing a controlled uncertainty. Starting from these data, the spectrum S ME (Θ), described by P degrees of freedom A(ω p ), is obtained through a Monte-Carlo sampling of Eq. ( 13) for a given value of Θ. Note that a well converged Monte-Carlo average will lead to a spectrum S ME (Θ) with an associated χ 2 ∼ O(Mε), where ε is a residual error, while the average χ 2 ∼ O(Mε + P/Θ). We denote CΘ (τ α ) the correlation function associated with this average spectrum.

In order to determine the optimal choice of Θ, therefore discriminating among different models for S(ω) (e.g., different finite discretizations on an ω-grid), 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] . We, therefore, generate P new sets of validation data, C val,i (τ α ) (i = 1, . . . , P ), by using the same technique (even not necessarily with the same accuracy) that we use to produce the original data set, and determine the associated,

χ 2 val = 1 P P ∑ i=1 M-1 ∑ α=0 [ CΘ (τ α ) -C val,i (τ α )] 2 . ( 15
)
Interestingly, this can be interpreted as a measure of the difference between the estimate CΘ (τ α ) and the exact correlation function, denoted by C exact (τ α ). Indeed, by writing

χ 2 val = 1 P P ∑ i=1 M-1 ∑ α=0 [ CΘ (τ α )-C exact (τ α )+C exact (τ α )-C val,i (τ α )] 2 , (16)
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 (τ α ). (17)
Here, the first term is the distance of the estimate to the exact data, while the second is the variance of the validation data leading to a background value χ 2 0 independent of Θ (or any other parameter entering the model description).The choice of Θ will therefore be eventually dictated by the behavior of the first term.

IV. IMPROVED ESTIMATORS FOR CURRENT CORRELATIONS

The computation of transport coefficients typically implies correlation functions involving the momentum operator, a prototypical one being C pp (τ) = p(τ)p(0) . In the path integral approach and within the primitive approximation of Eq. ( 5), the momentum operator is expressed as a difference of coordinates, so that the correlation function for τ = 0 takes the form

C pp (τ k ) = -m 2 h2 ∆τ 2 (x k+1 -x k )(x 1 -x 0 ) , where x k ≡ x(τ k ), and τ k = k∆τ ≡ k β
M is proportional to the discretized imaginary time. The MC evaluation of C pp (τ k ) is hampered by the fact that, when ∆τ gets 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 7 : instead of directly using the above expression for C pp (τ k ), the integrals entering the correlation function can be rearranged obtaining 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 ∆τ, and therefore does not suffer from the diverging variance associated with the "naive" estimator.

We now show that the strategy used to obtain the virial estimator can be generalized to any correlation function involving the momentum operator [START_REF] Carleo | Universal superfluid transition and transport properties of two-dimensional dirty bosons[END_REF] . Specifically, we consider correlation functions of the general form involved in calculation of transport coefficients, e. g., 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 and following this definition one obtains,

C pF (τ k ) = - m 2 h2 ∆τ 2 (x k+1 -x k )F(x k )(x 1 -x 0 )F(x 0 ) + m 2h∆τ (x k+1 -x k )F(x k )F (x 0 ) - + m 2h∆τ (x 1 -x 0 )F(x 0 )F (x k ) - 1 4 F (x k )F (x 0 ) , ( 18
)
This expression is valid for k ≥ 1, while the case k = 0 must be treated separately, along similar lines.

The MC calculation of Eq. ( 18) 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 have generalized the procedure originally used for the kinetic energy calculations (C pp (0)), and obtain a new estimator with reduced variance for general correlation functions. We start from the first term in Eq. ( 18), which has the strongest dependence on ∆τ, and can be expressed as,

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

m 2 h2 ∆τ 2 F(x k )(x k+1 -x k )F(x 0 )(x 1 -x 0 ) = = m 2 h2 ∆τ 2 Z dx 0 dy 0 • • • dy M-1 δ M-1 ∑ i=0 y i F k-1 ∑ i=0 y i + x 0 y k F(x 0 )y 0 ρ 0 (y 0 ; ∆τ) . . . ρ 0 (y M-1 ; ∆τ) exp[-∆τW ], (20)
with

W = M-1 ∑ j=0 V j ∑ i=0 y i + x 0 . (21)
By using the identity:

m h∆τ y k ρ 0 (y k ; ∆τ) = -∂ y k ρ 0 (y k , ∆τ), (22)
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(x k )(x k+1 -x k )F(x 0 )(x 1 -x 0 ) = m h F(x k )(x 1 -x 0 )F(x 0 ) 1 M M-1 ∑ j=1 jV (x j ) - M-1 ∑ j=k+1 V (x j ) - - mk (h∆τM) F (x k )(x 1 -x 0 )F(x 0 ) - m (h∆τM) F(x k )F(x 0 ) . (23)
For the special case F(x) ≡ 1, we can show that Eq. ( 23) reduces to a virial-like formula for the momenta correlations

C pp (τ k ) = x k V (x 0 ) (see App. A).
Repeating the procedure for the terms linear in 1 ∆τ , such as the second term in Eq. ( 23), we can write the correlation in a form that apparently does not depend on ∆τ (recall that M∆τ = β is a constant). The calculations, together with the expressions appropriate for the special case k = 0, are sketched in App. A.

In contrast with the initial expression Eq. ( 18), all terms are now well-defined as ∆τ → 0. We note, however, that the number of terms involved in the first part of Eq. ( 23) increases linearly with M = β /∆τ, so that the gain following our manipulation is not immediately obvious. The argument that Eq. ( 23) 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 |x 1x 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 explicitly verify the variance reduction numerically for the harmonic oscillator in the following section. We conclude this Section by emphasizing that the above derivation to improve generic estimators that involve momentum operators is by no means limited to the harmonic oscillator, but remains valid in general, in particular for the case of interacting particles. Also, note that the derivation of the improved estimator can be adapted beyond the use of the primitive approximation [START_REF] Ceperley | Path integrals in the theory of condensed helium[END_REF] . 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 [START_REF] Tuckerman | Statistical mechanics: theory and molecular simulation[END_REF] , as the variance of the estimator is entirely determined by the analytical form of the kinetic energy part of the action.

V. CASE STUDY I: THE SINGLE HARMONIC OSCILLATOR A. Computing correlation functions

We now apply the methods described above to our test cases. We start by considering the canonical example of a single quantum harmonic oscillator of frequency ω 0 in one dimension, with potential energy V = 1 2 mω 2 0 X 2 , and focus on the time correlation function of an operator with the structure of an energy current, e. g., C pV (τ) = (p(τ)V (τ)) s (p(0)V (0)) s . Note that, while the case of the harmonic oscillator could be considered as oversimplified, this choice of observables already leads to a non trivial structure of the correlation functions. Additional examples involving a continuous distribution of frequencies and a strongly anharmonic system are treated in Sect. VI and App. D, respectively.

The PIMC approach within the primitive approximation allows us to extract the values of the imaginary time correlation function C pV (τ k ), at M discrete time values, τ k = (k -1)β /M. Two main sources of inaccuracy are associated to this procedure: a systematic error, associated with the use of the primitive approximation for the density matrix, and the statistical uncertainty due to finite sampling. In the following we show how to control these issues.

For an 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 PIMC approach (which in this case can be obtained exactly) with the analytical expression for the correlation function C pV (τ), which corresponds to the continuous limit M → ∞. By applying the canonical formalism for the harmonic oscillator, we indeed obtain,

C exact pV (τ) = mh 3 ω 3 0 256 1 sinh 3 (hβ ω 0 /2) × 12 cosh 3hβ ω 0 2 -3hω 0 τ +2 4e -hβ ω 0 + e -2hβ ω 0 + 1 e hβ ω 0 cosh hβ ω 0 2 -hω 0 τ . (24)
In order to calculate the exact expression of the correlation function within the primitive approximation of the discretized path integral, we first note 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. Wick's theorem allows to recast such correlations x 1 . . . x 2n into products of pair correlation functions x i x j which are easily accessible as x i x j = A -1 i j . A is a symmetric M × M matrix, and we can write, x i x j = dXx i x j e -X T AX . We can therefore use numerical methods to calculate the matrix elements, as discussed in App. B. The relative difference between the two calculations is illustrated in Fig. 1. We observe that, for a sufficiently small value of β /M, the deviation is virtually not affected by a change of β .

In addition to this quantitative estimate, it is important to note that, for this system, the discretization preserves the qualitative shape of the correlation functions. One can show (see App. C) that the calculation using a finite but large M corresponds to the exact result (M → ∞) for slightly shifted oscillator strength and inverse temperature. The Trotter error therefore only introduces small quantitative deviations in the spectral density, but does not give rise to spurious qualitative features such as a broadening of the spectral lines. Also, our case studies below are performed by employing an imaginary time discretization τ = β /M = 0.01ω -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 cor- responding 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. 2, by comparing calculations for the energy current correlation function, C pV , using the naive estimator, Eq. ( 18), and the improved version of Eq. (23). The data of Fig. 2 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.

B. The inversion problem

We now use the reconstruction procedure outlined in Sect. III to extract the frequency spectrum for the correlation functions obtained in Sect. V A. In order to perform a reconstruction one needs both to define the set of parameters that expresses the spectral density in Eq. ( 12) and in the integration measure of Eq. ( 13 N ω delta-functions in the ω-space, see Eq. ( 9). 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. ( 24), 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 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 MC 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. 3 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. ( 15), 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. 5 we show χ 2 val as a function of the number of grid points. Clearly, increasing the number of coefficients A(ω i ) of Eq. ( 9) 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. 5 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. ( 24), 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 we show the behaviour of the spectral function for a chosen ω-grid at the indicated values of Θ. As expected from Eq. ( 9), by increasing Θ the result approaches the most probable configuration that describes the correlation function C pV (τ), reducing entropic effects. In Fig. 6 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 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 ≡ χ 2 valχ 2 0 as a function of the effective temperature Θ, at the indicated values of N ω . The value of χ 2 0 corresponds to χ 2 val for a shifted grid with N ω = 10 and Θ = 6. For the purpose of a better representation all values were shifted by 10 -4 . quality of the reconstruction.

We now consider the spectral for C pV (τ) at β = 3, again clarifying the influence of Θ and of the lattice discretization N ω . In Figs. 7 and8 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. 8. The data shown in Fig. 9 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 Θ does not further improve the result. We conclude the above discussion by observing 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 in App. D, 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 Θ.

VI. CASE STUDY II: CONTINUUM DISTRIBUTION OF OSCILLATORS

We now move to our second 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. ( 24) 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 ). (25)
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. Contrary to the single oscillator case, 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. (14).

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. 10, 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. ( 25), with the same variance σ k . The values of χ 2 val , shown in Fig. 12, indicate again a more statistically sound reconstruction corresponding to sparse grids. Unfortunately, none of the curves of Fig. 10, convincingly captures the sharp edges of the exact spectral density, which rather resemble two symmetrically broadened peaks. Considering shifted grids (Fig. 11), 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.

VII. DISCUSSION AND CONCLUSIONS

Here, we have examined the reconstruction of spectral functions for transport coefficients, starting from imaginary time correlation functions obtained by path integral Monte Carlo simulations. In particular, we have described a general strategy for wisely expressing improved estimators with reduced statistical variance for imaginary time correlation functions involving current or momentum operators. We have next introduced an inversion procedure 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 involved parameters, as we have illustrated in the case of the harmonic oscillator spectra employing different values for the effective inverse temperature, Θ, as well as different choices for the grid discretization, N ω , or offset, δ ω. Despite their apparent simplicity, the oscillator models studied here provide challenging benchmarks for the spectral reconstruction due to the sharp undamped deltafunctions they contain.

Pure Bayesian approaches suggest to eliminate the parameters dependence by using a flat prior with 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. In contrast, in our case studies 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.

Indeed, in practice, path integral Monte Carlo data are strongly correlated in imaginary time, undermining 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 (linear or logarithmic grids in regions where spectral densities are flat, for instance) in general strongly modify the results. The representation of a model must therefore be considered itself as a "parameter", making illusory in our view a "parameter-free" Bayesian inversion.

In this paper we have addressed exactly the above difficulties, and developed a validation procedure to quantitatively control any parameter dependence of the Bayesian inversion. Our proposal is based on the quantity χ 2 val constructed from independent data 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 χ 2 val . We have shown explicitly that the new validation step clearly identifies a discrete set of two delta functions in the case study of the single harmonic oscillator, and provides unambiguous indications towards the correct asymmetric sharp edges in the case of an underlying continuous frequency spectrum. In the case of the double well potential discussed in App. D, we demonstrate the utility of sparse frequency grids for describing discrete spectral functions. Furthermore, this example provides a possible recipe for reliable calculation of a general discrete spectrum -using sufficiently sparse grids and varying (non uniformly) the spacing between nodes to achieve an optimal reconstruction validated by the χ 2 val . In all cases, our validation procedure eventually selects models containing just a limited number of parameters, which intrinsically limits the resolution of the reconstruction.

Overall, combining in a consistent workflow Bayesian in-version together with an efficient validation procedure able to select model parameters and effective temperature dependence, indeed seems to offer promising perspectives for capturing qualitative and quantitative features in spectral reconstruction. We also stress that it is straightforward to integrate the proposed validation step in any existing flavor of maximum entropy reconstruction, possibly including additional prior information in the model [START_REF] Jarrell | Bayesian inference and the analytic continuation of imaginary-time quantum monte carlo data[END_REF] . In addition, the validation also provides an objective comparison of the reconstructed spectra with approximated ones, e.g., those obtained from the real-time centroid dynamics [START_REF] Perez | A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals[END_REF] . We conclude by noting that the Green-Kubo method, combined with the harmonic theory of solids and a numerical perturbative treatment of anharmonic effects, has recently proven to be remarkably effective for the determination of heat conductivity at low temperature in systems such as amorphous silicon [START_REF] Isaeva | Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach[END_REF][START_REF] Simoncelli | Unified theory of thermal transport in crystals and glasses[END_REF] . Our hope is to extend those works to arbitrary temperatures and stronger anharmonic effects, on one hand employing path integrals to relax the assumptions underlying the perturbative treatment of anharmonicity, and on the other hand using the strategies for the spectral reconstruction developed in the present paper.

we finally obtain,

m 2 h2 ∆τ 2 (x 1 -x 0 )F(x 0 ) 2 (x 1 -x 0 ) = = ( m h∆τ - m h∆τM ) F(x 0 )F(x 0 ) - m h M-1 ∑ k=0 F(x 0 ) 2 (x 1 -x 0 )V (x k ) + m hM F(x 0 )(x 1 -x 0 )F(x 0 ) M-1 ∑ j=1 jV (x j ) + MV (x 0 ) . (A8)
This expression (a similar argument applies to Eq. ( 23)) reduces to the usual virial formula for the kinetic energy when F = 1. Indeed, using cyclic invariance along the path, one obtains (setting for simplicity h and m to unity):

1 ∆τ 2 (x 1 -x 0 )(x 1 -x 0 ) = = 1 ∆τ - 1 ∆τM + 1 M ( 1 M M-1 ∑ j=0 x j )( M-1 ∑ j=0 V (x j )) -x 0 V (x 0 ) , (A9)
By employing the classical virial theorem for the center of mass of the path, the second and third terms cancel mutually, and we are left with the usual result for the kinetic energy estimator.

Appendix B: Exact computation of correlation for the harmonic oscillator with a discretized path integral

Path integral calculations are usually excessively involved to evaluate them analytically, even if the naive discretized version of the density matrix is used. For the case of a single harmonic oscillator, however, one can find the result explicitly. Within this approximation, the partition function, Z h.o. of the harmonic oscillator is written as a Gaussian integral,

Z h.o. = dXe -1 2 X T AX , (B1)
where X T is a short-hand form for the vector {x 0 , x 1 , ..., x M-1 }, and we have introduced the matrix,

A =     a b 0 . . . b b a b . . . 0 0 b a . . . 0 . . . . . . . . . . . . . . .     (B2)
with a = 2m h2 τ + mω 2 0 τ, and b = -m h2 τ . For a Gaussian weight, the correlation function x(t)x(0) ≡ x i x 0 can obtained as

x i x 0 = (A -1 ) i0 , (B3)
where A -1 is the inverse matrix of A. The formal expression for the partition function is therefore,

Z = (2π) M 2 [detA] -1/2 .
In general, the inversion of A is handled by numerical methods. The diagonal terms, however, can be obtained analytically in a straightforward manner. Indeed, we can first calculates the eigenvalues {λ } of the matrix A by posing, λ x j = ax j + bx j-1 + bx j+1 .

By looking for solutions of the form,

x (k) j = e 2iπ k j M with k ∈ [0, . . . M -1], the corresponding eigenvalues are, λ (k) = a + 2b cos 2πk M ,
and the determinant is det A = ∏ k λ (k) . The terms x 2 i are next obtained directly from

x 2 i = 1 M ∑ i x 2 i = 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 . (B4)
Appendix C: Influence of the primitive approximation for the density matrix on time correlations for the harmonic oscillator

In the case of the harmonic oscillator, the density matrix, ρ(X,Y, τ) = X|e -τ Ĥ |Y can be computed exactly, obtaining

ρ(X,Y, τ) = mω 2π h sinh(hωτ) exp - mω 4h f (X,Y ) , (C1) with f (X,Y ) = (X +Y ) 2 tanh hωτ 2 + (X -Y ) 2 coth hωτ 2 .
(C2) Comparing this expression with Eq. ( 5) specialized to the harmonic oscillator, we realize that the functional dependence on X and Y is identical for both the exact and the approximate expressions. As a consequence, a path integral simulation of an oscillator of frequency ω using the primitive approximation with an imaginary time step τ, will sample the same configurations as an exact calculation with a frequency ω and a step τ , provided that we meet the conditions, For a given value of hωτ and hω, these equations have a unique solution for ω and hω τ . As a consequence, a calculation only involving the configurations sampled by the path (and not the normalisation of the density matrix) will correspond to the exact result for the shifted frequency ω , and inverse temperature β = Mτ , with a relative shift of O(hωτ) 2 .

hω 2 coth hω τ 2 = 1 τ (1 + h2 ω 2 τ 2 4 ), (
Appendix D: Spectral function calculation for a double well potential

In order to illustrate that our discussion is not intrinsically limited to harmonic oscillators, we examine an additional classical benchmark example [START_REF] Perez | A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals[END_REF] with strongly anharmonic features. Let us consider the Hamiltonian describing a particle trapped in a double well potential,

Ĥ = p2 2m -α x2 + 1 2 γ x4 , (D1)
and we choose α = 1 and γ = 1. We can obtain the energy spectrum and the correlation functions of interest by numerical diagonalisation of Ĥ. In the following we consider the position correlation function, C β xx (τ) = x(τ) x(0) . 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 ). (D2)
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 ). We next employ the calculation scheme detailed in Sect. III, and apply the χ 2 val validation criteria in order to determine the optimal reconstruction. As discussed in the main text, we are interested in the effect of the modification of the simulation hyper-parameters on the 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-equallyspaced δ -peaks (see Table I). 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.(D2) 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. 13 we show the resulting spectral functions at the indicated values of N ω and δ ω, while Fig. 14 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 been 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. 15). We perform the reconstruction with the same parameters as in the case above.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. 16). 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 any qualitative effect on the result. Value of Θ, at which the best fit (minimum of χ 2 val ) is achieved, however, becomes smaller with less test data.

2 FIG. 1 .

 21 FIG. 1. 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.

FIG. 2 .

 2 FIG. 2. Difference between the exact correlation function C exactpV (τ) and the values obtained by Monte Carlo 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. IV. We show with line-points the primitive estimator and with the continuous line the improved estimator, both using the same Monte Carlo data.

FIG. 3 .

 3 FIG. 3. 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 deltafunctions 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.

FIG. 4 .

 4 FIG.4. 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.

25 FIG. 5 .

 255 FIG.5. N ω -dependence of ∆χ[START_REF] Ashcroft | Solid State Physics[END_REF] val ≡ χ 2 valχ 2 0 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, at the indicated values of N ω and δ ω. The value of χ 2 0 corresponds to χ 2 val for N ω = 5 and Θ = 50 (see Fig.6, for a better comparison all values were shifted by 10 -5 ).

FIG. 6 .FIG. 7 .

 67 FIG. 6. 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 ≡ χ 2 valχ 2 0 as a function of the effective temperature Θ, at the indicated values of N ω . The value of χ 2 0 corresponds to χ 2 val for N ω = 5 and Θ = 50. For the purpose of a better representation, all values were shifted by 10 -5 .

FIG. 8 .

 8 FIG.8. 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.

FIG. 9 .

 9 FIG. 9. 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 ≡ χ 2 valχ 2 0 as a function of the effective temperature Θ, at the indicated values of N ω . The value of χ 2 0 corresponds to χ 2 val for N ω = 5 and Θ = 40. For the purpose of a better representation all values were shifted by 10 -5 .

FIG. 10 .

 10 FIG. 10.Spectral reconstruction of C cont 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.

25 FIG. 11 .

 2511 FIG. 11.Spectral reconstruction of C cont pV (τ) for the continuous distribution of oscillators, for N ω = 10 and Θ = 1 and 10, 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.

FIG. 12 .

 12 FIG.12. 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 ≡ χ 2 valχ 2 0 as a function of the effective temperature Θ, at the indicated values of N ω . The value of χ 2 0 corresponds to χ 2 val for a shifted grid with N ω = 10 and Θ = 6. For the purpose of a better representation all values were shifted by 10 -4 .

FIG. 13 .

 13 FIG. 13. Spectral reconstruction of C β =8 xx (τ) for the double well potential, at the indicated values of the discretization N ω and Θ. Value of Θ corresponds to the minimum of χ 2 val for N ω = 25 model as shown in the Fig. 14. 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 frequency resolution ∆ω = 0.01.

FIG. 14 .FIG. 15 .

 1415 FIG. 14. Θ-dependence of ∆χ 2 val ≡ χ 2 valχ 2 0 calculated from the validation step of the reconstructed spectral functions for C β =8 xx (τ) for the frequency discretizations N ω = 10 and N ω = 25. The value of χ 2 0 is χ 2 val corresponding to N ω = 25 and Θ = 15. Spectral functions for the non-shifted N ω = 10 grid (shown at Θ = 15 as the blue spectrum in Fig. 13) provides considerably worse result and the corresponding χ 2 val are not shown here.

FIG. 16 .

 16 FIG. 16. Θ-dependence of ∆χ 2 val ≡ χ 2 valχ 2 0 calculated from the validation step of the reconstructed spectral functions for C β =1 xx (τ). The value of χ 2 0 corresponds to a reconstruction with N ω = 10 and a grid shifted by δ ω = 0.25, at Θ = 10. Here we compare shifted (δ ω = 0.25) and non-shifted lattices with the indicated values of N ω .

TABLE I .

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

	ω 1	0.72478
	ω 2	1.25495
	ω 3	1.46154
	ω 4	1.65143
	ω 5	1.80533
							C3)
		and,
		hω tanh	hω τ 2	=	(hω) 2 τ 2	.	(C4)

β xx (τ) for the double well potential.

ACKNOWLEDGMENTS

This work has been supported by the project Heatflow (ANR-18-CE30-0019-01) funded by the french "Agence Nationale de la Recherche". Computations were performed using the Froggy platform of the CIMENT infrastructure, which is supported by the Rhône-Alpes region (Grant No. CPER07-13 CIRA) and the project Equip@Meso (ANR-10-EQPX-29-01) of the ANR. The Authors thank Dr. Victor Hugo Purrello for helpful technical discussions.

Appendix A: Complete expression for virial-like estimators

As we mention in Sect. IV, when calculating currentcurrent correlations one needs to compute products of p and F. As these quantities involve terms proportional to 1/∆τ 2 and 1/∆τ, their variance is quite large. In the same Section, we have demonstrated how to re-express the first one in order to obtain a more accurate MC estimator. Quantities that are linear in 1/∆τ in Eq. ( 18) for C pF (τ k ) can be re-written as,

and,

When computing with Eq. ( 23), 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,

where,

By using the relation

we can re-write (A3) as,

y 0 F(x 0 )ρ 0 (y 0 ; ∆τ) . . . ρ 0 (y M-1 ; ∆τ) exp(-∆τW ). (A6)

By substituting,