Fast convergence of path integrals for many-body systems

# Fast Convergence of Path Integrals for Many-body Systems

## Abstract

We generalize a recently developed method for accelerated Monte Carlo calculation of path integrals to the physically relevant case of generic many-body systems. This is done by developing an analytic procedure for constructing a hierarchy of effective actions leading to improvements in convergence of -fold discretized many-body path integral expressions from to for generic . In this paper we present explicit solutions within this hierarchy up to level . Using this we calculate the low lying energy levels of a two particle model with quartic interactions for several values of coupling and demonstrate agreement with analytical results governing the increase in efficiency of the new method. The applicability of the developed scheme is further extended to the calculation of energy expectation values through the construction of associated energy estimators exhibiting the same speedup in convergence.

###### keywords:
Path integral, Effective action, Many-body system, Energy estimator
###### Pacs:
05.30.-d, 03.65.Db, 05.10.Ln

, , , and

## 1 Introduction

Originally introduced in quantum mechanics [1, 2] and later most widely used in high energy theory [3, 4] and condensed matter physics [5, 6], path integrals have become important tools throughout the physical sciences from atomic, molecular and nuclear physics, to chemistry and biophysics. Moreover, path integrals are starting to play important roles in several areas of mathematics and in modern finance [7], especially in option pricing applications [8, 9]. The downside of the formalism is that the mathematical properties of path integrals are not sufficiently well understood, and that an extremely small number of path integrals can be solved exactly [10]. An extensive overview of the path integral formalism and its various applications can be found in [11].

Beside their key position in analytical approaches to quantum theory, path integrals have an important role in direct numerical simulations of realistic many-body systems. Such simulations have made possible practical comparisons of the predictions of various theoretical models with the results of associated experiments. Further, they have led to a deeper understanding of the complex physical phenomena involved [12].

The starting point in all such calculations is the time-sliced expression for the general quantum-mechanical transition amplitude [2],

 AN(a,b;T)=1(2πε)MNd2∫dq1⋯dqN−1 e−SN. (1)

from initial state to final state , for time interval , where is the number of time slices, and is the naively discretized action for a system of non-relativistic particles in spatial dimensions. The limit of the above discretized amplitude gives the continuum amplitude . The evaluation of these types of discretized expressions is handled by a variety of well developed numerical integration methods, and this is one of the principle reasons for the popularity of the path integral approach in numerical simulations.

The performance of numerical algorithms for the calculation of path integrals is determined by the efficiency of the applied integration techniques, as well as by the speed of convergence of the discretized amplitudes to the continuum. The first aspect has been successfully dealt with by numerous Monte Carlo integration methods based on efficient sampling of trajectories. The second is usually addressed by the use of better discretizations, i.e by substituting better effective actions for the naively discetized action. The underlying idea here is to construct and use effective actions that lead to improved convergence of physical expressions to the (same) continuum limit.

Typically, physical expression calculated using the naively discretized action converge to the continuum as . A review of a variety of effective actions constructed by improving short-time propagation and using generalizations of the Trotter-Suzuki formula [13], together with their ranges of application and dependence on ordering prescriptions, are given in [14]. For a long time now, the state of the art result has been the convergence of discretized partition functions obtained by Takahashi and Imada [15], and Li and Broughton [16]. In these papers the authors used a generalized form [17] of the Trotter formula and the cyclic property of the trace to increase the speed of convergence of the discretized partition functions to . We stress that this increase in the speed of convergence only holds for partition functions and not for amplitudes. Some recent results regarding the efficient implementation of PIMC algorithms may be found in refs. [18, 19, 20, 21, 22]. In general, the main feature of the improved effective action – faster approach to the continuum limit, i.e with smaller number of time slices – translates into numerical speedup, simpler integration and smaller statistical fluctuations of the quantities calculated.

In a recent paper [23], we have introduced the concept of the ideal effective action . Discretized amplitudes calculated with the ideal effective action do not depend on the discretization coarseness , i.e. they are all equal to the sought-after continuum amplitude

 A∗N(a,b;T)=A(a,b;T). (2)

Note that for free particles is in fact just the naively discretized action. In [24] we derived the integral equation for the ideal action for a single particle moving in one spatial dimension in a general potential. We then solved this integral equation as an asymptotic expansion in . The obtained truncation of the ideal effective action to order systematically improves convergence to the continuum to :

 A(p)N(a,b;T)=A(a,b;T)+O(1/Np), (3)

where corresponds to the naive action. Up to now the work has focused on one particle theories in one dimension. Within this set of theories the outlined procedure has been used to determine explicit expressions for the set of effective discretized actions up to level . The ensuing speedup of several orders of magnitude has been verified by extensive numerical simulations [25]. The speedup holds for the calculation of all path integrals – for transition amplitudes, partition functions, energy expectation values (in combination with appropriate energy estimators), as well as for calculations of energy levels. In this paper we extend the above method to the physically relevant case of general many-particle non-relativistic systems in arbitrary number of spatial dimensions. This is done within a new analytical approach to the construction of the effective actions hierarchy. The new approach gives the same effective actions as the the old one when applied to one particle one dimensional theory, however its computational complexity is lower, making it possible to determine effective actions at higher levels. In particular, the new approach allows for a relatively straight-forward extension of the formalism to many-body theories and arbitrary dimensions.

The present paper gives the analytical derivation of these results, as well as a series of Monte Carlo simulations implementing the newly derived effective actions and explicitly displaying that the derived speedup indeed holds. The paper is organized as follows: Section 2 introduces the hierarchy of effective discretized actions and gives the new and extended (many particles, higher dimensions) analytical approach to the construction of these effective actions. Section 3 presents numerical results that demonstrate the validity of the analytically derived speedup in convergence. The efficiency of the new approach is further demonstrated in Section 4 through the calculation of energy spectra. A construction of virial energy estimators corresponding to newly derived effective actions and the results of numerical simulations that implement them are given in Section 5.

We conclude the paper with a brief summary of obtained results, indicating what we see to be the next steps in our line of research and its future applications. Appendix A gives a list of integration formulas needed for calculations up to level . Appendix B lists the corresponding effective actions representing the new state-of-the-art for path integral calculations of non-relativistic many-body systems in arbitrary number of dimensions. Higher level effective actions can be found on our web site [26].

## 2 New approach to the derivation of effective actions

We present a method for systematically increasing the efficiency of PIMC calculations of a non-relativistic quantum system consisting of distinguishable particles in spatial dimensions, with a Hamiltonian of the form , where is the usual kinetic energy operator, and is the potential describing inter-particle interactions and interactions with external fields. The above set of theories encompasses a large number of physically relevant models. However, as we shall see from the derivations, the method is in principle applicable to all quantum theories. For the considered class of theories, the naively discretized action is

 SN=N−1∑n=0ε⎛⎝M∑ℓ=112(δn,ℓε)2+V(¯qn)⎞⎠. (4)

Vectors with one index (e.g. ) represent the set of positions of all particles after time steps of length , while vectors with two indices (e.g. ) represent -dimensional positions of particle at time step . We have also introduced the associated discretized velocities and mid-point coordinates . To summarize, counts the time steps, the different particles. In the most compact notation, we may think of a configuration (trajectory) of the quantum system as a single vector whose individual components take on possible values, representing positions of all the particles at a given time step. In this notation the limit of the above discretized amplitude is symbolically written as the path integral

 A(a,b;T)=∫q(T)=bq(0)=a[dq] e−S[q(t)]. (5)

Note that we are using the mid-point ordering prescription and units in which and the particle masses have been set to unity.

The defining relation for path integrals as the continuum limit of discretized amplitudes given by equation (1) follows from the completeness relation (decomposition of unity)

 A(a,b;T)=∫dq1⋯dqN−1 A(a,q1;ε)⋯A(qN−1,b;ε), (6)

through the substitution of short-time amplitudes calculated to first order in time step , leading to the naive action in equation (1). This is what gives the convergence to standard path integral expressions. A faster converging result may be obtained by evaluating the amplitudes under the integral to higher orders in . From the above completeness relation, it follows that the ideal discretized action leads to exact propagation in time, and is given in terms of the exact amplitude, according to

 A(qn,qn+1;ε)=(2πε)−Md2 e−S∗n. (7)

The ideal discretized action is simply the sum of expressions :

 S∗N=N−1∑n=0S∗n(qn,qn+1;ε). (8)

Full knowledge of the ideal action is equivalent to knowing the exact expression for all the amplitudes . At first, this would seem to indicate that nothing new is to be gained by equation (7). This, however, is not the case. We will use equation (7) to input new analytical information into our numerical procedure by calculating amplitudes within some available analytical approximation scheme. We focus on the calculation of the short time propagation as a power series expansion in , which starts from the naive action (4). The details of this calculation have been inspired by a similar derivation given in [11]. One can also attack the problem of solving equation (7) using various other approximative schemes, especially in the case when short time approximation is not appropriate. The application of the Feynman and Kleinert variational approach [27, 28], further developed in [29, 30, 31, 32, 33], could prove useful is such cases.

Improved convergence of path integrals follows from calculating the general short-time amplitude up to the order . To this end, we first perform a simple shift of integration variable about a fixed referent trajectory ,

 A(qn,qn+1;ε)=e−Sn[ξ]∫x(ε/2)=0x(−ε/2)=0[dx] e−∫ε/2−ε/2ds(12˙x2+U(x;ξ)). (9)

We have also shifted the time from to . The referent trajectory satisfies the same boundary conditions as . As a result, the new integration variable vanishes at the boundaries. The action is defined as

 Sn[ξ]=∫ε/2−ε/2ds(12˙ξ2+V(ξ)), (10)

and , with dots representing derivatives over time . The amplitude may now be written as

 A(qn,qn+1;ε)=e−Sn[ξ](2πε)Md2 ⟨e−∫ε/2−ε/2ds U(x;ξ)⟩, (11)

where denotes the expectation value with respect to the free particle action. The above expression holds for any choice of referent trajectory . Equations (7) and (11) now give the expression for the ideal effective action

 S∗n=Sn[ξ]−ln⟨e−∫ε/2−ε/2ds U(x;ξ)⟩. (12)

The class of theories considered is free of ordering ambiguities, i.e. different ordering prescriptions yield the same continuum result. For concreteness we work in the mid-point prescription, calculating the above amplitude as a power series in time step . The free particle expectation value in equations (11) and (12) is calculated using a Taylor expansion in powers of :

 ⟨e−∫ds U(x;ξ)⟩ = 1−∫ds⟨U(x;ξ)⟩+ (13) +12∫∫dsds′⟨U(x;ξ)U(x′;ξ′)⟩+…

with shorthand notation . By expanding around the referent trajectory , we get

 U(x;ξ)=xi(∂iV(ξ)−¨ξi)+12xixj∂i∂jV(ξ)+… (14)

From now on we assume summation over repeated indices. The expectation values of products are now calculated in the standard way with the help of the free-particle generating functional given in terms of the propagator:

 Δ(s,s′)ij= δijεθ(s−s′)(ε2−s)(ε2+s′) (15) Missing or unrecognized delimiter for \left

From this follow the usual Wick’s theorem results , , etc. Note that the generating functional (and so all the expectation values) is independent of the specific choice of referent trajectory, i.e. the boundary conditions for the are the same for all choices of , and so the propagator is always given by (15). Different choices of simplify different approximation schemes: using the classical trajectory for is optimal for recovering semi-classical expansion, the choice of linear referent trajectories will turn out to be optimal for short-time expansion.

We next wish to perform the remaining integrations over . Because of the explicit dependence of the referent trajectory on , we first expand the potential and all its derivatives in (14) around some reference point. The choice of for that point corresponds to the mid-point ordering prescription. Once one chooses the referent trajectory , all expectation values in (11) are given in terms of quadratures. The choice of linear referent trajectories makes all of these integrals solvable in closed form, allowing us to determine the ideal effective action.

Before doing these explicit calculations, we first look at which terms need to be retained in order to get the sought-after convergence. It is easy to show that the ideal effective action is a sum of terms of the form

 εαδ2β(∂γ1V(¯qn))⋯(∂γkV(¯qn)), (16)

where , , and are nonnegative whole numbers constrained by simple dimensional analysis to satisfy

 α+β=k+12k∑r=1γr. (17)

Short time propagation satisfies the diffusion relation . Thus, convergence follows from keeping all the terms satisfying . An equivalent, but practically more useful, form of the criterion determining the relevant terms at level is

 k+12k∑r=1γr

We now continue with the explicit evaluation of the ideal effective action, illustrating the general procedure on calculations to level . The action is now

 Sn[ξ]=ε(12δ2nε2+V(¯qn)+δn,iδn,j24∂2ijV(¯qn))+O(ε3). (19)

Note that is the -th component of vector while is shorthand for . The first two terms in the above expression correspond to the naive action, while the third term gives contributions of order . The remaining contribution at level comes from the expectation value

 ⟨e−∫ds U(x;ξ)⟩ = Missing or unrecognized delimiter for \left = 1−∫ε2−ε2ds 12 ⟨xixj⟩ ∂i∂jV(ξ)+O(ε3)= = 1−12 δij∫ε2−ε2ds Δ(s,s)i,j ∂2ijV(ξ)+O(ε3)= = 1−12 ∂2V(¯qn)∫ε2−ε2ds Δ(s,s)i,i +O(ε3).

As promised, the remaining integral is easily evaluated in closed form to yield . Appendix A lists all the related integrals needed for higher level calculations. The expectation value now equals

 ⟨e−∫ds U(x;ξ)⟩=1−ε212∂2V(¯qn)+O(ε3)=e−ε212∂2V(¯qn)+O(ε3). (20)

Finally, the level discretized effective action is simply

 S(p=2)N=N−1∑n=0ε[12(δnε)2+V(¯qn)+ε12∂2V(¯qn)+δn,iδn,j24∂2ijV(¯qn)]. (21)

One similarly obtains the higher level effective actions. Appendix B gives the explicit analytical expressions for many-particle discretized effective actions at level . There are no obstacles in going to higher values of . The derived expressions become algebraically more complex, and calculations are best done using a package for symbolic calculus such as Mathematica. Higher level effective actions can be found on our web site [26].

## 3 Numerical results

In this section we present results of numerical PIMC simulations that confirm the analytically derived speedup in convergence of discretized path integrals. To do this, we have conducted a series of PIMC simulations of transition amplitudes for a two-dimensional system of two particles interacting through potential

 V(→r1,→r2)=12(→r1−→r2)2+g124(→r1−→r2)4+g22(→r1+→r2)2. (22)

Numerical simulations, based on our SPEEDUP [26] PIMC code, have been performed for different values of couplings and and for variety of initial and final states. The associated continuum limit amplitudes have been estimated by fitting polynomials in to the discretized values , according to the analytically derived relation (3):

 A(p)N=A(p)+B(p)Np+C(p)Np+1+…. (23)

For all values of the fitted continuum values agree within the error bars. The obtained dependence gives explicit verification of the analytically derived increase in convergence. The typical case is illustrated in Figure 1. The expected convergence may be most easily discerned from Figure 2, a plot of the deviations of discretized amplitudes from the continuum limit. As can be seen, the increase of level leads to an ever faster approach to the continuum.

As a result of the newly presented method, the usual simulations (in which one calculates specific physical quantities such as the one in Figure 1) proceed much faster than by using standard methods. On the other hand, Figure 2 is itself time consuming since it illustrates subdominant behavior. For this reason the figures contain only results obtained by effective actions to level . Note, however, that the curve corresponds to a precision of four decimal places even for an extremely coarse discretization such as .

In usual PIMC calculations one always chooses the number of MC samples so that the stochastic error of numerical results is of the order of deviations from the continuum limit. In Figure 2 the number of MC samples had to be much larger, in order for deviations from the continuum limit to be clearly visible.

## 4 Energy spectra

The improved convergence of path integral expressions for amplitudes leads directly to the same kind of improvement in the convergence of discretized partition functions, owing to the relation

 ZN(β)=∫dqAN(q,q;β), (24)

where the inverse temperature plays the role of the time of propagation . From the previous relation we directly obtain the path-integral presentation of the partition function:

 ZN(a,b;T)=1(2πε)MNd2∫dq1⋯dqN e−SN. (25)

The partition function is the central object for obtaining information about statistical systems. In addition, the partition function offers a straightforward way for extracting information about low lying energy levels of a system. This follows from evaluating the trace in the definition of the partition function in the energy eigen-basis

 Z(β)=∞∑n=0dne−βEn, (26)

where and denote corresponding energy levels and degeneracies. A detailed description of the procedure for extracting energy levels from the partition function is given in [34]. The free energy of the system, given by

 F(β)=−1βlnZ(β), (27)

tends to the ground state energy in the large limit. Similarly, we introduce auxiliary functions

 F(n)(β)=−1βlnZ(β)−∑n−1i=0di e−βEidn, (28)

which can be fitted for large to

 f(n)(β)=En−1βln(1+ae−βb). (29)

and which tend to the corresponding energy level .

In this way, by studying the large behavior of the free energy and the corresponding auxiliary functions, one obtains the low lying energy spectrum of the model. However, in numerical simulations we are inevitably limited to a finite range of inverse temperatures. In addition, the above procedure for the construction of the auxiliary functions is recursive, i.e. in order to construct we need to know all the energy levels below . This leads to the accumulation of errors as increases, practically limiting the number of energy levels that can be calculated. Note that the orders of magnitude increase in precision of the presented method reduces the magnitude of the accumulated error and thus allows us to extract a larger number of energy levels.

Figures 3 and 4 illustrate the calculation of low lying energy levels using the more efficient PIMC formalism presented in this paper on the case of a two-dimensional system of two distinguishable particles interacting through a quartic potential (22). Figure 3 demonstrates that the improved convergence of free energies is the same as in the case of amplitudes. Table 1 gives the calculated energy levels for quartic coupling from (free theory) to (strongly interacting theory).

Figure 4 presents a comparison of ground energy values calculated using the derived effective actions with those calculated by perturbative expansion. The graph gives the dependence of obtained values of ground energy on the coupling constant . The agreement with perturbative results is excellent for small values of the coupling. Around even the higher order perturbative results start deviating substantially from the exact result, while for larger couplings we are well in the non-perturbative regime in which such expansions are useless. Throughout, the calculations based on the use of the derived effective actions converge to the exact ground energy even for relatively rough discretizations. From Table 1 we see that this holds even in the case of strongly interacting theory.

## 5 Energy estimators

One of the basic properties of physical systems is the internal energy given as the energy expectation value. Various ways of evaluating energy expectation values in PIMC simulations are based on different types of energy estimators – expressions whose thermal expectation values are used to calculate the discretized internal energy. The practical choice of the best estimator usually depends on the model and on the specific values of physical parameters. It is determined through a comparison of the properties of different estimators (e.g. variance, convergence, etc.), both in terms of calculated numerical values and their functional dependencies. Details of the derivation and characteristics of several energy estimators are given in [12] and [35]. The virial energy estimator turns out to be most advantageous due to its roughly constant Monte Carlo variance with the increase in the number of time slices. In this section we construct a hierarchy of virial energy estimators with improved convergence to the continuum limit as , while keeping the feature of constant variance of the naive virial estimator.

First we briefly review the standard derivation of the virial estimator [35, 36]. Starting from the formula , the straightforward generalization to discretized expressions reads:

 UN(β)=−∂logZN(β)∂β. (30)

The above partial derivative can be rewritten in terms of as . In order to simplify the calculation of this derivative, we remove the dependence of the path integral measure and of the kinetic term of the discretized expression for by simply rescaling . Now the differentiation over only affects the rescaled potential term in the exponent of the expression for . After the differentiation, we reinstall the original variables , and obtain:

 UN(β)=1ZN(12πε)N2∫dq1…dqNEviriale−SN. (31)

Here, is the standard virial energy estimator given as

 Evirial=1NN−1∑n=0(V(¯qn)+12¯qn,i∂iV(¯qn)). (32)

This estimator yields a typical convergence.

As shown in [36], improved convergence of expectation values follows from using the appropriate effective actions and energy estimators. As we have seen, virial estimators follow directly from the form of the discretized action. The level estimator is obtained by substituting for in this procedure. Writing the ideal effective action and its associated virial energy estimator as

 S∗N = SN+∞∑p=2σ(p) (33) E∗virial = Evirial+∞∑p=2e(p) , (34)

where and are corresponding contributions proportional to . Each is a sum over all the time-slices, i.e. . The same relation holds between and . The outlined procedure now gives the following simple connection between ideal action and estimator:

 e(p)n=1T(p+12¯qn,i∂i)σ(p)n . (35)

The analytically derived improvement in convergence has been verified by a series of Monte Carlo simulations that have been performed on a two particle two dimensional system in a quartic potential (22) for a variety of different values of coupling constant and inverse temperature. Typical results are presented in Figures 5 and 6. Figure 5 shows convergence of discretized energy expectation values for different levels . As level increases, we see that the continuum limit is approached faster, with ever smaller values of . Numerical results conform precisely to the analytically derived increase in convergence as can be seen from the log-log plot of the deviations of discretized values from the continuum limit shown in Figure 6.

## 6 Conclusions

We have presented a derivation of discretized effective actions and energy estimators that lead to substantial, systematic speedup of numerical procedures for the calculation of path integrals of a generic many-particle non-relativistic theory. The derived speedup holds for all path integrals – for transition amplitudes, partition functions, expectation values, as well as for calculations of energy levels. The obtained analytical results have been numerically verified through simulations of path integrals for multi-particle model with quartic coupling. In the paper we present explicit expressions for the effective actions up to level . The developed calculation scheme has been completed and is ready for applications to relevant problems in condensed matter physics. Further analytical work will focus on the generalization of the outlined scheme to more complex bosonic and fermionic systems, e.g. field theories.

## Acknowledgements

This work was supported in part by the Ministry of Science of the Republic of Serbia, under project No. OI141035, and the European Commission under EU Centre of Excellence grant CX-CMCS. The numerical results were obtained on the AEGIS e-Infrastructure, which is supported in part by FP6 projects EGEE-II and SEE-GRID-2.

## Appendix A: Integrals of propagators

In this appendix we present the values of all integrals necessary for calculations of effective actions at levels .

 ∫ε2−ε2ds Δ(s,s)k sn−2k={0n odd2−1−n ε1−k+n B(1+k,1−2k+n2)n even ,

where is the Euler beta function. All the multiple integrals needed have the same boundaries of integration as the preceding case.

 ∫dsds′ Δ(s,s′)=ε312 ∫dsds′ Δ(s,s′)3=ε5560 ∫dsds′ Δ(s,s′)Δ(s′,s′)=ε460 ∫dsds′ Δ(s,s′)Δ(s,s′)=ε490 ∫dsds′ s2 Δ(s,s′)=ε5240 ∫dsds′ ss′ Δ(s,s′)=ε5720 ∫dsds′ ss′3 Δ(s,s′)=ε76720 ∫dsds′ s2s′2 Δ(s,s′)=ε74032 ∫dsds′ ss′ Δ(s,s′)2=ε65040 ∫dsds′ s4 Δ(s,s′)=ε72240 ∫dsds′ s′2 Δ(s,s′)2=ε62520 ∫dsds′ s2 Δ(s′,s′)Δ(s,s′)=ε61260 ∫dsds′ s2 Δ(s,s)Δ(s,s′)=ε61680 ∫dsds′ ss′ Δ(s,s′)Δ(s′,s′)=ε65040 ∫dsds′ Δ(s,s)Δ(s,s′)Δ(s′,s′)=17ε55040 ∫dsds′ Δ(s,s)Δ(s,s′)2=ε5420 ∫dsds′ Δ(s,s′)Δ(s,s)2=ε5280 ∫dsds′ds′′ Δ(s,s′′)Δ(s′,s′′)=ε580

## Appendix B: Effective actions for levels p≤5

The ideal effective action for a general non-relativistic multi-particle system in the mid-point prescription is given in (33) as a sum of terms containing all contributions of order . As we have seen, is the sum over all the time-slices, i.e. , where is the contribution of time-slice . In this appendix we give a list of the explicit expressions for up to in shorthand notation in which and . The expressions for higher levels can be found on our web site [26].

 σ(2) = ε212∂2V+εδiδj24∂2ijV σ(3) = −ε324∂iV∂iV+ε3240∂4V+ε2δiδj480∂2ij∂2V+εδiδjδkδl1920∂4ijklV σ(4) = ε46720∂6V−ε4120∂iV∂i∂2V−ε4360∂2ijV∂2ijV−ε3δiδj480∂kV∂3kijV+ + ε3δiδj13440∂2ij∂4V−ε3δiδj1440∂2ikV∂2kjV+ε2δiδjδkδl53760∂4ijkl∂2V+ + εδiδjδkδlδmδn322560∂6ijklmnV σ(5) = ε5241920∂8V−ε51680∂2ijV∂2ij∂2V−17ε540320∂i∂2V∂i∂2V− − ε52240∂iV∂i∂4V−ε56720∂3ijkV∂3ijkV+ε5240∂iV∂jV∂2ijV+ + ε4δiδj483840∂2ij∂6V−ε4δiδj6720∂kV∂3ijk∂2V−ε4δiδj10080∂2ikV∂2kj∂2V− − ε4δiδj10080∂2klV∂4ijklV−ε4δiδj5040∂k∂2V∂3ijkV−ε4δiδj20160∂3iklV∂3kljV+ + ε3δiδjδkδl1935360∂4ijkl∂4V−ε3δiδjδkδl53760∂mV∂5mijklV− − ε3δiδjδkδl40320∂2imV∂4mjklV−ε3δiδjδkδl32256∂3ijmV∂3mklV+ + ε2δiδjδkδlδmδn11612160∂6ijklmn∂2V+εδiδjδkδlδmδnδpδq92897280∂8ijklmnpqV

### References

1. R. P. Feynman, Rev. Mod. Phys. 20, 367 (1948).
2. R. P. Feynman, A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York, 1965.
3. C. Itzykson and J.-B. Zuber, Quantum Field Theory, McGraw-Hill, New York, 1980.
4. S. Coleman, Aspects of Symmetry, Cambridge University Press, 1985.
5. G. Parisi, Statistical Field Theory, Academic Press, 1988.
6. C. Itzykson and J-M. Drouffe, Statistical Field Theory, Cambridge Universityt Press, 1991
7. B. Baaquie, Quantum Finance, Cambridge University Press, 2004.
8. H. Kleinert, Physica A 311, 536 (2002).
9. H. Kleinert, Physica A 312, 217 (2002).
10. F. Steiner, C. Grosche, A Table of Feynman Path Integrals, Springer, New York, 1997.
11. H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets, 3rd edition, World Scientific, Singapore, 2004.
12. D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
13. M. Suzuki, Phys. Lett. A 146, 319 (1990).
14. A. Bogojević, A. Balaž, A. Belić, Phys. Rev. B 72, 064302 (2005).
15. M. Takahashi, M. Imada, J. Phys. Soc. Jpn. 53, 3765 (1984).
16. X. P. Li, J. Q. Broughton, J. Chem. Phys. 86, 5094 (1987).
17. H. De Raedt, B. De Raedt, Phys. Rev. A 28, 3575 (1983).
18. C. Predescu, Phys. Rev. E 69, 056701 (2004).
19. C. Predescu, Phys. Rev. E 71, 045701 (2005).
20. C. Predescu, J. Phys. Chem B 110, 667 (2006).
21. T. M. Yamamoto, J. Chem. Phys. 123, 104101 (2005).
22. S. Paulin, A. Alastuey, T. Dauxois, J. Stat. Phys. 128, 1391 (2007).
23. A. Bogojević, A. Balaž, A. Belić, Phys. Rev. E 72, 036128 (2005).
24. A. Bogojević, A. Balaž, A. Belić, Phys. Lett. A 344, 84 (2005).
25. A. Bogojević, A. Balaž, A. Belić, Phys. Rev. Lett. 94, 180403 (2005).
26. http://scl.phy.bg.ac.yu/speedup/
27. R. P. Feynman, H. Kleinert, Phys. Rev. A 34, 5080 (1986).
28. H. Kleinert, Phys. Lett. B 280, 251 (1992).
29. H. Kleinert, Phys. Lett. A 173, 332 (1993).
30. H. Kleinert, W. Kürzinger, A. Pelster, J. Phys. A 31, 8307 (1998).
31. F. Weissbach, A. Pelster, B. Hamprecht, Phys. Rev. E 66, 036129 (2002).
32. A. Pelster, H. Kleinert, M. Schanz, Phys. Rev. E 67, 016604 (2003).
33. S. F. Brandt, A. Pelster, J. Math. Phys. 46, 112105 (2005).
34. D. Stojiljković, A. Bogojević, A. Balaž, Phys. Lett. A 360, 217 (2006).
35. W. Janke, T. Sauer, J. Chem. Phys. 107 5821 (1997).
36. J. Grujić, A. Bogojević, A. Balaž, Phys. Lett. A 360, 205 (2006).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters