# The Schrödinger formulation of the Feynman path centroid density

## I INTRODUCTION

The thermodynamic properties of many quantum systems have been studied within the path integral formulation of statistical mechanics.  An essential aspect of this formulation is the mapping of the quantum system onto a classical model, whose equilibrium properties can be derived by classical molecular dynamics or Monte Carlo simulation techniques with arbitrarily accuracy. However, the efficient calculation of equilibrium dynamical properties of quantum many-body systems by this formulation, i.e., time dependent correlation functions, still remains as a challenging unsolved problem. 

Several applications of the path integral formulation have been developed around a quantity introduced by Feynman and Hibbs: , the effective classical potential (ECP). This concept allows the formulation of a variational principle for quantum systems in thermodynamic equilibrium, whose importance can be appreciated by quoting the conclusions of the Feynman and Hibbs’ book: ”\[This variational principle is\] the only example of a result obtained with path integrals which cannot be obtained in simple manner by more conventional methods”.  The original formulation of this principle used a free particle as a reference system, a fact that limits the range of validity of this approximation to temperatures where the system behaves nearly in a classical way. An essential improvement of this variational theory was formulated by Giachetti and Tognetti , and, independently, by Feynman and Kleinert , using a harmonic oscillator as reference system. With this improvement the equilibrium properties of quantum anharmonic systems can be realistically approximated even in the low temperature limit.  The name ”pure quantum self-consistent harmonic approximation” has been recently coined for the application of the variational principle at this level of approximation. 

Further applications of the ECP concept are related to kinetic and dynamical properties of quantum systems in thermodynamic equilibrium. Gillan formulated the basis of the quantum transition state theory (QTST), a kinetic approximation to calculate rate constants of thermal activated processes, that may be applied to anharmonic many-body systems at temperatures where quantum effects are important.  This approximation has been used in condensed matter and chemical physics investigations, e.g. in recent studies of reorientation rates of hydrogen around a boron atom in doped silicon, , or proton transfer reactions . The most important quantity in QTST is the so called centroid density, a function that carries exactly the same information as the ECP: the centroid density is an exponential function of the ECP. An interesting dynamical approximation, called centroid molecular dynamics (CMD), was formulated by Cao and Voth  to calculate real time correlation functions of quantum particles at finite temperatures. The ECP is an essential quantity for this approximation as the dynamical properties are derived from trajectories generated by classical equations of motion of particles moving in the ECP. A justification of the CMD approach has been given by showing that the centroid position correlation function agrees to second order, in a Taylor expansion in time, with the Kubo transformed position correlation function.  The latter correlation function is related to the dynamical response of the system to an external force. Unfortunately, the limits and capability of CMD for many-body systems at temperatures where quantum effects and anharmonicity are relevant have not been at present fully established. For recent applications of this dynamical approach see Refs.  and .

The centroid density or, equivalently, the ECP, are then important concepts in the theory of path integrals. The motivation of the present work is based on a simple idea: if the centroid density is the basis of several interesting physical applications within the path integral formulation (i.e., a variational principle, and the QTST and CMD approximations), then the definition of this quantity within the Schrödinger formulation may lead to new physical insight into these approximations. Our goal is to present the correspondence between path integral concepts, such as the centroid density or the ECP, and the Schrödinger formulation. In a recent contribution  we have analyzed the centroid density showing that is related to the quasi-static response of the quantum system to a constant external force. Some physical implications of this analysis have been explored in the zero temperature limit, e.g., the ECP corresponds to the mean energy of minimum energy wave packets (MEWP, to be defined below) and CMD is an approximate dynamics based on these wave packets. In this work we present a full account of these findings, with focus on static equilibrium properties.

This paper is organized as follows. The theoretical part is presented in Sec. II, which is divided into several subsections. Sections II A and II B review the Schrödinger and path integral formulation of the equilibrium density matrix. The influence of a constant external force on the Euclidean action of the system is treated in Sec. II C. Important path integral quantities, as the centroid density and the static-force response (SFR) density matrix, are introduced in Sec. II D. In Sec. II E the centroid density and the SFR density matrix are derived in the Schrödinger formulation by means of their moment generating functions. Sec. III presents some relevant physical consequences of the theoretical part. It includes two subsections, Sec. III A deals with the static isothermal susceptibility at finite temperatures, and Sec III B treats the zero temperature limit of the ECP and of the SFR density matrix. The concluding remarks are given in Sec. IV.

## II Theory

### A Schrödinger formulation of the density matrix

We consider the simple case of a quantum particle of mass $`m`$ having bound states in a one-dimensional potential $`V(x)`$. The extension to the many-particle case, for distinguishable particles, is presented in Appendix A. For the subsequent analysis of the Feynman path centroid density, we study the static response of the particle to a constant external force, $`f`$. Therefore, we consider the Hamiltonian of the particle, $`\widehat{H}(f)`$, as a function of the external force:

$$\widehat{H}(f)=\widehat{H}f\widehat{x},$$
(1)   
where $`\widehat{x}`$ is the position operator and $`\widehat{H}\widehat{H}(0)`$ is the Hamiltonian operator in the absence of the force:

$$\widehat{H}=\frac{\widehat{p}^2}{2m}+V(\widehat{x}),$$
(2)   
with $`\widehat{p}`$ being the momentum operator. The following convention is used throughout the text: operators (or functions) that share the same name, as $`\widehat{H}(f)`$ and $`\widehat{H}`$, but differ either by the presence or number of arguments, represent different operators (or functions). The eigenfunctions, $`|\psi _n(f)`$, of the Hamiltonian, $`\widehat{H}(f)`$, satisfy the Schrödinger equation:

$$\widehat{H}(f)|\psi _n(f)=E_n(f)|\psi _n(f),$$
(3)   
where $`E_n(f)`$ are the corresponding eigenvalues. In the position representation, the eigenfunctions are expressed as:

$$\psi _n(x;f)x|\psi _n(f),\text{with}n=0,1,2,\mathrm{}.$$
(4)   
If we consider a canonical ensemble of independent particles in thermal equilibrium at temperature $`T`$, the statistical state of the system in the presence of the external force, $`f`$, is described in terms of the unnormalised density operator:

$$\widehat{\rho }(f)=e^{\beta \widehat{H}(f)},$$
(5)   
where $`\beta `$ is the inverse temperature $`(k_BT)^1`$, with $`k_B`$ being the Boltzmann constant. The position representation of this operator is:

$$\rho (x,x^{};f)=x|\widehat{\rho }(f)|x^{}=\underset{n=0}{\overset{\mathrm{}}{}}e^{\beta E_n(f)}\psi _n(x;f)\psi _n^{}(x^{};f).$$
(6)   
The unnormalised particle’s probability density is given by the diagonal elements of the density matrix, $`\rho (x,x;f)`$. The canonical partition function in the presence of the external force is defined by the trace of the density matrix:

$$Z(f)=_{\mathrm{}}^{\mathrm{}}𝑑x\rho (x,x;f)=\underset{n=0}{\overset{\mathrm{}}{}}e^{\beta E_n(f)},$$
(7)   
where the second equality is obtained by substitution of $`\rho (x,x;f)`$ by the expression given in Eq. (6).

### B Path integral formulation of the density matrix

In the path integral formulation the elements of the unnormalised density matrix in the presence of the external force are given by: 

$$\rho (x,x^{};f)=_{xx(0)}^{x^{}x(\beta \mathrm{})}D[x(u)]e^{\frac{1}{\mathrm{}}S[x(u);f]}.$$
(8)   
$`u`$ is an imaginary time that varies between 0 and $`\beta \mathrm{}`$, and $`S[x(u);f]`$ is the functional of the Euclidean action of the path $`x(u)`$:

$$S[x(u);f]=_0^\beta \mathrm{}𝑑u\left(\frac{m}{2}\dot{x}(u)^2+V[x(u)]fx(u)\right).$$
(9)   
The function $`\dot{x}(u)`$ represents the derivative of $`x(u)`$ with respect to $`u`$. The integral measure of the path integral is given by:

$$D[x(u)]=\underset{N\mathrm{}}{lim}\underset{k=1}{\overset{N1}{}}dx_k\left(\frac{mN}{2\pi \beta \mathrm{}^2}\right)^{\frac{N}{2}},$$
(10)   
where the path $`x(u)`$ has been discretized as $`(x,x_1,x_2,\mathrm{},x_{N1},x^{})`$. 

### C Euclidean action under an external force

It is important to note the effect of a constant external force on the Euclidean action of a given path $`x(u)`$. As the external force is independent of the imaginary time $`u`$, it can be taken out of the time integral in Eq. (9) with the result:

$$S[x(u);f]=S[x(u)]f\beta \mathrm{}X,$$
(11)   
where $`S[x(u)]`$ is the Euclidean action in the absence of the force, and $`X`$ is the average point of the path, or path centroid:

$$X=\frac{1}{\beta \mathrm{}}_0^\beta \mathrm{}𝑑ux(u).$$
(12)   
For the set of paths contributing to the path integral in Eq. (8), the relation of having the same path centroid is an equivalence relation. Then, the sum over paths may be decomposed into a sum over equivalence classes.  A class of paths is the subset formed by all paths that have the same centroid coordinate $`X`$. We note that if two paths, say $`x(u)`$ and $`x^{}(u)`$, belong to the same class, then the contribution of the external force to the Euclidean action, \[i.e., the term $`f\beta \mathrm{}X`$ in Eq. (11)\], is identical for both paths. In other words, the contribution of a constant external force to the Euclidean action is identical for all paths that have the same centroid. This simple property is the origin of the interest of fixed centroid path integrals in statistical physics.

### D Static-force response density matrix and centroid density

The most important quantities to be defined for the class of paths with centroid at $`X`$ are the SFR density matrix  and the centroid density . We have introduced a new name, static-force response density matrix, because we will show that this matrix depends on the response of the system to constant external forces of arbitrary magnitude. The name reduced density matrix has been used before for this quantity.  However, this name is less convenient as it is often encountered in another context. The SFR density matrix is a constrained path integral over a given class of paths, defined by introducing a delta function in the integrand of Eq. (8):

$$\sigma (x,x^{};X;f)=_x^x^{}D[x(u)]\delta \left(X\frac{1}{\beta \mathrm{}}_0^\beta \mathrm{}𝑑ux(u)\right)e^{\frac{1}{\mathrm{}}S[x(u);f]}.$$
(13)   
The centroid density, $`C(X;f)`$, in the presence of the external force $`f`$, is defined as the trace of the SFR density matrix:

$$C(X;f)=_{\mathrm{}}^{\mathrm{}}𝑑x\sigma (x,x;X;f).$$
(14)   
By introducing the expression of the Euclidean action, \[Eq. (11)\], into Eq. (13), we derive an equation that relates the SFR density matrix in the presence of the external force $`f`$, with the same quantity in absence of the force:

$$\sigma (x,x^{};X;f)=\sigma (x,x^{};X)e^{\beta fX}.$$
(15)   
By setting $`x=x^{}`$ in the last equation and integrating over $`x`$, we get the following analogous relation for the centroid density:

$$C(X;f)=C(X)e^{\beta fX}.$$
(16)   
The last two equations show that if the centroid density, $`C(X)`$, and the SFR density matrix, $`\sigma (x,x^{};X)`$, are known for a given class of paths in absence of an external force, then the corresponding quantities in the presence of the force $`f`$ are easily obtained by multiplication with the constant factor $`e^{\beta fX}`$. Dividing the left and right hand sides of the last two equations, we get:

$$[C(X;f)]^1\sigma (x,x^{};X;f)=[C(X)]^1\sigma (x,x^{};X).$$
(17)   
This equation implies that the normalized SFR density matrix defined for a given centroid position $`X`$ is an invariant quantity with respect to the value of the force $`f`$. The normalization constant is the trace of the matrix, i.e., the centroid density. Setting $`x=x^{}`$, we get:

$$\varphi ^2(x;X)[C(X;f)]^1\sigma (x,x;X;f).$$
(18)   
$`\varphi ^2(x;X)`$ represents, as a function of $`x`$, a normalized particle’s probability density. The normalization of this probability density is easily checked by integrating Eq. (18) with respect to the variable $`x`$, and with the help of the definition given in Eq. (14):

$$_{\mathrm{}}^{\mathrm{}}𝑑x\varphi ^2(x;X)=1.$$
(19)   
This relation is valid irrespective of the value of the centroid coordinate, $`X`$.

The results of this Section may help to understand path integral results that are usually presented from a different point of view. For example, the CMD approximation  is an approximate dynamics for a normalized density matrix whose diagonal elements in the position representation are given by $`\varphi ^2(x;X)`$. An important property of this density matrix is its invariance with respect to the application of an external force to the particle. Associated to each centroid position, $`X`$, there is a different density matrix, and the approximate dynamics of CMD is formulated through a dynamical equation for the time evolution of the centroid coordinate.

We are now in position to relate the SFR density matrix and the centroid density to physical observables, an essential step for the understanding of important path integral concepts from the point of view of the Schrödinger formulation.

### E Moment generating function of the centroid density

After the decomposition of the path integral given in Eq. (8) into disjoint classes, we can recover the whole path integral, which is the quantity related to physical observables, by an integral over all classes. The elements of the unnormalised density matrix, in the presence of an external force $`f`$, are then obtained as:

$$\rho (x,x^{};f)=_{\mathrm{}}^{\mathrm{}}𝑑X\sigma (x,x^{};X;f).$$
(20)   
Introducing in this equation the expression given in Eq. (15) for the SFR density matrix, one gets:

$$\rho (x,x^{};f)=_{\mathrm{}}^{\mathrm{}}𝑑X\sigma (x,x^{};X)e^{\beta fX}.$$
(21)   
Setting $`x=x^{}`$ and integrating over the variable $`x`$ \[with the help of Eqs. (7) and (14)\], one arrives at:

$$Z(f)=_{\mathrm{}}^{\mathrm{}}𝑑XC(X)e^{\beta fX}.$$
(22)   
The last two equations provide the essential link for the correspondence between fixed centroid path integrals and the Schrödinger formulation. Although Eq. (22) can be found in Kleinert’s book,  the physical implications of this relation have not been explored until recently.  The physical content of Eqs. (21) and (22) can be derived by two complementary points of view, as suggested by the mathematical structure of these equations:

* The centroid density, $`C(X)`$, in absence of the external force, is related to the partition function, $`Z(f)`$, by an integral transformation defined by the kernel $`e^{\beta fX}`$. The same relation holds between the SFR density matrix, $`\sigma (x,x^{};X)`$, and the density matrix $`\rho (x,x^{};f)`$. This integral transformation is performed with respect to the centroid variable $`X`$, and the product $`\beta f`$ is the corresponding conjugate variable. Note that a general property of an integral transformation is that the physical information contained in the direct and transformed functions is identical. Therefore the centroid density $`C(X)`$, as a function of $`X`$, carries out the same physical information as the partition function $`Z(f)`$, as a function of $`f`$. The same relation holds for $`\sigma (x,x^{};X)`$ and $`\rho (x,x^{};f)`$. This integral transformation is called a two-sided Laplace transform, with properties similar to those of a Fourier transform. 
* The relation between the centroid density, $`C(X)`$, and the partition function in the presence of a force, $`Z(f)`$, can be interpreted as the relation between a probability density for the variable $`X`$ ant its moment generating function. The same relation holds for $`\sigma (x,x^{};X)`$ and $`\rho (x,x^{};f)`$. In the following, we develop this point of view.

The moments of the centroid density in the presence of an external force are defined as:

$$\left\{X^n\right\}_f=\frac{_{\mathrm{}}^{\mathrm{}}𝑑XC(X;f)X^n}{_{\mathrm{}}^{\mathrm{}}𝑑XC(X;f)}=[Z(f)]^1_{\mathrm{}}^{\mathrm{}}𝑑XC(X;f)X^n.$$
(23)   
The brackets $`\left\{\mathrm{}\right\}_f`$ indicate an average over the centroid density, $`C(X;f)`$, and the second equality is obtained from Eq. (22). The moments of $`X`$ in absence of the external force are:

$$\left\{X^n\right\}=Z^1_{\mathrm{}}^{\mathrm{}}𝑑XC(X)X^n,$$
(24)   
where the partition function is $`ZZ(0)`$. Note that the subindex ”$`f`$” is omitted from the brackets when $`f=0`$. The definition of the moment generating function,  $`M(\beta f)`$, of the normalized centroid density, $`Z^1C(X)`$, is:

$$M(\beta f)=\left\{e^{\beta fX}\right\}=Z^1_{\mathrm{}}^{\mathrm{}}𝑑XC(X)e^{\beta fX}.$$
(25)   
From Eqs. (22) and (25) we get:

$$M(\beta f)=\frac{Z(f)}{Z}.$$
(26)   
This result has a clear physical meaning: the ratio of partition functions $`Z(f)/Z`$ is the function generating the moments of the centroid density. Therefore, these moments may be defined as: 

$$\left\{X^n\right\}=\left[\frac{M(\beta f)}{(\beta f)}\right]_{f=0},$$
(27)   
and with the help of Eq. (26), we get:

$$\left\{X^n\right\}=\frac{1}{Z\beta ^n}\left[\frac{^nZ(f)}{f^n}\right]_{f=0}.$$
(28)   
The l.h.s. of the last equation, i.e., the moments of the centroid density, is a quantity typically defined by fixed centroid path integrals, while the r.h.s. contains physical quantities, that are defined within the Schrödinger formulation. It is more interesting to study the cumulant generating function, $`K(\beta f)`$, of the centroid density, which is defined as the logarithm of the moment generating function:

$$K(\beta f)=\mathrm{ln}\frac{Z(f)}{Z}=\beta [F(f)F],$$
(29)   
where $`F(f)`$ is the free energy derived from the partition function: $`Z(f)=\mathrm{exp}[\beta F(f)]`$ and $`FF(0)`$. The cumulants, $`\kappa _n`$, of the centroid density are obtained as:

$$\kappa _n=\left[\frac{K(\beta f)}{(\beta f)}\right]_{f=0}.$$
(30)   
With the help of Eq. (29) we get:

$$\kappa _n=\frac{1}{\beta ^{n1}}\left[\frac{^nF(f)}{f^n}\right]_{f=0}.$$
(31)   
The last equation shows that the cumulants of the centroid density are related to the change in the free energy as a constant external force is acting on the particle. The first cumulants of the centroid density correspond to the mean value and the dispersion of the centroid coordinate:

$$\kappa _1=\left\{X\right\}=\left[\frac{F(f)}{f}\right]_{f=0},$$
(32)   
$$\kappa _2=\delta X^2=\left\{X^2\right\}\left\{X\right\}^2=\frac{1}{\beta }\left[\frac{^2F(f)}{f^2}\right]_{f=0}.$$
(33)   
The centroid dispersion, $`\delta X^2`$, has been related with a ”classical delocalization” of the particle.  The precise physical meaning of this quantity will be presented in Sec. III.

As the mathematical structure of Eqs. (21) and (22) is identical, the derivation done for the moment and cumulant generating functions of the centroid density can be repeated step by step for the SFR density matrix. The most important result of this analysis may be enunciated as follows: the cumulant generating function, $`K_\sigma (\beta f)`$, of the SFR density matrix, $`\sigma (x,x^{};X)`$, is the logarithm of the ratio of the elements of the density matrix in the presence and in absence of the external force $`f`$:

$$K_\sigma (\beta f)=\mathrm{ln}\frac{\rho (x,x^{};f)}{\rho (x,x^{})}.$$
(34)   
As an application of the analysis presented so far we derive in Appendix B the centroid density and the SFR density matrix for a harmonic oscillator by means of their cumulant generating functions. The information needed for this task is the value of the partition function, $`Z(f)`$, and the density matrix, $`\rho (x,x^{};f)`$, as a function of the external force $`f`$. This example shows how results previously obtained by solving fixed centroid path integrals  can be easily derived in the Schrödinger formulation.

To close this theoretical Section, we explain the precise meaning of considering the centroid density as a classical-like density for the quantum particle, which is one of the suggestive pictures derived from the QTST and CMD approximations. Let $`\rho ^{cla}(x)`$ be the classical limit of the unnormalised particle’s probability density in absence of an external force. Then the classical partition function $`Z^{cla}(f)`$, in the presence of a constant external force, is given by:

$$Z^{cla}(f)=_{\mathrm{}}^{\mathrm{}}𝑑x\rho ^{cla}(x)e^{\beta fx}.$$
(35)   
This equation has formally the same structure as Eq. (22) for the centroid density, and illustrates in which sense the centroid density behaves as a classical density for the quantum particle. In the classical limit the ratio $`Z^{cla}(f)/Z^{cla}`$ is the moment generating function of the classical particle’s probability density, while in the quantum case, the ratio $`Z(f)/Z`$ is the moment generating function of the centroid density.

The unnormalised particle’s probability density is given in the classical case by a function of the potential energy:

$$\rho ^{cla}(x)=\left(\frac{m}{2\pi \mathrm{}^2\beta }\right)^{\frac{1}{2}}e^{\beta V(x)},$$
(36)   
In analogy to the classical result, the ECP, $`F_{ecp}(X)`$, for the quantum particle is defined from the centroid density as:

$$C(X)=\left(\frac{m}{2\pi \mathrm{}^2\beta }\right)^{\frac{1}{2}}e^{\beta F_{ecp}(X)}.$$
(37)   
We remark that the ECP is sometimes called the potential of mean force, , the centroid potential,  the effective centroid potential,  or the quantum effective potential.  However, these names refer all to the same physical quantity.

## III Physical implications

In this section we focus on some physical implications readily derived from the correspondence between fixed centroid path integrals and the Schrödinger formulation. Firstly we present some results valid at arbitrary temperatures.

### A Static susceptibility

The first moments of the centroid density have a clear physical meaning. We recall that the average position of the quantum particle in thermal equilibrium in the presence of an external force, $`f`$, is given by:

$$\overline{x}(f)=Z(f)^1_{\mathrm{}}^{\mathrm{}}𝑑xx\rho (x,x;f).$$
(38)   
From the definition of the SFR density matrix and of the centroid density \[Eqs. (13) and (14)\], it is easy to show that the average centroid position, defined from Eq. (23) by setting $`n=1`$, and the average particle position are identical quantities:

$$\overline{x}(f)=\left\{X\right\}_f=\frac{F(f)}{f},$$
(39)   
where the last equality is consequence of an exact generalization of Eq. (32), derived for $`f=0`$, to an arbitrary value of the external force. From Eqs. (33) and (39), we obtain the following expression for the dispersion of the centroid coordinate:

$$\delta X^2=\left\{X^2\right\}\left\{X\right\}^2=\frac{1}{\beta }\left[\frac{\overline{x}(f)}{f}\right]_{f=0}.$$
(40)   
The derivative that appears in the last equation is the static isothermal susceptibility, $`\chi _{xx}^T`$, a quantity that describes the static response of the average position of the quantum particle, with respect to the application of an external force.  We have then:

$$\delta X^2=\frac{1}{\beta }\chi _{xx}^T.$$
(41)   
This isothermal susceptibility, $`\chi _{xx}^T`$, is related to the canonical correlation of the position operator, $`\widehat{x};\widehat{x}`$ by the following relation: 

$$\chi _{xx}^T=\beta \left(\widehat{x};\widehat{x}\overline{x}^2\right),$$
(42)   
where $`\overline{x}\overline{x}(0)`$ and the canonical correlation is defined as: 

$$\widehat{x};\widehat{x}=Z^1\beta ^1_0^\beta 𝑑\lambda \text{Tr}\left[\mathrm{exp}(\beta \widehat{H})\mathrm{exp}(\lambda \widehat{H})\widehat{x}\mathrm{exp}(\lambda \widehat{H})\widehat{x}\right].$$
(43)   
Comparing the Eqs. (41) and (42), and noting that $`\overline{x}=\left\{X\right\}`$, we deduce that the second moment of the centroid density is identical to the canonical correlation of the position operator:

$$\left\{X^2\right\}=\widehat{x};\widehat{x}.$$
(44)   
This result is consistent with the analysis presented in Ref.  concerning the relation between the centroid time correlation function and the Kubo transformed position correlation function. Eq. (44) is interesting not only to understand the physical meaning of the centroid density, but also for practical applications. In the case of a harmonic oscillator of angular frequency $`\omega `$ the dispersion of the centroid coordinate, $`\delta X^2`$, is given in Eq. (B5). The static isothermal susceptibility for the harmonic case is then:

$$\chi _{xx}^T=\beta \delta X^2=\frac{1}{m\omega ^2}.$$
(45)   
For an arbitrary anharmonic potential, we expect that at low temperatures the static isothermal susceptibility will be determined by a frequency close to the first excitation energy, $`\mathrm{\Delta }E`$, of the system. Then, by substitution in the last equation of $`\omega `$ by $`\mathrm{\Delta }E_{app}/\mathrm{}`$, we get the following approximation to $`\mathrm{\Delta }E`$:

$$\mathrm{\Delta }E\mathrm{\Delta }E_{app}=\mathrm{}\left(\frac{k_BT}{m\delta X^2}\right)^{\frac{1}{2}}$$
(46)   
The capability of this approximation has been checked as a function of temperature for several one-dimensional model potentials, that are listed in Table I. The potentials $`V_2`$, $`V_4`$, and $`V_{10}`$ are power functions of $`x`$, whose coefficients were chosen so that a particle of mass 16 au displays the same value of the excitation energy in each of these potentials. $`V_{dw}`$ is a double-well potential where the first excitation energy corresponds to the tunnel splitting. By Monte Carlo path integral simulations  we have obtained the value of the centroid dispersion, $`\delta X^2`$, as a function of temperature, and the excitation energy has been estimated by Eq. (46). The results are shown in Fig. 1, where the value of $`\mathrm{\Delta }E_{app}`$ and the thermal energy, $`k_BT`$, are displayed in units of the $`\mathrm{\Delta }E`$ associated to each potential. For the harmonic potential, $`V_2`$, the approximation is exact at all temperatures. For the potentials $`V_4`$ and $`V_{10}`$ the approximation is remarkably good at temperatures where the thermal energy is lower than about $`1/4`$ of the first excitation energy. In this temperature range the tunnel splitting energy is approximated with an error of about 25 %, a value that at least is of the correct order of magnitude. This approximation can be applied to many-body problems, using the centroid dispersions resulting from the diagonalization of the tensor given in Eq. (A3).

### B $`T0`$ limit

In the zero temperature limit, the ECP and the normalized particle’s probability densities, $`\varphi ^2(x;X)`$, \[defined in Eq. (18)\] have a simple physical meaning: they are related to the eigenvalues and eigenfunctions of the Hamiltonian $`\widehat{H}(f)`$.

We are going to combine two pieces of information to obtain the low temperature limit of $`\varphi ^2(x;X)`$. Firstly, we know that, as $`T0`$, the normalized particle’s probability density derived from the path integral in Eq. (8) converges towards the probability density of the ground state of the Hamiltonian $`\widehat{H}(f)`$. The average value of the position operator $`\widehat{x}`$ for this ground state is $`\overline{x}_0(f)`$. Secondly, we have derived that the dispersion of the centroid coordinate goes to zero in this limit \[see Eq. (40)\]. This implies that the centroid density, $`C(X;f)`$, must be a delta function of $`X`$ centered at its average value $`\left\{X\right\}=\overline{x}_0(f)`$ \[see Eq. (39)\]. Therefore only the class of paths with centroid coordinate at $`X=\overline{x}_0(f)`$ contributes to the path integral in Eq. (8). As a consequence, the normalized probability density, $`\varphi ^2(x;X)`$, derived for the class of paths with centroid at $`X=\overline{x}_0(f)`$ should be identical to the ground state probability density of the particle in the presence of the external force $`f`$:

$$\underset{T0}{lim}\varphi ^2(x;X)=|\psi _0(x;f)|^2,\text{for}X\overline{x}_0(f).$$
(47)   
As illustration of this result we display in Fig. 2 the probability densities, $`\varphi ^2(x;X)`$, obtained by fixed centroid Monte Carlo simulations of a particle in the model potentials $`V_2`$, $`V_4`$, and $`V_{10}`$ at a temperature of $`k_BT=0.001`$ au, which is about 1/300 of the first excitation energy, and therefore a good approximation to the low temperature limit. For each model potential, the simulations were performed at three different centroid positions, $`X`$. The ground state densities $`|\psi _0(x;f)|^2`$ obtained by numerical solution of the time independent Schrödinger equation, for different values of the external force, $`f`$, are also shown in the figure. Both probability densities are identical, apart from tiny deviations (not visible in the scale of the figure) due to the finite temperature and the statistical uncertainty of the simulation. In the harmonic case the displayed probability densities, $`\varphi ^2(x;X)`$, differ only by a rigid displacement as a function of the centroid position $`X`$. We note that these curves are identical to the probability densities of the coherent states of the harmonic oscillator. 

In the following we show that the functions $`\varphi ^2(x;X)`$ are related to a variational principle. Suppose that we look for a quantum state, $`|\mathrm{\Psi }`$, of the particle with Hamiltonian $`\widehat{H}`$ whose mean energy,

$$E_{min}(X)=\mathrm{\Psi }|\widehat{H}|\mathrm{\Psi },$$
(48)   
is minimum against small variations of $`|\mathrm{\Psi }`$. Moreover, this state must satisfy two constraints: its mean position is fixed at an arbitrary value $`X`$, and it is normalized:

$$X=\mathrm{\Psi }|\widehat{x}|\mathrm{\Psi },$$
(49)   
$$1=\mathrm{\Psi }|\mathrm{\Psi }.$$
(50)   
By straightforward application of calculus of variations (see Appendix C), one finds that $`|\mathrm{\Psi }`$ must be the ground state of a Hamiltonian $`\widehat{H}(f)`$, i.e., $`|\mathrm{\Psi }|\psi _0(f)`$:

$$(\widehat{H}f\widehat{x})|\psi _0(f)=E_0(f)|\psi _0(f),$$
(51)   
where the force $`f`$ and the corresponding ground state energy, $`E_0(f)`$, of the Hamiltonian $`\widehat{H}(f)`$, appear as Lagrange multipliers. The value of $`f`$ must be chosen so that the constraint in Eq. (49) is satisfied. The fixed mean position $`X`$ corresponds then to the average position of the ground state $`|\psi _0(f)`$, i.e. $`X\overline{x}_0(f)`$. The average energy, $`E_{min}(X)`$, of these minimum energy states is derived as a function of $`X`$ with the help of Eqs. (51) and (49) as:

$$E_{min}(X)=E_0(f)+fX,\text{for}X\overline{x}_0(f).$$
(52)   
We call the states $`|\psi _0(f)`$ the MEWP’s of the unperturbed Hamiltonian $`\widehat{H}`$.  We can now reinterpretate Eq. (47) by saying that, in the zero temperature limit, the probability density $`\varphi ^2(x;X)`$ corresponds to the MEWP’s of the Hamiltonian $`\widehat{H}`$. A characteristic property of these states is that their average energy is stationary (i.e., minimum) with respect to any arbitrary change in their probability density that leaves constant their average position $`X`$. In the following we show that this average energy is identical to the ECP.

As $`T0`$ we have derived that the centroid density, $`C(X;f)`$, is a delta function centered at the average position of the ground state of $`\widehat{H}(f)`$, i.e. $`\left\{X\right\}_f=\overline{x}_0(f)`$. We also know that the integral of $`C(X;f)`$ with respect to $`X`$ is identical to the partition function $`Z(f)`$ \[see Eqs. (22) and (16)\]. Then, in the zero temperature limit we can write:

$$\underset{T0}{lim}C(X;f)=Z(f)\delta [X\overline{x}_0(f)].$$
(53)   
The asymptotic behavior of $`C(X;f)`$ as $`T0`$ is given by an exponential of the ECP $`\mathrm{exp}[\beta F_{ecp}(X;f)]`$, while the asymptotic behavior of $`Z(f)`$ is given by an exponential of the ground state energy $`\mathrm{exp}[\beta E_0(f)]`$. From Eq. (53) the asymptotic behavior of $`C(X;f)`$ and $`Z(f)`$ should be the same at $`X=\overline{x}_0(f)`$, that is:

$$\underset{T0}{lim}F_{ecp}(X;f)=E_0(f),\text{for}X\overline{x}_0(f).$$
(54)   
From the property given in Eq. (16) for the centroid density, and from the definition of the ECP in Eq. (37) it is easy to derive a relation valid at arbitrary temperature:

$$F_{ecp}(X;f)=F_{ecp}(X)fX.$$
(55)   
By substitution of the last expression into Eq. (54) we obtain the desired result:

$$\underset{\beta \mathrm{}}{lim}F_{ecp}(X)=E_0(f)+fX=E_{min}(X),\text{for}X\overline{x}(f),$$
(56)   
where the last equality corresponds to Eq. (52). Our final conclusion is that, in the limit $`T0`$, the value of the ECP, $`F_{ecp}(X)`$, is equal to the average energy of that MEWP whose average position coincides with the centroid coordinate $`X`$. The MEWP’s turn out to be the ground states of the Hamiltonian $`\widehat{H}(f)`$. The results obtained in this Subsection may be seen as a consequence of the original variational principle of Feynman and Hibbs , when it is applied to the zero temperature limit. However, the important role of the external force does not appear in the original path integral formulation of this variational principle.

The CMD approximation is, in the $`T0`$ limit, an approximate dynamics based on MEWP’s, which is accurate even for highly anharmonic potentials.  Cao and Voth have shown that CMD reproduces correctly the classical limit at high temperatures.  The classical limit of our formulation of CMD in terms of MEWP’s can be derived by the generalization of the results at $`T0`$ to arbitrary temperatures using the MEWP’s given by the excited states $`|\mathrm{\Psi }_n(f)`$ (see Appendix C). From this generalization we find that CMD can be formulated for harmonic systems at arbitrary temperatures as an exact MEWP’s dynamics, even in the classical limit.  The most important applications of CMD are related to condensed phase quantum dynamics of anharmonic systems.  From our analysis in terms of MEWP’s we find that for anharmonic potentials CMD provides accurate results in both the $`T0`$ and the classical limits, but further work is needed to clarify the capability of CMD at intermediate temperatures.

## IV Summary

We have presented the correspondence between fixed centroid path integrals and the Schrödinger formulation. This analysis shows that the Feynman path centroid density and the SFR density matrix depend on the static response of the quantum system to a constant external force. The path centroid density is related by a simple integral transformation to the partition function of the quantum system under the action of an external force. The same integral transformation is found between a classical phase space density and the classical partition function. Therefore, the interpretation of the centroid density as a classical-like density for the quantum system, which is one of the ideas suggested by the QTST and CMD approximations, has a precise physical meaning: the path centroid density behaves ”classically” in the sense that under the action of an external force it has the same static response properties as a classical phase space density. The same integral transformation is found between the normalized canonical density matrix in the presence of an external force and the SFR density matrix. Within the Schrödinger formulation it is essential to introduce a constant external force to define the centroid density and the elements of the SFR density matrix as transformed functions. However, in the path integral formulation the introduction of external forces is not needed, because one works directly with the transformed functions. This facility of the path integral formulation is the origin a large number of path integral applications, from a variational approximation of the thermodynamic properties of quantum system,  to the CMD approximation of time correlation functions of quantum particles in thermal equilibrium. 

In the present work, we have not tried to present the implications of our formulation in current applications based on fixed centroid path integrals. Nevertheless, our analysis has led to results that clarify the physical meaning of fixed centroid quantities. In particular, we have shown that the dispersion of the centroid coordinate is related to the static isothermal susceptibility, and that, in the zero temperature limit, the fixed centroid path integrals are related to the MEWP’s of the unperturbed system. At $`T=0`$, the normalized SFR density matrix is a pure state density matrix corresponding to a MEWP of the Hamiltonian, and the Feynman effective classical potential is the average energy of the MEWP’s.

###### Acknowledgements.

This work was supported by DGICYT (Spain) under contract PB96-0874. We thank E. Artacho for helpful discussions, and L.M. Sesé and M.C. Böhm for critically reading the manuscript.    

## A Extension to many-particles

The results presented for a particle moving in one dimension can be easily generalized to multidimensional $`n`$-body systems of distinguishable particles. We denote the set of $`n`$ particle coordinates as a 3$`n`$-dimensional vector $`𝐫=(𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐},\mathrm{},𝐫_𝐧)`$ where the $`i`$-particle position vector is $`𝐫_𝐢=(r_{ix},r_{iy},r_{iz})`$. The set of external forces is $`𝐟=(𝐟_\mathrm{𝟏},𝐟_\mathrm{𝟐},\mathrm{},𝐟_𝐧)`$. The external force $`𝐟_𝐢`$ acts on the particle $`𝐫_𝐢`$ through the linear term, $`𝐟_𝐢\widehat{𝐫}_𝐢`$, appearing in the potential energy of the Hamiltonian $`\widehat{H}(𝐟)`$. The centroid density associated to the $`n`$-body Hamiltonian in absence of external forces, $`\widehat{H}`$, can be defined by a generalization of Eq. (22) as:

$$Z(𝐟)=_{\mathrm{}}^{\mathrm{}}\mathrm{}_{\mathrm{}}^{\mathrm{}}\mathrm{𝐝𝐫}𝐂(𝐗)𝐞^{\beta \mathrm{𝐟𝐗}},$$
(A1)   
where $`𝐗=(𝐗_\mathrm{𝟏},𝐗_\mathrm{𝟐},\mathrm{},𝐗_𝐧)`$ is a 3$`n`$-dimensional vector formed by the centroid positions of each particle. The last equation represents a 3$`n`$-dimensional two-sided Laplace transform. The cumulants of the centroid density are obtained as the derivatives of the free energy $`F(𝐟)`$ associated to the partition function $`Z(𝐟)`$:

$$\left\{X_{ix}\right\}=\left[\frac{F(𝐟)}{f_{ix}}\right]_{𝐟=\mathrm{𝟎}},$$
(A2)   
$$\left\{X_{ix}X_{jy}\right\}\left\{X_{ix}\right\}\left\{X_{jy}\right\}=\frac{1}{\beta }\left[\frac{^2F(𝐟)}{f_{ix}f_{jy}}\right]_{𝐟=\mathrm{𝟎}}.$$
(A3)   
The last expression represents the components of a tensor that is related to the static isothermal susceptibity tensor by multiplication by the constant $`\beta `$.

## B Centroid density and static-force response density matrix for a linear harmonic oscillator

We consider a canonical ensemble of particles of mass $`m`$ moving in a one-dimensional harmonic potential $`V_{ho}(x)=(1/2)m\omega ^2x^2`$. In absence of an external force, the partition function and the free energy are: 

$$Z=\left[2\mathrm{sinh}\left(\frac{\beta \mathrm{}\omega }{2}\right)\right]^1,$$
(B1)   
$$F=k_BT\mathrm{ln}\left[2\mathrm{sinh}\left(\frac{\beta \mathrm{}\omega }{2}\right)\right].$$
(B2)   
By application of a constant external force, $`f`$, the potential changes to $`V_{ho}(x)fx`$. The new potential energy is also a quadratic function of $`x`$, and the corresponding free energy is easily derived as:

$$F(f)=F\frac{f^2}{2m\omega ^2}.$$
(B3)   
From Eqs. (32) and (33) we obtain the first two cumulants of the centroid density $`C(X)`$ as:

$$\left\{X\right\}=\left[\frac{f}{m\omega ^2}\right]_{f=0}=0,$$
(B4)   
$$\delta X^2=\frac{1}{\beta m\omega ^2}.$$
(B5)   
The higher order cumulants, which are proportional to successive derivatives of the last equation with respect to $`f`$, are zero. This result implies that the centroid density must be a Gaussian function of $`X`$:

$$C(X)=N_1G_X(\overline{X};\delta X^2),$$
(B6)   
where $`N_1`$ is a normalization constant and $`G_X(\overline{X};\delta X^2)`$ is a normalized Gaussian function of $`X`$:

$$G_X(\overline{X};\delta X^2)=\left(\frac{1}{2\pi \delta X^2}\right)^{\frac{1}{2}}\mathrm{exp}\left[\frac{(\overline{X}X)^2}{2\delta X^2}\right].$$
(B7)   
The constant $`N_1`$ can be obtained from Eq. (22) by setting $`f=0`$:

$$Z=_{\mathrm{}}^{\mathrm{}}𝑑XC(X)=N_1$$
(B8)   
The centroid density of the harmonic oscillator is then:

$$C(X)=ZG_X(0;\frac{1}{\beta m\omega ^2}).$$
(B9)   
The last result has been derived several times using fixed centroid path integrals.  The ECP for the harmonic oscillator is obtained from its definition in Eq. (37) and the last equation as:

$$F_{ecp}(X)=Fk_BT\mathrm{ln}\left(\beta \mathrm{}\omega \right)+\frac{1}{2}m\omega ^2X^2.$$
(B10)   
The SFR density matrix $`\sigma (x,x^{};X)`$ may be derived following the same scheme. The cumulant generating function in this case was given in Eq. (34). The elements of the density matrix for the harmonic oscillator in absence of an external force are: 

$$\rho (x,x^{})=\left[\frac{m\omega }{2\pi \mathrm{}\mathrm{sinh}(\beta \mathrm{}\omega )}\right]^{\frac{1}{2}}\mathrm{exp}\left\{\frac{m\omega }{2\mathrm{}\mathrm{sinh}(\beta \mathrm{}\omega )}\left[(x^2+x^2)\mathrm{cosh}(\beta \mathrm{}\omega )2xx^{}\right]\right\}.$$
(B11)   
While in the presence of an external force we obtain from the definition in Eq. (6):

$$\rho (x,x^{};f)=\mathrm{exp}\left(\frac{\beta f^2}{2m\omega ^2}\right)\rho (xx_{min},x^{}x_{min}),$$
(B12)   
where we have used the following relations:

$$E_n(f)=E_n\frac{f^2}{2m\omega ^2},$$
(B13)   
$$\psi _n(x;f)=\psi _n(xx_{min}).$$
(B14)   
$`x_{min}=f/(m\omega ^2)`$ is the position of minimum energy for the potential $`V_{ho}(x)fx`$. The cumulants of the centroid variable $`X`$, with $`\sigma (x,x^{};X)`$ as its probability density, are evaluated as the derivatives of the cumulant generating function \[Eq. (34)\] with respect to the force at $`f=0`$. We find that only the first and second cumulants $`(\kappa _1,\kappa _2)`$ are different from zero. Therefore $`\sigma (x,x^{};X)`$ must be a Gaussian function of the variable X, with $`\kappa _1`$ being the mean value of $`X`$, and $`\kappa _2`$ its dispersion:

$$\sigma (x,x^{};X)=N_2G_X(\kappa _1;\kappa _2).$$
(B15)   
The result for the first cumulant is obtained from the cumulant generating function after straightforward algebra as:

$$\kappa _1=\frac{\delta x_{cla}^2}{\delta x^2}\left(\frac{x+x^{}}{2}\right),$$
(B16)   
and the result for the second cumulant is:

$$\kappa _2=\frac{\delta x_{cla}^2}{\delta x^2}(\delta x^2\delta x_{cla}^2).$$
(B17)   
The constants $`\delta x^2`$ and $`\delta x_{cla}^2`$ correspond to the dispersion of the position coordinate for a canonical ensemble of harmonic oscillators in the quantum and classical case, respectively:

$$\delta x^2=\frac{\mathrm{}}{2m\omega }\mathrm{coth}\left(\frac{\beta \mathrm{}\omega }{2}\right),$$
(B18)   
$$\delta x_{cla}^2=\frac{1}{\beta m\omega ^2}.$$
(B19)   
The normalization constant, $`N_2`$, is determined from Eq. (21) by setting $`f=0`$. We get:

$$N_2=\rho (x,x^{}).$$
(B20)   
The final result is:

$$\sigma (x,x^{};X)=\rho (x,x^{})G_X[\frac{\delta x_{cla}^2}{\delta x^2}\left(\frac{x^{}+x}{2}\right);\frac{\delta x_{cla}^2}{\delta x^2}(\delta x^2\delta x_{cla}^2)].$$
(B21)   
The elements of the SFR density matrix are the product of the elements of the canonical density matrix by a normalized Gaussian function of the centroid coordinate $`X`$. This equation is identical to Eq. (B.4) of Ref.  (apart from a different grouping of factors).

## C Minimum energy wave packets

We want to find the MEWP’s associated to the Hamiltonian $`\widehat{H}`$. To simplify the notation we use the following abbreviations: $`\mathrm{\Psi }x|\mathrm{\Psi }`$, $`\mathrm{\Psi }^{}\mathrm{\Psi }/x`$, $`\mathrm{\Psi }^{\prime \prime }^2\mathrm{\Psi }/x^2`$, $`VV(x)`$, and $`D\mathrm{}^2/(2m)`$.

The MEWP’s minimize the functional:

$$E_{min}=_{\mathrm{}}^{\mathrm{}}𝑑x(D\mathrm{\Psi }\mathrm{\Psi }^{\prime \prime }+V\mathrm{\Psi }^2)$$
(C1)   
with the constraints given in Eqs. (49,50) for the wave function $`\mathrm{\Psi }`$. We call:

$$A=D\mathrm{\Psi }\mathrm{\Psi }^{\prime \prime }+V\mathrm{\Psi }^2fx\mathrm{\Psi }^2E\mathrm{\Psi }^2,$$
(C2)   
where $`f`$ and $`E`$ are two Lagrange multipliers. From the Euler equation:

$$\frac{A}{\mathrm{\Psi }}\frac{d}{dx}\frac{A}{\mathrm{\Psi }^{}}+\frac{d^2}{dx^2}\frac{A}{\mathrm{\Psi }^{\prime \prime }}=0,$$
(C3)   
one derives the differential equation that must be satisfied by the MEWP’s:

$$D\mathrm{\Psi }^{\prime \prime }+(Vfx)\mathrm{\Psi }=E\mathrm{\Psi },$$
(C4)   
which is the time independent Hamilton equation corresponding to the Hamiltonian $`\widehat{H}(f)`$. All functions satisfying this equation, i.e., the eigenfunction $`|\mathrm{\Psi }_n(f)`$, are MEWP’s in the sense that they minimize the functional in Eq. (C1) subject to the constrains given in Eqs. (49,50). The excited states ($`n>1`$) correspond to local minima of the functional. In the study of the $`T0`$ limit of the SFR density matrix we use only the states $`|\mathrm{\Psi }_0(f)`$ corresponding to the global minimum of the functional, but the excited states are necessary to generalize this study to arbitrary temperatures.