text:
# Cost of the Generalised Hybrid Monte Carlo Algorithm for Free Field Theory

## 1 Introduction

Hybrid Monte Carlo (HMC)  remains the most popular algorithm for simulation of Quantum Chromodynamics (QCD) including the dynamical effects of fermions. It is therefore imperative that we have a detailed understanding of how the computational cost of HMC varies as we change simulation parameters such as the lattice volume and the fermion mass. In this paper we present detailed analytical results for the computational cost of the generalised HMC algorithm for free field theory. We expect many of the results obtained to also be applicable to more general field theories where most of the “modes” are weakly interacting, for example the high momentum modes of asymptotically free theories such as QCD.

In this paper we give detailed descriptions of the techniques we have developed to enable us to calculate fictitious-time evolution operators, Metropolis acceptance probabilities, and autocorrelation functions for arbitrary polynomial operators in GHMC simulations of free field theory. This allows us to minimise the computational cost of such simulations, and to compare the efficiency of various popular limits of GHMC. This paper brings together and extends many results that have been presented previously . The techniques developed for this paper have also been used in several other papers .

The structure of this paper is as follows: in Section 2 we describe the Generalised Hybrid Monte Carlo (GHMC) Algorithm and discuss various limiting cases. Section 3 demonstrates that GHMC computations for free field theory are equivalent to GHMC computations for a set of uncoupled harmonic operators with judiciously chosen frequencies. In Section 4 we describe the generalised leapfrog discretisation schemes for the classical equations of motion in fictitious time, and introduce specific parameterisations of the evolution operators in order to facilitate the calculation of the Metropolis acceptance probability in Section 5. In Section 6 we introduce general techniques for calculating Laplace transforms of autocorrelation functions of polynomial operators in free field theory, and present explicit results for arbitrary linear and quadratic operators for several of the limiting cases of GHMC, both for fixed- and exponentially-distributed GHMC trajectory lengths. Section 7 presents a detailed analysis of the cost of HMC simulations in the approximation that the Metropolis acceptance rate is unity, whilst Section 8 analyses the relative cost of GHMC, HMC, and the second order Langevin algorithm (L2MC) for non-trivial acceptance rates under some very weak assumptions. Finally Section 9 discusses the somewhat related topic of how to optimise relative frequencies of updates and measurements in general Monte Carlo simulations. Some concluding remarks are contained in Section 10. Several detailed technical results are relegated to Appendices.

## 2 Generalised Hybrid Monte Carlo Algorithm

For simplicity, we shall describe the Generalised Hybrid Monte Carlo (GHMC) algorithm for a theory of scalar fields, denoted generically by $`\varphi `$, with action $`S(\varphi )`$.

In order to generate the distribution

$$\frac{1}{Z}e^{S(\varphi )}[d\varphi ]$$

we introduce a set of fictitious momenta $`\pi `$ and generate the fields $`\varphi `$ and $`\pi `$ according to the distribution

$$\frac{1}{Z^{}}e^{H(\varphi ,\pi )}[d\varphi ][d\pi ]$$
(1)
where

$$H(\varphi ,\pi )=\frac{1}{2}\underset{x}{}\pi _x^2+S(\varphi ),$$

and ignore $`\pi `$.

We begin by recalling  that a Markov Process will converge to some distribution of configurations if it is constructed out of update steps each of which has the desired distribution as a fixed point, and which taken together are ergodic. The generalised HMC algorithm  is constructed out of two such steps.

### 2.1 Molecular Dynamics Monte Carlo

This in turn consists of three parts:

1. *Molecular Dynamics (MD):* an approximate integration of Hamilton’s equations on phase space which is exactly area-preserving and reversible; this is a map $`U(\tau ):(\varphi ,\pi )(\varphi ^{},\pi ^{})`$ with $`det\frac{(\varphi ^{},\pi ^{})}{(\varphi ,\pi )}=1`$ and $`U(\tau )=U^1(\tau )`$.
2. A *momentum flip* $`F:\pi \pi `$.
3. *Monte Carlo (MC):* a Metropolis accept/reject test

$$(\varphi ^{},\pi ^{})=\{\begin{array}{cc}\hfill FU(\varphi ,\pi )& \text{with probability }\mathrm{min}(1,e^{\delta H})\hfill \\ \hfill (\varphi ,\pi )& \text{otherwise.}\hfill \end{array}$$

This may be implemented by generating a uniformly distributed random number $`y[0,1]`$, so that

$$_0^1𝑑y\theta \left(e^{\delta H}y\right)=\mathrm{min}(1,e^{\delta H}).$$

The composition of these gives the MDMC update step

$$\left(\begin{array}{c}\varphi ^{}\\ \pi ^{}\end{array}\right)=\left[FU(\tau )\theta (e^{\delta H}y)+\mathrm{𝕀}\theta (ye^{\delta H})\right]\left(\begin{array}{c}\varphi \\ \pi \end{array}\right).$$
(2)
This satisfies detailed balance because $`(FU)^2=\mathrm{𝕀}`$.

### 2.2 Partial Momentum Refreshment

This mixes the Gaussian-distributed momenta $`\pi `$ with Gaussian noise $`\xi `$

$$\left(\begin{array}{c}\pi ^{}\\ \xi ^{}\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\vartheta & \mathrm{sin}\vartheta \\ \mathrm{sin}\vartheta & \mathrm{cos}\vartheta \end{array}\right)F\left(\begin{array}{c}\pi \\ \xi \end{array}\right).$$
(3)
If $`\pi `$ and $`\xi `$ are Gaussian distributed, $`P_G(\pi )[d\pi ]e^{\pi ^2/2}[d\pi ]`$, then so is $`P_G(\pi ^{})`$:

$$P_G(\pi ^{})=_{\mathrm{}}^{\mathrm{}}𝑑\xi _{\mathrm{}}^{\mathrm{}}𝑑\pi P_G(\xi )P_G(\pi )\delta (\pi ^{}+\pi \mathrm{cos}\vartheta \xi \mathrm{sin}\vartheta ).$$

The extra momentum flip $`F`$ is included so that the trajectory is reversed upon an MC rejection instead of on an acceptance.

We may combine the MDMC update of equation (2) with the partial momentum refreshment of equation (3) to obtain the following expression for their combined effect

$$\left(\begin{array}{c}\varphi ^{}\\ \pi ^{}\\ \xi ^{}\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\vartheta & \mathrm{sin}\vartheta \\ 0& \mathrm{sin}\vartheta & \mathrm{cos}\vartheta \end{array}\right)\left[U(\tau )\theta (e^{\delta H}y)+F\theta (ye^{\delta H})\right]\left(\begin{array}{c}\varphi \\ \pi \\ \xi \end{array}\right).$$
(4)
### 2.3 Special Cases of GHMC

Several well-known algorithms are special cases of GHMC:

* The usual Hybrid Monte Carlo (HMC) algorithm is the special case where $`\vartheta =\pi /2`$, i.e., the fictitious momenta $`p`$ are replaced by the Gaussian distributed $`\xi `$, the old momenta being discarded completely. The momentum flips may be ignored in this case as long as MCMD and momentum refreshment steps are strictly alternated.
* $`\vartheta =0`$ corresponds to pure MDMC, which is an exact version of the MD or microcanonical algorithm. It is in general non-ergodic.
* The second order Langevin Monte Carlo (L2MC) algorithm of Horowitz  corresponds to choosing arbitrary $`\vartheta `$ but MDMC trajectories of a single leapfrog integration step, $`\tau =\delta \tau `$. This is an exact version of the second order Langevin algorithm .
* The Langevin Monte Carlo (LMC)  algorithm has $`\vartheta =\pi /2`$ and $`\tau =\delta \tau `$. This is an exact version of the Langevin algorithm.

### 2.4 Tunable Parameters

The GHMC algorithm has three free parameters, the trajectory length $`\tau `$ the momentum mixing angle $`\vartheta `$, and the integration step size $`\delta \tau `$. These may be chosen arbitrarily without affecting the validity of the method, except for some special values for which the algorithm ceases to be ergodic. We may adjust these parameters to minimise the cost of a Monte Carlo computation, and the main goal of this work is to carry out this optimisation procedure for free field theory.

Horowitz  pointed out that the L2MC algorithm has the advantage of having a higher acceptance rate than HMC for a given step size, but he did not take in to account that it also requires a higher acceptance rate to get the same autocorrelations because the trajectory is reversed at each MC rejection. It is not obvious a priori which of these effects dominates.

The parameters $`\tau `$ and $`\vartheta `$ may be chosen independently from some distributions $`P_R(\tau )`$ and $`P_M(\vartheta )`$ for each trajectory (but of course they cannot be chosen in a way which is correlated with the starting point in phase space). In the following we shall consider various choices for the momentum refreshment distribution $`P_M`$, but we shall always take a fixed value for $`\vartheta `$, $`P_M(\vartheta )=\delta (\vartheta \vartheta _0)`$. The generalisation of our results to the case of other distributions of $`\vartheta `$ values is trivial, but it is not immediately obvious that it is useful. We choose each trajectory $`\tau `$ length independently from some distribution $`P_R(\tau )`$, as this avoids the lack of ergodicity caused by choosing a fixed trajectory length which is a rational multiple of the period of any mode of the system . This is a disease of free field theory which in interacting models is removed to some extent by mode coupling.

## 3 Free Field Theory

### 3.1 Complex Fields on a Finite Lattice

Consider a $`d`$ dimensional free scalar field theory on a $`VL^d`$ lattice for which the expectation of an arbitrary operator $`\mathrm{\Omega }(\varphi )`$ is defined by the functional integral

$$\mathrm{\Omega }\frac{1}{Z}[d\varphi ]e^{S(\varphi )}\mathrm{\Omega }(\varphi )$$

where $`Z`$ is chosen such that $`1=1`$. For a complex field $`\varphi :\mathrm{}_L^d\mathrm{}`$ the action is

$$S(\varphi )=\frac{1}{2}\underset{x\mathrm{}_L^d}{}\varphi _x^{}\left(^2+m^2\right)\varphi _x$$

where $`^2`$ is the lattice Laplacian

$$^2\varphi _x\underset{\mu =1}{\overset{d}{}}(\varphi _{x+\widehat{\mu }}2\varphi _x+\varphi _{x\widehat{\mu }}).$$

For the Generalised Hybrid Monte Carlo (GHMC) algorithm discussed in section 2 we introduce a set of fictitious momentum fields $`\pi `$, and the corresponding Hamiltonian

$$H(\varphi ,\pi )\frac{1}{2}\underset{x\mathrm{}_L^d}{}\pi _x^{}\pi _x+S(\varphi ),$$

and we evaluate the functional integral

$$\mathrm{\Omega }\frac{1}{Z^{}}[d\varphi ][d\pi ]e^{H(\varphi ,\pi )}\mathrm{\Omega }(\varphi ).$$

For theoretical analysis it is useful to diagonalise the Hamiltonian by Fourier transformation to “real” (as opposed to “fictitious”) momentum space. From the identity of example (B.1.1) we have

$$\begin{array}{ccccc}\hfill \stackrel{~}{\varphi }_p& & (\varphi )_p& & L^{d/2}\underset{x\mathrm{}_L^d}{}e^{2\pi ipx/L}\varphi _x\hfill \\ \hfill \varphi _x& =& (^1\stackrel{~}{\varphi })_x& =& L^{d/2}\underset{p\mathrm{}_L^d}{}e^{2\pi ipx/L}\stackrel{~}{\varphi }_p\hfill \end{array}\}$$
(5)
and similarly for the fictitious momenta. We obtain

$$H(\varphi ,\pi )=\stackrel{~}{H}(\stackrel{~}{\varphi },\stackrel{~}{\pi })=\frac{1}{2}\underset{p\mathrm{}_L^d}{}\left(\stackrel{~}{\pi }_p^{}\stackrel{~}{\pi }_p+\omega _p^2\stackrel{~}{\varphi }_p^{}\stackrel{~}{\varphi }_p\right),$$
(6)
where the frequencies are

$$\omega _p^2m^2+4\underset{\mu =1}{\overset{d}{}}\mathrm{sin}^2\frac{\pi p_\mu }{L}.$$
(7)
The phase space configuration $`(\varphi ,\pi )`$ is generated by GHMC with probability density proportional to

$`e^{H(\varphi ,\pi )}[d\varphi ][d\pi ]`$ $`=`$ $`e^{\stackrel{~}{H}(\stackrel{~}{\varphi },\stackrel{~}{\pi })}det\left({\displaystyle \frac{\delta \varphi }{\delta \stackrel{~}{\varphi }}}\right)det\left({\displaystyle \frac{\delta \pi }{\delta \stackrel{~}{\pi }}}\right)[d\stackrel{~}{\varphi }][d\stackrel{~}{\pi }]`$
$``$ $`e^{\stackrel{~}{H}(\stackrel{~}{\varphi },\stackrel{~}{\pi })}[d\stackrel{~}{\varphi }][d\stackrel{~}{\pi }],`$
where the last relation follows because the Jacobian is a (real) constant. We thus see that free field theory corresponds to a set of harmonic oscillators with a specific choice of frequency spectrum.

### 3.2 Infinite Lattices

For an infinite lattice $`\mathrm{}^d`$ (instead of $`\mathrm{}_L^d`$) the corresponding momentum space is the torus $`S_1^d`$ (instead of $`\mathrm{}_L^d`$), and equations (5), (6), and (7) get replaced by

$$\begin{array}{ccccc}\hfill \stackrel{~}{\varphi }_p& & (\varphi )_p& & (2\pi )^{d/2}\underset{x\mathrm{}^d}{}e^{ipx}\varphi _x\hfill \\ \hfill \varphi _x& =& (^1\stackrel{~}{\varphi })_x& =& (2\pi )^{d/2}_{S_1^d}𝑑pe^{ipx}\stackrel{~}{\varphi }_p\hfill \end{array}\},$$

$$H(\varphi ,\pi )=\stackrel{~}{H}(\stackrel{~}{\varphi },\stackrel{~}{\pi })=\frac{1}{2}_{S_1^d}𝑑p\left(\stackrel{~}{\pi }_p^{}\stackrel{~}{\pi }_p+\omega _p^2\stackrel{~}{\varphi }_p^{}\stackrel{~}{\varphi }_p\right),$$

and

$$\omega (p)^2m^2+4\underset{\mu =1}{\overset{d}{}}\mathrm{sin}^2\frac{p_\mu }{2}.$$

where we made use of examples (B.1.2) and (B.1.3).

### 3.3 Real Fields

For real (as opposed to complex) fields $`\varphi :\mathrm{}_L^d\mathrm{}`$ we define

$`\stackrel{~}{\varphi }_p^{(e)}`$ $``$ $`Re\stackrel{~}{\varphi }_p`$
$`\stackrel{~}{\varphi }_p^{(o)}`$ $``$ $`Im\stackrel{~}{\varphi }_p`$
Since $`\stackrel{~}{\varphi }_p^{(e)}=\stackrel{~}{\varphi }_p^{(e)}`$ and $`\stackrel{~}{\varphi }_p^{(o)}=\stackrel{~}{\varphi }_p^{(o)}`$ these two fields are independent degrees of freedom only on a subset of the momentum space $`p\mathrm{}_L^d`$. We may choose to define the real momentum field $`\stackrel{~}{\varphi }:\mathrm{}_L^d\mathrm{}`$ as

$$\stackrel{~}{\varphi }_p=\sqrt{2}\{\begin{array}{ccc}\stackrel{~}{\varphi }_p^{(e)}\hfill & & \text{if }pp\hfill \\ \stackrel{~}{\varphi }_p^{(o)}\hfill & & \text{if }p<p\hfill \end{array}$$

where the ordering is arbitrary (*e.g.*, lexicographic). The Hamiltonian is then

$`\stackrel{~}{H}(\stackrel{~}{\varphi },\stackrel{~}{\pi })`$ $`=`$ $`\frac{1}{2}{\displaystyle \underset{p\mathrm{}_L^d}{}}\left(\left|\stackrel{~}{\pi }_p^{(e)}+i\stackrel{~}{\pi }_p^{(o)}\right|^2+\omega _p^2\left|\stackrel{~}{\varphi }_p^{(e)}+i\stackrel{~}{\varphi }_p^{(o)}\right|^2\right)`$
$`=`$ $`\frac{1}{2}{\displaystyle \underset{p\mathrm{}_L^d}{}}\left(\stackrel{~}{\pi }_p^{(e)\mathrm{\hspace{0.17em}2}}+\stackrel{~}{\pi }_p^{(o)\mathrm{\hspace{0.17em}2}}+\omega _p^2(\stackrel{~}{\pi }_p^{(e)\mathrm{\hspace{0.17em}2}}+\stackrel{~}{\pi }_p^{(o)\mathrm{\hspace{0.17em}2}})\right)`$
$`=`$ $`\frac{1}{2}{\displaystyle \underset{p\mathrm{}_L^d}{}}\left(\stackrel{~}{\pi }_p^2+\omega _p^2\stackrel{~}{\varphi }_p^2\right)`$
because $`\omega _p^2=\omega _p^2`$.

### 3.4 Spectral Averages

Results obtained for a general set of harmonic oscillators will be expressed in terms of “spectral averages,” which may be evaluated explicitly for free field theory. We use the notation $`𝜎(f)`$ to denote the spectral average of some function $`f:S_1^d\mathrm{}`$, $`𝜎(f)L^d_{p\mathrm{}_L^d}f\left(\frac{2\pi p}{L}\right)`$.

The finite lattice results may further be expanded about the infinite lattice results to obtain a large volume expansion by means of the Poisson resummation formula (92) of section (B.2),

$$𝜎(f)=\underset{k\mathrm{}}{}\left[\underset{\mu =1}{\overset{d}{}}\left(\frac{1}{2\pi }\right)_0^{2\pi }𝑑p_\mu \mathrm{cos}(p_\mu kL)\right]f(p).$$

###### Example 3.4.1

Consider the massless spectral average in one dimension (where the volume $`VL^d=L`$)

$`𝜎\left(\omega (p)^2\right)`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega \left({\displaystyle \frac{2\pi p}{V}}\right)^{2\alpha }={\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^{2\alpha }={\displaystyle \frac{4^\alpha }{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{sin}^{2\alpha }{\displaystyle \frac{\pi p}{V}}`$
$`=`$ $`{\displaystyle \frac{4^\alpha }{2\pi }}{\displaystyle _0^{2\pi }}𝑑p\mathrm{sin}^{2\alpha }{\displaystyle \frac{p}{2}}+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{4^\alpha }{\pi }}{\displaystyle _0^{2\pi }}𝑑p\mathrm{sin}^{2\alpha }{\displaystyle \frac{p}{2}}\mathrm{cos}pkV`$
$`=`$ $`{\displaystyle \frac{4^\alpha }{\pi }}B(\frac{1}{2},\alpha +\frac{1}{2})={\displaystyle \frac{1}{\pi }}\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right),`$
where all the higher order terms vanish for large $`V`$ because $`_0^{2\pi }𝑑p\mathrm{sin}^{2\alpha }\frac{p}{2}\mathrm{cos}pkV=0`$ for $`kV>\alpha `$ upon integration by parts. This result is to be compared with the exact answer given in equation (83) of Appendix B.

## 4 Leapfrog Evolution

### 4.1 Lowest Order Leapfrog

We wish to consider a system of $`V`$ uncoupled harmonic oscillators $`\{\varphi _p\}`$ for $`p\mathrm{}_V`$. The Hybrid Monte Carlo algorithm requires us to introduce a corresponding set of “fictitious” momenta $`\{\pi _p\}`$, and the dynamics on this “fictitious” phase space is described by the Hamiltonian

$$H=\frac{1}{2}\underset{p=1}{\overset{V}{}}\left(\pi _p^2+\omega _p^2\varphi _p^2\right).$$
(8)
The classical trajectory through phase space must obey Hamilton’s equations (using the symplectic 2-form $`_{i=1}^Nd\pi _pd\varphi _p`$)

$$\dot{\varphi }_p=\frac{H}{\pi _p}=\pi _p,\dot{\pi }_p=\frac{H}{\varphi _p}=\omega _p^2\varphi _p.$$

The leapfrog equations are the simplest discretisation of these which are exactly reversible and area-preserving,

$`\pi _p(\frac{1}{2}\delta \tau )=`$ $`\pi _p(0)+\dot{\pi }_p(0)\frac{1}{2}\delta \tau `$ $`=\pi _p(0)\omega _p^2\varphi _p(0)\frac{1}{2}\delta \tau `$
$`\varphi _p(\delta \tau )=`$ $`\varphi _p(0)+\dot{\varphi }_p(\frac{1}{2}\delta \tau )\delta \tau `$ $`=\varphi _p(0)+\pi _p(\frac{1}{2}\delta \tau )\delta \tau `$
$`\pi _p(\delta \tau )=`$ $`\pi _p(\frac{1}{2}\delta \tau )+\dot{\pi }_p(\delta \tau )\frac{1}{2}\delta \tau `$ $`=\pi _p(\frac{1}{2}\delta \tau )\omega _p^2\varphi _p(\delta \tau )\frac{1}{2}\delta \tau .`$
This leapfrog integration scheme is a linear mapping on phase space111For the rest of this section we shall consider only a single oscillator with $`\omega _p=1`$, as everything will be diagonal in $`p`$ and the frequency dependence can be recovered by dimensional analysis.

$$\left(\begin{array}{c}\varphi (\tau +\delta \tau )\\ \pi (\tau +\delta \tau )\end{array}\right)=U_0(\delta \tau )\left(\begin{array}{c}\varphi (\tau )\\ \pi (\tau )\end{array}\right)$$

where the matrix

$$U_0(\delta \tau )\left(\begin{array}{cc}1{\scriptscriptstyle \frac{1}{2}}\delta \tau ^2& \delta \tau \\ \delta \tau +{\scriptscriptstyle \frac{1}{4}}\delta \tau ^3& 1{\scriptscriptstyle \frac{1}{2}}\delta \tau ^2\end{array}\right)$$
(9)
satisfies $`detU_0=1`$, as it must because the mapping is area-preserving. This lowest-order leapfrog integration agrees with the exact Hamiltonian evolution up to errors of $`O(\delta \tau ^3)`$,

$$U_0(\delta \tau )=e^{\widehat{H}\delta \tau }\left[\mathrm{𝕀}+U_{0,1}\delta \tau ^3+O(\delta \tau ^5)\right]$$
(10)
where

$$\widehat{H}\left(\begin{array}{c}\varphi \\ \pi \end{array}\right)\left[\frac{H}{\pi }\frac{}{\varphi }\frac{H}{\varphi }\frac{}{\pi }\right]\left(\begin{array}{c}\varphi \\ \pi \end{array}\right)=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}\varphi \\ \pi \end{array}\right)$$
(11)
is the generator of a translation through fictitious time, and $`U_{0,1}`$ is some operator on phase space. The error must be an odd function of $`\delta \tau `$ because leapfrog integration is reversible to all orders in $`\delta \tau `$.

### 4.2 Higher Order Leapfrog

Following Campostrini *et al.*  we can easily construct a higher-order leapfrog integration schemes with errors of $`O(\delta \tau ^5)`$ by defining a “wiggle”

$$U_1(\delta \tau )U_0\left(\frac{\delta \tau }{\alpha _1}\right)U_0\left(\frac{\sigma _1\delta \tau }{\alpha _1}\right)U_0\left(\frac{\delta \tau }{\alpha _1}\right).$$

Clearly this is area-preserving and reversible because $`U_0`$ is. Using equation (10), we find that

$$U_1(\delta \tau )=\mathrm{exp}\left(\frac{2\sigma _1}{\alpha _1}\widehat{H}\delta \tau \right)\left[\mathrm{𝕀}+U_{0,1}\left(2\sigma _1^3\right)\left(\frac{\delta \tau }{\alpha _1}\right)^3+O(\delta \tau ^5)\right];$$

if we choose $`\sigma _1=\sqrt[3]{2}`$ then the coefficient of $`\delta \tau ^3`$ vanishes, and we can arrange the step size to equal that of the lowest-order leapfrog scheme by taking $`\alpha _1=2\sigma _1`$. The explicit form for the first-order wiggle $`U_1(\delta \tau )`$ is

$$\left(\begin{array}{cc}1\frac{1}{2}\delta \tau ^2+\frac{1}{24}\delta \tau ^4+\frac{6+5\sqrt[3]{2}+4\sqrt[3]{4}}{288}\delta \tau ^6& \delta \tau \frac{1}{6}\delta \tau ^3\frac{4+4\sqrt[3]{2}+3\sqrt[3]{4}}{144}\delta \tau ^5\\ \delta \tau +\frac{1}{6}\delta \tau ^3+\frac{\sqrt[3]{2}+\sqrt[3]{4}}{144}\delta \tau ^5\frac{25+20\sqrt[3]{2}+16\sqrt[3]{4}}{1728}\delta \tau ^7& 1\frac{1}{2}\delta \tau ^2+\frac{1}{24}\delta \tau ^4+\frac{6+5\sqrt[3]{2}+4\sqrt[3]{4}}{288}\delta \tau ^6\end{array}\right).$$

This construction can be iterated by defining

$$U_n(\delta \tau )U_{n1}\left(\frac{\delta \tau }{\alpha _n}\right)U_{n1}\left(\frac{\sigma _n\delta \tau }{\alpha _n}\right)U_{n1}\left(\frac{\delta \tau }{\alpha _n}\right),$$
(12)
which again clearly is area-preserving and reversible. This may be shown to have errors of the form

$$U_n(\delta \tau )=e^{\widehat{H}\delta \tau }\left[\mathrm{𝕀}+U_{n,1}\delta \tau ^{2n+3}+O(\delta \tau ^{2n+5})\right]$$
(13)
by induction on $`n`$: Assume equation (12) holds $`n<N`$, then from equations (12) and (13) we find that

$$U_N(\delta \tau )=\mathrm{exp}\left(\frac{2\sigma _N}{\alpha _N}\widehat{H}\delta \tau \right)\left[\mathrm{𝕀}+U_{N1,1}\left(2\sigma _N^{2N+1}\right)\left(\frac{\delta \tau }{\alpha _N}\right)^{2N+1}+O(\delta \tau ^{2N+3})\right],$$

which gives us equation (13) for $`n=N`$ upon setting $`\sigma _N=\sqrt[2N+1]{2}`$ and $`\alpha _N=2\sigma _N`$.

### 4.3 Parameterisation of Leapfrog Evolution Operators

In order to calculate the explicit form for a Molecular Dynamics trajectory consisting of $`\tau /\delta \tau `$ leapfrog steps it is useful to parameterise the leapfrog matrices $`U_n(\delta \tau )`$ in the following way. Reversibility requires that for any leapfrog matrix $`U_n(\delta \tau )=U_n(\delta \tau )^1`$ and area-preservation requires that $`detU_n(\delta \tau )=1`$, hence

$$U_n(\delta \tau )=\left(\begin{array}{cc}U_{n_{1,1}}(\delta \tau )& U_{n_{1,2}}(\delta \tau )\\ U_{n_{2,1}}(\delta \tau )& U_{n_{2,2}}(\delta \tau )\end{array}\right)=\left(\begin{array}{cc}U_{n_{2,2}}(\delta \tau )& U_{n_{1,2}}(\delta \tau )\\ U_{n_{2,1}}(\delta \tau )& U_{n_{1,1}}(\delta \tau )\end{array}\right).$$

For the lowest-order leapfrog matrix of equation (9) we observe that the diagonal elements are even functions of $`\delta \tau `$ and the off-diagonal elements are odd functions of $`\delta \tau `$; this property also holds for all the higher-order leapfrog matrices defined by equation (12). We may therefore parameterise $`U_n`$ in terms of two even functions

$$\begin{array}{ccc}\hfill \kappa _n(\delta \tau )& =& 1+\kappa _{n,1}\delta \tau ^{2n+2}+O(\delta \tau ^{2n+4})\hfill \\ \hfill \rho _n(\delta \tau )& =& 1+\rho _{n,1}\delta \tau ^{2n+2}+O(\delta \tau ^{2n+4}).\hfill \end{array}$$
(14)
as

$$U_n(\delta \tau )=\left(\begin{array}{cc}\mathrm{cos}[\kappa _n(\delta \tau )\delta \tau ]& \frac{\mathrm{sin}\left[\kappa _n\left(\delta \tau \right)\delta \tau \right]}{\rho _n\left(\delta \tau \right)}\\ \rho _n(\delta \tau )\mathrm{sin}[\kappa _n(\delta \tau )\delta \tau ]& \mathrm{cos}[\kappa _n(\delta \tau )\delta \tau ]\end{array}\right).$$
(15)
We may easily compute the leading terms in the Taylor expansions of these functions. For the lowest-order leapfrog scheme we obtain

$`\kappa _0`$ $`=`$ $`1+{\displaystyle \frac{1}{24}}\delta \tau ^2+{\displaystyle \frac{3}{640}}\delta \tau ^4+O(\delta \tau ^6)=1+0.041\dot{6}\delta \tau ^2+O(\delta \tau ^4)`$
$`\rho _0`$ $`=`$ $`1{\displaystyle \frac{1}{8}}\delta \tau ^2{\displaystyle \frac{1}{128}}\delta \tau ^4+O(\delta \tau ^6)=10.125\delta \tau ^2+O(\delta \tau ^4);`$
for the Campostrini “wiggle”

$`\kappa _1`$ $`=`$ $`1{\displaystyle \frac{\left(32+25\sqrt[3]{2}+20\sqrt[3]{4}\right)}{1440}}\delta \tau ^4+O(\delta \tau ^6)10.06614\delta \tau ^4+O(\delta \tau ^6)`$
$`\rho _1`$ $`=`$ $`1+{\displaystyle \frac{\left(4+3\sqrt[3]{2}+2\sqrt[3]{4}\right)}{288}}\delta \tau ^4+O(\delta \tau ^6)1+0.03804\delta \tau ^4+O(\delta \tau ^6);`$
and for the second-order “wiggle”

$`\kappa _2`$ $`=`$ $`1+{\displaystyle \frac{\left(\begin{array}{c}z\left(1321294+1039220\sqrt[3]{2}+831376\sqrt[3]{4}\right)\\ +\left(1146320+901600\sqrt[3]{2}+721280\sqrt[3]{4}\right)\sqrt[5]{2}\\ +\left(958174+753620\sqrt[3]{2}+602896\sqrt[3]{4}\right)\sqrt[5]{4}\\ +\left(811769+638470\sqrt[3]{2}+510776\sqrt[3]{4}\right)\sqrt[5]{8}\\ +\left(726418+571340\sqrt[3]{2}+457072\sqrt[3]{4}\right)\sqrt[5]{16}\end{array}\right)}{881798400}}\delta \tau ^6+O(\delta \tau ^8)1+0.02169\delta \tau ^6+O(\delta \tau ^8)`$
$`\rho _2`$ $`=`$ $`1{\displaystyle \frac{\left(\begin{array}{c}z\left(391582+312244\sqrt[3]{2}+247752\sqrt[3]{4}\right)\\ +\left(348752+277664\sqrt[3]{2}+219456\sqrt[3]{4}\right)\sqrt[5]{2}\\ +\left(293134+233308\sqrt[3]{2}+184248\sqrt[3]{4}\right)\sqrt[5]{4}\\ +\left(243521+194042\sqrt[3]{2}+153684\sqrt[3]{4}\right)\sqrt[5]{8}\\ +\left(211570+168880\sqrt[3]{2}+134352\sqrt[3]{4}\right)\sqrt[5]{16}\end{array}\right)}{125971200}}\delta \tau ^6+O(\delta \tau ^8)10.04573\delta \tau ^6+O(\delta \tau ^8).`$
Approximate values for $`\kappa _{n,1}`$ and $`\rho _{n,1}`$ are listed in Table 1. We note in passing that the magnitude of these leading non-trivial coefficients do not grow rapidly with increasing $`n`$.

### 4.4 Time Evolution Operators

The parameterisation given in equation (15) facilitates the computation of the time evolution operator $`U_n(\tau ,\delta \tau )U_n(\delta \tau )^{\tau /\delta \tau }`$ for trajectories of length $`\tau `$ where, in general, $`\tau \delta \tau `$. We assume as an induction hypothesis that

$$U_n(k\delta \tau ,\delta \tau )=U_n(\delta \tau )^k=\left(\begin{array}{cc}\mathrm{cos}[\kappa _n(\delta \tau )k\delta \tau ]& \frac{\mathrm{sin}\left[\kappa _n\left(\delta \tau \right)k\delta \tau \right]}{\rho _n\left(\delta \tau \right)}\\ \rho _n(\delta \tau )\mathrm{sin}[\kappa _n(\delta \tau )k\delta \tau ]& \mathrm{cos}[\kappa _n(\delta \tau )k\delta \tau ]\end{array}\right);$$

from this it immediately follows that $`U_n(\delta \tau )^jU_n(\delta \tau )^k=U_n(\delta \tau )^{j+k}`$ using simple trigonometric identities. Expressing the result in terms of the trajectory length rather than the number of integration steps, $`\tau k\delta \tau `$, we obtain

$$U_n(\tau ,\delta \tau )=\left(\begin{array}{cc}\mathrm{cos}[\kappa _n(\delta \tau )\tau ]& \frac{\mathrm{sin}\left[\kappa _n\left(\delta \tau \right)\tau \right]}{\rho _n\left(\delta \tau \right)}\\ \rho _n(\delta \tau )\mathrm{sin}[\kappa _n(\delta \tau )\tau ]& \mathrm{cos}[\kappa _n(\delta \tau )\tau ]\end{array}\right).$$
(18)
This time evolution matrix may be expanded about the exact Hamiltonian evolution as a Taylor series in $`\delta \tau `$,

$$U_n(\tau ,\delta \tau )=e^{\widehat{H}\tau }\left[\mathrm{𝕀}+U_{n,1}(\tau )\delta \tau ^{2n+2}+O(\delta \tau ^{2n+4})\right],$$
(19)
where from equation (11)

$$e^{\widehat{H}\tau }=\left(\begin{array}{cc}\mathrm{cos}\tau & \mathrm{sin}\tau \\ \mathrm{sin}\tau & \mathrm{cos}\tau \end{array}\right),$$
(20)
and

$$U_{n,1}(\tau )=\rho _{n,1}\mathrm{sin}\tau \left(\begin{array}{cc}\mathrm{sin}\tau & \mathrm{cos}\tau \\ \mathrm{cos}\tau & \mathrm{sin}\tau \end{array}\right)+\kappa _{n,1}\tau \left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$

It is instructive to compare equations (13) and (19).

## 5 Acceptance Rates

In this section we describe the calculation of the acceptance rate for MDMC (and thus for GHMC) for a system of $`V`$ uncoupled harmonic oscillators . The method of calculation is the same for both Langevin Monte Carlo and Hybrid Monte Carlo, and is independent of whether one uses lowest order leapfrog or a higher order discretisation scheme; the various algorithms differ only in the explicit form of the time evolution matrix $`U_n(\tau ,\delta \tau )`$ given in equation (18).

The Hamiltonian of equation (8) is a quadratic form

$$H\frac{1}{2}\underset{p=1}{\overset{V}{}}\left(\begin{array}{c}\omega _p\varphi _p\\ \pi _p\end{array}\right)^T\left(\begin{array}{c}\omega _p\varphi _p\\ \pi _p\end{array}\right),$$

so the change in “fictitious energy” over a trajectory is

$$\delta H_nH(\varphi (\tau ),\pi (\tau ))H(\varphi ,\pi )=\frac{1}{2}\underset{p=1}{\overset{V}{}}\left(\begin{array}{c}\omega _p\varphi _p\\ \pi _p\end{array}\right)^T\delta _n\left(\begin{array}{c}\omega _p\varphi _p\\ \pi _p\end{array}\right),$$

where $`\delta _nU_n^TU_n\mathrm{𝕀}`$, and we have abbreviated $`\varphi _p(0)`$ and $`\pi _p(0)`$ to $`\varphi _p`$ and $`\pi _p`$. Inserting the Taylor expansion of equation (19) and noting that $`(e^{\widehat{H}\tau })^Te^{\widehat{H}\tau }=1`$, corresponding to exact energy conservation for $`\delta \tau =0`$, we find that

$`\delta _n`$ $``$ $`\delta _{n,1}\delta \tau ^{2n+2}+\delta _{n,2}\delta \tau ^{2n+4}+O(\delta \tau ^{2n+6})`$
$`=`$ $`\left[\mathrm{𝕀}+U_{n,1}\delta \tau ^{2n+2}+O(\delta \tau ^{2n+4})\right]^T\left[\mathrm{𝕀}+U_{n,1}\delta \tau ^{2n+2}+O(\delta \tau ^{2n+4})\right]\mathrm{𝕀}`$
$`=`$ $`\left(U_{n,1}^T+U_{n,1}\right)\delta \tau ^{2n+2}+O(\delta \tau ^{2n+4})`$
The probability of the change in energy $`\delta H_n`$ having the value $`\xi `$ when averaged over the equilibrium distribution of starting points on phase space is

$$P_{\delta H_n}(\xi )=\frac{1}{Z}[d\varphi ][d\pi ]e^H\delta (\xi \delta H_n),$$

where as usual the “partition function” $`Z`$ is just the normalisation constant required to ensure that $`𝑑\xi P_{\delta H_n}(\xi )=1`$. It is convenient to choose an integral representation for the $`\delta `$-function as a superposition of plane waves, as then

$`P_{\delta H_n}(\xi )`$ $`=`$ $`{\displaystyle \frac{1}{Z}}{\displaystyle [d\varphi ][d\pi ]e^H_{\mathrm{}+i0}^{\mathrm{}+i0}\frac{d\eta }{2\pi }e^{i\eta (\xi \delta H_n)}}={\displaystyle \frac{d\eta }{2\pi }e^{i\eta \xi }\frac{1}{Z}[d\varphi ][d\pi ]e^{Hi\eta \delta H_n}}`$ (25)
$`=`$ $`{\displaystyle \frac{d\eta }{2\pi }e^{i\eta \xi }\frac{1}{Z^{}}[d\varphi ][d\pi ]\mathrm{exp}\left[\frac{1}{2}\underset{p,q=1}{\overset{V}{}}\left(\begin{array}{c}\omega _p\varphi _p\\ \pi _p\end{array}\right)^T\left\{\mathrm{𝕀}+i\eta \delta _n\right\}_{p,q}\left(\begin{array}{c}\omega _q\varphi _q\\ \pi _q\end{array}\right)\right]}`$
$`=`$ $`{\displaystyle \frac{d\eta }{2\pi }\frac{e^{i\eta \xi }}{\sqrt{det(\mathrm{𝕀}+i\eta \delta _n)}}}={\displaystyle \frac{d\eta }{2\pi }e^{i\eta \xi }e^{{\scriptscriptstyle \frac{1}{2}}tr\mathrm{ln}(\mathrm{𝕀}+i\eta \delta _n)}}.`$ (26)
It will prove useful to observe that the exact area-preservation property of the time evolution operator ensures that $`detU_n=1`$, and thus that

$$\mathrm{exp}tr\mathrm{ln}(\mathrm{𝕀}+\delta _n)=detU_n^TU_n=detU_n^TdetU_n=1,$$

and hence $`tr\mathrm{ln}(\mathrm{𝕀}+\delta _n)=0`$. This implies that the quantity $`tr\mathrm{ln}(\mathrm{𝕀}+i\eta \delta _n)`$ must vanish not only for $`i\eta =0`$ but also for $`i\eta =1`$. Expanding the logarithm, and making use of this fact, we find that

$$tr\mathrm{ln}(\mathrm{𝕀}+i\eta \delta _n)=i\eta (1i\eta )\frac{1}{2}tr\left(\delta _{n,1}^2\right)\delta \tau ^{4n+4}+O(\delta \tau ^{4n+6}).$$
(27)
We may perform an asymptotic expansion of the integral over $`\eta `$ using Laplace’s method. First we recast equation (27) into a form where the $`V`$ dependence is explicit

$$tr\mathrm{ln}(\mathrm{𝕀}+i\eta \delta _n)=i\eta (1i\eta )\frac{1}{2}𝜎\left(\delta _{n,1}^2\right)x+O\left(\frac{1,\eta ^3}{\sqrt[2n+2]{V}}\right),$$
(28)
where we have introduced the variable $`xV\delta \tau ^{4n+4}`$, and the spectral average $`𝜎(f)=\frac{1}{V}trf`$, which has a finite limit as $`V\mathrm{}`$. In equation (28) the correction terms are negligible if the only important contributions come from the regions where $`\eta ^3\sqrt[2n+2]{V}`$. Laplace’s method shows that such contributions are exponentially suppressed. Using equation (28) in equation (26) we obtain

$$P_{\delta H_n}(\xi )\frac{d\eta }{2\pi }e^{i\eta \xi }e^{\frac{1}{4}i\eta (1i\eta )𝜎\left(\delta _{n,1}^2\right)x}.$$

Completing the square gives

$$P_{\delta H_n}(\xi )\frac{1}{\sqrt{\pi 𝜎\left(\delta _{n,1}^2\right)x}}\mathrm{exp}\left[\frac{1}{𝜎\left(\delta _{n,1}^2\right)x}\left\{\xi \frac{𝜎\left(\delta _{n,1}^2\right)x}{4}\right\}^2\right],$$

and as $`\delta H_n=𝑑\xi P_{\delta H_n}(\xi )\xi =\frac{1}{4}𝜎\left(\delta _{n,1}^2\right)x`$ this may also be written as

$$P_{\delta H}(\xi )\frac{1}{\sqrt{4\pi \delta H}}\mathrm{exp}\left[\frac{(\xi \delta H)^2}{4\delta H}\right],$$

where we have supressed the index $`n`$.

The average Metropolis acceptance rate is now easily found,

$`P_{\text{acc}}`$ $``$ $`\mathrm{min}(1,e^{\delta H})={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\xi P_{\delta H}(\xi )\mathrm{min}(1,e^\xi )`$ (29)
$`=`$ $`{\displaystyle _{\mathrm{}}^0}𝑑\xi P_{\delta H}(\xi )+{\displaystyle _0^{\mathrm{}}}𝑑\xi P_{\delta H}(\xi )e^\xi `$
$`=`$ $`{\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle _{\frac{1}{2}\sqrt{\delta H}}^{\mathrm{}}}𝑑\xi e^{\xi ^2}=erfc\left({\displaystyle \frac{1}{2}}\sqrt{\delta H}\right)=erfc\left(\sqrt{{\displaystyle \frac{1}{8}}\delta H^2}\right).`$
The explicit form for

$`\delta H_n`$ $`=`$ $`{\displaystyle \frac{1}{4}}𝜎\left(\delta _{n,1}^2\right)x=2x\rho _{n,1}^2𝜎\left((\mathrm{sin}\omega \tau )^2\omega ^{4n+4}\right)`$ (30)
$`=`$ $`2x\rho _{n,1}^2\times {\displaystyle \frac{1}{V}}{\displaystyle \underset{p=1}{\overset{V}{}}}(\mathrm{sin}\omega _p\tau )^2\omega _p^{4n+4},`$
where

$$2\rho _{0,1}^2=\frac{1}{32}=0.03125$$

for the lowest-order leapfrog integration scheme, and

$`2\rho _{1,1}^2`$ $`=`$ $`{\displaystyle \frac{\left(40+32\sqrt[3]{2}+25\sqrt[3]{4}\right)}{41472}}0.002894,`$
$`2\rho _{2,1}^2`$ $`=`$ $`\begin{array}{c}\frac{\left(\begin{array}{c}z\left(221124374726+175505620552\sqrt[3]{2}+139299694184\sqrt[3]{4}\right)\\ +\left(192732972517+152971330676\sqrt[3]{2}+121414378108\sqrt[3]{4}\right)\sqrt[5]{2}\\ +\left(167855113064+133226114824\sqrt[3]{2}+105742133216\sqrt[3]{4}\right)\sqrt[5]{4}\\ +\left(145990198426+115872242816\sqrt[3]{2}+91967883784\sqrt[3]{4}\right)\sqrt[5]{8}\\ +\left(126961318646+100768983952\sqrt[3]{2}+79980474344\sqrt[3]{4}\right)\sqrt[5]{16}\end{array}\right)}{793437161472000}\hfill \end{array}0.004183`$
for higher-order “wiggles.”

### 5.1 One Dimensional Free Field Theory

For the case of free field theory we can compute the spectral averages using the explicit form for the frequency spectrum. Using the methods of appendix B we find that the spectral average

$`\overline{\sigma }_{m^2}^{(2)}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}(\mathrm{sin}\omega _p\tau )^2\omega _p^4`$ (33)
$`=`$ $`\left[33\mathrm{J}_0(4\tau )+4\mathrm{J}_2(4\tau )\mathrm{J}_4(4\tau )\right]`$
$`+`$ $`\left[22\mathrm{J}_0(4\tau )+4t\mathrm{J}_1(4\tau )+\mathrm{J}_2(4\tau )\right]m^2`$
$`+`$ $`{\displaystyle \frac{1}{2}}\left[1+(2t^21)\mathrm{J}_0(4\tau )+3t\mathrm{J}_1(4\tau )\right]m^4+O(m^6)`$
up to terms which vanish as $`V\mathrm{}`$. The finite volume corrections are given explicitly in appendix B. For the higher order integration schemes the corresponding spectral averages are also given in Appendix B.

The values for the logarithm of $`\delta H/x`$ are shown in Figure 1. The limiting values for $`\tau \mathrm{}`$ for the massless case ($`m=0`$) are given in Table 1; for $`m>0`$ the corresponding quantities diverge.

## 6 Autocorrelation Functions

### 6.1 Simple Markov Processes

Let $`(\varphi _0,\varphi _2,\mathrm{},\varphi _{N1})`$ be a sequence of field configurations generated by an equilibrated ergodic Markov process, and let $`\mathrm{\Omega }(\varphi )`$ denote the expectation value of some operator $`\mathrm{\Omega }`$ for $`\varphi `$ distributed according to the fixed point distribution of this Markov process. For simplicity we shall assume that $`\mathrm{\Omega }=0`$ in subsections 6.1 and 6.2. We may define an unbiased estimator $`\overline{\mathrm{\Omega }}`$ over the finite sequence of configurations by

$$\overline{\mathrm{\Omega }}\frac{1}{N}\underset{j=1}{\overset{N}{}}\mathrm{\Omega }(\varphi _j),$$

so $`\overline{\mathrm{\Omega }}=\frac{1}{N}_{j=1}^N\mathrm{\Omega }(\varphi _j)=\mathrm{\Omega }=0`$. The variance of this estimator is

$`\left(\overline{\mathrm{\Omega }}\overline{\mathrm{\Omega }}\right)^2`$ $`=`$ $`\overline{\mathrm{\Omega }}^2={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{j=0}{\overset{N1}{}}}{\displaystyle \underset{j^{}=0}{\overset{N1}{}}}\mathrm{\Omega }(\varphi _j^{})\mathrm{\Omega }(\varphi _j)`$
$`=`$ $`{\displaystyle \frac{1}{N^2}}\left\{{\displaystyle \underset{j=0}{\overset{N1}{}}}\mathrm{\Omega }(\varphi _j)^2+2{\displaystyle \underset{j=0}{\overset{N2}{}}}{\displaystyle \underset{j^{}=j+1}{\overset{N1}{}}}\mathrm{\Omega }(\varphi _j^{})\mathrm{\Omega }(\varphi _j)\right\}`$
$`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{\Omega }(\varphi )^2\left\{1+{\displaystyle \frac{2}{N}}{\displaystyle \underset{\mathrm{}=1}{\overset{N1}{}}}(N\mathrm{})C_\mathrm{\Omega }(\mathrm{})\right\},`$
where

$$C_\mathrm{\Omega }(\mathrm{})\frac{\mathrm{\Omega }(\varphi _{\mathrm{}})\mathrm{\Omega }(\varphi _0)}{\mathrm{\Omega }(\varphi )^2}$$
(34)
is the autocorrelation function for $`\mathrm{\Omega }`$. If the Markov process is ergodic, then for large $`\mathrm{}`$,

$$\left|C_\mathrm{\Omega }(\mathrm{})\right|\lambda _{\mathrm{max}}^{\mathrm{}}e^{\mathrm{}/N_{\text{exp}}},$$
(35)
where $`\lambda _{\mathrm{max}}`$ is the second-largest eigenvalue of the Markov matrix and $`N_{\text{exp}}`$ is the exponential autocorrelation time. If $`NN_{\text{exp}}`$ then

$`\overline{\mathrm{\Omega }}^2`$ $`=`$ $`\left\{1+2{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}C_\mathrm{\Omega }(\mathrm{})\right\}{\displaystyle \frac{\mathrm{\Omega }(\varphi )^2}{N}}\left[1+O\left({\displaystyle \frac{N_{\text{exp}}}{N}}\right)\right]`$ (36)
$``$ $`\{1+2A_\mathrm{\Omega }\}{\displaystyle \frac{\mathrm{\Omega }(\varphi )^2}{N}}\left[1+O\left({\displaystyle \frac{N_{\text{exp}}}{N}}\right)\right],`$
where $`A_\mathrm{\Omega }_{\mathrm{}=1}^{\mathrm{}}C_\mathrm{\Omega }(\mathrm{})`$ is the integrated autocorrelation function for the operator $`\mathrm{\Omega }`$.

This result tells us that on average $`1+2A_\mathrm{\Omega }`$ correlated measurements are needed to reduce the variance by the same amount as a single truly independent measurement.

### 6.2 Hybrid Stochastic Processes

Suppose that a sequence of configurations $`(\varphi (t_0),\varphi (t_1),\varphi (t_2),\mathrm{},\varphi (t_N))`$ is generated from $`\varphi (t_0)`$ as follows: the configuration $`\varphi (t_j)`$ is generated from $`\varphi (t_{j1})`$ by choosing momenta $`\pi (t_{j1})`$ as described in Section (2.2), and integrating Hamilton’s equations for a time interval $`\tau _jt_jt_{j1}`$, where each trajectory length $`\tau _j`$ is chosen randomly from the distribution $`P_R(\tau _j)`$. The autocorrelation function $`C_\mathrm{\Omega }`$ defined by equation (34) may be expressed in terms of the autocorrelation function

$$C_\mathrm{\Omega }(\tau _1,\mathrm{},\tau _{\mathrm{}})\frac{\mathrm{\Omega }\left(\varphi (t_{\mathrm{}})\right)\mathrm{\Omega }\left(\varphi (t_0)\right)}{\mathrm{\Omega }^2}$$

by averaging it over the refresh distribution

$$C_\mathrm{\Omega }(\mathrm{})=_0^{\mathrm{}}𝑑\tau _1\mathrm{}𝑑\tau _{\mathrm{}}P_R(\tau _1)\mathrm{}P_R(\tau _{\mathrm{}})C_\mathrm{\Omega }(\tau _1,\mathrm{},\tau _{\mathrm{}}).$$
(37)
The integrated autocorrelation function then becomes

$$A_\mathrm{\Omega }=\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}𝑑\tau _1\mathrm{}𝑑\tau _{\mathrm{}}P_R(\tau _1)\mathrm{}P_R(\tau _{\mathrm{}})C_\mathrm{\Omega }(\tau _1,\mathrm{},\tau _{\mathrm{}}).$$

If we wish to determine autocorrelations in terms of the total elapsed fictitious time $`t=t_{\mathrm{}}t_0=_{k=1}^{\mathrm{}}\tau _k`$ of the sequence of trajectories we may introduce yet another autocorrelation function by

$$C_\mathrm{\Omega }(t)\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}𝑑\tau _1\mathrm{}𝑑\tau _{\mathrm{}}P_R(\tau _1)\mathrm{}P_R(\tau _{\mathrm{}})\delta \left(t\underset{j=0}{\overset{\mathrm{}}{}}\tau _j\right)C_\mathrm{\Omega }(\tau _1,\mathrm{},\tau _{\mathrm{}}),$$
(38)
and in terms of this we find that

$$A_\mathrm{\Omega }=_0^{\mathrm{}}𝑑tC_\mathrm{\Omega }(t).$$
(39)
The function $`C_\mathrm{\Omega }(t)`$ has the advantage of giving the autocorrelations as a function of MD time which is approximately proportional to computer time.

If we make the reasonable assumption that the cost of the computation is proportional to the total fictitious (MD) time for which we have to integrate Hamilton’s equations, and to the volume $`V`$ of the lattice222For free field theory the volume $`V`$ is just the number of uncoupled harmonic oscillators., then the cost $``$ per independent configuration is proportional to $`(1+2A_\mathrm{\Omega })V\overline{\tau }/\delta \tau `$ with $`\overline{\tau }`$ denoting the average length of a trajectory. The meaning of “independent configuration” was discussed in section 6.1, and depends on the particular operator $`\mathrm{\Omega }`$ under consideration. The optimal trajectory length is obtained by minimising the cost, that is by choosing $`\overline{\tau }`$ so as to satisfy

$`{\displaystyle \frac{d}{d\overline{\tau }}}=0`$ $``$ $`1+2A_\mathrm{\Omega }+2\overline{\tau }{\displaystyle \frac{dA_\mathrm{\Omega }}{d\overline{\tau }}}=0`$
$`{\displaystyle \frac{d}{d\vartheta }}=0`$ $``$ $`{\displaystyle \frac{dA_\mathrm{\Omega }}{d\vartheta }}=0`$ (40)
### 6.3 Autocorrelation Functions for Polynomial Operators

In order to carry out these calculations we make a simplifying assumption:

###### Assumption 6.3.1

The acceptance probability $`\mathrm{min}(1,e^{\delta H})`$ for each trajectory may be replaced by its value averaged over phase space $`P_{\text{acc}}\mathrm{min}(1,e^{\delta H})`$; we neglect correlations between successive trajectories. Including such correlations leads to seemingly intractable complications. It is not obvious that our assumption corresponds to any systematic approximation except, of course, that it is valid when $`\delta H=0`$.

The action of a generalised HMC update on the fields $`\varphi `$, their conjugate fictitious momenta $`\pi `$, and the Gaussian noise $`\xi `$ for the trajectory under consideration from Equation (4) is

$$\left(\begin{array}{c}\omega \varphi ^{}\\ \pi ^{}\\ \xi ^{}\end{array}\right)=𝔄(\tau ,\vartheta ,y,\delta H)\left(\begin{array}{c}\omega \varphi \\ \pi \\ \xi \end{array}\right),$$

where the matrix $`𝔄`$ depends on the trajectory length $`\tau `$, the noise rotation angle $`\vartheta `$ (which can be chosen independently for each trajectory if we so wished), the uniform random number $`y`$ used in the Metropolis accept/reject test, and the value of $`\delta H`$.

We may ignore the corrections of non-leading order in $`\delta \tau `$ to the MD evolution operator because for any given value of $`P_{\text{acc}}`$ of $`O(1)`$ there is a corresponding value of $`x`$ which is also $`O(1)`$, and thus $`\delta \tau `$ is of order $`V^{1/(4n+4)}`$. These corrections therefore only contribute to the autocorrelations through the acceptance rate itself at leading order in the large volume expansion.

From the leading order contribution (20) we obtain

$$𝔄=\left(\begin{array}{ccc}\mathrm{cos}\omega \tau & \mathrm{sin}\omega \tau & 0\\ \mathrm{sin}\omega \tau \mathrm{cos}\vartheta & \mathrm{cos}\omega \tau \mathrm{cos}\vartheta & \mathrm{sin}\vartheta \\ \mathrm{sin}\omega \tau \mathrm{sin}\vartheta & \mathrm{cos}\omega \tau \mathrm{sin}\vartheta & \mathrm{cos}\vartheta \end{array}\right)\theta _++\left(\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\vartheta & \mathrm{sin}\vartheta \\ 0& \mathrm{sin}\vartheta & \mathrm{cos}\vartheta \end{array}\right)\theta _{}$$

with $`\theta _\pm \theta \left[\pm \left(e^{\delta H}y\right)\right]`$, and $`\omega `$ being the appropriate frequency for each mode.

#### 6.3.1 Linear Operators

We are interested in calculating the autocorrelation function for a general linear operator

$$\mathrm{\Omega }=\underset{p}{}\mathrm{\Omega }_p\varphi _p$$
(41)
for a set of uncoupled harmonic oscillators $`\{\varphi _p\}`$. Such an operator is of course connected, meaning that $`\mathrm{\Omega }=_p\mathrm{\Omega }_p\varphi _p=0`$. Its autocorrelation function may be expressed in terms of the autocorrelation functions for the individual harmonic oscillators themselves, for

$$C_\mathrm{\Omega }=\frac{\mathrm{\Omega }(t)\mathrm{\Omega }(0)}{\mathrm{\Omega }^2}=\underset{p,p^{}}{}\mathrm{\Omega }_p\mathrm{\Omega }_p^{}\frac{\varphi _p(t)\varphi _p^{}(0)}{\mathrm{\Omega }^2},$$

and as the oscillators are uncoupled

$$\varphi _p(t)\varphi _p^{}(0)=\delta _{p,p^{}}\varphi _p(t)\varphi _p(0)+(1\delta _{p,p^{}})\varphi _p\varphi _p^{}=\delta _{p,p^{}}\varphi _p(t)\varphi _p(0),$$

thus

$$C_\mathrm{\Omega }(t)=\frac{_p\mathrm{\Omega }_p^2C_{\varphi _p}(t)\varphi _p^2}{_p\mathrm{\Omega }_p^2\varphi _p^2},$$

where

$$C_{\varphi _p}(t)=\frac{\varphi _p(t)\varphi _p(0)}{\varphi _p^2}.$$

Let us proceed to calculate these single mode autocorrelation functions; while doing so we can drop the subscript $`p`$ for notational simplicity.

The average over the Gaussian distribution of $`\xi `$ gives $`\xi =0`$, so we can drop the last column of $`𝔄`$, and since a new $`\xi `$ is taken from a heatbath at the start of each trajectory we may drop the last row also. This leaves us with the basis $`\mathrm{\Phi }^{(1)}(\omega \varphi ,\pi )`$, which is updated by $`\mathrm{\Phi }_{}^{(1)}{}_{}{}^{}=𝔄^{(1)}\mathrm{\Phi }^{(1)}`$ where the matrix

$$𝔄^{(1)}\left(\begin{array}{cc}\mathrm{cos}\omega \tau & \mathrm{sin}\omega \tau \\ \mathrm{sin}\omega \tau \mathrm{cos}\vartheta & \mathrm{cos}\omega \tau \mathrm{cos}\vartheta \end{array}\right)\theta _++\left(\begin{array}{cc}1& 0\\ 0& \mathrm{cos}\vartheta \end{array}\right)\theta _{}.$$
(42)
#### 6.3.2 Quadratic Operators

Let $`\mathrm{\Omega }`$ be a generic quadratic operator for the set of uncoupled harmonic oscillators $`\{\varphi _p\}`$,

$$\mathrm{\Omega }\underset{pq}{}\mathrm{\Omega }_{pq}\varphi _p\varphi _q$$
(43)
whose connected part is $`\mathrm{\Omega }_c=\mathrm{\Omega }\mathrm{\Omega }`$. The autocorrelation function for $`\mathrm{\Omega }_c`$ is333The integrated autocorrelation function for a disconnected operator diverges in general.

$$C_{\mathrm{\Omega }_c}(t)=\frac{\mathrm{\Omega }_c(t)\mathrm{\Omega }_c(0)}{\mathrm{\Omega }_c^2}$$

If we define the elementary autocorrelation functions for linear and quadratic modes to be $`C_{\varphi _p}(t)`$ and $`C_{(\varphi _p^2)_c}(t)`$, then we may express $`C_{\mathrm{\Omega }_c}(t)`$ in terms of them:

$$C_{\mathrm{\Omega }_c}(t)=\frac{1}{\mathrm{\Omega }_c^2}\left\{\underset{p}{}\mathrm{\Omega }_{pp}^2C_{(\varphi _p^2)_c}(t)\left((\varphi _p^2)_c\right)^2+2\underset{p>q}{}\mathrm{\Omega }_{pq}^2C_{\varphi _p}(t)C_{\varphi _q}(t)(\varphi _p^2)_c(\varphi _q^2)_c\right\}$$

Higher degree polynomial operators may be treated in the same way. As before we calculate the single mode autocorrelation functions and again drop the subscript $`p`$ while doing so.

We cannot directly average the GHMC update matrix $`𝔄`$ over $`\xi `$ as we did in the linear case, but we can linearise the problem by considering a homogeneous quadratic operator as a linear combination of the quadratic monomials $`((\omega \varphi )^2,\pi \omega \varphi ,\pi ^2,\xi \omega \varphi ,\xi \pi ,\xi ^2)`$. These quadratic monomials are updated by the symmetrised tensor product of the update matrix $`𝔄`$,

$$\left(\begin{array}{c}(\omega \varphi ^{})^2\\ \pi ^{}\omega \varphi ^{}\\ \pi ^2\\ \xi ^{}\omega \varphi ^{}\\ \xi ^{}\pi ^{}\\ \xi ^2\end{array}\right)=𝔄_s𝔄\left(\begin{array}{c}(\omega \varphi )^2\\ \pi \omega \varphi \\ \pi ^2\\ \xi \omega \varphi \\ \xi \pi \\ \xi ^2\end{array}\right),$$

where the update matrix $`𝔄_s𝔄`$ is explicitly

$$\left(\begin{array}{cccccc}𝔄_{\varphi \varphi }^2& 2𝔄_{\varphi \pi }𝔄_{\varphi \varphi }& 𝔄_{\varphi \pi }^2& 2𝔄_{\varphi \xi }𝔄_{\varphi \varphi }& 2𝔄_{\varphi \xi }𝔄_{\varphi \pi }& 𝔄_{\varphi \xi }^2\\ 𝔄_{\pi \varphi }𝔄_{\varphi \varphi }& 𝔄_{\pi \pi }𝔄_{\varphi \varphi }+𝔄_{\pi \varphi }𝔄_{\varphi \pi }& 𝔄_{\pi \pi }𝔄_{\varphi \pi }& 𝔄_{\pi \xi }𝔄_{\varphi \varphi }+𝔄_{\pi \varphi }𝔄_{\varphi \xi }& 𝔄_{\pi \xi }𝔄_{\varphi \pi }+𝔄_{\pi \pi }𝔄_{\varphi \xi }& 𝔄_{\pi \xi }𝔄_{\varphi \xi }\\ 𝔄_{\pi \varphi }^2& 2𝔄_{\pi \pi }𝔄_{\pi \varphi }& 𝔄_{\pi \pi }^2& 2𝔄_{\pi \xi }𝔄_{\pi \varphi }& 2𝔄_{\pi \xi }𝔄_{\pi \pi }& 𝔄_{\pi \xi }^2\\ & & & & & \\ 𝔄_{\xi \varphi }𝔄_{\varphi \varphi }& 𝔄_{\xi \pi }𝔄_{\varphi \varphi }+𝔄_{\xi \varphi }𝔄_{\varphi \pi }& 𝔄_{\xi \pi }𝔄_{\varphi \pi }& 𝔄_{\xi \xi }𝔄_{\varphi \varphi }+𝔄_{\xi \varphi }𝔄_{\varphi \xi }& 𝔄_{\xi \xi }𝔄_{\varphi \pi }+𝔄_{\xi \pi }𝔄_{\varphi \xi }& 𝔄_{\xi \xi }𝔄_{\varphi \xi }\\ 𝔄_{\xi \varphi }𝔄_{\pi \varphi }& 𝔄_{\xi \pi }𝔄_{\pi \varphi }+𝔄_{\xi \varphi }𝔄_{\pi \pi }& 𝔄_{\xi \pi }𝔄_{\pi \pi }& 𝔄_{\xi \xi }𝔄_{\pi \varphi }+𝔄_{\xi \varphi }𝔄_{\pi \xi }& 𝔄_{\xi \xi }𝔄_{\pi \pi }+𝔄_{\xi \pi }𝔄_{\pi \xi }& 𝔄_{\xi \xi }𝔄_{\pi \xi }\\ & & & & & \\ 𝔄_{\xi \varphi }^2& 2𝔄_{\xi \pi }𝔄_{\xi \varphi }& 𝔄_{\xi \pi }^2& 2𝔄_{\xi \xi }𝔄_{\xi \varphi }& 2𝔄_{\xi \xi }𝔄_{\xi \pi }& 𝔄_{\xi \xi }^2\end{array}\right).$$

As the update is now linear we can average over the Gaussian distribution of $`\xi `$ as before. The basis monomials become $`((\omega \varphi )^2,\pi \omega \varphi ,\pi ^2,0,0,\xi ^2=1)`$, and this leads to three simplifications:

* The fourth and fifth columns are multiplied by zero, and can thus be dropped.
* The fourth and fifth rows are not of interest and may be dropped too, as $`\xi `$ is chosen from a heathbath before the next trajectory, and we already know how the linear monomials $`(\varphi ,\pi )`$ are updated.
* The last row may be replaced by $`(0,0,0,0,0,1)`$ as we know that $`\xi ^{}`$ is Gaussian distributed and thus $`\xi ^2=1`$.

We are thus led to consider the update of the inhomogenous monomials of even degree in $`\omega \varphi `$ and $`\pi `$ alone, namely $`\mathrm{\Phi }^{(2)}((\omega \varphi )^2,\pi \omega \varphi ,\pi ^2,1)`$, which is updated by $`\mathrm{\Phi }_{}^{(2)}{}_{}{}^{}=𝔄^{(2)}\mathrm{\Phi }^{(2)}`$ where

$`𝔄^{(2)}`$ $``$ $`\left(\begin{array}{cccc}𝔄_{\varphi \varphi }^2& 2𝔄_{\varphi \pi }𝔄_{\varphi \varphi }& 𝔄_{\varphi \pi }^2& 𝔄_{\varphi \xi }^2\\ 𝔄_{\pi \varphi }𝔄_{\varphi \varphi }& 𝔄_{\pi \pi }𝔄_{\varphi \varphi }+𝔄_{\pi \varphi }𝔄_{\varphi \pi }& 𝔄_{\pi \pi }𝔄_{\varphi \pi }& 𝔄_{\pi \xi }𝔄_{\varphi \xi }\\ 𝔄_{\pi \varphi }^2& 2𝔄_{\pi \pi }𝔄_{\pi \varphi }& 𝔄_{\pi \pi }^2& 𝔄_{\pi \xi }^2\\ 0& 0& 0& 1\end{array}\right)`$ (48)
$`=`$ $`\left(\begin{array}{cccc}\mathrm{cos}^2\omega \tau & 2\mathrm{cos}\omega \tau \mathrm{sin}\omega \tau & \mathrm{sin}^2\omega \tau & 0\\ \mathrm{cos}\omega \tau \mathrm{sin}\omega \tau \mathrm{cos}\vartheta & (\mathrm{cos}^2\omega \tau \mathrm{sin}^2\omega \tau )\mathrm{cos}\vartheta & \mathrm{cos}\omega \tau \mathrm{sin}\omega \tau \mathrm{cos}\vartheta & 0\\ \mathrm{sin}^2\omega \tau \mathrm{cos}^2\vartheta & 2\mathrm{cos}\omega \tau \mathrm{sin}\omega \tau \mathrm{cos}^2\vartheta & \mathrm{cos}^2\omega \tau \mathrm{cos}^2\vartheta & \mathrm{sin}^2\vartheta \\ 0& 0& 0& 1\end{array}\right)\theta _+`$ (58)
$`+\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \mathrm{cos}\vartheta & 0& 0\\ 0& 0& \mathrm{cos}^2\vartheta & \mathrm{sin}^2\vartheta \\ 0& 0& 0& 1\end{array}\right)\theta _{}.`$
### 6.4 Laplace Transforms of Autocorrelation Functions

We can extract more information about the autocorrelation function $`C_\mathrm{\Omega }(t)`$ than just the integrated autocorrelation function $`A_\mathrm{\Omega }`$. We shall discuss this further in Section 7. In order to do this it is convenient to compute the Laplace transform of the autocorrelation function

$`F_\mathrm{\Omega }(\beta )`$ $``$ $`C_\mathrm{\Omega }(t){\displaystyle _0^{\mathrm{}}}𝑑tC_\mathrm{\Omega }(t)e^{\beta t}`$
$`=`$ $`{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}𝑑\tau _1\mathrm{}𝑑\tau _{\mathrm{}}P_R(\tau _1)\mathrm{}P_R(\tau _{\mathrm{}})C_\mathrm{\Omega }(\tau _1,\mathrm{},\tau _{\mathrm{}})e^{\beta \tau _1}\mathrm{}e^{\beta \tau _{\mathrm{}}}.`$
We may generalise the results of sections 6.3.1 and 6.3.2 and observe that the update step for the vector of inhomogeneous monomials in $`\omega \varphi `$ and $`\pi `$, $`\mathrm{\Phi }^{(d)}`$, is of the form444After averaging over the distributions of $`\vartheta `$ and $`y`$ which are chosen independently for each trajectory.

$$\mathrm{\Phi }^{(d)}(\tau )=𝔄^{(d)}(\tau )\mathrm{\Phi }^{(d)}(0)$$

The matrix $`𝔄^{(d)}`$ is block upper triangular with with first block corresponding to the $`d+1`$ homogeneous monomials of degree $`d`$, the next block to the $`d1`$ homogeneous monomials of degree $`d2`$, and so forth. The dimension of the matrix, which is the number of even or odd inhomogeneous monomials of degree $`d`$ or less, is $`{\scriptscriptstyle \frac{1}{4}}(d+2)^2`$ for $`d`$ even and $`{\scriptscriptstyle \frac{1}{4}}(d+3)(d+1)`$ for $`d`$ odd.

The connected operator $`(\varphi ^d)_c\varphi ^d\varphi ^d=v^T\mathrm{\Phi }^{(d)}`$ where $`v=(1/\omega ^d,0,\mathrm{},0,\varphi ^d)`$, and the $`\mathrm{}`$-trajectory autocorrelation function is thus

$$C_{(\varphi ^d)_c}(\tau _1,\mathrm{},\tau _{\mathrm{}})=\frac{v^T[𝔄^{(d)}(\tau _1))\mathrm{}𝔄^{(d)}(\tau _{\mathrm{}})\mathrm{\Phi }^{(d)}(0)]v^T\mathrm{\Phi }^{(d)}(0)}{\left(v^T\mathrm{\Phi }^{(d)}\right)^2}.$$

By virtue of assumption 6.3.1 we may replace each matrix $`𝔄^{(d)}`$ by its phase space average,

$$C_{(\varphi ^d)_c}(\tau _1,\mathrm{},\tau _{\mathrm{}})=\frac{v^T_{j=1}^{\mathrm{}}𝔄^{(d)}(\tau _j)\mathrm{\Phi }^{(d)}\mathrm{\Phi }_{}^{(d)}{}_{}{}^{T}v}{v^T\mathrm{\Phi }^{(d)}\mathrm{\Phi }_{}^{(d)}{}_{}{}^{T}v},$$

and the Laplace transform of the degree $`d`$ single mode autocorrelation function is therefore

$$F_{(\varphi ^d)_c}(\beta )=C_{(\varphi ^d)_c}(t)=\frac{v^T\left[_{\mathrm{}=1}^{\mathrm{}}\overline{𝔄}_{}^{(d)}{}_{}{}^{\mathrm{}}\right]\mathrm{\Phi }^{(d)}\mathrm{\Phi }_{}^{(d)}{}_{}{}^{T}v}{v^T\mathrm{\Phi }^{(d)}\mathrm{\Phi }_{}^{(d)}{}_{}{}^{T}v},$$

where

$$\overline{𝔄}^{(d)}_0^{\mathrm{}}𝑑\tau P_R(\tau )e^{\beta \tau }𝔄^{(d)}(\tau ).$$
(60)
Summing the geometric series we obtain a simple expression for the Laplace transformed autocorrelation function,

$$1+F_{(\varphi ^d)_c}(\beta )=\frac{v^T\left[\mathrm{𝕀}\overline{𝔄}^{(d)}\right]^1\mathrm{\Phi }^{(d)}\mathrm{\Phi }_{}^{(d)}{}_{}{}^{T}v}{v^T\mathrm{\Phi }^{(d)}\mathrm{\Phi }_{}^{(d)}{}_{}{}^{T}v},$$
(61)
In order to evaluate the expectation values, we need the Gaussian averages for the equilibrium probability distribution of equation (1) for the Hamiltonian given in equation (8):

$$\begin{array}{ccc}\hfill \varphi _p^{2k+1}=& \pi _p^{2k+1}& =0\hfill \\ \hfill (\omega _p\varphi _p)^{2k}=& \pi _p^{2k}& =\frac{1}{\sqrt{\pi }}2^k\mathrm{\Gamma }(k+{\scriptscriptstyle \frac{1}{2}})\hfill \end{array}\}$$

and

$$\varphi _p^\alpha \varphi _q^\beta \pi _r^\gamma \pi _s^\delta =\varphi _p^\alpha \varphi _q^\beta \pi _r^\gamma \pi _s^\delta $$

when $`pq`$ and $`rs`$.

#### 6.4.1 Linear Operators

For linear operators the matrix of expectation values in equation (61) is

$$\mathrm{\Phi }^{(1)}\mathrm{\Phi }_{}^{(1)}{}_{}{}^{T}=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$

The explicit form of the Laplace transformed update matrix is obtained from equation (42) by first averaging over phase space, which replaces $`\theta _+`$ by $`P_{\text{acc}}(\omega \tau )`$ and $`\theta _{}`$ by $`1P_{\text{acc}}(\omega \tau )`$, and then equation (60) gives

$$\overline{𝔄}^{(1)}=\left(\begin{array}{cc}G_{1,0}^++G_{0,0}^{}& G_{0,1}^+\\ G_{0,1}^+\mathrm{cos}\vartheta & \left(G_{1,0}^+G_{0,0}^{}\right)\mathrm{cos}\vartheta \end{array}\right)$$
(62)
where

$$G_{j,k}^+_0^{\mathrm{}}𝑑\tau P_R(\tau )e^{\beta \tau }P_{\text{acc}}(\tau )(\mathrm{cos}\omega \tau )^j(\mathrm{sin}\omega \tau )^k$$

and

$$G_{j,k}^{}_0^{\mathrm{}}𝑑\tau P_R(\tau )e^{\beta \tau }\left(1P_{\text{acc}}(\tau )\right)(\mathrm{cos}\omega \tau )^j(\mathrm{sin}\omega \tau )^k.$$

The Laplace transform of the linear single mode autocorrelation function is

$$F_\varphi (\beta )=\frac{G_{1,0}^++G_{0,0}^{}\left(G_{1,0}^{+}{}_{}{}^{2}G_{0,0}^{}{}_{}{}^{2}+G_{0,1}^{+}{}_{}{}^{2}\right)\mathrm{cos}\vartheta }{1G_{1,0}^+G_{0,0}^{}+\left(G_{1,0}^++G_{0,0}^{}+G_{1,0}^{+}{}_{}{}^{2}G_{0,0}^{}{}_{}{}^{2}+G_{0,1}^{+}{}_{}{}^{2}\right)\mathrm{cos}\vartheta },$$
(63)
which for the special case of $`\vartheta =\pi /2`$ (HMC) simplifies to

$$F_\varphi (\beta ;\vartheta =\pi /2)=\frac{G_{1,0}^++G_{0,0}^{}}{1G_{1,0}^+G_{0,0}^{}}.$$

#### 6.4.2 Quadratic Operators

For quadratic operators the matrix of expectation values is

$$\mathrm{\Phi }^{(2)}\mathrm{\Phi }_{}^{(2)}{}_{}{}^{T}=\left(\begin{array}{cccc}3& 0& 1& 1\\ 0& 1& 0& 0\\ 1& 0& 3& 1\\ 1& 0& 1& 1\end{array}\right);$$

and from equations (58) and (60)

$$\overline{𝔄}^{(2)}=\left(\begin{array}{cccc}G_{2,0}^++G_{0,0}^{}& 2G_{1,1}^+& G_{0,2}^+& 0\\ G_{1,1}^+\mathrm{cos}\vartheta & \left(G_{2,0}^+G_{0,2}^+G_{0,0}^{}\right)\mathrm{cos}\vartheta & G_{1,1}^+\mathrm{cos}\vartheta & 0\\ G_{0,2}^+\mathrm{cos}^2\vartheta & 2G_{1,1}^+\mathrm{cos}^2\vartheta & \left(G_{2,0}^++G_{0,0}^{}\right)\mathrm{cos}^2\vartheta & G_{0,0}\mathrm{sin}^2\vartheta \\ 0& 0& 0& G_{0,0}\end{array}\right),$$

where

$$G_{j,k}G_{j,k}^++G_{j,k}^{}=_0^{\mathrm{}}𝑑\tau P_R(\tau )e^{\beta \tau }(\mathrm{cos}\omega \tau )^j(\mathrm{sin}\omega \tau )^k.$$
(64)
The Laplace transform of the quadratic single mode autocorrelation function $`F_{(\varphi ^2)_c}(\beta )`$ is

$$\frac{\left[\begin{array}{c}G_{2,0}^++G_{0,0}^{}+\left(G_{2,0}^{+}{}_{}{}^{2}2G_{1,1}^{+}{}_{}{}^{2}+G_{0,2}^+G_{0,0}^{}+G_{0,2}^+G_{2,0}^++G_{0,0}^{}{}_{}{}^{2}\right)\mathrm{cos}\vartheta \\ \left(G_{2,0}^+G_{0,2}^++G_{0,0}^{}\right)\left(G_{0,0}^{}+G_{2,0}^++G_{0,2}^+\right)(\mathrm{cos}\vartheta )^2\\ +\left(G_{0,0}^{}+G_{2,0}^++G_{0,2}^+\right)\left(G_{2,0}^{+}{}_{}{}^{2}2G_{0,2}^+G_{2,0}^++G_{0,2}^{+}{}_{}{}^{2}+4G_{1,1}^{+}{}_{}{}^{2}G_{0,0}^{}{}_{}{}^{2}\right)(\mathrm{cos}\vartheta )^3\end{array}\right]}{\left[\begin{array}{c}1G_{0,0}^{}G_{2,0}^+\\ \left(G_{0,2}^+G_{0,0}^{}2G_{1,1}^{+}{}_{}{}^{2}+G_{0,0}^{}{}_{}{}^{2}+G_{0,2}^+G_{2,0}^+G_{0,0}^{}+G_{2,0}^+G_{2,0}^{+}{}_{}{}^{2}G_{0,2}^+\right)\mathrm{cos}\vartheta \\ \left(G_{2,0}^{+}{}_{}{}^{2}+G_{0,0}^{}2G_{2,0}^+G_{0,0}^{}+G_{0,2}^{+}{}_{}{}^{2}G_{0,0}^{}{}_{}{}^{2}+G_{2,0}^+\right)(\mathrm{cos}\vartheta )^2\\ +\left(\begin{array}{c}2G_{1,1}^{+}{}_{}{}^{2}G_{0,2}^{+}{}_{}{}^{3}4G_{1,1}^{+}{}_{}{}^{2}G_{0,0}^{}+G_{0,2}^+G_{2,0}^{+}{}_{}{}^{2}G_{0,2}^+G_{2,0}^+\\ +G_{0,2}^{+}{}_{}{}^{2}G_{2,0}^+G_{0,2}^+G_{0,0}^{}+G_{2,0}^+G_{0,0}^{}{}_{}{}^{2}G_{2,0}^{+}{}_{}{}^{3}+G_{2,0}^{+}{}_{}{}^{2}\\ G_{0,0}^{}{}_{}{}^{2}4G_{2,0}^+G_{1,1}^{+}{}_{}{}^{2}G_{2,0}^{+}{}_{}{}^{2}G_{0,0}^{}4G_{1,1}^{+}{}_{}{}^{2}G_{0,2}^+\\ +G_{0,2}^+G_{0,0}^{}{}_{}{}^{2}G_{0,2}^{+}{}_{}{}^{2}G_{0,0}^{}+G_{0,0}^{}{}_{}{}^{3}+2G_{0,2}^+G_{2,0}^+G_{0,0}^{}\end{array}\right)(\mathrm{cos}\vartheta )^3\end{array}\right]},$$
(65)
which for the special case of $`\vartheta =\pi /2`$ (HMC) simplifies to

$$F_{(\varphi ^2)_c}(\beta ;\vartheta =\pi /2)=\frac{G_{2,0}^++G_{0,0}^{}}{1G_{2,0}^+G_{0,0}^{}}.$$

### 6.5 Exponentially Distributed Trajectory Lengths

To proceed further we need to choose a specific form for the momentum refresh distribution. In this section we will analyse the case of exponentially distributed trajectory lengths, $`P_R(\tau )=re^{r\tau }`$ where the parameter is just the inverse average trajectory length $`r=1/\overline{\tau }`$. In section 6.6 we will consider fixed length trajectories.

We make the approximation that $`\delta H`$ and thus $`P_{\text{acc}}`$ are independent of $`\tau `$ (c.f., Figure 1), so

$$G_{j,k}^+\overline{P}_{\text{acc}}G_{j,k},G_{j,k}^{}(1\overline{P}_{\text{acc}})G_{j,k};$$

where $`\overline{P}_{\text{acc}}P_{\text{acc}}(\overline{\tau })`$. This approximation is only made in order to avoid unpleasant integrals which cannot be evaluated in closed form. The integral in equation (64) may be evaluated, and we find

$$G_{j,k}^{\text{exp}}=\underset{\mu =0}{\overset{j}{}}\underset{\nu =0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{j}{\mu }\right)\left(\genfrac{}{}{0pt}{}{k}{\nu }\right)\frac{({\scriptscriptstyle \frac{1}{2}})^{j+k}(1)^{\nu +{\scriptscriptstyle \frac{1}{2}}k}B}{(\beta \overline{\tau }+1)^2+\phi ^2(j+k2\mu 2\nu )^2},$$
(66)
where we have introduced the dimensionless quantity $`\phi \omega \overline{\tau }`$ and

$$B=\{\begin{array}{cc}\beta \overline{\tau }+1& \text{if }{\scriptscriptstyle \frac{1}{2}}k\mathrm{},\hfill \\ (j+k2\mu 2\nu )\phi & \text{if }{\scriptscriptstyle \frac{1}{2}}k\mathrm{}\text{.}\hfill \end{array}$$

#### 6.5.1 Linear Operators

The Laplace transform of the linear single mode autocorrelation function for exponentially distributed trajectories is obtained by substituting the explicit form for $`G^\pm `$ into equation (63), and we obtain

$$F_\varphi ^{\text{exp}}(\beta )=\frac{\left[\begin{array}{c}r(\beta ^3+((12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +3)r\beta ^2\\ +\left(2(12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +(1\overline{P}_{\text{acc}})\phi ^2+3\right)r^2\beta \\ +(((1\overline{P}_{\text{acc}})^2\phi ^2+12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +(1\overline{P}_{\text{acc}})\phi ^2+1)r^3)\end{array}\right]}{\left[\begin{array}{c}\beta ^4+\left((12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +3\right)r\beta ^3+\left(2(12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +3+\phi ^2\right)r^2\beta ^2\\ +\left(\left((1\overline{P}_{\text{acc}})\phi ^2+12\overline{P}_{\text{acc}}\right)\mathrm{cos}\vartheta +(1+\overline{P}_{\text{acc}})\phi ^2+1\right)r^3\beta \\ +\left(\overline{P}_{\text{acc}}(1\overline{P}_{\text{acc}})\phi ^2\mathrm{cos}\vartheta +\overline{P}_{\text{acc}}\phi ^2\right)r^4\end{array}\right]}.$$

Evaluating this at $`\beta =0`$ gives the integrated autocorrelation function

$$A_\varphi ^{\text{exp}}=\frac{\left((1\overline{P}_{\text{acc}})^2\phi ^2+12\overline{P}_{\text{acc}}\right)\mathrm{cos}\vartheta +(1\overline{P}_{\text{acc}})\phi ^2+1}{\overline{P}_{\text{acc}}\phi ^2\left((1\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +1\right)}.$$

For HMC where $`\vartheta =\pi /2`$ we have

$$F_\varphi ^{\text{exp}}(\beta ;\vartheta =\pi /2)=\frac{\left(\beta ^2+2r\beta +\left((1\overline{P}_{\text{acc}})\phi ^2+1\right)r^2\right)r}{\beta ^3+2r\beta ^2+(\phi ^2+1)r^2\beta +\overline{P}_{\text{acc}}r^3\phi ^2},$$

with the corresponding integrated autocorrelation function being

$$A_\varphi ^{\text{exp}}(\vartheta =\pi /2)=\frac{1\overline{P}_{\text{acc}}}{\overline{P}_{\text{acc}}}+\frac{1}{\overline{P}_{\text{acc}}\phi ^2}.$$

In the limit of unit acceptance rate, $`\overline{P}_{\text{acc}}=1`$, we obtain

$$F_\varphi ^{\text{exp}}(\beta ;\overline{P}_{\text{acc}}=1)=\frac{\left(\beta +(1\mathrm{cos}\vartheta )r\right)r}{\beta ^2+(1\mathrm{cos}\vartheta )r\beta +r^2\phi ^2}$$

and

$$A_\varphi ^{\text{exp}}(\overline{P}_{\text{acc}}=1)=\frac{1\mathrm{cos}\vartheta }{\phi ^2}.$$

Finally, for HMC in the limit of unit acceptance rate

$$F_\varphi ^{\text{exp}}(\beta ;\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=\frac{(r+\beta )r}{\beta ^2+r\beta +r^2\phi ^2},$$
(67)
and

$$A_\varphi ^{\text{exp}}(\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=\frac{1}{\phi ^2}.$$

#### 6.5.2 Quadratic Operators

The Laplace transform of the quadratic single mode autocorrelation function for exponentially distributed trajectories is obtained by substituting the explicit form for $`G^\pm `$ into equation (65), and we obtain

$$F_{(\varphi ^2)_c}^{\text{exp}}(\beta )=\frac{r\left[\begin{array}{c}\beta ^4+\left((\mathrm{cos}\vartheta )^2+(12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +4\right)r\beta ^3\\ +\left(\begin{array}{c}(1+2\overline{P}_{\text{acc}})(\mathrm{cos}\vartheta )^33(\mathrm{cos}\vartheta )^2\\ +(36\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +(42\overline{P}_{\text{acc}})\phi ^2+6\end{array}\right)r^2\beta ^2\\ +\left(\begin{array}{c}(4\overline{P}_{\text{acc}}2)(\mathrm{cos}\vartheta )^3+\left((4\overline{P}_{\text{acc}}4)\phi ^23\right)(\mathrm{cos}\vartheta )^2\\ +\left(2(1+\overline{P}_{\text{acc}})(\overline{P}_{\text{acc}}2)\phi ^2+36\overline{P}_{\text{acc}}\right)\mathrm{cos}\vartheta \\ +(84\overline{P}_{\text{acc}})\phi ^2+4\end{array}\right)r^3\beta \\ +\left(\begin{array}{c}\left(4(1+\overline{P}_{\text{acc}})^2\phi ^21+2\overline{P}_{\text{acc}}\right)(\mathrm{cos}\vartheta )^3\\ +\left((4\overline{P}_{\text{acc}}4)\phi ^21\right)(\mathrm{cos}\vartheta )^2\\ +\left(2(1+\overline{P}_{\text{acc}})(\overline{P}_{\text{acc}}2)\phi ^2+12\overline{P}_{\text{acc}}\right)\mathrm{cos}\vartheta \\ +(42\overline{P}_{\text{acc}})\phi ^2+1\end{array}\right)r^4\end{array}\right]}{\left[\begin{array}{c}\beta ^5+\left((\mathrm{cos}\vartheta )^2+(12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +4\right)r\beta ^4\\ +\left(\begin{array}{c}(1+2\overline{P}_{\text{acc}})(\mathrm{cos}\vartheta )^33(\mathrm{cos}\vartheta )^2\\ +(36\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +6+4\phi ^2\end{array}\right)r^2\beta ^3\\ +\left(\begin{array}{c}(4\overline{P}_{\text{acc}}2)(\mathrm{cos}\vartheta )^3+\left((2\overline{P}_{\text{acc}}4)\phi ^23\right)(\mathrm{cos}\vartheta )^2\\ +\left((44\overline{P}_{\text{acc}})\phi ^2+36\overline{P}_{\text{acc}}\right)\mathrm{cos}\vartheta \\ +(8+2\overline{P}_{\text{acc}})\phi ^2+4\end{array}\right)r^3\beta ^2\\ +\left(\begin{array}{c}(2(1+\overline{P}_{\text{acc}})(\overline{P}_{\text{acc}}2)\phi ^21+2\overline{P}_{\text{acc}})(\mathrm{cos}\vartheta )^3\\ +(4\phi ^21)(\mathrm{cos}\vartheta )^2\\ +\left(2(\overline{P}_{\text{acc}}+2)(1+\overline{P}_{\text{acc}})\phi ^2+12\overline{P}_{\text{acc}}\right)\mathrm{cos}\vartheta \\ +(4\overline{P}_{\text{acc}}+4)\phi ^2+1\end{array}\right)r^4\beta \\ +\left(\begin{array}{c}2\overline{P}_{\text{acc}}(1+\overline{P}_{\text{acc}})\phi ^2(\mathrm{cos}\vartheta )^32(\mathrm{cos}\vartheta )^2\overline{P}_{\text{acc}}\phi ^2\\ 2\overline{P}_{\text{acc}}(1+\overline{P}_{\text{acc}})\phi ^2\mathrm{cos}\vartheta +2\overline{P}_{\text{acc}}\phi ^2\end{array}\right)r^5\end{array}\right]}.$$

Evaluating this at $`\beta =0`$ gives the integrated autocorrelation function

$$A_{(\varphi ^2)_c}^{\text{exp}}=\frac{\left[\begin{array}{c}\left(4(1\overline{P}_{\text{acc}})^2\phi ^2+12\overline{P}_{\text{acc}}\right)(\mathrm{cos}\vartheta )^3+\left(4(1\overline{P}_{\text{acc}})\phi ^2+1\right)(\mathrm{cos}\vartheta )^2\\ +\left(2(1\overline{P}_{\text{acc}})(\overline{P}_{\text{acc}}2)\phi ^21+2\overline{P}_{\text{acc}}\right)\mathrm{cos}\vartheta +2(\overline{P}_{\text{acc}}2)\phi ^21\end{array}\right]}{2\overline{P}_{\text{acc}}\phi ^2(\mathrm{cos}\vartheta 1)(\mathrm{cos}\vartheta +1)\left((1\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +1\right)}.$$

For HMC where $`\vartheta =\pi /2`$ these equations simplify to

$$F_{(\varphi ^2)_c}^{\text{exp}}(\beta ;\vartheta =\pi /2)=\frac{\left(\beta ^2+2r\beta +\left(2(2\overline{P}_{\text{acc}})\phi ^2+1\right)r^2\right)r}{\beta ^3+2r\beta ^2+(4\phi ^2+1)r^2\beta +2\overline{P}_{\text{acc}}r^3\phi ^2}$$

and

$$A_{(\varphi ^2)_c}^{\text{exp}}(\vartheta =\pi /2)=\frac{2\overline{P}_{\text{acc}}}{\overline{P}_{\text{acc}}}+\frac{1}{2\overline{P}_{\text{acc}}\phi ^2}.$$

In the limit of unit acceptance rate, $`\overline{P}_{\text{acc}}=1`$, we obtain

$$F_{(\varphi ^2)_c}^{\text{exp}}(\beta ;\overline{P}_{\text{acc}}=1)=\frac{\left[\begin{array}{c}\beta ^2+\left(\mathrm{cos}\vartheta +2(\mathrm{cos}\vartheta )^2\right)r\beta \\ +\left((\mathrm{cos}\vartheta )^3(\mathrm{cos}\vartheta )^2\mathrm{cos}\vartheta +2\phi ^2+1\right)r^2\end{array}\right]r}{\left[\begin{array}{c}\beta ^3+\left(\mathrm{cos}\vartheta +2(\mathrm{cos}\vartheta )^2\right)r\beta ^2\\ +\left((\mathrm{cos}\vartheta )^3\mathrm{cos}\vartheta +1+4\phi ^2(\mathrm{cos}\vartheta )^2\right)r^2\beta \\ +\left(2(\mathrm{cos}\vartheta )^2\phi ^2+2\phi ^2\right)r^3\end{array}\right]}$$

and

$$A_{(\varphi ^2)_c}^{\text{exp}}(\overline{P}_{\text{acc}}=1)=\frac{(\mathrm{cos}\vartheta )^3(\mathrm{cos}\vartheta )^2\mathrm{cos}\vartheta +2\phi ^2+1}{2\phi ^2\left((\mathrm{cos}\vartheta )^21\right)}.$$

Finally, for HMC in the limit of unit acceptance rate

$$F_{(\varphi ^2)_c}^{\text{exp}}(\beta ;\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=\frac{\left(\beta ^2+2r\beta +(1+2\phi ^2)r^2\right)r}{\beta ^3+2r\beta ^2+(4\phi ^2+1)r^2\beta +2r^3\phi ^2},$$
(68)
and

$$A_{(\varphi ^2)_c}^{\text{exp}}(\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=1+\frac{1}{2\phi ^2}.$$

### 6.6 Fixed Length Trajectories

In this section we consider the case of fixed length trajectories, $`P_R(\tau )=\delta (\tau \overline{\tau })`$. In this case we find that

$$G_{j,k}^+=\overline{P}_{\text{acc}}G_{j,k}\text{and}G_{j,k}^{}=(1\overline{P}_{\text{acc}})G_{j,k},$$

without making any approximations, and from equation (64) we obtain

$$G_{j,k}^{\text{fix}}=e^{\beta \overline{\tau }}(\mathrm{cos}\phi )^j(\mathrm{sin}\phi )^k.$$

#### 6.6.1 Linear Operators

The Laplace transform of the linear single mode autocorrelation function for fixed length trajectories $`F_\varphi ^{\text{fix}}(\beta )`$ is obtained by substituting the explicit form for $`G^\pm `$ into equation (63),

$$\frac{(1+2\overline{P}_{\text{acc}})\mathrm{cos}\vartheta e^{2\beta \overline{\tau }}(\overline{P}_{\text{acc}}\mathrm{cos}\phi +1\overline{P}_{\text{acc}})e^{\beta \overline{\tau }}}{(12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta e^{2\beta \overline{\tau }}+\left((\overline{P}_{\text{acc}}+\overline{P}_{\text{acc}}\mathrm{cos}\phi 1)\mathrm{cos}\vartheta +\overline{P}_{\text{acc}}\mathrm{cos}\phi +1\overline{P}_{\text{acc}}\right)e^{\beta \overline{\tau }}1}.$$

Evaluating this at $`\beta =0`$ gives the integrated autocorrelation function

$$A_\varphi ^{\text{fix}}=\frac{(12\overline{P}_{\text{acc}})\mathrm{cos}\vartheta +\overline{P}_{\text{acc}}\mathrm{cos}\phi +1\overline{P}_{\text{acc}}}{\overline{P}_{\text{acc}}(\mathrm{cos}\vartheta +1)(\mathrm{cos}\phi 1)}.$$

For HMC where $`\vartheta =\pi /2`$ we have

$$F_\varphi ^{\text{fix}}(\beta ;\vartheta =\pi /2)=\frac{(\overline{P}_{\text{acc}}\mathrm{cos}\phi +1\overline{P}_{\text{acc}})e^{\beta \overline{\tau }}}{1(\overline{P}_{\text{acc}}\mathrm{cos}\phi +1\overline{P}_{\text{acc}})e^{\beta \overline{\tau }}},$$

with the corresponding integrated autocorrelation function being

$$A_\varphi ^{\text{fix}}(\vartheta =\pi /2)=\frac{\overline{P}_{\text{acc}}\mathrm{cos}\phi +1\overline{P}_{\text{acc}}}{\overline{P}_{\text{acc}}(1\mathrm{cos}\phi )}.$$

In the limit of unit acceptance rate, $`\overline{P}_{\text{acc}}=1`$, we obtain

$$F_\varphi ^{\text{fix}}(\beta ;\overline{P}_{\text{acc}}=1)=\frac{\mathrm{cos}\phi e^{\beta \overline{\tau }}\mathrm{cos}\vartheta e^{2\beta \overline{\tau }}}{(\mathrm{cos}\vartheta \mathrm{cos}\phi +\mathrm{cos}\phi )e^{\beta \overline{\tau }}\mathrm{cos}\vartheta e^{2\beta \overline{\tau }}1},$$

and

$$A_\varphi ^{\text{fix}}(\overline{P}_{\text{acc}}=1)=\frac{\mathrm{cos}\vartheta \mathrm{cos}\phi }{\mathrm{cos}\vartheta \mathrm{cos}\phi +\mathrm{cos}\phi \mathrm{cos}\vartheta 1}.$$

Finally, for HMC in the limit of unit acceptance rate

$$F_\varphi ^{\text{fix}}(\beta ;\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=\frac{\mathrm{cos}\phi e^{\beta \overline{\tau }}}{1\mathrm{cos}\phi e^{\beta \overline{\tau }}},$$
(69)
and

$$A_\varphi ^{\text{fix}}(\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=\frac{\mathrm{cos}\phi }{1\mathrm{cos}\phi }.$$

#### 6.6.2 Quadratic Operators

The Laplace transform of the quadratic single mode autocorrelation function for fixed length trajectories is obtained by substituting the explicit form for $`G^\pm `$ into equation (65), and we obtain

$$F_{(\varphi ^2)_c}^{\text{fix}}(\beta )=\frac{\left[\begin{array}{c}(\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^21+\overline{P}_{\text{acc}})e^{\beta \overline{\tau }}\\ e^{3\beta \overline{\tau }}(1+2\overline{P}_{\text{acc}})(\mathrm{cos}\vartheta )^3\\ +e^{2\beta \overline{\tau }}(2\overline{P}_{\text{acc}}+1+2\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2)(\mathrm{cos}\vartheta )^2\\ +e^{2\beta \overline{\tau }}(\overline{P}_{\text{acc}}+\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^21)\mathrm{cos}\vartheta \end{array}\right]}{\left[\begin{array}{c}((\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^21+\overline{P}_{\text{acc}})(\mathrm{cos}\vartheta )^2\\ +(2\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2+1)\mathrm{cos}\vartheta \overline{P}_{\text{acc}}(\mathrm{cos}\phi )^21+\overline{P}_{\text{acc}})e^{\beta \overline{\tau }}\\ +\left(e^{2\beta \overline{\tau }}+e^{2\beta \overline{\tau }}\overline{P}_{\text{acc}}+e^{2\beta \overline{\tau }}\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2+e^{3\beta \overline{\tau }}2e^{3\beta \overline{\tau }}\overline{P}_{\text{acc}}\right)(\mathrm{cos}\vartheta )^3\\ +e^{2\beta \overline{\tau }}\left(2\overline{P}_{\text{acc}}+1+2\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2\right)(\mathrm{cos}\vartheta )^2\\ +e^{2\beta \overline{\tau }}\left(\overline{P}_{\text{acc}}+\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^21\right)\mathrm{cos}\vartheta +1\end{array}\right]}.$$

Evaluating this at $`\beta =0`$ gives the integrated autocorrelation function

$$A_{(\varphi ^2)_c}^{\text{fix}}=\frac{(12\overline{P}_{\text{acc}})(\mathrm{cos}\vartheta )^2+2\mathrm{cos}\vartheta \overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^21+\overline{P}_{\text{acc}}}{(\mathrm{cos}\vartheta 1)(\mathrm{cos}\vartheta +1)(\mathrm{cos}\phi 1)(\mathrm{cos}\phi +1)\overline{P}_{\text{acc}}}.$$

For HMC where $`\vartheta =\pi /2`$ these equations simplify to

$$F_{(\varphi ^2)_c}^{\text{fix}}(\beta ;\vartheta =\pi /2)=\frac{(\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2+1\overline{P}_{\text{acc}})e^{\beta \overline{\tau }}}{1\left(\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2+1\overline{P}_{\text{acc}}\right)e^{\beta \overline{\tau }}},$$

and

$$A_{(\varphi ^2)_c}^{\text{fix}}(\vartheta =\pi /2)=\frac{\overline{P}_{\text{acc}}(\mathrm{cos}\phi )^2+1\overline{P}_{\text{acc}}}{\overline{P}_{\text{acc}}(1\mathrm{cos}\phi )(1+\mathrm{cos}\phi )}.$$

In the limit of unit acceptance rate, $`\overline{P}_{\text{acc}}=1`$, we obtain

$$F_{(\varphi ^2)_c}^{\text{fix}}(\beta ;\overline{P}_{\text{acc}}=1)=\frac{\left[\begin{array}{c}(\mathrm{cos}\phi )^2e^{\beta \overline{\tau }}+(\mathrm{cos}\vartheta )^3e^{3\beta \overline{\tau }}\\ e^{2\beta \overline{\tau }}(1+2(\mathrm{cos}\phi )^2)(\mathrm{cos}\vartheta )^2\\ e^{2\beta \overline{\tau }}(\mathrm{cos}\phi )^2\mathrm{cos}\vartheta \end{array}\right]}{\left[\begin{array}{c}\left((\mathrm{cos}\vartheta )^2(\mathrm{cos}\phi )^2+(1+2(\mathrm{cos}\phi )^2)\mathrm{cos}\vartheta +(\mathrm{cos}\phi )^2\right)e^{\beta \overline{\tau }}\\ +\left(e^{2\beta \overline{\tau }}(\mathrm{cos}\phi )^2+e^{3\beta \overline{\tau }}\right)(\mathrm{cos}\vartheta )^3\\ e^{2\beta \overline{\tau }}(1+2(\mathrm{cos}\phi )^2)(\mathrm{cos}\vartheta )^2e^{2\beta \overline{\tau }}(\mathrm{cos}\phi )^2\mathrm{cos}\vartheta 1\end{array}\right]},$$

and

$$A_{(\varphi ^2)_c}^{\text{fix}}(\overline{P}_{\text{acc}}=1)=\frac{(\mathrm{cos}\vartheta )^22\mathrm{cos}\vartheta (\mathrm{cos}\phi )^2+(\mathrm{cos}\phi )^2}{(\mathrm{cos}\vartheta )^2(\mathrm{cos}\phi )^2(\mathrm{cos}\vartheta )^2(\mathrm{cos}\phi )^2+1}.$$

Finally, for HMC in the limit of unit acceptance rate

$$F_{(\varphi ^2)_c}^{\text{fix}}(\beta ;\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=\frac{(\mathrm{cos}\phi )^2e^{\beta \overline{\tau }}}{1(\mathrm{cos}\phi )^2e^{\beta \overline{\tau }}},$$

and

$$A_{(\varphi ^2)_c}^{\text{fix}}(\vartheta =\pi /2,\overline{P}_{\text{acc}}=1)=\frac{(\mathrm{cos}\phi )^2}{1(\mathrm{cos}\phi )^2}.$$

## 7 Autocorrelations and Costs for HMC with $`P_{\text{acc}}1`$

In this section we calculate autocorrelations for the special case $`\vartheta =0`$ (HMC) and in the approximation where the acceptance rate is unity, $`\overline{P}_{\text{acc}}=1`$. We shall consider more general cases in the following section.

We begin with some general observations:

* *Integrated autocorrelation time.* It is immediately obvious that the integrated autocorrelation time may be obtained from the Laplace transform (6.4) by evaluating it at $`\beta =0`$.
* *Exponential autocorrelation time.* For an ergodic Markov process we can write the typical asymptotic form of the autocorrelation function as

$$C_\mathrm{\Omega }(t)\text{constant}\times e^{t/t_{\text{exp}}}\text{for }t\mathrm{}.$$

This exponential autocorrelation time is a different quantity from $`N_{\text{exp}}`$ of Section 6.1, which was defined in terms of Markov steps. $`t_{\text{exp}}`$ can be extracted from $`F_\mathrm{\Omega }(\beta )`$ by considering its analytic structure in the complex $`\beta `$ plane. Since $`t_{\text{exp}}`$ governs the most slowly decaying exponential, it follows from the definition of the Laplace transform (6.4) that $`F_\mathrm{\Omega }(\beta )`$ will have its rightmost pole at $`Re\beta =1/t_{\text{exp}}`$.
* *Dynamical critical exponent.* One of the most relevant measures of the effectiveness of an algorithm for studying continuum physics is the exponent $`z`$ relating the cost $``$ to the correlation length $`\xi `$ of the system, $`\xi ^z`$.
* *Optimal choice for $`\vartheta `$.* For the GHMC algorithm we should minimise the cost by varying both $`\tau `$ and $`\vartheta `$. For the case of quadratic operators with exponentially distributed trajectory lengths, for instance, the optimal choice of parameters when $`P_{\text{acc}}=1`$ is to take $`\vartheta 0`$ and $`\tau ={\scriptscriptstyle \frac{1}{2}}\vartheta ^20`$. However, this ignores the fact that the cost does not decrease when we take $`\tau `$ smaller than the $`\delta \tau `$ required to obtain a reasonable Metropolis acceptance rate. If we choose $`\overline{\tau }_{\text{opt}}=\delta \tau `$ (L2MC) and the corresponding value for $`\vartheta _{\text{opt}}`$ we find that the cost is less than for the HMC case, but only by a constant factor. As the cost is only defined up to a implementation dependent constant factor anyhow we may conclude that generalised HMC does not appear to promise great improvements over HMC. The situation is more complex in the real world where $`P_{\text{acc}}1`$, which is explored in Section 8.

### 7.1 Exponentially Distributed Trajectory Lengths

#### 7.1.1 Linear Operators

The least uninteresting linear operator in free scalar field theory is the *magnetisation* $`M=_x\varphi (x)`$. From equations (5) and (7) this is simply the zero-momentum mode in Fourier space, $`\omega =m`$, and is thus expected to evolve most slowly in fictitious time.

###### Example 7.1.1

From equation (41), with $`\mathrm{\Omega }_{pq}=\delta _{p0}`$, and equation (67), the Laplace transform of the autocorrelation function for the magnetisation is

$$F_M(\beta )=\frac{r(r+\beta )}{\beta (r+\beta )+m^2}.$$

As explained previously, the exponential autocorrelation time $`t_{\text{exp}}=1/Re\beta _{\text{exp}}`$ corresponds to the rightmost pole in $`F_M`$, and this occurs when

$$\beta _{\text{exp}}^{}{}_{}{}^{\pm }=\frac{r}{2}\pm \sqrt{\left(\frac{r}{2}\right)^2m^2},$$

and therefore

$$t_{\text{exp}}=\{\begin{array}{cc}\frac{2\overline{\tau }}{1\sqrt{1(2m\overline{\tau })^2}}\hfill & \text{if }\overline{\tau }\frac{1}{2m}\hfill \\ 2\overline{\tau }\hfill & \text{if }\overline{\tau }\frac{1}{2m}\text{,}\hfill \end{array}$$

where we have used the fact that $`r=1/\overline{\tau }`$. We observe that the exponential autocorrelation time is minimised when the average trajectory length is chosen to be $`\overline{\tau }=1/2m`$. Note that $`M`$ only couples to the zero momentum mode, and its relaxation rate determines $`t_{\text{exp}}`$ in this case.

The integrated autocorrelation function $`A_M`$ is given by

$$A_M=F_M(0)=\left(\frac{1}{m\overline{\tau }}\right)^2,$$

and we can minimise the cost of computing the magnetisation by making use of equation (40). The optimal trajectory length is $`\overline{\tau }_{\text{opt}}=\sqrt{2}/m`$, which corresponds to the minimum integrated autocorrelation function value $`A_M(\overline{\tau }_{\text{opt}})=1/2`$. Note that the optimal trajectory length does not minimise the exponential autocorrelation time — they differ by a factor of $`2\sqrt{2}`$.

The correlation length for this system is $`\xi 1/m`$, and our result indicates that the optimal trajectory length is $`\overline{\tau }_{\text{opt}}\xi ^z`$ with a dynamical critical exponent $`z=1`$. Indeed, if we were to choose an average trajectory length $`\overline{\tau }\xi ^\alpha `$ then we would find that the cost per effectively independent configuration would grow as

$$\{\begin{array}{cc}\xi ^\alpha \hfill & \text{for }\alpha 1,\xi \mathrm{}\hfill \\ \xi ^{2\alpha }\hfill & \text{for }\alpha 1,\xi \mathrm{}\text{.}\hfill \end{array}$$

Keeping the average trajectory length fixed as the correlation length $`\xi `$ increases, i.e., choosing $`\alpha =0`$, thus leads to a dynamical critical exponent $`z=2`$.

We can calculate the autocorrelation function $`C_M(t)`$ in closed form by inverting the Laplace transform. If we expand $`F_M`$ in partial fractions

$$F_M(\beta )=\frac{1+\mathrm{\Delta }}{\mathrm{\Delta }(1\mathrm{\Delta }+2\overline{\tau }\beta )}\frac{1\mathrm{\Delta }}{\mathrm{\Delta }(1+\mathrm{\Delta }+2\overline{\tau }\beta )},$$

where $`\mathrm{\Delta }\sqrt{1(2m\overline{\tau })^2}`$, then it is easy to verify that the autocorrelation function is

$`C_M(t)`$ $`=`$ $`{\displaystyle \frac{1}{2\overline{\tau }\mathrm{\Delta }}}\left\{(1+\mathrm{\Delta })e^{(1\mathrm{\Delta })t/2\overline{\tau }}(1\mathrm{\Delta })e^{(1+\mathrm{\Delta })t/2\overline{\tau }}\right\}`$
$`=`$ $`{\displaystyle \frac{1}{\overline{\tau }}}e^{t/2\overline{\tau }}\left\{\mathrm{cosh}\left({\displaystyle \frac{\mathrm{\Delta }t}{2\overline{\tau }}}\right)+{\displaystyle \frac{1}{\mathrm{\Delta }}}\mathrm{sinh}\left({\displaystyle \frac{\mathrm{\Delta }t}{2\overline{\tau }}}\right)\right\}`$
when $`0<\mathrm{\Delta }<1`$. When $`\mathrm{\Delta }^2<0`$ we have

$$C_M(t)=\frac{1}{\overline{\tau }}e^{t/2\overline{\tau }}\left\{\mathrm{cos}\left(\frac{\mathrm{\Delta }^{}t}{2\overline{\tau }}\right)+\frac{1}{\mathrm{\Delta }^{}}\mathrm{sin}\left(\frac{\mathrm{\Delta }^{}t}{2\overline{\tau }}\right)\right\},$$

where $`\mathrm{\Delta }^{}i\mathrm{\Delta }=\sqrt{(2m\overline{\tau })^21}`$. We observe that there are oscillations in the autocorrelation function for long trajectories where $`\overline{\tau }>1/2m`$. For the critical case $`\mathrm{\Delta }=0`$ we have

$$C_M(t)=\frac{1}{\overline{\tau }}\left(1+\frac{t}{2\overline{\tau }}\right)e^{t/2\overline{\tau }}.$$

Finally, the autocorrelation function expressed in terms of Markov steps is $`C_M(n)=e^{n/T_0}`$, where the exponential autocorrelation time $`N_{\text{exp}}=T_0=1/\mathrm{ln}[1+(m\overline{\tau })^2]`$.

#### 7.1.2 Quadratic Operators

###### Example 7.1.2

We obtain the Laplace transform of the connected autocorrelation function555This is a synonym for the autocorrelation function for $`M_c^2`$, the connected part of $`M^2`$. for $`M^2`$ by setting $`\mathrm{\Omega }_{pq}=\delta _{p0}\delta _{q0}`$ in equation (43). From equation (68)

$$F_{M_c^2}(\beta )=\frac{r^3+2r^2\beta +r\beta ^2+2m^2r}{\beta ^3+2r\beta ^2+(r^2+4m^2)\beta +2m^2r}.$$
(70)
The integrated autocorrelation function is thus

$$A_{M_c^2}=F_{M_c^2}(0)=1+\frac{1}{2}\left(\frac{1}{m\overline{\tau }}\right)^2.$$

Minimising the cost by means of equation (40), we obtain $`\overline{\tau }_{\text{opt}}=1/\left(\sqrt{3}m\right)`$ and thence $`A_{M^2}(\overline{\tau }_{\text{opt}})=5/2`$. Again, the dynamical critical exponent is $`z=1`$. The optimum trajectory length is different from that for $`M`$, which is to be expected.

For exponentially distributed trajectory lengths the Laplace transform of the autocorrelation function is a rational function in $`\beta `$ with the numerator of lower degree than the denominator (see, eg, equations (67) and (68)), which implies that the autocorrelation function is a sum of exponentials. In this case the exponents are the roots of the cubic denominator, and they are either all real, or one is real and the other two are complex conjugates depending on the value of the mean trajectory length $`1/r`$. This is shown explicitly in Appendix A.

###### Example 7.1.3

Following equation (68), it is easy to write down the the Laplace transform of the connected autocorrelation function for the energy $`E=_p{\scriptscriptstyle \frac{1}{2}}\omega _p^2\varphi _p^2`$

$$F_E(\beta )=\frac{1}{N}\underset{p=1}{\overset{N}{}}\frac{r^3+2r^2\beta +r\beta ^2+2\omega _p^2r}{\beta ^3+2r\beta ^2+(r^2+4\omega _p^2)\beta +2\omega _p^2r}.$$
(71)
The integrated autocorrelation function for the energy is

$$A_E=F_E(0)=1+\frac{1}{2\overline{\tau }^2}\underset{p=1}{\overset{N}{}}\frac{1}{\omega _p^2}1+\frac{1}{2\overline{\tau }^2}\frac{1}{N}\sigma _{m^2}^{(1)},$$

and the optimal trajectory length is $`\overline{\tau }_{\text{opt}}=\sqrt{\sigma _{m^2}^{(1)}/3}`$, leading to an integrated autocorrelation function value of $`A_E(\overline{\tau }_{\text{opt}})=5/2`$.

For one dimensional free field theory it remains to evaluate the spectral sum $`\sigma _{m^2}^{(1)}`$, details of which are discussed in Appendix C. In the physically interesting limit of small $`m`$ and large $`V`$, we find $`\sigma _{m^2}^{(1)}{\scriptscriptstyle \frac{1}{2}}m`$, and thus $`\overline{\tau }_{\text{opt}}1/\sqrt{6m}`$. Hence the dynamical critical exponent for the energy is $`{\scriptscriptstyle \frac{1}{2}}`$.

For two dimensional free field theory we get $`z=0`$ up to logarithmic corrections (see Appendix D).

### 7.2 Fixed Length Trajectories

#### 7.2.1 Linear Operators

From Section 6.6 the Laplace transform of the autocorrelation function for the general linear operator of equation (41) is easily expressed in terms of equation (69). The exponential autocorrelation time $`t_{\text{exp}}`$ is related to the rightmost pole of the Laplace transform of the autocorrelation function,

$$t_{\text{exp}}=1/Re\beta _{\text{exp}}.$$

Equation (69) has poles in the complex $`\beta `$ plane for

$$\beta =\frac{1}{\overline{\tau }}\mathrm{ln}\mathrm{cos}\omega _p\overline{\tau },$$

and hence

$$Re\beta =\frac{1}{\overline{\tau }}\mathrm{ln}|\mathrm{cos}\omega _p\overline{\tau }|,$$

which leads to

$$t_{\text{exp}}=\underset{p}{\mathrm{max}}\frac{\overline{\tau }}{\left|\mathrm{ln}|\mathrm{cos}\omega _p\overline{\tau }|\right|}.$$

This diverges when $`\omega _p\overline{\tau }/\pi \mathrm{}`$ for any $`p`$, which just reflects the fact that the Hybrid Monte Carlo algorithm is not ergodic for these cases, as was first pointed out by Mackenzie . Perhaps this is most simply understood by considering the trajectory of the harmonic oscillator with frequency $`\omega _p`$ in the $`(\omega _p\varphi _p,\pi _p)`$ phase space. The Molecular Dynamics trajectories are circular arcs subtending an angle of $`\omega _p\overline{\tau }`$ about the origin, and the momentum refreshment corresponds to a change of the $`\pi _p`$ coordinate leaving the $`\omega _p\varphi _p`$ coordinate unchanged. If $`\omega _p\overline{\tau }`$ is an integer multiple of $`\pi `$ the value of $`\varphi _p`$ can at most change sign, and thus this mode certainly cannot explore the whole of its phase space.

###### Example 7.2.1

For the magnetisation the integrated autocorrelation function is

$$A_M=F_M(0)=\frac{\mathrm{cos}m\overline{\tau }}{1\mathrm{cos}m\overline{\tau }}.$$

If we minimise the cost we find that the optimal trajectory length corresponds to $`m\overline{\tau }`$ being an odd multiple of $`\pi `$, and the worst case occurs when $`m\overline{\tau }`$ is an even multiple of $`\pi `$. Both cases just reflect the non-ergodic nature of the updating scheme discussed above: the fact that the optimal “ergodic” update occurs when $`m\overline{\tau }`$ is taken very close to an odd multiple of $`\pi `$ is just an “accidental” consequence of the fact that $`M=0`$.

#### 7.2.2 Quadratic Operators

###### Example 7.2.2

In the case of the quadratic operator $`M_c^2`$ we find from equation (6.6.2) that

$$F_{M_c^2}(\beta )=\frac{\mathrm{cos}^2m\overline{\tau }}{e^{\beta \overline{\tau }}\mathrm{cos}^2m\overline{\tau }},$$

and thus $`A_{M_c^2}=F_{M_c^2}(0)=\mathrm{cot}^2m\overline{\tau }`$ for fixed length trajectories. The optimal trajectory length $`m\overline{\tau }_{\text{opt}}1.3866`$. As is to be expected, the non-ergodicity at $`m\overline{\tau }/\pi \mathrm{}`$ manifests itself as a divergence in $`A_{M_c^2}`$.

###### Example 7.2.3

For the energy of one dimensional free field theory we obtain

$$F_E(\beta )=\frac{1}{N}\underset{p}{}\frac{\mathrm{cos}^2\omega _p\overline{\tau }}{e^{\beta \overline{\tau }}\mathrm{cos}^2\omega _p\overline{\tau }},$$

and therefore $`A_E=F_E(0)=\frac{1}{N}_p\mathrm{cot}^2\omega _p\overline{\tau }`$ which diverges whenever $`\omega _p\overline{\tau }/\pi \mathrm{}`$ for any $`p`$.

## 8 Comparison of Costs for $`P_{\text{acc}}1`$

We wish to compare the performance of the HMC, L2MC and GHMC algorithms for one dimensional free field theory. To do this we shall compare the cost of generating a statistically independent measurement of the magnetisation $`M`$ and the magnetic susceptibility $`M_c^2`$, choosing the optimal values for the angle $`\vartheta `$ and the average trajectory length $`\overline{\tau }`$.

Following equation (40) we can minimise the cost with respect to $`\vartheta `$ without having to specify the form of the refresh distribution.

The next step is to minimise the cost with respect to the average trajectory length $`\phi =m\overline{\tau }`$. Strictly speaking we should note that the acceptance probability $`\overline{P}_{\text{acc}}`$ is a function of $`\overline{\tau }`$, but to a good approximation we may assume that $`P_{\text{acc}}`$ depends only upon the integration step size $`\delta \tau `$ *except* in the case of very short trajectories, such as Langevin-type algorithms (see Figure 1). The exact solution of the problem would clearly be highly transcendental. For this reason we shall find that although L2MC is a special case of GHMC and thus can never be cheaper, our “optimal” solution666*I.e.*, the solution obtained by neglecting the dependence of $`\overline{P}_{\text{acc}}`$ on $`\overline{\tau }`$ when optimising the parameters $`\vartheta `$ and $`\overline{\tau }`$. for very short trajectories (acceptance probabilities very close to unity) will in fact cost more than L2MC. Fortunately, this is irrelevant because the minimum cost for L2MC is far greater than the minimum for GHMC, the latter occurring for long trajectory lengths where our approximations are very good.

Another implicit approximation we make is that we treat $`\delta \tau `$ and $`\overline{\tau }`$ as independent parameters, although the trajectory length must be an integer multiple of the step size; again this is a very good approximation except for near to the Langevin limit. Of course we are also ignoring subtleties such as that a multistep leapfrog integration is cheaper than a series of single steps, that acceptance tests are a significant fraction of the cost for very short trajectories, and that HMC requires less memory than GHMC because the old momenta need not be saved. All of these facts would disfavour L2MC, but we shall see that it is not competitive even without these factors being taken into account.

### 8.1 Linear Operators

Applying equation (40) to equation (63) we find that the cost is minimised by choosing777Setting $`\vartheta =0`$ corresponds to never refreshing the momenta, and thus to a non-ergodic algorithm. This is just a peculiarity of linear operators in free field theory, and we can instead consider choosing $`\vartheta `$ to be very small but non-zero to circumvent this difficulty. $`\vartheta _{\text{opt}}=0`$, and at this value the Laplace transform of the single mode autocorrelation function becomes

$$F_M(\beta ;\vartheta _{\text{opt}})=\frac{G_{1,0}^{+}{}_{}{}^{2}+G_{0,0}^{}{}_{}{}^{2}G_{0,1}^{+}{}_{}{}^{2}+G_{1,0}^++G_{0,0}^{}}{1+2G_{1,0}^+G_{1,0}^{+}{}_{}{}^{2}+G_{0,0}^{}{}_{}{}^{2}G_{0,1}^{+}{}_{}{}^{2}}.$$

#### 8.1.1 Exponentially Distributed Trajectory Lengths

The integrated autocorrelation function for exponentially distributed trajectory lengths is obtained by substituting the results (66) for the integrals of equation (64) into $`F_M(\beta ;\vartheta =\vartheta _{\text{opt}})`$ and setting $`\beta =0`$. The cost becomes

$$_M^{\text{exp}}(\delta \tau ,\vartheta _{\text{opt}},\phi )=\frac{V}{m\delta \tau }\left[\frac{(\overline{P}_{\text{acc}}2)\phi }{\overline{P}_{\text{acc}}}+\frac{4(\overline{P}_{\text{acc}}1)}{\overline{P}_{\text{acc}}(\overline{P}_{\text{acc}}2)\phi }\right].$$

Minimising this with respect to $`\phi `$ we find that

$$\phi _{\text{opt}}=\frac{\sqrt{1\overline{P}_{\text{acc}}}}{1{\scriptscriptstyle \frac{1}{2}}\overline{P}_{\text{acc}}}$$

and that the cost at the optimal parameters is

$$_M^{\text{exp}}(\delta \tau ,\vartheta _{\text{opt}},\phi _{\text{opt}})=\frac{4V(1\overline{P}_{\text{acc}})}{m\delta \tau \overline{P}_{\text{acc}}\sqrt{1\overline{P}_{\text{acc}}}}$$

which corresponds to a dynamical critical exponent $`z=1`$.

For the lowest order leapfrog integration scheme where we must scale $`\delta \tau V^{{\scriptscriptstyle \frac{1}{4}}}`$ to keep $`\overline{P}_{\text{acc}}`$ constant we thus find $`_M^{\text{exp}}V^{{\scriptscriptstyle \frac{5}{4}}}`$ as expected.

### 8.2 Quadratic Operators

We can make the minimisation problem for quadratic operators manifestly algebraic by writing $`c\mathrm{cos}\vartheta `$ in equation (65). The condition for the cost to be a minimum is that $`c_{\text{opt}}`$ is a root of the polynomial

$`G_{1,1}^{+}{}_{}{}^{2}\left(G_{0,0}^{}+G_{0,2}^++G_{2,0}^+\right)^2c_{\text{opt}}^{}{}_{}{}^{4}`$
$`\left(2G_{1,1}^{+}{}_{}{}^{2}G_{0,2}^+G_{2,0}^++G_{0,2}^{+}{}_{}{}^{2}+G_{0,2}^+G_{0,0}^{}\right)^2c_{\text{opt}}^{}{}_{}{}^{3}`$
$`+\left(2G_{0,2}^{+}{}_{}{}^{2}G_{0,0}^{}2G_{0,2}^{+}{}_{}{}^{3}2G_{1,1}^{+}{}_{}{}^{2}G_{0,0}^{}2G_{2,0}^+G_{1,1}^{+}{}_{}{}^{2}+2G_{0,2}^{+}{}_{}{}^{2}G_{2,0}^+6G_{1,1}^{+}{}_{}{}^{2}G_{0,2}^+\right)c_{\text{opt}}^{}{}_{}{}^{2}`$
$`G_{0,2}^{+}{}_{}{}^{2}c_{\text{opt}}+G_{1,1}^{+}{}_{}{}^{2}`$ (72)
lying in the interval $`[1,1]`$.

#### 8.2.1 Exponentially Distributed Trajectory Lengths

Just as in equations (40) the extrema of the cost occur on the ideal defined by the polynomials $`d/dc`$ and $`d/d\phi `$. Finding the Gröbner basis with respect to the purely lexicographical ordering with $`c<\phi `$ we find the point $`(c_{\text{opt}},\phi _{\text{opt}})`$ at which the cost is minimal is defined by the equations

$`0`$ $`=`$ $`(4+3\overline{P}_{\text{acc}})c_{\text{opt}}^{}{}_{}{}^{4}+4(12\overline{P}_{\text{acc}})(1\overline{P}_{\text{acc}})c_{\text{opt}}^{}{}_{}{}^{3}`$ (73)
$`+16(1\overline{P}_{\text{acc}})c_{\text{opt}}^{}{}_{}{}^{2}+4c_{\text{opt}}+\overline{P}_{\text{acc}}4`$
$`0`$ $`=`$ $`4\overline{P}_{\text{acc}}^{}{}_{}{}^{2}(2\overline{P}_{\text{acc}})(4\overline{P}_{\text{acc}})(\overline{P}_{\text{acc}}^{}{}_{}{}^{2}6\overline{P}_{\text{acc}}+6)\phi _{\text{opt}}^{}{}_{}{}^{2}`$ (74)
$`+(4+3\overline{P}_{\text{acc}})(\overline{P}_{\text{acc}}^{}{}_{}{}^{3}+2\overline{P}_{\text{acc}}^{}{}_{}{}^{2}10\overline{P}_{\text{acc}}+8)c_{\text{opt}}^{}{}_{}{}^{5}`$
$`+2(1+\overline{P}_{\text{acc}})(4\overline{P}_{\text{acc}}^{}{}_{}{}^{4}+6\overline{P}_{\text{acc}}^{}{}_{}{}^{3}47\overline{P}_{\text{acc}}^{}{}_{}{}^{2}+68\overline{P}_{\text{acc}}32)c_{\text{opt}}^{}{}_{}{}^{4}`$
$`+(38\overline{P}_{\text{acc}}^{}{}_{}{}^{4}+108\overline{P}_{\text{acc}}^{}{}_{}{}^{3}16\overline{P}_{\text{acc}}^{}{}_{}{}^{2}152\overline{P}_{\text{acc}}+96)c_{\text{opt}}^{}{}_{}{}^{3}`$
$`+(16\overline{P}_{\text{acc}}^{}{}_{}{}^{5}+56\overline{P}_{\text{acc}}^{}{}_{}{}^{4}20\overline{P}_{\text{acc}}^{}{}_{}{}^{3}168\overline{P}_{\text{acc}}^{}{}_{}{}^{2}+248\overline{P}_{\text{acc}}96)c_{\text{opt}}^{}{}_{}{}^{2}`$
$`+(35\overline{P}_{\text{acc}}^{}{}_{}{}^{4}110\overline{P}_{\text{acc}}^{}{}_{}{}^{3}+54\overline{P}_{\text{acc}}^{}{}_{}{}^{2}+88\overline{P}_{\text{acc}}64)c_{\text{opt}}`$
$`+8\overline{P}_{\text{acc}}^{}{}_{}{}^{5}60\overline{P}_{\text{acc}}^{}{}_{}{}^{4}+126\overline{P}_{\text{acc}}^{}{}_{}{}^{3}62\overline{P}_{\text{acc}}^{}{}_{}{}^{2}48\overline{P}_{\text{acc}}+32`$
Using Sturm sequences we may easily show that equation (73) has exactly one real root in $`[0,1]`$ and none in $`[1,0]`$, and obviously equation (74) has exactly one positive solution for $`\phi _{\text{opt}}`$. The cost at the point $`(c_{\text{opt}},\phi _{\text{opt}})`$ is given by

$$_{M_c^2}^{\text{exp}}(\delta \tau ,\vartheta _{\text{opt}},\phi )=\frac{V\left(\begin{array}{c}(7\overline{P}_{\text{acc}}3\overline{P}_{\text{acc}}^{}{}_{}{}^{2}4)\phi _{\text{opt}}^{}{}_{}{}^{2}c_{\text{opt}}^{}{}_{}{}^{3}+(2\overline{P}_{\text{acc}}+1)c_{\text{opt}}+1\\ +(2\overline{P}_{\text{acc}}1)c_{\text{opt}}^{}{}_{}{}^{3}c_{\text{opt}}^{}{}_{}{}^{2}+(\overline{P}_{\text{acc}}+4)\phi _{\text{opt}}^{}{}_{}{}^{2}\\ +(\overline{P}_{\text{acc}}^{}{}_{}{}^{2}5\overline{P}_{\text{acc}}+4)c_{\text{opt}}\phi _{\text{opt}}^{}{}_{}{}^{2}+(4+3\overline{P}_{\text{acc}})c_{\text{opt}}^{}{}_{}{}^{2}\phi _{\text{opt}}^{}{}_{}{}^{2}\end{array}\right)}{\left(\begin{array}{c}\overline{P}_{\text{acc}}\delta \tau m(\overline{P}_{\text{acc}}1)c_{\text{opt}}^{}{}_{}{}^{3}\phi _{\text{opt}}+\phi _{\text{opt}}\overline{P}_{\text{acc}}\delta \tau m\\ \phi _{\text{opt}}\overline{P}_{\text{acc}}\delta \tau mc_{\text{opt}}^{}{}_{}{}^{2}\overline{P}_{\text{acc}}\delta \tau m(1+\overline{P}_{\text{acc}})c_{\text{opt}}\phi _{\text{opt}}\end{array}\right)}$$

This solution is a function of $`\delta \tau `$ and $`\overline{P}_{\text{acc}}`$ which are not independent variables, and using the results (29), (30) and (33) we can compute the cost as a function of $`\overline{P}_{\text{acc}}`$ as shown in Figure 2.

#### 8.2.2 HMC

The Hybrid Monte Carlo algorithm corresponds to setting $`\vartheta =\pi /2`$, and thus we find that the optimal trajectory length in this case is $`\phi _{\text{opt}}=1/\sqrt{4\overline{P}_{\text{acc}}}`$, corresponding to a cost

$$_{\text{opt}}=\frac{2V\sqrt{4\overline{P}_{\text{acc}}}}{\overline{P}_{\text{acc}}\delta \tau m}.$$

This is also shown in Figure 2.

### 8.3 Fixed Length Trajectories

For fixed length trajectories we shall only analyse the case of L2MC for which the trajectory length $`\phi =m\delta \tau `$. In this case $`c_{\text{opt}}`$ satisfies

$$\frac{\overline{P}_{\text{acc}}c_{\text{opt}}^{}{}_{}{}^{2}(\mathrm{cos}\phi _{\text{opt}})^2+3c_{\text{opt}}^{}{}_{}{}^{2}\overline{P}_{\text{acc}}+\overline{P}_{\text{acc}}(\mathrm{cos}\phi _{\text{opt}})^2\overline{P}_{\text{acc}}2c_{\text{opt}}^{}{}_{}{}^{2}4c_{\text{opt}}\overline{P}_{\text{acc}}(\mathrm{cos}\phi _{\text{opt}})^2+2}{\overline{P}_{\text{acc}}\left(c_{\text{opt}}^{}{}_{}{}^{2}(\mathrm{cos}\phi _{\text{opt}})^2c_{\text{opt}}^{}{}_{}{}^{2}(\mathrm{cos}\phi _{\text{opt}})^2+1\right)}=0.$$

This, together with the corresponding cost, is also plotted in Figure 2. From this figure it is clear that the minimum cost occurs for $`\overline{P}_{\text{acc}}`$ very close to unity, where the scaling variable $`x=V\delta \tau ^6`$ is very small. We may then express $`c_{\text{opt}}`$ as

$$c_{\text{opt}}=1m\left(\frac{x}{V}\right)^{1/6}+\frac{1}{2}m^2\left(\frac{x}{V}\right)^{1/3}\frac{11}{24}m^3\left(\frac{x}{V}\right)^{1/2}+\mathrm{},$$

and likewise

$$\overline{P}_{\text{acc}}=1\sqrt{\frac{5}{8\pi }}\left(1+\frac{9m^2}{20}+\frac{39m^4}{800}+O(m^6)\right)\sqrt{x}+O(x^{3/2}).$$

From these relations we find that the minimum cost for L2MC is

$$_{\text{opt}}=2\left(\frac{10}{\pi }\right)^{1/4}V^{5/4}m^{3/2}.$$

This result tells us that not only does the tuned L2MC algorithm have a dynamical critical exponent $`z=3/2`$, but also it has a volume dependence of exactly the same form as HMC. We may understand why this behaviour occurs rather than the naive  $`V^{7/6}m^1`$ by the following simple argument.

If $`\overline{P}_{\text{acc}}<1`$ then the system will carry out a random walk backwards and forwards along a trajectory because the momentum, and thus the direction of travel, must be reversed upon a Metropolis rejection. A simple minded analysis is that the average time between rejections must be $`O(1/m)`$ in order to achieve $`z=1`$. This time is approximately

$$\underset{n=0}{\overset{\mathrm{}}{}}\overline{P}_{\text{acc}}^{}{}_{}{}^{n}(1\overline{P}_{\text{acc}})n\delta \tau =\frac{\overline{P}_{\text{acc}}\delta \tau }{1\overline{P}_{\text{acc}}}=\frac{1}{m}.$$

For small $`\delta \tau `$ we have $`1\overline{P}_{\text{acc}}=erf\sqrt{kV\delta \tau ^6}\sqrt{V\delta \tau ^6}`$ where $`k`$ is a constant, and hence we must scale $`\delta \tau `$ so as to keep $`V\delta \tau ^4/m^2`$ fixed. Since the L2MC algorithm has a naive dynamical critical exponent $`z=1`$, this means that the cost should vary as $`V(V\delta \tau ^4m^2)^{1/4}/m\delta \tau =V^{5/4}m^{3/2}`$.

## 9 Autocorrelations and Frequency of Measurement

In this final section we perform an elementary analysis of the general problem of determining the optimal frequency for making measurements of observables on a Markov chain.

Suppose we wish to measure the expectation value $`\mathrm{\Omega }`$ of an operator $`\mathrm{\Omega }`$ by means of a Markov process (not necessarily HMC). If the cost888*I.e.*, the cost measured in units of computer time. of making one Markov step is $`_S`$, and the cost of making one measurement of $`\mathrm{\Omega }`$ is $`_M`$, how often should we make measurements in order to minimise the cost of the calculation? Due to the presence of correlations between successive measurements, the answer is not quite trivial.

Consider a sequence of $`N_u`$ uncorrelated measurements of $`\mathrm{\Omega }`$. The sample variance $`\overline{V}_\mathrm{\Omega }^u`$ is related to the intrinsic variance $`V_\mathrm{\Omega }`$ in the distribution of $`\mathrm{\Omega }`$ in the usual way

$$\overline{V}_\mathrm{\Omega }^u\left(\overline{\mathrm{\Omega }}\overline{\mathrm{\Omega }}\right)^2=\frac{\left(\mathrm{\Omega }\mathrm{\Omega }\right)^2}{N_u}\frac{V_\mathrm{\Omega }^u}{N_u}.$$
(75)
For the general case of $`N_c`$ successive correlated measurements, from equation (36), the sample variance $`\overline{V}_\mathrm{\Omega }`$ is

$$\overline{V}_\mathrm{\Omega }^u\left(\overline{\mathrm{\Omega }}\overline{\mathrm{\Omega }}\right)^2=(1+2A_\mathrm{\Omega })\frac{\left(\mathrm{\Omega }\mathrm{\Omega }\right)^2}{N_c}=(1+2A_\mathrm{\Omega })\frac{V_\mathrm{\Omega }}{N_c}$$
(76)
so that on average $`1+2A_\mathrm{\Omega }`$ correlated measurements are needed to reduce the variance by the same amount as a single truly independent measurement. If the cost of measuring $`\mathrm{\Omega }`$ is small (large), then it should be beneficial to make more (less) than one measurement per decorrelation time.

Let $`T`$ be the total number of Markov steps required to generate one independent sample, and let $`K`$ be the number of Markov steps between each measurement of $`\mathrm{\Omega }`$, so that the number of measurements performed per independent sample is $`T/K`$. If we wish to make sufficient measurements to generate $`N`$ independent samples, from equations (75) and (76) we have

$$\frac{1+2A_K}{T/K}=1,$$
(77)
where $`A_K`$ is the integrated autocorrelation function for measurements separated by $`K`$ Markov steps. The total cost $``$ of the entire process (steps plus measurements) is clearly

$$=N\left(_ST+_M\frac{T}{K}\right).$$
(78)
Eliminating $`T`$ from equations (77) and (78) we obtain

$$=N\left(_S+\frac{_M}{K}\right)\left(1+2A_K\right)K.$$

The integrated autocorrelation function $`A_K`$ can be written in terms of the autocorrelation function $`C_\mathrm{\Omega }(Kl)`$ of section (6.1). Assuming equation (35) we have

$$A_K=\underset{l=1}{\overset{\mathrm{}}{}}C_\mathrm{\Omega }(Kl)=\underset{l=1}{\overset{\mathrm{}}{}}e^{Kl/N_{\text{exp}}}=\frac{1}{e^{K/N_{\text{exp}}}1}$$

and so

$$=\left(K_S+_M\right)\left(\frac{e^{K/N_{\text{exp}}}+1}{e^{K/N_{\text{exp}}}1}\right).$$

Minimising with respect to $`K`$ gives

$$\mathrm{sinh}x=x+\alpha $$
(79)
where

$$x\frac{K}{N_{\text{exp}}}\text{and}\alpha \frac{_M}{K_S}.$$

For small $`\alpha `$ where measurements are “cheap” we find

$$K_{\mathrm{opt}}=N_{\text{exp}}\left(\frac{6_M}{_SN_{\text{exp}}}\right)^{1/3},$$

while for large $`\alpha `$ where they are “expensive”, we obtain

$$K_{\mathrm{opt}}=N_{\text{exp}}\mathrm{ln}\left(\frac{2_M}{_SN_{\text{exp}}}\right),$$

the slow logarithmic increase of $`K_{\mathrm{opt}}`$ being due to the exponential decay in the autocorrelation function $`C_\mathrm{\Omega }(l)`$. The crossover point, $`K=1`$, occurs when $`\alpha 0.1752`$.

## 10 Conclusions

We have introduced some powerful techniques for analysing a wide class of algorithms for free field theory. Whereas inventing algorithms which are efficacious for free field theory but useless for interesting ones is a popular but fruitless exercise, the analysis of generic algorithms for free fields is apparently very informative. The reason that this is so is because our algorithms are sufficiently dumb that they spent most of their time in dealing with the trivial almost free “modes” of the system.

As evidence of the more general applicability of our analysis we point to the excellent agreement with empirical Monte Carlo data of the $`erfc`$ form for the Metropolis acceptance probability as a function of integration step size . Futhermore preliminary results for simple models  indicate that their autocorrelations have at least the same form as expected from our free field theory results.

Perhaps the most surprising result is the failure of the L2MC (Kramers) algorithm to have superior performance to the HMC algorithm. Despite the fact that the noise is put into the Markov process more smoothly the desirable properties of the L2MC algorithm are upset by its “Zitterbewegung” caused by its rare Metropolis rejections.

It is also somewhat unexpected that the HMC algorithm performs nearly as well as the full GHMC algorithm, of which it is a special case, when the parameters of both are chosen optimally. The broad minimum of the costs of these algorithms as a function of acceptance rate has the pleasant consequence that no very fine tuning of parameters is required.

The result that the optimal HMC trajectory length is about the same as the correlation length of the underlying field theory is significant, in that it indicates that the temptation to use shorter, and hence cheaper, trajectories for systems near criticality should be avoided.

The results pertaining to higher order (Campostrini) integrators  are of theoretical interest, but in practice they do not seem to be very useful because the trajectories for interacting field theories are chaotic  and apparently limited by the intrinsic instability of the integrators .

## Acknowledgements

We gratefully acknowledge financial support from PPARC under grant number GR/L22744. This research was also supported in part by the U.S. Department of Energy through Contract Nos. DE–FG05–92ER40742 and DE–FC05–85ER250000.

We would like to thank David Daniel, Robert Edwards, Philippe de Forcrand, Alan Horowitz, Ivan Horváth, Alan Irving, Bálint Joó, Julius Kuti, Steffen Meyer, Hidetoshi Mino, Stephen Pickles, Jim Sexton, Stefan Sint, Alan Sokal, and Zbyszek Sroczynski for useful discussions, comments, and hospitality during the extremely prolonged gestation period of this paper.

## Appendix A Inverse Laplace Transforms
for Autocorrelation Functions

### A.1 Exponentially-Distributed Trajectory Lengths

In this case the Laplace transform of the autocorrelation function, $`F(\beta )=C(t)`$, is a rational function with the numerator of lower degree than the denominator. If the denominator is square free then we have the partial fraction expansion

$$F=\underset{k}{}\frac{A_k}{\beta B_k}.$$

Since

$$e^{\alpha t}_0^{\mathrm{}}𝑑te^{\beta t}e^{\alpha t}=\frac{1}{\alpha +\beta }(Re\alpha +\beta >0),$$

we have

$$^1F=\underset{k}{}A_ke^{B_kt}.$$

The roots $`B_k`$ of the denominator of $`F`$ come in complex conjugate pairs because the coefficients in $`F`$ are real. The autocorrelation function can therefore always be written as a sum of real exponentials (corresponding to the real roots) and of real exponentials multiplied by cosines (corresponding to pairs of complex conjugate roots).

If the denominator of $`F`$ is not square free then the repeated roots give terms of the form $`A/(\beta B)^n`$ in the partial fraction expansion. Consider then

$$t^{n1}e^{\alpha t}_0^{\mathrm{}}𝑑te^{\beta t}t^{n1}e^{\alpha t}=(\alpha +\beta )^n\mathrm{\Gamma }(n)(Re\alpha +\beta >0),$$

This serves to give the inverse Laplace transform in the general case.

### A.2 Fixed Length Trajectories

Here $`F`$ is a rational function of $`\zeta e^{\beta \tau }`$ with the numerator being of lower or equal degree than the denominator. For simplicity we first consider the case where the denominator is square free, whence by partial fractions

$$F=C+\underset{k}{}\frac{A_k}{\zeta B_k}.$$

Observe that

$$\delta (t)_0^{\mathrm{}}e^{\beta t}\delta (t)=1,$$

and, more generally,

$$\underset{j=0}{\overset{\mathrm{}}{}}\mu ^j\delta (t\tau j)=\underset{j=0}{\overset{\mathrm{}}{}}\mu ^je^{\beta \tau j}=\frac{1}{1\mu e^{\beta \tau }}(|\mu |<1),$$

so we have

$`^1F`$ $`=`$ $`C\delta (t){\displaystyle \underset{k}{}}{\displaystyle \frac{A_k}{B_k}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}B_k^j\delta (t\tau j)`$
$`=`$ $`C\delta (t){\displaystyle \underset{k}{}}A_k{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}e^{(j+1)\mathrm{ln}B_k}\delta (t\tau j).`$
All the the roots must satisfy $`|B_k|>1`$ for the geometric series to converge at $`\beta =0(\zeta =1)`$, which is necessary for the integrated autocorrelation function $`F(0)`$ to be finite.

In the general case where the denominator of $`F`$ has multiple roots the general form of the inverse Laplace transform may be obtained from the identity

$$\underset{j=0}{\overset{\mathrm{}}{}}\left(\genfrac{}{}{0pt}{}{n+j1}{j}\right)\mu ^j\delta (t\tau j)=(1\mu \zeta )^n(|\mu \zeta |<1).$$

### A.3 Example of Computation of Autocorrelation Function

To illustrate the computation of autocorrelation functions we shall explicitly evaluate the connected autocorrelation function for $`M^2`$ for the HMC algorithm with exponentially distributed trajectory lengths. The Laplace transform of the autocorrelation function is

$$F(\beta )=\frac{r^3+2r^2\beta +r\beta ^2+2m^2r}{\beta ^3+2r\beta ^2+\left(r^2+4m^2\right)\beta +2m^2r}.$$

If we write this in terms of two distinct roots $`\beta _1`$ and $`\beta _2`$ of the denominator

$$F(\beta )=\frac{r\left(r^2+2r\beta +\beta ^2+2m^2\right)}{\left(\beta \beta _1\right)\left(\beta \beta _2\right)\left(2r+\beta _1+\beta _2+\beta \right)}$$

we can expand this in terms of partial fractions to give

$`F(\beta )`$ $`=`$ $`{\displaystyle \frac{r\left(2r^4+4r^3\beta _1+\left(2\beta _{1}^{}{}_{}{}^{2}m^2\right)r^2+\beta _1rm^2+32m^4+4\beta _{1}^{}{}_{}{}^{2}m^2\right)}{\left(\beta \beta _1\right)\left(2r^413m^2r^2+64m^4\right)}}`$
$`+`$ $`{\displaystyle \frac{r\left(\begin{array}{c}4r^3\beta _1+\left(2\beta _{1}^{}{}_{}{}^{2}13m^22\beta _2\beta _1\right)r^2\\ +\left(7\beta _28\beta _1\right)m^2r+16m^4+\left(4\beta _{1}^{}{}_{}{}^{2}4\beta _2\beta _1\right)m^2\end{array}\right)}{\left(\beta \beta _2\right)\left(2r^413m^2r^2+64m^4\right)}}`$
$`+`$ $`{\displaystyle \frac{r\left(\left(2\beta _2\beta _1+m^2\right)r^2+\left(7\beta _2+7\beta _1\right)m^2r+16m^4+4\beta _2\beta _1m^2\right)}{\left(2r+\beta _1+\beta _2+\beta \right)\left(2r^413m^2r^2+64m^4\right)}}.`$
The inverse Laplace transform of this is immediately obvious

$`{\displaystyle \frac{r\left(2r^4+4r^3\beta _1+\left(2\beta _{1}^{}{}_{}{}^{2}m^2\right)r^2+\beta _1rm^2+32m^4+4\beta _{1}^{}{}_{}{}^{2}m^2\right)e^{\beta _1t}}{2r^413m^2r^2+64m^4}}`$
$`+{\displaystyle \frac{r\left(\begin{array}{c}4r^3\beta _1+\left(2\beta _{1}^{}{}_{}{}^{2}13m^22\beta _2\beta _1\right)r^2\\ +\left(7\beta _28\beta _1\right)m^2r+16m^4+\left(4\beta _{1}^{}{}_{}{}^{2}4\beta _2\beta _1\right)m^2\end{array}\right)e^{\beta _2t}}{2r^413m^2r^2+64m^4}}`$
$`+{\displaystyle \frac{r\left(\begin{array}{c}\left(2\beta _2\beta _1+m^2\right)r^2+\left(7\beta _2+7\beta _1\right)m^2r\\ +16m^4+4\beta _2\beta _1m^2\end{array}\right)e^{\left(2r+\beta _1+\beta _2\right)t}}{2r^413m^2r^2+64m^4}}.`$
In figure 3 we show this function for $`m=10^2`$ and the optimal value $`r_{\text{opt}}=\sqrt{3}m`$, together with the two neighbouring values of $`r=(3\sqrt{3}\pm \sqrt{15})m/2`$ for which the cost is 50% greater.

## Appendix B Evaluation of Spectral Averages

In this section we show how to evaluate the following spectral sums for one-dimensional free field theory:

$`\sigma _A^{(\alpha )}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{sin}^{2\alpha }\left({\displaystyle \frac{\pi p}{V}}\right)`$ (83)
$`=`$ $`{\displaystyle \frac{1}{2^{2\alpha }}}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right)+{\displaystyle \underset{k=1}{\overset{\alpha /V}{}}}()^{kV}\left[\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha kV}}\right)+\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha +kV}}\right)\right]\right\},`$ (84)
$`\sigma _{m^2}^{(\alpha )}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^{2\alpha }`$ (85)
$`=`$ $`{\displaystyle \underset{\beta =0}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha }{\beta }}\right)m^{2(\alpha \beta )}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\beta }{\beta }}\right)+2{\displaystyle \underset{k=1}{\overset{\beta /V}{}}}()^{kV}\left({\displaystyle \genfrac{}{}{0pt}{}{2\beta }{\beta kV}}\right)\right\},`$ (86)
$`\stackrel{~}{\sigma }_{m^2}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{cos}(2\tau \omega _p)`$ (87)
$`=`$ $`{\displaystyle \underset{\beta =0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2m^2)^\beta }{\beta !}}\left({\displaystyle \underset{\gamma =0}{\overset{\beta 1}{}}}\tau _\gamma {\displaystyle ^{\tau _\gamma }}𝑑\tau _{\gamma +1}\right){\displaystyle \underset{k\mathrm{}}{}}\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau ),`$ (88)
$`\overline{\sigma }_{m^2}^{(\alpha )}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^{2\alpha }\mathrm{sin}^2(\tau \omega _p)`$ (89)
$`=`$ $`{\displaystyle \frac{1}{2}}\sigma _{m^2}^{(\alpha )}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{4}}{\displaystyle \frac{d^2}{d\tau ^2}}\right)^\alpha \stackrel{~}{\sigma }_{m^2},`$ (90)
where

$$\omega _p^2m^2+4\mathrm{sin}^2\left(\frac{\pi p}{V}\right).$$
(91)
### B.1 Completeness relations

Let $`𝒢`$ be a group with a family of $`D`$-dimensional representations $`\varphi _q:𝒢(\mathrm{}^D\mathrm{}^D)`$, i.e., $`\varphi _q(g):\mathrm{}^D\mathrm{}^D`$, labelled by some parameter $`q`$; thus $`\varphi _q(gh)=\varphi _q(g)\varphi _q(h)`$ for $`g,h𝒢`$. Summing over group elements (or integrating with respect to Haar measure in the continuous case) gives

$$\sigma _q_𝒢𝑑\mu \varphi _q(g),$$

and multiplying by $`\varphi _q(h)`$ for any $`h𝒢`$ we find

$`\varphi _q(h)\sigma _q`$ $`=`$ $`\varphi _q(h){\displaystyle _𝒢}𝑑\mu \varphi _q(g)={\displaystyle _𝒢}𝑑\mu \varphi _q(h)\varphi _q(g)`$
$`=`$ $`{\displaystyle _𝒢}𝑑\mu \varphi _q(hg)={\displaystyle _𝒢}𝑑\mu \varphi _q(g^{})=\sigma _q`$
where we have made use of the (left) invariance of Haar measure $`d\mu `$. Hence either $`\varphi _q(h)=1h𝒢`$ or $`\sigma _q=0`$. If we let $`q=0`$ label the identity representation of $`𝒢`$, and all other $`\varphi _q`$ be non-identity representations, we have

$$\sigma _q=\delta (q).$$

###### Example B.1.1

Let $`𝒢=\mathrm{}_V`$, $`q\mathrm{}_V`$, and $`\varphi _q(p)=e^{2\pi ipq/V}`$:

$$\frac{1}{V}\underset{p\mathrm{}_V}{}e^{2\pi ipq/V}=\delta _{q,0};$$

###### Example B.1.2

Let $`𝒢=\mathrm{}`$, $`q\mathrm{𝕊}_1=[0,2\pi ]`$, and $`\varphi _q(p)=e^{ipq}`$:

$$\frac{1}{2\pi }\underset{p\mathrm{}}{}e^{ipq}=\delta (q);$$

###### Example B.1.3

Let $`𝒢=\mathrm{𝕊}_1`$, $`q\mathrm{}`$, and $`\varphi _q(p)=e^{ipq}`$:

$$\frac{1}{2\pi }_0^{2\pi }𝑑pe^{ipq}=\delta _{q,0}.$$

### B.2 A Poisson Resummation Formula

Let $`f:\mathrm{𝕊}_1\mathrm{𝕊}_1`$ be a periodic function. As it is periodic it obviously wants to be expanded in eigenfunctions of $`\mathrm{𝕊}_1`$:

$`f(p)`$ $`=`$ $`{\displaystyle _0^{2\pi }}𝑑p^{}\delta (p^{}p)f(p^{})[p\mathrm{𝕊}_1]`$
$`=`$ $`{\displaystyle _0^{2\pi }}𝑑p^{}{\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{q\mathrm{}}{}}e^{iq(p^{}p)}f(p^{})`$
$`=`$ $`{\displaystyle \underset{q\mathrm{}}{}}e^{iqp}{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑p^{}e^{iqp^{}}f(p^{}),`$
hence the “spectral average”

$`\sigma (f)`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}f\left({\displaystyle \frac{2\pi p}{V}}\right)`$
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}{\displaystyle \underset{q\mathrm{}}{}}e^{2\pi ipq/V}{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑p^{}e^{iqp^{}}f(p^{})`$
$`=`$ $`{\displaystyle \underset{q\mathrm{}}{}}\left\{{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}e^{2\pi ipq/V}\right\}\left\{{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑p^{}e^{iqp^{}}f(p^{})\right\}`$
$`=`$ $`{\displaystyle \underset{q\mathrm{}}{}}\left\{{\displaystyle \underset{k\mathrm{}}{}}\delta _{q,kV}\right\}\left\{{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑p^{}e^{iqp^{}}f(p^{})\right\}`$
$`=`$ $`{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑p^{}e^{ip^{}kV}f(p^{}).`$
As $`f`$ is real valued we may take the real part of this identity to obtain

$$\sigma (f)\frac{1}{V}\underset{p\mathrm{}_V}{}f\left(\frac{2\pi p}{V}\right)=\underset{k\mathrm{}}{}\frac{1}{2\pi }_0^{2\pi }𝑑p^{}\mathrm{cos}(p^{}kV)f(p^{}).$$
(92)
### B.3 Multiple Angle Expansion

The following simple identity is often useful

$`\mathrm{sin}^{2\alpha }\theta `$ $`=`$ $`\left[{\displaystyle \frac{1}{2i}}\left(e^{i\theta }e^{i\theta }\right)\right]^{2\alpha }(\alpha \mathrm{})`$ (93)
$`=`$ $`\left({\displaystyle \frac{1}{2i}}\right)^{2\alpha }{\displaystyle \underset{q=0}{\overset{2\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{q}}\right)()^qe^{2i(\alpha q)\theta }`$
$`=`$ $`\left({\displaystyle \frac{1}{2i}}\right)^{2\alpha }\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right)()^\alpha +{\displaystyle \underset{q=0}{\overset{\alpha 1}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{q}}\right)()^q2\mathrm{cos}[2(\alpha q)\theta ]\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2^{2\alpha }}}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right)+{\displaystyle \underset{q^{}=1}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha q^{}}}\right)()^q^{}2\mathrm{cos}(2q^{}\theta )\right\}.`$
### B.4 Free-field Spectral Sums

Let us first consider equation (83). From equation (92) we have

$`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{cos}\left({\displaystyle \frac{2\pi pq^{}}{V}}\right)={\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑p^{}\mathrm{cos}(p^{}kV)\mathrm{cos}(p^{}q^{})`$
$`={\displaystyle \underset{k_in\mathrm{}}{}}{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}{\displaystyle \frac{1}{2}}\left\{\mathrm{cos}[p^{}(kV+q^{})]+\mathrm{cos}[p^{}(kVq^{})]\right\}`$
$`={\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{1}{2}}\left\{\delta _{kV+q^{},0}+\delta _{kVq^{},0}\right\}.`$
Using the multiple-angle expansion (93) we obtain

$`\sigma _A^{(\alpha )}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{sin}^{2\alpha }\left({\displaystyle \frac{\pi p}{V}}\right)`$
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}{\displaystyle \frac{1}{2^{2\alpha }}}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right)+{\displaystyle \underset{q^{}=1}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha q^{}}}\right)()^q^{}2\mathrm{cos}\left({\displaystyle \frac{2\pi pq^{}}{V}}\right)\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2^{2\alpha }}}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right)+{\displaystyle \underset{q^{}=1}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha q^{}}}\right)()^q^{}{\displaystyle \frac{2}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{cos}\left({\displaystyle \frac{2\pi pq^{}}{V}}\right)\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2^{2\alpha }}}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right)+{\displaystyle \underset{q^{}=1}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha q^{}}}\right)()^q^{}{\displaystyle \underset{k\mathrm{}}{}}\left(\delta _{kV+q^{},0}+\delta _{kVq^{},0}\right)\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2^{2\alpha }}}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha }}\right)+{\displaystyle \underset{k=1}{\overset{\alpha /V}{}}}()^{kV}\left[\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha kV}}\right)+\left({\displaystyle \genfrac{}{}{0pt}{}{2\alpha }{\alpha +kV}}\right)\right]\right\}.`$
From this result we may obtain an expression for the quantity $`\sigma _{m^2}^{(\alpha )}`$ defined in equation (85)

$`\sigma _{m^2}^{(\alpha )}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^{2\alpha }={\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\left[m^2+4\mathrm{sin}^2\left({\displaystyle \frac{\pi p}{V}}\right)\right]^\alpha `$
$`=`$ $`{\displaystyle \underset{\beta =0}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha }{\beta }}\right){\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}m^{2(\alpha \beta )}4^\beta \mathrm{sin}^{2\beta }\left({\displaystyle \frac{\pi p}{V}}\right)`$
$`=`$ $`{\displaystyle \underset{\beta =0}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha }{\beta }}\right)m^{2(\alpha \beta )}4^\beta \sigma _A^{(\beta )}`$
$`=`$ $`{\displaystyle \underset{\beta =0}{\overset{\alpha }{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha }{\beta }}\right)m^{2(\alpha \beta )}\left\{\left({\displaystyle \genfrac{}{}{0pt}{}{2\beta }{\beta }}\right)+2{\displaystyle \underset{k=1}{\overset{\beta /V}{}}}()^{kV}\left({\displaystyle \genfrac{}{}{0pt}{}{2\beta }{\beta kV}}\right)\right\}.`$
###### Example B.4.1

$$\sigma _{m^2}^{(2)}=\frac{1}{V}\underset{p\mathrm{}_V}{}\omega _p^4=\left(6+4m^2+m^4\right)+\delta _{V,1}\left(64m^2\right)+\delta _{V,2}2;$$

###### Example B.4.2

$`\sigma _{m^2}^{(3)}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^6=\left(m^6+6m^4+18m^2+20\right)`$
$`+\delta _{V,1}\left(6m^418m^220\right)+\delta _{V,2}\left(6m^2+12\right)\delta _{V,3}2.`$
In order to evaluate Equation (87) we first investigate the simpler case where $`m=0`$

$`\stackrel{~}{\sigma }_0`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{cos}(2\tau \omega _p)`$
$`=`$ $`{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑p^{}\mathrm{cos}(p^{}kV)\mathrm{cos}\left[4\tau \mathrm{sin}\left({\displaystyle \frac{p^{}}{2}}\right)\right]`$
$`=`$ $`{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{1}{\pi }}{\displaystyle _0^\pi }𝑑\theta \mathrm{cos}(2kV\theta )\mathrm{cos}[4\tau \mathrm{sin}\theta ]`$
$`=`$ $`{\displaystyle \underset{k\mathrm{}}{}}\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau ).`$
For small $`m`$ we Taylor expand $`\stackrel{~}{\sigma }_{m^2}`$

$$\stackrel{~}{\sigma }_{m^2}=\underset{\beta =0}{\overset{\mathrm{}}{}}\frac{m^{2\beta }}{\beta !}\frac{d^\beta \stackrel{~}{\sigma }_{m^2}}{d(m^2)^\beta }|_{m^2=0}.$$

Consider the quantity

$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dm^2}}{\displaystyle \frac{d}{d\tau }}{\displaystyle \frac{\stackrel{~}{\sigma }_{m^2}}{\tau }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{d\tau }}\left({\displaystyle \frac{1}{\tau }}{\displaystyle \frac{d}{dm^2}}{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{cos}(2\tau \omega _p)\right)`$
$`=`$ $`{\displaystyle \frac{d}{d\tau }}\left({\displaystyle \frac{1}{\tau }}{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{sin}(2\tau \omega _p)\tau {\displaystyle \frac{d\omega _p}{dm^2}}\right)`$
$`=`$ $`{\displaystyle \frac{d}{d\tau }}\left({\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{sin}(2\tau \omega _p){\displaystyle \frac{1}{2\omega _p}}\right)`$
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\mathrm{cos}(2\tau \omega _p)=\stackrel{~}{\sigma }_{m^2};`$
upon integrating with respect to $`\tau `$ we obtain

$$\frac{d\stackrel{~}{\sigma }_{m^2}}{dm^2}=2\tau ^\tau 𝑑\tau ^{}\stackrel{~}{\sigma }_{m^2},$$

and hence

$$\frac{d^\beta \stackrel{~}{\sigma }_{m^2}}{d(m^2)^\beta }|_{m^2=0}=(2)^\beta \tau ^\tau 𝑑\tau _1\tau _1^{\tau _1}𝑑\tau _2\mathrm{}\tau _{\beta 1}^{\tau _{\beta 1}}𝑑\tau _\beta \stackrel{~}{\sigma }_0.$$

By this means we have obtained the desired result

$$\stackrel{~}{\sigma }_{m^2}=\underset{\beta =0}{\overset{\mathrm{}}{}}\frac{(2m^2)^\beta }{\beta !}\left(\underset{\gamma =0}{\overset{\beta 1}{}}\tau _\gamma ^{\tau _\gamma }𝑑\tau _{\gamma +1}\right)\stackrel{~}{\sigma }_0.$$

Equation (89) may be found by differentiating $`\stackrel{~}{\sigma }_{m^2}`$ with respect to $`\tau `$,

$`\overline{\sigma }_{m^2}^{(\alpha )}`$ $``$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^{2\alpha }\mathrm{sin}^2(\tau \omega _p)`$
$`=`$ $`{\displaystyle \frac{1}{2V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^{2\alpha }\left\{1\mathrm{cos}(2\tau \omega _p)\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2}}\sigma _{m^2}^{(\alpha )}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{4}}{\displaystyle \frac{d^2}{d\tau ^2}}\right)^\alpha \stackrel{~}{\sigma }_{m^2}.`$
###### Example B.4.3

$`\overline{\sigma }_{m^2}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{p\mathrm{}_V}{}}\omega _p^4\mathrm{sin}^2(\tau \omega _p)`$ (94)
$`=`$ $`{\displaystyle \frac{1}{2}}\sigma _{m^2}^{(2)}{\displaystyle \frac{1}{2}}{\displaystyle \underset{\beta =0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{m^{2\beta }}{\beta !}}\left({\displaystyle \frac{1}{4}}{\displaystyle \frac{d^2}{d\tau ^2}}\right)^2{\displaystyle \frac{d^\beta \stackrel{~}{\sigma }_{m^2}}{d(m^2)^\beta }}|_{m^2=0}`$
$`=`$ $`{\displaystyle \frac{1}{2}}\sigma _{m^2}^{(2)}{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \underset{\beta =0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2m^2)^\beta }{2\beta !}}({\displaystyle \frac{1}{4}}{\displaystyle \frac{d^2}{d\tau ^2}})^2\times `$
$`\times \left({\displaystyle \underset{\gamma =0}{\overset{\beta 1}{}}}\tau _\gamma {\displaystyle ^{\tau _\gamma }}𝑑\tau _{\gamma +1}\right)\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau ).`$
The $`\beta =0`$ term in the sum is

$$\frac{1}{2}\underset{k\mathrm{}}{}\left(\frac{1}{4}\frac{d^2}{d\tau ^2}\right)^2\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau )=8\underset{k\mathrm{}}{}\mathrm{J}_{2\mathrm{k}\mathrm{V}}^{}{}_{}{}^{\prime \prime \prime \prime }(4\tau ).$$

From the recurrence relation

$$\mathrm{J}_{\mathrm{n}}^{}{}_{}{}^{}=\frac{1}{2}(\mathrm{J}_{\mathrm{n}1}\mathrm{J}_{\mathrm{n}+1})$$

we find that

$$\frac{d^k\mathrm{J}_\mathrm{n}(z)}{dz^k}=\frac{1}{2^k}\underset{\gamma =0}{\overset{k}{}}\left(\genfrac{}{}{0pt}{}{k}{\gamma }\right)()^\gamma \mathrm{J}_{\mathrm{n}\mathrm{k}+2\gamma }(z),$$

which we may prove by induction:

$`{\displaystyle \frac{d^{k+1}\mathrm{J}_\mathrm{n}(z)}{dz^{k+1}}}={\displaystyle \frac{1}{2^k}}{\displaystyle \underset{\gamma =0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{\gamma }}\right)()^\gamma {\displaystyle \frac{d\mathrm{J}_{\mathrm{n}\mathrm{k}+2\gamma }(z)}{dz}}`$
$`={\displaystyle \frac{1}{2^{k+1}}}{\displaystyle \underset{\gamma =0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{\gamma }}\right)()^\gamma \left(\mathrm{J}_{\mathrm{n}\mathrm{k}+2\gamma 1}(z)\mathrm{J}_{\mathrm{n}\mathrm{k}+2\gamma +1}(z)\right)`$
$`={\displaystyle \frac{1}{2^{k+1}}}{\displaystyle \underset{\gamma =0}{\overset{k}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k}{\gamma }}\right)()^\gamma \left(\mathrm{J}_{\mathrm{n}(\mathrm{k}+1)+2\gamma }(z)\mathrm{J}_{\mathrm{n}(\mathrm{k}+1)+2(\gamma +1)}(z)\right)`$
$`={\displaystyle \frac{1}{2^{k+1}}}\{{\displaystyle \underset{\gamma =1}{\overset{k}{}}}[\left({\displaystyle \genfrac{}{}{0pt}{}{k}{\gamma }}\right)+\left({\displaystyle \genfrac{}{}{0pt}{}{k}{\gamma 1}}\right)]()^\gamma \mathrm{J}_{\mathrm{n}(\mathrm{k}+1)+2\gamma }(z)+`$
$`+\mathrm{J}_{\mathrm{n}(\mathrm{k}=1)}(z)\mathrm{J}_{\mathrm{n}+(\mathrm{k}+1)}(z)\}`$
$`={\displaystyle \frac{1}{2^{k+1}}}{\displaystyle \underset{\gamma =0}{\overset{k+1}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{k+1}{\gamma }}\right)()^\gamma \mathrm{J}_{\mathrm{n}(\mathrm{k}+1)+2\gamma }(z),`$
using $`\left(\genfrac{}{}{0pt}{}{k}{\gamma }\right)+\left(\genfrac{}{}{0pt}{}{k}{\gamma 1}\right)=\left(\genfrac{}{}{0pt}{}{k+1}{\gamma }\right)`$. We therefore find that

$`8{\displaystyle \underset{k\mathrm{}}{}}\mathrm{J}_{2\mathrm{k}\mathrm{V}}^{}{}_{}{}^{\prime \prime \prime \prime }(4\tau )=2{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \underset{\gamma =0}{\overset{4}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{4}{\gamma }}\right)()^\gamma \mathrm{J}_{2\mathrm{k}\mathrm{V}4+2\gamma }(4\tau )`$
$`=2{\displaystyle \underset{k\mathrm{}}{}}\{\mathrm{J}_{2\mathrm{k}\mathrm{V}4}(4\tau )4\mathrm{J}_{2\mathrm{k}\mathrm{V}2}(4\tau )+`$
$`+6\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau )4\mathrm{J}_{2\mathrm{k}\mathrm{V}+2}(4\tau )+\mathrm{J}_{2\mathrm{k}\mathrm{V}+4}(4\tau )\}`$
$`=3\mathrm{J}_0(4\tau )+4\mathrm{J}_2(4\tau )\mathrm{J}_4(4\tau )`$
$`2{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}\{\mathrm{J}_{2\mathrm{k}\mathrm{V}4}(4\tau )4\mathrm{J}_{2\mathrm{k}\mathrm{V}2}(4\tau )+`$
$`+6\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau )4\mathrm{J}_{2\mathrm{k}\mathrm{V}+2}(4\tau )+\mathrm{J}_{2\mathrm{k}\mathrm{V}+4}(4\tau )\}.`$
The $`\beta =1`$ term in the sum (94) is

$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k\mathrm{}}{}}2m^2\left({\displaystyle \frac{1}{4}}{\displaystyle \frac{d^2}{d\tau ^2}}\right)^2\tau {\displaystyle ^\tau }𝑑\tau _1\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau _1)`$
$`=m^2{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{d^4}{d\widehat{\tau }^4}}\widehat{\tau }{\displaystyle ^{\widehat{\tau }}}𝑑\tau ^{}\mathrm{J}_{2\mathrm{k}\mathrm{V}}(\tau ^{})(\mathrm{where}\widehat{\tau }4\tau )`$
$`=m^2{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{d^3}{d\widehat{\tau }^3}}\left\{{\displaystyle ^{\widehat{\tau }}}𝑑\tau ^{}\mathrm{J}_{2\mathrm{k}\mathrm{V}}(\tau ^{})+\widehat{\tau }\mathrm{J}_{2\mathrm{k}\mathrm{V}}(\widehat{\tau })\right\}`$
$`=m^2{\displaystyle \underset{k\mathrm{}}{}}\left\{4\mathrm{J}_{2\mathrm{k}\mathrm{V}}^{}{}_{}{}^{\prime \prime }(\widehat{\tau })+\widehat{\tau }\mathrm{J}_{2\mathrm{k}\mathrm{V}}^{}{}_{}{}^{\prime \prime \prime }(\widehat{\tau })\right\}`$
$`=m^2{\displaystyle \underset{k\mathrm{}}{}}\{{\displaystyle \underset{\gamma =0}{\overset{2}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2}{\gamma }}\right)()^\gamma \mathrm{J}_{2\mathrm{k}\mathrm{V}2+2\gamma }(\widehat{\tau })+`$
$`+{\displaystyle \frac{\widehat{\tau }}{8}}{\displaystyle \underset{\gamma =0}{\overset{3}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{3}{\gamma }}\right)()^\gamma \mathrm{J}_{2\mathrm{k}\mathrm{V}3+2\gamma }(\widehat{\tau })\}`$
$`=m^2{\displaystyle \underset{k\mathrm{}}{}}\{\mathrm{J}_{2\mathrm{k}\mathrm{V}2}(\widehat{\tau })2\mathrm{J}_{2\mathrm{k}\mathrm{V}}(\widehat{\tau })+\mathrm{J}_{2\mathrm{k}\mathrm{V}+2}(\widehat{\tau })+`$
$`+{\displaystyle \frac{\widehat{\tau }}{8}}\mathrm{J}_{2\mathrm{k}\mathrm{V}3}(\widehat{\tau }){\displaystyle \frac{3\widehat{\tau }}{8}}\mathrm{J}_{2\mathrm{k}\mathrm{V}1}(\widehat{\tau })+{\displaystyle \frac{3\widehat{\tau }}{8}}\mathrm{J}_{2\mathrm{k}\mathrm{V}+1}(\widehat{\tau }){\displaystyle \frac{\widehat{\tau }}{8}}\mathrm{J}_{2\mathrm{k}\mathrm{V}+3}(\widehat{\tau })\}`$
$`=m^2\{2\mathrm{J}_0(4\tau )+2\mathrm{J}_2(4\tau )+3\tau \mathrm{J}_1(4\tau )\tau \mathrm{J}_3(4\tau )+`$
$`+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}(2\mathrm{J}_{2\mathrm{k}\mathrm{V}2}(4\tau )4\mathrm{J}_{2\mathrm{k}\mathrm{V}}(4\tau )+2\mathrm{J}_{2\mathrm{k}\mathrm{V}+2}(4\tau )+`$
$`+\tau \mathrm{J}_{2\mathrm{k}\mathrm{V}3}(4\tau )3\tau \mathrm{J}_{2\mathrm{k}\mathrm{V}1}(4\tau )+`$
$`+3\tau \mathrm{J}_{2\mathrm{k}\mathrm{V}+1}(4\tau )\tau \mathrm{J}_{2\mathrm{k}\mathrm{V}+3}(4\tau ))\}.`$
### B.5 Results for Campostrini

$`\overline{\sigma }_{m^2}^{(4)}`$ $`=`$ $`\left[3535\mathrm{J}_0(4\tau )+56\mathrm{J}_2(4\tau )28\mathrm{J}_4(4\tau )+8\mathrm{J}_6(4\tau )\mathrm{J}_8(4\tau )\right]`$
$`+`$ $`\left[4040\mathrm{J}_0(4\tau )+64t\mathrm{J}_1(4\tau )+31\mathrm{J}_2(4\tau )8\mathrm{J}_4(4\tau )+\mathrm{J}_6(4\tau )\right]m^2`$
$`+`$ $`\left[18+(16t^218)\mathrm{J}_0(4\tau )+48t\mathrm{J}_1(4\tau )+8\mathrm{J}_2(4\tau )\mathrm{J}_4(4\tau )\right]m^4+O(m^6),`$
$`\overline{\sigma }_{m^2}^{(6)}`$ $`=`$ $`\left[\begin{array}{c}462462\mathrm{J}_0(4\tau )+792\mathrm{J}_2(4\tau )495\mathrm{J}_4(4\tau )\\ +220\mathrm{J}_6(4\tau )66\mathrm{J}_8(4\tau )+12\mathrm{J}_{10}(4\tau )\mathrm{J}_{12}(4\tau )\end{array}\right]`$
$`+`$ $`\left[\begin{array}{c}756756\mathrm{J}_0(4\tau )+1024t\mathrm{J}_1(4\tau )+698\mathrm{J}_2(4\tau )\\ 256\mathrm{J}_4(4\tau )+69\mathrm{J}_6(4\tau )12\mathrm{J}_8(4\tau )+\mathrm{J}_{10}(4\tau )\end{array}\right]m^2`$
$`+`$ $`\left[\begin{array}{c}525+(256t^2525)\mathrm{J}_0(4\tau )+1152t\mathrm{J}_1(4\tau )\\ +316\mathrm{J}_2(4\tau )74\mathrm{J}_4(4\tau )+12\mathrm{J}_6(4\tau )\mathrm{J}_8(4\tau )\end{array}\right]m^4+O(m^6),`$
$`\overline{\sigma }_{m^2}^{(8)}`$ $`=`$ $`\left[\begin{array}{c}64356435\mathrm{J}_0(4\tau )+11440\mathrm{J}_2(4\tau )8008\mathrm{J}_4(4\tau )\\ +4368\mathrm{J}_6(4\tau )1820\mathrm{J}_8(4\tau )+560\mathrm{J}_{10}(4\tau )\\ 120\mathrm{J}_{12}(4\tau )+16\mathrm{J}_{14}(4\tau )\mathrm{J}_{16}(4\tau )\end{array}\right]`$
$`+`$ $`\left[\begin{array}{c}1372813728\mathrm{J}_0(4\tau )+16384t\mathrm{J}_1(4\tau )\\ +14075\mathrm{J}_2(4\tau )6128\mathrm{J}_4(4\tau )+2185\mathrm{J}_6(4\tau )\\ 608\mathrm{J}_8(4\tau )+123\mathrm{J}_{10}(4\tau )16\mathrm{J}_{12}(4\tau )+\mathrm{J}_{14}(4\tau )\end{array}\right]m^2`$
$`+`$ $`\left[\begin{array}{c}12936+(4096t^212936)\mathrm{J}_0(4\tau )\\ +24576t\mathrm{J}_1(4\tau )+9280\mathrm{J}_2(4\tau )2807\mathrm{J}_4(4\tau )\\ +688\mathrm{J}_6(4\tau )128\mathrm{J}_8(4\tau )+16\mathrm{J}_{10}(4\tau )\mathrm{J}_{12}(4\tau )\end{array}\right]m^4+O(m^6),`$
$`\overline{\sigma }_{m^2}^{(10)}`$ $`=`$ $`\left[\begin{array}{c}9237892378\mathrm{J}_0(4\tau )+167960\mathrm{J}_2(4\tau )\\ 125970\mathrm{J}_4(4\tau )+77520\mathrm{J}_6(4\tau )38760\mathrm{J}_8(4\tau )\\ +15504\mathrm{J}_{10}(4\tau )4845\mathrm{J}_{12}(4\tau )+1140\mathrm{J}_{14}(4\tau )\\ 190\mathrm{J}_{16}(4\tau )+20\mathrm{J}_{18}(4\tau )\mathrm{J}_{20}(4\tau )\end{array}\right]`$
$`+`$ $`\left[\begin{array}{c}243100243100\mathrm{J}_0(4\tau )+262144t\mathrm{J}_1(4\tau )\\ +267814\mathrm{J}_2(4\tau )129872\mathrm{J}_4(4\tau )+54252\mathrm{J}_6(4\tau )\\ 19024\mathrm{J}_8(4\tau )+5420\mathrm{J}_{10}(4\tau )1200\mathrm{J}_{12}(4\tau )\\ +193\mathrm{J}_{14}(4\tau )20\mathrm{J}_{16}(4\tau )+\mathrm{J}_{18}(4\tau )\end{array}\right]m^2`$
$`+`$ $`\left[\begin{array}{c}289575+(65536t^2289575)\mathrm{J}_0(4\tau )\\ +491520t\mathrm{J}_1(4\tau )+233612\mathrm{J}_2(4\tau )82714\mathrm{J}_4(4\tau )\\ +25164\mathrm{J}_6(4\tau )6392\mathrm{J}_8(4\tau )+1300\mathrm{J}_{10}(4\tau )\\ 198\mathrm{J}_{12}(4\tau )+20\mathrm{J}_{14}(4\tau )\mathrm{J}_{16}(4\tau )\end{array}\right]m^4+O(m^6).`$
## Appendix C Spectral sum for one dimensional free field theory

We wish to evaluate the spectral sum

$$\sigma _{m^2}^{(1)}\frac{1}{V}\underset{p\mathrm{}_V}{}\omega _p^2$$

for large $`V`$. Using the Poisson resummation formula we obtain

$$\sigma _{m^2}^{(1)}=\underset{k\mathrm{}}{}\frac{1}{2\pi }_0^{2\pi }𝑑p^{}\mathrm{cos}(p^{}kV)\omega _{Vp^{}/2\pi }^2=\underset{k\mathrm{}}{}\frac{1}{2\pi }_0^{2\pi }𝑑p^{}\frac{\mathrm{cos}(p^{}kV)}{m^2+4\mathrm{sin}^2({\scriptscriptstyle \frac{1}{2}}p^{})}.$$

Upon substituting $`ze^{ip^{}}`$ we get

$`\sigma _{m^2}^{(1)}`$ $`=`$ $`{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{i}{4\pi }}{\displaystyle _{|z|=1}}𝑑z{\displaystyle \frac{z^{kV}+z^{kV}}{z^2(m^2+2)z+1}}`$
$`=`$ $`{\displaystyle \underset{k\mathrm{}}{}}{\displaystyle \frac{i}{2\pi }}{\displaystyle _{|z|=1}}𝑑z{\displaystyle \frac{z^{|k|V}}{(z\mu _+)(z\mu _{})}}={\displaystyle \frac{1}{m\sqrt{m^2+4}}}\left[{\displaystyle \frac{1+\mu _{}^V}{1\mu _{}^V}}\right],`$
where $`\mu _\pm 1+{\scriptscriptstyle \frac{1}{2}}m^2\pm {\scriptscriptstyle \frac{1}{2}}m\sqrt{m^2+4}`$. For $`m1`$ we have $`\mu _{}=1m+O(m^2)`$, and thus

$$\sigma _{m^2}^{(1)}=\frac{\mathrm{coth}{\scriptscriptstyle \frac{1}{2}}mV}{m\sqrt{m^2+4}}$$

which is correct to all orders in the large $`V`$ asymptotic expansion.

## Appendix D Spectral sum for two dimensional free field theory

We want to evaluate the spectral sum

$$\sigma \frac{1}{L_xL_y}\underset{\genfrac{}{}{0pt}{}{p_x\mathrm{}_{L_x}}{p_y\mathrm{}_{L_y}}}{}\omega _p^2$$

where

$$\omega _p^2m^2+4\left(\mathrm{sin}^2\frac{p_x\pi }{L_x}+\mathrm{sin}^2\frac{p_y\pi }{L_y}\right).$$

Using Poisson resummation we find that for $`L_x,L_y\mathrm{}`$

$$\sigma =\frac{1}{(2\pi )^2}_0^{2\pi }𝑑p_x^{}𝑑p_y^{}\left[m^2+4\left(\mathrm{sin}^2\frac{p_x^{}}{2}+\mathrm{sin}^2\frac{p_y^{}}{2}\right)\right]^1,$$

so upon substituting $`ze^{ip_x^{}}`$ we obtain

$`\sigma `$ $`=`$ $`{\displaystyle \frac{i}{4\pi ^2}}{\displaystyle _0^{2\pi }}𝑑p_y^{}{\displaystyle _{|z|=1}}{\displaystyle \frac{dz}{z^2\left(m^2+4\mathrm{sin}^2({\scriptscriptstyle \frac{1}{2}}p_y^{})\right)z+1}}`$
$`=`$ $`{\displaystyle \frac{i}{4\pi ^2}}{\displaystyle _0^{2\pi }}𝑑p_y^{}{\displaystyle _{|z|=1}}{\displaystyle \frac{dz}{(z\mu )(z1/\mu )}}`$
where

$$\mu \frac{m^2}{2}+2\mathrm{sin}^2\frac{p_y^{}}{2}+1\sqrt{\left(\frac{m^2}{2}+2\mathrm{sin}^2\frac{p_y^{}}{2}+1\right)^21}.$$

Evaluating the contour integral gives

$$\sigma =\frac{1}{2\pi }_0^{2\pi }𝑑p_y^{}\left(\mu \frac{1}{\mu }\right)^1=\frac{1}{4\pi }_0^{2\pi }𝑑p_y^{}\left[\left(\frac{m^2}{2}+2\mathrm{sin}^2\frac{p_y^{}}{2}+1\right)^21\right]^{{\scriptscriptstyle \frac{1}{2}}}.$$

We now make the substitution $`z^{}e^{ip_y^{}}`$, which leads to

$`\sigma `$ $`=`$ $`{\displaystyle \frac{i}{2\pi }}{\displaystyle _{|z^{}|=1}}{\displaystyle \frac{dz^{}}{\sqrt{\left[(m^2+4)z^{}z^21\right]^24z^2}}}`$
$`=`$ $`{\displaystyle \frac{i}{2\pi }}{\displaystyle _{|z^{}|=1}}{\displaystyle \frac{dz^{}}{\sqrt{\left[z^2(m^2+2)z^{}+1\right]\left[z^2(m^2+6)z^{}+1\right]}}}`$
$`=`$ $`{\displaystyle \frac{i}{2\pi }}{\displaystyle _{|z^{}|=1}}{\displaystyle \frac{dz^{}}{\sqrt{(z^{}a)(z^{}1/a)(z^{}b)(z^{}1/b)}}},`$
where

$`a`$ $``$ $`{\displaystyle \frac{m^2}{2}}+3\sqrt{\left({\displaystyle \frac{m^2}{2}}+3\right)^21}=(32\sqrt{2}){\displaystyle \frac{1}{8}}(3\sqrt{2}4)m^2+\mathrm{}`$
$`b`$ $``$ $`{\displaystyle \frac{m^2}{2}}+1\sqrt{\left({\displaystyle \frac{m^2}{2}}+1\right)^21}=1m+\frac{1}{2}m^2+\mathrm{}.`$
Since $`1/a>1/b>1>b>a>0`$ the contour integral may be shrunk to be the integral along the branch cut from $`a`$ to $`b`$,

$$\sigma =\frac{1}{\pi }_a^b\frac{dz^{}}{\sqrt{(z^{}a)(bz^{})(1/az^{})(1/bz^{})}}.$$

If we change variables to

$$x\frac{(a+1)(z^{}1)}{(a1)(z^{}+1)}$$

we obtain

$$\sigma =\frac{2\sqrt{ab}}{\pi (a+1)(b1)}_1^{1/k^{}}\frac{dx}{\sqrt{(x^21)(1k^2x^2)}}$$

where $`k^{}\frac{(a1)(b+1)}{(a+1)(b1)}`$. With the substitution $`\xi \frac{1}{k}\sqrt{1\frac{1}{x^2}}`$ we find that

$$_1^{1/k^{}}\frac{dx}{\sqrt{(x^21)(1k^2x^2)}}=_0^1\frac{d\xi }{\sqrt{(1\xi ^2)(1k^2\xi ^2)}}=K(k)$$

where $`K`$ is the complete elliptic integral and $`k^2+k^2=1`$; thus

$$\sigma =\frac{2\sqrt{ab}}{\pi (a+1)(b1)}K(k).$$

This may be further simplified using Landen’s transformation

$$K\left(\frac{2\sqrt{\eta }}{1+\eta }\right)=(1+\eta )K(\eta )$$

with $`\eta \frac{ab}{1ab}`$, which leads to

$$\sigma =\frac{2\sqrt{ab}}{\pi (1ab)}K\left(\frac{ba}{1ab}\right),$$

using the fact that $`K`$ is an even function. Expanding in powers of $`m`$ we obtain

$$\sigma =\frac{1}{\pi }\left(1+O(m)\right)K\left(1\sqrt{2}m+O(m^2)\right)=\frac{1}{4\pi }\mathrm{ln}\left(\frac{32}{m^2}\right)+O(m\mathrm{ln}m),$$

because

$$K(k)=K(\sqrt{1k^2})\mathrm{ln}\frac{4}{k^{}}\text{as }k^{}0.$$