Pathwise Differentiation of Worldline Path Integrals

In path integrals such as Eq. (4) for the Casimir–Polder potential, one can compute the Casimir–Polder force $\mathbf{F}_{\mathrm{\scriptscriptstyle CP}}=-\nabla E_{\mathrm{\scriptscriptstyle CP}}$ by differentiating with respect to the source point $\mathbf{x}_{0}$ of the paths. A simple, direct way to compute such a gradient numerically is the finite-difference method. In this technique the derivative is computed in terms of the contribution of each path $\mathbf{x}(t)$ and its shifted counterparts of the form $\mathbf{x}(t)+\hat{r}_{i}\delta$, where $\delta$ is a small number, and $\hat{r}_{i}$ form an orthonormal basis for the configuration space. The simplest finite difference for the first derivative, for example, yields the following path integral for the force, from Eq. (4):

	$\displaystyle\mathbf{F}(\mathbf{x}_{\mathrm{0}})=$	$\displaystyle-\frac{\hbar c\alpha_{0}}{4(2\pi)^{D/2}\epsilon_{0}}\sum_{i}\hat{r}_{i}\int_{0}^{\infty}\frac{d\mathcal{T}}{\mathcal{T}^{1+D/2}}$
		$\displaystyle\times{\bigg{\langle}\kern-3.50006pt\bigg{\langle}}\frac{\langle\epsilon_{\mathrm{r}}[\mathbf{x}(t)+\hat{r}_{i}\delta]\rangle^{-3/2}-\langle\epsilon_{\mathrm{r}}[\mathbf{x}(t)]\rangle^{-3/2}}{\delta}{\bigg{\rangle}\kern-3.50006pt\bigg{\rangle}}_{\mathbf{x}(t)}\!\!.$		(5)

Unfortunately, as a numerical method this path integral performs poorly. Specifically, in the case of interfaces between regions of otherwise uniform dielectric permittivity, the contribution to the force path integral for a given path $\mathbf{x}(t)$ vanishes except when $\smash{\!\left\langle{\epsilon_{\mathrm{r}}}\right\rangle}$ differs for the shifted and original paths $\mathbf{x}(t)+\hat{r}_{i}\delta$ and $\mathbf{x}(t)$, respectively. Depending on how $\smash{\!\left\langle{\epsilon_{\mathrm{r}}}\right\rangle}$ is computed numerically, this difference may vary substantially among the possible paths $\mathbf{x}(t)$ in the limit of small $\delta$. For example, in the limit of an interface between vacuum and a perfect conductor, the contribution for a given path $\mathbf{x}(t)$ vanishes except in the (unlikely) case that $\mathbf{x}(t)+\hat{r}_{i}\delta$ crosses the interface but $\mathbf{x}(t)$ does not. For small $\delta$, only a small subset of the paths make a nonvanishing contribution to the path integral, and the sample variance diverges as $\delta\longrightarrow 0$. This problem becomes progressively worse when computing additional derivatives. Thus, it is important to consider alternative methods for computing derivatives.

Because of the presence of the $\mathcal{T}$ integral in the path integrals (1) and (4), it is possible in particular geometries to implement finite differences in a way that avoids the convergence problems noted above. This follows from noting that, for example, in the case of a planar dielectric interface, a shift in the source point of a path is equivalent to a change in the path’s running time $\mathcal{T}$, as far as computing  $\smash{\!\left\langle{\epsilon_{\mathrm{r}}}\right\rangle}$ is concerned. Thus, the finite difference can effectively act on the $\mathcal{T}^{-1-D/2}$ factor instead of the path average. However, the rest of this paper will focus on methods that apply in general geometries and to path averages that lack the $\mathcal{T}$ integral.

To handle derivatives with respect to the atomic position in the Casimir–Polder energy (4), consider the general, schematic form of the path integral

	$\displaystyle I$	$\displaystyle={\big{\langle}\kern-2.29996pt\big{\langle}}\Phi[\mathbf{x}(t)]{\big{\rangle}\kern-2.29996pt\big{\rangle}}_{\mathbf{x}(t)}$
		$\displaystyle=\int\prod_{n=1}^{N-1}d\mathbf{x}_{n}\,P[\mathbf{x}(t)]\,\Phi[\mathbf{x}(t)]$		(6)

in terms of an arbitrary path functional $\Phi[\mathbf{x}(t)]$ and probability measure $P[\mathbf{x}(t)]$ for the paths. The relevant paths for Casimir–Polder energies are closed, Gaussian paths (Brownian bridges) starting and ending at $\mathbf{x}(0)=\mathbf{x}(\mathcal{T})=\mathbf{x}_{0}$; in discrete form with $N$ discrete path samples, the probability density for these paths in $D-1$ dimensions is

	$\displaystyle P[\mathbf{x}(t)]=\,$	$\displaystyle(2\pi\Delta\mathcal{T})^{-N(D-1)/2}\,e^{-(\mathbf{x}_{0}-\mathbf{x}_{N-1})^{2}/2\Delta\mathcal{T}}$
		$\displaystyle\times e^{-(\mathbf{x}_{0}-\mathbf{x}_{1})^{2}/2\Delta\mathcal{T}}\prod_{j=1}^{N-2}e^{-(\mathbf{x}_{j+1}-\mathbf{x}_{j})^{2}/2\Delta\mathcal{T}},$		(7)

where $\Delta\mathcal{T}=\mathcal{T}/N$. To proceed, we will assume the path functional $\Phi$ to be approximately independent of $\mathbf{x}_{0}$ ($\partial\Phi/\partial\mathbf{x}_{0}\approx 0$). This is appropriate, for example, in Casimir–Polder calculations where the dielectric permittivity is constant in the vicinity of the atom (e.g., as in the case of an atom in vacuum, at a finite distance from a dielectric interface). In this case, derivatives of the path integral (6) with respect to $\mathbf{x}_{0}$ act only on the probability density $P$. For the Gaussian density (7), the derivatives along a particular spatial direction have the form

$$\frac{\partial^{n}}{\partial x_{0}^{n}}P[{x}(t)]=(\Delta\mathcal{T})^{-n/2}H_{n}\bigg{(}\frac{\bar{x}_{1}-x_{0}}{\sqrt{\Delta\mathcal{T}}}\bigg{)}P[{x}(t)],$$

(8)

where we are now indicating the spatial dependence in only the relevant direction for notational simplicity, $\bar{x}_{1}=(x_{1}+x_{N-1})/2$, and the Hermite polynomials $H_{n}$ are defined according to the convention

$$\frac{d^{n}}{dx^{n}}e^{-x^{2}}=(-1)^{n}H_{n}(x)\,e^{-x^{2}}.$$

(9)

Note that we have written out the derivatives only for one-dimensional paths to simplify the expressions, but the results here generalize straightforwardly to multiple spatial dimensions (including multiple derivatives in one direction as well as cross-derivatives), in addition to non-Gaussian path measures (e.g., for paths corresponding to more general Lévy processes).

Combining Eq. (8) with the path integral (6) gives the expression

$\displaystyle\frac{\partial^{n}}{\partial x_{0}^{n}}I$

$\displaystyle=(\Delta\mathcal{T})^{-n/2}{\bigg{\langle}\kern-3.50006pt\bigg{\langle}}H_{n}\bigg{(}\frac{\bar{x}_{1}-x_{0}}{\sqrt{\Delta\mathcal{T}}}\bigg{)}\Phi[{x}(t)]{\bigg{\rangle}\kern-3.50006pt\bigg{\rangle}}_{{x}(t)}$

(10)

for the path integral derivatives. However, this expression is problematic as a numerical derivative estimate because of the factor $\Delta\mathcal{T}^{-n/2}$, which diverges in the continuum limit $\Delta\mathcal{T}\longrightarrow 0$. This factor indicates large sample variances and cancellations between paths in this limit, resulting in inefficient numerical computations.

One approach to ameliorating this problematic behavior comes from the observation that if the path functional $\Phi$ is approximately independent of $\mathbf{x}_{0}$, then in the limit of large $N$, $\Phi$ will also be approximately independent of neighboring points, such as $\mathbf{x}_{1}$, $\mathbf{x}_{N-1}$, and other nearby points. In this case, it is possible to carry out the integrals with respect to these points, to obtain a “partially averaged” path integral. (Note that this method is distinct from other partial-averaging methods [28, 29].) Carrying out the integrals over $x_{1},\ldots,x_{m-1}$ and $x_{N-k+1},\ldots,x_{N-1}$, Eq. (10) becomes

$\displaystyle\frac{\partial^{n}}{\partial x_{0}^{n}}I$

$\displaystyle\approx(2\nu_{mk})^{-n/2}{\bigg{\langle}\kern-3.50006pt\bigg{\langle}}H_{n}\bigg{(}\frac{\bar{x}_{mk}-x_{0}}{\sqrt{2\nu_{mk}}}\bigg{)}\Phi[\mathbf{x}(t)]{\bigg{\rangle}\kern-3.50006pt\bigg{\rangle}}_{{x}(t)},$

(11)

where $\nu_{mk}$ and $\bar{x}_{mk}$ are defined by

	$\displaystyle\nu_{mk}$	$\displaystyle:=\frac{mk\Delta\mathcal{T}}{m+k}$		(12)
	$\displaystyle\bar{x}_{mk}$	$\displaystyle:=\frac{kx_{m}+mx_{N-k}}{m+k},$		(13)

for $k,m\geq 1$. The probability density for the paths in the partially averaged expression (11) is

	$\displaystyle\bar{P}[{x}(t)]$	$\displaystyle:=\int dx_{1}\ldots dx_{m-1}dx_{N-k+1}\ldots dx_{N-1}\,P[x(t)]$
		$\displaystyle=e^{-({x}_{0}-{x}_{N-k})^{2}/2k\Delta\mathcal{T}}e^{-({x}_{0}-{x}_{m})^{2}/2m\Delta\mathcal{T}}$
		$\displaystyle\hskip 14.22636pt\times\prod_{j=m}^{N-k-1}e^{-({x}_{j+1}-{x}_{j})^{2}/2\Delta\mathcal{T}},$		(14)

which has “large steps” from $x_{0}$ to $x_{m}$ and $x_{N-k}$ to $x_{0}$. For closed paths, it is generally most sensible to choose $m=k$ to apply the same partial averaging to the beginning and end of the path, such that $\nu_{m}:=\nu_{mm}$ and $\bar{x}_{m}:=\bar{x}_{mm}$ simplify to

	$\displaystyle\nu_{m}$	$\displaystyle=\frac{m\Delta\mathcal{T}}{2}=\frac{m\mathcal{T}}{2N}$		(15)
	$\displaystyle\bar{x}_{m}$	$\displaystyle=\frac{x_{m}+x_{N-m}}{2},$		(16)

and the steps at the beginning and end of the path are equally large on average, as shown schematically (omitting the Brownian nature of the path for visual clarity) in Fig. 1.

Since the prefactor in the derivative (11) now scales as $\smash{\nu_{m}^{-n/2}}$, there is clearly some reduction in the sample variance in a numeric computation as $m$ increases. However, the critical point is this: when chosen appropriately, this prefactor becomes independent of $N$, curing the divergent sample variance in the continuum limit. To choose the number $m$ of points to partially average, the main consideration is that the averaging should not introduce significant error. Returning to the assumption that $\Phi$ is approximately independent of $\mathbf{x}_{0}$, at this point we will have to be somewhat more precise and assume that $\Phi$ depends on the path $\mathbf{x}(t)$ via a scalar function $V(\mathbf{x})$ in the form $\Phi[V(\mathbf{x})]$. Then the key assumption is that $V(\mathbf{x})$ is constant (or approximately constant) within a radius $d$ of the source point $\mathbf{x}_{0}$. The error in averaging over the point $\mathbf{x}_{m}$, while ignoring the corresponding dependence of $\Phi$, is small so long as there is a small probability for $\mathbf{x}_{m}$ to leave the constant region. Assuming $m$ will be chosen to be small compared to $N$, the path from $x_{0}$ to $x_{m}$ (and the corresponding path segment at the end) will behave approximately as an ordinary Wiener path. Since the probability for a one-dimensional Wiener path to cover a distance $d>0$ from $x_{0}$ in time $\tau$ is given by

$$P_{\text{cross}}=\mathrm{erfc}\left(\frac{d}{\sqrt{2\tau}}\right).$$

(17)

For $\tau\ll d^{2}$, this expression is bounded above by $e^{-d^{2}/2\tau}$. Then to keep the error acceptably small, the crossing probability is chosen to be below some small tolerance $\varepsilon>0$. This condition can be solved to find the total average step time $\tau=m\Delta\mathcal{T}$ such that $P_{\text{touch}}\leq\varepsilon$. The fraction of the path one can acceptably average over is then given by

$$\frac{m}{N}=\frac{d^{2}}{2\mathcal{T}\ln(1/\varepsilon)}.$$

(18)

(In practice we choose $\varepsilon$ to be much smaller than suggested by this formula.) Note that $m/N$ is constant even as the path resolution $N$ increases, so that the numerical fluctuations, governed by $\nu_{m}=m\mathcal{T}/2N$, are independent of $N$, as required. A bound for the numerical error then follows from multiplying the tolerance $\varepsilon$ by the amount over which $V(x)$ varies from $\mathbf{x}_{0}$ to the exterior of the constant region.

These apply readily to the computation of derivatives of the Casimir–Polder potential. For example, the path integral for the Casimir–Polder force becomes

	$\displaystyle\mathbf{F}=$	$\displaystyle\,\frac{\hbar c\alpha_{0}}{4(2\pi)^{D/2}}\int\frac{d\mathcal{T}}{\mathcal{T}^{1+D/2}}{\bigg{\langle}\kern-3.50006pt\bigg{\langle}}\frac{2\mathbf{x}_{0}-\mathbf{x}_{m}-\mathbf{x}_{N-m}}{m\Delta\mathcal{T}}$
		$\displaystyle\hskip 71.13188pt\times\Big{(}\langle\epsilon_{\mathrm{r}}\rangle^{-3/2}-[\epsilon_{\mathrm{r}}(\mathbf{x}_{0})]^{-3/2}\Big{)}{\bigg{\rangle}\kern-3.50006pt\bigg{\rangle}}_{\mathbf{x}(t)}$		(19)

after partial averaging, where an appropriate estimator for the path average of the permittivity is

	$\displaystyle\langle\epsilon_{\mathrm{r}}\rangle=$	$\displaystyle\,\frac{1}{N}\Bigg{[}m\epsilon_{\mathrm{r}}(\mathbf{x}_{0})+\frac{m}{2}\epsilon_{\mathrm{r}}(\mathbf{x}_{m})+\frac{m}{2}\epsilon_{\mathrm{r}}(\mathbf{x}_{N-m})$
		$\displaystyle\hskip 71.13188pt+\sum_{j=m+1}^{N-m-1}\epsilon_{\mathrm{r}}(\mathbf{x}_{j})\Bigg{]}.$		(20)

Note that the details of how the points $\mathbf{x}_{m}$ and $\mathbf{x}_{N-m}$ enter this estimator are not critical, since for small tolerances $\varepsilon$ the dielectric permittivity at these points should be equal to the source-point permittivity for the vast majority of paths. Fig. 1 shows an example of the physical setup for this path integral for an atom near a planar dielectric surface. In evaluating this path integral numerically here, the distance scale $d$ is simply the distance between the atom and the dielectric interface, and a small tolerance $\varepsilon$ guarantees that the points $x_{m}$ and $x_{N-m}$ are unlikely to cross over the interface.

In this differentiated path integral for the force, in principle no renormalization term is required—any renormalization term is independent of the position of the atom, and thus vanishes under the derivative. That is, the $\epsilon_{\mathrm{r}}(\mathbf{x}_{0})$ term in Eq. (19) does not ultimately contribute, because it vanishes under the ensemble average for any value of $\mathcal{T}$. However, on a path-wise basis this renormalization term is critical in a numerical implementation of this path integral, because it cuts off the divergence at $\mathcal{T}=0$ and dramatically reduces the sample variance of finite ensemble averages, particularly for small $\mathcal{T}$.

The partial-averaging approach is also useful in computing derivatives in a geometry with two or more material bodies, when the source point of the path integral is close to one of the bodies. This situation applies, for example, when computing the two-body stress-energy tensor in the vicinity of one of the bodies [25, 16, 18], or when extracting the nonadditive, multibody contribution to the Casimir–Polder potential.

As a concrete example, consider differentiating the Casimir–Polder potential of an atom between two dielectric half-spaces, where the atom is close to one interface (Fig. 2). In this problem there are two relevant scales for the path running time $\mathcal{T}$: the first time $\mathcal{T}_{1}$ that the path touches either body, and the first time $\mathcal{T}_{2}$ that the path touches both bodies. For the atom at position $x_{0}$ between interfaces at positions $d_{1}$ and $d_{2}$, these times are on the order of $\mathcal{T}_{1}\sim(d_{1}-x_{0})^{2}$ and $\mathcal{T}_{2}\sim(d_{2}-d_{1})^{2}$, respectively. If the atom is particularly close to one body, the dominant contribution to the worldline path integral comes from paths where $\mathcal{T}\sim\mathcal{T}_{2}\gg\mathcal{T}_{1}$. The approach of the previous section is limited in this case, because it relies on the dielectric permittivity being approximately constant in a neighborhood of the source point, with the effectiveness of the method increasing with the radius of this neighborhood. The method thus only effects partial averaging to a limited extent, for $\mathcal{T}_{m}:=m\Delta\mathcal{T}\ll\mathcal{T}_{1}$.

The partial-averaging method extends to this case if there are approximate analytical solutions available for the path integral. As in the case of the TE-polarization path integral [19], it is possible to find analytical solutions for path integrals corresponding to certain material geometries, and then to recast more general worldline path integrals in a form that leverages those solutions. In particular, the (un-renormalized) Casimir–Polder energy can be written as an ensemble average over paths divided into $N$ segments. The ensemble average in the path integral must then include an average over the positions of $N$ points that define a discrete path, as well as averages over all possible path segments connecting each pair of neighboring points:

	$\displaystyle V_{\mathrm{\scriptscriptstyle CP}}^{\mathrm{\scriptscriptstyle(TE)}}(\mathbf{x}_{\mathrm{0}})=\frac{\hbar c\alpha_{0}}{(2\pi)^{D/2}\sqrt{\pi}\epsilon_{0}}\int_{0}^{\infty}\!\!\,\frac{d\mathcal{T}}{\mathcal{T}^{1+D/2}}\,\int_{0}^{\infty}\!\!ds\,s^{2}\,e^{-s^{2}}\,$
	$\displaystyle\hskip 2.84526pt\times{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\prod_{j=0}^{N-1}{\kern 0.0pt\Bigg{\langle}\kern-6.60004pt\scalebox{0.66}[1.0]{$-$}\kern-2.20001pt\Bigg{\langle}\kern 0.0pt}\exp\!\Bigg{[}-\!\frac{s^{2}}{\mathcal{T}}\,\int_{\mathcal{T}_{j}}^{\mathcal{T}_{j+1}}\!\!d\tau\,\chi[\mathbf{x}(\tau)]\Bigg{]}{\kern 0.0pt\Bigg{\rangle}\kern-2.20001pt\scalebox{0.66}[1.0]{$-$}\kern-6.60004pt\Bigg{\rangle}\kern 0.0pt}_{\!\!\Delta x_{j}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}(\tau)}\!\!\!\!.$		(21)

Here, the connected-double-angle brackets $\smash{{\langle\kern-2.79999pt\scalebox{0.35}[1.0]{$-$}\kern-1.79993pt\langle}\;{\rangle\kern-1.79993pt\scalebox{0.35}[1.0]{$-$}\kern-2.79999pt\rangle}_{\Delta x_{j}}}$ denote an ensemble average over all bridges between $\smash{\mathbf{x}_{j}}$ and $\smash{\mathbf{x}_{j+1}}$, $\mathcal{T}_{j}=j\Delta\mathcal{T}=(j/N)\mathcal{T}$, and the dielectric susceptibility for the two bodies is $\chi[\mathbf{x}(\tau)]=\chi_{1}[\mathbf{x}(\tau)]+\chi_{2}[\mathbf{x}(\tau)]$. With this ordering of brackets, the averages over all path segments between points $\smash{\mathbf{x}_{j}}$ and $\smash{\mathbf{x}_{j+1}}$ occur for a fixed set of discrete points, and then the average must be taken over all possible discrete paths $\mathbf{x}_{j}$. The key point of this section is that for simple geometries, such as a single planar interface, it is possible to carry out the average over the path segments analytically. In the example of a single, planar interface, the solution corresponds to the generating function of the sojourn time of a Brownian bridge (see Ref [19], App. B).

For sufficiently small time steps ($\Delta\mathcal{T}\ll\mathcal{T}_{2}$), the probability for the path to interact with both bodies on a given discrete step $j$ is small, so the path integral over bridges is well-approximated by accounting only for the nearest body. For path increments $\Delta\mathbf{x}_{j}$ close to the $\chi_{2}$ body, the path-segment average is then

$${\kern 0.0pt\Big{\langle}\kern-4.90005pt\scalebox{0.38}[1.0]{{\hbox{-}}}\kern-1.79993pt\Big{\langle}\kern 0.0pt}e^{-(s^{2}/\mathcal{T})\int d\tau\,(\chi_{1}+\chi_{2})}{\kern 0.0pt\Big{\rangle}\kern-1.79993pt\scalebox{0.38}[1.0]{{\hbox{-}}}\kern-4.90005pt\Big{\rangle}\kern 0.0pt}_{\!\!\Delta\mathbf{x}_{j}}\approx{\kern 0.0pt\Big{\langle}\kern-4.90005pt\scalebox{0.38}[1.0]{{\hbox{-}}}\kern-1.79993pt\Big{\langle}\kern 0.0pt}e^{-(s^{2}/\mathcal{T})\int d\tau\,\chi_{2}}{\kern 0.0pt\Big{\rangle}\kern-1.79993pt\scalebox{0.38}[1.0]{{\hbox{-}}}\kern-4.90005pt\Big{\rangle}\kern 0.0pt}_{\!\!\Delta\mathbf{x}_{j}}.$$

(22)

The crucial point is that this path integral (including the Gaussian probability measure for the path increment), which we can write as

$\displaystyle U(\mathbf{x}_{j},\mathbf{x}_{j+1};t)$

$\displaystyle=\frac{e^{-(\mathbf{x}_{j}-\mathbf{x}_{j+1})^{2}/2t}}{(2\pi t)^{(D-1)/2}}{\kern 0.0pt\Big{\langle}\kern-4.90005pt\scalebox{0.38}[1.0]{$-$}\kern-1.79993pt\Big{\langle}\kern 0.0pt}e^{-s^{2}\!\int_{0}^{t}\!d\tau\,\chi[\mathbf{x}(\tau)]/\tau}{\kern 0.0pt\Big{\rangle}\kern-1.79993pt\scalebox{0.38}[1.0]{$-$}\kern-4.90005pt\Big{\rangle}\kern 0.0pt}_{\Delta\mathbf{x}_{j}},$

(23)

corresponds to a propagator for an associated diffusion equation [30, 31, 32, 33], so the solutions between adjacent steps can be composed into the solution for a larger, combined step. The partial averaging can then be carried out as in the method of the previous section by composing $m$ steps together. To ensure that the error associated with the partial averaging is small, the probability to touch the more-distant body in time $\mathcal{T}_{m}=m\mathcal{T}/N$ should remain below some small tolerance $\varepsilon$. In the example of two planar interfaces, this condition has the explicit form

$$\frac{m}{N}=\frac{(d_{2}-x_{0})^{2}}{2\mathcal{T}\ln(1/\varepsilon)},$$

(24)

using the same argument that led to Eq. (18). In more general material geometries, the analytical path integral must also approximate the shape of the nearest body. For example, a smooth but nonplanar interface (e.g., the surface of a spherical body) can be approximated locally by a plane for the purposes of the path-segment average in a given path increment. However, in this case, $\mathcal{T}_{m}$ is also constrained because the geometric approximation must hold over the typical extent of the path up to time $\mathcal{T}_{m}$. For example, in a local planar approximation, the distance that paths typically explore in time $\mathcal{T}_{m}$ must be small on the curvature scale of the interface.

After averaging, the derivatives are readily computed, although in general they will not yield the same Hermite-polynomial weighting factors as in Eq. (10):

	$\displaystyle\frac{\partial^{n}}{\partial x_{0}^{n}}V_{\mathrm{\scriptscriptstyle CP}}^{\mathrm{\scriptscriptstyle(TE)}}(\mathbf{x}_{\mathrm{0}})\!=\!\frac{\hbar c\alpha_{0}}{(2\pi)^{D/2}\sqrt{\pi}\epsilon_{0}}\int_{0}^{\infty}\!\!\frac{d\mathcal{T}}{\mathcal{T}^{1+D/2}}\int_{0}^{\infty}\!\!ds\,s^{2}\,e^{-s^{2}}\,$
	$\displaystyle\hskip 14.22636pt\times\int\prod_{k=m}^{N-m}d\mathbf{x}_{k}\prod_{j=m}^{N-m-1}U(\mathbf{x}_{j},\mathbf{x}_{j+1};\Delta\mathcal{T})$
	$\displaystyle\hskip 14.22636pt\times\frac{\partial^{n}}{\partial x_{0}^{n}}[U(\mathbf{x}_{0},\mathbf{x}_{m};m\Delta\mathcal{T})\,U(\mathbf{x}_{N-m},\mathbf{x}_{0};m\Delta\mathcal{T})].$		(25)

As discussed at the end of the previous section, this expression should be renormalized as necessary. Note that in the path integral here, the derivatives act on both the Gaussian path measure and the subaverages over path segments [via the bracketed factor in Eq. (23)]. In addition, the functional form of $U$ may vary depending on the local path position $\mathbf{x}_{j}$, because it should reflect only the presence of the interface nearest $\mathbf{x}_{j}$. Although the approximate solution near the source point $\mathbf{x}_{0}$ of the path is crucial for the partial averaging, it is not necessary to use such analytical results for the entire path. For example, one could also use the simpler trapezoidal estimator [19] for the average of the dielectric function around the remainder of the path.

The averaging approach described here helps to improve the sample variance in a numerical average when the source point of the paths is close to an interface, even though the functional integral cannot be treated as approximately constant, as assumed in the previous method of Section II.2. Both this and the previous methods of path-averaging described can be thought of as choosing non-uniform time steps in the discrete path. It is important to note, however, that in most path integrals, all of the time steps should be vanishingly small for good accuracy; but in the partially averaged path integrals here, the “large” first and last steps should be of fixed time step, even in the limit as the remaining step sizes become vanishingly small. Since the derivatives act on the first and last path steps, the sample variance decreases as their respective time steps increase. In order to improve the statistical convergence of the path average, the first and last steps should therefore be as large as possible, while controlling the errors associated with the simplified, approximate integrands. The methods here still allow a high resolution for the middle section of the path, in order to accurately capture the geometry dependence of the path integral.

The force on a body follows from the gradient of the Casimir energy, where the derivatives are taken with respect to the body’s position. For example, the force on body $1$ is given by differentiating the path integral in Eq. (26) with respect to the body position $\mathbf{R}_{1}$:

$$\mathbf{F}_{1}=-\nabla_{\mathbf{R}_{1}}E=-\,{\bigg{\langle}\kern-3.50006pt\bigg{\langle}}\frac{\partial\Phi}{\partial f}\,\nabla_{\mathbf{R}_{1}}f{\bigg{\rangle}\kern-3.50006pt\bigg{\rangle}}_{\mathbf{x}(t)}.$$

(28)

In this expression, $\partial\Phi/\partial f$ is independent of $f$ (since $\Phi$ is a linear function of $f$). In the case where $f$ is exponential, we can divide the path average into two steps: first, averaging over the $N$ discrete, uniform steps $\Delta\mathbf{x}_{j}:=\mathbf{x}_{j+1}-\mathbf{x}_{j}$, and second, averaging over all the subpaths between $\mathbf{x}_{j}$ and $\mathbf{x}_{j+1}$:

$$\mathbf{F}_{1}=-\,{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,\nabla_{\mathbf{R}_{1}}\prod_{j=0}^{N-1}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{j}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}}.$$

(29)

Here, the connected-double-angle brackets $\smash{{\langle\kern-2.79999pt\scalebox{0.35}[1.0]{$-$}\kern-1.79993pt\langle}~{}{\rangle\kern-1.79993pt\scalebox{0.35}[1.0]{$-$}\kern-2.79999pt\rangle}_{\Delta\mathbf{x}_{j}}}$ denotes an average over all such bridges between $\mathbf{x}_{j}$ and $\mathbf{x}_{j+1}$, as in Sec. II.3, and the subscript $\mathbf{x}_{n}$ indicates an average over all discrete paths of $N$ steps. Then applying the derivative to the product of subpaths,

$$\mathbf{F}_{1}=-\,{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,f\sum_{j=0}^{N-1}\frac{\nabla_{\mathbf{R}_{1}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{j}}}{{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{j}}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}}.$$

(30)

There are typically $N$ terms here, but if $\smash{{\langle\kern-2.79999pt\scalebox{0.35}[1.0]{$-$}\kern-1.79993pt\langle}f{\rangle\kern-1.79993pt\scalebox{0.35}[1.0]{$-$}\kern-2.79999pt\rangle}_{\Delta\mathbf{x}_{j}}}$ is associated with a sharp surface, only $\smash{O(N^{1/2})}$ of them will contribute to the sum. Again, this form of the path integral is specific to an exponential $f$; for example, Eq. (2) is of the form

	$\displaystyle f\big{(}\langle\chi\rangle\big{)}=e^{-s^{2}\mathcal{T}(1+\langle\chi\rangle)/2c^{2}}$		(31)
	$\displaystyle\Phi=-\frac{\hbar}{(2\pi)^{5/2}}\int_{0}^{\infty}\!\!ds\int_{0}^{\infty}\!\frac{d\mathcal{T}}{\mathcal{T}^{5/2}}\int d\mathbf{x}_{0}\,f(\langle\chi\rangle),$

and works with this method. Other path integrals, for example of the form of Eq. (1),

	$\displaystyle f(\langle\chi\rangle)=(1+\langle\chi\rangle)^{-1/2}$		(32)
	$\displaystyle\Phi=-\frac{\hbar c}{8\pi^{2}}\int_{0}^{\infty}\!\frac{d\mathcal{T}}{\mathcal{T}^{3}}\int d\mathbf{x}_{0}\,f(\langle\chi\rangle),$

can be handled using the alternate form for the force,

$$\mathbf{F}_{1}=-\,{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,f^{\prime}\sum_{j=0}^{N-1}\nabla_{\mathbf{R}_{1}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-3.80005pt\big{\langle}\kern 2.29996pt}\chi{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{j}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}},$$

(33)

where $\mathcal{T}_{j}=j\Delta\mathcal{T}=j\mathcal{T}/N$. Qualitatively, the Casimir force on a body arises from paths that pierce its surface, with a vector path weight in the direction of the surface normal at the path source point. Although the surface of an arbitrary body involves surface normals pointing in all directions, each surface normal obtains a different geometry-dependent weight via the path ensemble. The result is, in general, a nonzero net force.

As in Sec. II.2, a partial-averaging procedure accelerates convergence. If one surface is located at $d_{1}$, then the condition is

$$(x_{j+1}-d_{1})(x_{j}-d_{1})\leq\varepsilon$$

(34)

for a positive tolerance $\varepsilon$ [typically $10^{4}/(N\sqrt{\mathcal{T}})$; note that setting $\varepsilon=0$ simply detects a surface crossing]. When this condition is satisfied for any boundary, the solution advances forward by $m>1$ points, and the estimate of the function ${\langle\kern-2.79999pt\scalebox{0.35}[1.0]{$-$}\kern-1.79993pt\langle}f{\rangle\kern-1.79993pt\scalebox{0.35}[1.0]{$-$}\kern-2.79999pt\rangle}_{\Delta\mathbf{x}_{j}}$ or the susceptibility ${\langle\kern-2.79999pt\scalebox{0.35}[1.0]{$-$}\kern-1.79993pt\langle}\chi{\rangle\kern-1.79993pt\scalebox{0.35}[1.0]{$-$}\kern-2.79999pt\rangle}_{\Delta\mathbf{x}_{j}}$ is calculated with respect to this larger path length.

This method can be easily extended to the second derivative of the worldline energy, which yields the potential curvature,

$$C_{ij}:=(\hat{r}_{i}\cdot\nabla_{\mathbf{R}_{1}})(\hat{r}_{j}\cdot\nabla_{\mathbf{R}_{1}})E,$$

(35)

in basis $\hat{r}_{i}$. Applying this to the path integral (26) gives

$$C_{ij}={\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,(\hat{r}_{i}\cdot\nabla_{\mathbf{R}_{1}})(\hat{r}_{j}\cdot\nabla_{\mathbf{R}_{1}})\prod_{k=0}^{N-1}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}},$$

(36)

and applying the derivatives in the case where $f$ is exponential yields

	$\displaystyle C_{ij}=$	$\displaystyle{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,f\Bigg{[}\sum_{\begin{subarray}{c}\ell,k=0\\ \ell\neq k\end{subarray}}^{N-1}\frac{\hat{r}_{i}\cdot\nabla_{\mathbf{R}_{1}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{\ell}}}{{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{\ell}}}\frac{\hat{r}_{j}\cdot\nabla_{\mathbf{R}_{1}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}{{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}$		(37)
		$\displaystyle\hskip 28.45274pt+\sum_{k=0}^{N-1}\frac{(\hat{r}_{i}\cdot\nabla_{\mathbf{R}_{1}})(\hat{r}_{j}\cdot\nabla_{\mathbf{R}_{1}}){\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}{{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}\Bigg{]}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}}.$

However, it is possible to recast one derivative with respect to the other surface,

$$C_{ij}=-\,{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,(\hat{r}_{i}\cdot\nabla_{\mathbf{R}_{2}})(\hat{r}_{j}\cdot\nabla_{\mathbf{R}_{1}})\prod_{k=0}^{N-1}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{{\hbox{-}}}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}},$$

(38)

in which case, applying the derivatives, in the case where $f$ is exponential, and making the approximation $\Delta\mathbf{x}_{j}\ll d$ for well separated bodies, where $d$ characterizes the body separation, we obtain

	$\displaystyle C_{ij}\approx$	$\displaystyle-\,{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,f\Bigg{[}\sum_{k}\frac{\hat{r}_{i}\cdot\nabla_{\mathbf{R}_{2}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}{{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}\Bigg{]}$		(39)
		$\displaystyle\hskip 42.67912pt\times\Bigg{[}\sum_{k}\frac{\hat{r}_{j}\cdot\nabla_{\mathbf{R}_{1}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}{{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}f{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}}\Bigg{]}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}}.$

In the general case, where $f$ is any function, in the same approximation we obtain

	$\displaystyle C_{ij}\approx$	$\displaystyle-\,{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,f^{\prime\prime}\sum_{k}\hat{r}_{i}\cdot\nabla_{\mathbf{R}_{2}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}\chi{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{k}}$		(40)
		$\displaystyle\hskip 42.67912pt\times\sum_{\ell}\hat{r}_{j}\cdot\nabla_{\mathbf{R}_{1}}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}\chi{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{\ell}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}_{n}}.$

Partial averaging is then carried out as in Sec. III.1.

The Casimir torque on a body follows in a similar way by considering the energy shift due to a small rotation. Specifically, consider an infinitesimal rotation of body 1 about the point $\mathbf{R}_{1}$, with axis $\hat{m}$, and through angle $\phi$:

$$\chi(\mathbf{r})=\chi_{1}\Theta\big{\{}\sigma_{1}[\mathcal{R}(\phi)(\mathbf{r}-\mathbf{R}_{1})]\big{\}}+\chi_{2}\Theta[\sigma_{2}(\mathbf{r}-\mathbf{R}_{2})].$$

(41)

For a small rotation angle, the rotation matrix has the form

$$\mathcal{R}_{ij}(\phi)=\delta_{ij}-m_{k}\epsilon_{ijk}\phi+O(\phi^{2}),$$

(42)

where $\delta_{ij}$ is the Kronecker delta and $\epsilon_{ijk}$ is the Levi–Civita symbol; repeated indices are summed throughout this section. The torque for a rotation about axis $\hat{m}$ can be written as $K_{m}:=\hat{m}\cdot\mathbf{K}=-\partial_{\phi}E$. The $\phi$-derivative acts only on the path-averaged dielectric part of the energy integral,

	$\displaystyle\partial_{\phi}\langle\chi\rangle$	$\displaystyle=\chi_{1}\big{\langle}\partial_{\phi}\mathcal{R}_{ij}(\phi)[\mathbf{x}(t)-\mathbf{R}_{1}]_{j}[\hat{r}_{i}\cdot\nabla\Theta(\sigma_{1})]\big{\rangle}$
		$\displaystyle=-\chi_{1}\big{\langle}m_{k}\epsilon_{kij}[\mathbf{x}(t)-\mathbf{R}_{1}]_{j}[\hat{r}_{i}\cdot\nabla\Theta(\sigma_{1})]\big{\rangle}$
		$\displaystyle=\chi_{1}\hat{m}\cdot\big{\langle}[\mathbf{x}(t)-\mathbf{R}_{1}]\wedge\nabla\Theta(\sigma_{1})\big{\rangle},$		(43)

where in the last step we introduced the usual vector cross-product $(\mathbf{a}\wedge\mathbf{b})_{i}=\epsilon_{ijk}a_{j}b_{k}$. Substituting this into Eq. (26), the resulting path integral for the torque, after dropping the arbitrary axis $\hat{m}$, is

$\displaystyle\mathbf{K}$

$\displaystyle=\,{\Bigg{\langle}\kern-3.50006pt\Bigg{\langle}}\frac{\partial\Phi}{\partial f}\,f^{\prime}\chi_{1}\sum_{j=0}^{N-1}{\kern 2.29996pt\big{\langle}\kern-5.89996pt\scalebox{0.35}[1.0]{$-$}\kern-3.80005pt\big{\langle}\kern 2.29996pt}\nabla\Theta(\sigma_{1})\wedge[\mathbf{x}(t)-\mathbf{R}_{1}]{\kern 2.29996pt\big{\rangle}\kern-3.80005pt\scalebox{0.35}[1.0]{$-$}\kern-5.89996pt\big{\rangle}\kern 2.29996pt}_{\Delta\mathbf{x}_{j}}{\Bigg{\rangle}\kern-3.50006pt\Bigg{\rangle}}_{\mathbf{x}(t)}$

(44)

This expression for the torque has the intuitive interpretation where the contribution from each surface element is the cross-product of the displacement vector from the rotation center to the surface point with the Casimir force density associated with that surface element.

Fig. 4 demonstrates the utility of the Hermite-weighting method of Eq. (11) for computing derivatives of path integrals. The plot shows the results of the numerical evaluation of the path integral (4) for the Casimir–Polder potential of an atom near a dielectric half-space, along with the computed spatial derivatives up to the tenth. Since the potential and its derivatives are normalized to the potential in the perfect-conductor limit, the analytic result is distance-independent, and has the same form in each case, given by the expression

$$\eta_{\mathrm{\scriptscriptstyle TE}}(\chi)=\frac{1}{6}+\frac{1}{\chi}-\frac{\sqrt{1+\chi}}{2\chi}-\frac{\sinh^{-1}\!\sqrt{\chi}}{2\chi^{3/2}},$$

(45)

as discussed in Ref. [19] (while not completely realistic [34], this model suffices for a numerical test).

Some aspects of the calculations in this figure are worth discussing here. These computations nominally employed $N=10^{5}$ points per path. However, at finite $N$, the error analysis in Ref. [19] indicates that numerical accuracy suffers, particularly at large values of the susceptibility $\chi$, because the finitely sampled path does not have the same “extent” as the same path in the continuum limit. To help maintain accuracy at large $\chi$, the paths employ some limited, finer sampling: the two path intervals on either side of the point closest to the dielectric interface are each sampled with $10^{3}$ more points (see App. A for the algorithm to generate these sub-paths). This strategy substantially alleviates the major source of sampling error at large $\chi$, while negligibly impacting the computational effort. In terms of the number $m=k$ of partial averages, we make the simplistic choice of choosing a fixed value of $m$ for each calculation ($m=10$, 20, 200, 200, and 500 partial averages were employed for $n=1$, 2, 4, 7, and 10 derivatives, respectively), despite the indication of Eq. (18) that the number of partial averages should be chosen to be inversely proportional to $\mathcal{T}$ to maintain a controlled error. This is justified by noting that the partial averages can be chosen to have small errors for small values of $\mathcal{T}$, and the effects of the resulting larger error for large values of $\mathcal{T}$ will be suppressed by the factor $\mathcal{T}^{-3}$ in the path integral (4). (For high-accuracy evaluation of the path integral, the partial averages should instead be chosen in a $\mathcal{T}$-dependent way.)

Another major numerical challenge that arises in Fig. 4 in the case of high-order derivatives comes from the Hermite-polynomial factor in Eq. (11). One issue here is that the Hermite polynomial takes on large values for large values of $|\bar{x}_{m}-x_{0}|$. But the Gaussian path measure in the path integral only generates such large values rarely. This means that, with ordinary Gaussian paths, rare events make large contributions to the ensemble, leading to poor statistical accuracy in the numerical computations. In these calculations, we employed the direct solution of sampling $\bar{x}_{m}$ from the Hermite–Gaussian density

$$P_{H_{n}}(\bar{x}_{m})=\frac{\eta_{n}}{\sqrt{\pi m\Delta\mathcal{T}}}\left|H_{n}\!\left(\frac{\bar{x}_{m}-x_{0}}{\sqrt{m\Delta\mathcal{T}}}\right)\right|e^{-(\bar{x}_{m}-x_{0})^{2}/m\Delta\mathcal{T}},$$

(46)

where the factor

$$\eta_{n}^{-1}:=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\!\!dx\,\left|H_{n}(x)\right|\,e^{-x^{2}}$$

(47)

normalizes the probability distribution. Then a Gaussian deviate $\delta x_{m}$ of zero mean and variance $m(1-2m/N)\Delta\mathcal{T}/2$ allows the first and last path steps to be generated according to $x_{m}=\bar{x}_{m}+\delta x_{m}$ and $x_{N-m}=\bar{x}_{m}-\delta x_{m}$, and the remainder of the path is generated as a Brownian bridge running from $x_{m}$ to $x_{N-m}$ in time $(N-2m)\Delta\mathcal{T}$ (see App. A for the bridge-generating algorithm). With these modified paths, the path average (11) should be modified to read

$\displaystyle\frac{\partial^{n}}{\partial x_{0}^{n}}I$

$\displaystyle\approx\eta_{n}^{-1}(m\Delta\mathcal{T})^{-n/2}\,{\bigg{\langle}\kern-3.50006pt\bigg{\langle}}{\Phi[\mathbf{x}(t)]\,\mathrm{sgn}\big{[}H_{n}(\bar{z})\big{]}}{\bigg{\rangle}\kern-3.50006pt\bigg{\rangle}}_{H_{n}},$

(48)

where $\bar{z}:=(\bar{x}_{m}-x_{0})/\sqrt{m\Delta\mathcal{T}}$. Because the Hermite polynomial can take either sign, a great deal of cancellation occurs in Eq. (48) among the paths, which can again lead to poor statistical performance for high-order derivatives, when the Hermite polynomial has multiple regions of opposing sign. To combat this, for a given path, the path functional can be averaged over multiple values of $\bar{x}_{m}$ (holding the rest of the path fixed) to effect this cancellation on a pathwise basis. The calculations in Fig. 4, used averages over 1, 5, 100, 5000, and 5000 values of $\bar{x}_{m}$ per path for $n=1$, 2, 4, 7, and 10 derivatives, respectively, where the $\bar{x}_{m}$ were sampled from a sub-random (1D van der Corput) sequence to enhance the accuracy of this subaverage.

The numerical performance in Fig. 4 is generally quite good, although performance suffers for high-order derivatives. All of the data here represent statistical averages over $10^{8}$ paths (except for $n=10$ derivatives, which averages over $10^{9}$ paths). The points at a given derivative order $n$ were computed using the same random numbers, so the points within a curve are not statistically independent. The error bars (delimiting one standard error) are invisibly small on the plot, and amount to a relative statistical error of at most $0.2\%$, $0.5\%$, $0.7\%$, $0.8\%$, $2\%$, and $2\%$ for $n=0$, 1, 2, 4, 7, and 10 derivatives, respectively. The observed accuracy of the results is comparable to the standard error in each case except for the 10th derivative, where the observed error is at worst 14%, despite the extra computational effort. However, note that a comparable computation of the 10th derivative via the finite-difference method would be completely unfeasible. For an effective resolution of $N=10^{8}$ points per path (with $N=10^{5}$ and subsampling of $10^{3}$ points as described above), a length scale of the order $N^{-1/2}=10^{-4}$ is appropriate for the finite difference, but this leads to a standard deviation in the path values of the order of $10^{40}$ relative to the mean.

To study the numerical error of the partial-averaging method in more detail, we performed the same calculations as in the previous section for one and two spatial derivatives of the Casimir–Polder potential for an atom near a planar interface, while varying the number $m$ of partial averages (Figs. 5 and 6). To simplify the analysis, $m$ is held constant (independent of $\mathcal{T}$), ignoring the error-control criterion (18). For comparison, the plots also show the corresponding errors when computing the derivatives via finite differences. The first derivative employs the second-order, centered-difference formula $V^{\prime}_{\mathrm{\scriptscriptstyle CP}}(d)=[V_{\mathrm{\scriptscriptstyle CP}}(d+\delta)-V_{\mathrm{\scriptscriptstyle CP}}(d-\delta)]/2\delta+O(\delta^{2})$, where $\delta$ is a small spatial displacement. The second derivative exploits the analogous difference formula $V^{\prime\prime}_{\mathrm{\scriptscriptstyle CP}}(d)=[V_{\mathrm{\scriptscriptstyle CP}}(d+\delta)-2V_{\mathrm{\scriptscriptstyle CP}}(d)+V_{\mathrm{\scriptscriptstyle CP}}(d-\delta)]/\delta^{2}+O(\delta^{2})$. In implementing the finite-difference methods, it is important to apply these difference formulae on a path-wise basis (or equivalently, apply the difference formulae on ensemble averages over equivalent path samples) for this method to produce any kind of reasonable result. In applying the finite-difference method, the numerical details (path generation, calculation of path averages, etc.) otherwise match those of previous work (see Ref. [19], Section V). Although the error for both the partial-averaging and finite-difference methods behaves similarly in the plots, note that the quantities $m/N$ and $\delta/d$ from the two methods are not directly comparable. In the partial-averaging method, $m$ sets a length scale $\sqrt{m\Delta\mathcal{T}}=\sqrt{(m/N)\mathcal{T}}$ that varies with $\mathcal{T}$ in the integral in Eq. (4). However, it is customary practice for the finite-difference method to use a displacement $\delta$ independent of the details of the rest of the calculation, so $\delta$ represents a fixed length scale, independent of $\mathcal{T}$. Nevertheless, the plots provide a useful visual comparison of the performance of the two methods.