# Impenetrable Barriers in the Phase Space of a Particle Around a Kerr
Rotating Black Hole

Francisco Gonzalez Montoya<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Faculty of Physical Sciences and
Engineering, University of Leeds, Leeds, LS2 9JT, United Kingdom  Instituto de
Ciencias Físicas, Universidad Nacional Autónoma de México, Av. Universidad
s/n, Col. Chamilpa, Cuernavaca, Morelos, 62210, México  Facultad de Ciencias,
Universidad Nacional Autónoma de México, Av. Universidad 3000, Circuito
Exterior s/n, Coyoacán, CP 04510, Ciudad Universitaria, Ciudad de México,
México

###### Abstract

We study the phase space of a particle in the gravitational field of a
rotating black hole described by the Kerr metric from a geometrical
perceptive. In particular, we show the construction of a multidimensional
generalization of the unstable periodic orbits and their stable and unstable
invariant manifolds that direct the dynamics in the phase space. Those stable
and unstable invariant manifolds divide the phase space and are robust under
perturbations. To visualize the multidimensional invariant sets under the flow
in the phase space, we use a method based on the arclength of the trajectories
in phase space.

##  1 Introduction

The rotating black holes are the result of collisions of non-colinear masses
and are present in the centre of galaxies like our Milky Way. The study of the
dynamics of particles and light around rotating black holes has been important
to understanding diverse astronomical phenomena in the universe and detecting
them [1, 2]. Recently observations and new techniques in data data analysis
have made it possible to visualize a black hole using the images of the gas
around it [3].

The Kerr metric describes the gravitational field in a neighbourhood of a
rotating black hole [4]. The motion of a particle around a Kerr black hole is
Liouville integrable [5, 6], its Hamiltonian has 4 independent integrals of
motion: the Hamiltonian $H$, the energy of the particle $E$, the component of
the angular momentum around the rotation axis of the black hole $L$, and $Q$
the Carter integral. Using these integrals is possible to study the system in
detail. Recent analytical studies of the phase space of this system show: the
bifurcations of the main families of periodic orbits [7], homoclinic orbits in
the equatorial plane [8, 9], trajectories converging to photon sphere [11],
the classification of some bounded trajectories using angle-action variables
[10], the classification of radial motion [12], the spherical orbits [13, 14],
and other bounded trajectories relevant for astrophysics [15, 16].

In this work, we take a geometric approach to study and visualize the phase
space of the system. In particular, we consider invariant objects under the
dynamics in the phase space that generalize the unstable hyperbolic periodic
orbits and their stable and unstable invariant manifolds. Considering the hole
phase space, we can construct impenetrable barriers to understand the complete
dynamics.

The rotation of the black hole is characterized by a parameter $a\in[0,1]$
proportional to the angular momentum of the black hole. For $a=0$, the black
hole is not rotating and its gravitational field is described by the
Schwarzschild metric [17]. In this case, the spatial part of the metric is
spherically symmetric, it is natural to construct an effective potential
energy parameterized by the angular momentum $L$ to describe the radial motion
of the particle. Associated with the maximum value of this effective potential
energy, there is a family of unstable hyperbolic periodic orbits. Their stable
and unstable manifolds of unstable hyperbolic periodic orbits play a
fundamental role in the analysis of the phase space of the system. In
particular, for the construction of the impenetrable barriers in the phase
space in the non-integrable perturbed case.

Considering those geometrical properties of the system, we construct a
Normally Hyperbolic Invariant Manifold (NHIM) taking the union of all unstable
periodic orbits associated with the maximum of the effective potential for all
the different directions on the sphere. The stable and unstable manifolds of
this NHIM have coodimension one relative to the constant energy manifold.
Therefore, they are impenetrable barriers that direct the dynamics.

An important property of the NHIMs and their stable and unstable invariant
manifolds is that they are robust under generic perturbation of the
Hamiltonian [18, 19, 20, 21, 22, 23, 24]. Then, the internal NHIM of the
system and the impenetrable barriers in the phase space formed by their stable
and unstable manifolds persist when the symmetries of the system are broken
and the dynamics becomes chaotic.

The article is organized as follows: In section 2, we show the Hamiltonian
function of the system in the Boyer-Lindquist coordinates and the integrals of
motion of the system. The section 3, we construct the NHIM for the non-
rotating black hole, $a=0$, using the effective potential obtained from the
equation for the motion of the radial coordinate. Also, we show plots of the
stable and unstable manifolds of the NHIM using the technique of Lagrangian
descriptors constructed with the arclength of the trajectories in phase space
[25, 26, 27, 28]. In section 4, we break the spatial spherical symmetry with
the inclusion of the rotation of the black hole, $a>0$, and visualize the
changes in the invariant manifolds. Finally, we present the conclusions and
final remarks in the last section.

##  2  Hamiltonian Function and Integrals of Motion

A convenient set of coordinates to express the metric of the systems is the
Boyer-Lindquist coordinates $r,\theta,\phi,t$ [29]. These coordinates simplify
the expression of the Kerr metric and facilitate the analysis of the phase
space of the system [2, 5, 7, 8, 9, 10]. In these coordinates, the interval
has only one mixed term. The Kerr metric with natural units $G,c=1$, in the
Boyer-Lindquist coordinates is given by

$\displaystyle ds^{2}=-m^{2}d\tau^{2}=$ $\displaystyle-$
$\displaystyle\left(1-\frac{2Mr}{\Sigma}\right)dt^{2}-\frac{4Mar\sin^{2}\theta}{\Sigma}dtd\varphi$
$\displaystyle+$ $\displaystyle\frac{\Sigma}{\Delta}dr^{2}+\Sigma
d\theta^{2}+\sin^{2}\theta\left(r^{2}+a^{2}+\frac{2Ma^{2}r\sin^{2}\theta}{\Sigma}\right)d\varphi^{2}\;,$

where $\tau$ is the proper time of the particle, $m$ and $M$, are the mass of
the particle and the black hole respectively, $a$ is the spin parameter and
the functions $\Sigma$ and $\Delta$ are defined as

$\displaystyle\Sigma(r,\theta)$ $\displaystyle=$ $\displaystyle
r^{2}+a^{2}\cos^{2}\theta\;,$ $\displaystyle\Delta(r)$ $\displaystyle=$
$\displaystyle r^{2}-2Mr+a^{2}\;.$ (2)

In order to simplify the notation, we use the convention $M=1$. From his
metric is possible to calculate the Hamiltonian function of a particle around
a rotating black hole. In particular, we are going to use the Hamiltonian
function constructed in the reference [10]. The Hamiltonian function of the
systems in Boyer-Lindquist coordinates is given by

$H(r,\theta,\varphi,t,p_{r},p_{\theta},p_{\varphi},p_{t})=\frac{\Delta}{2\Sigma}p_{r}^{2}+\frac{1}{2\Sigma}p_{\theta}^{2}-\frac{R+\Delta\Theta}{2\Delta\Sigma}-\frac{1}{2}\;,$
(3)

where the functions $\Theta$, $R$ are defined as

$\displaystyle\Theta(\theta)$ $\displaystyle=$ $\displaystyle
Q-\cos^{2}\theta\;\left\\{a^{2}(1-E^{2})+\frac{L^{2}}{\sin^{2}\theta}\right\\}\;,$
$\displaystyle R(r)$ $\displaystyle=$ $\displaystyle
P^{2}-\Delta\left\\{r^{2}+(L-aE)^{2}+Q\right\\}\;,$ $\displaystyle P(r)$
$\displaystyle=$ $\displaystyle E(r^{2}+a^{2})-aL\;.\quad$ (4)

The corresponding Hamilton’s equations of motion parametrized by the proper
time $\tau$ of the particle are

$\displaystyle\dot{r}$ $\displaystyle=$
$\displaystyle\frac{\Delta}{\Sigma}p_{r},$ $\displaystyle\dot{p}_{r}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
r}\left(\frac{\Delta}{2\Sigma}\right)p^{2}_{r}-\frac{\partial}{\partial
r}\left(\frac{1}{2\Sigma}\right)p^{2}_{\theta}+\frac{\partial}{\partial
r}\left(\frac{R+\Delta\Theta}{2\Delta\Sigma}\right),$
$\displaystyle\dot{\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{\Sigma}p_{\theta},$ $\displaystyle\dot{p}_{\theta}$
$\displaystyle=$
$\displaystyle\frac{\partial}{\partial\theta}\left(\frac{\Delta}{2\Sigma}\right)p^{2}_{r}-\frac{\partial}{\partial\theta}\left(\frac{1}{2\Sigma}\right)p^{2}_{\theta}+\frac{\partial}{\partial\theta}\left(\frac{R+\Delta\Theta}{2\Delta\Sigma}\right),$
$\displaystyle\dot{\phi}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\Delta\Sigma}\frac{\partial}{\partial
E}\left(R+\Delta\Theta\right),$ $\displaystyle\dot{p}_{\phi}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\dot{t}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\Delta\Sigma}\frac{\partial}{\partial
L}\left(R+\Delta\Theta\right),$ $\displaystyle\dot{p}_{t}$ $\displaystyle=$
$\displaystyle 0.$ (5)

This Hamiltonian system is integrable, the system of differential equations
has 4 independent integrals of motion: the value of the Hamiltonian $H$, the
energy $E=-p_{t}$ of the particle, the angular momentum $L=p_{\varphi}$, and
the Carter integral $Q$. The first 3 integrals are immediate to obtain from
the structure of the Hamiltonian, however the last one, $Q$ was obtained by
Carter in [5, 15] proposing a separation of variables of solution of the
resulting Hamilton–Jacobi equation in other canonical coordinates. From that
approach is possible to obtain the first-order geodesic equations in a natural
way. Those differential equations in the Boyer-Lindquist coordinates are given
by

$\displaystyle m\Sigma\dot{r}$ $\displaystyle=$ $\displaystyle\pm\sqrt{R}\;,$
$\displaystyle m\Sigma\dot{\theta}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\Theta}\;,$ $\displaystyle m\Sigma\dot{\varphi}$
$\displaystyle=$ $\displaystyle\frac{a}{\Delta}(2rE-
aL)+\frac{L}{\sin{\theta}}\;,$ $\displaystyle m\Sigma\dot{t}$ $\displaystyle=$
$\displaystyle\frac{(r^{2}+a^{2})^{2}E-2arL}{\Delta}-a^{2}E\sin{\theta}\;.$
(6)

This set of differential equations is very convenient for analysing the
analytical important properties like bifurcations of the integrable system [7,
13, 15, 16]. However, for the numerical calculations in the present article,
we are going to use Hamilton’s equations 5 like in works [8, 9, 10]. We use
the equations of motion obtained directly from the Hamiltonian function 3
because we want to understand properties relevant to non-integrable cases,
where all the equations of motion are coupled. Also from the numerical point
of view the system 5 are easier to code because there are no change signs
related to the square roots like in 6.

Also, we can construct the reduced 6-dimensional phase space to describe the
dynamics of the system using the global time $t$ to parametrize the solutions
of the systems of differential equations. To construct the system of
differential equations for the reduced phase space, we only need to use the
differential equation for $\dot{t}=\frac{dt}{d\tau}$ and the chain rule
$\frac{d}{dt}=\frac{d\tau}{dt}\frac{d}{d\tau}$ to calculate the derivative
whit respect to $t$ of the other coordinates and momenta in the system, for
more details see the reference [9].

##  3 NHIM and its Stable and Unstable Manifolds for the Schwarzschild Black
Hole, $a=0$

The symmetries of the system facilitate the understanding of its
multidimensional phase space. In the present case, the symmetries allow us to
see the basic elements that constitute the invariant manifolds under the
dynamics. Some recent examples where is possible to find the Normally
Hyperbolic Invariant Manifold (NHIM) and their invariant manifolds using this
approach are in the works [22, 23, 32, 33, 26, 30, 31].

First, let us consider the simplest case, $a=0$, where the Kerr metric is
reduced to the Schwarzschild metric. In this case, the system has spatial
spherical symmetry it is easy to construct and visualize the generalization of
the hyperbolic periodic orbit and their invariant stable and unstable
manifolds. Without loss of generality, we set $Q=0$ that corresponds to the
motion on the equatorial plane, $\theta=\pi/2$. It is equivalent to any other
plane due to the spherical symmetry of the Schwarzschild metric. From the
geodesic equation for the radial coordinate 6, we can obtain the conservation
equation for the energy $E$.

$\frac{\dot{r}^{2}}{2}+V_{eff}(r;L)=E_{eff}=\frac{E^{2}}{2}\;,$ (7)

where

$V_{eff}(r;L)=\frac{m^{2}}{2}-\frac{m}{r}+\frac{L^{2}}{2r^{2}}-\frac{L^{2}}{r^{3}}\;$
(8)

is the effective potential parametrized by the conserved component of the
angular momentum $L$ and $E$ is the total energy of the particle.

Figure 1: Effective potential $V_{eff}(r;L)$ for different values of $L$. As
the value of $L$ is increased the maximum value of $V_{eff}(r;L)$ grows.
Associated with each maximum there is an unstable hyperbolic fixed point
$X_{L}$ in the phase space of the one-dimensional system with energy
$E=E_{h}$. This fixed point is located at $X_{L}=(r=r_{L},p_{r}=0)$ where
$r_{L}$ is the critical value of $V_{eff}(r;L)$.

Due to the existence of this effective potential $V_{e}ff(r;L)$ for the radial
motion, let us consider the reduced phase space only for the radial degree of
freedom. Associated with the maximum of the effective potential $V_{eff}(r;L)$
there is an unstable hyperbolic fixed point $X_{L}=(r=r_{L},p_{r}=0)$. There
are two special sets in the phase space that divide the reduced 2-dimensional
phase space in regions with different behaviour, the stable and unstable
manifolds of the unstable hyperbolic fixed point denoted by $W^{s/u}(X_{h})$.
The stable (unstable) manifold of the hyperbolic fixed point is the trajectory
that converges (diverges) when the time is increased.

When we include the angular coordinate $\phi$ in the analysis of the dynamics,
we have a 2-degree-of-freedom system. The unstable hyperbolic fixed point
$X_{L}$ corresponds to a circular unstable hyperbolic periodic orbit
$\gamma_{L}$ with radio $r_{L}$ in the configuration space. The unstable
hyperbolic periodic orbit $\gamma_{L}$ also have stable and unstable manifolds
$W^{s/u}(\gamma_{L})$. Those 2-dimensional invariant manifolds are formed by
all the trajectories that converge (diverge) to $\gamma_{L}$ when the time to
infinity. Also by construction, $W^{s/u}(\gamma_{L})$ divides the constant
energy manifold of the 2-degree-of-freedom system.

Considering the spatial spherical symmetry of the system, we construct the
generalization of this unstable hyperbolic periodic orbit in the constant
energy manifold, the Normally Hyperbolic Invariant Manifold (NHIM) denoted as
$\mathcal{M}_{E}$. Let us take the union of all the unstable circular periodic
orbits with the energy $E$, associated with the maximum of $V_{eff}(r;L)$, for
all the possible angular directions $\theta,\phi$ on the sphere in the
configuration space.

$\displaystyle\mathcal{M}_{E}$ $\displaystyle=$
$\displaystyle\bigcup_{\theta,\phi}\gamma_{L}=S^{2}\times\gamma_{L}$ (9)

By construction, the topology of the NHIM $\mathcal{M}_{E}$ is a 3-dimensional
sphere $S^{3}$. If we consider the degree of freedom corresponding to time $t$
in the construction of the NHIM, the result is a 4-dimensional unbounded NHIM
given by $\mathcal{M}_{E}\times\mathbb{R}$.

Analogously, we construct the stable and unstable manifolds of the NHIM
$\mathcal{M}_{E}$. Again, considering the spherical spacial symmetry of the
system and taking the union of all the stable or unstable manifolds of
$W^{s/u}(\gamma_{L})$ we obtain the invariant manifolds of the NHIM
$\mathcal{M}_{E}$,

$\displaystyle W^{s/u}(\mathcal{M}_{E})$ $\displaystyle=$
$\displaystyle\bigcup_{\theta,\phi}W^{s/u}(\gamma_{L}).$ (10)

By construction, the topology of $W^{s/u}(\mathcal{M}_{E})$ is a 4-dimensional
spherical cylinder $S^{3}\times\mathbb{R}$. This spherical cylinder divides
the 5-dimensional constant energy manifold in the reduced phase space that
considers only the spatial degrees of freedom. Also, if we consider the degree
of freedom corresponding to time $t$ in the construction of the invariant
stable and unstable manifolds of the NHIM, the result is a 5-dimensional
invariant manifold given by $W^{s/u}(\mathcal{M}_{E})\times\mathbb{R}$.

In order to visualize the invariant manifolds in multidimensional phase space,
we calculate the scalar field of arclengths of the trajectories. The arclength
of the trajectories is generated by their different behaviour give us
information on the invariant manifolds that intersect the set of initial
conditions where we calculate the scalar field. This scalar field is a type of
Lagrangian descriptor scalar field. Some recent examples where the Lagrangian
descriptors have been used to analyze the phase space of multidimensional
systems are in [25, 26, 34, 35]. The basic definition of the Lagrangian
descriptor based on arclength is in Appendix A and more illustrative examples
are in the reference therein.

Let us consider some relevant values of $L$ where the phase space changes
qualitatively. In this case, the changes in the phase space are related to
changes in the geometry of the effective potential energy $V_{eff}(r;L)$. In
Fig.4, we can see two important changes in the geometry of the effective
potential energy $V_{eff}(r;L)$ as the value of $L$ is increased. The first
important change in the effective potential $V_{eff}(r;L)$ is the collision
between its maximum with its minimum at $L=L_{sc}=\sqrt{12}$. This collision
generates a saddle-centre bifurcation in the phase space. To visualize the
changes in the phase space, we calculate the Lagrangian descriptor for initial
conditions in the canonical plane $r$–$p_{r}$ and two constant values of $L$,
one slightly smaller than the critical value and the other slightly bigger,
see Fig. 2 panels a) and b). We can see how the phase space grows abruptly and
a bounded region appears. In Fig. 2 b), the boundaries of phase space are the
intersection of invariant manifolds $W^{s/u}(\gamma_{L})$ with the set of
initial conditions and the intersection point is the intersection with the
unstable hyperbolic periodic orbit $\gamma_{L}$ .

The second important change in $V_{eff}(r;L)$ occurs when its maximum reach
the asymptotic value $1/2$. For that value of $L=4$, the trajectories can
escape to infinity and the phase space has a transition from bounded to
unbounded. Analogously, we calculate the Lagrangian descriptor for initial
conditions in the canonical plane $r$–$p_{r}$ and two constant values of $L$,
one value slightly smaller than the critical value and other chosen value
slightly bigger see Fig.2 panels c) and d). In Fig. 2 d), we can see the
external branches of invariant manifolds $W^{s/u}(\gamma_{L})$ go to infinity.

|
---
a) $L=3.4641<L_{sc}$
|
---
b) $L=3.6>L_{sc}$ , $E=E_{h}$
|
---
c) $L=3.65$, $E=E_{h}$
|
---
d) $L=4.0$, $E=E_{h}$
Figure 2:  Lagrangian descriptor evaluated on the plane $r$–$p_{r}$ for
different values of $E$, $L$, $p_{\theta}=0$, and $a=0$. The energy $E=E_{h}$
corresponds to the energy of the hyperbolic fixed point associated with the
maximum of the effective potential $V_{e}ff(r;L)$.

Now, we consider the Lagrangian descriptor evaluated on initial conditions in
the plane $r$–$p_{r}$ with the value of $E$ constant and the initial component
of the momentum $p_{\theta}=0$. Let us remark that for these calculations the
value of $L$ is determined by the other initial conditions and changes from
point to point in the next plots. The Fig. 3 shows the results for 4 different
values of $E$ or equivalently $E_{eff}$.

For $E_{eff}=0.44$, see Fig.3 a), the external branch of the stable manifold
$W^{s}(\mathcal{M}_{E})$ and the external branch of unstable manifolds
$W^{s/u}(\mathcal{M}_{E})$ coincide and divide the regions of the phase space
with different behaviour. We can differentiate 3 regions: the yellow region
corresponds to trajectories that are bounded, the dark blue region corresponds
to trajectories that cross the event horizon and go inside the black hole, and
the green region corresponds to trajectories that are unbounded and escape to
infinity.

For $E_{eff}=0.45$, see Fig.3 b), we have a similar phase space structure with
also 3 regions. However, the bounded yellow region does not reach the boundary
of the domain of the Lagrangian descriptor plot. The intersection of the
invariant manifolds of the NHIM $\mathcal{M}_{E}$ with the set of initial
conditions form a closed curve that defines the boundary of the trapped
region.

For $E_{eff}=0.49$, see Fig.3 c), shows that the trapped regions inside the
separatrix formed by $W^{s/u}(\mathcal{M}_{E})$ disappear and the region that
goes to the event horizon grows. The yellow-green regions correspond to
regions with a large value of final $r$.

For $E_{eff}=0.51$, see Fig.3 d), the domain of the Lagrangian descriptor plot
is unbounded, but the phase space with $r\in[0,30]$ has a similar structure
like the panel c).

|
---
a) $E_{eff}=0.44$
|
---
b) $E_{eff}=0.45$
|
---
c) $E_{eff}=0.49$
|
---
d) $E_{eff}=0.51$
Figure 3:  Lagrangian Descriptor evaluated on the plane $r$–$p_{r}$ for
different values of $E$, and $a=0$.

##  4 NHIM and its Invariant Stable and Unstable Manifolds for Kerr Black
Hole, $a>0$

For the rotating black hole, its angular momentum is non-zero, i.e $a>0$, the
system lost its spatial spherical symmetry, however, it is still an integrable
one [5]. A natural question is: What happens to the NHIM $\mathcal{M}_{E}$ and
their stable and unstable manifolds $W^{s/u}(\mathcal{M}_{E})$? In general,
the NHIMs and their stable and unstable manifolds are robust under
perturbations of the vector field that define the system of ordinary
differential equations [18, 21]. This fundamental geometrical property has
been essential to understanding the dynamics in the multidimensional phase
space for different systems where symmetry is broken, for example [32, 33, 26,
36, 37].

Let us denote the NHIM for $a>0$ as $\mathcal{M}_{E,a}$ and its stable and
unstable invariant manifolds as $W^{s/u}(\mathcal{M}_{E,a})$. Recent works
based on the integrability of the system and a change of variables [15, 13,
38] show analytical expressions for the trajectories contained in the NHIM
$\mathcal{M}_{E,a}$ and their stable and unstable manifolds
$W^{s/u}(\mathcal{M}_{E,a})$ for $a\in[0,1]$.

Moreover, we can understand the geometrical origin of the existence of the
NHIM and its stable and unstable manifolds analogously than in the non-
rotating case, $a=0$, considering again the equation 6 for the radial motion
and the dynamical systems theory. In this case, the equation for the radial
motion can be written as

$\frac{m^{2}\Sigma^{2}}{r^{4}}\frac{\dot{r}^{2}}{2}+V_{eff}(r;a,E,L,Q)=E_{eff}=\frac{E^{2}}{2}\;,$
(11)

where

$V_{eff}(r;a,E,L,Q)=\frac{m^{2}}{2}-\frac{m}{r}+\frac{L^{2}-a^{2}(E^{2}-m^{2})+Q}{2r^{2}}-\frac{(L-aE)^{2}+Q}{r^{3}}+\frac{a^{2}Q}{2r^{4}}\;.$
(12)

On the RHS of equation 11 we have the radial kinetic energy and the
generalized effective potential energy $V_{eff}(r;a,E,L,Q)$. On the LHS of
equation 11 we have a constant interpreted as the effective energy.

Let us notice that for $a=0$ and $Q=0$, the effective potential 12 and the
radial equation of motion 11 are reduced to corresponding equations 8 and 7 in
the previous section. Also, we obtain the same functional form for the
effective potential with the change $L^{2}+Q\to L^{2}$ for the case $a=0$ and
$Q\neq 0$. That represents the motion outside the equatorial plane in the non-
rotational case. For more information about the bifurcation diagrams for all
the parameters see the references [7, 13].

Figure 4: Effective potential $V_{eff}(r;a,E,L,Q)$ for different values of
$a,E,L,Q$. Associated with each maximum there is an unstable hyperbolic fixed
point in the phase space of the 1-dimensional system. This fixed point is
located at $(r,p_{r})=(r_{(a,E,L,Q)},0)$ where $r_{(a,E,L,Q)}$ is the critical
value of $V_{eff}(r;a,E,L,Q)$.

Let us remark that we can apply the theorem of persistence of NHIMs and their
invariant manifolds to $\mathcal{M}_{E,a}$ and $W^{s/u}(\mathcal{M}_{E,a})$
for a perturbed Hamiltonian originated by a perturbed metric. Typically, the
system becomes non-integrable and its dynamics is chaotic. Then the stable and
stable manifolds $W^{s/u}(\mathcal{M}_{E,a})$ intersect transversally an
infinite number of times and form a tangle in the phase space. This tangle is
the multidimensional generalisation of the Smale horseshoes for the Poincare
map of the 2-degree-of-freedom system, some examples are in the references
[32, 36, 39, 35].

The next figures show changes in the phase space for different values of the
effective energy of the particle $E_{eff}$ and different values of the
rotation parameter $a$. In general, we can appreciate how the trapped region
grows when the rotation of the black hole is increased.

For the effective energy $E_{eff}=0.44$, Fig. 5, the domain of the Lagrangian
descriptor plot is bounded and the separatrix formed by the invariant
manifolds $W^{s/u}(\mathcal{M}_{E,a})$. For $E_{eff}=0.45$, Fig. 6, the
trapped region is bigger than in the previous case but the phase space
structure is very similar. In Fig. 7, $E_{eff}=0.49$, the structure of the
phase space changes, and the trapped region does not reach the maximum value
of $r$ in the domain for moderate values of $a$. For $E_{eff}=0.51$, the phase
space becomes unbounded, see Fig. 8. However, the trapped region has a similar
behaviour to the previous case.

|
---
$a=0.1$
|
---
$a=0.2$
|
---
$a=0.3$
|
---
$a=1.0$
Figure 5:  Lagrangian Descriptor evaluated on the plane $r$–$p_{r}$ for $E_{eff}=0.44$, different values of $a$.  |
---
$a=0.1$
|
---
$a=0.2$
|
---
$a=0.3$
|
---
$a=1.0$
Figure 6:  Lagrangian Descriptor evaluated on the plane $r$–$p_{r}$ for $E_{eff}=0.45$, different values of $a$.  |
---
$a=0.1$
|
---
$a=0.2$
|
---
$a=0.3$
|
---
$a=1.0$
Figure 7:  Lagrangian Descriptor evaluated on the plane $r$–$p_{r}$ for $E_{eff}=0.49$, different values of $a$.  |
---
$a=0.1$
|
---
$a=0.2$
|
---
$a=0.3$
|
---
$a=1.0$
Figure 8:  Lagrangian Descriptor evaluated on the plane $r$–$p_{r}$ for
$E_{eff}=0.51$, different values of $a$.

##  5 Conclusions and Remarks

In this work, we study the phase space of a particle around a rotating Kerr
black hole using the Boyer-Lindquist coordinates. The dynamics of this system
is integrable and has 4 conserved quantities: the value of the Hamiltonian
$H$, the energy $E=-p_{t}$ of the particle, the angular momentum
$L=p_{\varphi}$, and the Carter integral $Q$. From the radial first integral
equation is possible to construct an effective potential $V_{eff}(r;a,E,L,Q)$
parametrized by the rotation parameter of the black hole $a$ and the conserved
quantities of the system. Using this radial effective potential, we construct
a family of unstable periodic orbits associated with the maximum value of the
effective potential for a fixed value of $E$. The union of those unstable
hyperbolic periodic orbits form a Normally Hyperbolic Invariant Manifold
(NHIM) $\mathcal{M}_{E,a}$ that has a central role in the dynamics in the
reduced phase space.

The stable and unstable manifolds $W^{s/u}(\mathcal{M}_{E,a})$ have
codimension one relative to the constant energy manifold and divide it into
regions with different behaviour, see Figs.3,5,6,7,8. Due to the integrability
of the system, one branch of the stable and unstable manifolds
$W^{s/u}(\mathcal{M}_{E,a})$ coincide and form a separatrix that encircles the
KAM tori associated with the minimal of the effective potential. The internal
branches of $W^{s/u}(\mathcal{M}_{E,a})$ form an impenetrable barrier in the
phase space that defines the set of trajectories that can cross the event
horizon. In one way, the phase space of this system is similar to other
multidimensional scattering systems [26, 32, 36, 39], but in this case the
open branches of $W^{s/u}(\mathcal{M}_{E,a})$ go to the central region instead
of the asymptotic region, $r\to\infty$.

We analyze the changes in the phase space of the system when the parameters
change using Lagrangian descriptor plots. On those plots, we can see the
intersections of the $W^{s/u}(\mathcal{M}_{E,a})$ with the set of initial
conditions. In particular, we appreciate that when the rotation of the black
hole grows the trapped region also grows. Therefore, when the angular moment
of the black hole grows more it drags more of the particles around it.

An important difference between the cases $a>0$ and $a=0$ is that the NHIM
constructed from the prograde orbits, $L>0$, and the NHIM constructed from the
retrograde orbits, $L<0$, have different spatial projections because the
effective potential energy $V_{e}ff(r;a,E,L,Q)$ is not symmetric respect the
change $L\to-L$.

All the results for bounded NHIMs in the reduced phase space defined by the
spatial degrees of freedom $r,\theta,\phi$ can be generalized to the case when
we include the time $t$ as a degree of freedom. In this case, we need to
consider the theory for unbounded NHIMs and their invariant stable and
unstable manifolds developed in [24].

In future works, we are going to study the phase space of a perturbed Kerr
Black hole. When the metric is perturbed, the symmetries are broken and the
dynamics of the system becomes chaotic. However, the NHIMs and their stable
and unstable invariant manifolds are robust under perturbations. This
fundamental property allows us to analyse the perturbed chaotic systems. For a
non-integrable perturbation, the stable and unstable manifolds of NHIM
intersect transversely and form a rich structure analogous to the Smale
horseshoes in more dimensions. The stable and unstable manifolds form a system
of tubes that determine the transport in the phase space in the region close
to the NHIM. Also, the internal dynamics of the NHIM becomes chaotic and give
us information about the bifurcations of the NHIM.

Also, the NHIMs $\mathcal{M}_{E,a}$ define a dividing surface in the phase
space like in Wiggner’s Transition State Theory for chemical reaction dynamics
in phase space [40, 41, 42, 43]. Using this surface in the phase space is
possible to calculate the total flux of trajectories that cross the saddle in
potential energy and can reach the external event horizon.

##  6 Acknowledgments

This research was funded by DGAPA UNAM grant number AG–101122,
CONAHCyT CF–2023–G–763, and CONAHCyT fronteras grant number 425854. The author
thanks the Virtual Institute of Physics at New York (VIP–NY), and Centro
Internacional de Ciencias AC–UNAM for their facilities during the early stages
of this work.

##  7 Appendix A: Lagrangian descriptors a tool to visualize the phase space

The Lagrangian descriptors are scalar fields to visualize invariant objects in
the multidimensional phase space. The main idea behind the method is based on
the different behaviour of the trajectories contained in different invariant
objects in the phase space [25, 26, 27, 28]. Some recent illustrative examples
of this phase space visualisation technique are in the references [44, 45,
46]. For instance, let us consider two trajectories with close initial
conditions. One trajectory is contained in a stable manifold of an unstable
periodic orbit $W^{s}(\gamma)$ and the other trajectory, with close initial
conditions, but not contained in $W^{s}(\gamma)$. The trajectory that belongs
to $W^{s}(\gamma)$ converges to $\gamma$, meanwhile the other trajectory is
close only for a finite interval of time. This different behaviour generates a
significant difference in the value of the scalar field.

A way to measure this different behaviour of the trajectories is the arclength
of the trajectories for a given evolution time. The Lagrangian descriptor
evaluated on the point $x_{0}$ in the phase space is defined as

$LD[x(t)]=\int^{t_{f}}_{t_{i}}\sqrt{\sum_{j}(\dot{x}_{j}(t))^{2}}dt\;,$ (13)

where $x(t)$ is the solution of the ODE system such that $x(t_{0})=x_{0}$,
$\dot{x}_{j}$ are the components of the vector field that define the ODE
system $\dot{x}=f(x)$, and $[t_{i},t_{f}]$ is the integration interval. In
this example, if the particle reaches the exterior event horizon, we stop the
integration.

##  8 Appendix B: Stability of Periodic Orbits on the NHIM and its relation
with the Effective Potential

The NHIM $\mathcal{M}_{E}$ is constructed with the periodic orbits
$\gamma_{L}$. Those periodic orbits and their stability on the radial
direction are defined by the equations

$\displaystyle\dot{r}$ $\displaystyle=$
$\displaystyle\pm\frac{\sqrt{R}}{m\Sigma}=0\;,$ (14)
$\displaystyle\frac{\partial\dot{r}}{\partial r}$ $\displaystyle>$
$\displaystyle 0\;.$ (15)

Substituting the definitions of $R$ and $\Sigma$, we can find that previous
equations are equivalent to

$\displaystyle R$ $\displaystyle=$ $\displaystyle 0\;,$ (16)
$\displaystyle\frac{dR}{dr}$ $\displaystyle=$ $\displaystyle 0\;,$ (17)
$\displaystyle\frac{d^{2}R}{dr^{2}}$ $\displaystyle>$ $\displaystyle 0\;.$
(18)

In terms of the effective potential energy, $V_{eff}$ those conditions are
related to its local maximum given by

$\displaystyle\frac{\partial V_{eff}}{\partial r}$ $\displaystyle=$
$\displaystyle 0\;,$ (19) $\displaystyle\frac{\partial^{2}V_{eff}}{\partial
r^{2}}$ $\displaystyle<$ $\displaystyle 0\;.$ (20)

For more information about the parameter space and the bifurcation diagrams of
the system, see the recent works [7, 13, 15, 16].

## References

  * [1] S. Chandrasekhar, The Mathematical Theory of Black Holes (Oxford: Clarendon Press, 1983)
  * [2] D. L. Wiltshire, M. Visser, S. M. Scott, The Kerr Spacetime: Rotating Black Holes in General Relativity (Cambridge University Press, 2009)
  * [3] The Event Horizon Telescope Collaboration et al, First M87 Event Horizon Telescope Results. IV, Imaging the Central Supermassive Black Hole, The Astrophysical Journal Letters. 875, 1, L1 (2019)
  * [4] R. P. Kerr, Gravitational field of a spinning mass as an example of algebraically special metrics. Phys. Rev. Lett. 11, 237 (1963)
  * [5] B. Carter, Global Structure of the Kerr Family of Gravitational Fields, Phys. Rev. 174, 1559 (1968)
  * [6] B. Carter, Axisymmetric Black Hole Has Only Two Degrees of Freedom, Phys. Rev. 26, 331 (1971)
  * [7] I. Bizyaev and I. Mamaev, Bifurcation diagram and a qualitative analysis of particle motion in a Kerr metric, Phys. Rev. D 105, 063003 (2022)
  * [8] J. Levin and G. Perez-Giz, Homoclinic Orbits around Spinning Black Holes I: Exact Solution for the Kerr Separatrix. Phys. Rev. D 79, 124013 (2009)
  * [9] J. Levin and G. Perez-Giz, Homoclinic Orbits around Spinning Black Holes II: The phase space portrait, Phys. Rev. D 79, 124014 (2009)
  * [10] J. Levin and G. Perez-Giz, A periodic table for hole orbits, Phys. Rev. D 77, 103005 (2008)
  * [11] A. Sneppen, Divergent reflections around the photon sphere of a black hole, Sci. Rep. 11, 14247 (2021).
  * [12] G. Compère, Y Liu, J Long, Classification of radial Kerr geodesic motion, Phys. Rev. D 105, 024075 (2022)
  * [13] E. Teo, Spherical orbits around a Kerr black hole. General Relativity and Gravitation, 53, 10 (2022)
  * [14] E. Teo, Spherical photon orbits around a Kerr black hole. General Relativity and Gravitation 35, 1909 (2003)
  * [15] P. Rana and A. Mangalam, Astrophysically relevant bound trajectories around a Kerr black hole. Classical and Quantum Gravity, 36, 045009 (2019)
  * [16] P. Rana and A. Mangalam, Bound orbit domains in the phase space of the Kerr geometry, The Fifteenth Marcel Grossmann Meeting, 858 (2022)
  * [17] R. Wald, General Relativity (Chicago, 1984)
  * [18] N. Fenichel, Persistence and smoothness of invariant manifolds for flows, Indiana Univ. Math. J., 21, 193 (1971)
  * [19] N. Fenichel, Asymptotic Stability With Rate Conditions. Indiana Univ. Math. J. 23, 12, 1109 (1974)
  * [20] N. Fenichel, Asymptotic Stability with Rate Conditions II, Indiana Univ. Math. J. 26, 1, 81 (1977)
  * [21] S. Wiggins, Normally Hyperbolic Invariant Manifolds in Dynamical Systems  (Springer Verlag, Berlin 1994)
  * [22] S. Wiggins, The role of normally hyperbolic invariant manifolds (NHIMS) in the context of the phase space setting for chemical reaction dynamics. Regul. Chaot. Dyn. 21, 621 (2016)
  * [23] S. Wiggins, L. Wiesenfeld, C. Jaffé, and T. Uzer, Impenetrable Barriers in Phase-Space, Phys. Rev. Lett. 86, 5478 (2001)
  * [24] J. Eldering, Normally Hyperbolic Invariant Manifolds. The noncompact case.  ( Atlantis Press Paris, 2013 )
  * [25] S. Wiggins and V. J. G. García-Garrido, Painting the Phase Portrait of a Dynamical System with the Computational Tool of Lagrangian Descriptors. Notices of the American Mathematical Society, 69, 936 (2022)
  * [26] F. Gonzalez Montoya, The Classical Action as a Tool to visualize the Phase Space of Hamiltonian Systems, Dynamics, 4, 3 (2023)
  * [27] A. M Mancho, S. Wiggins, J. Curbelo, C. Mendoza, Lagrangian Descriptors: A Method for Revealing Phase Space Structures of General Time Dependent Dynamical Systems. Commun. Nonlinear Sci. Numer. Simul., 18, 3530 (2013)
  * [28] C. Lopesino, F. Balibrea-Iniesta, V. J. García-Garrido, S. Wiggins, A. M. Mancho, A Theoretical Framework for Lagrangian Descriptors, International Journal of Bifurcation and Chaos, 27, 1730001 (2017)
  * [29] R. H. Boyer and R. W. Lindquist, Maximal Analytic Extension of the Kerr Metric, J. Math. Phys. 8, 265 (1967)
  * [30] M. Firmbach, A. Bäcker, R. Ketzmerick, Partial barriers to chaotic transport in 4D symplectic maps, Chaos 33, 1 013125 (2023)
  * [31] J. Stöber, Classical and quantum transport in 4d symplectic maps. Technische Universität Dresden, PhD tesis (2023)
  * [32] F. Gonzalez Montoya and C. Jung, The numerical search for the internal dynamics of NHIMs and their pictorial representation, Physica D 436, 133330 (2022)
  * [33] F. Gonzalez Montoya and C. Jung, Visualizing the perturbation of partial integrability, J. Phys. A: Math. Theor, 48 , 43 (2015)
  * [34] F. Gonzalez Montoya and S. Wiggins, Phase space structure and escape time dynamics in a van der Waals model for exothermic reactions, Phys. Rev. E, 102 062203 (2020)
  * [35] F. Gonzalez Montoya, F. Borondo, C. Jung, Atom scattering off a vibrating surface: an example of chaotic scattering with three degrees of freedom, Commun. Nonlinear Sci. Numer. Simul., 90, 105282 (2020)
  * [36] G. Drotos, F. Gonzalez, C. Jung, The decay of a normally hyperbolic invariant manifold to dust in a three degrees of freedom scattering system, J. Phys. A: Math. Theor, 47, 045101 (2014)
  * [37] G. Drótos, F. González Montoya, C. Jung, T. Tél, Asymptotic observability of low-dimensional powder chaos in a three-degrees-of-freedom scattering system, Phys. Rev. E, 90, 2, 22906 (2014)
  * [38] L.C. Stein, N. Warburton, Location of the last stable orbit in Kerr spacetime. Phys. Rev. D, 101, 064007 (2020)
  * [39] C. Jung, O. Merlo, T. H. Seligman and W. P. K. Zapfe, The chaotic set and the cross section for chaotic scattering in three degrees of freedom, New J. Phys. 12, 103021 (2010)
  * [40] H. Waalkens, R. Schubert, S. Wiggins, Wigner’s dynamical transition state theory in phase space: classical and quantum, Nonlinearity, 21, R1 (2007)
  * [41] M. Katsanikas, S. Wiggins, The Generalization of the Periodic Orbit Dividing Surface in Hamiltonian Systems with three or more degrees of freedom–I, International Journal of Bifurcation and Chaos, 31, 2130028 (2021)
  * [42] M. Katsanikas, S. Wiggins, The Generalization of the Periodic Orbit Dividing Surface in Hamiltonian Systems with three or more degrees of freedom–II, International Journal of Bifurcation and Chaos, 31, 2150188 (2021)
  * [43] M. Katsanikas, S. Wiggins, The Generalization of the Periodic Orbit Dividing Surface in Hamiltonian Systems with three or more degrees of freedom–III, International Journal of Bifurcation and Chaos, 33, 2350088 (2023)
  * [44] R. Crossley, M. Agaoglou, M Katsanikas, S. Wiggins, From Poincaré Maps to Lagrangian Descriptors: The Case of the Valley Ridge Inflection Point Potential, Regul. Chaot. Dyn. 26, 2 (2021)
  * [45] M. Hillebrand, S. Zimper, A. Ngapasare, M. Katsanikas, S. Wiggins, C. Skokos, Quantifying chaos using Lagrangian descriptors, Chaos, 32, 12 (2022)
  * [46] S. Zimper, A. Ngapasare, M. Hillebrand, M. Katsanikas, S. Wiggins, C. Skokos, Performance of chaos diagnostics based on Lagrangian descriptors. Application to the 4D standard map, Physica D, 453 133833 (2022)


# A single layer representation of the scattered field for multiple scattering
problems

Didier Felbacq, Anthony Gourdin and Emmanuel Rousseau  University of
Montpellier<EMAIL_ADDRESS>
###### Abstract

The scattering of scalar waves by a set of scatterers is considered. It is
proven that the scattered field can be represented as the integral of a
density over an arbitrary smooth surface enclosing the scatterers. This is a
generalization of the series expansion over spherical harmonics and spherical
Bessel functions for spherical geometries. It allows an extension of the Fast
Multiple Algorithm to non spherical domains.

††: J. Phys. A: Math. Gen.

Keywords: scattering theory, scalar waves, integral representation

##  1 Introduction and setting of the problem

We consider the scattering of scalar waves by a set of obstacles in
$\mathbb{R}^{p},\,p=2,3$ in the harmonic regime with a time dependence of
$e^{-i\omega t}$. When the number of obstacles is very large the solving of
the scattering problem requires the use of an efficient algorithm, such as the
Fast Multipole Method [1]. In the present work, we show that the scattered
field can be represented by an integral supported by a surface enclosing the
obstacles. This approach allows a drastic reduction of the number of unknowns
and an extension of the Fast Multipole Method. The theoretical results are
illustrated by simple numerical examples in the last section.

Figure 1: Sketch of the scattering problem under study.

Let us specify a few notations. The unit sphere of $\mathbb{R}^{p}$ is denoted
$S^{p-1}$. For $\mathbf{x}\in\mathbb{R}^{p}$, we denote
$x=\left|\mathbf{x}\right|,\,\hat{x}=\mathbf{x}/x$ and $k=\omega/c$. We denote
$\mathcal{H}_{a}u=\Delta u+k^{2}au$, where the potential $a$ belongs to
$L^{\infty}(\mathbb{R}^{p})$. The fundamental solution $g^{+}$ of the
Helmholtz equation: $\mathcal{H}_{1}g^{+}=\delta_{0}$ with outgoing wave
condition is: $g^{+}(\mathbf{x})=-\frac{1}{4\pi x}e^{ikx}$ for $p=3$ and
$g^{+}(\mathbf{x})=-\frac{i}{4}H_{0}^{(1)}(kx)$ for $p=2$. The Green function
with the incoming wave condition is denoted $g^{-}(\mathbf{x})$. Explicitly:
$g^{-}(\mathbf{x})=-\frac{1}{4\pi x}e^{-ikx}$ for $p=3$, and
$g^{-}(\mathbf{x})=-\frac{i}{4}H_{0}^{(2)}(kx)$ for $p=2$. The functions
$H_{0}^{(1)}$ and $H_{0}^{(2)}$ are the Hankel functions of first and second
type [2].

Let us consider the following time-harmonic scattering problem. Let $\Omega$
be a bounded domain of $\mathbb{R}^{p}$ with boundary $\partial\Omega=\Gamma$,
containing a collection of scatterers (see Figure 1). The scatterers are
characterized by a potential $a$ such that $a-1$ has a compact support
$K\subset\Omega$. For a given incident field $u^{\rm inc}(x)$ satisfying the
Helmholtz equation: $\mathcal{H}_{1}u^{\rm inc}=0$, the scattering problem
consists in finding the scattered field $u^{s}(\mathbf{x})$ such that the
total field $u=u^{\rm inc}+u^{s}$ satisfies:

$\mathcal{H}_{a}u=0\,,$

and $u^{s}$ satisfies a radiation condition at infinity:
$\partial_{n}u^{s}-iku^{s}=o\left(x^{-1}\right)$ and
$u^{s}(\mathbf{x})=O\left(x^{-1}\right)$ when $x\to\infty$.

This scattering problem has a unique solution, as stated in the following
lemma:

######  Lemma 1.

The scattered field $u^{s}$ exists and is unique. There is a linear operator
$\mathcal{T}$, the scattering amplitude, relating $u^{\rm inc}$ to $u^{s}$:
$u^{s}=\mathcal{T}(u^{\rm inc})$.

###### Proof.

$\mathcal{H}_{a}(u^{s})=(\mathcal{H}_{1}-\mathcal{H}_{a})(u^{\rm inc})$ and
$\mathcal{V}\equiv\mathcal{H}_{1}-\mathcal{H}_{a}$ is null outside the compact
region $K$. Then :
$\mathcal{H}_{1}(u^{s})=\mathcal{V}(u^{s})+\mathcal{V}(u^{\rm inc})$ and thus:
$u^{s}=(1-\mathcal{G}_{1}\mathcal{V})^{-1}\mathcal{G}_{1}\mathcal{V}(u^{\rm
inc})$ where the inverse operator $\mathcal{G}_{1}=\mathcal{H}_{1}^{-1}$ is an
integral convolution operator with kernel $g^{+}$. ∎

The existence of the resolvent operator is classical although rather subtle
(see for instance [3, 4]). This provides a decomposition of the total field in
the form:

$u=u^{\rm inc}+\mathcal{T}(u^{\rm inc}).$

Let $B_{e}=B(O,R_{e})$ be the smallest ball with center $O$ containing
$\Omega$ and $B_{i}=B(O,R_{i})$ the largest ball, with center $O$, contained
in $K$. A modal expansion for the scattered field is valid outside $B_{e}$:

$u^{s}(\mathbf{x})=\left\\{\begin{array}[]{l}\sum_{n,m}u^{s}_{nm}h^{(1)}_{n}(kx)Y_{n}^{m}(\hat{x})\hbox{
for }p=3\\\ \sum_{n}u^{s}_{n}H^{(1)}_{n}(kx)e^{in\theta}\hbox{ for
}p=2\end{array}\right.,\,x>R_{e}.$ (1)

Here, $h^{(1)}_{n}$ is the spherical Hankel function of first type and order
$n$ [2] and $\theta$ is the polar angle of $\mathbf{x}$ in $\mathbb{R}^{2}$.

Whether the functions defined by these series can be extended inside the ball
is a difficult problem known as Rayleigh hypothesis, it was essentially solved
in the 80’ [8, 9].

Note that, if there is only one scatterer, i.e. $a$ is constant inside $K$,
there is also a representation of the field inside $K$ by a series in the
following form:

$u(\mathbf{x})=\left\\{\begin{array}[]{l}\sum_{n,m}u^{s}_{nm}j_{n}(kx)Y_{n}^{m}(\hat{x})\hbox{
for }p=3\\\ \sum_{n}u^{s}_{n}J_{n}(kx)e^{in\theta}\hbox{ for
}p=2\end{array}\right.,\,x<R_{i}.$ (2)

Here, $J_{n}$ (resp. $j_{n})$ is the Bessel (resp. spherical Bessel) function
of order $n$ [2].

When $\overline{\Omega}=K$ and the boundary $\Gamma$ is a sphere, both series
can be matched on $\Gamma$, which leads to an explicit form of the scattering
coefficients. By considering the traces of the field on the boundary, one can
obtain a pseudo-differential operator relating the coefficient of the incident
field to that of the scattered field. In the case where $\Gamma$ is not a
sphere and $K$ is a proper subset of $\Omega$, this approach can be extended
by using an integral representation of the fields: this is the purpose of this
work. A pioneering work in that direction can be found in [5].

##  2 Integral representations of the incident and scattered fields

Our aim is to obtain a representation of the incident and scattered fields as
an integral supported by $\Gamma$. Let us first specify some notations.

For $u\in H^{1}(\Omega)$ (the Sobolev space of function of $L^{2}(\Omega)$
with gradient in $L^{2}(\Omega)^{2}$, see [3, chap. 2] for more results on
Sobolev spaces), the interior traces [3, chap. 2] of $u$ and its normal
derivative on $\Gamma$ are denoted by:

$\gamma^{-}(u)=\left.u\right|_{\Gamma},\,\gamma^{-}(\partial_{n}u)=\left.\partial_{n}u\right|_{\Gamma}$
(3)

For fields belonging to $H^{1}_{\rm loc}(\Omega\setminus\mathbb{R}^{p})$, we
denote the exterior traces by:

$\gamma^{+}(u)=\left.u\right|_{\Gamma},\,\gamma^{+}(\partial_{n}u)=\left.\partial_{n}u\right|_{\Gamma}.$
(4)

Given a field $u\in H^{1}_{\rm loc}(\mathbb{R}^{p})$, we denote $[f]_{\Gamma}$
the jump of $f$ across $\Gamma$, i.e.:

$[u]_{\Gamma}=\gamma^{+}(u)-\gamma^{-}(u)\hbox{ and
}[\partial_{n}u]_{\Gamma}=\gamma^{+}(\partial_{n}u)-\gamma^{-}(\partial_{n}u).$
(5)

In order to have an integral representation of the fields, we state a result
concerning the incident field. To do so, we first need a technical lemma. Note
that in the following, the proofs of the results are given for $p=3$ and can
be easily adapted for $p=2$ (or, in fact, any other dimension $>1$).

######  Lemma 2.

Given $\sigma\in H^{-1/2}(\Gamma)$, define:

$J[\sigma](\hat{x})=\int_{\Gamma}\sigma(\mathbf{x}^{\prime})e^{-ik\hat{x}\cdot\mathbf{x}^{\prime}}d\mathbf{x}^{\prime}.$
(6)

Assume that $k^{2}$ is not an eigenvalue of $-\Delta$ inside $\Omega$ with
Dirichlet boundary conditions on $\Gamma$. If $J[\sigma](\hat{x})=0$, then
$\sigma=0$. Moreover $J$ defines a bijection: $H^{-1/2}(\Gamma)\to
L^{2}(S^{2})$.

###### Proof.

Consider the unique function $v$ satisfying $\mathcal{H}_{1}v=0$ in
$\Omega\cup(\mathbb{R}^{3}\setminus\overline{\Omega})$ and the boundary
conditions:

$[v]_{\Gamma}=0,\,[\partial_{n}v]_{\Gamma}=\sigma.$

Then $v$ is represented in the following integral form:
$v(\mathbf{x})=\int_{\Gamma}g^{+}(\mathbf{x}-\mathbf{x}^{\prime})\sigma(\mathbf{x}^{\prime})d\mathbf{x}^{\prime}$.
Outside $B_{e}$, $v$ can be expanded in spherical harmonics:

$v(\mathbf{x})=\sum_{mn}v_{nm}h_{n}^{(1)}(kx)Y^{m}_{n}(\hat{x}).$

The spherical Hankel functions have the following asymptotic forms [2] as
$x\to\infty$:

$\displaystyle
h_{n}^{(1)}(kx)\equiv\sqrt{\frac{\pi}{2x}}H^{(1)}_{n+1/2}(x)\sim_{x\to\infty}\frac{e^{i\left(kx-(n+1)\pi/2\right)}}{kx},$
(7) $\displaystyle
h_{n}^{(2)}(kx)\equiv\sqrt{\frac{\pi}{2x}}H^{(2)}_{n+1/2}(x)\sim_{x\to\infty}\frac{e^{-i\left(kx-(n+1)\pi/2\right)}}{kx}.$
(8)

From these relations, we deduce that the following asymptotic behavior holds:
${v(\mathbf{x})\sim\frac{e^{ikx}}{kx}w(\hat{x})}$ with

$w(\hat{x})=\sum_{nm}v_{nm}e^{-i(n+1)\pi/2}Y_{n}^{m}(\hat{x}).$

Besides, using the asymptotic form of the Green function:

$g^{\pm}(\mathbf{x}-\mathbf{x}^{\prime})\sim_{x\to\infty}-\frac{e^{\pm
ikx}}{4\pi x}e^{\mp ik\mathbf{x}\cdot\mathbf{x}^{\prime}},$ (9)

we obtain:

$v(x)\sim_{x\to\infty}\frac{e^{ikx}}{kx}J[\sigma](\hat{x}),$

and thus $J[\sigma](\hat{x})=w(\hat{x})$. Consequently, the nullity of
$J[\sigma](\hat{x})$ implies that of $w$ and thus that of $v$ in
$\mathbb{R}^{3}\setminus\Omega$. Therefore $\gamma_{1}(v)=0$ and
$\gamma_{2}(v)$ is non zero iff $v$ is a solution of the Dirichlet problem,
therefore $v=0$ everywhere by the hypothesis on $k^{2}$, and thus $\sigma=0$.
The surjectivity is obtained as follows. Take a function $\phi\in
L^{2}(S^{2})$ and expand it in spherical harmonics:

$\phi(\hat{x})=\sum_{nm}\phi_{nm}Y_{n}^{m}(\hat{x}).$

Then construct a field:

$u(\mathbf{x})=\sum_{nm}\phi_{nm}h_{n}^{(1)}(kx)Y_{n}^{m}(\hat{x}).$

For a fixed $x>0$, this series converges in $L^{2}(S^{2})$, since
$(h_{n}^{(1)}(kx))_{n}$ is a bounded sequence for every $x>0$. Finally, define
the field $\tilde{u}$ equals to $u$ outside $\Omega$ and satisfying:

$\mathcal{H}_{1}\tilde{u}=0\hbox{ in
}\Omega,\gamma_{1}(\tilde{u})=\gamma_{1}(u)\hbox{ on }\Gamma.$

Then $\tilde{u}$ satisfies:
$\mathcal{H}_{1}\tilde{u}=[\partial_{n}\tilde{u}]_{\Gamma}\delta_{\Gamma}$ and
it holds:

$\tilde{u}(\mathbf{x})=\int_{\Gamma}[\partial_{n}\tilde{u}]_{\Gamma}g^{+}(\mathbf{x}-\mathbf{x}^{\prime})d\mathbf{x}^{\prime},\,\mathbf{x}\in\mathbb{R}^{3}\setminus\Omega.$

Therefore, we obtain the existence of
$\sigma=[\partial_{n}\tilde{u}]_{\Gamma}\in H^{-1/2}(\Gamma)$. ∎

Using this lemma, we are now in a position to prove the following result that
provides a representation of the incident field as an integral over $\Gamma$:

######  Theorem 1.

The incident field can be represented in the form:

$u^{\rm inc}(\mathbf{x})=i\int_{\Gamma}\sigma^{\rm
inc}(\mathbf{x}^{\prime})\Im(g^{+}(\mathbf{x}-\mathbf{x}^{\prime}))dr^{\prime},\,\mathbf{x}\in\mathbb{R}^{p}$

where $\sigma^{\rm inc}$ belongs to $H^{-1/2}(\Gamma)$.

###### Proof.

The incident field can be expanded in spherical harmonics in the form ${u^{\rm
inc}=u^{+}+u^{-}}$, where:

$u^{+}(\mathbf{x})=\frac{1}{2}\sum_{nm}i_{nm}h_{n}^{(1)}(kx)Y_{n}^{m}(\hat{x}),\,u^{-}(\mathbf{x})=\frac{1}{2}\sum_{nm}i_{nm}h_{n}^{(2)}(kx)Y_{n}^{m}(\hat{x}).$

Using the asymptotic forms of the spherical Hankel functions (7,8) we obtain
the existence of two functions $u^{\pm}_{\infty}(\hat{x})$ defined on $S^{2}$
and such that:

$u^{\rm
inc}(x)\sim\frac{e^{ikx}}{kx}u^{+}_{\infty}(\hat{x})+\frac{e^{-ikx}}{kx}u^{-}_{\infty}(\hat{x})$

Explicitly, these functions are given by:

$u^{+}_{\infty}(\hat{x})=\frac{1}{2}\sum_{nm}i_{nm}e^{-i(n+1)\pi/2}Y_{n}^{m}(\hat{x}),\,u^{-}_{\infty}(\hat{x})=\frac{1}{2}\sum_{nm}i_{nm}e^{i(n+1)\pi/2}Y_{n}^{m}(\hat{x}).$

Since:

$e^{-i(n+1)\pi/2}=(-1)^{n+1}e^{i(n+1)\pi/2}\hbox{ and
}Y_{n}^{m}(-\hat{x})=(-1)^{n}Y_{n}^{m}(\hat{x}),$

we have that:

$\displaystyle
u^{-}_{\infty}(-\hat{x})=\frac{1}{2}\sum_{nm}i_{nm}e^{i(n+1)\pi/2}Y_{n}^{m}(-\hat{x})=\frac{1}{2}\sum_{nm}i_{nm}(-1)^{n}e^{i(n+1)\pi/2}Y_{n}^{m}(\hat{x})$
$\displaystyle=-\frac{1}{2}\sum_{nm}i_{nm}e^{-i(n+1)\pi/2}Y_{n}^{m}(\hat{x})=-u^{+}_{\infty}(\hat{x}).$

Consider now the field $\tilde{u}^{\rm inc}$, defined by the following
integral:

$\tilde{u}^{\rm inc}(\mathbf{x})=\frac{1}{2}\int_{\Gamma}\left[\sigma^{\rm
inc}(\mathbf{x}^{\prime})g^{+}(\mathbf{x}-\mathbf{x}^{\prime})-\sigma^{\rm
inc}(\mathbf{x}^{\prime})g^{-}(\mathbf{x}-\mathbf{x}^{\prime})\right]d\mathbf{x}^{\prime}.$

Due to the continuity properties of a single layer potential [6], it satisfies
$\mathcal{H}_{1}\tilde{u}^{\rm inc}=0$ . Using (9), it is given asymptotically
as $x\to\infty$, by:

$\tilde{u}^{\rm inc}(\mathbf{x})\sim\frac{e^{ikx}}{kx}\int_{\Gamma}\sigma^{\rm
inc}(\mathbf{x}^{\prime})e^{-ik\hat{x}\cdot\mathbf{x}^{\prime}}d\mathbf{x}^{\prime}-\frac{e^{-ikx}}{kx}\int_{\Gamma}\sigma^{\rm
inc}(\mathbf{x}^{\prime})e^{ik\hat{x}\cdot\mathbf{x}^{\prime}}d\mathbf{x}^{\prime}$

Therefore, the existence of $\sigma^{\rm inc}$ satisfying:
$u^{+}_{\infty}(\hat{x})=\int_{\Gamma}\sigma^{\rm
inc}(\mathbf{x}^{\prime})e^{-ik\hat{x}\cdot\mathbf{x}^{\prime}}d\mathbf{x}^{\prime}$
follows from lemma (2) and the second relation:

$u^{-}_{\infty}(\hat{x})=-\frac{e^{-ikx}}{kx}\int_{\Gamma}\sigma^{\rm
inc}(\mathbf{x}^{\prime})e^{ik\hat{x}\cdot\mathbf{x}^{\prime}}d\mathbf{x}^{\prime}$

is fulfilled thanks to:

$u^{-}_{\infty}(\hat{x})=-u^{+}_{\infty}(-\hat{x})=-\int_{\Gamma}\sigma^{\rm
inc}(\mathbf{x}^{\prime})e^{ik\hat{x}\cdot\mathbf{x}^{\prime}}d\mathbf{x}^{\prime}.$

We conclude from Rellich lemma [3, p. 74] that $u^{\rm inc}=\tilde{u}^{\rm
inc}$, since both fields satisfy the same equation and the same asymptotic
behavior at infinity. The integral expression follows by noting that $g^{+}$
and $g^{-}$ are complex conjugated functions.

∎

Our next result is that the scattered field can also be represented by an
integral over $\Gamma$:

######  Theorem 2.

There exists $\sigma^{s}\in H^{-1/2}(\Gamma)$ such that:

$u^{s}(\mathbf{x})=\int_{\Gamma}\sigma^{s}(\mathbf{x}^{\prime})g^{+}(\mathbf{x}-\mathbf{x}^{\prime})d\mathbf{x}^{\prime},\mathbf{x}\in\mathbb{R}^{p}\setminus\Omega.$
(10)

###### Proof.

Consider the field $\tilde{u}^{s}$ that is equal to $u^{s}$ outside $\Omega$
and that satisfies the following problem inside $\Omega$:

$\mathcal{H}_{1}\tilde{u}^{s}=0\hbox{ inside
}\Omega,\,\left.\tilde{u}^{s}\right|_{\Gamma}=\left.u^{s}\right|_{\Gamma}.$

Over $\mathbb{R}^{p}$, it satisfies:
$\mathcal{H}_{1}\tilde{u}^{s}=[\partial_{n}\tilde{u}^{s}]_{\Gamma}\delta_{\Gamma}.$
Given the outgoing wave condition at infinity, this gives:

$u^{s}(\mathbf{x})=\int_{\Gamma}[\partial_{n}\tilde{u}^{s}]_{\Gamma}g^{+}(\mathbf{x}-\mathbf{x}^{\prime})d\mathbf{x}^{\prime},$

and consequently the existence of the density
$\sigma^{s}=[\partial_{n}\tilde{u}^{s}]_{\Gamma}$ belonging to
$H^{-1/2}(\Gamma)$. ∎

We can now deduce the following representation result for the total field:

######  Corollary 1.

The total field $u$ can be written in the form

$u(\mathbf{x})=u^{\rm tot,+}(\mathbf{x})+u^{\rm tot,-}(\mathbf{x}),$ (11)

where

$u^{\rm
tot,+}(\mathbf{x})=\int_{\Gamma}\sigma^{+}(\mathbf{x}^{\prime})g^{+}(\mathbf{x}-\mathbf{x}^{\prime})d\mathbf{x}^{\prime},\,u^{\rm
tot,-}(\mathbf{x})=\int_{\Gamma}\sigma^{-}(\mathbf{x}^{\prime})g^{-}(\mathbf{x}-\mathbf{x}^{\prime})dx^{\prime},$
(12)

and $\sigma^{+},\sigma^{-}$ belong to $H^{-1/2}(\Gamma)$.

###### Proof.

This is a direct consequence of theorem 1 and theorem 2. Putting the scattered
field and the incident field together, we obtain:

$\sigma^{+}=[\partial_{n}\tilde{u}^{s}]_{\Gamma}+\frac{1}{2}\sigma^{\rm
inc},\,\sigma^{-}=-\frac{1}{2}\sigma^{\rm inc}.$ (13)

∎

##  3 Discussion and a numerical example

###  3.1 A numerical example

Let us consider the scattering of an electromagnetic plane wave in $E_{||}$
polarization by a collection of cylinders contained in a domain whose cross
section $\Omega$ is bounded by an astroid $\Gamma$ (cf. Figure 2).

Figure 2: One quarter of the domain $\Omega$: it is a smooth astroid filled
with small dielectic rods.

For numerical purpose, the units are that defined by the wavelength and we
choose $\lambda=1$. The astroid contains $N=4657$ rods with relative
permittivity 12. We use the multiple scattering approach described in [7].

Figure 3: Normalized Discrete Fourier Transform coefficients of the sequence
of coefficients $(s_{k})$.

We assume that the rods, at positions $(\mathbf{x}_{k})$, are small enough
that they can each be characterized by only one scattering coefficient
$s^{0}_{q}$ [7]. We denote $\mathbf{x}=(x_{1},x_{2})$.

Figure 4: Map of the normalized modulus of the total field outside the
astroid, computed by a direct summation over all the scatteres contained
inside $\Omega$. Figure 5: Map of the normalized modulus of the total field
outside the astroid, computed by using the single layer representation of the
scattered field.

At this step the scattered field is represented, everywhere outside the
cylinders, by a sum over the rods in the form:

$u^{s}(\mathbf{x})=\sum_{q=1}^{N}s^{0}_{q}H^{(1)}_{0}(k|\mathbf{x}-\mathbf{x}_{q}|).$
(14)

The coefficients $\hat{s}=(s^{0}_{q})_{1,\ldots,N}$ are determined from the
multiple scattering theory [7]. The scattering coefficient $s^{0}_{q}$ is
related to the local incident field $u^{\rm inc,loc}_{q}$ through the
scattering amplitude $t^{0}_{q}$: $s^{0}_{q}=t^{0}_{q}u^{\rm inc,loc}_{q}$.
For a circular dielectric rod of radius $r_{q}$ and relative permittivity
$\varepsilon_{q}=\nu_{q}^{2}$, it holds:

$t^{0}_{q}=-\frac{J_{1}(kr_{q})J_{0}(k\nu_{q}r_{q})-\nu_{q}J_{0}(kr_{q})J_{1}(k\nu_{q}r_{q})}{H^{(1)}_{1}(kr_{q})J_{0}(k\nu_{q}r_{q})-\nu_{q}H^{(1)}_{0}(kr_{q})J_{1}(k\nu_{q}r_{q})}.$

The local incident field is:

$u^{\rm inc,loc}_{q}=u^{\rm inc}(\mathbf{x}_{q})+\sum_{j\neq
q}s^{0}_{j}H_{0}^{(1)}(|\mathbf{x}_{q}-\mathbf{x}_{j}|).$ (15)

Therefore, it holds:

$s^{0}_{q}=t^{0}_{q}\left(u^{\rm inc}(\mathbf{x}_{q})+\sum_{j\neq
q}s^{0}_{j}H_{0}^{(1)}(|\mathbf{x}_{q}-\mathbf{x}_{j}|)\right).$ (16)

Let us denote $\overline{\overline{t}}$ the diagonal matrix defined by
$\overline{\overline{t}}=diag(t^{0}_{1},\ldots,t^{0}_{N})$ and
$\overline{\overline{h}}$ the matrix with entries
$h_{ij}=H_{0}^{(1)}(|\mathbf{x}_{i}-\mathbf{x}_{j}|)$ for $i\neq j$ and
$h_{ii}=0$. The coefficients $\hat{s}$ are obtained by solving the system:

$\left((\overline{\overline{t}})^{-1}-\overline{\overline{h}}\right)\hat{s}=\hat{u}^{\rm
inc}$ (17)

where $\hat{u}^{\rm inc}=(u^{\rm inc}(\mathbf{x}_{q}))_{1,\ldots,N}$.

Figure 6: Modulus of the total field on a curve deduced from $\Gamma$ by a
homothety of ratio 1.3. The red curve corresponds to the total field
reconstructed by means of the single layer representation and the blue circles
correspond to the total field computed by a direct summation over all the
scatterers.

The point is to be able to represent the scattered field by means of a single
layer potential as explained in section (2). The unicity of $\sigma^{s}$ being
ensured by theorem (1), it can be determined by solving an integral equation
of the first kind by imposing

$\int_{\Gamma}\sigma^{s}(\mathbf{x}^{\prime})H^{(1)}_{0}(\mathbf{x}-\mathbf{x}^{\prime})d\mathbf{x}^{\prime}=u^{s}(\mathbf{x}),\mathbf{x}\in\Gamma.$

In order to do so numerically, that is, to obtain a discretized version of the
density $\sigma^{s}$, we simply write a discrete version of the integral:

$u^{s}_{P}(\mathbf{x})=\sum_{p=1}^{P}\sigma^{s}_{p}H^{(1)}_{0}(k|\mathbf{x}-\mathbf{y}_{p}|),$
(18)

and the points $(\mathbf{y}_{p})$ are put uniformly on $\Gamma$. On $\Gamma$,
the scattered field can be written:

$u^{s}_{p}\equiv
u^{s}(\mathbf{y}^{\prime}_{p})=\sum_{q=1}^{N}s^{0}_{q}H^{(1)}_{0}(k|\mathbf{y}^{\prime}_{p}-\mathbf{x}_{q}|),\,\mathbf{y}^{\prime}_{p}\in\Gamma,p=1,\ldots,P.$

The second set of points $(\mathbf{y}^{\prime}_{p})$ is different from
$(\mathbf{y}_{p})$ in order to avoid the $0$ singularity of the Hankel
function. In matrix form, this reads as:

$\overline{\overline{H}}\hat{s}=\hat{u}^{s}$ (19)

where $\hat{u}^{s}=(u^{s}_{p})$. Let us remark that $\overline{\overline{H}}$
is a $P\times N$ matrix. Then a square linear system is obtained by writing:

$\sum_{p=1}^{P}\sigma^{s}_{p}H^{(1)}_{0}(k|\mathbf{y}^{\prime}_{p^{\prime}}-\mathbf{y}_{p}|)=u^{s}_{p^{\prime}},\,p^{\prime}=1\ldots
P.$

In matrix form, this can be written:

$\overline{\overline{I}}\hat{\sigma}^{s}=\hat{u}^{s},$ (20)

where: $\hat{\sigma}^{s}=(\sigma^{s}_{p})$ and $\overline{\overline{I}}$ is a
$P\times P$ matrix with entries
$I_{pp^{\prime}}=H^{(1)}_{0}(k|\mathbf{y}^{\prime}_{p^{\prime}}-\mathbf{y}_{p}|)$.
Finally, we obtain a matrix $\overline{\overline{B}}$ relating $\hat{s}^{s}$
to $\hat{\sigma}^{s}$:

$\overline{\overline{B}}=(\overline{\overline{I}})^{-1}\overline{\overline{H}}.$
(21)

This matrix is a $P\times N$ matrix.

Heuristically, the number $P$ can be determined by recalling that the function
$\sigma^{s}$ is periodic, since it is defined on a bounded curve.
Consequently, the computation of the DFT of $(\sigma^{s}_{p})_{p\in\\{1,\ldots
P\\}}$ can indicate whether the approximation is good, by checking the
decreasing of the Fourier coefficients. This is examplified in Figure 3 where
we have computed the DFT of the finite sequence $(\sigma^{s}_{p})$. It is
important to have this criterium, since the discrete values of the density
$(\sigma^{s}_{p})$ do not have a decreasing behavior with $P$. The final value
is $P=160$.

We are able to reconstruct with a very good precision the diffracted field. In
Figure 4 we have plotted a map of the total field outside the region where the
scatterers are contained, obtained by summing the contributions of the
dielectric rods, and in Figure 5 it is the reconstructed field. Both fields
have been normalized so that their maximal value is equal to 1, in order to
have the same color scale. For a more direct comparison, in Figure 6 we have
plotted the total field on a curve deduced from $\Gamma$ by a homothety of
ratio 1.3. We stress that, thanks to this approach the representation of the
scattered field is now ensured by $160$ terms instead of 4657.

###  3.2 Extension of the Fast Multipole Method

It is important to remark that the single layer representation involves a
surface $\Gamma$, enclosing the scatterers, that can be chosen at will. By
this we mean that, given a set of scatterers and any smooth enough surface
$\Gamma$ enclosing this set, the field scattered by this set and the incident
field can be represented by an integral over $\Gamma$. This result is in fact
a generalization of the expansion over spherical harmonics and spherical
Bessel functions to an arbitrary surface. As a consequence, it is possible to
split a given set of scatterers into several subset, apply multiple scattering
theory to each smaller subset then use the representation by the single layers
to couple the subsets in between them. In order to do so, an iterative
algorithm is to be used. Let us be more specific: we put $p=2$ and consider
simply two subsets $O^{1}=\cup_{j=1}^{N^{1}}\Omega^{1}_{j}$ and
$O^{2}=\cup_{j=1}^{N^{2}}\Omega^{2}_{j}$ and two surfaces $\Gamma_{1}$ and
$\Gamma_{2}$ enclosing respectively $O^{1}$ and $O^{2}$ (see Figure 7).

Figure 7: Sketch of the scattering problem with two subsets.

The incident field $u^{\rm inc}$ illuminates $O^{1}$ and $O^{2}$. As in the
preceding section, we assume for simplicity that the wavelength is large
enough that the obstacle $\Omega^{\alpha}_{j}$ ($\alpha=1,2$) can be
considered to be a point at coordinate $\mathbf{x}_{j}^{\alpha}$, and that the
field scattered by $\Omega^{\alpha}_{j}$ reads as:

$u^{s,\alpha}_{j}(\mathbf{x})=s^{0,\alpha}_{j}H_{0}^{(1)}(k|\mathbf{x}-\mathbf{x}^{\alpha}_{j}|)$

For each scatterer $\Omega^{1}_{j}\in O^{1}$, the incident field is the sum of
the ”true” incident field $u^{\rm inc}(\mathbf{x}^{1}_{j})$ and the field
coming from the other subset $O^{2}$ and given by the discrete single layer
representation:

$u^{s,2}(\mathbf{x}^{1}_{j})=\sum_{p=1}^{P_{2}}\sigma_{p}^{s,2}(\mathbf{y}_{p})H_{0}^{(1)}(k|\mathbf{x}^{1}_{j}-\mathbf{y}_{p}|).$
(22)

Of course, there is the same set of relations obtained by making the switching
$1\leftrightarrow 2$. Here the local incident field is therefore:

$u^{\rm inc,loc}_{j}=u^{\rm
inc}(\mathbf{x}^{1}_{j})+u^{s,2}(\mathbf{x}^{1}_{j}).$ (23)

We denote the diagonal matrix
$\overline{\overline{t^{\alpha}}}=diag(t^{0,\alpha}_{1},\ldots,t^{0,\alpha}_{N_{\alpha}})$,
$\hat{u}^{s,2}=(u^{s,2}(x^{1}_{j}))$ (resp.
$\hat{u}^{s,1}=(u^{s,1}(\mathbf{x}^{2}_{j}))$) and $\hat{u}^{{\rm
inc},{\alpha}}=(u^{\rm inc}(\mathbf{x}^{\alpha}_{j}))$. As in the preceding
section, the operator that relates the scattering coefficients
$\hat{s}^{\alpha}=(s^{0,\alpha}_{j})$ to the discretized density
$\hat{\sigma}^{s,\alpha}=(\sigma^{s,\alpha}_{p})$ is denoted by
$\overline{\overline{B}}^{\alpha}$ (cf. (21)):

$\hat{\sigma}^{\alpha}=\overline{\overline{B}}^{\alpha}\,\hat{s}^{\alpha}.$
(24)

It is a $P_{\alpha}\times N_{\alpha}$ matrix. Finally,
$\overline{\overline{R}}^{\alpha}$ is the matrix that relates
$\hat{\sigma}^{s,\alpha}$ to $\hat{u}^{s,\alpha}$:

$\hat{u}^{s,\alpha}=\overline{\overline{R}}^{\alpha}\,\hat{\sigma}^{\alpha}.$
(25)

Explicitly (taking $\alpha=2$):

$\hat{u}^{s,2}(\mathbf{x}^{1}_{j})=\sum_{p=1}^{P_{2}}H_{0}^{(1)}(|\mathbf{x}^{1}_{j}-\mathbf{y}^{2}_{p}|)\,\sigma^{s,2}_{p},\,j=1,\ldots
N_{1}.$ (26)

Therefore $\overline{\overline{R}}^{2}$ as entries:
$R^{2}_{ij}=H_{0}^{(1)}(|\mathbf{x}^{1}_{i}-\mathbf{y}^{2}_{j}|),\,i=1,\ldots
N_{1},\,j=1,\ldots P_{2}$. It is a $N_{1}\times P_{2}$ matrix. Finally, the
linear system to be solved is:

$\left[\begin{pmatrix}(\overline{\overline{t^{1}}})^{-1}&0\\\
0&(\overline{\overline{t^{2}}})^{-2}\end{pmatrix}-\begin{pmatrix}\overline{\overline{h}}^{1}&\overline{\overline{R}}^{2}\overline{\overline{B}}^{2}\\\
\overline{\overline{R}}^{1}\overline{\overline{B}}^{1}&\overline{\overline{h}}^{2}\end{pmatrix}\right]\begin{pmatrix}\hat{s}^{1}\\\
\hat{s}^{2}\end{pmatrix}=\begin{pmatrix}\hat{u}^{\rm inc,1}\\\ \hat{u}^{\rm
inc,2}\end{pmatrix}.$ (27)

The numerical gain here lies in the diagonal terms
$\overline{\overline{R}}^{\alpha}\overline{\overline{B}}^{\alpha}$. If we were
to use directly the multiple scattering theory, we would have to compute these
terms by coupling directly each cylinder to every other cylinder, which would
involve the direct computation of the matrix with entries
$H_{0}^{(1)}(|\mathbf{x}^{1}_{i}-\mathbf{x}^{2}_{j}|)$ of size $N_{1}\times
N_{2}$. Here, we have to compute directly the matrix
$\overline{\overline{h}}^{1}$ (resp. $\overline{\overline{h}}^{2}$) of size
$N_{1}\times N_{1}$ (resp. $N_{2}\times N_{2}$) and then the coupling between
the two subsets is ensured by the terms
$\overline{\overline{R}}^{1}\overline{\overline{B}}^{1}$ and
$\overline{\overline{R}}^{2}\overline{\overline{B}}^{2}$. While these matrices
are of course of size $N_{2}\times N_{1}$ (resp. $N_{1}\times N_{2}$) they are
obtained as the product of matrices of size $N_{2}\times P_{1}$ and
$P_{1}\times N_{1}$ (resp. $N_{1}\times P_{2}$ and $P_{2}\times N_{1}$ ).

Figure 8: Two elliptic subsets containing each $3234$ dielectric cylinders.
The red ellipsis indicates where the field is computed in Figure 10.

Let us give a numerical example. We consider two subsets of dielectric rods
contained in two elliptical domains, as depicted in Figure 8. The incident
field is a plane wave $u^{\rm inc}(\mathbf{x})=e^{-ikx_{2}}$ . As in the first
numerical example, the length unit is that of the wavelength and we choose
$\lambda=10$. There are $3234$ rods with radius $0.015$ and relative
permittivity $\varepsilon=10$ in each domain. The map of the field is given in
Figure 9. We have computed the maps by using the single layer representation
(left panel) and by using the multiple scattering theory for the entire set of
rods in the right panel. We have used $P_{1}=P_{2}=100$ points to compute the
single layer representations (22). The maps and the scattering coefficients
agree to a precision below $0.3\%$ (in $L^{2}$ norm for the entire region
covered by the map).

Figure 9: Maps of the modulus field. On the left (a) panel, the field is
computed by using the extended Fast Multipole Method, on the right (b) panel,
it is computed by using directly the multiple scattering theory for the entire
set of cylinders.

In Figure 10, we have plotted the modulus of the scattered field on the red
ellipsis plotted in Figure 8. The fields coincide to $0.6\%$.

Figure 10: The curve in solid line corresponds to the field computed by using
directly the single layer representations of the fields scattered by each
subset (cf. (18)). The circles corresponds to the field computed by summing
over the contribution of each cylinder (cf. (14)).

The calculation time on a laptop for the single layer approach is around $2$
to $3$ times faster than for the direct multiple scattering approach. It is to
be noted that the number of points $P_{1,2}$ plays a negligible role in the
total calculation time: reducing $P_{1}=P_{2}$ to $30$ does not change the
calculation time up to the time fluctuations.

In conclusion, we have established a new way of representing the field
scattered by a large collection of object by using a single layer
representation. The scattered field is characterized by a density supported by
the boundary of a domain containing the scatterers. From a numerical point of
view, the gain lies in the number of parameters needed to represent the field.
Since the sources are supported by a region of codimension 1, much less
information is needed, as compared to a volumetric representation (of
codimension 0). This result in a drastic reduction of the number of values
required for representing the scattered field with a given precision. This
result is a generalization of the representation of the field by spherical
harmonics used in the Fast Multiple Method and extend this algorithm beyond
the spherical geometry.

## References

## References

  * [1] Gumerov N and Duraiswami R 2004 Fast Multipole Methods for the Helmholtz Equation in Three Dimensions (Elsevier: Amsterdam).
  * [2] Abramowitz M and Stegun I 1964 Handbook of Mathematical Functions (Dover: London).
  * [3] Cessenat M 1996 Mathematical Methods in Electromagnetism(World Scientific: Singapore).
  * [4] Melrose R 1995 Geometric scattering theory (Cambridge University Press: Cambridge).
  * [5] Maystre D 2006 Electromagnetic scattering by a set of objects: an integral methods based on scattering properties PIER 57 55-84.
  * [6] Colton D and Kress R 2013 Integral Equation Method in Scattering Theory (SIAM: Philadelphia).
  * [7] Felbacq D, Tayeb G and Maystre D 1995 Scattering by a random set of parallel cylinders J. Opt. Soc. Am. A 11 2526-2538.
  * [8] Millar R F 1980 The Analytic Continuation of Solutions to Elliptic Boundary Value Problems in Two Independent Variables J. Math. Anal. Appl. 76 498-515.
  * [9] Maystre D and Cadilhac M 1985 Singularities of the continuation of fields and validity of Rayleigh’s hypothesis J. Math. Phys. 26 2201-2204.


# A concise proof of Benford’s law

Luohan Wang  Bo-Qiang Ma

###### Abstract

This article presents a concise proof of the famous Benford’s law when the
distribution has a Riemann integrable probability density function and
provides a criterion to judge whether a distribution obeys the law. The proof
is intuitive and elegant, accessible to anyone with basic knowledge of
calculus, revealing that the law originates from the basic property of the
human number system. The criterion can bring great convenience to the field of
fraud detection.

###### keywords:

Benford’s law , First-digit law , Significant digit law , Proof , Criterion

††journal: Fundamental Research, 4 (2024) 842-845, by L. Wang & B.-Q. Ma,
doi:10.1016/j.fmre.2023.01.002

[inst1]organization=School of Physics,addressline=Peking University,
city=Beijing, postcode=100871, country=China

[inst2]organization=Center for High Energy Physics,addressline=Peking
University, city=Beijing, postcode=100871, country=China

[inst3]organization=Collaborative Innovation Center of Quantum
Matter,city=Beijing, country=China

##  1 Introduction

One may instinctively assume that the occurrence of the first digit of a
randomly selected number is uniformly distributed among 1 to 9. Nonetheless,
the Newcomb-Benford law, also known as the first-digit law or Benford’s law,
reveals a completely different fact. In many natural and human phenomena, the
distribution of the first significant digit of a random number follows a
logarithm-type law.

Early in 1881, Newcomb [1] noticed that the first pages of the logarithmic
table tend to wear out faster than the last ones. Consequently, he found that
the first significant figure is oftener 1 than any other digit, with the
frequency diminishing from 1 to 9, and figured out the logarithm distribution
of the first two significant digits of natural numbers. In 1938, by analyzing
over 20,000 numbers collected from 20 widely divergent sources, Benford [2]
found that these data show unbelievable adherence to a logarithm distribution,
and rediscovered the empirical law

$P(d)=\log_{10}\left(1+\frac{1}{d}\right),\quad d=1,2,\cdots,9,$ (1)

where $P(d)$ is the probability of a randomly selected number with a first
digit of $d$. However, no convincing proof was proposed by him at that time.

This bias in the frequency of occurrence of the leading digits of numbers is
evident in everything from the lengths of rivers, the areas of countries,
stock prices, and population numbers to death rates. In terms of mathematics,
Benford’s law appears in integer sequences, such as the celebrated Fibonacci
sequence [3] and the factorial sequence [4], and the distribution of prime
numbers follows a general Benford’s law [5]. The law also emerges in the field
of physics. The collection of many physics constants [6] conforms to this law.
The half-lives of unstable nuclei [7, 8, 9] and many other experimental
quantities [10] in nuclear physics also agree with the logarithm distribution.
Its applications in physics involve astronomy [11], dynamic system [12],
statistical physics [13, 14], nuclear physics [7, 8, 9, 10] and particle
physics [15]. For example, the Bose-Einstein distribution conforms to it
strictly whatever the temperature is, while the Fermi-Dirac distribution
fluctuates slightly in a periodic manner depending on the temperature of the
system [13].

The detection of frauds is one of the most prominent applications of the
Newcomb–Benford law for significant digits. In 1972, a seminal letter by
Varian [16] suggested the notion that Benford’s law can be applied in
detecting fraud data. In 1992, Nigrini proposed in his doctoral dissertation
to use the law to check for false accounting. Since then, the law has become
rather popular in auditing financial data, testing the fairness of electoral
processes [17], while a study anticipating the validity of this law in
detecting frauds in international trades was carried out in 2019 [18]. The law
also succeeds in application to many other domains, including social science
[19] and environmental science [20].

Nevertheless, Benford’s law has its scope of application. It works in a great
many natural date sets and some human phenomena, but most artificial data sets
and distributions that do not span several orders of magnitude (telephone
number, height, weight) or are artificially truncated do not follow Benford’s
law.

With regard to the underlying mechanism of Benford’s law, there have been
explanations from different views. However, due to the intriguing feature of
Benford’s law, the available explanations are not completely satisfactory.
They are either under a restricted condition or mathematically not rigorous.
The inchoate research focused on the cases where the samples are number sets
defined on $\mathbb{R^{+}}$, such as the set of natural integers. In 1976,
Cohen [21] raised the condition under which a measure on $\mathbb{N}$
satisfies Benford’s law. In 1992, Jech [22] stated the necessary and
sufficient conditions for a probability measure with a finitely additive
function defined on $\mathbb{R^{+}}$ to observe the law. The law also has
proof of physical versions due to its frequent appearance in physics. Lemon
[23] and Iafrate [24] derived this law using the maximum entropy model in
thermodynamics.

An explanation widely considered as rigorous is by Hill [25] in 1995. He
illustrated that when the distributions and the sample data are both selected
at random, the first digits of the combinations follow Benford’s law. However,
this work did not answer the question of why data from a specific distribution
obey this law and did not pay much attention to the internal reasons and
conditions of the law.

In addition, little attention has been paid to the circumstances where the
random variable set is continuous and has a continuous probability density
function defined on it. Studying this topic helps to explain why a continuous
distribution follows or violates the law. For example, among the common
distributions, the exponential distribution conforms well to the law, while
the uniform distribution and the normal distribution do not obey the law.
Research conducted by Engel and Leuenberger [26] proved that exponential
functions obey this law within an allowed error range. In 2019, a study at
Peking University gave a proof of the law when the probability density
function has an inverse Laplace transform [27, 28], revealing that the first
digit law originates from the basic property of the human number system.

In this article, first, we aim to derive this law for distributions with
Riemann integrable probability density functions on the positive real number
set, and the basic idea of the proof is inspired by Refs. [27, 28]. Our proof
is concise, intuitive, and accessible to anyone with basic knowledge of
calculus. Second, we estimate the potential error caused by a hypothesis in
the proof and discuss the criterion to judge the validity of this law in a
distribution.

##  2 Proof for the continuous case

###  2.1 The general proof of Benford’s law

Let the sample space of the variable $x$ be the positive real number set
$\mathbb{R^{+}}$, and let $f(x)$ be a continuous probability density function
of an arbitrary distribution defined on $\mathbb{R^{+}}$, with $f(x)$
satisfying the normalization condition. The set of positive numbers with the
first significant of $d\ (d=1,2,\cdots,9)$ is

$S(d)=\bigcup\limits_{n=-\infty}^{+\infty}[d\cdot 10^{n},(d+1)\cdot 10^{n}).$
(2)

Therefore, it is obvious that

$P(d)=\sum_{n=-\infty}^{+\infty}\int_{d\cdot 10^{n}}^{(d+1)\cdot
10^{n}}f(x)\,\text{d}x,$ (3)

where $P(d)$ is the probability of a number (base 10) chosen randomly from the
distribution $f(x)$ beginning with $d$. Evidently, $P(d)$ can also be written
as

$P(d)=\sum_{n=-\infty}^{+\infty}\int_{d\cdot 10^{n}}^{(d+1)\cdot
10^{n}}f(x)\,\text{d}x\,\Delta n,\quad(\Delta n=1).$ (4)

Then we will make an important assumption that the sum in the above equation
is approximately equal to the integral as follows, which is similar to Eq.
(11) in Refs. [27, 28]

$\displaystyle P(d)\approx\int_{-\infty}^{+\infty}\left[\int_{d\cdot
10^{n}}^{(d+1)\cdot 10^{n}}f(x)\,\text{d}x\right]\text{d}n,$ (5)

since $\Delta n$ is not exactly the differential d$n$, this equation does not
always hold. This explains why some $f(x)$ break Benford’s law when the
difference between the two equations cannot be ignored. We discuss this
situation in the next section.

Substitute the variable $n$ with $t=d\cdot 10^{n}$, $x$ with $x=ty$, therefore
$\text{d}n=\frac{1}{t\ln 10}\,\text{d}t$, and the same $P(d)$ in Eq. (5) can
be written as

$\displaystyle P(d)=\frac{1}{\ln
10}\int_{0}^{+\infty}\text{d}t\int_{1}^{1+\frac{1}{d}}f(ty)\,\text{d}y.$ (6)

######  Definition 1.

$\phi(y,t)=f(ty)$, and $\phi(y,t)$ is continuous on
$\left[1,1+\frac{1}{d}\right]\times[0,+\infty)$.

Apparently, the integral $\int_{0}^{+\infty}\phi(y,t)\,\text{d}t$ is uniform
convergent on $\left[1,1+\frac{1}{d}\right]$.

######  Lemma 1.

If $\phi(y,t)$ is continuous on $[a,b]\times[c,+\infty)$,
$a,b,c\in\mathbb{R^{+}}$, and the integral
$\int_{c}^{+\infty}\phi(y,t)\,\mathrm{d}t$ is uniform convergent on $[a,b]$,
then there is

$\int_{c}^{+\infty}\mathrm{d}t\int_{a}^{b}\phi(y,t)\,\mathrm{d}y=\int_{a}^{b}\mathrm{d}y\int_{c}^{+\infty}\phi(y,t)\,\mathrm{d}t,$
(7)

which is a basic theorem in calculus.

Therefore, according to Lemma 1,

$\displaystyle P(d)$ $\displaystyle=\frac{1}{\ln
10}\int_{0}^{+\infty}\text{d}t\int_{1}^{1+\frac{1}{d}}\phi(y,t)\,\text{d}y$
$\displaystyle=\frac{1}{\ln
10}\int_{1}^{1+\frac{1}{d}}\text{d}y\int_{0}^{+\infty}\phi(y,t)\,\text{d}t$
$\displaystyle=\frac{1}{\ln
10}\int_{1}^{1+\frac{1}{d}}\frac{1}{y}\,\text{d}y\int_{0}^{+\infty}f(x)\,\text{d}x$
$\displaystyle=\log_{10}\left(1+\frac{1}{d}\right)$ (8)

is exactly the logarithm relation in Benford’s law. The proof above reveals
that Benford’s law is essentially a mathematical principle. It arises simply
because of humans’ method of counting.

###  2.2 The extensional versions of Benford’s law

Apart from the commonest first digit law in the decimal system, Benford’s law
has other properties and extensional versions. The law is both scale-invariant
and base-invariant. The scale-invariance is the direct conclusion of the proof
above. By replacing $x$ with $ax$ ($a$ is a positive real constant) in Eq.
(3), we can find that the result is the same since the integrating range is
from zero to infinity.

In 1995, Hill [25] derived a general $i$th-significant digit law: For all
positive integers $k$, all $d_{1}\in\left\\{1,2,\cdots,9\right\\}$ and all
$d_{j}\in\left\\{0,1,\cdots,9\right\\}$, $j=1,2,\cdots,k$, the probability of
a number beginning with first $k$ digit, $D_{j}=d_{j}$ can be expressed as

$\displaystyle P$ $\displaystyle\left(D_{1}=d_{1},\cdots,D_{k}=d_{k}\right)$
$\displaystyle=\log_{10}\left[1+\left(\sum_{i=1}^{k}d_{i}\cdot
10^{k-i}\right)^{-1}\right].$ (9)

Here is an extension of Eq. (2.2) to a number with continuous coefficient in
base $m$. Suppose $x=k\cdot m^{n},\ m\in\mathbb{N^{+}},\ k\in[1,m)$, and $n$
is an integer. Specially, this is the scientific notation of $x$ when $m=10$.
Then use the same method mentioned above, we can calculate that the
probability of a number $x$ with its $k\in[a,b),\ a\geq 1,\ b<m$ is

$\displaystyle P_{m}(a,b)$
$\displaystyle=\sum_{n=-\infty}^{+\infty}\int_{a\cdot m^{n}}^{b\cdot
m^{n}}f(x)\,\text{d}x$ $\displaystyle=\log_{m}b/a.$ (10)

Clearly, this general version of Benford’s law is base-invariant and Eq. (2.2)
is a special case when $m=10,\ a=\sum\limits_{i=1}\limits^{k}d_{i}\cdot
10^{k-i},\ b=a+1$ of Eq. (10).

##  3 Discussions over the errors

In the proof above, we use an approximation as Eq. (5). Apparently, replacing
the sum with an integration is not universally feasible. When $n$ gets large,
the integral interval of $f(x)$ for each $n$ grows wider, and the integration
is more possibly to disagree with the original sum. Therefore, we should
discuss carefully the rationality of this replacement, as well as what kind of
$f(x)$ obeys Benford’s law.

######  Definition 2.

$g_{d}(x)$ is a function satisfying

$\displaystyle g_{d}(x)=\sum_{n=-\infty}^{+\infty}\left[\eta\left(x-d\cdot
10^{n}\right)-\eta\left(x-(d+1)\cdot 10^{n}\right)\right],$ (11)

which is the same function in Refs. [27, 28], and $\eta(x)$ is the Heaviside
step function

$\displaystyle\eta(x)=\left\\{\begin{aligned} 0\quad x<0\\\ 1\quad x\geq
0\end{aligned}\right..$ (12)

Therefore, $P(d)$ can also be rewritten as

$P(d)=\int_{0}^{+\infty}g_{d}(x)f(x)\,\text{d}x.$ (13)

######  Definition 3.

Define a function $\Delta(x)$ as follows:

$\Delta(x)=g_{d}(x)-\log_{10}\left(1+\frac{1}{d}\right).$ (14)

Denote the Error function $\mathrm{Er}(f)$ a functional of $f(x)$ which
represents the difference between the original sum and the integration in Eq.
(5), therefore

$\displaystyle\mathrm{Er}(f)$
$\displaystyle=\int_{0}^{+\infty}g_{d}(x)f(x)\,\text{d}x-\log_{10}\left(1+\frac{1}{d}\right)$
$\displaystyle=\int_{0}^{+\infty}f(x)\left[g_{d}(x)-\log_{10}\left(1+\frac{1}{d}\right)\right]\text{d}x$
$\displaystyle=\int_{0}^{+\infty}f(x)\Delta(x)\,\text{d}x,$ (15)

$\Delta(x)$ oscillates within
$\left[-\log_{10}\left(1+\frac{1}{d}\right),1-\log_{10}\left(1+\frac{1}{d}\right)\right]$
and the interval is proportional to $10^{n},\ n\in\mathbb{Z}$. Qualitatively,
if $f(x)$ strides across many orders of magnitude, the positive and negative
part of $\Delta(x)$ may be offset. However, if $f(x)$ is defined on a narrow
interval, $\mathrm{Er}(f)$ could be significant, that’s why Benford’s law is
broken when the distribution is narrow. It can also explain why artificial
data sets always contradict this law, since humans tend to make a uniform
distribution while uniform distributions violate Benford’s law.

To judge whether the replacement used in Eq. (5) is allowed, we only have to
calculate the integration in Eq. (3). If $\mathrm{Er}(f)$ is minute, $f(x)$
obeys Benford’s law. The examples in Refs. [27, 28] and Ref. [26] both verify
the effectiveness of this criterion. Here we take a simple case where
$f(x)=\lambda e^{-\lambda x},\ \lambda>0$ for example:

$\text{Er}(f)=\sum\limits_{n=-\infty}\limits^{+\infty}e^{-\lambda d\cdot
10^{n}}\left(1-e^{-\lambda\cdot
10^{n}}\right)-\log_{10}\left(1+\frac{1}{d}\right).$ (16)

According to the calculation in Ref. [26], $|\text{Er}(f)|\leq 0.03$ for any
$\lambda$ and $d$. To give an intuitive illustration, we plot a figure of the
difference between Benford’s distribution and the exact distributions of
exponential functions as shown in Fig.1.

Figure 1: Graph of $P(d)$ and the range of error. The red part is the maximum
$|\text{Er}(f)|$ of $f(x)$ with different $\lambda$ for each $d$, and the
black dots are Benford’s distributions.

Apparently, Benford’s law is a good estimate for exponential distributions.
This also explains the prevalence of Benford’s law in nature. A number of
distributions in nature are exponential, or can be written as superposition of
exponential functions, such as $f(x)=\int_{0}^{+\infty}c(t)e^{-tx}\text{d}t$.
The normalized probability density function $f(x)$ can be seen as the
superposition of $e^{-tx}$ with coefficient $c(t)$ in this expression, meaning
$f(x)$ has an inverse Laplace transformation. Therefore, it is intuitive to
guess that such $f(x)$ with positive $c(t)$ converges to Benford’s law, which
is exactly the conclusion of Refs. [27, 28].

What else, though the proof above is under the condition that $f(x)$ is
continuous, it can be extended to discontinuous $f(x)$. If $f(x)$ can be
approximated by a continuous function, in other words, for all
${\displaystyle\varepsilon>0}$, there exists a continuous function $h(x)$ on
$[0,+\infty]$ satisfying

$\int_{0}^{+\infty}\left|f(x)-h(x)\right|\text{d}x<\varepsilon,$ (17)

then the difference between $\mathrm{Er}(f)$ and $\mathrm{Er}(h)$ is also
infinitesimal. As long as $\mathrm{Er}(h)$ is a small number, $f(x)$ converges
to the logarithm distribution. $h(x)$ can be the Fourier expansion of $f(x)$
or a piecewise linear function, which is easy to calculate its
$\mathrm{Er}(h)$. Therefore, the condition of the proof is weakened as $f(x)$
is Riemann integrable on $\mathbb{R^{+}}$.

We also have to mention that our methodology does not work for all the
distributions. The method loses efficacy when $f(x)$ is not Riemann integrable
or the distribution does not have an asymptotic density. Related works in
these circumstances can be found in Refs. [21, 29, 30].

##  4 Summary

Over more than one hundred years, Benford’s law has drawn much attention in
different disciplines. Despite its long history, the law still lacks a
satisfactory explanation and is even considered as a mystery of nature. Does
nature favor logarithms? Does nature count at the base of natural exponential
functions instead of our natural number scale? Benford’s law has inspired so
many conjectures.

In this article, we prove Benford’s law for Riemann integrable probability
density functions defined on $\mathbb{R^{+}}$ and give the criterion to
estimate the accuracy of the law. Our proof is concise, intuitive, and
accessible to anyone with basic knowledge of calculus. With the criterion we
can easily determine whether a distribution with a continuous variable obeys
Benford’s law, thus such criterion can bring great convenience to the field of
fraud detection.

## Declaration of competing interest

The authors declare that they have no conflict of interests.

## Acknowledgments

This work is supported by National Natural Science Foundation of China (Grant
No. 12075003).

## References

  * [1] S. Newcomb, Note on the frequency of use of the different digits in natural numbers, Amer. J. Math. 4 (1-4) (1881) 39–40.  doi:10.2307/2369148.
  * [2] F. Benford, The law of anomalous numbers, Proceedings of the American Philosophical Society 78 (4) (1938) 551–572, http://www.jstor.org/stable/984802.
  * [3] A. Berger, T. P. Hill, An introduction to Benford’s law, Princeton University Press, 2015, https://www.jstor.org/stable/j.ctt1dr35m0.
  * [4] P. Diaconis, The distribution of leading digits and uniform distribution mod 1, The Annals of Probability 5 (1) (1977) 72 – 81.  doi:10.1214/aop/1176995891.
  * [5] B. Luque, L. Lacasa, The first-digit frequencies of prime numbers and riemann zeta zeros, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 465 (2107) (2009) 2197–2216.  doi:10.1098/rspa.2009.0126.
  * [6] J. Burke, E. Kincanon, Benford’s law and physical constants: The distribution of initial digits, American Journal of Physics 59 (10) (1991) 952–952.  doi:10.1119/1.16838.
  * [7] D. Ni, Z. Ren, Benford’s law and half-lives of unstable nuclei, The European Physical Journal A 38 (2008) 251–255.  doi:10.1140/epja/i2008-10680-8.
  * [8] D. Ni, L. Wei, Z. Ren, Benford’s law and $\beta$-decay half-lives, Communications in Theoretical Physics 51 (4) (2009) 713–716.  doi:10.1088/0253-6102/51/4/25.
  * [9] X. J. Liu, X. P. Zhang, D. D. Ni, Z. Z. Ren, Benford’s law and cross-sections of A(n,$\alpha$)B reactions, The European Physical Journal A 47 (2011) 78.  doi:10.1140/epja/i2011-11078-3.
  * [10] H. Jiang, J. Shen, Y. Zhao, Benford's law in nuclear structure physics, Chinese Physics Letters 28 (3) (2011) 032101.  doi:10.1088/0256-307x/28/3/032101.
  * [11] L. Shao, B.-Q. Ma, Empirical mantissa distributions of pulsars, Astroparticle Physics 33 (4) (2010) 255–262.  doi:10.1016/j.astropartphys.2010.02.003.
  * [12] M. A. Snyder, J. H. Curry, A. M. Dougherty, Stochastic aspects of one-dimensional discrete dynamical systems: Benford’s law, Phys. Rev. E 64 (2001) 026222.  doi:10.1103/PhysRevE.64.026222.
  * [13] L. Shao, B.-Q. Ma, The significant digit law in statistical physics, Physica A: Statistical Mechanics and its Applications 389 (16) (2010) 3109–3116.  doi:10.1016/j.physa.2010.04.021.
  * [14] L. Shao, B.-Q. Ma, First-digit law in nonextensive statistics, Phys. Rev. E 82 (2010) 041110.  doi:10.1103/PhysRevE.82.041110.
  * [15] L. Shao, B.-Q. Ma, First digit distribution of hadron full width, Modern Physics Letters A 24 (40) (2009) 3275–3282.  doi:10.1142/S0217732309031223.
  * [16] HR.Varian, Letters to the editor, The American Statistician 26 (3) (1972) 62–65.  doi:10.1080/00031305.1972.10478934.
  * [17] L. Pericchi, D. Torres, Quick anomaly detection by the Newcomb–Benford law, with applications to electoral processes data from the USA, Puerto Rico and Venezuela, Statistical Science 26 (4) (2011) 502 – 516.  doi:10.1214/09-STS296.
  * [18] A. Cerioli, L. Barabesi, A. Cerasa, M. Menegatti, D. Perrotta, Newcomb–Benford law and the detection of frauds in international trade, Proceedings of the National Academy of Sciences 116 (1) (2019) 106–115.  doi:10.1073/pnas.1806617115.
  * [19] L. M. Leemis, B. W. Schmeiser, D. L. Evans, Survival distributions satisfying Benford’s law, American Statistician 54 (4) (2000) 236–241.  doi:10.1080/00031305.2000.10474554.
  * [20] R. J. C. Brown, Benford’s law and the screening of analytical data: the case of pollutant concentrations in ambient air, Analyst 130 (2005) 1280–1285.  doi:10.1039/B504462F.
  * [21] D. I. Cohen, An explanation of the first digit phenomenon, Journal of Combinatorial Theory, Series A 20 (3) (1976) 367–370.  doi:10.1016/0097-3165(76)90032-7.
  * [22] T. Jech, The logarithmic distribution of leading digits and finitely additive measures, Discrete Mathematics 108 (1) (1992) 53–57.  doi:10.1016/0012-365X(92)90659-4.
  * [23] D. S. Lemons, Thermodynamics of Benford’s first digit law, American Journal of Physics 87 (10) (2019) 787–790.  doi:10.1119/1.5116005.
  * [24] J. R. Iafrate, S. J. Miller, F. W. Strauch, Equipartitions and a distribution for numbers: A statistical model for Benford’s law, Phys. Rev. E 91 (2015) 062138\.  doi:10.1103/PhysRevE.91.062138.
  * [25] T. P. Hill, A statistical derivation of the significant-digit law, Statistical Science 10 (4) (1995) 354–363.  doi:10.1214/ss/1177009869.
  * [26] H.-A. Engel, C. Leuenberger, Benford’s law for exponential random variables, Statistics & Probability Letters 63 (4) (2003) 361–365.  doi:10.1016/S0167-7152(03)00101-9.
  * [27] M. Cong, B.-Q. Ma, A proof of first digit law from Laplace transform, Chinese Physics Letters 36 (7) (2019) 070201.  doi:10.1088/0256-307X/36/7/070201.
  * [28] M. Cong, C. Li, B.-Q. Ma, First digit law from Laplace transform, Physics Letters A 383 (16) (2019) 1836–1844.  doi:10.1016/j.physleta.2019.03.017.
  * [29] J. Cigler, Methods of summability and uniform distribution mod 1, Compositio Mathematica 16 (1964) 44–45, http://www.numdam.org/item?id=CM_1964__16__44_0.
  * [30] B. J. Flehinger, On the probability that a random integer has initial digit $A$, Amer. Math. Monthly 73 (1966) 1056–1061.  doi:10.2307/2314636.

Figure 2: The graphical abstract of this article, published in Fundamental
Research, 4 (2024) 842-845, by L. Wang & B.-Q. Ma,
doi:10.1016/j.fmre.2023.01.002.