# Construction of explicit symplectic integrators in general relativity. II.
Reissner-Nordström black holes

Ying Wang1,2, Wei Sun1, Fuyao Liu1, Xin Wu1,2,3,† 1\. School of Mathematics,
Physics and Statistics, Shanghai University of Engineering Science, Shanghai
201620, China
2\. Center of Application and Research of Computational Physics, Shanghai
University of Engineering Science, Shanghai 201620, China
3\. Guangxi Key Laboratory for Relativistic Astrophysics, Guangxi University,
Nanning 530004, China  Emails<EMAIL_ADDRESS>(Y. W.),
<EMAIL_ADDRESS>(W. S<EMAIL_ADDRESS>(F. L.); ${\dagger}$
Corresponding Author<EMAIL_ADDRESS>(X. W.)

###### Abstract

In a previous paper, second- and fourth-order explicit symplectic integrators
were designed for a Hamiltonian of the Schwarzschild black hole. Following
this work, we continue to trace the possibility of the construction of
explicit symplectic integrators for a Hamiltonian of charged particles moving
around a Reissner-Nordström black hole with an external magnetic field. Such
explicit symplectic methods are still available when the Hamiltonian is
separated into five independently integrable parts with analytical solutions
as explicit functions of proper time. Numerical tests show that the proposed
algorithms share the desirable properties in their long-term stability,
precision and efficiency for appropriate choices of step sizes. For the
applicability of one of the new algorithms, the effects of the black hole’s
charge, the Coulomb part of the electromagnetic potential and the magnetic
parameter on the dynamical behavior are surveyed. Under some circumstances,
the extent of chaos gets strong with an increase of the magnetic parameter
from a global phase-space structure. No the variation of the black hole’s
charge but the variation of the Coulomb part is considerably sensitive to
affect the regular and chaotic dynamics of particles’ orbits. A positive
Coulomb part is easier to induce chaos than a negative one.

_Unified Astronomy Thesaurus concepts_ : Black hole physics (159);
Computational methods (1965); Computational astronomy (293); Celestial
mechanics (211)

##  1 Introduction

Based on Einstein’s theory of general relativity, black holes are a group of
solutions of the Einstein’s field equations. Usually, black holes have
singularities covered by event horizon surfaces. The first black hole’s
solution for the description of a static and spherically symmetric
gravitational field around a point-like mass was given by Schwarzschild
(1916). When the mass source in the origin is charged, the Reissner-Nordström
(RN) metric (Reissner 1916) became available. A rigorous solution representing
the gravitational field around a rotating central mass is the Kerr metric
(Kerr 1963). On the other hand, evidence from observations demonstrated the
existence of supermassive black holes with masses from millions to tens of
billions of solar masses in centers of nearly all galaxies. In terms of images
of M87, a central Kerr black hole is estimated to have mass $M=(6.5\pm
0.7)\times 10^{9}M_{\odot}$ ($M_{\odot}$ being the Sun’s mass), which is
consistent with the result predicted by the general theory of relativity (EHT
Collaboration et al. 2019a, b, c). Successful gravitational-wave measurements
(Abbott et al. 2016) also provided powerful evidence for the presence of black
holes.

Although the relativistic spacetimes like the RN or Kerr metric are highly
nonlinear, they are integrable and have analytical solutions because of the
presence of enough constants of motion. The solutions have only formal
expressions in terms of quadratures, but cannot be expressed as elementary
functions or explicit functions of time. To know detailed information on the
solutions how to evolve with time, one had better employ a numerical
integration technique to solve the integrable problems. If the central bodies
are suffered from perturbations, such as external magnetic fields, the
spacetimes become non-integrable in most cases. Under some circumstances,
chaos occurs (Takahashi $\&$ Koyama 2009; Kopáček et al. 2010; Kopáček $\&$
Karas 2014; Stuchlík $\&$ Kološ 2016; Pánis et al. 2019; Li $\&$ Wu 2019;
Stuchlík et al. 2020; Yi $\&$ Wu 2020). This chaoticity indicates that a
dynamical system is exponentially sensitive dependence on initial conditions
(Lichtenberg $\&$ Lieberman 1983). In this case, the numerical technique is
more indispensable to study the non-integrable systems.

Reliable results from numerical integrators with a good behavior are required,
especially in the case of long-term integration of chaotic orbits. The most
appropriate solvers are geometric or structure preserving algorithms (Hairer
et al. 1999; Seyrich $\&$ Lukes-Gerakopoulos 2012; Bacchini et al. 2018a,
2018b; Hu et al. 2019), such as symplectic methods for Hamiltonian systems
(Ruth 1983; Wisdom $\&$ Holman 1991). They have several advantages over
standard explicit integrators, such as the family of explicit Runge-Kutta
solvers. The integrals of motion (e.g., energy integral) along the trajectory
are nearly conserved for the structure preserving integrators, but their
errors increase linearly with time for the standard integration schemes. In
addition, the overall phase error only grows linearly with time for the former
algorithms, whereas it is normally proportional to the square of the length of
the integration interval for the latter schemes (Deng et al. 2020). The above-
mentioned curved spacetimes can be expressed in terms of Hamiltonian systems,
and thus symplectic integrators are naturally chosen. Standard explicit
symplectic integrators, such as a second-order Verlet integrator (Swope et al.
1982), become useless. However, completely implicit symplectic methods
(Kopáček et al. 2010; Seyrich $\&$ Lukes-Gerakopoulos 2012; Tsang et al.
2015), such as the implicit midpoint scheme (Feng 1986; Brown 2006), or
implicit and explicit combined symplectic methods (Liao 1997; Preto $\&$ Saha
2009; Lubich et al. 2010; Zhong et al. 2010; Mei et al. 2013a, 2013b) are
used. This is because the Hamiltonians have no separable forms of variables or
can be split into two integrable parts without analytical solutions as
explicit functions of time. Unfortunately, such implicit integrators are more
computationally demanding at the expense of computational time than the same
order standard explicit methods. Extended phase-space methods (Pihajoki 2015;
Liu et al. 2016; Luo et al. 2017; Li $\&$ Wu 2017) are explicit and have good
long-term stable behavior in energy errors, but are not symplectic. In a
previous work (Wang et al. 2021), we overcame the difficulty in the
construction of explicit symplectic integrators for the Schwarzschild metric.
In our construction, the Hamiltonian for the Schwarzschild spacetime can be
separated into four integrable parts with analytical solutions as explicit
functions of proper time. Then, these explicit solvable operators
symmetrically composed second- and fourth-order explicit symplectic
integrators.

Following the previous work (Wang et al. 2021), we design explicit symplectic
integrators for the RN black hole immersed into an external magnetic field.
This is one of the main aims in the present paper. Another aim is to know how
an increase of the black hole’s charge, the Coulomb paramter of the
electromagnetic potential or the magnetic parameter exerts an influence on the
dynamical transition of orbits of charged particles around the RN black hole.
For the sake of these purposes, we introduce a dynamical model of charged
particles moving around the RN black hole surrounded with an external magnetic
field in Section 2. Second- and fourth-order explicit symplectic integrators
are designed for the magnetized RN spacetime in Section 3. In Section 4, we
evaluate the numerical performance of the proposed algorithms, and apply a new
integrator to address the question of how the related parameters affect the
orbital dynamics of order and chaos. Finally, the main results are concluded
in Section 5.

##  2 Reissner-Nordström black holes

The Schwarzschild black hole with charge $Q$ is the RN black hole. In
dimensionless spherical-like coordinates $(t,r,\theta,\phi)$, the RN spacetime
(Reissner 1916) takes the following metric

$\displaystyle-\tau^{2}$ $\displaystyle=$ $\displaystyle
ds^{2}=g_{\alpha\beta}dx^{\alpha}dx^{\beta}$ (1) $\displaystyle=$
$\displaystyle-(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})dt^{2}+(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})^{-1}dr^{2}$
$\displaystyle+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\phi^{2}.$

The speed of light $c$ and the constant of gravity $G$ use geometrized units,
$c=G=1$. $M$ is the mass of black hole, and also takes one unit, $M=1$. Proper
time $\tau$, coordinate time $t$, radial separation $r$ and charge $Q$ are
dimensionless. In practice, the dimensionless operations are obtained via
scale transformations: $\tau\rightarrow M\tau$, $t\rightarrow Mt$,
$r\rightarrow Mr$ and $Q\rightarrow MQ$. When $|Q|<1$, this spacetime
corresponds to black holes with two event horizons
$r_{\pm}=1\pm\sqrt{1-Q^{2}}$. The spacetime is still a black hole with an
event horizon $r=1$ for $Q=\pm 1$. It has no event horizon but has naked
singularities if $|Q|>1$. Hereafter, the case of black holes with $|Q|\leq 1$
is considered.

The motion of a test particle around the black hole is described by the
Lagrangian system

$\ell=\frac{1}{2}(\frac{ds}{d\tau})^{2}=\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu},$
(2)

where $\dot{x}^{\mu}=\mathbf{U}$ is a four-velocity satisfying the relation

$\mathbf{U}\cdot\mathbf{U}=U^{\alpha}U_{\alpha}=g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}=-1.$
(3)

A covariant generalized momentum $\mathbf{p}$ is defined as

$p_{\mu}=\frac{\partial\ell}{\partial\dot{x}^{\mu}}=g_{\mu\nu}\dot{x}^{\nu}.$
(4)

Because $t$ and $\phi$ do not explicitly appear in the Lagrangian, there are
two constant momentum components

$\displaystyle p_{t}$ $\displaystyle=$
$\displaystyle-(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})\dot{t}=-\mathcal{E},$ (5)
$\displaystyle p_{\phi}$ $\displaystyle=$ $\displaystyle
r^{2}\sin^{2}\theta\dot{\phi}=\mathfrak{L}.$ (6)

$\mathcal{E}$ and $\mathfrak{L}$ denote the particle’s energy and angular
momentum, respectively.

The Lagrangian corresponds to the Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\mathbf{U}\cdot\mathbf{p}-\ell=\frac{1}{2}g^{\mu\nu}p_{\mu}p_{\nu}$
(7) $\displaystyle=$
$\displaystyle-\frac{\mathcal{E}^{2}}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})^{-1}+\frac{p^{2}_{r}}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})$
$\displaystyle+\frac{1}{2}\frac{p^{2}_{\theta}}{r^{2}}+\frac{1}{2}\frac{\mathfrak{L}^{2}}{r^{2}\sin^{2}\theta}.$

Because of the four-velocity relation (3), this Hamiltonian is always
identical to -1/2,

$H=-\frac{1}{2}.$ (8)

By separating the variables in the Hamilton-Jacobi equation, one can find a
second integral excluding the two integrals $\mathcal{E}$ and $\mathfrak{L}$
in the Hamiltonian system (Carter 1968). Thus, this system is integrable and
has formally analytical solutions.

Now, suppose the black hole surrounded by an external magnetic field whose
four-vector potential has two nonzero covariant components

$A_{t}=-\frac{Q}{r},~{}~{}~{}A_{\phi}=\frac{B}{2}g_{\phi\phi}=\frac{B}{2}r^{2}\sin^{2}\theta,$
(9)

where $A_{t}$ represents the Coulomb part of the electromagnetic potential
(Kopaccek $\&$ Karas 2014), and $B$ is the strength of the magnetic field
parallel to the $z$ axis (Felice $\&$ Sorge 2003). The motion of a particle
with charge $q$ under the interactions of the black hole’s gravity and
electromagnetic force is described by the Hamiltonian

$\displaystyle K$ $\displaystyle=$
$\displaystyle\frac{1}{2}g^{\mu\nu}(p_{\mu}-qA_{\mu})(p_{\nu}-qA_{\nu})$ (10)
$\displaystyle=$
$\displaystyle-\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})^{-1}(E-\frac{Q^{*}}{r})^{2}$
$\displaystyle+\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})p^{2}_{r}+\frac{1}{2}\frac{p^{2}_{\theta}}{r^{2}}$
$\displaystyle+\frac{1}{2r^{2}\sin^{2}\theta}(L-\frac{1}{2}\beta
r^{2}\sin^{2}\theta)^{2},$

where $Q^{*}=qQ$ is a Coulomb parameter of the electromagnetic potential, and
$\beta=qB$. To make the system (10) be dimensionless, we take $K\rightarrow
m^{2}K$, $E\rightarrow mE$, $p_{r}\rightarrow mp_{r}$, $p_{\theta}\rightarrow
mMp_{\theta}$, $L\rightarrow mML$, $q\rightarrow mq$ and $B\rightarrow B/M$,
where $m$ is the particle’s mass. The expressions of energy $E$ and angular
momentum $L$ of the charged particle become

$\displaystyle E$ $\displaystyle=$
$\displaystyle(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})\dot{t}+\frac{Q^{*}}{r},$
(11) $\displaystyle L$ $\displaystyle=$ $\displaystyle
r^{2}\sin^{2}\theta\dot{\phi}+\frac{1}{2}\beta r^{2}\sin^{2}\theta.$ (12)

Similar to $H$, $K$ always satisfies the constraint

$K=-\frac{1}{2}.$ (13)

However, $K$ unlike $H$ has no second integral. Thus, it is non-integrable and
has no formally analytical solutions. In this case, a numerical integration
method is a convenient tool to work out such a non-integrable system.

##  3 Construction of explicit symplectic integrators

In view of a symplectic integrator with good geometric and physical
properties, it is naturally a prior choice of numerical integrator for the
description of long-term qualitative evolution of the Hamiltonian system (7).
An explicit symplectic method becomes useless without doubt if this
Hamiltonian is separated into two analytically solvable parts. This is because
not all analytical solutions of the two splitting parts are explicit functions
of proper time $\tau$. An explicit symplectic algorithm fails to be built if
the Hamiltonian is split into four analytically integrable parts, as the
Hamiltonian of Schwarzschild black hole is in our previous paper (Wang et al.
2021). Thus, it may be necessary that more analytically integrable splitting
parts should be given to the Hamiltonian for the construction of explicit
symplectic schemes.

Let the Hamiltonian of RN black hole be separated into five separable parts

$H=H_{1}+H_{2}+H_{3}+H_{4}+H_{5},$ (14)

where the five sub-Hamiltonians are written as follows:

$\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{\mathfrak{L}^{2}}{r^{2}\sin^{2}\theta}-\frac{\mathcal{E}^{2}}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})^{-1},$
(15) $\displaystyle H_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}p^{2}_{r},$ (16) $\displaystyle H_{3}$
$\displaystyle=$ $\displaystyle-\frac{1}{r}p^{2}_{r},$ (17) $\displaystyle
H_{4}$ $\displaystyle=$ $\displaystyle\frac{p^{2}_{\theta}}{2r^{2}},$ (18)
$\displaystyle H_{5}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{Q^{2}}{r^{2}}p^{2}_{r}.$ (19)

$H_{2}$, $H_{3}$ and $H_{4}$ are the same as those in the Hamiltonian
splitting of Schwarzschild black hole in the previous work (Wang et al. 2021).

$H_{1}$ has its canonical equations $\dot{r}=\dot{\theta}=0$ and

$\displaystyle\frac{dp_{r}}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{\partial H_{1}}{\partial
r}=\frac{\mathfrak{L}^{2}}{r^{3}\sin^{2}\theta}-\frac{\mathcal{E}^{2}}{r^{2}}(1-\frac{Q^{2}}{r})$
(20)
$\displaystyle\cdot(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})^{-2}=\Re(r,\theta),$
$\displaystyle\frac{dp_{\theta}}{d\tau}$ $\displaystyle=$
$\displaystyle-\frac{\partial
H_{1}}{\partial\theta}=\frac{\mathfrak{L}^{2}\cos\theta}{r^{2}\sin^{3}\theta}=\Theta(r,\theta).$
(21)

If $\mathcal{A}$ is taken as a differential operator

$\displaystyle\mathcal{A}=\Re\frac{\partial}{\partial
p_{r}}+\Theta\frac{\partial}{\partial p_{\theta}},$ (22)

then $\mathcal{A}p_{r}=\dot{p}_{r}=\Re$ and
$\mathcal{A}p_{\theta}=\dot{p}_{\theta}=\Theta$. Because $r$ and $\theta$ are
constants, $p_{r}$ and $p_{\theta}$ are easily solved. From proper time
$\tau_{0}$ over a proper time step $h$ to proper time $\tau=\tau_{0}+h$, the
solutions are expressed as

$\displaystyle p_{r}$ $\displaystyle=$ $\displaystyle
p_{r0}+h\Re(r_{0},\theta_{0}),$ (23) $\displaystyle p_{\theta}$
$\displaystyle=$ $\displaystyle p_{\theta 0}+h\Theta(r_{0},\theta_{0}),$ (24)

where $\textbf{z}(0)=(r_{0},\theta_{0},p_{r0},p_{\theta 0})$ are the solutions
at the beginning of the step of length $h$. We use an exponential operator
$e^{h\mathcal{A}}$ to represent the analytical solutions (23) and (24), i.e.,
$(p_{r},p_{\theta})=e^{h\mathcal{A}}\textbf{z}(0)$.

Set $\mathcal{B}$, $\mathcal{C}$, $\mathcal{D}$ and $\mathcal{F}$ as
differential operators of $H_{2}$, $H_{3}$, $H_{4}$ and $H_{5}$, respectively.
They are of the following expressions

$\displaystyle\mathcal{B}$ $\displaystyle=$ $\displaystyle
p_{r}\frac{\partial}{\partial r},$ (25) $\displaystyle\mathcal{C}$
$\displaystyle=$ $\displaystyle-\frac{2}{r}p_{r}\frac{\partial}{\partial
r}-\frac{p^{2}_{r}}{r^{2}}\frac{\partial}{\partial p_{r}},$ (26)
$\displaystyle\mathcal{D}$ $\displaystyle=$
$\displaystyle\frac{p_{\theta}}{r^{2}}\frac{\partial}{\partial\theta}+\frac{p^{2}_{\theta}}{r^{3}}\frac{\partial}{\partial
p_{r}},$ (27) $\displaystyle\mathcal{F}$ $\displaystyle=$
$\displaystyle\frac{Q^{2}}{r^{2}}p_{r}\frac{\partial}{\partial
r}+Q^{2}\frac{p^{2}_{r}}{r^{3}}\frac{\partial}{\partial p_{r}}.$ (28)

The four sub-Hamiltonians have their analytical solutions

$\displaystyle e^{h\mathcal{B}}:~{}r$ $\displaystyle=$ $\displaystyle
r_{0}+hp_{r0};$ (29) $\displaystyle e^{h\mathcal{C}}:~{}r$ $\displaystyle=$
$\displaystyle[(r^{2}_{0}-3hp_{r0})^{2}/r_{0}]^{1/3},$ $\displaystyle p_{r}$
$\displaystyle=$ $\displaystyle
p_{r0}[(r^{2}_{0}-3hp_{r0})^{2}/r^{2}_{0}]^{1/3};$ (30) $\displaystyle
e^{h\mathcal{D}}:~{}\theta$ $\displaystyle=$
$\displaystyle\theta_{0}+hp_{\theta 0}/r^{2}_{0},$ $\displaystyle p_{r}$
$\displaystyle=$ $\displaystyle p_{r0}+hp^{2}_{\theta 0}/r^{3}_{0};$ (31)
$\displaystyle e^{h\mathcal{F}}:~{}r$ $\displaystyle=$
$\displaystyle\sqrt{r^{2}_{0}+2hQ^{2}p_{r0}/r_{0}},$ $\displaystyle p_{r}$
$\displaystyle=$
$\displaystyle\frac{p_{r0}}{r_{0}}\sqrt{r^{2}_{0}+2hQ^{2}p_{r0}/r_{0}}.$ (32)

It is clear that all the analytical solutions in Equations (23), (24) and
(29)-(32) are explicit functions of proper time $\tau$ or step size $h$.
Although the compositions $H_{2}+H_{3}$, $H_{2}+H_{3}+H_{4}$ and
$H_{2}+H_{3}+H_{4}+H_{5}$ can be solved analytically, their solutions are not
expressed in terms of explicit functions of $\tau$. The present splitting form
(14) of the Hamiltonian $H$ is one possible choice to satisfy the need.

The solutions of the Hamiltonian (7) over the time step $h$ can be obtained
approximately by a second order explicit symplectic integrator, namely,
symmetric products of these exponential operators

$\displaystyle S^{H}_{2}(h)$ $\displaystyle=$ $\displaystyle
e^{\frac{h}{2}\mathcal{F}}e^{\frac{h}{2}\mathcal{D}}e^{\frac{h}{2}\mathcal{C}}e^{\frac{h}{2}\mathcal{B}}e^{h\mathcal{A}}$
(33) $\displaystyle\otimes
e^{\frac{h}{2}\mathcal{B}}e^{\frac{h}{2}\mathcal{C}}e^{\frac{h}{2}\mathcal{D}}e^{\frac{h}{2}\mathcal{F}}.$

It can compose a fourth-order symplectic scheme of Yoshida (1990)

$S^{H}_{4}(h)=S^{H}_{2}(\gamma h)\circ S^{H}_{2}(\delta h)\circ
S^{H}_{2}(\gamma h),$ (34)

where $\delta=1-2\gamma$ and $\gamma=1/(2-\sqrt[3]{2})$.

The two explicit symplectic algorithms for $H$ are also suitable for $K$. The
only one difference is $H_{1}$ replaced with

$\displaystyle K_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2r^{2}\sin^{2}\theta}(L-\frac{1}{2}\beta
r^{2}\sin^{2}\theta)^{2}$ (35)
$\displaystyle-\frac{1}{2}(1-\frac{2}{r}+\frac{Q^{2}}{r^{2}})^{-1}(E-\frac{Q^{*}}{r})^{2}.$

Then, we obtain two explicit symplectic methods $S^{K}_{2}$ and $S^{K}_{4}$
for the Hamiltonian $K$.

Table 1: Performance of algorithms S2 and S4 with different time steps $h$. In the brackets, e.g., ($10^{-8}$, U, $11^{\prime}40^{\prime\prime}$), $10^{-8}$ denotes the order of Hamiltonian error, U (or B) indicates the unboundedness (or boundedness) of Hamiltonian error, and $11^{\prime}40^{\prime\prime}$ corresponds to CPU time (minute ′, second ′′). The integration time reaches $\tau=10^{8}$ for each step size. h | 0.1 | 1 | 4 | 10
---|---|---|---|---
S2 | ($10^{-8}$, U, $11^{\prime}40^{\prime\prime}$) | ($10^{-6}$, B, $1^{\prime}12^{\prime\prime}$) | ($10^{-5}$, B, $20^{\prime\prime}$) | ($10^{-4}$, B, $7^{\prime\prime}$)
S4 | no tested | ($10^{-9}$, U, $3^{\prime}34^{\prime\prime}$) | ($10^{-8}$, B, $56^{\prime\prime}$) | ($10^{-6}$, B, $23^{\prime\prime}$)

##  4 Numerical simulations

At first, let us check the numerical performance of the proposed explicit
symplectic integration algorithms for solving the system (10). Then, one of
the new methods is selected to explore the orbital dynamics of charged
massless particles in the system.

###  4.1 Evaluations of the new algorithms

In the previous work (Wang et al. 2021), the established explicit symplectic
integrators for the Schwarzschild black hole surrounded by an external
magnetic field were compared with a conventional fourth-order Runge-Kutta
integrator, second- and fourth-order explicit and implicit mixed symplectic
algorithms (Mei et al. 2013b) and second- and fourth-order extended phase-
space explicit symplectic-like methods (Luo et al. 2017). It was shown that
the Runge-Kutta method has a secular drift in Hamiltonian errors and performs
the poorest performance. The other algorithms at same order can exhibit good
long-term stable error behavior for appropriate time steps and have no
explicit differences among their Hamiltonian errors. Therefore, only the newly
proposed explicit symplectic integrators in the present paper are considered
to work out the Hamiltonian $K$.

Taking proper time step $h=1$, we consider the parameters to be $E=0.995$,
$L=4.6$, $\beta=6.4\times 10^{-4}$, $Q=0.1$ and $Q^{*}=10^{-4}$. The new
second-order explicit symplectic integrator S2 (or the new fourth-order method
S4) is used to integrate an orbit with initial conditions $r=25$, $p_{r}=0$
and $\theta=\pi/2$. The starting value of $p_{\theta}>0$ is determined by
Equation (13). In Figure 1(a), Hamiltonian errors $\Delta K=-1/2-K$ in
Equation (13) can remain bounded in an order of $\mathcal{O}(10^{-6})$ for S2
when the number of integration steps is $10^{8}$. S4 gives a higher accuracy
with an order of $\mathcal{O}(10^{-9})$, but its errors grow linearly with
time due to roundoff errors. Here are some analysis to these results. The test
orbit has an approximate average period $T\approx 5000$. Truncation
Hamiltonian error is $(h/T)^{2}\sim 10^{-8}$ for S2 and $(h/T)^{4}\sim
10^{-16}$ for S4. In fact, the error outputted at the end of the first step is
$1.98\times 10^{-15}$ for S4. In addition, the machine yields a roundoff error
in per computation, e.g. $\epsilon=10^{-16}$ in a double-precision level. The
roundoff errors grow in a rough estimation $n\epsilon$, where $n$ is a number
of computations. The roundoff errors are more important than the truncation
errors when $n$ is large enough. In an integration time $t=200$, the
Hamiltonian errors for S4 fast grow to $9\times 10^{-12}$. When $t$ spans this
time and is less than $10^{7}$, the errors much slowly grow and basically
remain stable at an order of $\mathcal{O}(10^{-9})$. With the integration
continuing, the boundness of the Hamiltonian errors is destroyed by the
roundoff errors. If the step size gets larger, e.g. $h$=10, the Hamiltonian
errors are stabilized at an order of $\mathcal{O}(10^{-6})$. For $h$=4, the
Hamiltonian errors remain bounded in an order of $\mathcal{O}(10^{-8})$. These
results roughly indicate that a symplectic integrator can stabilize the
Hamiltonian errors at the values larger than $\mathcal{O}(10^{-8})$ for
$10^{8}$ integration steps, which yield roundoff errors in the order of
$\mathcal{O}(10^{-8})$. When the time step is $h=1$ and the number of
integration steps is $10^{8}$, the main error source for S2 is the truncation
errors and therefore the Hamiltonian errors can remain stable at the order of
$\mathcal{O}(10^{-6})$. However, the roundoff errors for S2 with time step
$h=0.1$ reach an order of $\mathcal{O}(10^{-7})$ after $10^{9}$ integration
steps. This forces the Hamiltonian errors with an order of
$\mathcal{O}(10^{-8})$ to grow linearly. To clearly show the dependence of the
magnitude and boundness of Hamiltonian errors and computational efficiency for
algorithms S2 and S4 on the time step $h$, we give Table 1. Obviously, S2 with
$h=1$ or S4 with $h=4$ is an optimal choice in the present cases.

The test orbit in Figure 1(a) is Orbit 1 colored red in Figure 1(b). Because
this orbit is a single Kolmogorov-Arnold-Moser (KAM) torus on the Poincaré
section, it is a regular quasi-periodic orbit. Black Orbit 2 with the initial
value $r=90$, consisting of 11 small islands, is also a regular many-islands
KAM torus. The occurrence of resonance and chaos will become easy for such an
orbit with many islands. However, blue Orbit 3 with the initial value $r=50$
has many discrete points distributed in a small area on the Poincaré section.
This kind of phase space structure indicates the chaoticity of Orbit 3. The
purple orbit with the initial value $r=110$ is also chaotic. In spite of the
onset of chaos, the possibility of “islands of regularity” like Orbits 1 and 2
is still existent. Based on KAM theorem, the minima of the effective potential
in the equatorial plane correspond to stable circular orbits, which are
related to regular harmonic oscillatory motions for the description of
Keplerian accretion disks of stellar mass black holes. Moreover, all
trajectories that are bounded in the vicinity of the equatorial plane are also
regular (Kološ et al. 2015). These regular motions can successfully explain
the quasi-periodic oscillations of X-ray flux from several microquasars (Kološ
et al. 2017; Tursunov $\&$ Kološ 2018). There is another island of regularity
related to the motion along the magnetic field lines (Tursunov et al. 2020a).
In some cases, a hot spot can exhibit quasi-circular motion along a single
orbit (Tursunov et al. 2020b). However, the quasicircular motion may become
chaotic because the axial symmetry of the system is broken so that the hot
spot’s angular momentum is not conserved. The absence of the angular momentum
is caused by the inclination angle of the hot-spot orbit from the equatorial
plane or of the magnetic field lines with respect to the black hole’s spin
axis.

If regular single-island Orbit 1 is replaced with regular many-islands Orbit 2
or chaotic Orbit 3, the numerical performance of the two algorithms S2 and S4
has no explicit differences. In other words, no dynamical behavior of orbits
but a step size mainly affects the quality of the proposed algorithms. In the
later discussions, we employ S4 with the appropriate time step $h=4$ to
investigate the related dynamical features of the Hamiltonian $K$ when charge
parameters $Q$ and $Q^{*}$, and magnetic parameter $\beta$ are varied.

###  4.2 Applications

To show the dependence of the orbital dynamics of order and chaos on the black
hole’s charge $Q$ or the Coulomb parameter $Q^{*}$, we fix the parameters $E$,
$L$ and $\beta$ in Figure 1. Of course, different values of $Q$ and $Q^{*}$
are adopted.

In fact, the phase-space structures for the case of $Q=Q^{*}=0$ are similar to
those for the case of $Q=0.1$ and $Q^{*}=10^{-4}$ in Figure 1(b). To clearly
show how the orbital dynamics of order and chaos depends on the charge $Q$, we
consider the choice of $Q\gg Q^{*}$. Fixing $Q^{*}=10^{-4}$, we give $Q$
different larger values. The result for $Q=0.1$ is also suitable for the case
of $Q=0.3$. However, the phase-space structures for $Q=0.6$ in Figure 2(a) are
somewhat different from those for $Q=0.1$. Orbit 2 is ordered in Figure 1(b),
but becomes chaotic in Figure 2(a). Orbit 3 that is chaotic in Figure 1(b) is
a regular orbit with many loops in Figure 2(a). Compared with those for the
case of $Q=0.6$ in Figure 2(a), the orbits exist some differences for the case
of $Q=0.8$ in Figure 2(b). The blue islands in Figure 2(a) become an ordered
single torus in Figure 2(b). The green ordered single torus in Figure 2(a) is
weakly chaotic in Figure 2(b). When $Q=1$, the extent of chaos in Figure 2(c)
is not typically strengthened. A result seems to be concluded from Figures
1(b) and 2(a)-2(c). An increase of the black hole’s charge $Q$ may exert some
influence on the phase-space structures, but does not bring an apparent
dynamical transition from order to chaos. Namely, it is not considerably
sensitive to alter the dynamical orbital properties. It does not typically
enhance the extent of chaos, either.

What about the dynamical transition with an increase of the Coulomb parameter
$Q^{*}$ for a given smaller value (e.g., $Q=10^{-4}$)? Figures 2(d) and 2(e)
describe that the chaotic behavior existing in the case of $Q=Q^{*}=0$
gradually dies out when $Q^{*}$ increases, such as $Q^{*}=$0.1 and 0.3. As $Q$
runs from a smaller value $Q=10^{-4}$ to a larger value $Q=0.3$, the phase-
space structures have no dramatic differences between Figures 2(f) and 2(e).
In fact, the phase-space structures for the case of $Q=Q^{*}=0.1$ are
basically similar to those for the case of $Q=10^{-4}$ and $Q^{*}=0.1$ in
Figure 2(d). The orbits for the case of $Q=0.1$ and $Q^{*}=0.3$ are also the
same as those for the case of $Q=10^{-4}$ and $Q^{*}=0.3$ in Figure 2(e). This
result shows again that no $Q$ but $Q^{*}$ mainly affects the regular and
angular dynamics of orbits. In particular, a positive value of $Q^{*}$ weakens
the strength of chaos. On the other hand, a negative value of $Q^{*}$ can
easily induce the occurrence of chaos, and the extent of chaos is drastically
strengthened when the magnitude of negative Coulomb parameter $Q^{*}$
increases, as shown in Figures 2(g)-2(i). Notice that the orbits between the
case of $Q=0.1$ and $Q^{*}=-0.1$ and the case of $Q=10^{-4}$ and $Q^{*}=-0.1$
in Figure 2(g) are almost the same. So are the orbits between the case of
$Q=0.1$ and $Q^{*}=-0.3$ and the case of $Q=10^{-4}$ and $Q^{*}=-0.3$ in
Figure 2(i). When Coulomb parameter $Q^{*}$ is negative, the Coulomb part of
the electromagnetic potential $A_{t}$ in Equation (9) is positive.

What will happen if magnetic parameter $\beta$ increases but parameters $E$,
$L$, $Q$ and $Q^{*}$ are fixed? Red Orbit 1 with parameters $Q=0.1$ and
$Q^{*}=10^{-4}$ in Figure 1(b) is tested. The orbit is twisted for
$\beta=9.7\times 10^{-4}$ in Figure 3(a), becomes a three-islands orbit for
$\beta=9.9\times 10^{-4}$ in Figure 3(b), and is finally evolved to a strong
chaotic orbit for $\beta=1.1\times 10^{-3}$ in Figure 3(c). Given $Q=10^{-4}$
and $Q^{*}=0.1$, the orbit becomes many-islands, weakly chaotic and strong
chaotic orbits as $\beta$ increases from $9.7\times 10^{-4}$ to $1.1\times
10^{-3}$ in Figures 3(d)-3(f). For $Q=Q^{*}=0.1$, the orbit is evolved to a
twisted single torus, a many-islands orbit and a strong chaotic orbit with an
increase of $\beta$ in Figures 3(g)-3(i). Particular for $Q^{*}=$-0.1, -0.3
and $\beta=9.7\times 10^{-4}$, $9.9\times 10^{-4}$, $1.1\times 10^{-3}$,
strong chaos (not plotted) always occurs. All the results prove that an
increase of magnetic parameter $\beta$ gives rise to enhancing the chaotic
effect.

Why do the two charge parameters $Q$ and $Q^{*}$ have completely different
effects on the dynamical behavior of orbits? Why does an increase of negative
parameter $Q^{*}$ or magnetic parameter $\beta$ lead to strengthening the
extent of chaos? To answer these questions, we expand the term $(~{})^{-1}$ in
Equation (35) and rewrite Equation (35) as follows:

$\displaystyle K_{1}$ $\displaystyle\approx$ $\displaystyle-\frac{1}{2}(\beta
L+E^{2})+\frac{\beta^{2}}{8}r^{2}\sin^{2}\theta-\frac{E^{2}}{r}$ (36)
$\displaystyle+\frac{EQ^{*}}{r}+\frac{L^{2}}{2r^{2}\sin^{2}\theta}+\frac{Q^{2}E^{2}}{2r^{2}}$
$\displaystyle+\frac{Q^{*}}{2r^{2}}(4E-Q^{*})+\cdots.$

The second term $V_{1}=\beta^{2}r^{2}\sin^{2}\theta/8$ in Equation (36) is a
magnetic field force acting as a gravitational effect to the particle. The
third term $V_{2}=-E^{2}/r$ is the gravity of the black hole to the particle.
The Coulomb term $V_{3}=EQ^{*}/r$ acts as a repulsive force effect to the
particle for $Q^{*}>0$, but a gravitational force effect for $Q^{*}<0$. The
fifth term $V_{4}=L^{2}/(2r^{2}\sin^{2}\theta)$ is an inertial centrifugal
force caused by the particle’s angular momentum $L$. The sixth term
$V_{5}=Q^{2}E^{2}/(2r^{2})$ is an electric field force, which acts as a
repulsive force effect to the particle. The magnetic field part $V_{1}$ is a
fundamental source for causing the nonintegrability and chaoticity of the
system (10). For $\beta=0$, the system (10) is integrable and nonchaotic. When
$\beta$ is extremely small in the case of $Q=0$ corresponding to
$V_{3}=\emph{V}_{5}=0$, the black hole’s gravity $V_{2}$ is a dominant force
and therefore chaos does not possibly occur, either. With $\beta$ increasing,
the magnetic field force increases. Only when $V_{1}$ appropriately matches
with $V_{2}$, may chaos occur. The extent of chaos can be strengthened with an
increase of the magnetic parameter from the global phase-space structure. As
to $V_{3}$ and $V_{5}$ to the contributions of particle’s dynamics, $V_{3}\sim
1/r$ is a primary part, and $V_{5}\sim 1/r^{2}$ is a secondary part for $r\gg
2$ and $|Q|\leq 1$. This can explain why the variation of $Q^{*}$ rather than
the variation of $Q$ is considerably sensitive to affect the regular and
chaotic dynamics of particles’ orbits. For $Q^{*}>0$, the Coulomb term, as a
repulsive force, reduces the gravitational effect from the magnetic field. On
the contrary, the Coulomb term, as a gravitational force, enhances the
magnetic field gravitational force effect. Thus, an increase of the magnitude
of negative Coulomb parameter $Q^{*}$ leads to strengthening the extent of
chaos, whereas an increase of positive Coulomb parameter $Q^{*}$ does not.

##  5 Conclusions

In this paper, we are devoted to designing explicit symplectic integrators for
a Hamiltonian system of charged test particles moving around the Reissner-
Nordström black hole immersed into an external magnetic field. In our
construction, the Hamiltonian is separated into five independently integrable
parts with analytical solutions as explicit functions of proper time. These
analytical solutions are used to yield second- and fourth-order explicit
symplectic integrators in symmetric combinations.

The proposed algorithms are shown to exhibit good long term numerical
performance in the Hamiltonian conservation, numerical accuracy and
computational efficiency. Such good numerical performance does not mainly
depend on the regular and chaotic dynamical behavior of orbits but a step
size. Thus, an optimal step size is necessary.

The fourth-order explicit symplectic integrator with an optimal step size is
applied to well explore the dynamics of charged particles around the Reissner-
Nordström black hole with an external magnetic field. We focus on the
influences of the black hole’s charge, the Coulomb part of the electromagnetic
potential and the magnetic parameter on the dynamical behavior. The magnetic
parameter plays an important role in causing the nonintegrability and
chaoticity of the system. Under some circumstances, the extent of chaos is
strengthened from the global phase-space structure as the magnetic parameter
increases. No the variation of the black hole’s charge but the variation of
the Coulomb part is considerably sensitive to affect the regular and chaotic
dynamics of particles’ orbits. A positive Coulomb part is easier to induce
chaos than a negative one.

## Acknowledgments

The authors are very grateful to a referee for useful suggestions. This
research has been supported by the National Natural Science Foundation of
China [Grant Nos. 11533004, 11973020 (C0035736), 11803020, 41807437, U2031145]
and the Natural Science Foundation of Guangxi (Grant Nos. 2018GXNSFGA281007
and 2019JJD110006).

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phy. Rev. Lett., 116, 061102
  * Bacchini et al. (2018a) Bacchini, F., Ripperda, B., Chen, A. Y., Sironi, L. 2018a, Astropys. J. Suppl., 237, 6
  * Bacchini et al. (2018b) Bacchini, F., Ripperda, B., Chen, A. Y.,Sironi, L. 2018b, Astropys. J. Suppl., 240, 40
  * Brown (2006) Brown, J. D. 2006, Phys. Rev. D, 73, 024001
  * Carter (1968) Carter, B. 1968, Phy. Rev., 174, 1559
  * Deng et al. (2020) Deng, C., Wu, X., $\&$ Liang, E. 2020, MNRAS, 496, 2946
  * EHT et al. (2019a) EHT Collaboration, et al. 2019a, ApJL, 875, L1 (Paper I)
  * EHT et al. (2019b) EHT Collaboration et al. 2019b, ApJL, 875, L4 (Paper IV)
  * EHT et al. (2019c) EHT Collaboration et al. 2019c, ApJL, 875, L6 (Paper VI)
  * Felice & Sorge (2003) Felice, D. d, $\&$ Sorge, F. 2003, Class. Quantum Grav., 20, 469
  * Feng (1986) Feng, K. 1986, Journal of Computational Mathematics, 44, 279
  * Hairer et al. (1999) Hairer, E., Lubich, C. $\&$ Wanner, G. 1999, Geometric Numerical Integration, Springer-Verlag, Berlin
  * Hu et al. (2019) Hu, S., Wu, X., Huang, G., $\&$ Liang, E. 2019, ApJ, 887, 191
  * Kerr (1963) Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237
  * Kološ et al. (2015) Kološ, M., Stuchlík, Z., $\&$ Tursunov, A. 2015, Class. Quantum Grav., 32, 165009
  * Kološ et al. (2017) Kološ, M., Tursunov, A., $\&$ Stuchlík, Z. 2017, Eur. Phys. J. C, 77, 860
  * Kopáček et al. (2010) Kopáček, O., Karas, V., Kovář, J., $\&$ Stuchlík, Z. 2010, ApJ, 722, 1240
  * Kopáček & Karas (2014) Kopáček, O., $\&$ Karas, V. 2014, ApJ, 787, 117
  * Li & Wu (2017) Li, D., $\&$ Wu, X. 2017, Mon. Not. R. Astron. Soc., 469, 3031
  * Li & Wu (2019) Li, D., $\&$ Wu, X. 2019, Eur. Phys. J. Plus, 134, 96
  * Liao (1997) Liao, X. H. 1997, Celest. Mech. Dyn. Astron., 66, 243
  * Lichtenberg & Lieberman (1983) Lichtenberg, A. J., $\&$ Lieberman, M. A. 1983, Regular and Chaotic Dynamics (Springer-Verlag, New York)
  * Liu et al. (2016) Liu, L., Wu, X., Huang, G. Q., $\&$ Liu, F. 2016, Mon. Not. R. Astron. Soc., 459, 1968
  * Lubich et al. (2010) Lubich, C., Walther, B., $\&$ Brügmann, B. 2010, Phys. Rev. D, 81, 104025
  * Luo et al. (2017) Luo, J., Wu, X., Huang, G., $\&$ Liu, F. 2017, ApJ, 834, 64
  * Mei et al. (2013a) Mei, L., Ju, M., Wu, X., $\&$ Liu, S. 2013a, Mon. Not. R. Astron. Soc., 435, 2246
  * Mei et al. (2013b) Mei, L., Wu, X., $\&$ Liu, F. 2013b, Eur. Phys. J. C, 73, 2413
  * Pánis et al. (2019) Pánis, R., Kološ, M., Stuchlík, Z. 2019, Eur. Phys. J. C, 79, 479
  * Pihajoki (2015) Pihajoki, P. 2015, Celest. Mech. Dyn. Astron., 121, 211
  * Preto & Saha (2009) Preto, M., $\&$ Saha, P. 2009, ApJ, 703, 1743
  * Reissner (1916) Reissner, H. 1916, Ann. Phys., 50, 106
  * Ruth (1983) Ruth, R. D, 1983, IEEE Trans. Nucl. Sci. NS 30, 2669-2671
  * Seyrich & Lukes-Gerakopoulos (2012) Seyrich, J., $\&$ Lukes-Gerakopoulos, G. 2012, Phys. Rev. D, 86, 124013
  * Schwarzschild (1916) Schwarzschild, K. 1916, Stizber. Deut. Akad. Wiss., Berlin, K1. Math.-Phys. Tech. s., 189
  * Stuchlík & Kološ (2016) Stuchlík, Z., $\&$ Kološ, M. 2016, Eur. Phys. J. C, 76, 32
  * Stuchlík et.al (2020) Stuchlík, Z., Kološ, M., Kovář, J., $\&$ Tursunov, A. 2020, Universe, 6, 26
  * Swope et al. (1982) Swope, W. C., Andersen, H. C., Berens, P. H., $\&$ Wilson, K. R. 1982, J.Chem. Phys. 76, 637
  * Takahashi & Koyama (2009) Takahashi, M., $\&$ Koyama, H. 2009, ApJ, 693, 472
  * Tsang et al. (2015) Tsang, D., Galley, C. R., Stein, L. C., $\&$ Turner, A. 2015, ApJL, 809, L9
  * Tursunov et al. (2016) Tursunov, A. A., $\&$ Kološ, M. 2018, Physics of Atomic Nuclei, 81, 279
  * Tursunov et al. (2020a) Tursunov, A., Stuchlík, Z., Kološ, M., Dadhich, N., $\&$ Ahmedov, B. 2020a, ApJ, 895, 14
  * Tursunov et al. (2020b) Tursunov, A., Zajaček, M., Eckart, A., Britzen, S., $\&$ Stuchlík, Z. 2020b, ApJ, 897, 99
  * Wang et al. (2021) Wang Y., Sun W., Liu F., Wu X. 2021, ApJ (Paper I), 907, 66
  * Wisdom & Holman (1991) Wisdom, J., $\&$ Holman, M. 1991, AJ, 102, 1528
  * Yi & Wu (2020) Yi, M., $\&$ Wu, X. 2020, Phys. Scr., Phys. Scr., 95, 085008
  * Yoshida (1990) Yoshida, H. 1990, Phys. Lett. A, 150, 262
  * Zhong et al. (2010) Zhong, S. Y., Wu, X., Liu, S. Q., $\&$ Deng, X. F. 2010, Phys. Rev. D, 82, 124040

Figure 1: (a) Hamiltonian errors $\Delta K=-1/2-K$ in Equation (13) for the
proposed algorithms solving the system (10). The parameters are $E=0.995$,
$L=4.6$, $\beta=6.4\times 10^{-4}$, $Q=0.1$ and $Q^{*}=1\times 10^{-4}$. A
test orbit has initial conditions $r=25$, $p_{r}=0$ and $\theta=\pi/2$. The
starting value of $p_{\theta}>0$ is determined by Equation (13). The new
second-order and fourth-order explicit symplectic integrators S2 and S4 take
proper time step $h=1$. The realistic errors for S4 are 1000 times smaller
than the plotted errors. The errors for S2 remain bounded and stable in an
order of $\mathcal{O}(10^{-6})$, whereas do not for S4 due to roundoff errors.
However, such a secular drift in Hamiltonian errors is missing when large
proper time step $h=10$ is used in the fourth-order method Y4. These facts
show that the new algorithms with appropriate time steps can share good
properties of a standard symplectic integrator in long-term stabilized error
behavior. (b) Poincaré sections on the plane $\theta=\pi/2$ and
$p_{\theta}>0$, given by algorithm S2 with proper time step $h=1$. The test
orbit in panel (a) is regular Orbit 1 colored red in panel (b).

Figure 2: Same as Figure 1(b), but different values of parameters $Q$ and
$Q^{*}$ are given.

Figure 3: Same as Orbit 1 colored red in Figure 1(b), but different
combinations of parameters $Q$, $Q^{*}$ and $\beta$ are given.


aainstitutetext: Grup de Física Teòrica, Departament de Física, Universitat
Autònoma de Barcelona,
08193 Bellaterra, Spainbbinstitutetext: Institut de Física d’Altes Energies
(IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193
Bellaterra , Spainccinstitutetext: Center for Cosmology and AstroParticle
Physics (CCAPP), Ohio State University, Columbus, OH 43210ddinstitutetext:
Department of Physics, Ohio State University, Columbus, OH
43210eeinstitutetext: Departament de Física Quàntica i Astrofísica and
Institut de Ciències del Cosmos, Universitat de Barcelona, Diagonal 647,
E-08028 Barcelona, Spainffinstitutetext: Institució Catalana de Recerca i
Estudis Avançats (ICREA), Barcelona, Spainffinstitutetext: C.N. Yang Institute
for Theoretical Physics, Stony Brook University, Stony Brook, USA

# On neutrino-mediated potentials in a neutrino background

Diego Blas c,d Ivan Esteban e,f,g M.C. Gonzalez-Garcia e Jordi Salvado
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The exchange of a pair of neutrinos with Standard Model weak interactions
generates a long-range force between fermions. The associated potential is
extremely feeble, $\propto G_{F}^{2}/r^{5}$ for massless neutrinos, which
renders it far from observable even in the most sensitive experiments testing
fifth forces. The presence of a neutrino background has been argued to induce
a correction to the neutrino propagator that enhances the potential by orders
of magnitude. In this brief note, we point out that such modified propagators
are invalid if the background neutrino wavepackets have a finite width. By
reevaluating the 2–$\nu$ exchange potential in the presence of a neutrino
background including finite width effects, we find that the background-induced
enhancement is reduced by several orders of magnitude. This confirms the
intuitive result that coherent enhancements to the potential can only be large
for distances of the order of the background wavepacket width. Unfortunately,
this pushes the resulting 2–$\nu$ exchange potential away from present and
near-future sensitivity of tests of new long-range forces.

##  1 Introduction: short review of neutrino-mediated forces

Soon after the existence of neutrinos was proposed, it was pointed out that
the exchange of neutrino pairs would induce a long-range force between
nucleons Bethe:1936zz ; gamow-teller ; rusos . Most of these early works
estimated the potential to behave as $1/r$, with $r$ the distance between two
test particles, and even suggested that it could be responsible for
gravitation gamow-teller ; gamow ; corben ; bocchieri . However, as early as
1936, Iwanenko and Sokolow rusos  obtained that the potential due to the
exchange of two massless fermions behaves as $1/r^{5}$. This did not prevent
Feynman to consider a $1/r^{3}$ dependence in Ref. Feynman:1963ax ;
Feynman:1996kb , reinforcing the idea that neutrinos could be mediators of a
long-range interaction sourced by mass.

The confusion in the literature concerning the $r$-dependence of the potential
was finally resolved by Feinberg and Sucher Feinberg:1968zz , who carried out
the first calculation of the potential using the effective four-fermion
interaction of weak charged currents of the Standard Model (SM) via dispersion
techniques (the calculation was later complemented to include the neutral
current contribution Feinberg:1989ps ). They found that, for the exchange of
massless neutrino-antineutrino pairs, the interaction potential between two
electrons separated by a distance $r$ is given by

$V(r)=-\frac{G_{F}^{2}}{4\pi^{3}r^{5}},$ (1)

where $G_{F}$ is the Fermi constant, thus confirming the $1/r^{5}$ behaviour.
Hsu and Sikivie Hsu:1992tg  obtained the same result using Feynman
diagrammatic methods.

The effect of the neutrino mass on the potential was first introduced in Ref.
Fischbach:1996qf  for a single massive Dirac neutrino species, and in Ref.
Grifols:1996fk  for the Majorana case; obtaining a potential with a finite
range ${\cal O}(m_{\nu}^{-1})$. The effects from several neutrino masses and
the associated flavor mixing was recently addressed in Refs. Lusignoli:2010gw
; LeThien:2019lxh ; Segarra:2020rah ; Costantino:2020bei . In Ref.
Bolton:2020xsm  the general form of the neutrino forces generated by SM and
beyond the Standard Model (BSM) interactions in the framework of effective
field theories was derived. Finally, Ref. Xu:2021daf  has included
interactions beyond four-fermion contact interactions.

To sum up, the existence of the 2-$\nu$ mediated long-range force is a well-
grounded prediction of the SM, but this force is extremely weak. The
$G_{F}^{2}$ suppression together with the $r^{-5}$ dependence implies that it
is already much weaker than gravity at distances ${\cal O}$(nm). This renders
the effect orders of magnitude below present and near-future sensitivity of
experiments testing the gravitational inverse-square law Hoskins:1985tn ;
Tan:2020vpf ; Lee:2020zjt  and the weak equivalence principle Smith:1999cr ;
Schlamminger:2007ht .

However, if the interaction takes place in the presence of a background of
neutrinos, the 2-$\nu$ mediated potential could be significantly modified.
This was first discussed by Horowitz and Pantaleone Horowitz:1993kw  and later
by Ferrer et. al Ferrer:1998ju ; Ferrer:1999ad , who evaluated the 2-$\nu$
mediated potential in the presence of the cosmic neutrino background
(C$\nu$B). For this case they found a modest enhancement. Most recently, the
effect of neutrino backgrounds has been restudied and extended in Ref.
Ghosh:2022nzo . This work claims that the intense neutrino fluxes from the
Sun, supernovae (SN), or nuclear reactors, may substantially enhance the
potential (by up to 20 orders of magnitude), in particular in the direction of
the incoming background neutrino flux. Such an enhancement could render the
2-$\nu$ potential close to the sensitivity of current and near future fifth
force experiments, which motivated the revision of the effect performed in the
rest of this note.

##  2 Neutrino-mediated forces in a neutrino background: general remarks

The key technical ingredient to evaluate neutrino background effects in Refs.
Horowitz:1993kw ; Ferrer:1998ju ; Ferrer:1999ad ; Ghosh:2022nzo  is the use of
the background-modified propagator obtained in finite temperature field theory
(TQFT) (for reviews of TQFT, see for example Refs. Kapusta:2006pm ;
Laine:2016hma ),

$S_{F}(p)=)(/\\!\\!\\!p+m)\left[\frac{i}{p^{2}-m^{2}+i\epsilon}-2\,\pi\,\delta(p^{2}-m^{2})\left(\theta(p^{0})f_{\nu}(\vec{p})+\theta(-p^{0})f_{\overline{\nu}}(\vec{p})\right)\right]\,,$
(2)

where $f_{\nu}(\vec{p})$ and $f_{\overline{\nu}}(\vec{p})$ are the momentum
distributions of neutrinos and antineutrinos, respectively, normalized so that
$\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}f_{\nu(\overline{\nu})}(\vec{p})$ gives
their number density. In what follows we show that such a modified propagator
does not _always_ correctly quantify the effect of the neutrino background in
the evaluation of the 2-$\nu$-exchange potential. This is true, in particular,
when the finite width of the wavepackets of the neutrinos of the medium are
considered.

Our first observation is that the derivation of the background-modified
fermion propagator in TQFT (Eq. 2) relies on the assumption that the fermions
in the background are in thermal equilibrium, while neither the present
C$\nu$B, nor the neutrino fluxes from the Sun, SN, or nuclear reactors are in
thermal equilibrium close to Earth. Interestingly, one can still use QFT
techniques to show that Eq. (2) may be used for the neutrino propagator in the
presence of _a class of_ backgrounds, see Ghosh:2022nzo . This derivation,
however, implicitly assumes that the neutrino background is well-described by
an incoherent superposition of delocalized wavepackets, i.e., plane waves.
Neither the present C$\nu$B, nor the neutrino fluxes from the Sun, SN, or
nuclear reactors can be well-described as plane waves at _all distance
scales_. Even more, on general grounds, one can show that the Pauli exclusion
principle forbids a superposition of delocalized fermionic wavepackets over
scales smaller than the size of the wavepacket.

We proceed now to evaluating the $2$-$\nu$-exchange in the presence of
realistic neutrino backgrounds, accounting for the finite width of the
wavepackets describing the background state.

##  3 Neutrino-mediated forces in a neutrino background: formalism

The static interaction potential plays a central role in non-relativistic
scattering in quantum mechanics as well as in describing classical forces. To
obtain it given a relativistic QFT, recall that for an interaction between two
fermions $f_{1}$ and $f_{2}$ located at positions $\vec{r}_{1}$ and
$\vec{r}_{2}$, the aforementioned potential reads, Itzykson:1980rh

$V(\vec{r}\equiv\vec{r}_{1}-\vec{r}_{2})=-\int
d^{3}\vec{q}\,e^{i\vec{q}\cdot\vec{r}}\,T^{f_{1}f_{2}}_{\rm NR}(\vec{q})\;,$
(3)

where $T_{\rm NR}^{ff^{\prime}}(\vec{q})$ can be obtained from the fully
relativistic scattering amplitude starting with the ${\cal S}$ matrix for the
transition $\ket{I}\rightarrow\ket{F}$ computed in QFT in momentum space

${\cal
S}_{IF}=\delta_{IF}+i(2\pi)^{4}\delta^{4}(p_{F}-p_{I})T_{IF}^{f_{1}f_{2}}\,,$
(4)

as

$T_{\rm
NR}^{f_{1}f_{2}}(\vec{q})\equiv\frac{T_{IF}^{f_{1}f_{2}}}{4\,m_{f}\,m_{f^{\prime}}}\,,$
(5)

with the four-momentum transfer between the fermions given by
$q\equiv(0,\vec{q})\equiv p^{\prime}_{1}-p_{1}\equiv p_{2}-p^{\prime}_{2}$,
and the helicities of $f_{1}$ and $f_{2}$ not changing in the process.

In the SM the relevant weak interaction Lagrangian leading to the 2–$\nu$
mediated potential can be written in the effective four-fermion interaction
approximation as

$\mathcal{L_{I}}=-\frac{G_{F}}{2\sqrt{2}}\bar{\psi}_{\nu_{i}}\gamma^{\mu}(1-\gamma_{5})\psi_{\nu_{j}}\bar{\psi}_{f}\gamma_{\mu}(g^{f}_{V_{ij}}-g^{f}_{A_{ij}}\gamma_{5})\psi_{f}\,,$
(6)

where $\nu_{i}$ are the neutrino fields with masses $m_{\nu_{i}}$ and $f$
indicates the fermions in the external legs. The effective couplings are

$\displaystyle g_{V_{ij}}^{e}$
$\displaystyle=2U_{ei}U_{ej}-(1-4s_{w}^{2})\delta_{ij}\,,$ (7) $\displaystyle
g_{V_{ij}}^{p}$ $\displaystyle=(1-4s_{w}^{2})\delta_{ij}\,,$ (8)
$\displaystyle g_{V_{ij}}^{n}$ $\displaystyle=-\delta_{ij}\,,$ (9)
$\displaystyle g_{A_{ij}}^{e}$
$\displaystyle=g_{A_{ij}}^{e}=-g_{A_{ij}}^{p}=-\delta_{ij}\,.$ (10)

In what follows, for simplicity we neglect leptonic mixing and assume Dirac
neutrinos. Mixing effects and Majorana neutrinos do not affect our conclusions
and can be easily accounted for.

Figure 1: Feynman diagram contributing to the relevant scattering amplitude in
vacuum

The leading contribution to the 2-$\nu$ exchange potential generated in vacuum
can be obtained from the amplitude represented by the Feynman diagram in Fig.
1. In the above notation, it corresponds to initial and final states
$\ket{I}\equiv\ket{{p}_{1},s_{1}}_{f_{1}}\ket{{p}_{2},s_{2}}_{f_{2}}$ and
$\ket{F}\equiv\ket{{p^{\prime}}_{1},s_{1}}_{f_{1}}\ket{{p^{\prime}}_{2},s_{2}}_{f_{2}}$,
where $p_{i}$ are the initial fermion 4-momenta (with
$p^{0}_{i}=\sqrt{\vec{p}_{i}+m_{i}^{2}}$) and $s_{i}$ their helicities. For
massless neutrinos this leads to the well-known result

$V(r)=-\frac{G_{F}^{2}\,g_{V}^{f_{1}}\,g_{V}^{f_{2}}}{16\,\pi^{3}}\frac{1}{r^{5}}\;.$
(11)

In the presence of a neutrino background new diagrams contribute to the
amplitude at the same order in perturbation theory. More concretely, there is
a process capturing a neutrino from the background with momentum $\vec{k}$,
exchanging a virtual neutrino, and returning another neutrino with momentum
$\vec{k}$ as represented in Fig. 2 (below, we justify that both background
neutrinos have the same momentum). The contributions from an antineutrino
background can be trivially obtained by inverting the neutrino lines.111Notice
that a diagram where both a neutrino and an antineutrino are absorbed or
emitted would modify the fermion energy and hence it does not contribute to
the static potential. This is also true for other scattering processes beyond
those we consider in Fig. 2, and we ignore them.

Figure 2: Feynman diagrams contributing to the relevant scattering amplitudes
in the presence of a neutrino background.

Being of the same order in perturbation theory, these contributions should be
taken into account for computing the potential in any environment where
neutrinos are present. As mentioned above, the neutrino backgrounds that we
shall consider are far from being in thermal equilibrium, and may even entail
momentum distributions far from thermal. This is why the use of TQFT
techniques is not justified _a priori_ , and we will study the influence of
the background explicitly.

The basic assumption in what follows is that the background can be well
described by an incoherent superposition of neutrino wavepackets of the
form222One can add a random phase at each momentum state, but it does not play
any role for the final observable and we ignore it.

$\ket{\psi(\vec{x},\vec{k}^{\prime},s)}=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}\sqrt{2E_{k}}}\omega(\vec{k}^{\prime},\vec{k})e^{-i\vec{k}\cdot\vec{x}}\ket{\vec{k},s}_{\nu}\,,$
(12)

where $\omega(\vec{k}^{\prime},\vec{k})$ are neutrino wavepackets centered at
momentum of $\vec{k}^{\prime}$ and with helicity $s$. The exponential factor
centers the packet at position $\vec{x}$. In what follows, for concreteness in
our quantifications we use Gaussian wavepackets,

$|\omega(\vec{k}^{\prime},\vec{k})|^{2}=\frac{1}{\pi^{3/2}\sigma^{3}}e^{-\frac{|\vec{k}^{\prime}-\vec{k}|^{2}}{\sigma^{2}}}\,,$
(13)

with $\sigma$ the width in momentum space (unlike the position-space width,
that increases with time for free particles, this stays constant when the
particles evolve freely).

To account for the stochastic properties of the medium, we parametrize it in
terms of a density matrix. Assuming that the background is homogeneous (this
can be realized by large enough densities over the probed distances or by
averaging the force over a long enough time),

$\begin{split}\rho^{bkg}&=\sum_{s}\int
d^{3}\vec{k}^{\prime}f_{\nu}(\vec{k}^{\prime})\int
d^{3}\vec{x}\ket{\psi(\vec{x},\vec{k}^{\prime},s)}\bra{\psi(\vec{x},\vec{k}^{\prime},s)}\\\
&=\sum_{s}\int
d^{3}\vec{k}^{\prime}f_{\nu}(\vec{k}^{\prime})\int\frac{d^{3}\vec{k}}{(2\pi)^{3}\sqrt{2E_{k}}}|\omega(\vec{k}^{\prime},\vec{k})|^{2}\ket{\vec{k},s}_{\nu}\bra{\vec{k},s}_{\nu}\end{split}$
(14)

with $f_{\nu}(\vec{k}^{\prime})$ the momentum distribution of the neutrino
ensemble which we assume to be spin-independent. For polarized
ultrarelativistic backgrounds, such as solar or reactor neutrinos, only one of
the helicities contributes, and hence this density matrix is equally
applicable. As we see, when computing averages both background neutrinos will
have the same momentum.

Before proceeding to the calculation of the potential, it is worth pointing
out that the Pauli exclusion principle imposes relevant constraints on the
posible form of a state composed of identical fermions. For a thermal neutrino
background for which $f_{\nu}(\vec{k}^{\prime})$ is a Dirac-Fermi
distribution, Eq. (14) is fully consistent for any width of the wave packet.
On the contrary, for a non-thermal neutrino background, characterized by a
general momentum distribution _with arbitrary normalization_ , the previous
description of the background as an incoherent sum of _single-particle states_
is only physical for interparticle distances larger that the spatial extent of
the wavepackets, $\sim 1/\sigma$. For interparticle distances shorter than the
extent of the wavepacket, the Pauli exclusion principle renders the assumed
momentum distribution inconsistent. That said, in the scenarios we consider
the interparticle distances are always larger than wavepacket extents.

In this approach, we first compute the amplitudes in Fig. 2 and we later trace
out the result with the state that describes the background, Eq. (14). This
corresponds to the following initial and final states in the evaluation of the
amplitude

$\displaystyle\ket{I}$
$\displaystyle=\ket{p_{1},s_{1}}_{f_{1}}\ket{p_{2},s_{2}}_{f_{2}}\ket{k,s}_{\nu}\,,$
(15) $\displaystyle\ket{F}$
$\displaystyle=\ket{p^{\prime}_{1},s_{1}}_{f_{1}}\ket{p^{\prime}_{2},s_{2}}_{f_{2}}\ket{k,s}_{\nu}\,.$
(16)

After some rearrangements, the corresponding ${\cal S}$ matrix element can be
written as,

$\begin{split}{\cal S}_{IF}=-\frac{G_{F}^{2}}{8}\int
d^{4}x_{1}\,d^{4}x_{2}\;{}_{f_{1}}\\!\\!\bra{p^{\prime}_{1},s_{1}}\bar{\psi}_{f}({x}_{1})\gamma_{\mu}(g^{f}_{V}-g^{f}_{A}\gamma_{5})\psi_{f}({x}_{1})\ket{p_{1},s_{1}}_{f_{1}}\times\\\
{}_{f_{2}}\\!\\!\bra{p^{\prime}_{2},s_{2}}\bar{\psi}_{f}({x}_{2})\gamma_{\nu}(g^{f}_{V}-g^{f}_{A}\gamma_{5})\psi_{f}({x}_{2})\ket{p_{2},s_{2}}_{f_{2}}\times\\\
\Bigl{[}{}_{\nu}\\!\\!\bra{{k},s}\bar{\psi}_{\nu}({x}_{1})\gamma^{\mu}(1-\gamma_{5})\,iS_{F}(x_{1}-x_{2})\,\gamma^{\nu}(1-\gamma_{5})\psi_{\nu}({x}_{2})\ket{k,s}_{\nu}+\\\
{}_{\nu}\\!\\!\bra{k,s}\bar{\psi}_{\nu_{i}}({x}_{2})\gamma^{\nu}(1-\gamma_{5})\,iS_{F}(x_{1}-x_{2})\,\gamma^{\mu}(1-\gamma_{5})\psi_{\nu_{j}}({x}_{1})\ket{{k},s}_{\nu}\Bigr{]}\,,\end{split}$
(17)

where $S_{F}(x_{1}-x_{2})$ is the vacuum Feynman neutrino propagator in
position space.

Fourier-expanding the fermion fields in terms of creation and annihilation
operators with well-determined momenta and using the Fourier-space neutrino
propagator, one can easily compute the expectation values. After factorizing
out the global energy-momentum conservation Dirac delta,
$i\,(2\pi)^{4}\delta^{4}(p_{F}-p_{I})=i\,(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-p^{\prime}_{1}-p^{\prime}_{2})$,
and taking the non-relativistic limit we get,

$\begin{split}\sum_{s}T^{f_{1}f_{2}\,\nu}_{NR}=-G_{F}^{2}g^{f_{1}}_{V}g^{f_{2}}_{V}&\Biggl{(}\int
d^{4}x\int\frac{d^{4}q^{\prime}}{(2\pi)^{4}}\frac{[2k^{0}q^{0}-\vec{k}\cdot\vec{q^{\prime}}]}{q^{\prime
2}-m^{2}}e^{-ix(q-q^{\prime}+k)}\\\ &~{}~{}~{}~{}~{}+\int
d^{4}x\int\frac{d^{4}q^{\prime}}{(2\pi)^{4}}\frac{[2k^{0}q^{0}-\vec{k}\cdot\vec{q^{\prime}}]}{{q^{\prime
2}-m^{2}}}e^{-ix(q+q^{\prime}-k)}\Biggr{)}\,,\end{split}$ (18)

where the superindex $\nu$ indicates that the background neutrino state is
still not traced out. Using Eq. 3 and computing the trace with the density
matrix in Eq. 14, we arrive at the following form for the potential between
the two fermions,

$\begin{split}V^{\rm
bkg}_{\sigma}(r)=&\,G_{F}^{2}g^{f_{1}}_{V}g^{f_{2}}_{V}\int
d^{3}\vec{q}\,e^{i\vec{q}\cdot\vec{r}}\int\frac{d^{3}\vec{k}^{\prime}}{(2\pi)^{3}}\,f_{\nu}(\vec{k}^{\prime})\,\int\frac{d^{3}\vec{k}}{(2\pi)^{3}2E_{k}}\,|\omega(\vec{k}^{\prime},\vec{k})|^{2}\\\
&\hskip
113.81102pt\times\Biggl{[}\frac{2|\vec{k}|^{2}+m^{2}+\vec{k}\cdot\vec{q}}{|\vec{k}|^{2}+2\vec{k}\cdot\vec{q}}+\frac{2|\vec{k}|^{2}+m^{2}-\vec{k}\cdot\vec{q}}{|\vec{k}|^{2}-2\vec{k}\cdot\vec{q}}\Biggr{]}\,.\end{split}$
(19)

We find that the same result could have been obtained starting with the
amplitude corresponding to the diagram in Fig. 1 and using an effective
neutrino propagator in the neutrino medium given by

$\displaystyle S_{F}(p)$ $\displaystyle=$
$\displaystyle(/\\!\\!\\!p+m)\left\\{\frac{i}{p^{2}-m^{2}+i\epsilon}-2\,\pi\,\delta(p^{2}-m^{2})\right.$
(20) $\displaystyle\left.\times\left[\theta(p^{0})\int
d^{3}{\vec{p}\,}^{\prime}|\omega_{\nu}(\vec{p},{\vec{p}\,}^{\prime})|^{2}f_{\nu}({\vec{p}\,}^{\prime})+\theta(-p^{0})\int
d^{3}{\vec{p}\,}^{\prime}|\omega_{\nu}(\vec{p},{\vec{p}\,}^{\prime})|^{2}f_{\overline{\nu}}({\vec{p}\,}^{\prime})\right]\right\\}\,\,,$

which the thermal propagator reproduces Eq. 2 in the $\sigma\rightarrow 0$
limit.

In fact, we can formally recover the results in the literature from Eq. 19 by
taking the limit of neutrino wavepackets to be plane waves

$\lim_{\sigma\rightarrow
0}\,|\omega(\vec{k}^{\prime},\vec{k})|^{2}=\delta^{3}(\vec{k}^{\prime}-\vec{k})\;.$
(21)

In particular, if $f(\vec{k}^{\prime})$ is a thermal distribution we recover
the results obtained for the C$\nu$B Horowitz:1993kw ; Ferrer:1998ju ;
Ferrer:1999ad ; Ghosh:2022nzo . Notice however, that as discussed above, the
Pauli exclusion principle makes this limit unphysical for the non-thermal
neutrino backgrounds from the Sun, nuclear reactors, or SN.

##  4 Results and conclusions

From the previous expressions, one can extract the corrections of the width of
wavepackets in the 2-$\nu$ exchange potential. Although the integrals in Eq.
19 cannot be performed analytically for most $f_{\nu}(\vec{p})$ distributions,
the basic lessons can be extracted from the case of a monochromatic
directional flux $\Phi_{0}$

$f_{\nu}(\vec{p})=(2\pi)^{3}\delta^{3}(\vec{p}-\vec{p}_{0})\Phi_{0}\,,$ (22)

with $k_{0}$ parallel to $\vec{r}$, which allows for an analytic treatment.
This example can be used to estimate the modification of the potential due to
solar, SN, and reactor neutrinos. Furthermore, this is the scenario for which
the largest background-induced enhancement is obtained in Ref. Ghosh:2022nzo .
For this case, and neglecting the neutrino mass, we find for
$E_{\nu}=|\vec{k}_{0}|$

$\displaystyle V^{\rm bkg}_{\sigma}(r)$
$\displaystyle=\,-\frac{G_{F}^{2}\,g_{V}^{f_{1}}\,g_{V}^{f_{2}}}{4\
\pi^{3/2}\,r^{2}}\Phi_{0}\frac{1}{(4\hat{E}^{2}+\hat{\sigma}^{4})^{2}}\left\\{-2e^{-\frac{\hat{E}_{\nu}^{2}}{\hat{\sigma}^{2}}}\hat{\sigma}(\hat{\sigma}^{2}-2\hat{E}^{2})(4\hat{E}^{2}+\hat{\sigma}^{4})\right.$
$\displaystyle+\sqrt{\pi}\left[2\hat{E}\left(8\hat{E}^{4}+2(1+\hat{E})^{2}\hat{\sigma}^{4}+\hat{\sigma}^{6}\right){\rm
Erf}\left[\frac{\hat{E}}{\hat{\sigma}}\right]\right.$
$\displaystyle\left.+e^{-\hat{\sigma}^{2}}\left(8\hat{E}^{4}+\hat{\sigma}^{6}+2\hat{\sigma}^{8}+\hat{E}^{2}(6\hat{\sigma}^{4}-4\hat{\sigma}^{2})\right)\left(\cos(2\hat{E}){\rm
Im}(X)+\sin(2\hat{E}){\rm Re}(X)\right)\right.$
$\displaystyle\left.\left.+e^{-\hat{\sigma}^{2}}\left(4\hat{E}\hat{\sigma}^{2}(2\hat{E}^{2}+\hat{\sigma}^{2}+\hat{\sigma}^{4})\right)\left(\cos(2\hat{E}){\rm
Im}(X)-\sin(2\hat{E}){\rm Re}(X)\right)\right]\right\\}\;,$ (23)

where for convenience we have introduced the dimensionless variables
$\hat{E}\equiv E_{\nu}\,r$ and $\hat{\sigma}\equiv\sigma\,r$ (with
$r=|\vec{r}|$), and $X\equiv{\rm
Erf}\left[\frac{\hat{E}^{2}}{\hat{\sigma}^{2}}+i\,\hat{\sigma}\right]$.

The comparison of the potential in Eq. 23 with the expression obtained for
$\sigma\rightarrow 0$

$V^{\rm bkg}_{\sigma\rightarrow
0}(r)=-\frac{G_{F}^{2}\,g_{V}^{f_{1}}\,g_{V}^{f_{2}}}{4\pi\,r}\Phi_{0}\left[E_{\nu}+\frac{\sin(2E_{\nu}r)}{2r}\right]\,,$
(24)

is shown in Fig. 3, where we plot $V^{\rm bkg}_{\sigma}(r)/V^{\rm
bkg}_{\sigma\rightarrow 0}(r)$ as a function of $r$ for energies around
$1\,\mathrm{MeV}$, characteristic of solar, SN, and reactor neutrinos; and
different wavepacket widths. As we see from the figure, including the
wavepacket width suppresses the effect of the background by orders of
magnitude for distances probed by the most sensitive current experiments
(which range from $\mu m$ to $\sim 10^{4}$ km).

Figure 3: Suppression of the background-induced 2-$\nu$-mediated potential
between two fermions due to the inclusion of a finite width packet of the
background states as a function of the distance between the fermions. The
figures illustrate the dependence of the effect on the neutrino energy and the
width of the wavepacket.

For widths $\sigma\ll E_{\nu}$ and distances $r\gg\frac{1}{E_{\nu}}$, Eq. 23
is well approximated by

$V^{\rm bkg}_{\sigma}(r)\simeq\frac{V^{\rm bkg}_{\sigma\rightarrow
0}(r)}{1+r^{2}\sigma^{4}/(4E_{\nu}^{2})}\,.$ (25)

Although in Fig. 3 we use the full expression, this approximation is
excellent. To have an unsuppressed potential due to a background of neutrinos
of energies ${\cal O}$(MeV) over a distance of $\mu m$, the background
neutrino states should have momentum widths $\sigma\ll 10$ keV. For longer
distances, $\sim 10^{4}\,\mathrm{km}$, the width has to be even smaller,
$\sigma\ll\,\mathrm{meV}$.

Estimations on the typical wavepacket sizes for neutrino backgrounds differ in
the literature. For thermal systems with temperature $T$, collisions act as
quantum-mechanical measurements that localize the particles, and $\sigma\sim
T$ Akhmedov:2017mcc . This would lead to $\sigma\sim\mathrm{keV}$ for solar
neutrinos (assuming that neutrinos inherit the wavepacket width of their
parent nuclei), and $\sigma\sim\mathrm{MeV}$ for SN neutrinos Akhmedov:2017mcc
. For reactor neutrinos, estimates are more uncertain, ranging from $\sim
20\,\mathrm{MeV}$ to $\sim 0.1\,\mathrm{keV}$ Akhmedov:2017mcc ;
deGouvea:2021uvg ; Arguelles:2022bvt ; Akhmedov:2022bjs ; Jones:2022hme . In
all cases, the estimated wavepacket widths are far from the required values to
yield enhanced long-ranged effects at large distances.

Our results can be intuitively understood. In quantum mechanics, a potential
induces a _coherent_ phase shift of the wavefunction. For such a shift to be
induced by a background, the latter should be coherent over the considered
distances. It seems reasonable to find a suppression for the effects of the
background for distances longer than the typical wave-packet, corresponding to
the inverse of the width of the background wave-packet in momentum-space,
$\sim 1/\sigma$.

Overall, we find that the long-discussed 2-$\nu$ mediated force remains
undetectable for fifth-force experiments and the discussed neutrino
backgrounds. The effects only seem to be unsuppressed at extremely short
distances (c.f. Fig. 3), making it difficult to envision physical scenarios
able to probe neutrino-induced long-range forces.

###### Acknowledgements.

We thank L. Badurina, O. Eboli, J. Ghiglieri, C. Manuel, N. Sabti, J. Taron
and R. Zukanovich for valuable discussions. We thank Y. Grossman for comments
on an earlier draft. This project has received funding/support from the
European Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie grant agreement No 860881-HIDDeN, as well as from grants
PID2020-115845GBI00 and PID2019-105614GB-C21, “Unit of Excellence Maria de
Maeztu 2020-2023" award to the ICC-UB CEX2019-000918-M, funded by
MCIN/AEI/10.13039/501100011033. MCGG is also supported by USA-NSF grant
PHY-1915093. DB is supported by a ‘Ayuda Beatriz Galindo Senior’ from the
Spanish ‘Ministerio de Universidades’, grant BG20/00228. IFAE is partially
funded by the CERCA program of the Generalitat de Catalunya.

## References

  * (1) H. A. Bethe and R. F. Bacher, _Nuclear Physics A. Stationary States of Nuclei_ , _Rev. Mod. Phys._ 8 (1936) 82.
  * (2) G. Gamow and E. Teller, _Some Generalizations of the $\beta$ Transformation Theory_, _Phys. Rev._ 51 (1937) 289.
  * (3) D. Iwanenko and Z. Sokolow, _Zur Wechselwirkung der schweren Teilchen_ , _Z. Physik_ 102 (1936) 119.
  * (4) G. Gamow, _Probability of Nuclear Meson-Absorption_ , _Phys. Rev._ 71 (1947) 550.
  * (5) H. Corben, _Long Range Nuclear Forces_ , _Nuovo Cimento_ 10 (1953) 1485.
  * (6) P. Bocchieri and P. Gulmanelli, _Gravitational Forces and Quantum Theory_ , _Nuovo Cimento_ 4 (1955) 1016.
  * (7) R. P. Feynman, _Quantum theory of gravitation_ , _Acta Phys. Polon._ 24 (1963) 697.
  * (8) R. P. Feynman, _Feynman lectures on gravitation_. 1996.
  * (9) G. Feinberg and J. Sucher, _Long-Range Forces from Neutrino-Pair Exchange_ , _Phys. Rev._ 166 (1968) 1638.
  * (10) G. Feinberg, J. Sucher and C. K. Au, _The Dispersion Theory of Dispersion Forces_ , _Phys. Rept._ 180 (1989) 83.
  * (11) S. D. H. Hsu and P. Sikivie, _Long range forces from two neutrino exchange revisited_ , _Phys. Rev. D_ 49 (1994) 4951 [hep-ph/9211301].
  * (12) E. Fischbach, _Long range forces and neutrino mass_ , _Annals Phys._ 247 (1996) 213 [hep-ph/9603396].
  * (13) J. A. Grifols, E. Masso and R. Toldra, _Majorana neutrinos and long range forces_ , _Phys. Lett. B_ 389 (1996) 563 [hep-ph/9606377].
  * (14) M. Lusignoli and S. Petrarca, _Remarks on the forces generated by two-neutrino exchange_ , _Eur. Phys. J. C_ 71 (2011) 1568 [1010.3872].
  * (15) Q. Le Thien and D. E. Krause, _Spin-Independent Two-Neutrino Exchange Potential with Mixing and $CP$-Violation_, _Phys. Rev. D_ 99 (2019) 116006 [1901.05345].
  * (16) A. Segarra and J. Bernabéu, _Absolute neutrino mass and the Dirac/Majorana distinction from the weak interaction of aggregate matter_ , _Phys. Rev. D_ 101 (2020) 093004 [2001.05900].
  * (17) A. Costantino and S. Fichet, _The Neutrino Casimir Force_ , _JHEP_ 09 (2020) 122 [2003.11032].
  * (18) P. D. Bolton, F. F. Deppisch and C. Hati, _Probing new physics with long-range neutrino interactions: an effective field theory approach_ , _JHEP_ 07 (2020) 013 [2004.08328].
  * (19) X.-j. Xu and B. Yu, _On the short-range behavior of neutrino forces beyond the Standard Model: from 1/r 5 to 1/r4, 1/r2, and 1/r_, _JHEP_ 02 (2022) 008 [2112.03060].
  * (20) J. K. Hoskins, R. D. Newman, R. Spero and J. Schultz, _Experimental tests of the gravitational inverse square law for mass separations from 2-cm to 105-cm_ , _Phys. Rev. D_ 32 (1985) 3084.
  * (21) W.-H. Tan et al., _Improvement for Testing the Gravitational Inverse-Square Law at the Submillimeter Range_ , _Phys. Rev. Lett._ 124 (2020) 051301.
  * (22) J. G. Lee, E. G. Adelberger, T. S. Cook, S. M. Fleischer and B. R. Heckel, _New Test of the Gravitational $1/r^{2}$ Law at Separations down to 52 $\mu$m_, _Phys. Rev. Lett._ 124 (2020) 101101 [2002.11761].
  * (23) G. L. Smith, C. D. Hoyle, J. H. Gundlach, E. G. Adelberger, B. R. Heckel and H. E. Swanson, _Short range tests of the equivalence principle_ , _Phys. Rev. D_ 61 (2000) 022001.
  * (24) S. Schlamminger, K. Y. Choi, T. A. Wagner, J. H. Gundlach and E. G. Adelberger, _Test of the equivalence principle using a rotating torsion balance_ , _Phys. Rev. Lett._ 100 (2008) 041101 [0712.0607].
  * (25) C. J. Horowitz and J. T. Pantaleone, _Long range forces from the cosmological neutrinos background_ , _Phys. Lett. B_ 319 (1993) 186 [hep-ph/9306222].
  * (26) F. Ferrer, J. A. Grifols and M. Nowakowski, _Long range forces induced by neutrinos at finite temperature_ , _Phys. Lett. B_ 446 (1999) 111 [hep-ph/9806438].
  * (27) F. Ferrer, J. A. Grifols and M. Nowakowski, _Long range neutrino forces in the cosmic relic neutrino background_ , _Phys. Rev. D_ 61 (2000) 057304 [hep-ph/9906463].
  * (28) M. Ghosh, Y. Grossman, W. Tangarife, X.-J. Xu and B. Yu, _Neutrino forces in neutrino backgrounds_ , 2209.07082.
  * (29) J. I. Kapusta and C. Gale, _Finite-temperature field theory: Principles and applications_ , Cambridge Monographs on Mathematical Physics. Cambridge University Press, 2011, 10.1017/CBO9780511535130.
  * (30) M. Laine and A. Vuorinen, _Basics of Thermal Field Theory_ , vol. 925. Springer, 2016, 10.1007/978-3-319-31933-9, [1701.01554].
  * (31) C. Itzykson and J. B. Zuber, _Quantum Field Theory_ , International Series In Pure and Applied Physics. McGraw-Hill, New York, 1980.
  * (32) E. Akhmedov, J. Kopp and M. Lindner, _Collective neutrino oscillations and neutrino wave packets_ , _JCAP_ 09 (2017) 017 [1702.08338].
  * (33) A. de Gouvêa, V. De Romeri and C. A. Ternes, _Combined analysis of neutrino decoherence at reactor experiments_ , _JHEP_ 06 (2021) 042 [2104.05806].
  * (34) C. A. Argüelles, T. Bertólez-Martínez and J. Salvado, _Impact of Wave Packet Separation in Low-Energy Sterile Neutrino Searches_ , 2201.05108.
  * (35) E. Akhmedov and A. Y. Smirnov, _Damping of neutrino oscillations, decoherence and the lengths of neutrino wave packets_ , _JHEP_ 11 (2022) 082 [2208.03736].
  * (36) B. J. P. Jones, E. Marzec and J. Spitz, _The Width of a Beta-decay-induced Antineutrino Wavepacket_ , 2211.00026.