A new method for the calculation of functional and path integrals

This paper addresses a disconnect between the pivotal role of functional (path) integrals in modern theories, such as quantum mechanics and statistical thermodynamics, and the currently limited ability to perform the actual calculation. We present a new method for calculating functional integrals, based on a finite-element formulation, which solves all limitations of existing methods. This approach is far more robust, versatile, and powerful than the prevailing methods, thus allowing for more sophisticated computations and the study of problems that could not previously be tackled. Importantly, existing procedures, element libraries and shape functions, which have been developed throughout the years in the context of engineering analysis and partial differential equations, may be directly employed for this purpose.

Functional (or path) integrals, are ubiquitous in a wide range of physical and mathematical problems, ranging from quantum mechanics to statistical thermodynamics through biology, chemistry, engineering and finance [1][2][3][4][5][6][7][8][9].Functional integrals generalize the notion of integration over vector spaces to integration over function spaces.Much like vector integration, f (v) dv, where the integral is the sum of all volume elements in the integration space each weighted by an integrand function, functional integration, g[u]Du, is the summation process over all admissible functions each weighted by an integrand functional.
Currently, the prevailing method for evaluating path-integrals is a "slicing method" which involves a naive discretization of the spatial (or temporal) space followed by summation over the function values at the discrete points [1][2][3].In that approach, for example, derivatives of the state function with respect to the independent coordinates are approximated using finite-differences.The "slicing method" is straight forward and useful, however it is limited to 1-D systems in time or space, or to very simple and regular spatial geometries, such as rectangles.Further, even in 1-D, the application of certain boundary conditions, as well as constraints inside the domain, can be cumbersome.Thus, there is a clear need in a new, more sophisticated, method for the calculation of functional integrals.
In this paper, we propose a new approach for the calculation of functional integrals that is based on the finite-element (FE) formulation rather than the aforementioned naïve discretization.The FE method is a numerical approximation method predominantly used for solving partial differential equations (PDEs) in all fields of engineering, ranging from stress analysis and heat transfer in structures to electromagnetic scattering of objects and the design of photonic crystals [10][11][12][13][14]. Driven by the growing need to tackle more complicated engineering problems, the FE method has become the standard numerical tool for engineering design and mechanical analysis, thanks to its generality, robustness and versatility.Still, the use of the FE approach for evaluating functional integrals has been largely overlooked, practically limiting the computation of path integrals to 1-D or to very simple geometries.
As we show below, applying the FE approach to functional integrals requires some technical care; however, the powerful formulation opens the door for more sophisticated computations and for the study of problems that could not previously be tackled.In that sense, it may be reminiscent of the revolution brought by the FE approach to partial differential equations, which has enabled solving complex engineering problems with complicated geometries and all types of boundary conditions or constraints.Further, although the method we propose is new, the underlying foundations are mature and well-established: Countless papers and textbooks have been published on the theory of FE [10,11,[15][16][17][18], and highly sophisticated computational schemes and software have been developed; all may readily be repurposed for our needs.For example, one may directly use the large established libraries of finite elements and associated shape-functions [19], or apply well-established meshing procedures and related software [20].Perhaps the most important advantage that the FE formulation introduces to the computation of functional integrals is the preservation of the state functions as functional entities; this is unlike the existing approaches where the state functions are reduced to a finite set of discrete points.This property of the method allows for a more rigorous treatment and valuable insights.Moreover, even in simple cases, where the formulation does not provide fundamentally new results, we may still find ourselves appreciating the refreshing interpretations the method gives, or by the words of Richard P.
Feynman there is a pleasure in recognizing old things from a new point of view [2].
The Method.-Consider a stochastic system whose (micro)state is described by the real valued function u(x) ∈ V , where V is the space of all possible states, or admissible functions, of the system, and x is a point in a parameter space Ω.At the moment don't make any assumption on the class of continuity of V .Let the functional p[u], with p : V → R, be the probability density corresponding to state u(x); and let g[u] be some functional, g : V → R, dependent on the system's state.The average of the state-dependent functional g[u] is defined as the sum of g[u] over all possible states of the system weighed by the respective probability p[u]Du, hence expressed as a functional integral The mathematical rigor for such operation is quite subtle, but the concept is well established, e.g., the Feynman path integral [1][2][3]; in that context, equation ( 1) should be regarded as a formal way to express the functional averaging process.
We hereby present a new method for the calculation of functional integrals, based on a FE formulation.As a first step, we approximate the function space V by introducing a subspace V h ⊂ V .The superscript h is called the mesh parameter and it implies that the functions u h ∈ V h are associated with a mesh, or a discretization, of the domain Ω.The mesh parameter, h, is a measure for the size of the largest element in the mesh, therefore Here η is the set of open nodes, i.e. nodes where u is variable, and η u is the set of closed nodes, i.e. nodes where u is prescribed.The functions φ A (x) are called FE shape functions; these are continuous functions that are related to the mesh, in the sense that each function φ A (x) gets the value of one at node A and vanishes in all elements that don't contain that node.
It is therefore apparent that the coefficients d A and ūA are the values of u at the respective node.It is shown in FE theory that as h → 0 then u h → u and the approximation error is given by u − u h = Ch α , where α is the rate of convergence which depends on the type of FE shape functions used [10].
Once we are able to represent u by a finite number of degrees of freedom (DOFs), we may substitute the approximation u h into the probability density and obtain p h (d).Note that p h is a function of the variables d = {d A }, rather than a functional, and it is normalized such that p h (d) dd = 1.Similarly, we express g h (d) = g[u h ] as a function of d.Since d uniquely determines u h and g h , then the problem of finding g[u h ] reduces to computing g h (d) , where Above, N is the number of DOFs (the cardinal number of η) and R N is the space of all real vectors of length N .Equation ( 3) is thus the FE approximation of the functional integral (1), and the finite-dimensional integration can be carried out analytically (where possible) or numerically.The Markov-Chain Monte-Carlo method is especially appropriate for this task, because its rate of convergence is independent on N and it is particularly suitable for finding the statistical moments of complicated distributions [21].
Note that the method described here can be readily generalized for vector functions u = u(x) by writing equation ( 2) separately for each component of u.Further, the domain Ω may be a simple 1-D interval, as in the case of the path integral in quantum mechanics where the coordinate x above represents time, or be a multi-dimensional domain of any geometry, such as in the case of statistical thermodynamics of a 3-D body.Unfortunately, the "slicing method", commonly adopted for computing path integrals, cannot be applied to the latter.In terms of the proposed method, however, the only difference between these two cases is the use of different finite elements and corresponding shape functions; while in the 1-D case the elements are lines (or 1-D segments) and the shape functions are described in terms of one coordinate, in 2-D the elements have a 2-D geometry, such as triangles or quadrilaterals, and the corresponding shape functions are described using two coordinates.
Similarly, if the domain Ω is three-dimensional, 3-D elements are used, etc.It is noted that the use of a non-uniform mesh, where the domain Ω is divided into elements of different sizes, is a standard practice of the FE method as illustrated in figure 1.This allows, for example, to use of a finer mesh in regions where high accuracy is needed.This feature is another important attribute of the versatile and powerful finite-element formulation.Finally, the physics of the problem dictates the number of DOFs at each node of the element.This is exemplified in the two examples below, where the first example involves one DOF at each node, while for the second example two DOFs are used at each node.
Example I: string.-Consider a string of length L and uniform tension σ with both ends held fixed at a horizontal level.A lateral force f (x) is distributed along the string, and the entire system is submerged in a heat reservoir of temperature T .Let u(x) ∈ V describe the transverse displacement of the string at x ∈ Ω = [0, L] and regard u as the state of the system.We would like to find, for example, the average state of the system.The space is the set of all square-integrable functions over Ω with square-integrable first derivatives (Sobolev space) that admit the boundary conditions FIG. 1. Example of a non-uniform mesh over a 2D domain in the shape of a cat's silhouette.The mesh was generated using the open-source Gmsh [20].
u(0) = u(L) = ū = 0.The probability density corresponding to micro-state u is [22] p where the partition function, Z, is a normalization constant, β = (k B T ) −1 , and E[u] is the energy functional with u ,x = ∂u/∂x .Accordingly, the average state of the system is given by the functional integral Following the method described above, we write the finite-element approximation as where, for simplicity, φ A (x) are linear FE shape functions [10], and thus u is approximated by u h in a continuous piece-wise-linear manner.Next, we substitute the FE approximation u h into the energy functional E[u] so it becomes a function E h (d) = E[u h ] of the variable d.This energy integral is calculated at the element level.Thus, define the vectors φ(x) = {φ e 1 (x), φ e 2 (x)} T of the element shape functions and d e = {d e 1 , d e 2 } T of the element state values, such that the state in each element is given by u h = φ T d e .Then, rewrite the integral (5) as a sum of integrals over Ω e ⊂ Ω, the domain of the e-th element.The energy of that element is then where f e = f e 1 , f e 2 T is a vector whose components are the values of f at the element nodes, such that f ≈ f h = φ T f e describes the force inside the element.The symmetric matrices k e = σ Ω e φ ,x φ T ,x dΩ and m e = Ω e φφ T dΩ are respectively called the element stiffness and mass matrices and they are calculated in each element separately; however, for many elements, including the 1D element considered here, a formula for these matrices is readily found in the literature [10,19].Next, define the matrix K, vector v and scalar S such that the expression for the global energy becomes Here, K is the global stiffness matrix, the vectors v and F are related to the contribution of prescribed displacements and of external loads, respectively, and the scalar S corresponds solely to the contribution of boundary condition to the energy.This quadratic energy may now be substituted into equation ( 4) so it becomes an off-centered Gaussian and thus d may be calculated analytically from equation ( 3) by considering the particular case of g h (d) = d .
Example II: Beam.-Consider a stochastic system that is modeled by an Euler-Bernoulli beam of length L and uniform bending stiffness K B with one end fixed and the other free.A lateral load f (x) is distributed along the beam, and the entire system is submerged in a heat reservoir of temperature T .Let u(x) ∈ V describe the transverse deflection of the beam at x ∈ Ω = [0, L] and regard u as the state of the system.The space is the set of all square-integrable functions over Ω with square-integrable first and second derivatives that admit the fixed boundary conditions u(0) = 0, u x (0) = 0. Similar to the previous example, we want to find the average state of the system.The statistical distribution is given by equation ( 4), with the energy functional The fundamental difference compared to the previous example is that for the energy (10) to be well defined, we demand stronger requirements on the class of continuity of u, namely that V ⊂ H 2 .Accordingly, u h ∈ V h must satisfy these continuity requirements.To this end, we approximate u using the Hermite cubic shape functions φ (Ai) where i = 1, 2 [10].
These shape functions and their derivative vanish everywhere, except in the elements that share node A. Moreover, at node A, the shape functions satisfy φ (A1) = φ (A2),x = 1 and φ (A2) = φ (A1),x = 0.This property allows dictating separately the displacements (u) and rotations (u ,x ) at the element nodes, thus Here, N ndof is the number of nodal DOFs (in our case, N ndof = 2).The values of u at the nodes are d (A1) and the values of u ,x at the nodes are d (A2) , the sets η (i) and η T of the element state values, such that the state in each element is given by the cubic function u h = φ T d e .The energy expression for the element is the same as in equation ( 8), other than k e = K B Ω e φ ,xx φ T ,xx dΩ; therefore the global energy also has the quadratic form of equation ( 9) and d may be calculated analytically from equation (3).
Adhesion of elastic body to a rigid substrate: a numerical example.-Inwhat follows, we present numerical results obtained using the proposed finite-element formulation.The model considered is prototypical to phenomena such as detachment of biological cells, peeling of a thin film from a substrate, etc., and demonstrates how the formulation can be conveniently applied to complex systems composed of coupled linear and non-linear elements.We emphasize that while the model may be suitable for describing real phenomena, such as those mentioned above, it is presented here merely for demonstrating the proposed method; thus justification of the model and its assumptions are not further discussed.Consider an Euler-Bernoulli beam of the sort described in the previous example, but instead of a lateral distributed force, the non-fixed end of the beam is supported at a height ū.Thus, The beam is adhered to a rigid substrate as illustrated in the inset of figure 2. The adhesion is modeled by a set of N bonds connected at points x A ∈ Ω (A = 1, ..., N ) along the beam.
When the A-th bond is connected it acts as a linear spring of stiffness k A , and when it is broken it exerts no force.Hence we use the following potential function to describe the adhesion Here U A is constant of units length that describes the broken state potential in terms of elongation of the spring.Note that due to the stochastic nature of our system, each bond may break and reconnect randomly.Accordingly, the state of each bond, either connected or broken, is identified by a two-state spin variable.A similar, yet simpler, adhesion-decohesion model was introduced and discussed by Florio et al. [23].There, it was suggested to introduce a single N -state spin variable, ξ, for the entire array of bonds.This is based on the assumption that due to the one-sided decohesion process we have ξ connected bond at the fixed-end-side and N − ξ broken bonds at the supported-end-side.The potential energy of the system is therefore The next step is to approximate u using Hermite cubic shape functions as was done in the previous example.In principle, one may use any mesh as long as it has nodes at all the {x A } points where the beam is attached to the substrate through a breakable bond.For simplicity, we consider here a mesh with N unknowns located where the springs are connected.The approximate energy function is then where K, v and S are defined as usual, and K conc. is a matrix full of zeros except for kU Discussion and conclusions.-Wepresented a new method for calculating functional integrals based on finite-elements formulation.The proposed method is far more robust, versatile and powerful than any prevailing method, as it allows the calculation of functional integrals over any domain subjected to any boundary conditions or constraints, while not limited to 1-D domains like the "slicing method".Due to the nature of the discretization, a finer mesh may be used in regions where high accuracy is needed.Moreover, by employing the FE formulation, the functional identity of the state-function is naturally maintained throughout the calculation, enabling insightful perspectives even in 1-D.Just as importantly, existing finite-element routines, elements libraries and shape functions, which have been developed throughout the years for solving PDEs, can be directly employed for as in equation (2) but they include only nodes with open or closed i-th DOF.Define the vector φ(x) = φ e (11) (x), φ e (12) (x), φ e (21) (x), φ e (22) (x) T of the element shape functions and the vector d e = d e (11) , d e (12) , d e (21) , d e

FIG. 2 .
FIG. 2. Relation between the end-displacement ū and (a) the mean force f , or (b) the mean number of attached bonds ξ , for various values of β.Inset: schematic illustration of the model.

2
at ξ entries along the diagonal corresponding to connected bonds.The statistical distribution of the system is p[u; ξ] = exp(−βE[u; ξ])/Z, and the discretized version is the well-studied Gaussian distribution.Once we have p h (d, ξ) we essentially know everything about the system, and we can calculate some interesting statistical properties and study how they are influenced by temperature.For example, figure 2(a) shows the effect of temperature on the force-displacement relation of the mean force, f = −(∂Z h ∂ ū)/β, applied by the support with respect to the prescribed displacement ū.We may also calculate ξ , the mean number of connected bonds, as a function of ū, as shown in figure 2(b).These results are given in non-dimensional form after the energy was rescaled by E 0 = K B U 2 /L 3 and lengths were rescaled by U ; the calculations were carried out with the values: kU 2 /E 0 = 5, β = {15, 10, 6, 4, 3, 2, 1}E 0 and N = 6.