We present a straightforward implementation of contour integration by setting options for and , taking advantage of powerful results in complex analysis. As such, this article can be viewed as documentation to perform numerical contour integration with the existing built-in tools. We provide examples of how this method can be used when integrating analytically and numerically some commonly used distributions, such as Wightman functions in quantum field theory. We also provide an approximating technique when time-ordering is involved, a commonly encountered scenario in quantum field theory for computing second-order terms in Dyson series expansion and Feynman propagators. We believe our implementation will be useful for more general calculations involving advanced or retarded Green’s functions, propagators, kernels and so on.

Introduction

It is well known that we can integrate analytically a large class of functions with known anti-derivatives via ; otherwise, we can use for numerical results. There are various settings that one can use to evaluate integrals, depending on the task at hand. Crucially, this can be performed even if the integrand is complex valued, such as the function . In some cases, when the integrand is a distribution, this integration can also be done analytically by setting the option , for example, or by imposing a regulator, which in physics is often called a UV cutoff.

In this article, we are interested in a contour integral of the form

(1)

where is a complex-valued function or distribution and is an integration contour, with possibly higher-dimensional generalizations. For a large class of functions and a contour , this can be done in various ways, such as using powerful techniques in complex analysis like Cauchy’s integral formula or the residue theorem. For many of these functions, can provide answers immediately, sometimes involving special functions.

The motivation for this article is based on the observation that there is no good documentation on how to perform explicit contour integration in Mathematica. It turns out that the ingredients necessary to perform this task are already available within the software. We believe that some of these ingredients are already used elsewhere for different purposes. Our task is to synthesize these components in a coherent way to document how to perform numerical contour integration, using only built-in commands and their options. We will show by way of examples that in many cases this method proves superior to direct integration with prescription, especially when the integrand is complicated or when the integral is multidimensional.

Example 1: Simple Textbook Example—One-Dimensional Contour Integral with Closed-Form Solutions

Consider the integral

(2)

where , , and is the entire real line. The first equality comes from Sokhotsky’s formula [1] and the second equality uses principal value integration. This is a commonly used “epsilon-regularization” in physics when evaluating certain apparently divergent integrals, which has to do with the distributional nature of the integrand in question (see Appendix A in [1] for examples).

This integral can be solved using complex analysis instead of Sokhotsky’s formula, based on the observation that can be computed using the following contour.

The contour is the horizontal blue segment and the semi-circular arc . The integral over the semi-circular arc vanishes as by Jordan’s lemma [2, 3]. The contribution from the closed loop can be calculated in various ways, such as the residue theorem. The residue theorem says that since the closed loop encloses the pole at (using the standard counterclockwise convention [2, 3]), . Taking the limit as , we get the same result as before.

There is indeed a built-in function for calculating residues.

Another way to evaluate is based on the fact that the “ prescription” is an instruction to perform a specific contour integration (without an -regulator) in the sense that

(3)

where is a contour chosen to go from to but deformed near the pole to the lower complex plane.

The prescription in (2) shifts the pole at to in the upper half-plane, so that we can integrate over the real line and then take the limit as . Consequently, if we were to remove the -regulator, the deformation theorem [2] would require that we deform the contour to the lower complex plane near in order to obtain the same value of the integral. (The deformation theorem [2] in complex analysis states that a contour integral remains constant under deformation of the contour that does not cross any poles or branch cuts.) For this contour, the residue is instead given by the following.

If we choose instead, the integral is zero. Here is an explicit computation for fixed .

Case 1

The zero values have a simple explanation in terms of contour integration: since the pole is now shifted to the lower half-plane, the contour does not enclose any poles, which by Cauchy’s integral theorem guarantees that the integral is zero.

Case 2

This is a good place to introduce the numerical contour integration scheme. As a simple example, we can use a rectangular contour.

Case 3

For any and , this calculation gives the correct answer, as expected from direct calculation using residues or textbook pen-and-paper calculations (see e.g. [3]). Varying does not change the integral, as a consequence of the deformation theorem. In particular, other contours also work, such as, for example, a triangular contour.

Case 4

For completeness, let us show that in this simple case the prescription (which is commonly used in physics) can also be done. We add the condition .

The final result involves , which tends to 1 as . Therefore, in this simple example, all the methods yield good results. It is straightforward to show that many common textbook examples (see e.g. [3]) can be performed in this way.

Case 5

Finally, our main goal is to show that can also perform the task reliably. First do the prescription.

Here we see the first instance where the prescription starts to fail: when performed numerically, the most basic setting does not work, and increasing the working precision does not help much.

Remark: The command silences all error messages and is used for readability. In practice, in order to discover the issues and possibly infer the sources of problems, it is useful to not use on a first try.

Case 6: Numerical Contour Integration

Now do the contour integral using a triangular contour.

The contour integration works as is without further modification. The size of the contour does not matter, as long as it is not so small that it leads to numerical instabilities associated with small numbers.

We refer to the numerical computation using that involves contour deformation in the complex plane as numerical contour integration. We will see recurring examples where prescriptions, albeit simplest to implement, almost never work when needed in spite of optimizing various settings. Thus numerical contour integration is a highly competitive technique to use.

Example 2: Two-Dimensional Integral on over a Distribution

(Remark: Some of the discussion here was briefly outlined in [4, 5], which is based on what this article contains.)

The main example we consider in this article involves the integral over a distribution or bi-distribution of the form

(4)

where , , and is a distribution or bi-distribution given by

(5)

In physics, more specifically quantum field theory, this distribution is the two-point vacuum correlation function (also called the two-point Wightman function) of a massless scalar field in ()-dimensional Minkowski space [6].

(6)

The arguments of the two-point function are the coordinates of two spacetime events and , where and in Cartesian coordinates and . The integral is evaluated at and is thus independent of the spatial coordinates. The prescription tells us that the Wightman function is a distribution or bi-distribution [6].

Our task is to evaluate the integral in (4). This is a two-dimensional complex integral over and with a continuum of poles at for any fixed . In the literature, gives the transition probability (divided by a small coupling constant ) of an Unruh–DeWitt (UDW) detector consisting of a two-level quantum system (qubit) with energy gap interacting with the massless scalar field for a finite duration prescribed by the Gaussian function with timescale set by (see e.g. [7] and references therein). The Gaussian functions in the integrand ensure that the interaction between the detector and the quantum field is adiabatically switched on and off smoothly; that is, the Fourier transform of the Wightman function has polynomial tails and is strongly suppressed by the exponential tails of the Fourier transform of the Gaussians. This guarantees that there is no spurious divergence in the transition probability due to switching on the detector very sharply or suddenly, such as with a switching function with discontinuity at an endpoint, like a rectangular function, which is unphysical. An analogy is switching on a lightbulb: in realistic scenarios the current in the circuit does not simply increase from zero to a constant value instantaneously, but rather smoothly over some short time interval.

The example in equation (4) is important for three reasons.

1. It is one of the simplest examples in quantum field theory; many other examples in physics involve distributions that take similar and sometimes more complicated forms (e.g. different power laws where or for some function ).

2. It is one of the simplest examples of a two-dimensional integral with a continuum of poles, thus going beyond standard textbook examples (e.g. in [3] all the examples are one dimensional).

3. Despite its somewhat complicated appearance, it admits a closed-form expression [7].

Benchmark Value for

In [7] it is shown that the exact closed-form expression is given by

(7)

where is the complementary error function. We want to show that numerically this result can be obtained in a satisfactory way. For this purpose, let us focus on a particular numerical value by setting , so that

(8)

We can then use this value as a benchmark for our methods.

Direct Integration

Direct integration for is given by the following command.

The most modest settings yield poor results.

Clearly, although the real parts are all positive (hence physically at least sensible for a transition probability), the results are very different from the benchmark value. The third value is wrong because the probabilities should be real (hence numerically the imaginary part should be much closer to zero).

Increasing as a first remedy does not help much, if at all.

Restricting the Domain

Restricting the integration domain to (, ), where the Gaussian is strongly supported, improves the result, but not enough. (The negligible part of the integral associated to Gaussian tails was ignored.) To see that the behavior scales badly with , let us choose a particular setting.

The first two entries, while not accurate enough, are again at least physically valid results because the imaginary part is negligible and hence can represent transition probabilities. In contrast, the last three values have negative real parts and are thus invalid by default.

We leave it to the reader to verify that other possible settings, such as increasing or , do not help, nor does changing the setting from to another. There may be some combination of settings that makes this integration work, but if they exist then they are for our purposes not worth the time. We will in fact show that the exact same setting that yields the bad results above does work with numerical contour integration.

A possible reason for these issues is that prescription changes the asymptotic behavior at . Numerically, we are evaluating the integral from to , so has to deal with points at infinity that are shifted with a very small imaginary component.

Numerical Contour Integration

Here is transformed into a numerical contour integral, . In contrast to the previous subsection, the -regulator now serves as a contour deformation near a reasonably shaped contour. Hence the asymptotic behavior of the integral is not altered by varying .

(9)

where is a contour along the real line but deformed around the poles of the Wightman function. Here we have in the Wightman function because we have converted the prescription into an instruction to perform contour integration, so the is now part of the definition of the contour .

Here is the contour , where is the location of the continuum of poles for every fixed .

Since many contours work, let us first choose a particularly simple one: a rectangular deformation near the poles. If we integrate over first, then the continuum of poles is located in the complex plane at . Therefore, the contour should be deformed to the upper complex plane.

Based on the previous subsection, we see that using yields a better answer, since more than 99.99% of the area of the Gaussian in the integrand is contained in the interval . We will not set in what follows. Consequently, the parts of the contour from to and to have been truncated and do not yield significant error. This gives the following results.

This result may look undesirable, but the first value is numerically the same as the benchmark value up to very small imaginary part. This suggests that in the neighborhood of there is a good range that we can use to get a numerically accurate result; we will show that the result is invariant under smooth deformation of the contour. For example, if we zoom in between and , we get the following results.

That the results are unchanged as we vary (only the negligible error changes a little) is what we expect from the deformation theorem in complex analysis, since the deformation does not cross any singularities. This invariance under contour deformation is very useful because it is one of the consistency checks we can do (or perhaps need to do) in the absence of exact closed-form expressions and other cross-validation. Direct integration does not allow this, because every yields different values (because the asymptotics are also shifted by ). Furthermore, numerical contour integration is also remarkably fast.

A Different Contour

Here is a different contour. The width is fixed but the height toward the imaginary axis varies.

This choice of contour highlights that while contour integration is much more powerful than the direct integration method, one should avoid numerical issues associated with values of that are either too large or too small; sometimes different contours yield different stability. Nonetheless, the main takeaway is that we have found that there is a large enough range of for which the contour can be varied but still yields the correct value (in this case, ). In this range, we can vary and the answer is invariant. The correct range of that provides reliable results must however be found empirically and depends strongly on the integrand. For completeness, we show the results for zooming into the range of (except that this time we need this to be ).

A significant advantage of numerical contour integration over direct integration is that it allows for robust consistency check: if we can find a range for in which the integral is constant (up to small numerical errors like a small imaginary part), then the integral is likely to be correct. This is an extremely important necessary condition, especially when a closed-form expression is not available. Of course, if there are other ways to check the correctness of the value (as in this example where we know the exact closed-form expression or if there are physical arguments to back it up), we should always also perform these other checks.

Example 3: Two-Dimensional Integral over a Distribution Containing the Heaviside Step Function

Next, we consider a more complicated integral,

(10)

where , , , is the Heaviside step function and the two-point Wightman function is given in equation (5). This integral differs conceptually from Example 2 in two respects:

(1) the additional Heaviside step function in the integrand

(2) we also allow for and changed the sign on the phase

The Heaviside step function appearing in the integrand naturally arises in physics calculations where a Dyson series expansion is involved (i.e. where the notion of time ordering is necessary). In particular, this occurs in various time-dependent perturbation theory calculations within quantum mechanics and quantum field theory. The specific integral in (10) appears as part of the calculation of the nonlocal part of the joint two-detector density matrix of two qubits interacting with a massless scalar field in the so-called entanglement harvesting protocol [7, 8].

Perhaps somewhat surprisingly, this integral admits a closed-form solution [7], given by

(11)

where . As before, we can assume that the coordinate system is aligned so that .

Benchmark Value for

For benchmarking, assume and set . This gives

(12)

Here is the calculation.

This is an alternative if one prefers arbitrary in the expression.

This number is complex, so we no longer have the benefit of verifying our numerical calculation by requiring that the imaginary part be small (which we used in Example 2).

Let us now evaluate numerically. We will not bother with the prescription anymore because it does not work; this is not surprising since the Heaviside step function only complicates the integral compared to Example 2. However, we will evaluate this integral in two ways, both using numerical contour integration, with the two methods differing in how they handle the Heaviside step function.

Absorbing the Heaviside Function into the Upper Limit

We rewrite the integral as:

(13)

By suitable relabeling, the two integrals are equal (this property is used in [7, 8]), but we will not need this fact, since in more general situations the two integrals are different. (See e.g. [5, 9], where the upper limits differ slightly by a “relative redshift” factor.) Although equation (13) is written in the notation of an prescription, this is merely a shorthand for performing the numerical contour integration, since we already know that direct integration does not work.

Pick a rectangular contour and set it so that . For this choice, the two integrals have a continuum of poles at and , respectively. The contour will be chosen so that it does not cross any poles as varies.

We can then write the following integral command.

The most modest settings yield good results.

The contour parameter is good in the range and starts to deviate at around . Remarkably, we can compute this result very well despite the extra complications, with nearly no additional settings beyond .

Approximate Heaviside Step Function

The integral in equation (9) is actually written in a more useful form than the upper limit form in equation (12), in the sense that it accommodates a more general scenario where the two-point functions may not be simple functions of and . In general, the poles of the Wightman functions for the problem at hand can be located at , where is some function that does not have a simple inverse. In such a situation, the upper limit of the integrals is not simply or and may have to be worked out numerically.

More specifically, in the context of quantum field theory in curved spacetimes, the Wightman function in equation (5) is given in flat spacetime in terms of the Minkowski coordinates, and the integrals we computed thus far correspond to two Unruh–DeWitt (UDW) detectors separated by proper distance , both of which are at rest relative to this coordinate system. When the detectors are in motion, we have to replace with , where is Minkowski coordinate time and and are the proper times of each detector, which take the roles of and in the previous sections. Similarly, we have to replace with , where now takes the role of Minkowski spatial coordinates. In this case, in general it is not true that is a simple function of , so it may not always be possible to find an upper limit version similar to equation (12). The static detectors (at rest relative to Minkowski coordinates) correspond to the special case where , which can then be rewritten as . This complexity occurs more generally when the spacetime is not flat, such as when we consider two-dimensional truncation of Schwarzschild spacetime [5].

For this reason, it is useful to try to integrate equation (10) using the Heaviside step function directly. Since the built-in function is not defined for complex numbers, we approximate it with a smooth function. For example, the Heaviside step function is a limit of a logistic function:

(14)

provided that we define . We can thus define an approximate Heaviside step function for any fixed by

(15)

Any other function used in computing a cumulative probability distribution would work.

Our integral now becomes

(16)

again with the rectangular contour of Figure 2. The task now is to manually find and that yield good results by trial and error. Again we set so that .

We define and .

The results are promising, since there is a range of that gives a result very close to the benchmark , but we can do better.

The key is to set so that . For example, for , a good choice of is 1/5 but not 1 or 1/12. The only difference in the next command from the last is in the fourth argument of , where is instead of 5. The results are slightly better.

We can study this relationship better for a smaller set of values of and judge the quality.

Thus we have obtained a good range of where the deformation theorem clearly works. We need to work with because the smooth approximation of the Heaviside step functions introduces infinitely many new poles in the complex plane, since has (countably) infinitely many poles along the imaginary axis. This can be easily seen by rewriting as:

(17)

Therefore the choice of must depend on so that the contour does not cross extra poles, making the integral gain spurious contributions. Nonetheless, finding the right that works is an empirical process, especially when the Wightman function is very complicated. Our examples thus far are simple enough that they can be easily cross-validated and evaluated without much effort.

Further Applications

The above examples have been used by the author to calculate the density matrix elements of Unruh–DeWitt detectors interacting with massless scalar fields for various background spacetimes, where the calculations require time-dependent perturbation theory involving Dyson series expansions. Starting from the calculations in [4] where a primitive version was first used to study entanglement harvesting in two-dimensional spacetimes with accelerating mirror (see the Appendix in [4]), the technique was subsequently improved (using a full rectangular contour that is not close to the poles) for the entanglement harvesting problem in Schwarzschild and collapsing shell spacetimes [5]. In [5], the derivative-coupling functions (strictly speaking they are not Wightman functions but the proper time derivatives of the Wightman functions) are so complicated that it is somewhat surprising the upper-limit technique outlined above worked very well across large parameter spaces. The author has verified that the approximation of the Heaviside step function works well for calculations in [5], although with the caveat: when the derivative-coupling Wightman functions can be expressed as a sum of distinct terms, the part of the derivative Wightman function that contains and requires relatively small but the part that contains or may need very large . This highlights the need for empirically testing the choice of before proceeding with physical calculations.

We extend the results of [5] for harvesting of correlations by two quantum-mechanical detectors interacting with a quantum field in [9], where we consider more complicated detector trajectories. In this context, the authors verified that the upper-limit form analogous to equation (12) proved numerically unstable in physically relevant regimes, thus making indispensable the approximate Heaviside step functions as a way to calculate the nonlocal correlation term in entanglement harvesting contexts. In that setting, the calculation of the nonlocal term involves an integral similar to in equation (10, but with much more complicated numerical calculations than in [5]; in particular the upper limits are not simply affine functions of or . We will explain the technicalities specific to the problem at hand when invoking the techniques outlined in [9].

In quantum field theory in generic spacetimes (flat or curved), there are many situations of computing integrals over distributions. Our techniques should be general enough to use in such contexts, where correlation functions take a form similar to Wightman functions. For example, we have verified that our techniques work for computing field commutators [5] and in ()-dimensional rotating BTZ black holes [10, 11]. We are also able to reproduce the computations of transition rates as done for example in [12] (also see Appendix of [4]), which are fundamentally simpler because they are one-dimensional integrals over the Wightman functions.

Conclusion

In this article we have presented a straightforward implementation of contour integration using the and commands, taking advantage of powerful results in complex analysis. The goal is to provide a systematic example-based implementation of contour integration that does not require any user-defined functions and only requires setting options. We provide explicit examples of how this can be used when integrating analytically and numerically some commonly used distributions, such as Wightman functions in quantum field theory. We also provide efficient approximation schemes to compute quantities involving time-ordering operations where Heaviside step functions are required. Our examples are geared toward research in relativistic quantum information, where finite-dimensional quantum systems interact with a quantum field on a possible curved background spacetime. We believe our implementation will be useful for calculations involving various Green’s functions, propagators, integral kernels and so on for both students and researchers in fields requiring numerical evaluations of integrals over distributions without resorting to unnecessary or excessive approximation schemes.

Acknowledgments

The author thanks Kensuke Gallock-Yoshimura for help in testing the numerical schemes proposed here in [9], and Nicholas Funai for reading through the draft and providing useful feedback. The author acknowledges support from the Mike-Ophelia Lazaridis Fellowship from the Institute for Quantum Computing (IQC), University of Waterloo, Canada. The author also would like to thank both the reviewer and editor for very useful and clarifying comments and recommendations.

References

[1]	J. M. Howie, Complex Analysis, London: Springer Science & Business Media, 2003.
[2]	J. W. Brown and R. V. Churchill, Complex Variables and Applications, 8th ed., New York: McGraw-Hill Education, 2009.
[3]	V. Mukhanov and S. Winitzki, Introduction to Quantum Effects in Gravity, New York: Cambridge University Press, 2007.
[4]	E. Tjoa. “Aspects of Quantum Field Theory with Boundary Conditions,” UWSpace, 2019. hdl.handle.net/10012/14843.
[5]	E. Tjoa and R. B. Mann, “Harvesting Correlations in Schwarzschild and Collapsing Shell Spacetimes,” Journal of High Energy Physics, 08(155), 2020. doi.org:10.1007/JHEP08(2020)155.
[6]	N. D. Birrell and P. C. W. Davies, Quantum Fields in Curved Space, Cambridge: Cambridge University Press, 1982.
[7]	A. R. H. Smith. “Detectors, Reference Frames, and Time,” UWSpace, 2017. hdl.handle.net/10012/12618.
[8]	G. Ver Steeg and N. C. Menicucci, “Entangling Power of an Expanding Universe,” Physical Review D, 79(4), 2009 044027. doi.org:10.1103/PhysRevD.79.044027.
[9]	K. Gallock-Yoshimura, E. Tjoa and R. B. Mann, “Harvesting Entanglement with Detectors Freely Falling into a Black Hole.” arxiv.org/abs/2102.09573.
[10]	M. P. G. Robbins, L. J. Henderson and R. B. Mann, “Entanglement Amplification from Rotating Black Holes.” arxiv.org/abs/2010.14517.
[11]	L. J. Henderson, R. A. Hennigar, R. B. Mann, A. R. H. Smith and J. Zhang, “Harvesting Entanglement from the Black Hole Vacuum,” Classical and Quantum Gravity, 35(21) 21LT02, 2018. doi.org/10.1088/1361-6382/aae27e.
[12]	J. Louko and A. Satz, “Transition Rate of the Unruh–DeWitt Detector in Curved Spacetime,” Classical and Quantum Gravity, 25(5) 055012, 2008. doi.org/10.1088/0264-9381/25/5/055012.
Erickson Tjoa, “Numerical Contour Integration,” The Mathematica Journal, 2021. doi.org/10.3888/tmj.23–3.

About the Author

Erickson Tjoa was born in Jakarta, Indonesia. He moved to Singapore at the age of 15 and studied physics and mathematics as a double major at Nanyang Technological University from 2013-2017. He then completed an MSc at the University of Waterloo, where he is currently a second-year PhD student in relativistic quantum information.

Erickson Tjoa
Department of Physics and Astronomy
Institute for Quantum Computing (IQC)
University of Waterloo, ON, Canada
200 University Avenue West
Waterloo, ON, Canada N2L 3G1
e2tjoa@uwaterloo.ca