A Monte Carlo Approach to the Worldline Formalism in Curved Space

We present a numerical method to evaluate worldline (WL) path integrals defined on a curved Euclidean space, sampled with Monte Carlo (MC) techniques. In particular, we adopt an algorithm known as YLOOPS with a slight modification due to the introduction of a quadratic term which has the function of stabilizing and speeding up the convergence. Our method, as the perturbative counterparts, treats the non-trivial measure and deviation of the kinetic term from flat, as interaction terms. Moreover, the numerical discretization adopted in the present WLMC is realized with respect to the proper time of the associated bosonic point-particle, hence such procedure may be seen as an analogue of the time-slicing (TS) discretization already introduced to construct quantum path integrals in curved space. As a result, a TS counter-term is taken into account during the computation. The method is tested against existing analytic calculations of the heat kernel for a free bosonic point-particle in a D-dimensional maximally symmetric space.


Introduction
The worldline formalism applied to Quantum Field Theory is a powerful tool to compute relevant physical quantities. The guidelines of such method were already introduced in the 1950 by Feynman, who proposed a quantum mechanical path integral model representing the dressed scalar propagator in scalar QED [1]. In the 80's the method started to be systematically used as an alternative to standard secondquantized approaches to QFT, in particular for the computation of anomalies [2][3][4][5][6], effective actions and Feynman diagrams [7][8][9]. This formalism was then extended to various quantum field theories in curved spacetime (see Ref. [10] for a review) and several applications have so far been considered such as, oneloop effective actions [11][12][13], the one-loop graviton photon mixing in an electromagnetic field [14], gravitational corrections to Euler-Heisenberg lagrangians [15], worldline representations of quantum gravity [16,17], supergravity [18] and higher spin field theories [19,20], and simplified methods for anomaly computations [21,22], just to name a few.
The effort to move towards non-perturbative worldline calculations led to the development of numerical methods for the evaluation of quantum mechanical path integrals in flat space. The need for such a generalization was in the several interesting applications that rely on non-perturbative physics, such as the chiral symmetry breaking. The main issue was the numerical generation of the worldlines in an Euclidean space according to their relative probability distribution. For such a purpose, Gies and Langfeld [23] proposed a first numerical algorithm, adopting a Monte Carlo sampling of the coordinate space as an improvement of a previous work by Nieuwenhuis and Tjon [24]. Such methods are now known as Worldline Monte Carlo (WLMC) and have been used for several application in QFT, such as the Casimir effect [25], Schwinger pair production in inhomogeneous fields [26], quantum effective actions [27] and strongly-coupled, large-N fermion models [28]. In these works, all the proposed algorithms are based on a Monte Carlo sampling which models the free Brownian motion. Thus, instead of using the full (kinetic + potential) action as a weight to select the points on the worldline, only the kinetic part is used. As explained in the recent work by Edwards et al. [29], the reason for such a choice is universality, which then allows the application of the method to a variety of different cases, regardless of the specific potential involved. Moreover, therein they propose an efficient algorithm, called YLOOPS, to generate closed worldlines around a point, 1 which is a modified version of the alghoritm VLOOPS originally designed in [25]. Below, due to its efficiency and easiness of application, we find it convenient to adopt the YLOOPS algorithm for our calculation: more precisely, we present a slight modification of such algorithm, where a fictitious quadratic term is inserted in the kinetic term, in order to make the evaluation of the numerical path integral faster and more stable.
The main goal of the present manuscript, is to extend WLMC techniques to curved spaces, i.e. to non-linear sigma models representing the propagation of a scalar particle in curved space. This nonperturbative numerical method would presumably help tackling a large class of problems in quantum field theory in curved space such as, for example, the Casimir effect and the Schwinger effect on a non trivial background, and the computation of effective actions and anomalies in curved space. In order to test our construction we restrict ourselves to maximally symmetric spaces. In such a context indeed, recently it was shown that the heat kernel expansion for a scalar particle [30,31] (later generalized to a spinor particle [32]) can be efficiently reproduced using an effective model with a flat kinetic term, where the effects of the curvatures are taken into account through a suitable effective potential. Therefore, as we show, such a model can be easily simulated with the conventional (flat space) WLMC methods, mentioned above. Moreover, this constitutes a perfect benchmark to verify our extension of the WLMC techniques to a genuine non-linear sigma model representing a particle in curved space.
The paper is organized as follows. In Section 2 we review the Worldline Monte Carlo theory in flat space, providing an example of the calculation of the propagator associated to a 4-dimensional harmonic oscillator, to get familiar with the method. In Section 3 we build the WLMC setup in curved space and propose our strategy to compute numerical path integrals on non-trivial backgrounds. In Section 4 we describe the theoretical basis of the case study which we will use to test our method, i.e. the diagonal part of the heat kernel of a free scalar point particle moving on a D-sphere and D-hyperboloid, i.e. maximally symmetric spaces. The worldline realization of such quantity is given in terms of a quantum mechanical path integrals with periodic boundary conditions (defining closed trajectories or loops). After that, we will report our numerical simulations, studying their convergence with respect to the discretization parameters and the curvature of the sphere. A possible extension to the case of open trajectories is discussed in Section 5. In Section 6 we sum up our conclusions, with possible outlook for future applications. Finally, in Appendix A we report the details of the YLOOPS algorithm generalized to the case of a non-vanishing quadratic term, whereas in Appendix B we see the effect of inserting such mass term in the WL sampling discussed in Section 2.

Worldline Monte Carlo in flat space
In order to introduce our computation, let us first review the main ingredients of worldline Monte Carlo in flat space. Let us consider a D-dimensional flat space heat kernel computed through an Euclidean quantum mechanical path integral with and Expression (1) can be seen as the heat kernel of a unit mass non-relativistic point-particle x(τ ) in Ddimensional Euclidean space, or conversely the (Schwinger integrand of the) propagator of a bosonic relativistic particle, where the affine parameter τ is the proper time. Its Hamiltonian reads The formal expression (1) involves the sum over all the possible worldlines joining x and x in time β, weighted by the action (3). To provide a numerical realization of the worldline path integral (1), some comments are needed. First of all, it is not possible to span numerically the entire configuration space, hence a selection on the worldlines must be taken into account: let us denote N W L the number of worldlines which are considered. Secondly, each worldline has to be discretized: this is usually done with respect to the affine parameter, taking a number N of points per worldline. As pointed out in [29], here discretization is not realized on the points x(τ ) of the manifold, 2 but directly on the affine parameter τ . To fulfill the above two requirements, what is usually done in WLMC approaches [23,[25][26][27]29] is to compute worldline averages like x(β)=x with S KIN [x] and S P OT [x] being the kinetic and potential terms in (3) respectively. In this way, the exponential e −S KIN characterizing each worldline can be used as a weight function to build the worldlines used for the calculation. To be more precise, worldlines are constructed point by point using Monte Carlo algorithms: each trajectory produced, then, satisfies a Monte Carlo-selection criterion which allows to consider all the trajectories obtained as a set of almost equally meaningful trajectories to simulate the quantum propagation of the point-particle in spacetime. The aforementioned choice of the WL weight function is the most popular in literature, as it has an evident feature of universality, in fact in such case the selection of the worldline is independent of the model considered. Once the configuration space is sampled with the points x where all the worldlines start at x and end at x . At this point, expression (5) can be approximated by the arithmetic average for which, clearly, a good estimate of the path integral is obtained when N and N W L are large enough. Expression (7) is the main prescription for our WLMC computations. Note that, unlike (5), the latter does not occur in the form of a weighted average. However, the weighing procedure takes place beforehand, upon the selection of the different worldlines present in the ensemble. As mentioned before, different algorithms have been developed in order to generate worldlines. However, historically, the most commonly used algorithm is the VLOOPS developed by Gies [23]. Here we find it convenient to use the YLOOPS routine developed by Edwards et al. [29], with a slight modification due to the inclusion of a regularizing quadratic contribution in the kinetic term that will be taken into account for the computations in curved space. Such algorithm is developed to diagonalize the discretized version of the flat kinetic term in (3), having vanishing Dirichlet boundary conditions x 1 = x N = 0. The details are discussed in Appendix A.

A simple flat-space test
In order to conclude this preliminary section let us review a simple example where the WLMC setup can be applied, i.e. the propagator of a bosonic harmonic oscillator with null endpoints in a D-dimensional Euclidean space. Such quantity can be written as where we have assumed ω = 1 = m. The analytical result is well known and reads On the other hand, the propagator can be written as where we have multiplied and divided by the free path integral, which provides the Feynman factor x(0)=0 . Expression (10) thus reads and, combining (11) with (9), we obtain which is accessible to WLMC calculation. Indeed, in Fig. 1 we report a WLMC calculation of G(β) in four Euclidean dimensions, with N = 1000 points per loop and N W L = 1000 worldlines. It shows a satisfactory agreement with the analytical result in (12). In order to remove correlation between the data, this calculation-as well as those presented henceforth-has been performed using a different set of worldlines for each β-value [29]. As we have seen, WLMC in flat space is relatively easy to implement and, even for not very large values of N and N W L , it faithfully reproduces the theoretical curve. We will show that the flat space WLMC setup can be fruitfully exploited also for numerical applications in curved space problems: the price to pay is the introduction of further potential-like terms that incorporates curvature effects.

Worldline Monte Carlo in curved space
Similarly to what we did in the previous section for the flat space case, let us consider the classical Hamiltonian of a D-dimensional bosonic scalar point-particle in Euclidean curved space At the quantum level the Einstein invariant Hamiltonian operator reads (see [10] for a complete treatment of particle path integrals in curved space) with g(x) = det g µν (x), and the heat kernel associated to (14) can be expressed as the following particle path integral where the Einstein-invariant formal measure reads and the resulting particle action is a one-dimensional non-linear sigma model. With respect to the flat space case, there are three main differences: the root factor in (16), the spacetime dependent metric tensor g µν (x) and the counter-term V CT (x) appearing in (17). A few comments are in order. In perturbative computations, one needs a regularization scheme to treat divergent Feynman diagrams and V CT (x) is the suitable, finite, counterterm needed, which depends on the adopted regularization scheme. On the other hand, in non-perturbative computations, no divergences are expected to occur since these Feynman path integrals represent quantum mechanical transition amplitudes which are thus finite. However, in curved spaces, there are ordering ambiguities as the kinetic hamiltonian mixes coordinates and momenta. Yet, a well known scheme that allows to obtain Feynman path integrals from the associated quantum mechanical transition amplitudes exists: it is the Time Slicing approach, which requires to Weyl-order the Einstein-inviariant hamiltonian, in order to employ the so-called "mid-point rule" in the costruction of the path integral. Such reordering produces an order-potential, V CT (x), which is the same potential needed at the perturbative level, since the Time Slicing also encodes a regularization procedure. For historical reasons we will refer to such potential as the "counterterm potential". Finally, the root factor in (16) renders the formal measure Einstein invariant: for perturbative calculations it is often customary to exponentiate the » g(x) in terms of a particle path integrals over Lee-Yang ghost fields. However, for our non-perturbative, numeric computation we find more convenient to keep it as an external factor, evaluated at each point of the various worldlines involved in the sums. Now, the key idea in order to perform a numerical WLMC average in curved space (like the one in (5)) is to extract the flat space kinetic term from the curved one and use it to sample the worldlines, and treat the remaining part as a potential contribution, so that it is possible to bring the curved space problem back to a flat space one-provided one takes into account the suitable counterterm and measure factors. Namely, after a flat space Monte Carlo sampling of the worldlines, the path integral average will thus be given by where and is the middle point between adjacent points of the worldline s.
In particular, the first potential term appearing in (20) is the one arising from (18) and reads It is easy to verify that, as discussed in the Appendix A, definitions (20) and (22) respect the backward convention on the first derivative. Focusing on the counter-term potential V CT , it ought to be observed that the proper time discretization of the path integral which we performed is equivalent to the Time Slicing (TS) procedure adopted to regularize particle path integrals in curved space from first principles [10]. Thus, the discretized counter-term is 4 Now we have all the ingredients to test our setup.
In the following section we consider a case study which comes particularly handy for our purpose: it is the study of the free heat kernel of a scalar point particle constrained on a D-sphere. Recently, in Ref. [30] (see also [31]), it was shown that, using Riemann normal coordinates (RNC), it is possible to map the problem of computing the heat kernel on spheres, to the evaluation of a flat space heat kernel with a suitable potential which takes into account the curvature effects. Moreover, as found in [21], for a maximally symmetric geometry, the metric tensor, in RNC, can be neatly written in closed form as the flat space metric plus a curvature-dependent contribution. Hence, our numerical problem can be studied both from a purely flat perspective, by mapping it into the effective flat space model of [30]-thus using the conventional flat space Worldline Monte Carlo reviewed in Section 2-, and using the curved space WLMC setup discussed in Section 3. In other words, we use the effective flat space model as a benchmark test for our curved construction.
4 Case study: free scalar heat kernel on maximally symmetric spaces Let us briefly review the model [30] which we are going to use to implement our numerical WLMC method in curved space. The Hamiltonian operator of a free scalar point-particle in curved space can be used to build its heat kernel satisfying the operatorial equation The matrix elements of the kernel operator satisfy the heat equation DefiningK and using RNC, the curved space heat equation (28) can be turned into the flat space heat equation if the curved space is a space with maximally symmetry (a D-sphere for definiteness), for which In such a case we have Moreover, the effective potential appearing in (30) can be expressed in closed form as and can be used in the associated path integral representation of the problem Hence, the heat kernel for a particle on a maximally symmetric space can be obtained through a linear sigma model with an effective potential that encodes the curvature effects. We use such effective representation of the heat kernel to test our numerical method.
Effective  Specifically, we summarize in Table 1 the details of the two numerical calculations that we have performed. The first one, which we refer to as "Effective Potential Method" (EPM) makes use of the potential given in Eq. (33) to effectively reproduce the heat kernel of a scalar particle on a maximally symmetric space through a linear sigma model (flat metric). In the second calculation, we compute the same heat kernel making use of a "genuine" non-linear sigma model action, which amounts to consider the two potentials (22) and (23) within the WLMC technique-we refer to this case as the "Non-linear Sigma Model Method" (NSMM). Note that, for the EPM computation we have adopted the original YLOOPS algorithm, developed in [29], whereas for the NSMM we found it convenient to use a generalized version of YLOOPS, where a tiny quadratic term is included in the kinetic part of the action. Namely, Let us stress that such mass term is used solely at the stage of the construction of the worldline ensemble, to improve the convergence, and does not contribute to the potential term, which is left untouched. The details of the generalized YLOOPS algorithm are discussed in Appendix A, whereas its numerical outcomes are shown in Appendix B.
(a) (b) Figure 2: Calculations of the heat kernel of a free scalar particle in a D=4 sphere with endpoints x = x = 0 for various parameters (N W L , N ) both using the effective potential method (green data) with known potential and using the non-linear sigma model method (blue data). The red curve is the analytical and perturbative computation for small β performed in [30].
In Figure 2, we report the results of our calculations for a D = 4 sphere with M = 1. Figures 2(a) and 2(b) show a comparison between the data of the Non-linear Sigma Model Method (NSMM, blue) and the Effective Potential Method (EPM, green), together with the perturbative, analytical computation performed in [30], for couples of values (N W L , N ) = (100, 100) and (3000, 3000) respectively. As we can see, the EPM is considerably more precise than the NSMM, presumably because of the form of the associated potential: in the first case it is expressed in regular form (33), whilst in the second case by means of numerical derivatives (22), and this contributes to reduce precision. Thus, we take the green points as a benchmark to test the blue ones. We notice that the agreement between the two calculations gets better as the discretized approximation improves, i.e. when N W L and N become large. In fact, the absolute value of the maximum relative error passes from ∼ 35% with (N W L , N ) = (100, 100) ( fig.  2(a)) to ∼ 4% with (N W L , N ) = (3000, 3000) ( fig. 2(b)). In fact, it was noticed that WLMC calculations generally lose accuracy at large time values due to the numerical phenomenon of undersampling [29] of worldlines. A priori, it would seem that our results confirm this feature. However, in the present case the lose of accuracy is to be ascribed to minor efficiency in representing the derivative interactions. In the attempt to increase the efficiency of the code, one may think of using a higher order discrete derivative to insert into (22), like the two-point method or, say, the five-point method However, the specific form of the counter-term (23) is compatible only with a first order backward derivative, as it was derived [10] in this way, and numerical simulations confirmed that. What we observe is that, locally, the decrease of curvature improves the agreement between NSMM and EPM-the precision has gained roughly two orders of magnitude (the full-scale changes from 4% to 0.04%) passing from M = 1 to M = 0.1. This is due to the fact that, magnifying the radius of the sphere, the space in the vicinity of x = x = 0 gets closer to be flat, reducing the stronger effects which affect the precision of the NSMM.
In general, the different degree of precision between the two methods considered, can be ascribed to the fact that-as already observed in perturbative computations [30]-the flat effective model (EPM) takes into account of the curvature effects in a much more efficient way than the non-linear sigma model, as in the former the curvature effects are all encoded in the effective potential. In fact, in the non-linear sigma-model, derivative interactions do appear, which are less efficiently represented in the heat-kernel expansion. This can be understood, by rewriting the action (17) in terms of a rescaled time s = τ /β which yields in which, the potentials appear with an order-β 2 higher than the derivative terms. Thus, in a perturbative expansion about a fixed point, one needs higher order terms from the Taylor expansion of g µν to match the corresponding orders in the potential. On the other hand, the flat effective model-unlike the nonlinear sigma model-has, so far, been shown to work only in maximally symmetric backgrounds. Finally, we present the calculation of the free particle heat kernel constrained in a D = 4 maximally symmetric Anti-de Sitter space, i.e. a 4-hyperboloid. Following [30], the sectional curvature is negative M 2 < 0 and, defining |M | = √ −M 2 , the metric is written as withf All other quantities appearing in (32) are defined in the same way as before, whereas in the potential (33) we have the replacement M → iM . Figure 5 shows a satisfying agreement between the curved method and the flat one for (N W L , N ) = (1000, 1000) for a range of β-values which is now extended to 10. The reason for such agreement at large β may be reasonably ascribed to the exponential suppression of the heat kernel on this scale, implying a sensitive reduction of Monte Carlo fluctuations too. Figure 5: Calculations of the heat kernel for the free scalar particle on a constant curvature-AdS space with D = 4 and |M | = 1.

An extension to open worldlines
So far we have seen an efficient way to numerically compute particle path integrals in curved space, with a direct application to the case of the heat kernel of a free scalar particle on a maximally symmetric space with Dirichlet boundary conditions at the origin, x(0) = x(β) = 0. In this section we present a possible way to easily extend the above computation to the case of non-zero boundary conditions. First of all, let us denote the endpoints with x(0) = y and x(β) = z for the particle x(τ ) propagating from time τ = 0 to time τ = β. The path integral we are interested in, is with where the potentialṼ (x) includes the counter-term (23) and a possible external potential. Now, we perform a splitting of the particle path into the classical one x cl (τ ) (i.e. a straight line, being the solution of the classical equation of motion for the free particle) and quantum fluctuations q(τ ), In particular, the fluctuations satisfy q(0) = q(β) = 0. Hence, the average of the path integral (41) takes the form with As equations (44)-(47) show, it is possible to calculate the averaged path integral I(y, z; β) with nonvanishing endpoints in terms of another one with null Dirichlet boundary conditions (like those computed in Section 4), provided a few replacements/adjustments, namely: • a factor which depends upon the squared difference of the endpoints appears in front of the average (44)-which corresponds to the free flat action; • the sampled points of the worldlines around zero (which here are denoted by q) must be displaced by the corresponding points of the classical path x cl when we evaluate the potentials (46) and (47).
In Figure 6 we report the calculation of the heat kernel of a free scalar particle propagating on a maximally symmetric space with D = 4, from point y = (0, 0, 0, 0) to point z = (5,5,5,5), done using both numerical methods introduced above, and showing very good agreement.

Conclusions and Outlook
In the past couple of decades the numerical Monte Carlo approach to the Worldline formalism in flat space has seen a variety of applications. In the meantime, worldline techniques have been extended to QFT in curved spaces, but restricted to perturbative, analytic computations. The present manuscript intends to be a first attempt to link these two avenues, in order to be able to extend the Worldline approach to QFT in curved space beyond analytic perturbation theory. Applications can be numerous. Firstly, the study of quantum field theory effective actions can be extended beyond the regime where the inclusion of gravity is treated perturbatively, as for instance in the derivation of gravitational corrections to the Euler-Heisenberg lagrangians [15]. On the other hand, the effect of curvatures in strongly-coupled fermionic models, such as the Gross-Neveu model, which describe the low-energy limit of several physical systems, can be studied with great efficiency with the present setup. In fact, in flat space, the Worldline Monte Carlo approach to large N fermionic models was already successfully considered [28,33], whereas, at the perturbative level, the heat kernel expansion of the Gross-Neveu model in 3d curved space with constant curvature was studied, for example in Ref. [34]. The results that we have presented here clearly show the possibility of performing worldline path integrals associated to scalar point particles moving in curved space by means of Monte Carlo numerical procedures.
In the present manuscript, we have considered maximally symmetric geometries. By using a recently studied effective flat space model-which encodes the curvature effects inside a suitable potential-as a benchmark, we have described a possible extension of the Worldline Monte Carlo to curved spaces with maximal symmetry. We have found it convenient to add a mass term in the kinetic action that serves to construct the worldline ensemble. Let us stress that, although the computation presented above focuses on maximally symmetric spaces, in the numerical NSMM approach there is a priori no particular constraint on the curved space background, nor is there any limitation on the chosen set of coordinates. However, as in other numerical methods, the main drawback of the present approach is the necessity of having a specific background, dependent on a finite set of real parameters.
The construction described in the present manuscript focuses on the computation of the diagonal part of a particle heat kernel which is linked to the computation of one-loop effective actions. However, a possible strategy to simulate open worldlines has been presented and might also be used to study dressed particle propagators, for instance. Finally, the inclusion of spinorial degrees of freedom would certainly be a welcome generalization. This problem can be tackled either with the inclusion of suitable matrix-valued potentials (spin factors) or in terms of spinning particle models, which involve Grassmann odd coordinates, along with the geometric, Grassmann even, coordinates. Numerical approaches in the presence of Grassmann odd variables were already considered in [35], and it would be helpful to apply such construction to the present worldline Monte Carlo method.

A Modified YLOOPS algorithm
To diagonalize the flat kinetic term (3), we first define the discretized derivative as the backward differenceẋ The index k parameterizes the discretized x-points in the same way as the continuous parameter τ parameterizes the continuous x-points. Now, the full kinetic term will be proportional to the summation which we generalize to Following the same steps of [29], and The algorithm finally reads as follows.
3. Construct the loop according to It is easy to see that coefficients C (α) k reproduce those of [29] for vanishing α, i.e. k+1 k .

B WLMC convergence with tiny mass term
The effect of adding a tiny mass term to the sampling algorithm of the worldlines (cfr. Eq. (50)) has an important effect on WLMC computations. In Figure 7, we report the case of two calculations of the heat kernel introduced in Section 4. Both calculations have been performed with the same parameters, namely N W L = 1000, N = 1000, D = 4 and M = 1: however, for the one on the left a tiny fictitious mass has been introduced in the sampling, whilst for the other it has not. It can be easily noticed that 7(a) follows with good accuracy the expected behaviour, while 7(b) not only has a lot more broadened data, but they also qualitatively deviate from the green reference points. In this context, the addition of a mass term has the role of a regulator for the sampling of the worldlines. We stress that this fictitious mass has not been taken into account during the evaluation of the potential; rather, its only role is to modify the diagonalization coefficients (53) appearing in (51) and (52). Finally, we investigated a possible qualitative dependence of the fictitious mass which minimizes the error between the non-linear sigma model method and the effective potential method, with respect to the curvature of the maximally symmetric space. We considered the hyperboloid case with |M | ranging from 1 to 0.01. However, nothing would forbid to consider the positive curvature case, and the results are expected to hold unchanged.
In Figure 8 we report the cases of |M | = 1 and |M | = 0.01 (intermediate curvatures exhibit the same behaviour). Each point of the plot represents the average error between a simulation with the nonlinear sigma model method and the associated effective potential method for different values of the α parameter. It can be seen that the optimization value α min does not change over two orders of magnitude of the curvature of the space. Hence we conclude that the fictitious mass parameter α is substantially curvature-independent.  Figure 8: Comparison between heat kernel calculations of a free scalar particle on the hyperboloid for |M | = 1, 0.01. Each point denotes a single simulation: on the horizontal axis we have the value of the fictitious mass parameter α adopted, whereas vertically we have reported an arithmetic average of distance between the data of the non-linear sigma model method and the ones of the effective potential method (assumed as a benchmark). The α-value which minimizes the error does not qualitatively depend on the curvature of the space.