Functional density matrix formulation of quantum statistics

A. Bessa,* C. A. A. de Carvalho, and E. S. Fraga Escola de Ciências e Tecnologia, Universidade Federal do Rio Grande do Norte, Caixa Postal 1524, 59072-970 Natal, RN, Brazil Instituto de Física, Universidade de São Paulo, Caixa Postal 66318, 05315-970 São Paulo, SP, Brazil Instituto de Física, Universidade Federal do Rio de Janeiro, Caixa Postal 68528, 21941-972 Rio de Janeiro, RJ, Brazil Received 23 December 2008; revised manuscript received 8 November 2009; published 5 January 2010


I. INTRODUCTION
Quantum statistical physics provides the conceptual and computational framework for the treatment of interacting many-body systems in thermal equilibrium with a heat bath.It describes a great variety of systems and phenomena over a wide range of scales, thus covering practically all areas of physics, from cosmology to particle physics, with an extensive number of examples in condensed matter.
Its broad spectrum of applications includes specialrelativistic systems.For those, the quantum statistical treatment is known as finite-temperature field theory ͓1͔, which is synonymous to quantum statistical field theory.In fact, just as quantum statistical mechanics adds stochastic probabilities to the quantum probabilities of quantum mechanics, finite-temperature field theory does likewise with respect to quantum-field theory, the natural combination of special relativity, and quantum mechanics.
Quantum statistical physics can be formulated in the language of Euclidean ͑imaginary-time͒ functional integrals, quite natural for finite-temperature field theories but also applicable in nonrelativistic contexts since quantum mechanics can be viewed as a field theory in zero spatial dimensions.Functional integrals thus furnish a powerful unifying formalism to compute correlations.
Indeed, the formalism expresses the partition function of a given system as a generating functional of euclidean Green's functions ͑correlations͒, similarly to zero-temperature field theory, but with the time direction made compact, and specific boundary conditions that constrain the domain of field configurations in the functional integral.
In order to best exploit the unifying aspect of the formalism, we propose to view the partition function as an integral of the diagonal density matrix element of the theory.The integral is performed over the stochastic variable which characterizes the Schrödinger representation of the matrix element ͑position representation for quantum mechanics, field representation for second-quantized quantum mechanics or field theory͒.The density matrix is just the Boltzmann operator, whose Hamiltonian operator may either come from ͑second-quantized͒ quantum mechanics or from field theory.
The density matrix element may be written as a functional integral in euclidean time which starts from a given configuration ͑a point in quantum mechanics, a field configuration in field theory͒ at = 0 and returns to that same configuration at = ␤ប.Clearly, there is a sum to be performed over configurations which coincide at the end points of the Euclidean time interval.Those boundary configurations are identified as the stochastic variable of the remaining integral.
This density matrix method for deriving an effective stochastic problem was already used in ͓2͔, where it led to the construction of dimensionally reduced effective actions.Likewise, in Ref. ͓3͔ we used the density matrix method, and a semiclassical approximation, to investigate the thermodynamics of scalar fields.In the latter reference, the boundary configuration of the field ͑considered as a fluctuation around a homogeneous background͒ played a very nontrivial role in the semiclassical formulation for the partition function of scalar fields.
Different approximation schemes can be used in order to obtain an effective theory for the stochastic field.Our prior uses of the density matrix method drew upon generalizations of a thermal semiclassical treatment previously developed and applied to quantum statistical mechanics, which produced excellent results for the ground-state energy and the specific heat of the anharmonic ͑quartic͒ oscillator ͓4-6͔.In those prior uses, however, the connection with the stochastic variable was indirect, through classical solutions of equations of motion.In this paper, we have built a simple effective stochastic theory by expanding the action directly around the stochastic variable.Besides providing a realization of new conceptual ideas, we believe that this approach may be physically relevant to the discussion of infrared problems in several applications of quantum statistical physics.
We perform an explicit calculation of the free energy and specific heat for the anharmonic oscillator ͑for references in the literature, see ͓7-18͔͒ from the partition function of an interacting one-dimensional quantum-mechanical system at finite temperature from a path-integral representation for the density matrix.The method of calculation that we propose provides an alternative to the usual sum over periodic trajectories: it sums over paths with coincident end points, and includes the contribution of nonvanishing boundary terms.
The paper is organized as follows: in Sec.II, we introduce our formulation for quantum statistical field theory, in order to clarify the multiple role of our boundary fields, which are at the same time field eigenstates for the subjacent field theory problem, and stochastic variables for the ensuing statistical integral; in Sec.III, we particularize the discussion to quantum mechanics, presenting the contribution of nonperiodic trajectories in the density matrix formalism which lead us to a modified series expansion in Matsubara frequencies; besides, we review the high-temperature limit of the quantum theory as a classical effective theory for the boundary variable; in Sec.IV, using a simple effective theory for the stochastic variable, we compute several thermodynamic quantities for the quadratic and the quartic potentials, and compare the results to the semiclassical findings of ͓4͔, and to some exact results; and Sec.V contains our conclusions and outlook.

II. DENSITY MATRIX FOR SCALAR THEORIES
In statistical mechanics, the partition function for a system in contact with a thermal reservoir at temperature 1 T ͑␤ =1/ T͒ is expressed as a sum ͑integral͒ over a stochastic variable, whose probabilistic weight is given by the density matrix.
For a system described by a ͑self-interacting͒ scalar field theory, the stochastic variable is the field eigenvalue in the functional Schrödinger picture, given by the function 0 ͑x͒, which satisfies ˆ͉ 0 ͑x͒͘ = 0 ͑x͉͒ 0 ͑x͒͘.
The partition function is then a functional integral over the field eigenvalue of the diagonal element of the density matrix functional ␤ , The Hamiltonian in the preceding formula specifies the underlying dynamics.In the present case, it is the dynamics of a ͑self-interacting͒ scalar field ͑t , x͒.
As is well known ͓19͔, ␤ can also be expressed as a functional integral over dynamical fields ͑ , x͒ defined for Euclidean time , subject to the boundary conditions ͑0,x͒ = ͑␤ , x͒ = 0 ͑x͒, The boundary conditions relate the stochastic variable of the integral in Eq. ͑1͒ to the dynamical fields that are integrated over in Eq. ͑3͒.The Hamiltonian which involves the time-dependent field ͑t , x͒ and the conjugate momentum ⌸͑t , x͒, leads to the Euclidean action ͑we assume d spatial dimensions͒, The density matrix is a functional of 0 ͑x͒ only.The remaining integral over 0 ͑x͒ which is required to obtain the partition function is unrestricted ͑except for the vacuum boundary conditions that must be imposed at spatial infinity͒.
We may write = e −S d , S d being a certain T-dependent dimensionally reduced action.The field 0 ͑x͒ which is the argument of S d depends only on the d spatial coordinates; all the dependence of the original d + 1 theory has been eliminated through the integration.
The fields 0 ͑x͒ are the natural degrees of freedom of the reduced theory.Any thermal observable can be constructed by integrating over the fields 0 ͑x͒ an appropriate functional of 0 ͑x͒ weighed with the corresponding diagonal element of the density matrix.The Euclidean time evolution can be viewed as an intermediate step which calculates the weights.
Notice that the density matrix is not, in general, of the form e −␤H d with H d being independent of ␤ as in ordinary statistical mechanics; its ␤ dependence is far more complicated.This had already been pointed out in the analogous discussion of the transfer matrix carried out in Ref. ͓20͔.The density matrix provides a direct but alternative way of deriving a dimensionally reduced theory.
In previous works, we calculated the density matrix element functional given by the integral in Eq. ͑3͒ either from a perturbative expansion ͓2͔ or from a semiclassical expansion ͓3͔ around the extremum of the integrand for the special case of a homogeneous ͑x-independent͒ field 0 .
The semiclassical approximation corresponds to the usual saddle-point approximation with Gaussian fluctuations.For the general case, one has to obtain the extrema of the euclidean action in Eq. ͑3͒; in other words, classical solutions subject to the boundary conditions ͑0,x͒ = ͑␤ , x͒ = 0 ͑x͒.Since the density matrix can only depend on 0 ͑x͒, it should be clear that the classical solution is itself a functional of 0 ͑x͒.Thus, the semiclassical expansion would lead to a somewhat indirect dependence of the density matrix on 0 ͑x͒.Furthermore, the general problem of finding classical solutions for arbitrary boundary conditions 0 ͑x͒ does not have a general solution.To circumvent these difficulties, in the 1 In order to simplify the notation, we will adopt natural units where ប =1, c = 1, and k B =1.Using 0 ͑x͒ as a background eliminates the dependence that would be present in the classical solution, and provides expressions that depend directly on 0 ͑x͒.On the other hand, we are forced to deal with linear terms when we resort to a quadratic approximation to the integral over ͑ , x͒ which would be absent in a quadratic approximation to the integral over ͑ , x͒.We emphasize that our alternative expansion is not a saddle-point approximation.However, when T is much bigger than any other scale of the theory, cl ͓ 0 ͔ is approximately given by 0 ͓x͔ Dealing with a linear term is not a problem within a quadratic approximation, and can be handled by the usual shift in the integration variable.As for the validity of the expansion, it should be clear that it should be quite reliable at high temperatures, where ͑ , x͒ cannot differ appreciably from zero since the Euclidean time interval ␤ is too small for deviations from the boundary values to become large.In fact, in a high-temperature regime, the classical solution is not very different from 0 ͑x͒.
For the sake of simplicity, we have chosen to test our strategy for computing the density matrix element in quantum mechanics, which can be viewed as a field theory in zero spatial dimension.By doing so, we do not have to deal with the subtleties of renormalization.We also need not worry about the infrared problems that plague various treatments of field theories at finite temperature ͓1,21,22͔.Furthermore, we can compare the outcome our proposal with previous calculations done using ordinary perturbation theory, as well as to a semiclassical treatment of the integral.
In the sequel, we calculate the partition function over a range of temperatures where we believe that our approximation will hold.In fact, we even estimate a lower bound on temperature above which we expect our computations to be valid and verify that this is indeed the case through numerical calculations.
Before proceeding in the calculation, we discuss some important features of the present approach: the inclusion of nonperiodic configurations and its exact high-temperature limit.

A. Inclusion of nonperiodic configurations
If we repeat the steps of the previous section, the partition function Z for a one-dimensional quantum-mechanical system in contact with a thermal reservoir at temperature T =1/ ␤ will be written as

͑9͒
where ␤ ͑␤ ; x 0 , x 0 ͒ is the diagonal element of the density matrix.Note that the role of the field 0 ͑x͒ of the previous section is now played by a mere point x 0 , and that accordingly functional integrals become ordinary integrals.
In the path-integral formulation, the density matrix ␤ is obtained from the imaginary-time evolution of the trajectory x͑͒ : ͓0,␤͔ → R ͓which replaces ͑ , x͔͒ determined by the Euclidean action S, where V is the potential, and the dot refers to differentiating in .As we have emphasized for 0 ͑x͒, it is natural to regard x 0 as the effective stochastic degree of freedom of the theory, and to compute the distribution of the variable x 0 by integrating over the dynamical variable x͑͒, with boundary conditions determined by x 0 .
Integrating action ͑11͒ by parts, we obtain,

͑12͒
Among the trajectories with boundary value x 0 at = 0 and ␤, the one which contributes the most for the partition function is the ͑classical͒ solution of In the special case of the harmonic oscillator ͑V͑x͒ = m 2 x 2 / 2͒, it is a simple matter to show that which is not a periodic function of .As a consequence, one is led to deal with boundary terms in the action ͓see Eq. ͑12͔͒.We stress that, as long as the only condition over paths in the path integral is to have coincident values at 0 and ␤, nonperiodicity is allowed.
In order to see that nonperiodic trajectories are a general feature of the present approach, let us now decompose the path x as where Clearly, fluctuations obeying Eq. ͑16͒ can be expanded in a sine-Fourier series, x n sin ˆn, ͑17͒ where ˆn = n / ␤.The notation with a dot is a reminder that the r.h.s. of Eq. ͑17͒ is a Fourier approximation of the l.h.s.The complete path is, then, given by x n sin ˆn, ͑18͒ which does not correspond to a periodic function in ͓0,␤͔ ͑see Fig. 1͒.Series ͑18͒ does not impose conditions on the derivative of x.It is easy to convince oneself that evaluating the euclidean action S with the r.h.s. of Eq. ͑18͒ one obtains exactly S͓x͔͑͒.In particular, we have for the classical solution, sin ˆn ͑19͒ so that the quantity x ˙c͑␤͒ − x ˙c͑0͒ is given by leading to the correct result.Obviously, the series defined by Eq. ͑18͒ does not correspond to the most general trajectory in Z.In particular, that series imposes a serious restriction on the second derivative of x by forcing it to vanish at 0 and ␤.However, at least for the calculation of usual action functionals, the boundary terms do not involve second derivatives and the aforementioned restriction is actually unimportant.
A remarkable property of decomposition ͑18͒ is that the boundary value x 0 plays the role of the zero, static component of x͑͒. 2 Likewise, 0 ͑x͒ is the static component of ͑ , x͒ in the context of thermal field theories, which thus suggests that it plays a special role in the infrared ͑low-momentum͒ limit of the theory.

B. Exact high-temperature limit
When the typical spacing between quantum levels, ប, is such that ͑or ␤ Ӷ 1, in natural units͒, we are in the classical limit of the theory.It is a known fact that the classical regime is correctly described by the effective theory for the static modes.Indeed, let us consider x c , the solution of the equation of motion ͑13͒.In the high-temperature regime ͑␤ Ӷ 1͒, thermal fluctuations dominate.In this limit, x c → x 0 and S͑x c ͒ → ␤V͑x 0 ͒, and we obtain for the partition function

͑22͒
where N is a normalization factor which incorporates quantum fluctuations.A proper normalization of Z leads to
As discussed in ͓23͔, for ␤ Ӷ 1 the correlation ͗x 2 ͑͒͘ follows the classical linear scaling with T. This behavior is entirely associated with the static component of x, and represents a problem for the high-temperature limit of the usual perturbative expansion.In contrast, the subtracted thermal propagator ͑without the contribution from the static mode͒ goes to zero as T increases so that one should strongly improve the convergence of the perturbation expansion by calculating diagrams with that subtracted propagator.
From the aforementioned reasons, we conclude that a separate treatment of the boundary ͑static͒ modes is convenient.The main advantage of the density matrix approach is to provide this separation in a natural fashion.In terms of y defined in Eq. ͑16͒, one can express Z conveniently as and different strategies can be used to handle S͓x 0 + y͔ in order to obtain an effective theory 3 for x 0 , e −␤V eff ͑x 0 ͒ .͑25͒

IV. HIGH-TEMPERATURE EXPANSION
In the last section, we have shown that the hightemperature limit of the theory is an effective theory for the 2 Rigorously speaking, x 0 is not a true Fourier mode since it is not orthogonal to sin ˆn in ͓0,␤͔ but in ͓−␤ , ␤͔.With this caveat in mind, we will continue to refer to x 0 as a mode of x.

3
With the normalization factor used in Eq. ͑25͒, one has the property: V eff ͑x 0 ͒ → V͑x 0 ͒ when T → ϱ, as follows from Eq. ͑23͒.boundary static variable x 0 .In the sequel, we use a very simple expansion around x 0 in two cases: the harmonic potential, as a consistency check, and the single-well quartic potential, in order to verify the quality of our alternative approach when we extrapolate the high-temperature regime.

A. Quadratic case
The action in this case is given by Without any loss of generality, one can take m = 1 in the previous equation by rescaling x by a factor 1 / ͱ m.Now, we use the decomposition of x given by Eq. ͑16͒ to obtain

͑27͒
where G y is the solution of G y ͑0,Ј͒ = G y ͑␤,Ј͒ = 0. ͑28b͒ A simple calculation gives ͓where Ͼ = max͑ , Ј͒ and Ͻ = min͑ , Ј͔͒, and Shifting the variable of integration to which obeys ͑0͒ = ͑␤͒ = 0, and performing the Gaussian integration, we obtain Integrating over x 0 and taking the logarithm, one reproduces the known result for the free partition function.Therefore, for quadratic theories the present high-temperature expansion is exact.

B. Anharmonic oscillator
Let us now consider the following action: It is convenient to define the dimensionless variables: q = ͱ / ͑m 2 ͒x, = , ⌰ = ␤, and g = / ͑m 2 3 ͒.The asso- ciated action for the variable q is Using the decomposition q͑͒ = q 0 + ͱ g͑͒, we write for the action S͓q͔ = S 2 ͓q 0 ,͔ + S I ͓q 0 ,͔, ͑35͒ where S 2 contains S͓q 0 ͔, the kinetic term, and those which are linear or quadratic in .The remaining contributions to S͓q͔ are collected in S I .It is a simple matter to show that where ¯0 2 ͑q 0 ͒ =3q 0 2 + 1 and ␣ 0 ͑q 0 ͒ = q 0 3 + q 0 ͑here and in the sequel, the subscript is a reminder of the dependence on q 0 ͒.One obtains a quadratic approximation for the partition function by considering Proceeding as in the previous section, we perform the Gaussian integration over the variable where the Green's function G 0 ͑again, the subscript is a reminder of the dependence of this propagator on q 0 ͒, defined by is given by where the q 0 dependence comes through ¯0͑q 0 ͒.The result for Z 2 is then where det G 0 = ¯0 2 sinh͑ ¯0⌰͒ ͑43͒ and From Eq. ͑42͒, one directly obtains an expression for ␤ ͑␤ ; q 0 , q 0 ͒.The positive function ␤ is peaked at q 0 = 0 and decays to zero around q 0 =2 ͱ g.
To go beyond the quadratic order in , we use the following approximation:

͑47͒
As usual, we can introduce a linear coupling of with an external current j: and obtain the correlations ͓Eq.͑47͔͒ as functional derivatives of Z 2 ͓j͔ with respect to j.With minor modifications on the recipe for Z 2 , one obtains where and 0 ͑͒ = ␣ 0 / ͱ g − j͑͒.Finally, we use

͑52͒
where I 0 ͑͒ plays the role of an expectation value of ͑͒, The remaining integrations over can be done analytically, providing all the ingredients for the calculation of the first correction to Z 2 .From Eq. ͑45͒, we read that corrections to the quadratic approximation become important when the expectation value of interaction ͑37͒ is about 1.As the integrand of Eq. ͑46͒ is dominated by the vicinity of q 0 = 0, where the quadratic action reaches its minimum value, one can estimate the importance of corrections using ͗S I ͑0,͒͘.We have The temperature ⌰ min where ͗S I ͑0,͒͘ goes to 1 provides a rough estimate for the validity of the quadratic approximation.Figure 2 displays ⌰ min as a function of the dimensionless coupling constant g for the high-temperature and the usual saddle-point approximation ͑see Ref. ͓4͔͒.From Fig. 2, we conclude that the two methods have almost the same range of applicability.In particular, we see that ⌰ min is of order 1 for g as large as 50.
The results for the free energy when g = 0.4 are displayed in two plots with different ranges of temperature: Fig. 3 shows the high-temperature region, while Fig. 4 is devoted to low temperatures, according to the criterium used in Fig. 2. The results for the specific heat are displayed in Fig. 5.In the figures, classical corresponds to calculation using Eq.͑23͒, FIG. 2. ͑Color online͒ Plot of ⌰ min as a function of the dimensionless coupling constant g.Each approximation is applicable for temperatures larger than the associated ⌰ min .The inset is the same plot for smaller values of g. high-temperature uses Eq. ͑42͒, improved includes the correction corresponding to Eq. ͑46͒, semiclassical refers to the semiclassical result obtained in Ref. ͓4͔, perturbative stands for the usual first-order perturbative result, and exact corresponds to results from Refs.͓13,14͔ combined with WKB estimates.
The first remarkable thing in the plots is the failure of the first-order perturbative result in the high-temperature region ͑T ӷ ͒.As discussed in ͓23͔, the static component of the usual thermal propagator spoils the convergence of the perturbative series for temperatures large compared with .As seen in Ref. ͓4͔, the high-temperature region is well described by the semiclassical curve, and asymptotically, by the classical one.This is evident from the plots.It is interesting that, for g = 0.4, the agreement extends to temperatures down to ⌰ Ϸ 1 / 2, a region where the character of the system is far from classical.However, it is surprising that our simple calculation using Eq.͑42͒ works as well as the semiclassical one in that region.
Indeed, the cited semiclassical calculation is based on a quadratic expansion around exact solutions of Eq. ͑13͒ for the quartic potential.In practice, one has to deal with rather involved Jacobi elliptic functions ͓4,24͔.In contrast, using expansion ͑15͒ no knowledge of classical solutions is demanded.This simplification represents an economic alternative which can be crucial contexts where exact classical solutions are not available.From the plots, one also concludes that, in the region where the approximation is supposed to be good, the improved calculation exhibits a stronger convergence, indicating that the approximation is consistent.
A detailed look at the quadratic curve in the free-energy plot ͑Fig.4͒ reveals a nonmonotonic behavior for low temperatures: the free energy passes by a maximum and decreases toward the free value 0.5.This is a general feature regardless of the value of g, and it is a clear signal of the breakdown of the approximation below ⌰ min ͑see Fig. 2͒.The improved free energy diverges at ⌰ min because at this value the sum Z 2 + ␦ ͑1͒ Z vanishes.The exact curve for the free energy reaches its maximum value at the ground-state energy, where it rests down to ⌰ = 0. Using the maximum value assumed by the quadratic curve as an approximate value for the ground-state energy, we obtain unexpectedly reasonable results, as shown in Table I.A path-integral general method for deriving the ground-state wave function and energy can be found in Ref. ͓25͔.
Figure 5 compares different predictions for the specific heat for g = 0.4.Notice that all curves, except for the perturbative one, have the correct high-temperature limit.The anomalous behavior of the quadratic approximation near ⌰ min is even more evident in this plot.The improvement obtained in the low-temperature regime when one corrects the quadratic calculation with Eq. ͑46͒ is also clearly shown.

V. CONCLUSIONS AND OUTLOOK
The theoretical description of the thermodynamics of a given interacting many-body system kept in thermal equilibrium with a heat bath is an issue of capital importance in most realms of physics, from cosmology to condensedmatter settings.A very convenient framework is provided by finitetemperature field theory, where one can use the imaginarytime formalism and make a direct connection to the powerful formulation of thermal averages in terms of functional integrals.
Rigorously speaking, the paths need not be periodic in the path-integral representation of the partition function.Indeed, the necessity of paths which are not periodic is particularly clear in the formulation of the partition function as an integral of the diagonal density matrix element of the theory.
In this paper, we have followed our alternative strategy in the case of quantum statistical mechanics, viewed not only as a toy model for the case of finite-temperature field theory, but also as a prototype for various relevant systems in statistical physics.More specifically, we have explored an alternative way of computing the partition function in the pathintegral formalism that includes nonperiodic trajectories.As was shown above, one can properly incorporate the contribution of nonperiodic paths by using a modified Matsubara series expansion in which the static mode is identified with the boundary value of the path.The latter turns out to be the stochastic variable that survives in the final effective theory which has the character of a high-temperature approximation.
We built a very simple effective theory in the nontrivial problem of the anharmonic oscillator, obtaining very precise results for practically the whole range of temperatures without any information about classical solutions.We expect that other potentials can profit from the simplicity of the present method in quantum-mechanical systems.
When more than one classical solution is present ͑such as in the case of double-well potentials͒, the method probably would still apply in the high-temperature regime down to the first bifurcation point ͑a discussion on caustics in the semiclassical approximation is given in ͓6͔͒.However, we do not expect the method can correctly describe zero-temperature features such as tunneling driven by instantonlike solutions.
As mentioned previously, the infrared physics of bosonic theories is not accessible to plain perturbation theory, and resummation techniques are required in order to produce sensible results.In common, all such techniques render a special treatment to the static mode.Therefore, we expect that in field theory the physics of 0 ͑x͒ will play a prominent infrared role, although this is still under investigation.
In the case of finite-temperature field theory there are, of course, subtle issues of renormalization, which are absent in the quantum-mechanical setting, that will have to be addressed.Nevertheless, the results obtained in quantum statistical mechanics are encouraging, and we hope that this new perspective may shed light on the problem of infrared divergences.Results in this direction will be presented elsewhere ͓26͔.
FIG.3.͑Color online͒ Free energy of the anharmonic oscillator as a function of the dimensionless temperature for g = 0.4 in the high-temperature region ͑see text for details͒.

TABLE I .
Ground-state energy of the quartic oscillator for different values of the coupling ͑m = =1͒.Exact data were extracted from Ref. ͓14͔, whereas E 0 ͑quadratic͒ was obtained using Eq.͑42͒.