Worldline Numerics for EnergyMomentum Tensors in Casimir Geometries
Abstract
We develop the worldline formalism for computations of composite operators such as the fluctuation induced energymomentum tensor. As an example, we use a fluctuating real scalar field subject to Dirichlet boundary conditions. The resulting worldline representation can be evaluated by worldline MonteCarlo methods in continuous spacetime. We benchmark this worldline numerical algorithm with the aid of analytically accessible singleplate and parallelplate Casimir configurations, providing a detailed analysis of statistical and systematic errors. The method generalizes straightforwardly to arbitrary Casimir geometries and general background potentials.
I Introduction
The worldline formalism Feynman (1950); Halpern and Siegel (1977); Bern and Kosower (1991); Strassler (1992) is a field theory technique with various computational advantages such as the reduction of the number of diagrams in perturbative expansions. It is particularly powerful for amplitude computations in external backgrounds Schmidt and Schubert (1993); Reuter et al. (1997); Shaisultanov (1996); Schubert (2001), as stringinspired techniques become highly efficient in this case. Moreover, fluctuationinduced quantities such as quantum actions, energies, and forces can also be efficiently computed by means of Monte Carlo simulations of the worldline path integral. Even fully nonperturbative problems can be tackled within the worldline formalism Affleck et al. (1982); Nieuwenhuis and Tjon (1996); Gies et al. (2005); Dietrich (2014); Bastianelli et al. (2014); Dietrich (2015).
The numerical method has given unprecedented access to background problems with nontrivial spacetime dependencies in QED Gies and Langfeld (2001, 2002); Langfeld et al. (2002); Schmidt and Stamatescu (2003); Mazur and Heyl (2015), QCD Deschner (2008); Epelbaum et al. (2015) fermionic systems Dunne et al. (2009), as well as the physics of Casimir forces Gies et al. (2003); Gies and Klingmuller (2006a, b, c); Weber and Gies (2009); Schaden (2009a, b); Aehlig et al. (2011); analogous analytical results for QED have also come within reach using worldline instanton approximations Dunne and Schubert (2005); Dunne et al. (2006); Dietrich and Dunne (2007); Strobel and Xue (2015); Linder et al. (2015); Ilderton et al. (2015); Dumlu (2015). Important advantages of the MonteCarlo method are that spacetime is not discretized but remains continuous in the algorithm, and renormalization can be performed on the level of finite and numerically controllable quantities. In addition the worldline formalism offers an intuitive picture of quantum fluctuations, as the worldlines can be viewed as spacetime trajectories of the fluctuating particles. The numerical method has also been applied to Minkowskivalued 2point functions for external legs Gies and Roessler (2011).
In the present work, we demonstrate that the technique can also be extended to composite operators, using the induced energymomentum tensor (EMT) as an example. The quest for such a method was first initiated in Schäfer et al. (2012). The new difficulty arises from the fact that standard tools such as pointsplitting lead to expressions involving the 2point function of the fluctuating field (internal line). As shown below, this can nevertheless be dealt with in a comparatively straightforward fashion by using open in addition to closed worldlines.
The present conceptual and in large parts technical work is motivated by the computation of energymomentum tensors for Casimir configurations. This problem is paradigmatic as the EMT sources the gravitational field. Since induced Casimir energies are typically negative (if associated with an attractive force) Casimir (1948); Bordag et al. (2001); Milton (2001, 2009); Bordag et al. (2009), information about the EMTs for Casimir configurations is a testbed for conceptual considerations of the interplay of quantum field theory and gravity, see, e.g., Fulling et al. (2007); Milton et al. (2007); Shajesh et al. (2008) for a discussion of how the Casimir energy “falls”.
Generically, many fundamental investigations of dynamical spacetime properties start with assumptions imposed on the properties of the EMT on the righthand side of Einstein’s equations. Most prominently, in order to exclude exotic phenomena Morris and Thorne (1988); Morris et al. (1988); Visser (1989); Kar et al. (2004); Fewster and Roman (2005); Olum (1998), energy conditions (ECs) on the properties of the EMT are imposed. For instance, certain “mechanisms” for superluminal travel can be excluded if an energy condition of the form is imposed, where is the tangent vector of a geodesic . If is a null vector this condition is called the null energy condition (NEC). The Casimir effect violates some energy conditions such as the NEC Graham and Olum (2003); Olum and Graham (2003); Graham and Olum (2005); Graham (2006). It can therefore not be used to rule out superluminal travel. In fact, superluminal (but still causality preserving) phase velocities are known to occur in the parallelplate Casimir configuration Scharnhorst (1990); Barton (1990); Barton and Scharnhorst (1993); Dittrich and Gies (1998, 2000). The Casimir configurations studied so far at least obey the weaker averaged NEC (ANEC). For a further discussion of this topic, see Graham and Olum (2005, 2007); Visser (1996a, b, c, 1997); Fewster et al. (2007); Graham and Olum (2007). In the present work, we use the violation of the NEC by the Casimir effect as an illustration for our new computational method. For conventional computational strategies to determine the induced EMT, for example, by means of mode summation or expansion, image charge methods, or similar techniques, see, e.g., Brown and Maclay (1969); Milton (2003, 2004); Scardicchio and Jaffe (2006).
This paper is intended to be a manual for the numerical computation of composite operators using the worldline formalism. It is organized as follows: in Sec. II we apply the worldline formalism to composite operators, specifically the EMT of a real scalar field obeying Dirichlet boundary conditions (BCs). We discuss the differences between the formalism for composite and noncomposite operators or functionals like the effective action. In Sec. III we test our numerical worldline algorithm by calculating components of the energymomentum tensor for a single plate with Dirichlet BCs analytically and numerically. The numerical calculation of the EMT and of the NEC is presented for the parallel plate configuration in Sec. IV. We compare our numerical results with known analytical results and discuss the arising systematic and statistical errors. We conclude with a summary of our results and an outlook on future worldline calculations.
Ii Worldline Formalism for the energymomentum tensor
Composite, or local, operators in quantum field theory are local products of field operators and their derivatives. Being distributionvalued, their product at the same spacetime point may be illdefined and give rise to divergences. Such operators can nevertheless be used for calculations after regularizing the divergences, for example, by pointsplitting, function, or dimensional regularization. In the following, we compute the vacuum expectation value of the energymomentum tensor operator of a scalar field; a composite operator constructed from the scalar field operator and its derivatives.
Our presentation of the basic formalism follows the one of Scardicchio and Jaffe (2006) with only a few differences. We study a quantum scalar field that is a map from a domain of dimensional Minkowski spacetime onto the linear space of selfadjoint operators on the Fock space
(1) 
has mass and is minimally coupled to a static classical background potential . Later on will impose boundary conditions on the fluctuations of . We use the dimensional Minkowski metric with the ”mostly minus” signature . From the Lagrangian density operator
(2) 
we derive the equation of motion for ,
(3) 
as well as its canonical energymomentum tensor operator ()
(4) 
The EMT operator in Eq. (4) contains products of field operators at the same spacetime point that can lead to divergences. We regularize them by a pointsplitting procedure in the spatial components , i.e.,
(5) 
The pointsplit EMT operator is then ()
(6) 
Our final goal is to compute the effects of boundary conditions on the energymomentum tensor. For the sake of convenience, we restrict ourselves to computing those EMT components that contain information relevant for the null energy condition (NEC) in the direction, where the coordinate is always the th coordinate of our spatial vectors . This null energy condition is given by the vacuum expectation value of the projection of the EMT onto a null curve with the null tangent vector ,
(7) 
where denotes the vacuum expectation value. All our calculations generalize straightforwardly to other EMT components. In order to compute the vacuum expectation values of Eq. (6), we expand in terms of spatial eigenmodes and corresponding creation and annihilation operators and ,
(8) 
where a discrete notation is used for both discrete and continuous parts of the spectrum for simplicity. The are defined as normalized eigenmodes of the Laplacian in the presence of the potential
(9) 
It is convenient to parametrize the energy eigenvalues also in terms of momenta according to . Then, the eigenvalue equation reads , which is a Helmholtztype equation. The corresponding Green’s function equation,
(10) 
is solved by the spectral representation
(11) 
We aim at expressing Eq. (7) in terms of . We therefore insert unity in the form
(12) 
and exchange the integration and the summation. This yields
(13) 
The sum inside the imaginary part is just . In general, instead of unity a decaying exponential must be inserted in order to construct local counterterms for renormalization (cf. Scardicchio and Jaffe (2006)). then serves as a cutoff for large momenta . This is, however, not necessary in our calculations because we are going to evaluate the EMT only at points for which , so that all local counterterms, that is, all counterterms that depend on , are automatically zero. Furthermore, since the term is in general complex, pulling it into the argument of the imaginary part generates an additional term that is proportional to . This term vanishes in the limit and hence will be ignored. The effects of boundary conditions imposed on the field fluctuations by are described by the difference between the EMT with nonvanishing potential and the EMT with . So far, we have left the potential arbitrary to emphasize that our calculations are independent of its specific properties. We can therefore repeat all the above steps with a vanishing potential. The only changes that occur are the mode functions in Eq. (9). The corresponding Green’s function is , defined by Eq. (10) for , rather than . We switch to the common point variables
(14) 
and with , evaluated at and we have
(15) 
Equation (15) is independent of any specific mode expansion of . It only depends on the Green’s function , which is representationindependent by definition. Therefore, any method for computing can be used at this stage, e.g., an optical approach in Scardicchio and Jaffe (2006). The subtraction of also removes the divergent terms that are independent of the potential and normalizes the EMT such that it vanishes for .
ii.1 The worldline representation of
For simplicity, we consider here only the massless case . Our calculations can straightforwardly be carried over to the massive case, for details, see App. A. In order to express the Green’s function in the worldline formalism, we interpret as the matrix element of an operator and Eq. (9) as a quantum mechanical Schrödinger problem whose Hamiltonian is . Then corresponds to a quantum mechanical propagator, Fourier transformed to energy space, from which the free motion has been subtracted. Hence, it can be written in position space in terms of a propertime representation:
(16) 
The matrix element in Eq. (16) is the quantum mechanical transition amplitude of a fictitious particle moving from at the fictitious time to at with a Hamiltonian . The corresponding Feynman path integral in position space is then straightforward to find. Next, we perform formal Wick rotations in both the and planes, such that and , which are consistent with causality. This casts into its doubly Wickrotated form
(17) 
The variables and are called Minkowskian and Euclidean propertime, respectively. Both describe the time evolution of the fictitious Schrödinger problem, but neither is a physical, measurable time. An analogous propertime representation for the effective action has been used in previous calculations of effective interaction energies for the Casimir effect and similar boundary configurations Gies et al. (2005); Gies and Klingmuller (2006c, b). The worldline representation of the effective action contains, however, a path integral over closed loops, whereas for , open worldlines running from to must be computed. The implicitly normalized path integral in Eq. (17) gives, for the free path integral, the standard free propagator,
(18) 
From this normalization one derives the shorthand notation for the path integral expectation value, the worldline average of an arbitrary operator :
(19) 
Pulling out the imaginary part prescription in front of the integral in Eq. (15), we can perform both Wick rotations for the components of the EMT. Details can be found in App. A. We finally arrive at the manifestly real worldline expressions for the normalized vacuum expectation values of and ,
(20) 
The expectation value of the path integral in Eq. (20) is defined by Eq. (19) with the operator . It is an average over worldlines that start at and end at , and that are weighted with a Gaußian velocity distribution. It is convenient to switch to coordinates such that for all paths are fixed with respect to one common point (cf. Fig. 1). We therefore call these worldlines common point loops or lines. The remaining path is written in terms of a dimensionless unit worldline with and .
The unit line itself can be written as the sum of the classical path from to and a deviation that obeys . We choose the origin of the coordinate system to be the point ,
(21) 
In the limit the worldline closes, it becomes a loop . We have set above
(22) 
The dependence in this relation is crucial for the computation of the derivatives in Eq. (20). Rescaling the worldlines, , makes the weight factor of the resulting path integral independent of . This allows us to compute the expectation value by generating only one ensemble of unit worldlines , which are all defined with respect to the same coordinate system. The worldline average is now calculated by replacing the path integral with a sum over a finite number of unit worldlines that are themselves approximated by a finite number of points per line with . The points are random numbers distributed according to the Gaußian weight factor , being a discretized realization of an open line with the vector connecting its endpoints. As a result, the expectation value of an operator is written as an average,
(23) 
The common point paths are conveniently generated with the loop algorithm Gies et al. (2005), which also works for open worldlines.
ii.2 Dirichlet constraints on
The entire formalism that we have outlined so far works, of course, for arbitrary static background field configurations . In our calculations, we use
(24) 
where is a surface element on . We recover Dirichlet BCs in the limit (cf. Graham et al. (2004, 2003)). The subtraction of the vacuum Green’s function in Eq. (15) already removed all divergences independent of . Therefore, can only diverge at points for which and the energymomentum tensor in Eq. (20) is finite on , where vanishes. All remaining divergences are located on the boundary and can be related to the infinite amount of energy necessary to constrain on all momentum scales to fulfill the Dirichlet boundary condition on .
The parametrization Eq. (24) and the subsequent Dirichlet limit greatly simplify the worldline average because now we have
(25) 
Equation (25) states that only paths which violate the boundary conditions lead to deviations from the trivial vacuum and thus contribute to the expectation value. We call the function the intersection condition. It gives a geometric description of how and for which values of the worldlines intersect the boundary. As long as the start and end point of a given worldline are not on the boundary (), the intersection condition will always determine a minimal nonzero value for the propertime for which the worldline first intersects . denotes the minimum propertime that is necessary for a worldline to propagate from the start to the end point and intersect a boundary in between. For any deviations from the straight line between and are typically strongly suppressed. Only for sufficiently large propertimes does the diffusive Brownian motion process, described by the path integral, create sufficiently large random detours that can intersect . In the propertime integral, serves as a lower bound and removes the divergence for . From a physics point of view, acts as an ultraviolet cutoff because small propertimes correspond to large momenta. The intersection condition may also provide an upper bound . This value would then mark the propertime for which the worldlines no longer intersect the boundary. This needs to be considered especially for configurations where the boundary consists of compact objects (see e.g. Gies and Klingmuller (2006c)).
ii.3 Compact expressions of and for worldline numerics
The EMT components are completely finite on by virtue of Eq. (25). We may therefore decompose and further and study the resulting, more compact, terms one by one. Toward this end, we parametrize the intersection condition as , where the function has inverse length dimension and describes the geometrical conditions for a worldline to intersect the boundary . It depends on and , and through the latter also on . This dependence is not explicitly known but it vanishes as . In general, is not a smooth function of these parameters. In fact, it can be nondifferentiable in either variable. Since depends on the shape of the boundary for the setup under consideration, its functional properties must be investigated each time anew. In addition, there can be several ways in which the intersection condition may be written. Whereas different parametrizations are equivalent, they can exhibit very different properties during numerical evaluation. Thus, the following calculations should be understood as primarily formal. First, we assume the function to be differentiable twice in both and . We additionally assume that all limits can be evaluated, and that only determines a lower bound on the propertime integral. In this case, we have .
We now insert Eq. (25) into Eqs. (20) and find
(26) 
In order to evaluate the propertime integral, we need to interchange several limits, which requires justification. First, we note that all expressions we deal with are finite by construction due to the Dirichlet constraint Eq. (25). Since this is valid for all , the limit may be interchanged with all other limits except . Furthermore, since we approximate the worldline average by a finite sum, this average can be exchanged with other limits and be conveniently computed last. We may also formally interchange propertime integration and differentiations because , and are independent variables. There is, however, the intersection condition , which depends on all three variables. The function also depends on the combination . As a consequence, the derivatives of , or more specifically of , must be computed before the propertime integration (cf. Sec. B). Picturing the worldlines as paths in space helps examine the situation: we need not only compute the intersection condition itself but also its derivatives, that is, we need to determine how the intersection is altered if or are changed. A derivative with respect to can be viewed geometrically as moving the complete worldline through space without changing its shape. By contrast, the derivative corresponds to opening and closing the worldline at a fixed point in space, changing its shape in the process. The required order in which these manipulations of the worldline expressions should be performed is then:

compute the derivatives of with respect to and ,

let , that is, ,

perform the propertime integration, and

average the expression over all worldlines in the ensemble.
With these consideration at hand, we write where we define
(27a)  
(27b) 
With we proceed accordingly and write it as a sum of four terms:
(28a)  
(28b)  
(28c)  
(28d) 
The fourth term, , comes from acting with on the exponential . Another term, the mixed term vanishes as .
We note that, with the help of the above decomposition into compact terms, the null energy condition along the axis Eq. (7) now reduces to computing
(29) 
Iii The energymomentum tensor for a single plate
The numerical computation of and in and space dimensions in the case that is a single dimensional surface, i.e., a plate, is our first proofofprinciple example. The plate imposes Dirichlet BCs on the fluctuations of . It is placed at such that its normal is the axis. The single Dirichlet plate configuration is also sometimes referred to as the perfect mirror.
iii.1 Analytic calculation for a single plate
Before we use worldline numerics, we compute the EMT for the single Dirichlet plate analytically. For that, we use Eq. (15) and compute by solving the equation of motion Eq. (9) for different boundary conditions. Denoting the BCs by a superscript , we must solve
(30) 
Toward that end, we decompose any in the dimensional vector parallel to and the component of . The solution of Eq. (30) in the half space with Dirichlet boundary conditions at is then
(31) 
where and . We assumed that there are only outgoing waves at spatial infinity (Sommerfeld radiation condition). In the same manner, the free solution without boundary conditions is found to be . Both solutions are normalized for and , . The Green’s functions and are now computed according to the spectral representation Eq. (11). We perform the momentum integration in polar coordinates and find for (see also Sommerfeld (1958); Jackson (1982); Polyanin and Zaitsev (1996))
(32) 
The are modified Bessel functions of the second kind (Macdonald functions).
We can now solve Eq. (15) analytically and even use the decomposition of and that we developed in Sec. II.3. Denoting the distance from the plate with , we find for the EMT in
(33a)  
(33b) 
And in the case of 3 spatial dimensions the values for and are
(34a)  
(34b) 
We note that and cannot be calculated separately in a direct manner since we used the functional structure of the worldline representation of to define these functions. Despite that, we can always compute from the sum using the fact that . We thus have
(35) 
According to Eq. (29) the NEC along the axis is then violated,
(36) 
These are the same values for the NEC which were derived in Graham and Olum (2005).
iii.2 Worldline calculation for one Dirichlet plate
The first step in all our worldline calculations is determining the intersection condition . For the single plate setup, this is easily done from Fig. 2.
The worldlines start at the point and intersect the plate at for all that fulfill . We call the component of the point on the loop that is closest to the plate . It is negative for our choice of coordinates in the setup Fig. 2. Thus we find for the function
(37) 
This corresponds to a minimal propertime