Computational modeling of the condensed phase based on classical statistical mechanics has been rapidly developing over the last few decades and has yielded important information on various systems containing up to millions of atoms. However, if a system of interest contains important quantum effects, well-developed classical techniques cannot be used. One way of treating finite temperature quantum systems at equilibrium has been based on Feynman’s imaginary time path integral approach and the ensuing quantum-classical isomorphism. This isomorphism is exact only in the limit of infinitely many classical quasiparticles representing each physical quantum particle. In this work, we present a reductionist perspective on this problem based on the emerging methodology of coarse-graining. This perspective allows for the representations of one quantum particle with only two classical-like quasiparticles and their conjugate momenta. One of these coupled quasiparticles is the centroid particle of the quantum path integral quasiparticle distribution. Only this quasiparticle feels the potential energy function. The other quasiparticle directly provides the observable averages of quantum mechanical operators. The theory offers a simplified perspective on quantum statistical mechanics, revealing its most reductionist connection to classical statistical physics. By doing so, it can facilitate a simpler representation of certain quantum effects in complex molecular environments.

## I. INTRODUCTION

Quantum statistical mechanics^{1} relates macroscopic properties of quantum systems to a microscopic description of particles constituting such systems. It is vital for the correct description of hydrogen tunneling and nuclear quantum effects,^{2–11} electron transfer,^{3,12–14} infrared and Raman spectra,^{5,15,16} etc., in addition to the paradigmatic examples of condensed phase helium and hydrogen.^{17–21} In comparison to classical statistical mechanics, the quantum theory is arguably more difficult to understand, since it involves particle delocalization, therefore effectively more degrees of freedom, and uses a more advanced mathematical apparatus. From the viewpoint of numerical simulations, quantum statistical mechanics is also more complex than its classical counterpart. For these reasons, classical-like approximations of the exact quantum statistical theory, providing a reductionist picture and facilitating conceptual modeling and numerical simulation of quantum systems, are important.

Feynman’s concept of the imaginary time path integral, with the ensuing quantum-classical isomorphism, can serve as a basis for such reductionist formulations.^{1,22} By replacing each quantum particle with infinitely many classical quasiparticles via a discretization of the path integral, this theory expresses any observable property of a quantum system at equilibrium and at finite temperature in terms of certain properties of an isomorphic classical system.^{23–28} A number of computational approaches based on combining the discretized path integral representation with molecular dynamics (MD), Monte Carlo (MC), or other classical algorithms have been developed.^{18,25–29} These methods would provide numerically exact quantum mechanical results in the limit of an infinite number of classical quasiparticles per quantum particle, the latter being usually denoted by the integer *P*.

In this paper, the most essential, minimalist picture of the manifestation of quantum effects in systems at finite temperature and at equilibrium is presented, based on a coarse-grained (CG) representation of the intermediate isomorphic quasiparticles into a single CG quasiparticle. The CG methodology is commonly used for dealing with classical atomistic models of molecules (see, e.g., Ref. 30) and, to the best of our knowledge, has not yet been applied to quantum problems, except in the case of the imaginary time Feynman path integral “centroid” variable.^{22,31} The resulting coarse-grained path integral (CG-PI) theory provides a new reductionist perspective on quantum statistical mechanics, allowing for a simple and intuitive representation of quantum systems at equilibrium.

## II. THEORY

Starting from a standard imaginary time path integral representation of the thermal density matrix in the position basis (Fig. 1(a)), we first demonstrate that the thermal density matrix can be written as

where $ H \u02c6 $ is the Hamiltonian of the system, *m* is the mass of the quantum particle, *β* is 1/*k _{B}T*,

*Q*is the coordinate of the effective particle introduced by coarse-graining the (infinitely) many intermediate isomorphic classical quasiparticles, and

*V*is the effective CG potential (Fig. 1(c)). It is to be noted that in this paper, the word “quasiparticle” refers only to the classical-like particles that are utilized in the discretized imaginary time Feynman path integral. We are not referring to other kinds of quasiparticles, such as excitations in a Fermi liquid. For notational convenience, Eq. (1) is written for a one-dimensional case; generalizations to more dimensions are straightforward, as shown in Sec. S7 in the supplementary material.

_{CG}^{32}It is informative to compare Eq. (1) to the following path integral approximation for the density matrix with only one intermediate imaginary time slice (Fig. 1(b)):

where *V*_{2} is Feynman’s isomorphic potential^{22} *V _{P}* at

*P*= 2,

where *V*(*q*) is the potential in which the quantum particle is placed. Thus, the proposed CG approach can be considered as a renormalization of Feynman’s isomorphic potential for *P* = 2, leading to an, in principle, exact expression for the density matrix, Eq. (1). The formal definition of the effective potential *V _{CG}* is given below.

This approach is based on the extension of the CG methodology to the field of quantum statistical mechanics. In the original multiscale coarse-graining (MS-CG) theory,^{33,34} coarse-graining is performed for systems formed by *n* interacting particles obeying the laws of classical mechanics, for example, for (bio)molecules described by a set of classical particles (“atoms”) governed by classical force fields. A projection from the space of the coordinates of all atoms in the system **r**^{n} onto the space of *N* CG variables **R**^{N}, where *N* ≪ *n*, is defined by means of *N* linear mapping operators $ M R N = M R 1 r n , \u2026 , M R N r n $,

where *c _{Ii}* are coefficients. One possible interpretation of Eq. (4) is that the system of interest is divided into

*N*domains, and CG variables

**R**

^{N}designate positions of the centers of mass of each domain.

^{33}Other interpretations are also possible, for example, when only some atoms (e.g., C

_{α}in each amino acid) are used as representative atoms, or when the CG variables are defined in terms of vibrational normal modes.

^{35}The CG potential $U R N $ in the MS-CG theory is defined as

^{33}

where *u*(**r**^{n}) is the atomistic potential.

The extension of this CG methodology to quantum systems is possible due to the quantum-classical isomorphism. Coarse-graining the intermediate beads or quasiparticles that correspond to the intermediate states in a path integral (Fig. 1(a)) allows not only for a formally exact description of a quantum system with classical-like quasiparticles but also for a major reduction in the number of quasiparticles used. The first and the last quasiparticles, however, should not be coarse-grained, since the resulting expression for the thermal density matrix must keep an explicit dependence on *q* and *q*′. This background motivates the definition of the CG variable in the present theory of CG-PI as the center of mass of all intermediate quasiparticles in the path integral,

In the limit of *P* → ∞ and if *q* = *q*′, the CG variable *Q* coincides with the centroid defined as^{22,31}

where $q \tau $ is a closed imaginary time Feynman path. Here, $q \tau $, with *τ* running from 0 to *ħ β*, is a continuous analogue of *q _{i}*, with

*i*running from 1 to

*P*+ 1, in the limit of

*P*→ ∞.

The effective potential *V _{CG}* is defined here as

where

The proposed definition, Eqs. (8) and (9), is in general similar to the definition of the CG potential in the MS-CG theory, Eq. (5), but differs from the latter in two aspects. First, integration over variables *q* and *q*′ is not performed, and the resulting CG potential *V _{CG}* depends on a set of mixed resolution variables, namely, fine-grained (

*q*,

*q*′) and coarse-grained (

*Q*) ones. Second, Eqs. (8) and (9) do not include an unspecified constant; therefore, in principle, exact values of the partition function and its correct dependence on temperature can be obtained.

In the case of a free particle, the CG potential assumes the following form:

where

and

^{32}

The first term on the right-hand side of Eq. (11) penalizes large deviations of the CG particle from the average position between the initial and final points, while the second term penalizes large distances between the initial and final points. In the classical or high-temperature limit, the initial, final, and CG coordinates all collapse into the same point in space, while in the low-temperature limit, corresponding to strong quantum effects, the particles tend to be equidistributed over the whole space. A number of other analytical results for the CG potential are provided in the Appendix and the supplementary material.^{32}

The CG potential can be used for computing the thermal density matrix, as well as for the analysis of quantum delocalization in the system of interest at a given temperature. In particular, maximum contributions to the thermal density matrix, according to Eq. (1), come from positions of the CG particle *Q* close to the value *Q _{min}* defined as follows:

For example, for the harmonic oscillator with the potential energy $V q =m \Omega 2 q 2 /2$, $ Q min = tanh f / f q \u0304 $, where *f* = *ħ β* Ω/2. In the classical limit, this $ Q min \u2192 q \u0304 $, while in the limit of strong quantum effects, *Q _{min}* → 0 regardless of the $ q \u0304 $ value. This analysis can be generalized to the case of an arbitrary molecular system. For any two geometrical configurations of the system of interest, one may compute the multi-dimensional

*Q*and interpret it as the most important geometry of the system as it propagates from the first given geometrical configuration to the second one.

_{min}Series expansions of *V _{CG}* in terms of ℏ as a small parameter can be used to study quantum systems near the classical limit. As demonstrated in Sec. S1 in the supplementary material,

^{32}the lowest order nontrivial terms in

*V*have the following form:

_{CG}where

where $ R Q q \u0304 $ and *R*_{Δq2} are renormalization factors,

where *f _{eff}* is defined as

Equation (15) demonstrates that the most essential part of the three-body interactions between the initial, final, and CG intermediate beads reduces in the limit of *ħ* → 0 to a sum of three simple less-than-three-body terms.

The CG particle feels the potential

*V*in a classical way. This is a one-body potential.The CG particle is bound by a harmonic spring to the middle point between the initial and final beads. Since this middle point is, in general, “empty” which is not coinciding with

*q*or*q*′, this term can be considered as a quasi-two-body potential. In the classical limit, the spring stiffness $ k Q q \u0304 $ reduces to the corresponding spring constant for a free particle, while in the limit of strong quantum effects, it tends to $ k Q q \u0304 \u2009 free particle \u22c5 f eff ( q \u0304 ) /3$.The initial and final quasiparticles attract to each other via a two-body effective harmonic potential. In the classical limit, the spring stiffness

*k*_{Δq2}is the same as that of a free particle, while in the regime of strong quantum effects, it tends to $ k \Delta q 2 \u2009 free particle \u22c5 f eff ( q \u0304 ) $.

Two general strategies of improving the accuracy of Eq. (15) are possible. First, one may analytically derive explicit expressions for higher order terms [see Eq. (A15) and Secs. S1.2-1.4 in the supplementary material^{32}]. Second, one may choose an *Ansatz* expression for the CG potential, for example,

and find the exact values of the renormalization factors $ R Q , R Q q \u0304 , R \Delta q 2 , R 0 $ using one of the variational methods described in the Appendix (see also Sec. S4 in the supplementary material^{32}). Both strategies (the use of higher order terms or an *Ansatz* expression) can employ the potential energy *V*(*q*) obtained from electronic structure methods, for example *ab initio* molecular dynamics (AIMD), or approximated by empirical force fields.

It is valuable to further manipulate the preceding expressions to cast them in a more classical-like form. An expectation value of a quantum mechanical operator $ A \u02c6 p \u02c6 , q \u02c6 $ at equilibrium at finite temperature is given by

and can be computed by expressing the thermal density matrix in terms of the CG potential, Eq. (1), and the matrix elements of the operator $ A \u02c6 p \u02c6 , q \u02c6 $ in terms of the corresponding Weyl map^{36} *W _{A}*, such that

(Recall that the Weyl map is a strict way to get a classical-like function corresponding to a given quantum operator. For example, if $ A \u02c6 = q \u02c6 n $, where *n* is an integer, then *W _{A}* =

*q*, and if $ A \u02c6 = p \u02c6 n $, then

^{n}*W*=

_{A}*p*.) The final result (see Sec. S5.1 in the supplementary material for a derivation

^{n}^{32}) is given by

where $ H eff P Q , p , Q , q \u0304 $ has the form of a classical two-particle Hamiltonian. The exact expression for *H _{eff}* is given by

where *M _{Q}* and

*C*are arbitrary positive constants. Based on Eq. (15), the following approximation for

*H*can be obtained:

_{eff} where *m _{eff}* is a positive-valued function of $ q \u0304 $ interpreted as an effective mass. For the derivation of Eq. (26) and explicit expressions for

*m*, see Sec. S5.2 in the supplementary material.

_{eff}^{32}In particular, when $ V \u2033 ( q \u0304 ) >\u2212 \pi 2 m \u0127 \u2212 2 \beta \u2212 2 $, the effective mass is given by

*m*=

_{eff}*m*⋅

*R*

_{Δq2}with

*R*

_{Δq2}defined by Eq. (19).

A simple classical-like picture for the most essential part of the quantum effects emerges then from Eqs. (24) and (26). A physical quantum particle is represented here by two classical quasiparticles that we will call a centroid particle and an observed particle (Fig. 2). The centroid particle, but not the observed particle, feels the physical potential. On the other hand, only the coordinates and momenta of the observed particle, and not those of the centroid particle, are used to compute the values of the function $ W A p , q \u0304 $ to be averaged in order to calculate the observable $ A \u02c6 p \u02c6 , q \u02c6 $. The centroid and observed particles are connected by an effective harmonic spring. An algorithmic implementation of this approach seems straightforward based on the existing codes for classical MD or MC simulations and should not add significantly to the cost of simulations of molecular and biomolecular systems, especially if only some of the atoms in the system (e.g., hydrogen) exhibit quantum behavior.

In the present scheme, the spring stiffness $ k Q q \u0304 $ depends on the current position of the observable particle $ q \u0304 $. Note that $ k Q q \u0304 $ is well-defined (real) not only for positive, but also for negative $ V \u2033 ( q \u0304 ) $. In the latter case, the value of $ f eff ( q \u0304 ) $ is purely imaginary, and the right hand side of Eq. (18) is real-valued. The negative values of $ k Q q \u0304 $ achieved when $ V \u2033 ( q \u0304 ) <\u22124 \pi 2 m \u0127 \u2212 2 \beta \u2212 2 $ [e.g., near the top of a sufficiently high and narrow barrier in the case of a double-well potential *V*(*q*)] correspond to an effective repulsion between the centroid and observed quasiparticles due to the term $ k Q q \u0304 q \u0304 Q \u2212 q \u0304 2 /2$ in Eq. (26). In this case, the two quasiparticles (those with the coordinates *Q* and $ q \u0304 $) tend to move away from each other, falling from the top of the barrier toward two different local minima on the potential energy surface. However, as soon as the observed quasiparticle leaves the vicinity of the top of the barrier, implying that $ V \u2033 ( q \u0304 ) $ increases to − 4*π*^{2}*mħ*^{−2}*β*^{−2} and above, the spring stiffness becomes positive again, and the two quasiparticles start attracting each other, tending to return to the vicinity of the same local minimum. Events when the observed and centroid quasiparticles in the CG-PI model sporadically appear in basins of two different local minima on the potential energy surface reflect the ability of a quantum particle to tunnel. Therefore, possible negative values of $ k Q q \u0304 $ do not contradict to the physical sense and should not lead to numerical divergences in MD or MC simulations with the use of the effective potential *H _{eff}* given by Eq. (26).

The effective mass of the observed quasiparticle *m _{eff}* in Eq. (26) also depends on $ q \u0304 $. In this case,

*m*is real and positive for any value of $ q \u0304 $, as discussed in Sec. S5.2 in the supplementary material.

_{eff}^{32}This definition of the effective mass does not seem to lead to any technical complications in MC simulations. However, in the case of MD, the varying effective mass causes the dependence of the force on the momentum, making common MD integrators, such as the leapfrog algorithm, inapplicable. The problem of how to reliably and efficiently run MD simulations for particles with changing masses in a thermostat goes beyond the scope of the present paper. To deal with this issue, we introduce here one more approximation for the effective Hamiltonian

*H*, namely,

_{eff} where $ m eff const $ is a positive constant. Importantly, for momentum-independent quantum mechanical operators $ A \u02c6 = A \u02c6 q \u02c6 $, the error introduced by the approximation of going from Eqs. (25) to (26) and the error of going from Eqs. (26) to (27) appear to exactly cancel each other (for the proof, see Sec. S5.4 in the supplementary material^{32}). Thus, for a wide class of practically important observables, including average potential energies, average interatomic distances and functions thereof, radial distribution functions, etc., the CG-PI model can employ quasiparticles with constant masses without introducing any errors beyond those coming from approximating the CG potential by Eq. (15) or (21) or (A15).

Equation (24) with *H _{eff}* with constant

*m*and $ k Q q \u0304 $ is exact in the case of a quantum harmonic oscillator for arbitrarily strong quantum effects (that is, to all powers of ℏ) and for arbitrary operators $ A \u02c6 p \u02c6 , q \u02c6 $. Renormalization of $ k Q q \u0304 $—and, for momentum-dependent observables, of

_{eff}*m*—accounts for the major part of quantum corrections in the case of a general potential

_{eff}*V*(

*q*). The remaining corrections follow from higher order terms in the CG potential, Eq. (15), resulting in additional corrections to the spring stiffness $ k Q q \u0304 $, the effective mass

*m*and additional terms in

_{eff}*H*anharmonically coupling

_{eff}*Q*and $ q \u0304 $, and higher order powers of

*p*. In particular, the $O \beta V ( I V ) q \u0304 \Delta q 4 $ term in the CG potential [see Eq. (S79) in the supplementary material

^{32}] results in relativistic-like $O p 4 $ and higher order in

*p*terms in

*H*. If necessary, these terms can be analytically derived in a systematic manner or accounted for by means of a variational treatment with the

_{eff}*Ansatz*representations of the CG potential, as in Eq. (21). The second approach leads to the following expression for the effective Hamiltonian

*H*:

_{eff} where the renormalization factors *R*_{Δq2}, *R _{Q}*, and $ R Q q \u0304 $ are the same as in Eq. (21).

Finally, note that the CG-PI formulas discussed in this section can be generalized to multidimensional cases (see Sec. S7 in the supplementary material^{32}). In particular, the three-dimensional generalization of Eq. (15) has the following form:

where $q, q \u2032 ,Q, q \u0304 $, and Δ**q** are the corresponding three-dimensional vectors, $ K Q q \u0304 $ and **K**_{Δq2} are 3 × 3 matrices given by

and $ \Lambda Q q \u0304 ( q \u0304 ) , \Lambda \Delta q 2 ( q \u0304 ) $ and $ \Lambda H q \u0304 $ are diagonal 3 × 3 matrices with the elements obeying the following conditions:

with *k*_{0} given by

The multidimensional analogue of the expression for the effective Hamiltonian, Eq. (26), assumes the following form:

where the effective mass matrix $ M eff =mU ( q \u0304 ) \Lambda \Delta q 2 ( q \u0304 ) U ( q \u0304 ) \u2020 $.

## III. NUMERICAL EXAMPLES

To illustrate the behavior of the CG-PI method, we consider here two model problems that represent, though in a simplified form, typical characteristics of physical potentials in realistic molecular systems. The CG-PI results are based on Eq. (27) with the physical particle mass used for $ m eff const $. The first model problem is intended to imitate nuclear quantum effects in a condensed phase at low temperatures. In this problem, a particle is placed in an asymmetric anharmonic one-dimensional potential

where the coordinate *q* is in angströms and *V* is in kcal/mol. The values of *a*, *b*, and *c* are 0.252, −0.570, and 0.814, respectively. This potential *V*(*q*) is essentially anharmonic, as demonstrated in Fig. 3(a) by a comparison to the corresponding harmonic potential

The temperature is taken to be 5 K, and the mass of the particle is 20.0 g/mol. (For more details on the motivation of this problem and the numerical techniques, see Sec. S6.1 in the supplementary material.^{32})

The CG-PI approach—more specifically, the simplest version of it, with quadratic coupling between the centroid and observed quasiparticles, and in the assumption of constant masses of both quasiparticles, as in Eq. (27)—yields in this case an asymmetric non-Gaussian shape of the thermal density matrix as a function of the coordinate *q* (Fig. 3(b), black dotted line). This curve is much closer to the curve for the exact thermal density matrix (Fig. 3(b), thick blue line) than to the probability density for the same potential *V*(*q*) in classical statistical mechanics (thin red line), suggesting that the simplest version of the CG-PI method nearly completely accounts for the quantum effects in this model system. In particular, the CG-PI method precisely reflects broadening of the particle distribution in the quantum thermodynamic ensemble due to tunneling of the physical particle. It is important to note that the CG-PI result for this system is closer to the exact thermal density matrix for the potential *V*(*q*) than to the exact thermal density matrix for the potential *V _{HO}*(

*q*) (blue dashed line) over the most part of the

*q*axis, except for

*q*> 0.5 Å, confirming the statement made in Sec. II that a renormalization of the effective interaction between the centroid and observed quasiparticles allows for the major part of quantum corrections to anharmonic contributions to a thermal density matrix.

The second model problem is designed to imitate tunneling of a hydrogen atom through a potential energy barrier at temperatures typical of living organisms. A particle is placed in an asymmetric double-well one-dimensional potential

where *q* is in angströms and *V* is in kcal/mol (Fig. 4(a)). The constants *c*_{0}, *c*_{1}, *c*_{2}, *c*_{3}, and *c*_{4} are 6, −0.5, −44, 2, and 88, respectively. The temperature is 310 K and the particle mass is 1.008 g/mol. (For more details, see Sec. S6.2 in the supplementary material.^{32})

Under these circumstances, the particle exhibits quantum behavior, despite the high temperature and the presence of a reasonably high barrier (5 kcal/mol) to tunnel through. The diagonal values of the thermal density matrix predicted with the CG-PI method (Fig. 4(b), black dotted line) are much closer to the values of the exact thermal density matrix (Fig. 4(b), thick blue line) than to the probability density for the same potential *V*(*q*) in classical statistical mechanics (thin red line). Broadening of the particle distribution, which is particularly noticeable on the sides of the curves closer to the potential energy barrier, is quantitatively accounted for by the CG-PI approach.

## IV. DISCUSSION AND CONCLUSIONS

The present CG-PI approach to imaginary time path integrals and thermal density matrices provides a new reductionist perspective on quantum statistical mechanics. This approach is based on a coarse-graining of Feynman’s imaginary time path integral, providing a simple way to analyze the most essential nature of quantum effects. As is demonstrated in this paper, the minimalist description employs only two quasiparticle coordinates: *Q*, the coordinate of the centroid particle that feels the physical potential, and $ q \u0304 $, the coordinate of a second quasiparticle that is used to evaluate the values of observables. The expectation value of a quantum operator in this representation equals an average of the classical counterpart of this operator over a classical thermodynamic ensemble. Unlike theories aiming at a representation of a quantum particle with only the centroid variable,^{31} the present model naturally avoids the problem of evaluating quantum operators in the centroid-only representation.

The main advantages of the CG-PI theory are the following: First, it is important from the conceptual viewpoint, since it defines a simple framework and provides a toolkit for studying quantum systems (especially near the classical limit) in a way maximally close to the well-developed and computationally efficient tools of classical statistical mechanics. For example, the use of two quasiparticles (the centroid and observed ones) to describe thermal quantum delocalization of an atom in a molecular system is a way to provide a minimalist approach for understanding quantum systems in a classical framework, including the necessary rules arising from the non-commuting position and momentum operators.

Next, using Eq. (26), one can directly run MD or MC simulations for quantum systems, employing only two quasiparticles for each physical quantum particle. Even in the simplest version of the method, namely, with the use of constant spring stiffness $ k Q q \u0304 $ and constant effective mass *m _{eff}*, the results for a harmonic oscillator are exact at an arbitrary strength of quantum effects (that is, to all powers of ℏ). Renormalization of $ k Q q \u0304 $ and, for momentum-dependent observables, of

*m*with changes in $ V \u2033 ( q \u0304 ) $ accounts in a compact, minimalistic way for the major part of quantum corrections especially near classical limit in the case of a general physical potential. The CG-PI approach also allows for adding higher order in ℏ corrections in a systematic manner by explicitly deriving higher order terms not shown explicitly in Eq. (15), as illustrated in the Appendix, Eq. (A15), and Sec. S1 in the supplementary material.

_{eff}^{32}These corrections result in anharmonic terms in the effective potential of interaction between the centroid and observed quasiparticles,

*p*

^{4}and higher order terms in the effective kinetic energy in the Hamiltonian in Eq. (26), and additional terms in the expression for the effective mass

*m*. Another way of adding higher order correction may proceed via variational determination of the CG potential or by implementing various

_{eff}*Ansätze*to approximate and parameterize these terms. Likewise, the CG-PI theory can serve as a basis for the development of various approximate semi-phenomenological models of specific physical phenomena. In this case, an

*Ansatz*expression for the CG potential can be chosen based on a qualitative analysis of the problem, while the specific values of the coefficients in such an expression can be related to the microscopic description of the system via a variational procedure. The computational cost of the CG-PI method is comparable to that of path integral molecular dynamics (PIMD) with only

*P*= 2 quasiparticles per quantum particle, while the accuracy of the former is much higher, as discussed above. The use of the variable spring stiffness $ k Q q \u0304 q \u0304 $, and possibly the effective mass $ m eff q \u0304 $, does not significantly increase the cost of computations, at least in the case of analytically defined potentials

*V*(

*q*), since the right hand sides of Eqs. (18) and (19) can be quickly computed with the use of standard computational techniques, such as look-up tables and piecewise interpolations.

Finally, the CG-PI approach can be used in a combination with the existing PIMD, path integral Monte Carlo (PIMC), etc., methods as a tool for processing and analyzing the data generated by these more exact numerical methods. Here, the advantages of coarse-graining are exhibited at the stage of the analysis and interpretation of exact numerical data. This mode of using the CG-PI theory is analogous to the use of CG models to analyze all-atom trajectory data of biomolecules (see, e.g., Ref. 37). Another potentially useful strategy of combining the CG-PI approach and the traditional PIMD and PIMC methods would be to use approximate CG-PI models, e.g., based on Eq. (15), to enhance the sampling in the PIMD or PIMC simulations.

Possible directions for the future development of the CG-PI theory will address the following questions: computations for benchmark quantum systems need to be performed with a subsequent comparison of the results with numerically exact PIMD results. Generalizations of the CG-PI approach to treat many-particle systems of bosons and fermions might also be worked out. There are also possible modifications of the CG-PI method based on various other CG approaches; the theory of ultra-coarse-graining (UCG)^{38,39} may prove useful for extending this approach, for example, in the description of tunneling of quantum particles between two or more deep local minima in terms of hops between different discrete states. The UCG theory is developed precisely for the situation in which the coarse-graining is very “aggressive” in going from many particles to just a few at a low CG resolution, as in the CG-PI.

To conclude, in this work, we suggest the most essential, minimalist picture of how quantum effects are manifested in a system at finite temperature and at equilibrium. The proposed CG-PI methodology may offer a new spectrum of tools for qualitative and quantitative investigation of quantum systems, especially near the classical limit such as molecular systems at room temperature.

## Acknowledgments

This research was supported by the National Science Foundation (NSF Grant No. CHE-1465248). The authors are grateful to Professor David M. Ceperley, James F. Dama, Stephan Köhler, Dr. Yuxing Peng, and Dr. Frank X. Vazquez for discussing this project on different stages of its completion.

### APPENDIX: ADDITIONAL ANALYTICAL RESULTS IN THE CG-PI THEORY

#### 1. Exact CG potential for a harmonic oscillator

In the case of a one-dimensional harmonic oscillator with the potential energy *V*(*q*) given by

where *m* is the particle mass and *Ω* is the angular frequency of the oscillator, the exact CG potential *V _{CG}* assumes the following form:

#### 2. General perturbative expansion for the CG potential

If the physical potential *V*(*q*) can be represented as

such that the CG potential $ V CG , 0 q , Q , q \u2032 $ for the reference potential $ V 0 q $ is known, then the CG potential for *V*(*q*) can be found in a form of a perturbative series in the small parameter *λ*. One variant of the perturbative series is given by

where the averaging is defined by

Alternatively, the perturbative series expansion for the CG potential can be written as

where $\sigma =\u0127 \beta / m $, the averaging is defined by

and

#### 3. Variational approach to determining the CG potential based on force-matching

In the case of an arbitrary physical potential *V*(*q*), the CG potential (up to an additive constant) can be found variationally,

where the variational functional *χ*^{2} is defined as

#### 4. Variational approach to determining the CG potential based on Gibbs-Bogolyubov-Feynman (GBF) inequality

An alternative variational approach can be based on a generalization of the GBF inequality to the context of the CG-PI theory. This generalization leads to

where *g* is an arbitrary function. Thus, the GBF approach can be used in the cases when the dependence of *V _{CG}* on Δ

*q*is not important (e.g., for the computation of diagonal elements of the density matrix or expectation values of momentum-independent operators). The variational functional

*χ*is defined by

_{GBF} where *F _{trial}* is the free energy corresponding to the trial CG potential, specifically

and the averaging in Eq. (A12) is defined by

#### 5. Third order terms in the small *ħ* expansion of the CG potential

The series expansion of *V _{CG}* in terms of ℏ as a small parameter [see Eq. (15)] with the

*O*(

*ħ*

^{3}) terms explicitly shown assumes the following form:

where $\sigma =\u0127 \beta / m $,

Note that, in practice, the values of *κ*_{3}, *φ*_{312}, and *φ*_{310} can be computed as fast as the value of $tanh f eff ( q \u0304 ) $ by means of using reference tables and asymptotic expansions. For the derivation of these and some higher order corrections, see Sec. S1 in the supplementary material.^{32}