# Test Problems for Potential Field Source Surface Extrapolations of Solar and
Stellar Magnetic Fields

David Stansby Mullard Space Science Laboratory, University College London,
Holmbury St. Mary, Surrey RH5 6NT, UK  Daniel Verscharen Mullard Space Science
Laboratory, University College London, Holmbury St. Mary, Surrey RH5 6NT, UK

###### Abstract

The potential field source surface (PFSS) equations are commonly used to model
the coronal magnetic field of the Sun and other stars. As with any
computational model, solving equations using a numerical scheme introduces
errors due to discretisation. We present a set of tests for quantifying these
errors by taking advantage of analytic solutions to the PFSS equations when
the input field is proportional to a single spherical harmonic. From the
spherical harmonic solutions we derive analytic equations for magnetic field
lines traced through the three dimensional magnetic field solution. We propose
these as a set of standard analytic solutions that all PFSS solvers should be
tested against to quantify their inherent errors. We apply these tests to the
pfsspy software package, showing that it reproduces spherical harmonic
solutions well with a slight overestimation of the unsigned open magnetic
flux. It is also successful at reproducing analytic field line equations, with
errors in field line footpoints typically much less than one degree.

##  1 Introduction

The potential field source surface (PFSS) set of equations (Altschuler &
Newkirk, 1969; Schatten et al., 1969) are commonly used to model the coronal
magnetic field of the Sun (e.g. Badman et al., 2020; Stansby et al., 2021;
Fargette et al., 2021) and other stars (e.g. Jardine et al., 2017; Saikia et
al., 2020; Kochukhov, 2020). The key assumption of PFSS models is the absence
of electric current within the domain of interest. Under this assumption the
electromagnetic equations for the magnetic field ($\mathbf{B}$) reduce to

$\displaystyle\nabla\cdot\mathbf{B}$ $\displaystyle=0$ (1)
$\displaystyle\nabla\times\mathbf{B}$ $\displaystyle=\mathbf{0}$ (2)

These equations are equivalent to the Laplace equation,

$\nabla^{2}\Phi=0,$ (3)

where the scalar potential is defined as $\mathbf{B}=\nabla\Phi$. The solution
is derived between the stellar surface (at a radius r⊙) and a specified source
surface (at a radius $r_{ss}>r_{\odot}$), and three boundary conditions are
required for a unique solution. One of these is given by a single component of
the magnetic field at the inner boundary ($r=r_{\odot}$), and the other two
are set by the assumption that the field is purely radial (ie. the two
transverse components are zero) on the outer boundary ($r=r_{ss}$). The second
assumption is motivated by accelerating stellar wind plasma forcing the field
to be near-radial at the source surface.

Different numerical methods are available for solving the above equations,
including spherical harmonic expansions (e.g. Altschuler & Newkirk, 1969;
Hakamada, 1995; Tóth et al., 2011) and finite difference methods (e.g. Jiang &
Feng, 2012; Tóth et al., 2011; Caplan et al., 2021). To understand the
topology of the magnetic field (e.g. van Driel-Gesztelyi et al., 2012; Boe et
al., 2020; Baker et al., 2021) and how stellar wind flows through a corona
(e.g. Neugebauer et al., 1998; Stansby et al., 2020a), magnetic field line
connectivities between the outer and inner boundaries are also needed, which
requires tracing magnetic field lines through the three-dimensional magnetic
field solution using a field line tracer.

In order to check whether numerical methods work as expected, and to quantify
any errors inherent in the numerical scheme employed, it is helpful to compare
their output with exact analytical solutions. While in some areas such as
hydrodynamic modelling test problems are common (Sod, 1978), to our knowledge
no such test problems have been published for PFSS solvers and field line
tracers.

In this paper we provide a set of analytical solutions for the PFSS equations
(Section 2.1) and use these to derive solutions for the unsigned open flux at
the source surface (Section 2.2) and field lines traced through these
solutions (Section 2.3). We then compare these to the numerical solutions
computed by the pfsspy solver (Stansby et al., 2020b) (Section 4), to
demonstrate its usefulness, accuracy, and limitations.

##  2 Analytical solutions

In this section we recount analytical solutions to the PFSS equations using
spherical harmonics (Section 2.1). Using these solutions we then derive
equations for the total open magnetic flux (Section 2.2) and magnetic field
lines traced through the spherical harmonic solutions (Section 2.3).

###  2.1 PFSS solutions

Following Wang & Sheeley (1992); Mackay & Yeates (2012), when the radial
component of the magnetic field on the stellar surface ($r_{\odot}$) is
specified the general solution to the PFSS equations in a spherical coordinate
system is given by the spherical harmonic decomposition

$\displaystyle B_{r}$
$\displaystyle=\sum_{\ell=1}^{\infty}\sum_{m=-\ell}^{m=\ell}a_{\ell
m}c_{\ell}\left(r\right)Y_{\ell m}\left(\theta,\phi\right)$ (4) $\displaystyle
B_{\theta}$ $\displaystyle=\sum_{l=1}^{\infty}\sum_{m=-\ell}^{m=\ell}a_{\ell
m}d_{\ell}\left(r\right)\frac{\partial Y_{\ell
m}\left(\theta,\phi\right)}{\partial\theta}$ (5) $\displaystyle B_{\phi}$
$\displaystyle=\sum_{\ell=1}^{\infty}\sum_{m=-\ell}^{m=\ell}a_{\ell
m}d_{\ell}\left(r\right)\frac{1}{\sin\theta}\frac{\partial Y_{\ell
m}\left(\theta,\phi\right)}{\partial\phi}$ (6)

where

$Y_{\ell m}\left(\theta,\phi\right)=A_{\ell m}P_{\ell
m}\left(\cos\theta\right)\begin{Bmatrix}\cos m\phi\\\ \sin
m\phi\end{Bmatrix}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\begin{matrix}m>0\\\ m\leq
0\end{matrix}$ (7)

are real spherical harmonics, $A_{\ell m}$ are normalisation coefficients and
$P_{\ell m}$ are the associated Legendre polynomials.

$a_{\ell m}$ are constant coefficients derived from the input magnetic field
via.

$a_{\ell m}=\int_{0}^{2\pi}d\phi\int_{0}^{\pi}\sin\theta
d\theta~{}B_{r\odot}\left(\theta,\phi\right)Y_{\ell
m}\left(\theta,\phi\right)$ (8)

and the functions $c_{\ell}$ and $d_{\ell}$ are

$\displaystyle c_{\ell}\left(z\right)$
$\displaystyle=\left(\frac{r}{r_{\odot}}\right)^{-\ell-2}\left[\frac{\ell+1+\ell\left(r/r_{ss}\right)^{2\ell+1}}{\ell+1+\ell(r_{\odot}/r_{ss})^{2\ell+1}}\right]$
(9) $\displaystyle d_{\ell}\left(z\right)$
$\displaystyle=\left(\frac{r}{r_{\odot}}\right)^{-\ell-2}\left[\frac{1-\left(r/r_{ss}\right)^{2\ell+1}}{\ell+1+\ell(r_{\odot}/r_{ss})^{2\ell+1}}\right]$
(10)

The polar coordinate range is $[0,\pi]$, with $\theta=0$ the north pole.

If the input magnetic field is directly proportional to a single spherical
harmonic, $B_{r}=B_{0}Y_{\ell^{\prime}m^{\prime}}$ where $B_{0}$ is a
constant, and the $a_{\ell m}$ coefficients simplify greatly to

$a_{\ell m}=B_{0}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}$ (11)

and only a single term in each sum is non-zero. For the remainder of this
paper we set the input magnetic field proportional to a single spherical
harmonic, so for brevity drop the apostrophes that denote a specific choice of
$\ell,m$. Under these assumptions the solutions are

$\displaystyle B_{r}$ $\displaystyle=B_{0}c_{\ell}\left(r\right)Y_{\ell
m}\left(\theta,\phi\right)$ (12) $\displaystyle B_{\theta}$
$\displaystyle=B_{0}d_{\ell}\left(r\right)\frac{\partial Y_{\ell
m}\left(\theta,\phi\right)}{\partial\theta}$ (13) $\displaystyle B_{\phi}$
$\displaystyle=B_{0}d_{\ell}\left(r\right)\frac{1}{\sin\theta}\frac{\partial
Y_{\ell m}\left(\theta,\phi\right)}{\partial\phi}$ (14)

In the next sections we provide equations for calculating the total unsigned
open flux (Section 2.2) and magnetic field line equations (Section 2.3) from
these solutions.

###  2.2 Open flux

The total unsigned open flux is defined by integrating the radial component of
the magnetic field on the source surface

$\Phi_{open}=r^{2}_{ss}\int_{0}^{2\pi}\int_{0}^{\pi}\left|B_{r}\left(\theta,\phi,r_{ss}\right)\right|\sin\theta
d\theta d\phi$ (15)

For a single harmonic this simplifies to

$\displaystyle\Phi_{open}=r^{2}_{ss}B_{0}c_{\ell}\left(r_{ss}\right)\int_{0}^{2\pi}\int_{0}^{\pi}\left|Y_{\ell
m}\left(\theta,\phi\right)\right|\sin\theta d\theta d\phi$ (16)

We numerically evaluated the double integral to get a number for the analytic
open flux when comparing to values computed from the PFSS solver.

###  2.3 Magnetic field lines

In spherical coordinates the magnetic field tracing equations are

$\displaystyle\frac{dr}{ds}$ $\displaystyle=\hat{B}_{r}$ (17)
$\displaystyle\frac{d\theta}{ds}$ $\displaystyle=\frac{\hat{B}_{\theta}}{r}$
(18) $\displaystyle\frac{d\phi}{ds}$
$\displaystyle=\frac{\hat{B}_{\phi}}{r\sin\theta}$ (19)

where the paramater $s$ is the physical distance along the field line, and
$\hat{B}_{i}:=B_{i}/\left|\mathbf{B}\right|$ are components of a unit vector
pointing in the direction of the magnetic field. These form a set of three
coupled equations that can be integrated from an initial seed point to
evaluate coordinates along a magnetic field line.

In the case where the input field is proportional to a single spherical
harmonic, equations 12 – 14 can be substituted in and $ds$ eliminated to give

$\displaystyle\frac{d\theta}{dr}$
$\displaystyle=\left[\frac{1}{r}\frac{d_{\ell}\left(r\right)}{c_{\ell}\left(r\right)}\right]\left[\frac{1}{Y_{\ell
m}}\frac{\partial Y_{\ell m}\left(\theta,\phi\right)}{\partial\theta}\right]$
(20) $\displaystyle\frac{d\phi}{dr}$
$\displaystyle=\frac{1}{r\sin^{2}\theta}\frac{d_{\ell}\left(r\right)}{c_{\ell}\left(r\right)}\left[\frac{1}{Y_{\ell
m}}\frac{\partial Y_{\ell m}\left(\theta,\phi\right)}{\partial\phi}\right]$
(21) $\displaystyle\frac{d\theta}{d\phi}$
$\displaystyle=\sin^{2}\theta\left[\frac{\partial Y_{\ell
m}\left(\theta,\phi\right)}{\partial\theta}\right]\left[\frac{\partial Y_{\ell
m}\left(\theta,\phi\right)}{\partial\phi}\right]^{-1}$ (22)

Because the spherical harmonics are separable in $\theta,\phi$, equation 20 is
a function of only $r$ and $\theta$ and equation 22 is only a function of
$\phi$ and $\theta$. This decouples the three field line tracing equations,
allowing two of them to be integrated independently to give the longitude and
latitude as functions of radius. In the following subsections we integrate
these two equations and give analytic solutions for field lines in a spherical
harmonic solution.

####  2.3.1 The $\theta$ field line equation

To find $\theta\left(r\right)$ we start by separating variables and
integrating Equation 20

$\int_{\theta_{ss}}^{\theta}P_{\ell
m}\left(\bar{\theta}\right)\left[\frac{\partial P_{\ell
m}\left(\bar{\theta}\right)}{\partial\bar{\theta}}\right]^{-1}d\bar{\theta}=\int_{r_{ss}}^{r}\frac{1}{\bar{r}}\frac{d_{\ell}\left(\bar{r}\right)}{c_{\ell}\left(\bar{r}\right)}d\bar{r}$
(23)

Barred symbols denote dummy integration variables. To trace field lines down
from the source surface to the solar surface, the radial integration limits
are set to the source surface radius, $r_{ss}$ and the radial coordinate along
the field line, $r$. The $\theta$ integration limits set to a fixed initial
latitude $\theta_{ss}$, and the latitude along the field line, $\theta$.
Defining $\rho=r/r_{ss}$, the radial integration is given in equation 23 of
Gregory (2011) as

$\int_{r_{ss}}^{r}\frac{1}{\bar{r}}\frac{d_{\ell}\left(\bar{r}\right)}{c_{\ell}\left(\bar{r}\right)}d\bar{r}=\frac{1}{\ell}\frac{1}{\ell+1}\ln\left[\frac{\rho^{\ell}(2\ell+1)}{l\rho^{2\ell+1}+(\ell+1)}\right]$
(24)

The $\theta$ integration on the left hand side of equation 23 is more
complicated – Gregory (2011) give a general solution for the left hand side of
23 $\forall\ell$ and $m=0$, but for $\ell>1$ the solution is not available in
closed form. Instead we rearrange the integral on the left hand side of
equation 23 to define the function

$f_{\ell m}\left(\theta\right):=\exp\left[2l\int_{\theta_{ss}}^{\theta}P_{\ell
m}\left(\bar{\theta}\right)\left[\frac{\partial P_{\ell
m}\left(\bar{\theta}\right)}{\partial\bar{\theta}}\right]^{-1}d\bar{\theta}\right]$
(25)

With this definition, equations 24 and 25 are substituted into equation 23 to
give the equation for $\theta\left(r\right)$ along a field line

$f_{\ell m}\left(\theta\right)=f_{\ell
m}\left(\theta_{ss}\right)\left[\frac{\rho^{\ell}(2\ell+1)}{\ell\rho^{2\ell+1}+(\ell+1)}\right]^{\frac{2}{\ell+1}}$
(26)

The functions $f_{\ell m}\left(\theta\right)$ for low order spherical
harmonics are tabulated in Table 1. Solving this equation analytically
requires that the inverse of $f_{\ell m}$ exist in a closed form.

Table 1: Low order associated Legendre polynomials, and related functions. See Equation 25 for the definition of $f_{\ell m}$ and Equation 29 for the definition of $g_{\ell m}$. The functions listed are the same for $m\rightarrow-m$. l, m | $P_{\ell m}$ | $f_{\ell m}\left(\theta\right)$ | $g_{\ell m}\left(\theta\right)$
---|---|---|---
1, 0 | $\cos\theta$ | $\sin^{2}\theta$ |
1, 1 | $\sin\theta$ | $\cos^{2}\theta$ | $\frac{\sin\theta}{\cos\theta}$
2, 0 | $3\cos^{2}\theta-1$ |  |
2, 1 | $\cos\theta\sin\theta$ | $\cos 2\theta$ | $\frac{\sin\theta}{\sqrt{cos2\theta}}$
2, 2 | $\sin^{2}\theta$ | $\cos^{2}\theta$ | $\left(\frac{\sin\theta}{\cos\theta}\right)^{2}$
3, 0 | $5\cos^{3}\theta-3\cos\theta$ |  |
3, 1 | $(5\cos^{2}\theta-1)\sin\theta$ |  |
3, 2 | $\cos\theta\sin^{2}\theta$ | $3\cos(2\theta)+1$ | $\frac{\sin^{2}\theta}{2-3\sin^{2}\theta}$
3, 3 | $\sin^{3}\theta$ | $\cos^{2}\theta$ | $\frac{\sin\theta}{\cos\theta}$

####  2.3.2 The $\phi$ field line equation

Equation 22 is the field line equation relating $\theta$ and $\phi$. For
$m=0$, $B_{\phi}=0$ and this has the trivial solution $\phi=const$, so we only
consider the $m\neq 0$ case.

Separating variables and integrating gives

$\int_{\theta_{ss}}^{\theta}\frac{P_{\ell
m}}{\sin\bar{\theta}^{2}}\left[\frac{\partial P_{\ell
m}}{\partial\bar{\theta}}\right]^{-1}d\bar{\theta}=\int_{\phi_{ss}}^{\phi}Y_{\ell
m}\left[\frac{\partial Y_{\ell m}}{\partial\bar{\phi}}\right]^{-1}d\bar{\phi}$
(27)

The $\phi$ integration on the right hand side evaluates as

$\int_{\phi_{ss}}^{\phi}Y_{\ell m}\left[\frac{\partial Y_{\ell
m}}{\partial\bar{\phi}}\right]^{-1}d\bar{\phi}=\frac{1}{m^{2}}\begin{Bmatrix}\ln\frac{\sin
m\phi_{ss}}{\sin m\phi}\\\ \ln\frac{\cos m\phi_{ss}}{\cos
m\phi}\end{Bmatrix}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\begin{matrix}m>0\\\
m<0\end{matrix}$ (28)

For low order spherical harmonics the integral on the left hand side of 27 is
rearranged to define the function

$g_{\ell
m}\left(\theta\right):=\exp\left[m^{2}\int_{\theta_{ss}}^{\theta}\frac{P_{\ell
m}}{\sin\bar{\theta}^{2}}\left[\frac{\partial P_{\ell
m}}{\partial\bar{\theta}}\right]^{-1}d\bar{\theta}\right]$ (29)

The functions $g_{\ell m}$ for low order spherical harmonics are given in
Table 1. Substituting equations 28 and 29 into 27 gives the the equation for
$\phi$ as a function of $r$ along a field line

$\begin{Bmatrix}\frac{\sin m\phi}{\sin m\phi_{ss}}\\\ \frac{\cos m\phi}{\cos
m\phi_{ss}}\end{Bmatrix}=\frac{g_{\ell m}\left(\theta\right)}{g_{\ell
m}\left(\theta_{ss}\right)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\begin{matrix}m>0\\\
m<0\end{matrix}$ (30)

Once $\theta(r)$ is known from solving equation 26, it is used in Equation 30
to solve for $\phi(\theta(r))$. Unlike $f_{\ell m}$, the inverse of $g_{\ell
m}$ does not need to exist in closed form to derive $\phi$ along a field line.

##  3 The pfsspy solver

We briefly recount how pfsspy calculates the PFSS solution, in order to
understand the user configurable options that affect the accuracy of the
solver. Details on the numerical scheme are available in the numerical methods
document archived alongside the software (Stansby et al., 2022).

pfsspy uses a finite difference method (van Ballegooijen et al., 2000,
Appendix B) to calculate the magnetic vector potential on a grid regularly
spaced in $z=\log\left(r/r_{\odot}\right)$, $\lambda=\cos\theta$, $\phi$. The
number of grid points in the angular dimensions, $n_{\phi}$ and $n_{\theta}$,
are fixed by the resolution of the input $B_{r\odot}$ grid which must span the
full sphere. Since the smallest magnetograms widely used from the Sun have a
grid size of of 360 $\times$ 180 in (longitude, latitude), for simplicity we
keep this as the fixed angular grid size throughout the tests. The number of
radial grid points, $n_{r}$ is user configurable.

To trace magnetic field lines pfsspy offers two different field line tracers.
Here we perform comparisons with the FORTRAN
implementation111https://github.com/dstansby/streamtracer, which is the
fastest of the two.

The field line equations in the coordinates that pfsspy uses are

$\displaystyle\frac{d\rho}{ds}$ $\displaystyle=\hat{B}_{\rho}$ (31)
$\displaystyle\frac{d\lambda}{ds}$
$\displaystyle=\sqrt{1-\lambda^{2}}\hat{B}_{\lambda}$ (32)
$\displaystyle\frac{d\phi}{ds}$
$\displaystyle=\frac{1}{\sqrt{1-\lambda^{2}}}\hat{B}_{\phi}$ (33)

In the field line tracer three equations are integrated numerically using a
4th order Runge-Kutta method with a fixed discrete step size, $\Delta s$,
which is specified by the user as a fraction of the grid cell size in the
radial direction, $\Delta r=\log(r_{ss}/r_{\odot})/n_{r}$. Because the step
size $\Delta s$ is kept constant relative to the log-scaled grid, larger
physical steps are taken further away from the solar surface at larger values
of $r$.

Table 2: Configurable parameters in the pfsspy solver and their values for each figure. | $n_{\phi},n_{\theta}$ | $n_{r}$ | $\Delta s$
---|---|---|---
Figure 1 | (360, 180) | 40 | 1
Figure 2 | (360, 180) | Varied | 1
Figure 3 | (360, 180) | 40 | 1
Figure 4 | (360, 180) | 40 | Varied

##  4 Analytic solutions as test problems

In this section we compare the magnetic field on the source surface (Section
4.1), the open flux (Section 4.2), and the field line connectivity (Section
4.3) between analytic and numerical pfsspy solutions. For clarity analysis is
limited to harmonics with $\ell\leq 5$. For the Sun these are the dominant
harmonics at all times in the solar cycle (DeRosa et al., 2012).

A range of software is used for creating comparisons, including pfsspy
(Stansby et al., 2020b), numpy (Harris et al., 2020), scipy (Virtanen et al.,
2020), pandas (Reback et al., 2021), sympy (Meurer et al., 2017) and astropy
(The Astropy Collaboration et al., 2018). Code for producing the comparisons
is available in version 1.1.0 of the pfsspy
repository222https://github.com/dstansby/pfsspy and is archived at Stansby et
al. (2022).

In the numerical solutions we set default the grid sizes to $n_{\phi}=360$,
$n_{\theta}=180$, $n_{r}=40$, and the default field line tracing step size to
$\Delta s=1$. The following subsections show a series of different tests,
where one of these parameters may be varied. In all tests the source surface
height is set to $r_{ss}=2r_{\odot}$. Each test has one or more corresponding
figures. For a summary of which parameters are fixed or varied for each
figure, see Table 2.

###  4.1 $\mathbf{B}$ on the source surface

Figure 1: Comparison between analytic (orange) and pfsspy (blue) spherical
harmonic solutions to the PFSS equations as functions of latitude at
$r=r_{ss},\phi=15^{\circ}$. Each panel represents a different spherical
harmonic number, and plotted is the radial component of magnetic field. Dashed
horizontal lines show $B_{r}=0$.

The two transverse magnetic field components, $B_{\phi,\theta}$ must always be
zero on the source surface within a PFSS solution. This is always the case in
in pfsspy which forces $B_{\phi,\theta}=0$ at $r=r_{ss}$.

To compare $B_{r}$ at $r=r_{ss}$, equation 14 is evaluated at the same points
as the numerical pfsspy solution. 1D cuts of both the analytic and pfsspy
solutions at a constant longitude of $\phi=15^{\circ}$ are compared in Figure
1. pfsspy reproduces the analytic solutions well, with the magnitude of
solutions slightly larger than the analytic solutions. This suggests that
pfsspy systematically over-estimates the total unsigned open flux, which we
investigate quantitively in the next section.

###  4.2 Open flux

Figure 2: Ratio of unsigned open flux in numerical pfsspy solutions to
analytic solutions as a function of number of radial cells within the pfsspy
solution for various spherical harmonics.

Equation 16 gives an analytic integral for the unsigned open magnetic flux,
which we evaluate using the nquad integration function in scipy. To check the
accuracy of the integration we successfully verify the result against the
analytical solution to the integral for $\ell,m=1,0$:

$\displaystyle\Phi_{open}/B_{0}$
$\displaystyle=c_{1}\left(r_{ss}\right)\int_{0}^{2\pi}d\phi\int_{0}^{\pi}\left|\cos\theta\right|\sin\theta
d\theta$ (34) $\displaystyle=\frac{3}{2z_{ss}^{3}+1}\cdot 2\pi\cdot
2\int_{0}^{\pi/2}\cos\theta\sin\theta d\theta$ (35)
$\displaystyle=\frac{6\pi}{2z_{ss}^{3}+1}$ (36)

To evaluate the total unsigned open flux within a pfsspy result the magnetic
field values on the source surface is multiplied by the area of their
associated cell and summed over the full sphere.

Figure 2 shows the ratio of open flux in pfsspy results to the analytic value
as a function of the number of radial grid points for spherical harmonics up
to $\ell=5$. As hinted from Figure 1, the unsigned open fluxes within pfsspy
are systematically larger than the analytic solution. This difference
decreases with an increasing number of radial grid cells, approaching
asymptotic values of $\leq 5\%$ at around $n_{r}=40$. This motivates our
choice of $n_{r}=40$ for the default number of radial grid cells – with more
radial cells there is not a significant reduction in the open flux error.

At larger values of $\ell$ the error in unsigned open flux increases from
$\sim$1% at $\ell=1$ to $\sim$5% at $\ell=5$. The contribution of each
spherical harmonic to the open flux on the source surface reduces as
$(r_{ss}/r_{\odot})^{-\ell-2}$ (Equation 10). As long as the spherical
harmonics of the input $B_{r}$ map do not scale exponentially with $\ell$, the
increasing open flux error with $\ell$ within pfsspy will be suppressed by the
$(r_{ss}/r_{\odot})^{-\ell-2}$ factor, preventing the total open flux error
growing when summing over multiple harmonics.

###  4.3 Field line connectivity

From a regularly spaced grid of field line seed points on the source surface
analytic solar surface field line footpoints are calculated using Equations 26
and 30 and the numerical field line footpoints computed using pfsspy . Figure
3 shows the difference between the analytic and traced solar surface
footpoints in latitude (top panel) and longitude (lower panel) for
$\ell,m=3,3$ a step size of $\Delta s=1$, and a radial grid size of
$n_{r}=40$. Overall the errors are small in both latitude and longitude at
$\leq 0.5^{\circ}$.

Figure 3: The difference between numerically traced and analytic solar surface
footpoint latitude (top panel) and longitude (bottom panel), for $l=3,m=3$ and
a tracer step size of $\Delta s=1$. Field lines are traced downwards from an
evenly spaced grid on the source surface.

There are two situations where pfsspy fails to correctly trace the field
lines. The first is near polarity inversion lines, which occur at
$30^{\circ},90^{\circ},150^{\circ},...$ in Figure 3. Traced field lines
started near polarity inversion lines on the source surface turn around and
return immediately to the source surface. This is due to finite resolution
effects within the model. pfsspy automatically tags these field lines as
incorrectly traced, and they show up as thin white strips in Figure 3. This
occurs at all spherical harmonic numbers, but only for a limited number of
field lines near polarity inversion lines.

The second situation occurs near polarity inversion lines at high latitudes,
when the inversion line in the model is slightly displaced from the expected
inversion line. This causes field lines to deviate strongly in longitude from
their analytic solution, as seen in the bottom panel of Figure 3 where a few
points have large errors outside the range of the errorbars. This is due to a
combination of finite resolution effects and no attempt by the tracer to
handle the spherical coordinate singularity at the poles. Up to $\ell=3$ it
only occurs for $\ell=3$ and $m=-3,3$.

Figure 4: Maximum error in field line footpoints as a function of tracer step
size for various spherical harmonics. Error in longitude is shown in blue and
error in latitude in orange. Errors with absolute values $>10^{\circ}$ are
excluded in these plots - this can occur close to polarity inversion lines on
the source surface.

To quantify how errors in field line tracing vary with tracing step size we
produce error maps like Figure 3 for a range of spherical harmonic numbers.
Figure 4 shows how the maximum error in field line footpoint across a whole
map varies with step size for different spherical harmonics. Errors with
absolute values $>10^{\circ}$ are excluded from these plots, to exclude any
points that are incorrectly traced near the polarity inversion line (see
previous paragraph).

With decreasing magnetic field tracing step size, the maximum error drops off
until reaching a plateau at a step size of $\sim$ 4\. Because the field line
integrator uses a 4th order Runge-Kutta method, it samples the field at 1/4 of
the step size, so a levelling off of the error at a step size of around 4 is
expected. Below this the integrator samples at sub-grid resolution, and the
error is dominated by the error in the magnetic field solution itself. This
result justifies our default choice of $\Delta s=1$ for the other tests.
Making the integration step size any smaller than $\Delta s=1$ increases
computation and memory use, but does not increase the accuracy at which the
field lines are traced.

##  5 Conclusions

We have derived a set of analytical closed solutions for the PFSS equations
when the input magnetic field is a single spherical harmonic (Section 2.1),
along with analytical equations for magnetic field lines traced through these
solutions (Section 2.3). These solutions have then been used to test and
quantify the accuracy of the numerical pfsspy solver (Section 4).

This set of tests reveal both the accuracy and limitations of the pfsspy
solver. The total open magnetic flux is systematically overestimated within
the pfsspy solver, but when the number of radial grid cells exceeds 40 is
accurate to within 6% of the true value for individual spherical harmonics
(Figure 2). When tracing field lines from the source surface to the solar
surface, pfsspy is always accurate to within 0.5∘ with a small enough tracer
step size (Figure 4), but large errors in field line tracing can occur near
polarity inversion lines and at the poles (Figure 3). Note that our tests use
spherical harmonic solutions which are relatively smooth functions, so
different tests are needed to understand how accurate the field line tracing
is in areas of rapidly changing magnetic field, e.g., surrounding active
regions.

These test problems form a solid benchmark for PFSS solvers and their ability
to reproduce stellar potential magnetic fields accurately. We recommend that
all PFSS solvers be tested against the analytical test problems presented
here.

D.S. and D.V. are supported by STFC Consolidated Grant ST/S000240/1. D.V. is
supported by STFC Ernest Rutherford Fellowship ST/P003826/1.

## References

  * Altschuler & Newkirk (1969) Altschuler, M. D., & Newkirk, G. 1969, Solar Physics, 9, 131, doi: 10.1007/BF00145734
  * Badman et al. (2020) Badman, S. T., Bale, S. D., Oliveros, J. C. M., et al. 2020, The Astrophysical Journal Supplement Series, 246, 23, doi: 10.3847/1538-4365/ab4da7
  * Baker et al. (2021) Baker, D., Mihailescu, T., Démoulin, P., et al. 2021, Solar Physics, 296, 103, doi: 10.1007/s11207-021-01849-7
  * Boe et al. (2020) Boe, B., Habbal, S., & Druckmüller, M. 2020, The Astrophysical Journal, 895, 123, doi: 10.3847/1538-4357/ab8ae6
  * Caplan et al. (2021) Caplan, R. M., Downs, C., Linker, J. A., & Mikic, Z. 2021, The Astrophysical Journal, 915, 44, doi: 10.3847/1538-4357/abfd2f
  * DeRosa et al. (2012) DeRosa, M. L., Brun, A. S., & Hoeksema, J. T. 2012, The Astrophysical Journal, 757, 96, doi: 10.1088/0004-637X/757/1/96
  * Fargette et al. (2021) Fargette, N., Lavraud, B., Rouillard, A. P., et al. 2021, The Astrophysical Journal, 919, 96, doi: 10.3847/1538-4357/ac1112
  * Gregory (2011) Gregory, S. G. 2011, American Journal of Physics, 79, 461, doi: 10.1119/1.3549206
  * Hakamada (1995) Hakamada, K. 1995, Solar Physics, 159, 89, doi: 10.1007/BF00733033
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Jardine et al. (2017) Jardine, M., Vidotto, A. A., & See, V. 2017, Monthly Notices of the Royal Astronomical Society: Letters, 465, L25, doi: 10.1093/mnrasl/slw206
  * Jiang & Feng (2012) Jiang, C., & Feng, X. 2012, Solar Physics, 281, 621, doi: 10.1007/s11207-012-0074-x
  * Kochukhov (2020) Kochukhov, O. 2020, The Astronomy and Astrophysics Review, 29, 1, doi: 10.1007/s00159-020-00130-3
  * Mackay & Yeates (2012) Mackay, D. H., & Yeates, A. R. 2012, Living Reviews in Solar Physics, 9, 6, doi: 10.12942/lrsp-2012-6
  * Meurer et al. (2017) Meurer, A., Smith, C. P., Paprocki, M., et al. 2017, PeerJ Computer Science, 3, e103, doi: 10.7717/peerj-cs.103
  * Neugebauer et al. (1998) Neugebauer, M., Forsyth, R. J., Galvin, A. B., et al. 1998, Journal of Geophysical Research: Space Physics, 103, 14587, doi: 10.1029/98JA00798
  * Reback et al. (2021) Reback, J., jbrockmendel, McKinney, W., et al. 2021, Pandas-Dev/Pandas: Pandas 1.3.4, Zenodo, doi: 10.5281/zenodo.5574486
  * Saikia et al. (2020) Saikia, S. B., Jin, M., Johnstone, C. P., et al. 2020, Astronomy & Astrophysics, 635, A178, doi: 10.1051/0004-6361/201937107
  * Schatten et al. (1969) Schatten, K. H., Wilcox, J. M., & Ness, N. F. 1969, Solar Physics, 6, 442, doi: 10.1007/BF00146478
  * Sod (1978) Sod, G. A. 1978, Journal of Computational Physics, 27, 1, doi: 10.1016/0021-9991(78)90023-2
  * Stansby et al. (2022) Stansby, D., Badman, S., Ancellin, M., & Barnes, W. 2022, Dstansby/Pfsspy: Pfsspy 1.1.0, Zenodo, doi: 10.5281/zenodo.5879440
  * Stansby et al. (2020a) Stansby, D., Berčič, L., Matteini, L., et al. 2020a, Astronomy & Astrophysics, doi: 10.1051/0004-6361/202039789
  * Stansby et al. (2021) Stansby, D., Green, L. M., van Driel-Gesztelyi, L., & Horbury, T. S. 2021, Solar Physics, 296, 116, doi: 10.1007/s11207-021-01861-x
  * Stansby et al. (2020b) Stansby, D., Yeates, A., & Badman, S. T. 2020b, Journal of Open Source Software, 5, 2732, doi: 10.21105/joss.02732
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, The Astronomical Journal, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Tóth et al. (2011) Tóth, G., van der Holst, B., & Huang, Z. 2011, The Astrophysical Journal, 732, 102, doi: 10.1088/0004-637X/732/2/102
  * van Ballegooijen et al. (2000) van Ballegooijen, A. A., Priest, E. R., & Mackay, D. H. 2000, The Astrophysical Journal, 539, 983, doi: 10.1086/309265
  * van Driel-Gesztelyi et al. (2012) van Driel-Gesztelyi, L., Culhane, J. L., Baker, D., et al. 2012, Solar Physics, 281, 237, doi: 10.1007/s11207-012-0076-8
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang & Sheeley (1992) Wang, Y.-M., & Sheeley, Jr., N. R. 1992, The Astrophysical Journal, 392, 310, doi: 10.1086/171430


# Finite perturbation theory for the relativistic Coulomb problem

Scott E. Hoffmann  School of Mathematics and Physics,
The University of Queensland,
Brisbane, QLD 4072
Australia<EMAIL_ADDRESS>
###### Abstract

We present a novel form of relativistic quantum mechanics and demonstrate how
to solve it using a recently derived unitary perturbation theory, within
partial wave analysis. The theory is tested on a relativistic problem, with
two spinless, equal mass particles, in which the interaction is entirely given
by a Coulomb potential. As such, it is not meant to reproduce experimental
results for the scattering of two electrons, but is intended as a test of our
calculation methods. We find that this perturbation theory gives finite
results at second order. This is unlike other versions of perturbation theory,
which find divergent results at second and all higher orders. We calculate
differential cross sections in the nonrelativistic regime, where we find
excellent agreement with the Rutherford formula. Then, well into the
relativistic regime, we find differential cross sections with similar shapes
to the Møller formula and differing from that formula by less than an order of
magnitude.

##  I Introduction

There are two aims of this paper. The first is to apply a novel perturbation
theory to the problem of the Coulomb scattering of two (spinless) charged
particles. This problem has been treated before with other forms of
perturbation theory, with the well-known result that all corrections at second
and higher order diverge (Dalitz [1]). The reason, in this context, is that
certain integrals over particle separation do not converge because the Coulomb
potential falls off insufficiently fast at large particle separation. We note
that even the first-order result is problematic (Collas [2]).

The new method used here (Hoffmann [3]) involves working in momentum space and
a prescription to use principal part integration on the momentum integrals. We
will find that this procedure produces finite results in agreement with
expectations when applied to the Coulomb potential.

The other aim of this paper is to show that some features of relativity can be
easily incorporated into our model. For our purposes, relativity gives a
different dependence of single-particle energy on momentum from the
nonrelativistic form, and can add new terms to the interaction. For the
latter, an example is the emergence of the spin-orbit coupling term when the
Foldy-Wouthuysen transformation is applied to the Dirac equation (Foldy and
Wouthuysen [4]). In our model, we will introduce the relativistic form of the
single-particle energy. We will consider an entirely Coulomb potential in this
spinless case.

We note that relativity constrains measurements made in different frames, but
does not tell us about measurements in one frame. The position profile of a
porcupine depends on how that animal has grown, and varies from animal to
animal. Relativity tells us only how a moving frame would measure a changed
porcupine profile. We will propose that there must be a four-component
potential operator that transforms as a four-vector to have a covariant
relativistic quantum mechanics. We will find that in a centre of mass (CM)
frame, the spatial components must vanish. The form of the potential in a CM
frame is arbitrary for our considerations, but is determined by the model of
microscopic interactions in quantum electrodynamics (Weinberg [5]). Future
work will define these potentials and explore the consequences.

Throughout this paper, we use Heaviside-Lorentz units, in which
$\hbar=c=\epsilon_{0}=\mu_{0}=1.$

##  II Finite, relativistic, perturbation theory

###  II.1 Relativistic Schrödinger equations

We are motivated by the success of nonrelativistic quantum mechanics, in which
the total Hamiltonian is written as the sum of a free part and a part
responsible for interaction (a potential $\hat{V}$):

$\hat{H}_{i}=\hat{H}_{f}+\hat{V}.$ (1)

Here $f$ indicates free, $i$ indicates interacting, hats indicate operators.

Thus we propose for relativistic quantum mechanics the four equations
($\mu=0,1,2,3$)

$\hat{P}_{i}^{\mu}=\hat{P}_{f}^{\mu}+\hat{V}^{\mu}.$ (2)

Here the potential of nonrelativistic quantum mechanics becomes the zero
component of what must transform as a four-vector, as do the components of
$P_{f}^{\mu}$ and $P_{i}^{\mu},$ to guarantee relativistic covariance. This is

$U^{\dagger}(L)\,V^{\mu}\,U(L)=L_{\phantom{\mu}\nu}^{\mu}\,V^{\nu},$ (3)

where $U(L)$ is the unitary representation of a general Lorentz
transformation, $L,$ and $L_{\phantom{\mu}\nu}^{\mu}$ is the corresponding
coordinate transformation matrix. For processes involving photons, we assert
that $V^{\mu}$ must be gauge invariant to arrive at physical results
independent of the gauge.

For the scattering of two spinless “electrons”, we choose as free basis
vectors $|\,P,l,m;f\,\rangle.$ In the CM frame, these are defined in terms of
the individual momentum eigenvectors as

$|\,(2\omega^{\mathrm{\scriptscriptstyle
CM}},\boldsymbol{0}),l,m;f\,\rangle=\int
d^{2}\hat{p}^{\mathrm{\scriptscriptstyle
CM}}\,|\,+\boldsymbol{p}^{\mathrm{\scriptscriptstyle
CM}},-\boldsymbol{p}^{\mathrm{\scriptscriptstyle
CM}};f\,\rangle\,p^{\mathrm{\scriptscriptstyle
CM}}\,Y_{lm}(\hat{\boldsymbol{p}}^{\mathrm{\scriptscriptstyle CM}}).$ (4)

Then a boost takes them to general four-momentum

$|\,P,l,m;f\,\rangle=U(\Lambda(\frac{\boldsymbol{P}}{E}))\,|\,(2\omega^{\mathrm{\scriptscriptstyle
CM}},0),l,m;f\,\rangle.$ (5)

We can ignore the total linear momentum, $\boldsymbol{P},$ by using the
normalization

$\langle\,+\boldsymbol{p}_{1}^{\mathrm{\scriptscriptstyle
CM}},-\boldsymbol{p}_{1}^{\mathrm{\scriptscriptstyle
CM}};f\,|\,+\boldsymbol{p}_{2}^{\mathrm{\scriptscriptstyle
CM}},-\boldsymbol{p}_{2}^{\mathrm{\scriptscriptstyle
CM}};f\,\rangle=\delta^{3}(\boldsymbol{p}_{1}^{\mathrm{\scriptscriptstyle
CM}}-\boldsymbol{p}_{2}^{\mathrm{\scriptscriptstyle CM}}),$ (6)

which has no consequences.

These are eigenvectors of total four-momentum
$P^{\mu}=(E,\boldsymbol{P})^{\mu}.$ They are also eigenvectors of
$\boldsymbol{L}^{2}$ and $L_{z},$ where $\boldsymbol{L}$ is the total angular
momentum (which in this case is entirely orbital angular momentum) in the CM
frame (where $\boldsymbol{P}=0$). In the more general case, where the two
particles can be any combination of electrons, positrons and photons, the
basis vectors take the form $|P,J,M,\lambda_{1},\lambda_{2};f\,\rangle,$ where
$J$ and $M$ are the quantum numbers associated with the total CM-frame angular
momentum, $\boldsymbol{J},$ and $\lambda_{1}$ and $\lambda_{2}$ are the CM-
frame helicities of the two particles.

It is advantageous to work in momentum space, where the square roots in the
single-particle energies,
$\omega(\boldsymbol{p})=\sqrt{\boldsymbol{p}^{2}+m_{e}^{2}},$ pose no problem.
At this point there is no need to introduce spacetime coordinates for the two
particles. Note that to see a potential with spacetime dependence requires
taking matrix elements between two state vectors, at least one of which has
spacetime dependence.

Corresponding to each free basis vector is an interacting basis vector that
solves the Schrödinger equations (eq. (2)). We will explore this
correspondence shortly.

We take matrix elements

$\langle\,P_{f},l_{f},m_{f};f\,|\,\hat{P}_{i}^{\mu}\,|\,P_{i},l_{i},m_{i};i\,\rangle=\langle\,P_{f},l_{f},m_{f};f\,|\,\hat{P}_{f}^{\mu}+\hat{V}^{\mu}\,|\,P_{i},l_{i},m_{i};i\,\rangle,$
(7)

which reduce to

$\displaystyle(P_{f}^{\mu}-P_{i}^{\mu})\,\langle\,P_{f},l_{f},m_{f};f\,|\,P_{i},l_{i},m_{i};i\,\rangle$
$\displaystyle=-\langle\,P_{f},l_{f},m_{f};f\,|\,\hat{V}^{\mu}\,|\,P_{i},l_{i},m_{i};i\,\rangle$
$\displaystyle=-\int
d^{4}P_{f}^{\prime}\sum_{l_{f}^{\prime}=0}^{\infty}\sum_{m_{f}^{\prime}=-l_{f}^{\prime}}^{l_{f}^{\prime}}\langle\,P_{f},l_{f},m_{f};f\,|\,\hat{V}^{\mu}\,|\,P_{f}^{\prime},l_{f}^{\prime},m_{f}^{\prime};f\,\rangle\langle\,P_{f}^{\prime},l_{f}^{\prime},m_{f}^{\prime};f\,|\,P_{i},l_{i},m_{i};i\,\rangle,$
(8)

a set of integral equations. We ensure that the free basis vectors are
orthonormalized and complete so that we can insert the completeness relation
as we have just done.

###  II.2 Simplification in the centre of mass frame

A simplification is possible. We will find below that it is possible to define
a unitary transformation that relates free and interacting basis vectors with
the same quantum numbers:

$|\,P,l,m;i\,\rangle=U_{if}\,|\,P,l,m;f\,\rangle.$ (9)

Note that in the quantum electrodynamics formalism (Weinberg [5]), the
electron is given a free mass different from its mass after including
electromagnetic interactions. But the difference of these masses is calculated
to diverge and only the interacting mass is observable. As we will see, the
correspondence proposed here produces no divergences.

The transformation of operators that gives the same expectation values is

$\langle\,P,l,m;i\,|\,P_{i}^{\mu}\,|\,P,l,m;i\,\rangle=\langle\,P,l,m;f\,|\,P_{f}^{\mu}\,|\,P,l,m;f\,\rangle\Rightarrow
P_{i}^{\mu}=U_{if}\,P_{f}^{\mu}\,U_{if}^{\dagger}.$ (10)

Note that we will have similar relations for the transformation of rotation
generators ($\boldsymbol{J}$) and boost generators ($\boldsymbol{K}$):

$J_{i}^{k}=U_{if}\,J_{f}^{k}\,U_{if}^{\dagger}\quad\mathrm{and}\quad
K_{i}^{k}=U_{if}\,K_{f}^{k}\,U_{if}^{\dagger}\quad\mathrm{for}\ k=1,2,3.$ (11)

Consequently, if

$\hat{\boldsymbol{P}}_{f}\,|\,P,l,m;f\,\rangle=\boldsymbol{0},$ (12)

then

$\hat{\boldsymbol{P}}_{i}\,|\,P,l,m;i\,\rangle=U_{if}\,\hat{\boldsymbol{P}}_{f}\,U_{if}^{\dagger}\,U_{if}\,|\,P,l,m;f\,\rangle=U_{if}\,\hat{\boldsymbol{P}}_{f}\,|\,P,l,m;f\,\rangle=\boldsymbol{0}.$
(13)

So, from eq. (), we see the vanishing of the matrix elements of the three-
vector part of the potential, $\hat{\boldsymbol{V}},$ in the one CM frame that
is the same in the free and interacting theories:

$\langle\,(E_{f},\boldsymbol{0}),l_{f},m_{f};f\,|\,\hat{\boldsymbol{V}}\,|\,(E_{i},\boldsymbol{0}),l_{i},m_{i};i\,\rangle=\boldsymbol{0}.$
(14)

In what follows, we will perform calculations only in the CM frame, knowing
that we could transform differential cross sections to other frames, if
necessary.

###  II.3 Translational and rotational invariance

We define $\hat{V}=\hat{V}^{0}.$ We know

$\langle\,P_{f},l_{f},m_{f};f\,|\,P_{i},l_{i},m_{i};i\,\rangle\propto\delta^{3}(\boldsymbol{P}_{f}-\boldsymbol{P}_{i})$
(15)

by translational invariance, so we can remove all reference to the three-
momenta and arrive at the simplified equation in the CM frame (understanding
$\boldsymbol{P}_{f}=\boldsymbol{P}_{i}=0$ throughout)

$(E_{f}-E_{i})\,\langle\,E_{f},l_{f},m_{f};f\,|\,E_{i},l_{i},m_{i};i\,\rangle=-\int_{2m_{e}}^{\infty}dE_{f}^{\prime}\sum_{l_{f}^{\prime}=0}^{\infty}\sum_{m_{f}^{\prime}=-l_{f}^{\prime}}^{l_{f}^{\prime}}\langle\,E_{f},l_{f},m_{f};f\,|\,\hat{V}\,|\,E_{f}^{\prime},l_{f}^{\prime},m_{f}^{\prime};f\,\rangle\langle\,E_{f}^{\prime},l_{f}^{\prime},m_{f}^{\prime};f\,|\,E_{i},l_{i},m_{i};i\,\rangle.$
(16)

We are considering the case of a potential that is rotationally invariant
($[\boldsymbol{J}_{i},\hat{V}]=0$). Thus, by the Wigner-Eckhart theorem
(Messiah [6]), the matrix elements take the form

$\langle\,E_{f},l_{f},m_{f};f\,|\,\hat{V}\,|\,E_{f}^{\prime},l_{f}^{\prime},m_{f}^{\prime};f\,\rangle\propto\delta_{l_{f}l_{f}^{\prime}}\delta_{m_{f}m_{f}^{\prime}}.$
(17)

Furthermore, if the potential is rotationally invariant, we expect that the
free to interacting transformation will also be rotationally invariant,
preserving the quantum numbers $l$ and $m$. A further reduction, ignoring $m,$
gives

$(E_{f}-E_{i})\,\langle\,E_{f},l;f\,|\,E_{i},l;i\,\rangle=-\int_{2m_{e}}^{\infty}dE_{f}^{\prime}\langle\,E_{f},l;f\,|\,\hat{V}\,|\,E_{f}^{\prime},l;f\,\rangle\langle\,E_{f}^{\prime},l;f\,|\,E_{i},l;i\,\rangle.$
(18)

###  II.4 Eigenvectors of particle separation

The phase shifts that we want to calculate come from the asymptotic behaviour
of interacting position wavefunctions as the separation, $r,$ of the two
particles increases without bound. Thus we need to construct eigenvectors of
the operator associated with this separation. Note that we cannot here use the
techniques of nonrelativistic quantum mechanics, where the kinetic energy
depends quadratically on the momentum, to form the total linear momentum and
the relative momentum.

We consider two different individual particle translations applied to a free
individual momentum eigenvector:

$U(T_{1}(+\boldsymbol{r}/2))\,U(T_{2}(-\boldsymbol{r}/2))\,|\,\frac{1}{2}\boldsymbol{P}+\boldsymbol{p},\frac{1}{2}\boldsymbol{P}-\boldsymbol{p};f\,\rangle=|\,\frac{1}{2}\boldsymbol{P}+\boldsymbol{p},\frac{1}{2}\boldsymbol{P}-\boldsymbol{p}\,\rangle\,e^{-i\boldsymbol{p}\cdot\boldsymbol{r}}.$
(19)

This transformation changes the separation of the two particles by
$\boldsymbol{r}$ but leaves the average position unchanged. So we define
simultaneous eigenvectors of separation and total three-momentum as

$|\,\boldsymbol{r},\boldsymbol{P};f\,\rangle=\int\frac{d^{3}p}{(2\pi)^{\frac{3}{2}}}\,|\,\frac{1}{2}\boldsymbol{P}+\boldsymbol{p},\frac{1}{2}\boldsymbol{P}-\boldsymbol{p};f\,\rangle\,e^{-i\boldsymbol{p}\cdot\boldsymbol{r}},$
(20)

noting that

$\hat{\boldsymbol{r}}=i\frac{\partial}{\partial\boldsymbol{p}}=i\frac{\partial}{\partial\boldsymbol{p}_{1}}-i\frac{\partial}{\partial\boldsymbol{p}_{2}}$
(21)

and $\hat{\boldsymbol{P}}$ commute.

Then we form eigenvectors of the magnitude of the separation and of orbital
angular momentum,

$|\,r,l,m,\boldsymbol{P};f\,\rangle=\int
d^{2}\hat{r}\,|\,\boldsymbol{r},\boldsymbol{P};f\,\rangle\,r\,Y_{lm}(\hat{\boldsymbol{r}}),$
(22)

where the argument of the spherical harmonic is here a unit vector, not an
operator. We note that
$\hat{\boldsymbol{L}}=\hat{\boldsymbol{r}}\times\hat{\boldsymbol{p}}$ defined
in the CM frame also commutes with $\boldsymbol{P}.$

Note that we will only be requiring these basis vectors for the special case
$\boldsymbol{P}=\boldsymbol{0},$ where $r$ is the separation and $l,m$ label
the orbital angular momentum, both in the CM frame. Note also that these basis
vectors will have calculable but very complicated Lorentz transformation
properties, as with any position eigenvectors.

If we choose the orthonormalization

$\langle\,\frac{1}{2}\boldsymbol{P}_{a}+\boldsymbol{p}_{a},\frac{1}{2}\boldsymbol{P}_{a}-\boldsymbol{p}_{a};f\,|\,\frac{1}{2}\boldsymbol{P}_{b}+\boldsymbol{p}_{b},\frac{1}{2}\boldsymbol{P}_{b}-\boldsymbol{p}_{b};f\,\rangle=\delta^{3}(\boldsymbol{P}_{a}-\boldsymbol{P}_{b})\,\delta^{3}(\boldsymbol{p}_{a}-\boldsymbol{p}_{b}),$
(23)

then we find

$\langle\,r_{a},l_{a},m_{a},\boldsymbol{P}_{a};f\,|\,r_{b},l_{b},m_{b},\boldsymbol{P}_{b};f\,\rangle=\delta^{3}(\boldsymbol{P}_{a}-\boldsymbol{P}_{b})\,\delta(r_{a}-r_{b})\,\delta_{l_{a}l_{b}}\delta_{m_{a}m_{b}}.$
(24)

We find the relation between eigenvectors of momentum magnitude and
eigenvectors of separation magnitude:

$\langle\,r,l_{a},m_{a},\boldsymbol{P}_{a};f\,|\,p,l_{b},m_{b},\boldsymbol{P}_{b};f\,\rangle=\delta^{3}(\boldsymbol{P}_{a}-\boldsymbol{P}_{b})\,\delta_{l_{a}l_{b}}\delta_{m_{a}m_{b}}\,\sqrt{\frac{2}{\pi}}\,pr\,j_{l_{a}}(pr)=\delta^{3}(\boldsymbol{P}_{a}-\boldsymbol{P}_{b})\,\delta_{l_{a}l_{b}}\delta_{m_{a}m_{b}}\,y_{l_{a}}^{(f)}(r,p),$
(25)

in terms of the free spherical waves, $y_{l_{a}}^{(f)}(r,p),$ familiar from
partial wave analysis (Messiah [6]). This motivates our definition for the
interacting case,

$\langle\,r,l_{a},m_{a},P_{a};f\,|\,p,l_{b},m_{b},P_{b};i\,\rangle\equiv\delta^{3}(\boldsymbol{P}_{a}-\boldsymbol{P}_{b})\,\delta_{l_{a}l_{b}}\delta_{m_{a}m_{b}}\,y_{l_{a}}^{(i)}(r,p),$
(26)

in terms of the interacting wavefunctions, $y_{l_{a}}^{(i)}(r,p).$

Note that the position eigenvectors are common to the free and interacting
theories. We do not define interacting position eigenvectors.

###  II.5 Unitary perturbation theory

A unitary transformation connecting the free and interacting theories, in the
nonrelativistic case, was introduced in Hoffmann [3], and a perturbation
series was developed to second order in the coupling strength. Only minor
changes are needed to treat the relativistic case.

In the CM frame, the total energies are

$E(p)=2\sqrt{p^{2}+m_{e}^{2}},$ (27)

where $p$ is the magnitude of both single-particle momenta. Then a change of
quantum numbers to eq. (18) gives

$(E(k)-E(p))\,\langle\,k,l;f\,|\,p,l;i\,\rangle=-\int_{0}^{\infty}dk^{\prime}\,\langle\,k,l;f\,|\,\hat{V}\,|\,k^{\prime},l;f\,\rangle\langle\,k^{\prime},l;f\,|\,p,l;i\,\rangle.$
(28)

The normalizations are

$\langle\,E_{1},l;f/i\,|\,E_{2},l;f/i\,\rangle=\delta(E_{1}-E_{2})\Leftrightarrow\langle\,p_{1},l;f/i\,|\,p_{2},l;f/i\,\rangle=\delta(p_{1}-p_{2}).$
(29)

We expand the unitary transformation to second order in $\alpha$:

$U_{if}=e^{-i\Theta}\cong
1-i\alpha\,\Theta^{(1)}-\frac{i}{2}\alpha^{2}\Theta^{(2)}-\frac{1}{2}\alpha^{2}\Theta^{(1)2},$
(30)

where

$\Theta\cong\alpha\,\Theta^{(1)}+\frac{1}{2}\alpha^{2}\,\Theta^{(2)}.$ (31)

Then we solve the equation

$H_{i}+V=U_{if}\,H_{f}\,U_{if}^{\dagger}$ (32)

order by order, with $V=\mathcal{O}(\alpha)$. The results are

$i[H_{f},\alpha\,\Theta^{(1)}]=V\quad\mathrm{and}\quad[H_{f},\alpha^{2}\,\Theta^{(2)}]=[\alpha\,\Theta^{(1)},V].$
(33)

The solutions for the matrix elements of $\Theta^{(1)}$ and $\Theta^{(2)}$ are

$\langle\,k_{1},l;f\,|\,\alpha\,\Theta^{(1)}\,|\,k_{2},l;f\,\rangle=-iP\frac{\mathcal{V}_{l}(k_{1},k_{2})}{E(k_{1})-E(k_{2})}$
(34)

and

$\langle\,k_{1},l;f\,|\,\alpha^{2}\,\Theta^{(2)}\,|\,k_{2},l;f\,\rangle=-iP\,\frac{1}{E(k_{1})-E(k_{2})}\int_{0}^{\infty}dk^{\prime}\,\left\\{\frac{\mathcal{V}_{l}(k_{1},k^{\prime})\mathcal{V}_{l}(k^{\prime},k_{2})}{E(k_{1})-E(k^{\prime})}-\frac{\mathcal{V}_{l}(k_{1},k^{\prime})\mathcal{V}_{l}(k^{\prime},k_{2})}{E(k^{\prime})-E(k_{2})}\right\\},$
(35)

where

$\mathcal{V}_{l}(k_{1},k_{2})=\langle\,k_{1},l;f\,|\,V\,|\,k_{2},l;f\,\rangle.$
(36)

Also

$\langle\,k_{1},l;f\,|\,\alpha^{2}\,\Theta^{(1)2}\,|\,k_{2},l;f\,\rangle=-P\int_{0}^{\infty}dk^{\prime}\,\frac{\mathcal{V}_{l}(k_{1},k^{\prime})}{E(k_{1})-E(k^{\prime})}\frac{\mathcal{V}_{l}(k^{\prime},k_{2})}{E(k^{\prime})-E(k_{2})}.$
(37)

The symbol $P$ indicates that when these expressions are integrated over
momentum, the principal part prescription (also known as finding the Cauchy
principal value) must be used. We have seen in (Hoffmann [3]) that use of this
prescription eliminates divergences and give results in agreement with an
exactly solvable model.

Then the interacting position separation wavefunction to second order is

$\displaystyle y_{l}^{(2)}(r,p)$
$\displaystyle=\langle\,r,l;f\,|\,U_{if}^{(2)}\,|\,p,l;f\,\rangle$
$\displaystyle=y_{l}^{(f)}(r,p)-P\int_{0}^{\infty}dk\,y_{l}^{(f)}(r,k)\,\frac{\mathcal{V}_{l}(k,p)}{E(k)-E(p)}+\frac{1}{2}P\int_{0}^{\infty}dk\,y_{l}^{(f)}(r,k)\int_{0}^{\infty}dk^{\prime}\,\frac{\mathcal{V}_{l}(k,k^{\prime})\mathcal{V}_{l}(k^{\prime},p)}{(E(k)-E(k^{\prime}))(E(k^{\prime})-E(p))}$
$\displaystyle\quad-\frac{1}{2}P\int_{0}^{\infty}dk\,y_{l}^{(f)}(r,k)\int_{0}^{\infty}dk^{\prime}\,\left\\{\frac{\mathcal{V}_{l}(k,k^{\prime})\mathcal{V}_{l}(k^{\prime},p)}{(E(k)-E(p))(E(k)-E(k^{\prime}))}-\frac{\mathcal{V}_{l}(k,k^{\prime})\mathcal{V}_{l}(k^{\prime},p)}{(E(k)-E(p))(E(k^{\prime})-E(p))}\right\\}.$
(38)

###  II.6 Calculation of phase shifts

For the Coulomb potential in the CM frame, the matrix elements are

$\displaystyle\mathcal{V}_{l}(k_{1},k_{2})$
$\displaystyle=\int_{0}^{\infty}dr\,y_{l}^{(f)}(r,k_{1})\,\frac{\alpha}{r}\,y_{l}^{(f)}(r,k_{2})$
$\displaystyle=\frac{\alpha}{\sqrt{\pi}}\,\frac{l!}{\Gamma(l+\frac{3}{2})}\,\begin{cases}\left(\rho\right)^{l+1}\,_{2}F_{1}(\frac{1}{2},l+1;l+\frac{3}{2};\rho^{2})&\rho<1,\\\
\left(\frac{1}{\rho}\right)^{l+1}\,_{2}F_{1}(\frac{1}{2},l+1;l+\frac{3}{2};\frac{1}{\rho^{2}})&\rho>1,\end{cases}$
(39)

with $\rho=k_{1}/k_{2}.$ Both sections of this function are singular at
$k_{1}=k_{2}.$ Close to the singularity, with $k_{1}=k_{2}(1+x)$ for small
$x,$ we find the approximation (Olver  _et al._ [7], their eq. (15.8.10))

$\mathcal{V}_{l}(k_{2}(1+x),k_{2})\rightarrow\frac{\alpha}{\pi}(C_{l}-\ln|x|)+\mathcal{O}(x),$
(40)

with

$C_{l}=\ln 2-\boldsymbol{C}-\psi(l+1)$ (41)

and

$\boldsymbol{C}=0.577215\dots$ (42)

known as Euler’s constant. Here $\psi(z)$ is the logarithmic derivative of the
gamma function (Gradsteyn and Ryzhik [8], their section 8.36).

These functions $\mathcal{V}_{l}(k_{1},k_{2})$ have the property that they
have a singular point but are integrable across the singularity. This will be
sufficient for the principal part integrals to converge to finite values.

To find the asymptotic form of $y_{l}^{(2)}(r,p)$ in eq. (38), we first use
the asymptotic forms of the free spherical waves (Messiah [6])

$y_{l}^{(f)}(r,k)\rightarrow\sqrt{\frac{2}{\pi}}\,\sin(kr-l\frac{\pi}{2}).$
(43)

For the first-order term (with $k=p(1+x)$)

$I^{(1)}\rightarrow-P\int_{-1}^{\infty}dx\,\sqrt{\frac{2}{\pi}}\,\\{\sin(pr-l\frac{\pi}{2})\cos(prx)+\cos(pr-l\frac{\pi}{2})\sin(prx)\\}\,\frac{\mathcal{V}_{l}(p(1+x),p)}{x}\left\\{\frac{E(p(1+x))+E(p)}{4p(2+x)}\right\\}.$
(44)

We review the rules we learned from the similar case in Hoffmann [3]:

  * •

For the $\cos(prx)$ and $\sin(prx)$ integrals, the contributions from $x>1$
will vanish asymptotically like $1/pr,$ since the integrand has no
singularities on that region and the sinusoids oscillate rapidly.

  * •

The factor $\sin(prx)/x\rightarrow\pi\,\delta(x)$ acts as a delta function.
When combined with the singularity of $\mathcal{V}_{l}(p(1+x),p),$ we have
(with $z=prx$)

$-\frac{\alpha}{\pi}\int_{-\infty}^{\infty}dz\,\mathrm{sinc}(z)(C_{l}-\ln|z/pr|)=\alpha(-C_{l}-\ln(pr)-\boldsymbol{C})=\alpha(-\ln(2pr)+\psi(l+1)).$
(45)

  * •

There is no delta function in the integral containing $\cos(prx).$ On
$x\in[-1,1],$ only integrands even in $x$ will survive the principal part
integration. In this integrand, there will be factors of the form
$[A]_{-}[B]_{+}/x+[A]_{+}[B]_{-}/x,$ with $A=\mathcal{V}_{l}(p(1+x),p)$ and
$B=\\{E(p(1+x))+E(p)\\}/p(2+x)$ and $[f(x)]_{\pm}=\frac{1}{2}\\{f(x)\pm
f(-x)\\}.$ We note that the first of these terms and $[B]_{-}/x$ are finite at
the origin. So we consider integrals of the form

$\int dx\,\ln x\,\cos(prx)=x\,\ln
x\,\mathrm{sinc}(prx)-\frac{\mathrm{Si}(prx)}{pr},$ (46)

where $\mathrm{Si}(x)=\int_{0}^{x}dt\,\mathrm{sinc}\,t.$ This integral
vanishes at $x=0$ and decreases like $1/pr$ for $x\neq 0.$

Thus we find (with $B(0)=E(p)/2p=1/2\beta,$ where
$\beta=p/\sqrt{p^{2}+m_{e}^{2}}$ is the relativistic velocity magnitude for
either of the particles)

$I^{(1)}=\sqrt{\frac{2}{\pi}}\,\cos(pr-l\frac{\pi}{2})\,\eta\,\\{-\ln(2pr)+\psi(l+1)\\}\equiv\sqrt{\frac{2}{\pi}}\,\cos(pr-l\frac{\pi}{2})\,\\{\delta_{l}^{(1)}(p)\\},$
(47)

with

$\delta_{l}^{(1)}(p)=-\eta\ln(2pr)+\eta\,\psi(l+1).$ (48)

Here $\eta=\alpha/2\beta$ corresponds to the parameter defined in the
nonrelativistic case, $\eta_{\mathrm{NR}}=\alpha/(p/\mu),$ with $\mu=m_{e}/2$
the reduced mass in this case. We will find that our perturbation series is in
powers of $\eta$ rather than of $\alpha.$

We use similar methods for the second-order terms, employing the four rules

$\displaystyle\int_{0}^{\infty}dk\,\frac{\sin(kr-l\frac{\pi}{2})}{E(k)-E(p)}\,f(k)$
$\displaystyle\rightarrow\cos(pr-l\frac{\pi}{2})\,\frac{\pi}{2\beta}\,f(p),$
$\displaystyle\int_{0}^{\infty}dx\,\frac{\cos(kr-l\frac{\pi}{2})}{E(k)-E(p)}\,f(k)$
$\displaystyle\rightarrow-\sin(pr-l\frac{\pi}{2})\,\frac{\pi}{2\beta}\,f(p),$
$\displaystyle\int_{-1}^{\infty}dx\,\frac{\sin(kr-l\frac{\pi}{2})}{E(k)-E(p)}\,\mathcal{V}_{l}(p(1+x),p)\,f(x)$
$\displaystyle\rightarrow\cos(pr-l\frac{\pi}{2})\,\delta_{l}^{(1)}(p)\,f(0)$
$\displaystyle\int_{-1}^{\infty}dx\,\frac{\cos(kr-l\frac{\pi}{2})}{E(k)-E(p)}\,\mathcal{V}_{l}(p(1+x),p)\,f(x)$
$\displaystyle\rightarrow-\sin(pr-l\frac{\pi}{2})\,\delta_{l}^{(1)}(p)\,f(0),$
(49)

all for $f(k)$ having no singularities. We evaluate the terms in eq. (38) as
written, performing each $k^{\prime}$ integral before the $k$ integral. We
find that changing the order of the double momentum integrals changes some of
the terms, but the total is invariant. An ambiguity remains. Using partial
fractions to cancel some terms gives a different final result. We suspect this
is because doing this changes the set of poles in the integrands. This matter
is under further investigation.

We find the result

$y_{l}^{(2)}(r,p)\rightarrow\sqrt{\frac{2}{\pi}}\,\sin(pr-l\frac{\pi}{2})\,\left\\{-\frac{1}{2}\delta_{l}^{(1)}(p)^{2}\right\\}+\sqrt{\frac{2}{\pi}}\,\cos(pr-l\frac{\pi}{2})\,\left\\{\delta_{l}^{(1)}(p)+\frac{1}{2}\,\frac{\pi}{2\beta}\,P\int_{0}^{\infty}dk^{\prime}\,\frac{\mathcal{V}_{l}(k^{\prime},p)^{2}}{E(k^{\prime})-E(p)}\right\\}.$
(50)

We note that the term proportional to $\sin(pr-l\frac{\pi}{2})$ has the
correct form required by unitarity at this order. The last term is

$\delta_{l}^{(2)}(p)=\frac{1}{2}\,\frac{\pi}{2\beta}\,P\int_{0}^{\infty}dk^{\prime}\,\frac{\mathcal{V}_{l}(k^{\prime},p)^{2}}{E(k^{\prime})-E(p)}=\frac{1}{2}\,\frac{\pi}{(2\beta)^{2}}\,P\int_{-1}^{\infty}dx\,\frac{\mathcal{V}_{l}(p(1+x),p)^{2}}{x}\,g(p,x),$
(51)

where

$g(p,x)=2\beta\,\frac{E(p(1+x))+E(p)}{4p(2+x)}\rightarrow 1$ (52)

at $x=0.$

##  III Calculation of differential cross sections

We use the results of Hoffmann [9], where we considered the scattering of a
single-particle wavepacket from a Coulomb potential. A slight modification of
the differential cross section formula, here in the CM frame, with the
resulting form shown here, was necessary to treat the scattering of two
particles,

$\frac{d\sigma}{d\Omega}=\frac{1}{4p^{2}}\left|\sum_{l=0}^{\infty}(2l+1)\,e^{-\epsilon^{2}(l+\frac{1}{2})^{2}}\,e^{i2(\bar{\delta}_{l}^{(1)}(p)+\delta_{l}^{(2)}(p))}\,e^{-(\delta-\Delta_{l}(p))^{2}/4}\,P_{l}(\cos\theta)\right|^{2}.$
(53)

The wavepacket treatment makes the sum in eq. (53) converge, even for
$\theta=0.$ It was previously thought that partial wave analysis was not
applicable to the Coulomb potential, as the sum over $l$ in eq. (53) diverges
for all scattering angles if the convergence factor
$\exp(-\epsilon^{2}(l+\frac{1}{2})^{2})$ is not present. This treatment
introduces another physical parameter,

$\epsilon=\frac{\sigma_{p}}{p},$ (54)

where $\sigma_{p}$ is the standard deviation of the Gaussian probability
density in momentum space of all the initial and final wavepackets. The
standard deviation in position space at $t=0$ is $\sigma_{x},$ with
$\sigma_{x}\sigma_{p}=\frac{1}{2}.$ In practice, we choose $\epsilon=0.001.$
In the nonrelativistic case, we found generally excellent agreement with the
Rutherford formula (Rutherford [10]),

$\frac{d\sigma}{d\Omega}_{\mathrm{Rutherford}}=\frac{m_{e}^{2}\alpha^{2}}{4p^{4}\sin^{4}(\frac{\theta}{2})},$

but deviations close to $\theta=0$ that are influenced by the choice of
$\epsilon.$

In eq. (53), the Coulomb phase shifts, $\sigma_{l}(p),$ have been replaced by
the phase shifts of our model to second order. Here

$\delta=\frac{\beta T-2R}{\sigma_{x}}$ (55)

measures the time delay (positive) or advancement (negative) over the time,
$T,$ of the experiment, relative to free evolution. The wavepacket centres are
separated by $2R$ at the start of the experiment. As discussed in Hoffmann
[9], the choice $R/\sigma_{x}=\sqrt{\epsilon}$ gives an initial separation
that grows as $\epsilon$ is decreased, while wavepacket spreading is minimal
over the course of the experiment. The quantities

$\Delta_{l}(p)=4\epsilon\eta\,\\{\ln(4pR/e)-\psi(l+1)+\mathcal{O}(\eta)\\},$
(56)

that influence this time shift involve the logarithmic phase, $-\eta\ln(2pr),$
from $\delta_{l}^{(1)}(p),$ which has been removed in the definition

$\bar{\delta}_{l}^{(1)}(p)=\delta_{l}^{(1)}(p)+\eta\ln(2pr)=\eta\,\psi(l+1).$
(57)

The approximation in eq. (56) will be sufficient for our purposes, as we will
mostly be considering small $\eta.$

In what follows, we evaluate the phase shifts numerically and insert them into
the cross section formula, sum the series numerically and take the modulus-
squared.

##  IV Numerical results

First we calculate the phase shifts and compare them to the Coulomb phase
shifts in the nonrelativistic regime ($p=0.02\,\mathrm{MeV}$, $\beta=0.04$,
$\eta=0.09$), where we expect agreement. We find the results plotted in fig.
1. The second-order corrections are seen to be much smaller than the Coulomb
phase shifts. Note that the exact and approximate phase shifts are negative at
$l=0$, so are not included in this logarithmic plot.

That the second-order corrections were calculable and finite achieves the
first goal of this paper. That the corrections were small is as it should be,
since the first-order result gave a good approximation to the Coulomb phase
shifts.

Figure 1: Comparison of exact and computed phase shifts.

As could be expected from this result, the plots of the model differential
cross section and the Rutherford formula were visually indistinguishable,
except close to $\theta=0,$ for this nonrelativistic choice of momentum. We
took $\delta_{l}^{(2)}(p),$ small but slowly varying, up to $l=50.$ The
profile in $\delta$ at $\theta=90\text{\textdegree}$ peaked very near
$\delta=0.$

Into the relativistic regime, we plotted the model differential cross section
(just to first order) at $\theta=90\text{\textdegree}$ for a range of momenta
from about $0.03\,\mathrm{MeV}$ to $19\,\mathrm{MeV}$ in fig. (2), compared
with the Rutherford differential cross section.

Figure 2: Model differential cross section at $\theta=90\text{\textdegree}$
for a range of nonrelativistic and relativistic momenta.

We see agreement with the Rutherford differential cross section at low momenta
and a distinct deviation at higher momenta. One reason that the relativistic
result would be larger than the Rutherford result is that the coupling
strength in the former can be no lower than $\eta>\alpha/2,$ since $\beta<1,$
while the nonrelativistic form $\eta_{\mathrm{NR}}=\alpha m_{e}/p$ has a lower
bound of zero.

At $p=5\,\mathrm{MeV}$ ($\eta=0.004$) we plotted the differential cross
section as a function of angle, again comparing with the Rutherford result.
This is shown in fig. (3). At $\theta=90\text{\textdegree},$ the $\delta$
profile peaked very near $\delta=0.$ We note that the model cross section was
seen to take a large but finite value at $\theta=0,$ while the Rutherford
formula diverges there. This finiteness was discussed in Hoffmann [9].

Figure 3: Angular dependence at $p=5\,\mathrm{MeV}.$

At higher momenta, a more meaningful comparison is with the Møller formula
(Itzykson and Zuber [11]). This is the differential cross section derived from
tree-level quantum electrodynamics for the scattering of two electrons. To
imitate a spinless system, it has been averaged over initial spins and summed
over final spins. The formula is

$\frac{d\sigma}{d\Omega}_{\mathrm{M\not{o}ller}}=\frac{1}{m_{e}^{2}}\frac{\alpha^{2}(1+\beta^{2})(1-\beta^{4})}{4\beta^{4}}\left[\frac{4}{\sin^{4}\theta}-\frac{3}{\sin^{2}\theta}\right]+\frac{1}{m_{e}^{2}}\frac{\alpha^{2}(1-\beta^{2})}{4}\left[1+\frac{4}{\sin^{2}\theta}\right].$

To make the comparison, we could not antisymmetrize these two fictitious
bosons. Instead, to have a definite particle exchange symmetry, we symmetrized
our model scattering amplitude, treating the two particles as identical. This
involved inserting a factor $1+(-)^{l}$ in the sum over $l$ that defines the
scattering amplitude.

For $p=5\,\mathrm{MeV},$ we plotted the angular dependence of the model cross
section along with that of the Møller formula and the Rutherford formula, in
fig. (4).

Figure 4: Angular dependence at $p=5\,\mathrm{MeV},$ compared with the Møller
formula.

For the scattering angle $\theta=90\text{\textdegree},$ we plotted the
dependence on momentum of the model cross section and the Møller formula, in
fig. (5).

Figure 5: Momentum dependence at $\theta=90\text{\textdegree},$ compared with
the Møller formula.

It is remarkable that our model lacks many of the essential features of the
electron-electron scattering problem, other than the Coulomb potential, yet
shows quite good agreement with the Møller formula. Note that the symmetrized
model cross section no longer converges to the Rutherford formula at
nonrelativistic energies. Instead, it must converge to a symmetric cross
section that evidently has some features in common with the Rutherford
formula. We recall that the unsymmetrized cross section was seen to converge
to the unsymmetrized Rutherford formula. It is then puzzling why the Møller
formula converges to a cross section that has no definite particle exchange
symmetry.

##  V Conclusions

There were two aims of this paper. One was to use a novel perturbation theory
on the Coulomb potential and show that finite results could be obtained. In
comparison, other forms of perturbation theory find infinite results at second
and all higher orders. This aim was achieved using only one rule that was
imposed rather than derived: the requirement that principal part integration
be used on the momentum integrals.

The other aim was to show that it is possible to construct a theory that is
relativistic and quantum-mechanical, but not quantum field-theoretic and not
the wave mechanics of a single particle in a potential. This was achieved by
working in momentum space, where the relativistic dependence of single-
particle energy on momentum can be easily introduced. We found that the three-
vector part of the four-potential must vanish in the centre of mass frame.
That leaves an arbitrary zero component, at first sight not constrained by
relativity. The further success of this method depends on being able to choose
a CM-frame potential that reproduces experimental results.

The model of a Coulomb potential in the centre of mass frame was not intended
to provide a full description of the scattering of two electrons. It was
intended only as a proof of principle for the method considered. There are
other terms that contribute to the scattering of two electrons. In a future
work, we will include those terms, with the expectation of finite results.

Also in a future work, we will analyze the forms of the terms in our
perturbation theory, to see if it is possible to prove finiteness at higher
orders.

## References

  * Dalitz [1951] R. H. Dalitz, Proc. Roy. Soc. Lond. A: Math., Phys. and Eng. Sci. 206, 509 (1951).
  * Collas [2021] P. Collas, arXiv:2102.13105 (2021).
  * Hoffmann [2021] S. E. Hoffmann, J. Math. Phys. 62, 032105 (2021).
  * Foldy and Wouthuysen [1950] L. L. Foldy and S. A. Wouthuysen, Phys. Rev. 78, 29 (1950).
  * Weinberg [2012] S. Weinberg,  _The Quantum Theory of Fields_ , Vol. 1 (Cambridge University Press, N.Y., 2012).
  * Messiah [1961] A. Messiah,  _Quantum Mechanics_ , Vol. 1 and 2 (North-Holland, Amsterdam and John Wiley and Sons, N.Y., 1961).
  * Olver  _et al._ [2020] F. W. J. Olver _et al._ , “NIST Digital Library of Mathematical Functions,” http://dlmf.nist.gov/, Release 1.0.25 of 2019-12-15 (2020).
  * Gradsteyn and Ryzhik [1980] I. S. Gradsteyn and I. M. Ryzhik,  _Tables of Integrals, Series and Products_ , corrected and enlarged ed. (Academic Press, Inc., San Diego, CA, 1980).
  * Hoffmann [2017] S. E. Hoffmann, J. Phys. B: At. Mol. Opt. Phys. 50, 215302 (2017).
  * Rutherford [1911] E. Rutherford, Phil. Mag. Series 6, vol. 21, 669 (1911).
  * Itzykson and Zuber [1980] C. Itzykson and J.-B. Zuber,  _Quantum Field Theory_ , 1st ed. (McGraw-Hill Inc., 1980).