# Relativistic levels of mesic atoms

R. Giachetti  Physics Department, Università di Firenze, Italy
Istituto Nazionale di Fisica Nucleare, Sezione di Firenze, Italy  E. Sorace
Istituto Nazionale di Fisica Nucleare, Sezione di Firenze, Italy

###### Abstract

We revisit the derivation of the covariant two-body scalar-fermion equation
with a Coulomb interaction, presented in a previous paper. We show that it can
be given the formal aspect of a Dirac equation, but for the fact that the
eigenvalue is also contained in one of the coefficients and thus it is not
linearly included. The discussion of the boundary value problem is therefore
different, although some properties of the Dirac equation can be recovered in
an approximation bringing back to the concept of reduced mass. We discuss a
mixed analytic-numerical method of solution which allows to obtain very
accurate results and we calculate the lowest levels, states and QED
corrections for pionic and kaonic atoms.

PACS numbers: 12.39.Pn, 03.65.Ge, 03.65.Pm, 36.10.-k
Keywords Covariance; Two-body; Wave-equation; Pion; Kaon.

##  1 Introduction

Some years ago we began a research program on the two body covariant wave
equations and we developed a relativistic wave equation for two interacting
fermions [1, 2, 3, 4, 5]. We studied the covariance, the different limits and
the spectral properties of such an equation and the first applications were
made on the hyperfine levels of the Hydrogenic atoms. We then addressed to
high energy physics, calculating the mass of the mesons and their radiative
decays. In all cases the results were in excellent agreement with experimental
results, even when dealing with the masses of light mesons for which potential
models had generally failed. The essential reason for the accuracy of the
results can be traced back to the complete covariance of the treatment, which
includes all the relativistic and the two-body effects. In more recent times,
in order to broaden the applicability of the method, we have extended our
analysis by considering relativistic objects of different nature. We have thus
determined the wave equations for interacting scalar-scalar and scalar-fermion
particles [6], enabling thus the study of further elementary two-body systems.
Again we have proved the covariance, the different limits of such equations
and the general properties of their spectrum, giving them a clear physical
interpretation.

In this paper we revisit and we give a further development of the relativistic
scalar-fermion equation together with its application to the pionic and kaonic
atoms. The physical interest for such systems is due to the fact that meson
masses are two-three order of magnitude larger than the electron mass, their
atomic orbits are of the order of hundred fm and their binding energies of
some keV. Therefore they provide a particularly useful opportunity for
studying nuclear interactions in QCD low energy regime [7, 8, 9]. The effects
of the strong interactions on the spectrum is given by a shift of the purely
electromagnetic levels and an absorption width, both of them experimentally
measured through X-ray spectroscopy [10]. In fact the size of the atom is much
larger than the range of strong interactions, so that QED itself provides
almost exact values of the atomic levels and the fundamental state is actually
the only one for which the nuclear force plays a role. Hence the difference
between the measured transition energies from excited states to the
fundamental one and the corresponding energies calculated by pure QED are
essentially the same as the shift of the fundamental level. A history of the
determination of the pion-nucleon coupling constant can be found in [11]. Up
to the most recent papers [12, 13], the mass of the $\pi$ meson has been
evaluated by matching the measured values of the photon energy emitted in the
transitions between two rather high level of pionic atoms. The energy of those
same decays were calculated by the Schrödinger equation with a perturbation
treatment of the relativistic and QED effects [14]. Different corrections of a
more phenomenological nature have also been considered, such as, for instance,
the finite dimensions of the particles, usually treated as if they were
uniformly charged spheres with a non vanishing radius [15, 16, 17]. A
different measure involving the $\mu$-neutrino mass has been presented in
[18].

We aim here at a more essential treatment which maintains the covariance, such
to take exactly into account the relativistic effects. We therefore will add
the lowest order QED corrections, due to vacuum polarization, while assuming a
basic model formed by two particles which are point-like and interact by a
Coulomb potential. Of course this does not mean that we believe that
phenomenological corrections are negligible, but only that we defer their
inclusion until we find a formulation fitting the covariant two body
framework. For instance, we know that it is comparatively easy to introduce a
term mimicking a finite radius of an attracting center for the case of a
single particle in an external potential. For two interacting points, after
the separation of the global coordinate, the relative part becomes an
entangled system and the definition of a potential producing a finite radius
effect in terms of the relative coordinate only seems tough. Starting with
extended particles would directly lead to the formulation of a field theory
producing the appropriate form factors. This is outside of our present
purposes. On the other hand we notice that the final reduced form of the
scalar-fermion equation develops a radial term clearly related to relativistic
corrections of the the classical reduced mass. Thus, besides the results, we
also take the opportunity to illustrate some properties of the fermion-scalar
equation, as well as to explain the ideas of the mixed analytical-numerical
treatment we have used in order to obtain a very accurate solution of the
spectral problem. High precision appears to be necessary, as seen by the
numerical results, since the spectrum presents pairs of very close levels,
whose difference becomes fainter and fainter with increasing energy. These
states appear as a direct consequence of the covariant treatment which takes
into account the fermionic nature of the proton. For the sake of completeness,
we begin Section 2 with a sketchy derivation of the scalar-fermion equation
relegating to Appendix I the summary of the kinematical variables introduced
in our previous papers. We then show that the scalar-fermion equation can be
cast into the the form of a Dirac equation in a central potential, but for the
fact that the reduction procedure makes a coefficient depending upon the
eigenvalue. Unfortunately this non-linear inclusion of the eigenvalue forbids
the application of most results holding for the Dirac equation [19] and
produces, in particular, the close pairs of levels we mentioned above. In
Section 3 some details of the discussion of the boundary value problem are
given. In particular we work out the series solutions and the recurrence
relations at the two singular points, the origin and infinity, as well as at
any regular point of the positive real line. We then present the results for
the pure Coulomb relativistic levels of the pionic and kaonic atoms. In
Section 4 we state the form of the potentials we have used for the one- and
two-loop vacuum polarization potentials. The detailed results, obtained with
the use of our eigenfunctions, are presented in the table given in Appendix
II. It can be seen that the one-loop and irreducible two-loop corrections are
very regular, namely monotonically decreasing for increasing energy of the
levels. The apparently more bizarre behavior of the corrections due to
reducible two-loop graph comes from a change of sign of the corresponding
potential. As a consequence, the overall results are very sensitive to the
wave-function density distribution and a little bit longer to be calculated
with the desired accuracy. In particular, for some levels, the different
ratios of the meson to proton mass of pionic and kaonic atoms may produce
effects with opposite sign. In any case these last corrections are one-two
orders less than the previous ones. A brief comment is also in order about the
seven decimal figures which we have presented the results with. All of them
are numerically meaningful, although supernumerary for a comparison with the
experimental results. Their inclusion is motivated in order to give a clear
distinction of the levels and to show the regular trend of the behavior both
of spectrum and corrections. Some brief conclusions are given in the last
Section 5. In the final Appendix III we show an intriguing feature of our
scalar-fermion equation, in which it is possible to single out an interaction
term related to the corrections connected with the classical ’reduced mass’.
Indeed, although not directly deduced from the covariant framework, this
quantity gives an approximation of the levels to the first order of the fine
structure constant. Relativistic corrections to the reduced mass improve the
accuracy of the levels [20].

##  2 The scalar-fermion equation and its spectrum

In this section we recall the basic steps leading to the relativistic wave
equation for a fermion and a scalar particle interacting through a Coulomb
potential [6]. We then make explicit the dimensionless form of the second
order system which defines the boundary value problem leading to the spectrum.
An improvement with respect to the original derivation is the Dirac form the
reduced equation can be given, which makes somewhat more familiar the
treatment. The solution has been discussed analytically as far as possible in
order to obtain a high precision for the spectral levels. In the next Section
we will therefore develop the recurrence relations which make it possible a
rapid achievement of very accurate numerical results. As we said in the
Introduction, the physical systems we have considered are the pionic atom
$\,({H}^{+}\\!,\,\pi^{-})\,$ and the kaonic atom $\,({H}^{+}\\!,\,{K}^{-}),$
whose lowest Coulomb levels are given in Tables 1 and 2.

Let us recall some elementary properties of the Klein-Gordon (KG) and Dirac
equations in order to establish the two-body relativistic equation for a
scalar and a fermion. The covariance of the method has been largely discussed
in the series of our previous papers [1, 2, 3, 4, 5, 6], so that, for the sake
of completeness, we just report in Appendix I the definition of the phase-
space we are using and some few words to describe its physical and geometrical
properties.

A Hamiltonian formulation for the Klein-Gordon equation can be given in the
form of a 2-dim system:

$\displaystyle i\,{\partial\Phi}/{\partial
t}-H_{S}\,\Phi=0\,,\quad\Phi={}^{T}\bigl{(}\phi_{1},\phi_{2}\bigr{)}$ (2.1)

In the standard notation $\,\sigma_{z},\,\sigma_{\pm}\,$ for the Pauli
matrices, the Hamiltonian $H_{S}$ has the expression

$\displaystyle
H_{S}=-\,\sigma_{-}\,\,\bigl{(}\,{{\nabla}^{2}}/{2m_{S}}\,\bigr{)}+(\sigma_{z}+2\sigma_{+})\,m_{S}=\left(\begin{matrix}m_{S}&2m_{S}\vrule
height=1.99997pt,width=0.0pt,depth=1.99997pt\cr-{{\nabla}^{2}}/{2m_{S}}&-m_{S}\end{matrix}\,\right)$
(2.2)

where $\,m_{S}\,$ is the mass of the scalar. It is obtained from the Feshbach-
Villars [21] representation of the KG equation with a transformation generated
by $\sigma_{+}+\sigma_{z}$.

The free Dirac Hamiltonian $H_{F}$ in spherical coordinates reads

$\displaystyle H_{F}=\left(\begin{matrix}m_{F}&0&q_{0}&\sqrt{2}\,q_{-}\vrule
height=1.99997pt,width=0.0pt,depth=1.99997pt\\\
0&m_{F}&-\sqrt{2}\,q_{+}&q_{0}\vrule
height=1.99997pt,width=0.0pt,depth=1.99997pt\\\
q_{0}&\sqrt{2}\,q_{-}&-m_{F}&0\vrule
height=1.99997pt,width=0.0pt,depth=1.99997pt\\\
-\sqrt{2}\,q_{+}&-q_{0}&0&-m_{F}\end{matrix}\right)$ (2.3)

where $\,m_{F}\,$ is the mass of the fermion and the spherical derivatives
$q_{+}$, $q_{-}$ and $q_{0}$ are defined by

$\displaystyle q_{\pm}=-(\,\pm{q_{x}}+i\,{q_{y}}\,)/\sqrt{2}\,,\qquad\qquad
q_{0}={q_{z}}\,,\qquad\qquad\,q_{k}\rightarrow-i\partial/\partial r_{k}\,$
(2.4)

Denoting the fermion radial coordinate by $\xi=(x^{2}+y^{2}+z^{2})^{1/2}$, the
use of the spherical spinors

$\Omega_{\,\ell+\frac{1}{2},\,\ell,\,m}(\theta,\phi)=\begin{pmatrix}\bigl{(}{j+m}\bigr{)}^{1/2}\,\bigl{(}{2j}\bigr{)}^{-1/2}\,Y_{\ell,\,m-\frac{1}{2}}(\theta,\phi)\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
\bigl{(}{j-m}\bigr{)}^{1/2}\,\bigl{(}{2j}\bigr{)}^{-1/2}\,Y_{\ell,\,m+\frac{1}{2}}(\theta,\phi)\end{pmatrix},~{}~{}\Omega_{\,\ell-\frac{1}{2},\,\ell,\,m}(\theta,\phi)=\begin{pmatrix}-\bigl{(}{j-m+1}\bigr{)}^{1/2}\,\bigl{(}{2j+2}\bigr{)}^{-1/2}\,Y_{\ell,\,m-\frac{1}{2}}(\theta,\phi)\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
\phantom{-}\bigl{(}{j+m+1}\bigr{)}^{1/2}\,\bigl{(}{2j+2}\bigr{)}^{-1/2}\,Y_{\ell,\,m+\frac{1}{2}}(\theta,\phi)\end{pmatrix}$
(2.5)

allows for the construction of a Dirac state where angular momentum and parity
are diagonal:

$\displaystyle\psi(\xi,\theta,\phi)=\begin{pmatrix}\varphi(\xi,\theta,\phi)\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
\chi(\xi,\theta,\phi)\end{pmatrix}=\begin{pmatrix}a(r)\,\,\Omega_{j\ell
m}(\theta,\phi)\vrule height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
b(r)\,\,\Omega_{j\ell^{\prime}m}(\theta,\phi)\end{pmatrix}\,,\qquad\ell=j-\frac{1}{2},\quad\ell^{\prime}=j+\frac{1}{2}$
(2.6)

We now write the relativistic Hamiltonian for a scalar and a fermion
interacting through a Coulomb potential expressed in the coordinates
introduced in Appendix I:

$\displaystyle H=H_{S}\otimes\mathbf{I}_{4}+\mathbf{I}_{2}\otimes
H_{F}+\mathbf{I}_{4}\,V(r)\,,\qquad V(r)=-\alpha/r$ (2.7)

It generates the stationary eigenvalue equation

$\displaystyle H\,\Phi=M\,\Phi\,,\qquad M=m_{F}+m_{S}+E$ (2.8)

where the states of opposite parity have the form

$\displaystyle\Phi_{\mathrm{I\phantom{I}}}(\bm{r})=\begin{pmatrix}a_{1}(r)\,\,\Omega_{j\ell
m}(\theta,\phi)\vrule height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
a_{2}(r)\,\,\Omega_{j\ell^{\prime}m}(\theta,\phi)\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\ a_{3}(r)\,\,\Omega_{j\ell
m}(\theta,\phi),\vrule height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
a_{4}(r)\,\,\Omega_{j\ell^{\prime}m}(\theta,\phi)\end{pmatrix}\,,\qquad\Phi_{\mathrm{II}}(\bm{r})=\begin{pmatrix}a_{1}(r)\,\,\Omega_{j\ell^{\prime}m}(\theta,\phi)\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\ a_{2}(r)\,\,\Omega_{j\ell
m}(\theta,\phi)\vrule height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
a_{3}(r)\,\,\Omega_{j\ell^{\prime}m}(\theta,\phi)\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\ a_{4}(r)\,\,\Omega_{j\ell
m}(\theta,\phi)\end{pmatrix}$ (2.9)

and $a_{i}(r)$ are coefficients depending upon the relative radial variable
only. We redefine two of the unknown functions by letting
$\,\,a_{2}(r)\rightarrow-ia_{2}(r)\,,\quad$ $\,\,a_{4}(r)\rightarrow-
ia_{4}(r)\,\,~{}$ and we introduce the dimensionless variables

$\displaystyle\bm{s}=m_{F}\,\bm{r},\quad~{}\sigma=m_{S}/m_{F},\quad~{}\epsilon=E/m_{F},\quad~{}\lambda=M/m_{F}=1+\sigma+\epsilon,\quad$
(2.10)

By substituting (2.9) into (2.7) we get a system of eight equations equal in
pairs, and the final system to be solved is

$\displaystyle\biggl{(}\,\frac{d}{ds}-\frac{j-1/2}{s}\,\biggr{)}\,a_{1}(s)+\Bigl{(}\,\lambda+\frac{\alpha}{s}-\sigma+P\,\Bigr{)}\,a_{2}(s)-2\sigma\,a_{4}(s)=0\vrule
height=11.99998pt,width=0.0pt,depth=11.99998pt$ (2.11)
$\displaystyle\biggl{(}\,\frac{d}{ds}+\frac{j+3/2}{s}\,\biggr{)}\,a_{2}(s)-\Bigl{(}\,\lambda+\frac{\alpha}{s}-\sigma-P\,\Bigr{)}\,a_{1}(s)+2\sigma\,a_{3}(s)=0\vrule
height=11.99998pt,width=0.0pt,depth=11.99998pt$ (2.12)
$\displaystyle\biggl{(}\,\frac{d}{ds}+\frac{j+3/2}{s}\,\biggr{)}\,a_{4}(s)-\Bigl{(}\,\lambda+\frac{\alpha}{s}+\sigma-P\,\Bigr{)}\,a_{3}(s)-\frac{\nabla^{2}}{2\sigma}\,a_{1}(s)=0\vrule
height=11.99998pt,width=0.0pt,depth=11.99998pt$ (2.13)
$\displaystyle\biggl{(}\,\frac{d}{ds}-\frac{j-1/2}{s}\,\biggr{)}\,a_{3}(s)+\Bigl{(}\,\lambda+\frac{\alpha}{s}+\sigma+P\,\Bigr{)}\,a_{4}(s)+\frac{\nabla^{2}}{2\sigma}\,a_{2}(s)=0$
(2.14)

where $\,P=\pm 1\,$, the positive sign referring to the state
$\,\Phi_{\mathrm{I}}(\bm{s})\,$, the negative sign to
$\,\Phi_{\mathrm{II}}(\bm{s})\,$. Notice that the Laplace operator of the
third equation contains an angular momentum contribution
$r^{-2}\,(j-1/2)(j+1/2)$, while in the last equation the angular momentum
contribution is $r^{-2}\,(j+1/2)(j+3/2)$. Actually, from (2.14) it turns out
that $a_{3}(r)$ and $a_{4}(r)$ can be expressed in terms of $a_{1}(r)$ and
$a_{2}(r)$ by the algebraic relations, equal for both parities:

$a_{3}(s)=\frac{\Bigl{(}\bigl{(}(\lambda-\sigma)s+\alpha\bigr{)}^{2}-s^{2}\Bigr{)}\,a_{1}(s)-\alpha\,a_{2}(s)}{4s\,(\lambda
s+\alpha)\,\sigma},\qquad
a_{4}(s)=\frac{\alpha\,a_{1}(s)+\Bigl{(}\bigl{(}(\lambda-\sigma)s+\alpha\bigr{)}^{2}-s^{2}\Bigr{)}\,a_{2}(s)}{4s\,(\lambda
s+\alpha)\,\sigma}\,.$ (2.15)

Denoting for simplicity $f(s)=\,a_{1}(s)\,$, $g(s)=-\,a_{2}(s)\,$ and
introducing the parameters

$\displaystyle\Lambda=\frac{\lambda}{2}+\frac{1-\sigma^{2}}{2\,\lambda},\qquad
B=\frac{\alpha\,(1-\sigma^{2})}{2\,\lambda^{2}},\qquad\rho=\frac{\alpha}{\lambda}$
(2.16)

the differential problem to be solved reduces to the second order boundary
condition problem with eigenvalue $\epsilon$ (or $\lambda$)

$\displaystyle\frac{d}{ds}f(s)+\Bigl{(}\,-\frac{j}{s}+\frac{1}{2\,(s+\rho)}\,\Bigr{)}\,f(s)-\Bigl{(}\,\Lambda+1+\frac{\alpha}{2\,s}-\frac{B}{s+\rho}\,\Bigr{)}\,g(s)=0\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (2.17)
$\displaystyle\frac{d}{ds}g(s)+\Bigl{(}\,\frac{j+1}{s}+\frac{1}{\,2\,(s+\rho)}\,\Bigr{)}\,g(s)+\Bigl{(}\,\Lambda-1+\frac{\alpha}{2\,s}-\frac{B}{s+\rho}\,\Bigr{)}\,f(s)=0$
(2.18)

The system (2.18) can be given a more convenient form by substituting the
angular momentum $\,j\,$ with $\,\kappa=-j-1/2\,$ as in [19] and by changing
the unknown functions $\,f(s),\,g(s)\,$ with $\,\phi(s),\,\chi(s)\,$ according
to

$\displaystyle f(s)=\frac{\phi(s)}{\sqrt{s\,(s+\rho)}},\qquad
g(s)=\frac{\chi(s)}{\sqrt{s\,(s+\rho)}}.$ (2.19)

We get

$\displaystyle\frac{d}{ds}\phi(s)+\frac{\kappa}{s}\,\phi(s)-\Bigl{(}\,\Lambda+1+\frac{\alpha}{2\,s}-\frac{B}{s+\rho}\,\Bigr{)}\,\chi(s)=0\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (2.20)
$\displaystyle\frac{d}{ds}\chi(s)-\frac{\kappa}{s}\,\,\chi(s)+\Bigl{(}\,\Lambda-1+\frac{\alpha}{2\,s}-\frac{B}{s+\rho}\,\Bigr{)}\,\phi(s)=0$
(2.21)

The system (2.21) is formally identical to a Dirac equation in a central
potential [19]. However, besides the fact that for a given value of the
parameter $\,\Lambda\,$ there are two determinations of the eigenvalue
$\,\lambda\,$ we are concerned with, i.e.
$\,\lambda=\Lambda\pm\bigl{(}\Lambda^{2}+\sigma^{2}-1)^{1/2}\,$, the major
point that forbids the system (2.21) to be treated as a genuine Dirac equation
is the dependence of the parameter $\,B\,$ on $\,\lambda\,$ and hence on
$\,\Lambda\,$. In particular, for (2.21) the spectrum is no more dependent
only upon the absolute value of $\,\kappa\,$ as for the Dirac equation in a
central potential [19]. Thus the level degeneracy is removed and the states
with $\,\pm\kappa\,$ appear separately, although their difference is really
small, as we can see from the numerical results. Finally notice that an actual
Dirac equation in a Coulomb potential with coupling constant $\,\alpha/2\,$ is
obtained for equal fermion and scalar masses: in this case
$\,\sigma=1,~{}B=0\,$ and $\,\Lambda=\lambda/2,$.

For general masses, the solution of (2.21) cannot be completely analytical. It
is a boundary value problem with singularities at the origin and at infinity
and we use a double shooting framework for obtaining the eigenvalues and the
eigenfunctions. However, in order to keep a high standard of accuracy in the
results, we can push far enough the analytical development and determine the
exact recurrence relations for the approximate solutions needed at the two
singular points. Since the system is regular in $(0,\infty)$, we calculate two
more holomorphic series solutions, each of which centered sufficiently close
to one of the singular points. We then match them to the approximate solution
coming from the corresponding singular point. The continuity of components of
these two series at an arbitrarily chosen crossing point gives, as usual, the
spectral condition.

We conclude this section with an observation on the limits of (2.14) for large
values of the mass of each one of the two components. Since
$\lambda+\alpha/s-\sigma=1+\epsilon+\alpha/s$, when the mass of the scalar
tends to infinity, namely $\sigma\rightarrow\infty$, from (2.10) we see that
the last two equations of (2.14) are identically vanishing and the two first
ones give the Dirac equation [19]

$\displaystyle\Bigl{(}\,\frac{d}{ds}+\frac{1+\kappa}{s}\,\Bigr{)}\,a_{1}(s)+\Bigl{(}\,1+\epsilon+\frac{\alpha}{s}+1\,\Bigr{)}\,a_{2}(s)=0\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (2.22)
$\displaystyle\Bigl{(}\,\frac{d}{ds}+\frac{1-\kappa}{s}\,\Bigr{)}\,a_{2}(s)-\Bigl{(}\,1+\epsilon+\frac{\alpha}{s}-1\,\Bigr{)}\,a_{1}(s)=0$
(2.23)

When the mass of the fermion tends to infinity we recover the KG equation for
the scalar. This is easily seen by first making the change of variable
$\,{\bm{s}}=(1/\sigma)\,{\bm{u}}\,$ and taking the mass of the lighter
particle as mass scale. In order to normalize the physical parameters to the
mass of the scalar, we define $\,\eta=\epsilon/\sigma$. Developing
$\,a_{j}(u)=(1/\sigma)\,a_{j}^{(1)}(u)+a_{j}^{(0)}(u),\,j=1..4\,$, at the
order $\,\sigma^{-2}\,$ we find the vanishing of $\,a_{2}^{(1)}(u)\,$ and
$\,a_{4}^{(1)}(u)\,$. At the order $\,\sigma^{-1}\,$ we then see that
$\,a_{2}^{(0)}(u),\,a_{3}^{(1)}(u)\,$ and $\,a_{4}^{(0)}(u)\,$ can be
expressed in terms of $\,a_{1}^{(1)}(u)$, which, in turn, satisfies

$\nabla^{2}\,a_{1}^{(1)}(u)+\biggl{(}\,\Bigl{(}1+\eta+\frac{\alpha}{u}\Bigr{)}^{2}-1\,\biggr{)}\,a_{1}^{(1)}(u)=0$
(2.24)

namely the KG equation in a Coulomb field.

##  3 The solution of the boundary value problem.

In this section we give the details of our method for calculating the spectrum
of the system (2.18), together with the results of the pure Coulomb levels of
the $\,({H}^{+}\\!,\,\pi^{-})\,$ and $\,({H}^{+}\\!,\,{K}^{-})\,$ atoms. A
global exact analytical solution is not available. Therefore we look for a
result in the form of piecewise approximate analytical solutions expressed by
series and, in in order to accelerate the convergence, by Padé approximants.
This allows to obtain any required precision and to have an effective control
of it. Clearly, any improvement of the accuracy occurs to the detriment of the
computation time. However, numerical values of very large numbers of
coefficients of the series involved are efficiently and rapidly found by
recurrence procedures. Only the calculations of the Padé are somewhat longer,
still completely remaining within acceptable time limits even for a normal
workstation. Here below we give, thus, the explicit form of the recurrence
relations of (2.21) for the series solution at the origin, the asymptotic
solution at infinity and the holomorphic solution at any regular point of the
open positive line $(0,\infty)$.

### The solution at the origin.

For the solution at the origin we make the usual expansion of the unknown
functions $\phi(s)$ and $\chi(s)$:

$\displaystyle\phi(s)=s^{\nu}\,\sum_{n=0}^{\infty}\phi_{n}\,s^{n},\qquad\chi(s)=s^{\nu}\,\sum_{n=0}^{\infty}\chi_{n}\,s^{n}$
(3.1)

The index $\nu$ is obtained by the relations with $n=0$

$\displaystyle(\nu+\kappa)\,\phi_{0}-({\alpha}/2)\,\chi_{0}=0\vrule
height=6.00006pt,width=0.0pt,depth=6.00006pt$ (3.2)
$\displaystyle({\alpha}/2)\,\phi_{0}+(\nu-\kappa)\,\chi_{0}=0$ (3.3)

which actually are an eigenvalue system in $\nu$, yielding

$\displaystyle\nu=\sqrt{\kappa^{2}-\alpha^{2}/4},\qquad\phi_{0}=1,\qquad\chi_{0}=(2/{\alpha})\,(\nu+\kappa)$
(3.4)

By substituting (3.1) in the system (2.21), expanded in the neighborhood of
$s=0$, and letting

$\displaystyle\Sigma_{\phi,\,n}=\sum_{p=0}^{n}(-\rho)^{-p}\,\phi_{n-p},\qquad\Sigma_{\chi,\,n}=\sum_{p=0}^{n}(-\rho)^{-p}\,\chi_{n-p}$
(3.5)

we obtain the recurrence relations for $n\geq 1$. They read

$\displaystyle\begin{pmatrix}\nu+\kappa+n&-\alpha/2\vrule
height=8.00003pt,width=0.0pt,depth=8.00003pt\\\
\alpha/2&\nu-\kappa+n\end{pmatrix}\,\begin{pmatrix}\phi_{n}\vrule
height=8.00003pt,width=0.0pt,depth=8.00003pt\\\
\chi_{n}\end{pmatrix}=\begin{pmatrix}\phantom{-}(\Lambda+1)\,\chi_{n-1}-(B/\rho)\,\Sigma_{\chi,n-1}\vrule
height=8.00003pt,width=0.0pt,depth=8.00003pt\\\
-(\Lambda-1)\,\phi_{n-1}+(B/\rho)\,\Sigma_{\phi,n-1}\end{pmatrix}$ (3.6)

Solving in $\phi_{n}$ and $\chi_{n}$, we get the expressions to be implemented
by a numerical code:

$\displaystyle\phi_{n}={\displaystyle{\frac{1}{2n\,(n+2\nu)}}}\,\biggl{(}\,\,2\,(\nu+n-\kappa)\,\Bigl{(}\,(\Lambda+1)\,\chi_{n-1}-(B/\rho)\,\,\Sigma_{\chi,n-1}\,\Bigr{)}-\alpha\,\Bigl{(}\,(\Lambda-1)\,\phi_{n-1}-(B/\rho)\,\,\Sigma_{\phi,n-1}\,\Bigr{)}\,\biggr{)}\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (3.7)
$\displaystyle\chi_{n}={\displaystyle{\frac{1}{2n\,(n+2\nu)}}}\,\biggl{(}\,-\alpha\,\Bigl{(}\,(\Lambda+1)\,\chi_{n-1}-(B/\rho)\,\,\Sigma_{\chi,n-1}\,\Bigr{)}-(\nu+\kappa+n)\,\Bigl{(}\,(\Lambda-1)\,\phi_{n-1}-(B/\rho)\,\,\Sigma_{\phi,n-1}\,\Bigr{)}\,\biggr{)}$
(3.8)

The corresponding series have a finite radius of convergence $s_{0}$. A point
in $(0,s_{0})$ will be chosen to match this solution with the holomorphic
series described here below.

### The solution at a regular point.

We first shift the system (2.21) at a point $\,u_{0}>0\,$ by letting
$\,s=u+u_{0}\,$, so that $\,u=0\,$ is a regular point. We have

$\displaystyle\frac{d}{du}\phi(u)+\frac{\kappa}{u+u_{0}}\,\phi(u)-\Bigl{(}\,\Lambda+\frac{\alpha}{2\,(u+u_{0})}-\frac{B}{u+u_{0}+\rho}+1\,\Bigr{)}\,\chi(u)=0\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (3.9)
$\displaystyle\frac{d}{du}\chi(u)-\frac{\kappa}{u+u_{0}}\,\,\chi(u)+\Bigl{(}\,\Lambda+\frac{\alpha}{2\,(u+u_{0})}-\frac{B}{u+u_{0}+\rho}-1\,\Bigr{)}\,\phi(u)=0$
(3.10)

The expansion of the solutions at the regular point $\,u=0\,$ are now the
holomorphic series

$\displaystyle\phi(u)=\sum_{n=0}^{\infty}\phi_{n}\,u^{n},\qquad\chi(u)=\sum_{n=0}^{\infty}\chi_{n}\,u^{n}$
(3.11)

The coefficients $\,(\phi_{0},\chi_{0})\,$ are arbitrary integration constants
that will be chosen to be the pairs $\,(1,0)\,$ and $\,(0,1)\,$ in order to
obtain a set of fundamental solutions, necessary for the matching with the
expansions at the origin and at infinity. For $\,n\geq 1\,$, assuming that the
coefficients with negative index are vanishing, a substitution into (3.10)
yields the following expressions of $\phi_{n}$ and $\chi_{n}$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\phi_{n}={\displaystyle{\frac{1}{n\,(u_{0}+\rho)\,u_{0}}}}\,\biggl{(}\,-\Bigl{(}(\rho+2\,u_{0})\,(n-1)+(\rho+u_{0})\,\kappa\Bigr{)}\,\phi_{n-1}+\Bigl{(}(\Lambda+1)\,u_{0}^{2}+\Bigl{(}(\Lambda+1)\,\rho-B+\alpha/2\Bigr{)}\,u_{0}+\rho\,\alpha/2\Bigr{)}\,\chi_{n-1}$
(3.12)
$\displaystyle\phantom{\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\phi_{n}={\displaystyle{\frac{1}{n\,(u_{0}+\rho)\,u_{0}}}}\,\biggl{(}}-\Bigl{(}\kappa+n-2\Bigr{)}\,\phi_{n-2}+\Bigl{(}(\rho+2\,u_{0})\,(\Lambda+1)-B+\alpha/2\Bigr{)}\,\chi_{n-2}+(\Lambda+1)\,\chi_{n-3}\,\biggr{)}$
(3.13)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\chi_{n}={\displaystyle{\frac{1}{n\,(u_{0}+\rho)\,u_{0}}}}\,\biggl{(}\,-\Bigl{(}(\Lambda-1)\,u_{0}^{2}+\Bigl{(}(\Lambda-1)\,\rho-B+\alpha/2\Bigr{)}\,u_{0}+\,\rho\,\alpha/2\Bigr{)}\,\phi_{n-1}-\Bigl{(}(\rho+2\,u_{0})\,(n-1)-(\rho+u_{0})\,\kappa\Bigr{)}\,\chi_{n-1}$
(3.14)
$\displaystyle\phantom{\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\phi_{n}={\displaystyle{\frac{1}{n\,(u_{0}+\rho)\,u_{0}}}}\,\biggl{(}}-\Bigl{(}(\Lambda-1)\,(\rho+2\,u_{0})-B+\alpha/2\Bigr{)}\,\phi_{n-2}+\Bigl{(}\kappa-n+2\Bigr{)}\,\chi_{n-2}-(\Lambda-1)\,\phi_{n-3}\,\biggr{)}\,,$
(3.15)

by which accurate numerical values of the solutions are rapidly obtained.

### The solution at infinity.

General theorems [22] state that asymptotic solutions at infinity of (2.21)
are of the form

$\displaystyle\phi(s)=e^{-\omega
s}\,s^{\delta}\,\Phi(s),\qquad\chi(s)=e^{-\omega s}\,s^{\delta}\,X(s)$ (3.16)

The unknown functions $\,(\Phi(s),X(s))\,$ are solutions of the system

$\displaystyle\frac{d}{ds}\Phi(s)+\Bigl{(}-\omega+\frac{\delta+\kappa}{s}\Bigr{)}\,\Phi(s)-\Bigl{(}\,\Lambda+\frac{\alpha}{2\,s}-\frac{B}{s+\rho}+1\,\Bigr{)}\,X(s)=0\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (3.17)
$\displaystyle\frac{d}{ds}X(s)+\Bigl{(}-\omega+\frac{\delta-\kappa}{s}\Bigr{)}\,X(s)+\Bigl{(}\,\Lambda+\frac{\alpha}{2\,s}-\frac{B}{s+\rho}-1\,\Bigr{)}\,\Phi(s)=0$
(3.18)

By substituting into (3.18) the expansions

$\displaystyle\Phi(s)=\sum_{n=0}^{\infty}\Phi_{n}\,s^{-n},\qquad
X(s)=\sum_{n=0}^{\infty}X_{n}\,s^{-n}\qquad(s+\rho)^{-1}=\sum_{n=0}^{\infty}(-\rho)^{n}\,s^{-(n+1)}$
(3.19)

at the zero order we find the system

$\displaystyle\begin{pmatrix}\,-\omega&-\Lambda-1\,\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
\,\Lambda-1&-\omega\,\end{pmatrix}\,\begin{pmatrix}\,\Phi_{0}\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\ X_{0}\,\end{pmatrix}=0$ (3.20)

which yields

$\displaystyle\omega=\sqrt{1-\Lambda^{2}},\qquad\Phi_{0}=1,\qquad
X_{0}=-\omega/(1+\Lambda)$ (3.21)

Letting

$\displaystyle\Sigma_{\Phi,\,n}=\sum_{p=0}^{n}(-\rho)^{p}\,\Phi_{n-p},\qquad\Sigma_{X,\,n}=\sum_{p=0}^{n}(-\rho)^{p}\,X_{n-p}$
(3.22)

at order $\,-n\,$, with $\,n\geq 1\,$ we find the relations

$\displaystyle\begin{pmatrix}\,-\omega&-\Lambda-1\,\vrule
height=3.99994pt,width=0.0pt,depth=3.99994pt\\\
\,\Lambda-1&-\omega\,\end{pmatrix}\,\begin{pmatrix}\,\Phi_{n}\vrule
height=8.00003pt,width=0.0pt,depth=8.00003pt\\\
X_{n}\,\end{pmatrix}=\begin{pmatrix}\,-(\delta+\kappa-n+1)\,\Phi_{n-1}+(\alpha/2)\,X_{n-1}-B\,\Sigma_{X,\,n-1}\vrule
height=8.00003pt,width=0.0pt,depth=8.00003pt\\\
-(\delta-\kappa-n+1)\,X_{n-1}-(\alpha/2)\,\Phi_{n-1}+B\,\Sigma_{\Phi,\,n-1}\end{pmatrix}\,$
(3.23)

The matrix in (3.23) is singular in view of (3.20), with eigenvalues
$\,(\,0,\,-2\omega\,)\,$. Therefore we transform the system to the basis where
the matrix is diagonal. Obviously, in such a basis, one of the two relations
obtained from (3.23) is homogeneous. For $\,n=1\,$ this homogeneous equation
gives the value of $\,\delta\,$, namely

$\displaystyle\delta=-(2\,\omega)^{-1}\,\Bigl{(}\Lambda\,(2B-\alpha)\Bigr{)}$
(3.24)

In order to obtain a solvable set of recurrence equations allowing to
determine the coefficients $\,\Phi_{n},\,X_{n}\,$ for $\,n\geq 1\,$, we have
considered a first equation by subtracting the homogeneous relation from the
inhomogeneous one, and a second equation obtained by rescaling the homogeneous
equation from $\,n\,$ to $\,n+1\,$, which thus becomes an inhomogeneous
equation in $\,\Phi_{n}\,$ and $\,X_{n}\,$. The result is

$\displaystyle\Phi_{n}=\Delta^{-1}\,\biggl{(}\,-A_{X}\,C_{\Phi}~{}\Phi_{n-1}-A_{X}\,C_{X}~{}X_{n-1}+\Upsilon\,\Bigl{(}2\,\Upsilon\,\omega\,B\,\rho-
A_{X}\,B\,\Bigr{)}~{}\Sigma_{\Phi,n-1}+\Bigl{(}2\,\Upsilon\,\omega\,B\,\rho+A_{X}\,B\Bigr{)}~{}\Sigma_{X,n-1}\,\biggr{)}\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (3.25) $\displaystyle
X_{n}=\Delta^{-1}\,\biggl{(}\,\phantom{-}A_{\Phi}\,C_{\Phi}~{}\Phi_{n-1}+A_{\Phi}\,C_{X}~{}X_{n-1}+\Upsilon\,\Bigl{(}-2\,\omega\,B\,\rho+A_{\Phi}\,B\,\Bigr{)}~{}\Sigma_{\Phi,n-1}-(2\,\omega\,B\,\rho+A_{\Phi}\,B\Bigr{)}~{}\Sigma_{X,n-1}\,\biggr{)}$
(3.26)

where, for simplicity in writing, we have defined

$\displaystyle
A_{\Phi}=n-\delta-\kappa-(2\,B-\alpha)\,(\Lambda+1)/(2\,\omega),\qquad
A_{X}=(1/2)\,\alpha-(\Lambda+1)\,(n-\delta+\kappa)/\omega-B,\vrule
height=8.00003pt,width=0.0pt,depth=8.00003pt$ (3.27)
$\displaystyle\phantom{2\omega)}C_{\Phi}=n-1-\kappa-\delta-(\Lambda+1)\,\alpha/(2\,\omega),\qquad
C_{X}=(1/2)\,\alpha+(\Lambda+1)\,(n-1-\delta+\kappa)/\omega,\vrule
height=8.00003pt,width=0.0pt,depth=8.00003pt$ (3.28)
$\displaystyle\phantom{XXXXXXXXXXXX}\phantom{(2\omega),}\Upsilon=(\Lambda+1)/\omega,\qquad\phantom{{}_{X}}\Delta=(-2\,A_{\Phi}\,\Upsilon+2\,A_{X})\,\omega$
(3.29)

It is well known that the expansion of the solution at infinity produces an
asymptotic series. Thus, for convergence, in the numerical treatment we have
applied to this series the Padé approximant technique.

The results for the pure Coulomb levels of the pionic and kaonic atoms are
reported in the following Tables 1 and 2.

State | $\kappa$ |  $n=1$9011 |  $n=2$901. |  $n=3$901. |  $n=4$901. |  $n=5$901.
---|---|---|---|---|---|---
$s_{1/2}$ | $-1$ | -3235.0931859 | -808.7611897 | -359.4465684 | -202.1877820 | -129.3998072
$p_{1/2}$ | $+1$ |  -09319 | -808.7426715 | -359.4410815 | -202.1854673 | -129.3986221
$p_{3/2}$ | $-2$ |  -09319 | -808.7424909 | -359.4410280 | -202.1854447 | -129.3986105
$d_{3/2}$ | $+2$ |  -09319 |  -09319 | -359.4399305 | -202.1849817 | -129.3983734
$d_{5/2}$ | $-3$ |  -09319 |  -09319 | -359.4399127 | -202.1849742 | -129.3983696
$f_{5/2}$ | $+3$ |  -09319 |  -09319 |  -09319 | -202.1847757 | -129.3982680
$f_{7/2}$ | $-4$ |  -09319 |  -09319 |  -09319 | -202.1847720 | -129.3982661
$g_{7/2}$ | $+4$ |  -09319 |  -09319 |  -09319 |  -09319 | -129.3982096
$g_{9/2}$ | $-5$ |  -09319 |  -09319 |  -09319 |  -09319 | -129.3982085
Table 1: The pure Coulomb spectrum of the pionic atom $\,({H}^{+}\\!,\,\pi^{-})\,$ expressed in eV. We have taken the proton mass $m_{P}$=938.2720813 MeV, the pion mass $m_{\pi}$=139.57061 MeV, the fine structure constant $\alpha$=0.0072973525698163 [23]. In the first row $n=n_{r}+|\kappa|$ is the principal quantum number, where $n_{r}$ is the radial quantum number [19]. State | $\kappa$ |  $n=1$9011 |  $n=2$901. |  $n=3$901. |  $n=4$901. |  $n=5$901.
---|---|---|---|---|---|---
$s_{1/2}$ | $-1$ | -8612.9384350 | -2153.2252598 | -956.9862506 | -538.3038471 | -344.5140798
$p_{1/2}$ | $+1$ |  -09319 | -2153.2097136 | -956.9816444 | -538.3019039 | -344.5130848
$p_{3/2}$ | $-2$ |  -09319 | -2153.2063064 | -956.9806348 | -538.3014780 | -344.5128669
$d_{3/2}$ | $+2$ |  -09319 |  -09319 | -956.9797133 | -538.3010892 | -344.5126677
$d_{5/2}$ | $-3$ |  -09319 |  -09319 | -956.9793768 | -538.3009472 | -344.5125950
$f_{5/2}$ | $+3$ |  -09319 |  -09319 |  -09319 | -538.3007806 | -344.5125097
$f_{7/2}$ | $-4$ |  -09319 |  -09319 |  -09319 | -538.3007096 | -344.5124734
$g_{7/2}$ | $+4$ |  -09319 |  -09319 |  -09319 |  -09319 | -344.5124260
$g_{9/2}$ | $-5$ |  -09319 |  -09319 |  -09319 |  -09319 | -344.5124042
Table 2: The pure Coulomb spectrum of the kaonic atom
$\,({H}^{+}\\!,\,K^{-})\,$ expressed in eV. $m_{P}$ and $\alpha$ are as in
Table 1. We have taken the mass of the $K$ meson $m_{K}$=493.677 MeV [23]. In
the first row $n=n_{r}+|\kappa|$ is the principal quantum number, where
$n_{r}$ is the radial quantum number [19].

##  4 The lowest order radiative corrections

The lowest order QED corrections to the spectrum is given by the one-loop
electron vacuum polarization or Uehling potential. It takes into account the
vacuum polarization produced by the creation of electron-positron pairs, since
the contributions due to other particles, such as the lepton $\,\mu\,$, are
much smaller. Its form, in terms of the dimensionless variable $s$ and
expressed in proton mass units, is

$V_{U}(s)=-\frac{2\alpha^{2}}{3\pi
s}\,\int\limits_{1}^{\infty}{\rm{e}}^{\displaystyle{-\Bigl{(}\,\frac{2m_{e}}{m_{P}}\,sx\,\Bigr{)}}}\,\sqrt{x^{2}-1}~{}\frac{2x^{2}+1}{2x^{4}}~{}dx$
(4.1)

where $m_{P}$ is the proton mass and $m_{e}$ the electron mass. It can be
expressed in terms of special functions

$V_{U}(s)=-\frac{2\alpha^{2}}{3\pi
s}\,\left(\,\frac{m_{e}s}{4m_{P}}\,G^{30}_{13}\,\biggl{(}{{1}\atop{0,-1/2,-1/2}}\,\biggl{|}\,\Bigl{(}\frac{m_{e}}{m_{P}}\,s\Bigr{)}^{2}\biggr{)}+\Bigl{(}\frac{m_{e}s}{4m_{P}}\Bigr{)}^{3}\,G^{30}_{13}\,\biggl{(}{{1}\atop{-1/2,-1,-3/2}}\biggr{)}\,\biggl{|}\,\Bigl{(}\frac{m_{e}}{m_{P}}\,s\Bigr{)}^{2}\,\right)$
(4.2)

where
$\,G^{mn}_{pq}\bigl{(}{{a_{1},...,a_{p}}\atop{b_{1},...,b_{q}}}\,|\,z\bigr{)}$
is the Meijer G function [24].

The next order corrections are given by the two-loop irreducible part of the
the vacuum polarization, also known as Källen-Sabry potential

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!V_{KS}(s)=-\frac{2\alpha^{3}}{3\pi^{2}s}\int\limits_{0}^{1}\frac{{x\,\rm{e}}^{\displaystyle{-\Bigl{(}\,\frac{2m_{e}}{m_{P}}\frac{s}{\sqrt{1-x^{2}}}\,\Bigr{)}}}}{(1-x^{2})}\,\biggl{[}x\,(3-x^{2})\,\biggl{(}\frac{3}{2}\ln\Bigl{(}\frac{1-x^{2}}{4}\Bigr{)}-2\,\ln(x)\biggr{)}+\ln\Bigl{(}\frac{1+x}{1-x}\Bigr{)}\,\Bigl{(}\frac{11}{16}\,(3-x^{2})\,(1+x^{2})+\frac{x^{4}}{4}\Bigr{)}+\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (4.3)
$\displaystyle\phantom{xxxxx}(3-x^{2})(1+x^{2})\,\biggl{(}{\rm{Li_{2}}}\Bigl{(}-\frac{1-x}{1+x}\Bigr{)}+2{\rm{Li_{2}}}\Bigl{(}\frac{1-x}{1+x}\Bigr{)}+\ln\Bigl{(}\frac{1+x}{1-x}\Bigr{)}\,\Bigl{(}\frac{3}{2}\ln\Bigl{(}\frac{1+x}{2}\Bigr{)}-\ln(x)\Bigr{)}\biggr{)}+\frac{3}{8}x\,(5-3\,x^{2})\biggr{]}\,dx$
(4.4)

where $\,{\rm{Li_{2}}}(z)\,$ is the Dilogarithm function.

We finally add the contribution of the reducible part of the two-loop vacuum
polarization potential, namely

$V_{R}(s)=-\frac{\alpha^{3}}{9\pi^{2}s}\int\limits_{0}^{1}\frac{{x^{2}\,\rm{e}}^{\displaystyle{-\Bigl{(}\,\frac{2m_{e}}{m_{P}}\frac{s}{\sqrt{1-x^{2}}}\,\Bigr{)}}}}{(1-x^{2})}\,\Bigl{(}1-\frac{x^{2}}{3}\Bigl{)}\,\biggl{(}16-6x^{2}+3x\,(3-x^{2})\,\ln\Bigl{(}\frac{1-x}{1+x}\Bigr{)}\biggr{)}$
(4.5)

The detailed results of each one of the corrections on each atomic level are
shown in Appendix II, together with some comments. Here below we give the QED-
corrected transitions from the excited to the ground state for
$(H^{+},\pi^{-})$ and $(H^{+},K^{-})$

Transition | $n=2\phantom{110}$ | $n=3\phantom{110}$ | $n=4\phantom{110}$ | $n=5\phantom{110}$
---|---|---|---|---
${}^{n}s_{1/2}\rightarrow{}^{1}s_{1/2}$ | 2429.2202909 | 2878.7969532 | 3036.1184511 | 3108.9286213
${}^{n}p_{3/2}\rightarrow{}^{1}s_{1/2}$ | 2429.5729174 | 2878.8990487 | 3036.1611897 | 3108.9504262
${}^{n}d_{5/2}\rightarrow{}^{1}s_{1/2}$ |  -09319 | 2878.9112441 | 3036.1663831 | 3108.9530963
${}^{n}f_{7/2}\rightarrow{}^{1}s_{1/2}$ |  -09319 |  -09319 | 3036.1668317 | 3108.9533385
${}^{n}g_{9/2}\rightarrow{}^{1}s_{1/2}$ |  -09319 |  -09319 |  -09319 | 3108.9533999
Table 3: The values in eV of the transitions from the excited to the ground state of $(H^{+},\pi^{-})$. The QED corrected value of the ground level is $\epsilon_{0,\,(H^{+},\pi^{-})}=-3238.3591996$ eV.  Transition | $n=2\phantom{110}$ | $n=3\phantom{110}$ | $n=4\phantom{110}$ | $n=5\phantom{110}$
---|---|---|---|---
${}^{n}s_{1/2}\rightarrow{}^{1}s_{1/2}$ | 6479.1535386 | 7677.1203701 | 8096.2083183 | 8290.1408666
${}^{n}p_{3/2}\rightarrow{}^{1}s_{1/2}$ | 6480.8162040 | 7677.6072875 | 8096.4123161 | 8290.2449542
${}^{n}d_{5/2}\rightarrow{}^{1}s_{1/2}$ |  -09319 | 7677.7851662 | 8096.4842761 | 8290.2811709
${}^{n}f_{7/2}\rightarrow{}^{1}s_{1/2}$ |  -09319 |  -09319 | 8096.5009525 | 8290.2897706
${}^{n}g_{9/2}\rightarrow{}^{1}s_{1/2}$ |  -09319 |  -09319 |  -09319 | 8290.2910456
Table 4: The values in eV of the transitions from the excited to the ground
state of $(H^{+},K^{-})$. The QED corrected value of the ground level is
$\epsilon_{0,\,(H^{+},K^{-})}=-8634.8934971$ eV.

We have done a comparison with experimental evaluations of the transitions
between higher levels of the pionic atom and we have found complete agreement
within the experimental errors [25, 26]. The experimental results for kaonic
atom are more rare and affected by larger errors [27]. There is however
agreement between our numerical values and those reported in [27] taken from
[28].

##  5 Conclusions.

In this paper we have calculated levels and states, up to $n=5$, of the pionic
and kaonic Hydrogen atoms in a two body quantum relativistic framework for a
scalar and a fermion. The original four dimensional system is reduced to a
pair of first order differential equations formally identical to a Dirac
equation in a central potential, but for the eigenvalue. This fact, together
with the not purely Coulomb interaction, implies the non conservation of the
Johnson-Lippmann operator [29] and thus the splitting of the levels having the
same $n$ and $\kappa$ of equal absolute value but opposite sign. This is quite
evident for $n=2$, much less for increasing values of $n$ and $|\kappa|$, as
shown in Tables 1 and 2. The radiative corrections, Table 5, share the same
behavior, with the hints concerning the contributions due to $V_{R}$. The
numerical calculations have been carried on with high accuracy, in order to
give a clearcut distinction and enumeration of the levels although, nowadays,
they still cannot be experimentally observed. In any case the great physical
relevance of the mesic atoms is well known: for instance many decades ago
Dalitz underlined the great importance of the “definitive determination of the
energy level shifts in the $K$-Hydrogen and Deuterium, because of their direct
connection with the physics of $K$-nucleon interaction and their complete
independence from all other kinds of measurements which bear on this
interaction” (quoted by [9]). This justifies the efforts in order to obtain
the electromagnetic spectrum as precise as possible. The experiments now in
progress mainly in kaonic atoms [9, 30] will hopefully give a breakthrough in
this direction.

## Appendix I. Variables, kinematics, dynamics.

In previous papers we have explained in detail the kinematic setting and the
covariance of the treatment. Referring to those papers for details, we report
here, for the sake of completeness, a brief synthesis useful to contextualize
the definitions of dimensional and dimensionless variables which appear in our
equations. We are given the phase space of two relativistic points, spanned by
the Minkowski coordinates $\,x_{1}^{\mu},\,x_{2}^{\mu}$ and conjugate momenta
$\,p_{1}^{\mu},\,p_{2}^{\mu}$, the Lorentz metric tensor being
$\eta^{\mu}_{\nu}={{\rm diag}}(1,-1,-1,-1).$ It is first instrumental to
introduce a canonical transformation to ‘relative’ and ‘global’ coordinates
$\,\tilde{r}^{\mu}=x_{1}^{\mu}-x_{2}^{\mu},\,\,X^{\mu}=(1/2)\,(x_{1}^{\mu}+x_{2}^{\mu}),$
with conjugate momenta
$\,\tilde{q}^{\mu}=(1/2)\,(p_{1}^{\mu}-p_{2}^{\mu}),\,\,P^{\mu}=p_{1}^{\mu}+p_{2}^{\mu}$.
We then make a second canonical transformation defined by its action on the
momenta: it acts as the identity on $P^{\mu}$ and transforms $\tilde{q}^{\mu}$
by boosting it to the $\bm{P}=0$ frame, thus defining
$\,q_{\alpha}=\varepsilon_{\alpha}^{\mu}(P)\,\tilde{q}_{\alpha}$, where
$\,\alpha=(0,A),$ $A=1,2,3\,$ and
$\,\varepsilon_{\alpha}^{\mu}(P)=L^{-1}(\bm{P})_{\alpha}^{\mu}$ is the matrix
of the Lorentz boost. This is a point-like transformation easily completed by
the generating function yielding the corresponding coordinates conjugate to
$q_{\alpha}$ and $P^{\mu}$, namely:
$\,r_{\alpha}=\varepsilon_{\alpha}^{\mu}(P)\,\tilde{r}_{\alpha}\,$ and
$\,Z^{\mu}=X^{\mu}+P^{-2}\,\bigl{(}P^{\nu}\tilde{q}_{\nu}\,\tilde{r}^{\mu}-P^{\nu}\tilde{r}_{\nu}\,\tilde{q}^{\mu}+(P_{0}+P)^{-1}\,W^{0\mu}\bigr{)},\,$
where $\,W^{\mu\nu}\,$ is the Pauli-Lubanski tensor. It is then obvious that
$\,q_{0}\,$ and $\,r_{0}\,$ are Lorentz invariant such as
$\,q=\bigl{(}q_{A}\,q_{A}\bigr{)}^{1/2},~{}r=\bigl{(}r_{A}\,r_{A}\bigr{)}^{1/2},\,$
since $\,\bm{q}=\\{q_{A}\\}\,$ and $\,\bm{r}=\\{r_{A}\\}\,$ are Wigner
vectors. Moreover Z is a position vector of the Newton-Wigner type and
$\,Z^{0}\,$ has the covariance of a Lorentz time. The essential point is that
the mass shell conditions $\,p_{1}^{2}=m_{1}^{2}\,$ and
$\,p_{2}^{2}=m_{2}^{2}\,$ expressed in these coordinates, amount to the
elimination of $\,q_{0}\,$ from the resulting energy (and Hamiltonian) of the
relative system:
$\,(m_{1}^{2}+\bm{q}^{2})^{1/2}+(m_{2}^{2}+\bm{q}^{2})^{1/2},\,$ which plays
the role of the invariant mass. We then see that $\,Z^{0}\,$ remains the only
one time in the system, since the relative time $\,r_{0}\,$ is a cyclic
variable and disappears: we have therefore realized a canonical reduction of
the phase space where $\,m_{1}\,$ and $\,m_{2}\,$ appear as fixed parameters,
overcoming the problem of the relative energy which becomes an irrelevant
quantity. The cyclic time $\,r_{0}\,$ assumes the role of a gauge-type
variable, arbitrarily fixed in order to recover separate descriptions of the
two particles world lines. The physical mass spectrum is determined by the
invariant Hamiltonian completed with invariant interaction terms depending on
$\,\bm{q}\,$ and $\,\bm{r}.$ Thus the action of the Poincaré group is well
defined and the dynamics is covariant. Finally, the quantization is done in
the tensor product of the particle Hilbert spaces, using gamma matrices for
fermions and an adapted Feshbach-Villars representation for the scalars, so
that the square roots are eliminated from $\,H\,$ in any case.

## Appendix II. The lowest order radiative correction.

We report here the detailed list of the QED corrections to the levels of
Tables 1 and 2.

  State  $\kappa$  $V_{U}~{}{(H^{+},\,\pi^{-})}$
$V_{KS}~{}{(H^{+},\,\pi^{-})}$ $V_{R}~{}{(H^{+},\,\pi^{-})}$
$V_{U}~{}{(H^{+},\,K^{-})}$ $V_{KS}~{}{(H^{+},\,K^{-})}$
$V_{R}~{}{(H^{+},\,K^{-})}$90 $1s_{1/2}$  $-1$
-3.2413174$\cdot{10^{-1}}$  -2.0900666$\cdot{10^{-2}}$
3.7955745$\cdot{10^{-3}}$ -21.7953398$\cdot{10^{-3}}$
-1.1473759$\cdot{10^{-1}}$  4.4984666$\cdot{10^{-2}}$ $2s_{1/2}$  $-1$
 -3.6831740$\cdot{10^{-1}}$ -2.2911091$\cdot{10^{-3}}$
4.8071080$\cdot{10^{-4}}$ -2.4169328$\cdot{10^{-3}}$
-1.2919109$\cdot{10^{-2}}$  5.1225963$\cdot{10^{-3}}$ $2p_{1/2}$  $+1$
 -3.5796665$\cdot{10^{-2}}$ -3.7779853$\cdot{10^{-4}}$
-2.5718595$\cdot{10^{-5}}$ -7.7448229$\cdot{10^{-1}}$
-6.5377131$\cdot{10^{-3}}$ -8.4622213$\cdot{10^{-6}}$ $2p_{3/2}$  $-2$
 -3.5796506$\cdot{10^{-2}}$ -3.7779729$\cdot{10^{-4}}$
-2.5718666$\cdot{10^{-5}}$ -7.7447115$\cdot{10^{-1}}$
-6.5376417$\cdot{10^{-3}}$ -8.4748167$\cdot{10^{-6}}$ $3s_{1/2}$  $-1$
 -1.0756349$\cdot{10^{-1}}$ -6.6578532$\cdot{10^{-4}}$
1.4242775$\cdot{10^{-4}}$ -6.9476605$\cdot{10^{-1}}$
-3.6514009$\cdot{10^{-3}}$  1.5104177$\cdot{10^{-3}}$ $3p_{1/2}$  $+1$
 -1.1407162$\cdot{10^{-2}}$ -1.1722417$\cdot{10^{-4}}$
-7.3772433$\cdot{10^{-6}}$ -2.1392561$\cdot{10^{-1}}$
-1.7166683$\cdot{10^{-3}}$  3.3703260$\cdot{10^{-5}}$ $3p_{3/2}$  $-2$
 -1.1407110$\cdot{10^{-2}}$ -1.1722378$\cdot{10^{-4}}$
-7.3772720$\cdot{10^{-6}}$ -2.1392247$\cdot{10^{-1}}$
-1.7166496$\cdot{10^{-3}}$  3.3698833$\cdot{10^{-5}}$ $3d_{3/2}$  $+2$
 -4.4463586$\cdot{10^{-4}}$ -6.3140024$\cdot{10^{-6}}$
-7.0773072$\cdot{10^{-7}}$ -3.8501657$\cdot{10^{-2}}$
-4.4371663$\cdot{10^{-4}}$ -3.9638974$\cdot{10^{-5}}$ $3d_{5/2}$  $-3$
 -4.4463503$\cdot{10^{-4}}$ -6.3139928$\cdot{10^{-6}}$
-7.0772989$\cdot{10^{-7}}$ -3.8501463$\cdot{10^{-2}}$
-4.4371479$\cdot{10^{-4}}$ -3.9638894$\cdot{10^{-5}}$ $4s_{1/2}$  $-1$
 -4.5156315$\cdot{10^{-2}}$ -2.7904635$\cdot{10^{-4}}$
6.0079123$\cdot{10^{-5}}$ -2.9047738$\cdot{10^{-1}}$
-1.5204611$\cdot{10^{-3}}$ 6.3556694$\cdot{10^{-4}}$ $4p_{1/2}$  $+1$
 -4.9207863$\cdot{10^{-3}}$ -5.0167795$\cdot{10^{-5}}$
-3.0743111$\cdot{10^{-6}}$ -8.9047505$\cdot{10^{-2}}$
-7.0538006$\cdot{10^{-4}}$  1.7877271$\cdot{10^{-5}}$ $4p_{3/2}$  $-2$
 -4.9207642$\cdot{10^{-3}}$ -5.0167631$\cdot{10^{-5}}$
-3.0743244$\cdot{10^{-6}}$ -8.9046201$\cdot{10^{-2}}$
-7.0537249$\cdot{10^{-4}}$  1.7875322$\cdot{10^{-5}}$ $4d_{3/2}$  $+2$
 -2.4729903$\cdot{10^{-4}}$ -3.4598119$\cdot{10^{-6}}$
-3.8498620$\cdot{10^{-7}}$ -1.8088916$\cdot{10^{-2}}$
-1.9921591$\cdot{10^{-4}}$ -1.6422068$\cdot{10^{-5}}$ $4d_{5/2}$  $-3$
 -2.4729857$\cdot{10^{-4}}$ -3.4598067$\cdot{10^{-6}}$
-3.8498576$\cdot{10^{-7}}$ -1.8088823$\cdot{10^{-2}}$
-1.9921508$\cdot{10^{-4}}$ -1.6422045$\cdot{10^{-5}}$ $4f_{5/2}$  $+3$
 -4.5984890$\cdot{10^{-6}}$ -6.4383418$\cdot{10^{-8}}$
-9.7052084$\cdot{10^{-9}}$ -1.8361008$\cdot{10^{-3}}$
-2.6515938$\cdot{10^{-5}}$ -3.0478845$\cdot{10^{-6}}$ $4f_{7/2}$  $-4$
 -4.5984835$\cdot{10^{-6}}$ -6.4383351$\cdot{10^{-8}}$
-9.7051988$\cdot{10^{-6}}$ -1.8360956$\cdot{10^{-3}}$
-2.6515873$\cdot{10^{-5}}$ -3.0478779$\cdot{10^{-6}}$ $5s_{1/2}$  $-1$
 -2.3068171$\cdot{10^{-2}}$ -2.7904635$\cdot{10^{-4}}$
3.0758120$\cdot{10^{-5}}$ -1.4813228$\cdot{10^{-1}}$
-7.7410848$\cdot{10^{-4}}$  3.2499175$\cdot{10^{-4}}$ $5p_{1/2}$  $+1$
 -2.5443058$\cdot{10^{-3}}$ -5.0167795$\cdot{10^{-5}}$
-1.5648296$\cdot{10^{-6}}$ -4.5359931$\cdot{10^{-2}}$
-3.5741716$\cdot{10^{-4}}$  9.9266091$\cdot{10^{-6}}$ $5p_{3/2}$  $-2$
 -2.5442943$\cdot{10^{-3}}$ -5.0167631$\cdot{10^{-5}}$
-1.5648367$\cdot{10^{-5}}$ -4.5359271$\cdot{10^{-2}}$
-3.5741336$\cdot{10^{-4}}$  9.9255934$\cdot{10^{-6}}$ $5d_{3/2}$  $+2$
 -2.4729903$\cdot{10^{-4}}$ -3.4598119$\cdot{10^{-6}}$
-2.1636864$\cdot{10^{-7}}$ -9.6488853$\cdot{10^{-3}}$
-1.0462533$\cdot{10^{-4}}$ -8.3402250$\cdot{10^{-6}}$ $5d_{5/2}$  $-3$
 -2.4729857$\cdot{10^{-4}}$ -3.4598067$\cdot{10^{-6}}$
-2.1636840$\cdot{10^{-7}}$ -9.6488355$\cdot{10^{-3}}$
-1.0462489$\cdot{10^{-4}}$ -8.3402159$\cdot{10^{-6}}$ $5f_{5/2}$  $+3$
 -4.5984890$\cdot{10^{-6}}$ -8.0260480$\cdot{10^{-8}}$
-7.7854603$\cdot{10^{-9}}$ -1.2641487$\cdot{10^{-3}}$
-1.7661783$\cdot{10^{-5}}$ -1.9992710$\cdot{10^{-6}}$ $5f_{7/2}$  $-4$
 -4.5984835$\cdot{10^{-6}}$ -8.0260397$\cdot{10^{-8}}$
-7.7854526$\cdot{10^{-9}}$ -1.2641450$\cdot{10^{-3}}$
-1.7661739$\cdot{10^{-5}}$ -1.9992667$\cdot{10^{-6}}$ $5g_{7/2}$  $+4$
 -3.7750656$\cdot{10^{-8}}$ -0.7700092$\cdot{10^{-9}}$
-0.0927066$\cdot{10^{-9}}$ -7.6482739$\cdot{10^{-5}}$
-1.3143021$\cdot{10^{-6}}$ -1.5986779$\cdot{10^{-7}}$ $5g_{9/2}$  $-5$
 -3.7750621$\cdot{10^{-8}}$ -0.7700085$\cdot{10^{-9}}$
-0.0927065$\cdot{10^{-9}}$ -7.6482586$\cdot{10^{-5}}$
-1.3142998$\cdot{10^{-6}}$ -1.5986751$\cdot{10^{-7}}$

Table 5: QED corrections to the $\,({H}^{+}\\!,\,\pi^{-})\,$ and
$\,({H}^{+}\\!,\,K^{-})\,$ levels. In the first two columns we indicate the
state and the corresponding $\kappa$. In the last six columns of the first
line we specify the atomic systems and the potentials responsible for the
energy corrections. The results are expressed in eV.

## Appendix III. How the reduced mass can be brought to bear in a two-body
relativistic framework.

The notion of ‘reduced mass’ comes from the reduction of the Newtonian two-
body dynamics and is indeed artificial in a relativistic formulation of a two-
body problem. Despite this, however, this quantity may produce some benefit
even in a relativistic treatment, giving an approximation of the levels to the
order $\alpha$. The accuracy could possibly be increased by carrying on series
expansions to higher orders. Indeed we are able to single out the interaction
term of (2.21) leading to the reduced mass and we give some arguments to show
the connection.

We recall that the levels of the Dirac equation in a Coulomb field with
coupling constant $\alpha$, in units of the fermion mass have the well known
analytic expression

$\displaystyle\epsilon_{D}(\alpha,\kappa,n_{r})=\biggl{(}\,1+\frac{\alpha^{2}}{\Bigl{(}\,n_{r}+\sqrt{\kappa^{2}-\alpha^{2}}\,\Bigr{)}^{2}}\,\biggr{)}^{-1/2}-1$
(III.1)

where $\,n_{r}\,$ is the radial quantum number [19]. For future convenience we
observe that, for constant $c$,

$\displaystyle\epsilon_{D}(c\,\alpha,\kappa,n_{r})=c^{2}\,\epsilon_{D}(\alpha,\kappa,n_{r})+O(\alpha^{4})$
(III.2)

Taking into account the definitions of the parameters, if we develop the two
equations (2.21) to the first order in $\alpha$ and the subsequent result to
the first order in $\epsilon$ we find a new Dirac equation in the Coulomb
field with interaction constant $\mu_{R}\,\alpha$ and eigenvalue
$\mu_{R}\,\epsilon$, namely

$\displaystyle\frac{d}{ds}\phi(s)+\frac{\kappa}{s}\,\phi(s)-\Bigl{(}\,\mu_{R}\,\frac{\alpha}{s}+\mu_{R}\,\epsilon+2\,\Bigr{)}\,\chi(s)=0\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (III.3)
$\displaystyle\frac{d}{ds}\chi(s)-\frac{\kappa}{s}\,\,\chi(s)+\Bigl{(}\,\mu_{R}\,\frac{\alpha}{s}+\mu_{R}\,\epsilon\,\Bigr{)}\,\phi(s)=0$
(III.4)

where $\mu_{R}=\sigma/(1+\sigma)$. Therefore, using the relation (III.2) with
$c=\mu_{R}$, the solution to the first order in $\alpha$ of the original
system (2.21) is given by

$\displaystyle\epsilon=\mu_{R}\,\epsilon_{D}(\alpha,\kappa,n_{r})$ (III.5)

Multiplying both sides by the fermion mass $m_{F}$, we finally see that

$E=m_{F}\,\epsilon=m_{R}\,\epsilon_{D}(\alpha,\kappa,n_{r})$ (III.6)

where $m_{R}=m_{S}\,m_{F}/(m_{S}+m_{F})$ is the classical reduced mass. Notice
that for $\sigma=1$, namely when the scalar and fermion masses are equal, the
coefficient $B$ is vanishing and (2.21) is exactly the Dirac equation (III.4)
with $\mu_{R}=1/2$. In this case the Coulomb coupling constant is $\alpha/2$
and $\Lambda=\lambda/2$. The perturbation treatment of (2.21) could start from
the eigenvalues (III.1) with $\alpha$ substituted by $\mu_{R}\,\alpha$. The
corresponding eigenfunctions of (III.4)

$\displaystyle\phi=\phantom{-}e^{\displaystyle{-as}}\,(as)^{\displaystyle{-1+\gamma}}\,(\,1+\mu_{R}\epsilon\,)^{1/2}\,\Bigl{(}\,-N_{1}\,L_{n_{r}-1}^{2\gamma}(2as)+N_{2}\,L_{n_{r}}^{2\gamma}(2as)\,\Bigr{)}\vrule
height=10.00002pt,width=0.0pt,depth=10.00002pt$ (III.7)
$\displaystyle\chi=-e^{\displaystyle{-as}}\,(as)^{\displaystyle{-1+\gamma}}\,(\,1-\mu_{R}\epsilon\,)^{1/2}\,\Bigl{(}\,\phantom{-}N_{1}\,L_{n_{r}-1}^{2\gamma}(2as)+N_{2}\,L_{n_{r}}^{2\gamma}(2as)\,\Bigr{)}$
(III.8)

where $L_{n}^{b}$ are the associated Legendre polynomials [31],
$\gamma=(\kappa^{2}+\mu_{R}^{2}\alpha^{2})^{1/2}$,
$a=\mu_{R}\alpha\,(\kappa^{2}+n_{r}^{2}+2n_{r}\gamma)^{-1/2}$ and the
normalization constants $N_{1},\,N_{2}$ can be recovered from [19]. We finally
observe that, being a genuine Dirac equation, (III.4) conserves the Johnson-
Lippmann operator which is not conserved in (2.21). This suggests that indeed
the breaking has to be very small, as confirmed by the numerical results.

## References

  * [1] R. Giachetti, E. Sorace, J. Phys. A, 38, 1345 (2005).
  * [2] R. Giachetti, E. Sorace, J. Phys. A 39, 15207 (2006).
  * [3] R. Giachetti, E. Sorace, Phys. Rev. D, 87, 034021 (2013).
  * [4] A. Barducci, R. Giachetti, E. Sorace, J. Phys. B: At. Mol. Opt. Phys. 48, 085002 (2015).
  * [5] A. Barducci, R. Giachetti, E. Sorace, Phys. Rev. D, 95, 054022 (2017).
  * [6] R. Giachetti, E. Sorace, Annals of Physics 401, 202 (2019).
  * [7] M. Iwasaki et al., Nucl. Phys. A 639, 501 (1998).
  * [8] H.Ch. Schröder et al., Phys. Lett. B 469, 25 (1999).
  * [9] C. Curceanu, C. Guaraldo et al., Symmetry, 12, 547 (2020).
  * [10] M. Bazzi et al., Phys. Lett. B 704, 113 (2011).
  * [11] E. Matsinos, A brief history of the pion-nucleon coupling constant, e-Print: arXiv:1901.01204 (2019).
  * [12] M. Trassinelli, P. Indelicato, Phys. Rev. A 76 , 012510 (2007).
  * [13] M. Trassinelli, Phys. Lett. B 759, 583 (2016).
  * [14] S. Schlesser, E.O. Le Bigot, P. Indelicato, K. Pachucki, Phys. Rev. C 84 , 015211 (2011).
  * [15] G. Ron, The proton radius, Talk 11/01/2018, Open Access at Springerlink.com,
https://www.npl.washington.edu/sites/default/files/speaker-materials/Proton

  * [16] A.S.M. Patoary, N.S. Oreshkina, Eur. Phys. J. D72 , 54 (2018).
  * [17] W. Xiong, A. Gasparian, H.Gao et al., Nature 575 , 147 (2019).
  * [18] M. Daum, R. Frosch, P.R. Kettle, Phys. Lett. B 796, 11 (2019).
  * [19] V.B. Berestetskii, E.M. Lifshitz, L.P. Pitaevskii, Quantum Electrodynamics (Pergamon Press, 1982).
  * [20] A.P. Martynenko, R.N. Faustov, Theor. Math. Phys. 64, 765 (1985).
  * [21] H. Feshbach, F. Villars, Rev. Mod. Phys. 30 (1958) 1.
  * [22] F.T. Cope, Am. J. Math. 56, 411 (1934) and 58 130 (1936).
  * [23] Review of Particle Physics, Particle Data Group, Chinese Physics C, 38, Number 9, (2014).
  * [24] R. Beals, J. Szmigielski, Notices of the AMS 60, 866 (2013),
  * [25] M. Hennebach et al., Eur. Phys. J. A50, 329 (2014), Erratum: Eur .Phys. J. A55, 24 (2019).
  * [26] D. Gotta et al., Physics Procedia 17, 69 (2011).
  * [27] C. Corceanu et al., Rev. Mod Phys. 91, 1 (2019).
  * [28] M. Iwasaki et al., Phys. Rev. Lett. 78, 3067 (1997).
  * [29] M.H. Johnson, B.A. Lippmann, Phys. Rev 78 , 329 (1950).
  * [30] J. Marton et al., EPJ Web of Conference 199, 03004 (2019).
  * [31] F.B. Pidduck, London J. of London Math. Soc. s1-4, 163 (1929).


# Viscous effects on plasmoid formation from nonlinear resistive tearing
growth in a Harris sheet

Nisar AHMAD1, Ping ZHU(朱平)2,3∗, Ahmad ALI4, and Shiyong ZENG(曾市勇)5 1 CAS Key
Laboratory of Geospace Environment and Department of Engineering and Applied
Physics, University of Science and Technology of China, Hefei 230026, People’s
Republic of China 2 International Joint Research Laboratory of Magnetic
Confinement Fusion and Plasma Physics, State Key Laboratory of Advanced
Electromagnetic Engineering and Technology, School of Electrical and
Electronic Engineering, Huazhong University of Science and Technology, Wuhan,
430074, People’s Republic of China  3 Department of Engineering Physics,
University of Wisconsin-Madison, Madison, Wisconsin 53706, USA  4 Pakistan
Tokamak Plasma Research Institute, Islamabad 3329, Pakistan  5 Department of
Plasma Physics and Fusion Engineering, University of Science and Technology of
China, Hefei 230026, People’s Republic of China<EMAIL_ADDRESS>
###### Abstract

In this study, the evolution of a highly unstable ${m=1}$ resistive tearing
mode, leading to plasmoid formation in a Harris sheet is studied in the
framework of full MHD model using the NIMROD simulation. Following the initial
nonlinear growth of the primary $m=1$ island, the X-point develops into a
secondary elongated current sheet that eventually breaks into plasmoids. Two
distinctive viscous regimes are found for the plasmoid formation and
saturation. In the low viscosity regime (i.e. $P_{r}\lesssim 1$), the plasmoid
width increases sharply with viscosity, whereas in the viscosity dominant
regime (i.e. ${P_{r}\gtrsim 1}$ ), the plasmoid size gradually decreases with
viscosity. Such a finding quantifies the role of viscosity in modulating the
plasmoid formation process through its effects on the plasma flow and the
reconnection itself.

Keywords: viscosity, reconnection, plasmoids, Prandtl number

(Some Figures may appear in colour only in the online journal)

##  1 Introduction

Plasmoid instability (PI) is known to develop on the elongated current sheet
formed during the externally driven Sweet-Parker (SP) reconnection, or from
the intrinsically growing nonlinear kink or tearing mode. In general, when the
aspect ratio of the elongated current sheet becomes sufficiently large,
unstable secondary tearing can lead to the formation of plasmoids [1, 2]. The
problem of the transition from the laminar reconnection during the early
nonlinear stage, to the subsequent highly unstable one, characterized by
sporadic production of plasmoids inside the sheet itself, with faster average
reconnection rates, has been addressed by a number of past numerical and
theoretical studies [3, 4], in the context of PIs following the externally
driven SP reconnection [5, 6, 7, 8, 9], or the intrinsically nonlinear tearing
mode [2, 10, 11, 12], on the scaling and dynamics of plasmoid formation with
different Lundquist numbers.

Previous studies have found the critical roles of plasma flow in the processes
of reconnection in general and plasmoid formation in particular [13, 14, 15,
16, 17, 18]. Whereas the plasma outflow is stabilizing on the primary tearing
mode or reconnection process [13, 14, 15, 16], the effects of plasma flow
itself, including both inflow and outflow, may contribute to the initial onset
of plasmoid instability [5]. The plasma viscosity can affect the properties
and topologies of plasma flow close to the thin current sheet as well as the
reconnection rate. Because of the narrowness of the current sheet, viscosity
can influence the non-linear regime. In fact, viscosity increases the
possibility of local changes in the flow topology. The robustness of the flow
cells around the sheet might be weakened or even unstable due to the existence
of strong flow gradients in the current sheet region [19]. Finite viscosity
inserts dissipation to the flow patterns that in turn interact with the island
evolution and reconnection [19]. Thus the plasma viscosity, both collisional
and collisionless, is expected to be one of the key parameters that determine
the onset and saturation conditions for the plasmoid instability.

The effects of viscosity on linear and nonlinear resistive tearing mode as
well as plasmoid instability have been studied by many [9, 10, 20, 21, 22, 19,
23, 24, 25, 26, 27, 28]. In this paper we focus on exploring the impact of
viscosity on the onset and saturation of plasmoid instability. Most of the
past studies on visco-resistive tearing and kink modes were made in a 2D
reduced MHD model, however in this study, we use the complete resistive MHD
equations implemented in NIMROD code [29]. Both the onset and the dynamics of
plasmoid differ greatly from those found in previous reduced MHD simulations
[18, 19].

The rest of the paper is organized as follows. In section 2, we briefly
describe our simulation model. In section 3, both linear and nonlinear
simulation results are reported. At the end in section 4, summary and
discussion are presented.

##  2 Simulation model and equilibrium

Our simulations are based on the single-fluid full MHD model implemented in
the NIMROD (Non-Ideal Magnetohydrodynamics with Rotation, Open Discussion)
code [29].

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0$ (1)

$\rho\left(\frac{\partial}{\partial
t}+\mathbf{v}\cdot\nabla\right)\mathbf{v}=\mathbf{J}\times\mathbf{B}-\nabla
p+\rho\nu\nabla^{2}\mathbf{v}$ (2)

$\frac{N}{\gamma-1}\left(\frac{\partial}{\partial
t}+\mathbf{v}\cdot\nabla\right)\mathbf{T}=-{\frac{p}{2}}\nabla\cdot\mathbf{v}-\nabla\cdot\mathbf{q}$
(3)

$\frac{\partial\mathbf{B}}{\partial
t}=\nabla\times(\mathbf{v}\times\mathbf{B}-\eta\mathbf{J})$ (4)

$\nabla\times\mathbf{B}=-\mu_{0}\mathbf{J}$ (5)

${\bf{q}}=-N\left(\chi_{\parallel}\nabla_{\parallel}{T}+\chi_{\perp}\nabla_{\perp}{T}\right)$
(6)

where
${\rho,N,p,\mathbf{J},\mathbf{v},\mathbf{B},{\bf{q}},\eta,\gamma,\nu,\chi_{\parallel}}$
and ${\chi_{\perp}}$ are the plasma mass density, number density, pressure,
current density, velocity, magnetic field, heat flux, resistivity, specific
heat ratio, viscosity, parallel and perpendicular thermal conductivity,
respectively. The Lundquist number ${S=\tau_{R}/\tau_{A}}$, where
${\tau_{R}=\mu_{0}a^{2}/\eta}$ is the resistive time and ${\tau_{A}=a/v_{A}}$
is the Alfvénic time (${a}$ represents the current sheet width), the Alfvén
speed ${v_{A}=B_{0}/\sqrt{\mu_{0}\rho}}$ (${B_{0}}$ is the magnitude of
magnetic field at the edge of the current sheet), and ${P_{r}=\nu/\eta}$ is
the Prandtl number. The Harris current sheet model is adopted for the
equilibrium magnetic field [30]

$\mathbf{B_{0}}{(x)}=B_{0}\mathrm{tanh}\left(\frac{x}{a}\right)\hat{\mathbf{y}}$
(7)

The corresponding pressure profile from the static MHD force balance is
determined as

${P_{0}}{(x)}=\frac{B_{0}^{2}}{2\mu_{0}}\mathrm{sech^{2}}\left(\frac{x}{a}\right)$
(8)

The resistive MHD equations (1)–(6) are numerically solved in a rectangular
domain ${[-L_{x},L_{x}]\times{[-L_{y},L_{y}]}}$. The periodic boundary
conditions are imposed at the ${y}$-boundaries, and the solid, perfect
conducting walls are assumed at the ${x}$-boundaries. For the Harris current
sheet ${\Delta^{\prime}a=2[(ka)^{-1}-ka]}$ [31], so that the unstable modes
have wave vector satisfying ${ka<1}$. Here ${\Delta^{\prime}}$ is the
discontinuity of logarithmic derivative of the outer flux function when
approaching the singular layer at ${x=0}$, which is a measure of the free
energy of the system. In our simulations ${\Delta^{\prime}=49.66}$, ${k}$ is
the mode wave number along ${y}$, and ${L_{y}=2\pi m/k}$, with ${m}$ being the
mode number. Simulations are performed for a uniform plasma resistivity
${\eta=2.8\times 10^{-4}}$ and a wide range of viscosity ( ${P_{r}}$ = 0.33 to
10), and the equilibrium plasma number density ${n=3\times 10^{19}m^{-3}}$ .

(a)

(b)

Figure 1: (a) Harris sheet equilibrium magnetic field and (b) pressure
profiles.

(a)

(b)

Figure 2: (a) Linear growth rates as functions of the resistivity for the
fixed value of ${\Delta^{\prime}=49.66}$. The blue and red lines, represent
the theoretical scaling of the linear growth rate for low
(${\gamma\sim\eta^{0.33}}$) and large (${\gamma\sim\eta^{0.66}}$) Prandtl
number regimes respectively, whereas, the green line represents the simulation
results. (b) Linear growth rates as functions of the viscosity for the fixed
value of instability parameter ${\Delta^{\prime}=49.66}$ and various values of
resistivity ${\eta}$.

##  3 Simulation results

###  3.1 Linear scaling

The plasmoid instability tends to develop from the primary tearing growth in
the large ${\Delta^{\prime}}$ regime [2]. One such case,
${\Delta^{\prime}=49.66}$ is examined first in simulations for its linear
scaling in comparison with theory. The linear growth rate of the ${m=1}$
resistive tearing mode obtained from our NIMROD simulations scales with the
resistivity ${\gamma}$ as ${\gamma\sim\eta^{0.30}}$, which is close to the
resistive tearing scaling of ${\gamma\sim\eta^{0.33}}$ in the large
${\Delta^{\prime}}$ regime previously derived in theory [33, 21] (Figure 2a).
Viscosity in general introduces dissipation that reduce the linear growth of
resistive tearing mode. This viscous dissipation is stronger in the
${P_{r}>1}$ regime, where the growth rate ${\gamma}$ of the ${m=1}$ resistive
tearing scales with the viscosity as ${\gamma\sim\nu^{-0.33}}$ , similar to
the viscosity scaling obtained in previous reduced MHD simulations [26] and
theory [21] (Figure 2b).

###  3.2 Nonlinear results

####  3.2.1 Critical $\Delta^{\prime}$ for the X-point collapse and plasmoid
instability

Our nonlinear simulations find that the onset of secondary tearing instability
and plasmoid formation occur only when the ${\Delta^{\prime}}$ is above a
certain threshold value. In our nonlinear simulation, we mostly employ
${64\times 48}$ 2D ﬁnite elements with a polynomial degree of 5, which ensures
the numerical convergence (Figure 3).

Figure 3: Kinetic energy evolution for different numerical resolutions at
${\Delta^{\prime}=49.66}$. Figure 4: Kinetic energy evolution for
$\Delta^{\prime}$ = 9.306, 12.057, 13.462, 14.05, 14.53, 15.6419, 17.31,
24.5007, 39.35 and 49.66.

(a)

(b)

Figure 5: 2D contours of the current density in z-direction and the 2D
magnetic field lines at (a) time = 450 sec and (b) time = 570 sec for
${\Delta^{\prime}=13.462}$.

The evolution of kinetic energy reaches its maximum sooner as we increase the
value of ${\Delta^{\prime}}$ (Figure 4). The minimum value of
${\Delta^{\prime}}$ at which the X-point evolves into the Y-type is ${14.04}$.
For ${\Delta^{\prime}=13.42}$, the current sheet remains the shape of X-point
over the entire time (Figure 5). As the value of ${\Delta^{\prime}}$ increases
to above 14.05, the X-point evolves into a Y-type current sheet as shown in
Figures 6 and 7. Waelbroeck [34] first predicted the criterion for the
collapse of X-point into Y-type current sheet to be ${W>W_{c}\approx
25/\Delta^{\prime}}$, where ${W}$ represents the width of primary island and
${W_{c}}$ represents the critical width at which X-point collapses. This
conversion of X-point into Y-type current sheet is termed as the secondary
instability (SI)[34]. The subsequent collapse of the Y-type current sheet into
plasmoids is known as the plasmoid instability (PI) [35].

For ${\eta=2.8\times 10^{-4}}$ and ${P_{r}=1}$ the critical
${\Delta^{\prime}}$ value for the onset of plasmoid instability (PI) is
${\Delta^{\prime}=17.31}$, which is in agreement with the previous reduced MHD
simulation [18]. In this case, at time 200 Sec, Figure 8(a), an X-point is
formed as shown. At time 340 Sec, Figure 8(b), the Y-type current sheet
develops during the nonlinear stage [3, 6, 36]. At time 360 Sec, Figure 8(c)
the current density reaches maximum and the current sheet becomes more
stretched and thinner. Finally, at time 365 Sec, Figure 8(d), formation of
plasmoid takes place.

(a)

(b)

Figure 6: 2D contours of the current density in z-direction and the 2D
magnetic field lines at (a) time = 300 sec and (b) time = 470 sec for
${\Delta^{\prime}=14.05}$. Figure 7: 2D contours of the current density in
z-direction and the 2D magnetic field lines at (a) time = 300 sec, (b) time =
400 sec, (c) time = 430 sec and (d) time = 500 sec for
${\Delta^{\prime}=15.6419}$. Figure 8: 2D contours of the current density in
z-direction and the 2D magnetic field lines at (a) time = 200 sec, (b) time =
340 sec, (c) time = 360 sec and (d) time = 365 sec for
${\Delta^{\prime}=17.31}$. Figure 9: Kinetic energy evolution for
$\Delta^{\prime}$ = 49.66 at $P_{r}$ = 0.33, 0.5, 1 and $P_{r}$ = 10. Figure
10: Growth rate evolution for all ${P_{r}=0.33,0.5,1}$ and ${10}$ cases with
$\Delta^{\prime}$=49.66. Figure 11: 2D contours of the current density in
z-direction and the 2D magnetic field lines at (a) time = 70 sec, (b) time =
100 sec, (c) time = 152 sec, (d) time = 168 sec, (e) time = 172 sec and (f)
time = 176 sec for $\Delta^{\prime}$=49.66, ${P_{r}=0.33}$. Figure 12: 2D
contours of the current density in z-direction and the 2D magnetic field lines
at (a) time = 100 sec, (b) time = 156 sec, (c) time = 180 sec, (d) time = 196
sec, (e) time = 204 sec and (f) time = 208 sec for $\Delta^{\prime}$=49.66,
${P_{r}=0.5}$. Figure 13: 2D contours of the current density in z-direction
and the 2D magnetic field lines at (a) time = 100 sec, (b) time = 150 sec, (c)
time = 190 sec, (d) time = 200 sec, (e) time = 205 sec and (f) time = 215 sec
for $\Delta^{\prime}$=49.66, ${P_{r}=1}$. Figure 14: 2D contours of the
current density in z-direction and the 2D magnetic field lines at (a) time =
150 sec, (b) time = 200 sec, (c) time = 280 sec, (d) time = 292 sec, (e) time
= 296 sec and (f) time = 304 sec for $\Delta^{\prime}$=49.66, ${P_{r}=10}$.

(a)

(b)

Figure 15: (a) Saturated island width and (b) plasmoid width as functions of
the Prandtl number.

####  3.2.2 Effect of viscosity on the non-linear evolution of resistive
tearing mode for highly unstable system

In our cases, we only vary the viscosity in terms of ${P_{r}}$ by keeping the
resistivity and instability parameter constant, which are ${\eta=2.8\times
10^{-4}}$ and ${\Delta^{\prime}=49.66}$ respectively. To study the effect of
viscosity, we choose four different cases with ${P_{r}=0.33,0.5,1}$ and ${10}$
(Figure 9). In Figure 10, the dynamics of the visco-resistive tearing mode
growth are divided in 5-stages. The first stage is the initial transient stage
when the linear instability starts to grow. The second stage is the ${FKR}$
stage, during which both the reconnection rate and the magnetic island width
grow exponentially. In the third stage, the so-called Rutherford stage, the
island evolves toward saturation and subsequent decay. The X-point collapse
and the Y-type, SP-like current sheet forms during the fourth stage (Figure 11
(a)-(c)). The transition from the X-type geometry to the Y-type current sheet
is known as the secondary instability. The first peak that appears in the
kinetic energy evolution is due to the X-point collapse with the onset of
secondary instability at t = 160 sec as shown in Figure 9. After the X-point
collapse, the SP-like current sheet starts to become elongated in the poloidal
direction. During the fifth stage, the Y-type current sheet starts to become
more elongated and SP-like. After the collapse of SP-like current sheet, the
secondary island appears along with two X-points on both ends (Figure 11(e)).
After t = 160 sec, the growth rate starts to decrease up to t = 170 sec, and a
significant change in the growth rate occurs due to the collapse of the Y-type
current sheet and the formation of small plasmoid chain. The size of the
secondary island increases up to some extent and both X-points collapse to
form two tertiary current sheets with the passage of time. The second large
peaks in the growth rate and the kinetic energy plots represent the collapse
of these two tertiary current sheets. As the width of the secondary island
approaches some critical value, the ejection of the secondary island takes
place. At this point kinetic energy increases abruptly and the secondary
island coalesces with the primary one (Figure 12(d)). As the Y-type current
sheet collapse and the plasmoid instability (PI) appears, a drastic increase
in the growth rate takes place (Figure 10) which is much larger as compared to
the growth rate of the secondary instability.

The nonlinear stages of our simulation results described above are quite
similar to the nonlinear secondary island evolution reported by N. F. Loureiro
et al. [2]. But the collapse of the secondary island and direction of ejection
are different in our cases. For example in our simulations for the
${P_{r}=0.5}$ and the ${P_{r}=5}$ cases, the direction of plasmoid ejection is
upward, whereas for the ${P_{r}=1}$ and the ${P_{r}=10}$ cases the direction
is downward. The directions of plasmoid ejection are different for different
${P_{r}}$ cases.

Our simulations for various ${P_{r}}$ numbers show viscous effects on both the
timing and the spatial structure of plasmoid instability. In the ${P_{r}}$ =
0.33 (i.e. very low viscosity), we find secondary islands that saturate early
at smaller size (Figure 11(f)). As we increase the viscosity further
(${P_{r}}$ = 0.5), the appearance and ejection of the secondary island becomes
more prominent, along with the larger saturated island and plasmoid (Figure
12). A clear transition occurs at ${P_{r}=1}$, when the size of the primary
plasmoid becomes the maximum, and the onsets of the secondary instability, the
PI and the island saturation are significantly postponed. The direction of
plasmoid ejection also switches to the opposite. The second peak of kinetic
energy increases with viscosity and reaches maximum at ${P_{r}=1}$ too. For
higher Prandtl number (${P_{r}}$ = 10, with ${\eta=2.8\times 10^{-4}}$ and
viscosity = 0.0028) the size of the plasmoid becomes smaller (Figure 14).

Figure 15 summarizes the relationships among plasmoid width, saturated island
width, and ${P_{r}}$ number. At lower ${P_{r}}$, the width of saturated island
is small, but as we increase viscosity, the width of saturated island
increases sharply up to ${P_{r}=1}$, beyond which the saturated island width
becomes almost independent of the ${P_{r}}$ number. The ${P_{r}=1}$ number
also separates two regimes for the plasmoid width. In the ${P_{r}<1}$ regime,
the width of plasmoid increases drastically with viscosity, whereas in the
${P_{r}>1}$ regime, the width of plasmoid slowly decreases with the viscosity.

##  4 Summary

The key objective of this study is to explore the viscous effects on secondary
instability, plasmoid formation, their merging and ejection process during the
nonlinear evolution of a resistive tearing mode in the large
${\Delta^{\prime}}$ regime. For our equilibrium, we find the critical
instability parameter for the onset of the secondary instability, and the
minimum value of ${\Delta^{\prime}}$ at which PI can take place. Two
distinctive regimes of the ${P_{r}}$ number are found for the plasmoid
instability, which are separated by the value of ${P_{r}=1}$. In the
${P_{r}\lesssim 1}$ regime, the amplitude of the second peak for kinetic
energy increases up to ${P_{r}=1}$, whereas in the ${P_{r}\gtrsim 1}$ regime,
this amplitude decreases with the viscosity. Both the saturated island width
and the plasmoid size increase sharply with the viscosity in the
${P_{r}\lesssim 1}$ regime, however, the former slowly increases whereas the
later decreases with the viscosity in the ${P_{r}\gtrsim 1}$ regime. In
another word, the plasmoid size reaches maximum at ${P_{r}\simeq 1}$. We plan
to further explore the significance of such a finding in future work.

This research was supported by the National Magnetic Confinement Fusion
Science Program of China (No.2019YFE03050004), National Natural Science
Foundation of China (Nos.11775221 and 51821005), U.S.DOE (Nos.DE-
FG02-86ER53218 and DESC0018001), and the Fundamental Research Funds for the
Central Universities at Huazhong University of Science and Technology
(No.2019kfyXJJS193). We are grateful for the support from NIMROD team. This
research used the computing resources from the Supercomputing Center of
University of Science and Technology of China. The author Nisar Ahmad
acknowledges the support from the Chinese Government Scholarship.

## References

## References

  * [1] Marala F et al 1992 Phys. Fluids 4 3072
  * [2] Loureiro N F et al 2005 Phys. Rev. Lett. 95 235003
  * [3] Lapenta G 2008 Phys. Rev. Lett. 100 235001
  * [4] Bhattacharjee A et al 2009 Phys. Plasmas 16 112102
  * [5] Loureiro N F et al 2007 Phys. Plasmas 14 100703
  * [6] Huang Y-M and Bhattacharjee A 2010 Phys. Plasmas 17 062104
  * [7] Huang Y-M and Bhattacharjee A 2013 Phys. Plasmas 20 055702
  * [8] Huang Y-M et al 2011 Phys. Plasmas 18 072109
  * [9] Loureiro N F et al 2013 Phys. Rev. E 87 013102
  * [10] Ali A et al 2014 Phys. Plasmas 21 052312
  * [11] Uzdensky D A and Loureiro N F 2016 Phys. Rev. Lett. 116 105003
  * [12] Loureiro N F et al 2012 Phys. Plasmas 19 042303
  * [13] Dobrott D et al 1977 Phys. Fluids 20 1850
  * [14] Bulanov S V et al 1978 J. Exp. Theor. Phys. Lett. 28 177
  * [15] Bulanov S V et al 1979 Sov. J. Plasma Phys. 5 157
  * [16] Biskamp D 1986 Phys. Fluids 29 1520
  * [17] Tenerani A et al 2015b Astrophys. J. Lett. 813 L32
  * [18] Ali A et al 2015 Phys. Plasmas 22 042102
  * [19] Takeda K et al 2008 Phys. Plasmas 15 022502
  * [20] Bondeson M and Sobel J R 1984 Phys. Fluids 27 2028
  * [21] Porcelli F 1987 Phys. Fluids 30 1734
  * [22] Ofman L et al 1991 Phys. Fluids B 3(6) 1364
  * [23] Grasso D et al 2008 Phys. Plasmas 15 072113
  * [24] Militello F et al 2011 Phys. Plasmas 18 112108
  * [25] Tenerani A et al 2015a Astrophys. J. 801 1
  * [26] Betar H et al 2020 Phys. Plasmas 27 102106
  * [27] Comisso L and Grasso D 2016 Phys. Plasmas 23 032111
  * [28] Comisso L et al 2017 Astrophys. J. 850 16
  * [29] Sovinec R C et al 2004 J. Comput. Phys. 195 355
  * [30] Harris E G 1962 Nuo. Sim. 23 115
  * [31] Tenerani A et al 2016 J. Plasma Phys. 5 82
  * [32] Furth H P et al 1963 Phys. Fluids 6 459
  * [33] Coppi B et al 1976 Sov. J. Plasma Phys. 2 533
  * [34] Waelbroeck F L 1993 Phys. Rev. Lett. 70 3259
  * [35] Biskamp D 2000 Magnetic Reconnection in Plasmas (Cambridge University Press)
  * [36] Uzdensky D et al 2010 Phys. Rev. Lett. 105 235002