# Geodesic Motions near an improved Schwarzschild black hole

Surajit Mandal

###### Abstract

In this paper, we studied the geodesics of timelike and null like particles
near an improved Schwarzschild black hole. The lapse function has been plotted
and was found that only one horizon is possible. The equation of motion and
effective potential of test particle have been calculated. This equation has
an importance in studying the radial free fall and also in studying the
stability of radial orbits (trajectories). The energy and angular momentum
have also been calculated to analysis the cicrular motion and stability of
circular orbits. Moreover, Innermost stable circular orbit radius has
determined. To get a deeper insight of the nature of these trajectories, we
have studied the timelike and null geodesics with the help of the dynamical
systems approach. This analysis help us to determine the stability as well as
the fixed point of phase space trajectories.

##  1 Introduction

Schwarzschild solution is necessarily a vacuum solution of Einstein’s field
equations, first developed by Schwarzschild in 1916 [1]. This vacuum solution
has importance in general relativity to describe a black hole which is
spherically symmetric in nature. The solution of improved Schwarzschild was
first derived by Bonanno and Reuter [21]. Nowadays it is very challenging
issues from theoretical prespective when quantum effect of geometry are taken
into consideration. It is believed that those quantum effects will make a key
role in the evaporation process of black hole(Planck size) as well as during
the late stages of the gravitational collapse. According to semiclassical
point of view, at a temperature a black hole emits Hawking radiation which is
inversely proportional to the mass of the black hole. In addition to this
process of Hawking radiation having continuous energy, a negative energy flux
through the horizon can produced. As a result, mass of the black hole will be
lowered and the temperature will be increased. It is yet unknown to us whether
this process will continues until the whole mass of the black hole is
converted to radiation or whether it will stops when the temperature of
radiative black hole reaches to the Planck temperature limit where the
semiclassical scenario are likely to break down.

One could speculate that the evaporation process comes to an end where the
Schwarzschild radius is very close to the Planck length where the
semiclassical results are no longer valid. Consequently the final state of the
Hawking evaporation could be some kind of “cold” remnant has a mass not too
far above the Planck mass. It is evident that the problem of the final state
might be addressed within the framework of quantum gravity theory. The matter
field is quantize only in semiclassical derivation of the Hawking temperature
and it helps us to consider the spacetime metric as a fixed classical
background. Moreover, quantum fluctuations of the metric has an important role
to inquire the black holes with a radius close to Planck length. The
perturbative quantization of Einstein gravity has a little importance here
since it provides a non-renormalizable theory. The fundamental theory of
quantum gravity such as string theory, loop quantum gravity, etc., do not
gives us a satisfactory answer yet [2]. In order to study quantum effects in
Schwarzschild spacetime, the idea of the Wilsonian renormalization group [3]
comes into the picture and this essentially leads to an improved Schwarzschild
black hole. Also, RG improvement of coupling effect on Schwarzschild black
hole using action improvement has been investigated by R. Moti and A. Shojai
[4, 5].

Geodesic motion is not only a problem of theoretical interest. This has far-
reaching consequences for the observations of these objects, as the motion of
light and gravitational waves on geodesic trajectories have consequences for
the observation of them, via shadows and mergers. In order to check the
viability of geodesic motion near black hole, it is important to study the
geodesic motion in four-dimensions as well as in higher-dimensions. For
n-dimensional rotating black holes the equations comes from geodesic study are
completely integrable, have studied by Don N. Page and others [7]. Geodesic
study in Schwarzschild–Ads space–time can ben found in [8]. The geodesic
structure for the Schwarzschild-Ads black hole has been studied by Norman Cruz
and others [9]. The motion of test particle around a black hole in a
braneworld can be found in [10]. Moreover, the motion of test particle in
different spacetimes was found in [11, 12, 13, 14, 15]. Stuchlik and Calvani
have studied on the effective potentials in case of radial null geodesics in
Kerr-de sitter and RN-de sitter spacetimes [16]. The study on bound orbits has
been found in [17, 18]. Hackmann and others [19] have been investigated the
behavior of timelike particles in four-dimensions. Recently, Wajiha Javed and
others [6] have studied on deflection angle and greybody bound of an improved
Schwarzschild black hole. Their analysis did not address the geodesic study
around this black hole. In this paper we tried to do this analysis using
effective Newtonian potential method. This paper also includes dynamical
approach to understand the phase-space trajectories near this black hole.

The paper is organized in six parts. In Section, 2, the line element and the
lapse function has been defined and the lapse function has been plotted to
understand the nature of event horizon. In Section 3, the equations for the
geodesics motion in 4-dimension around the black hole are found. In Section 4,
to understand the nature of orbits of test particles the method of effective
Newtonian orbit calculation has been used. The equation of motion of the test
particle of unit mass is estimated and the test particle moves under effective
potential with a given energy and angular momentum has found to depend on the
variuos parameter like free parameter, fixed parameter, mass of the black hole
and also on the radial coordinate. This equation helps us to study radial
motion of test particle and the corresponding stability of radial orbits. To
do analysis the circular motion of test particles and the stability of
circular orbits, we have used the method of transformation and determined the
energy and angular momentum of massive particles. The radius of the innermost
stable circular orbit (ISCO) of timelike (massive) particles is completely
depends on their angular momentum and fixed parameter for a vanishing free
parameter. To get a deeper insight of the nature of these trajectories, we
have studied the timelike and null geodesics with the help of the dynamical
systems approach in Section 5. The fixed points are found and the nature of
the phase space trajectories have been studied. Finally we make the summary
and conclusions in Section 6.

##  2  Preliminaries

The line element for an improved Schwarzschild black hole spacetime is defined
as,

$ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Omega_{2}^{2}$ (2.1)

where r is the radial coordinate and

$d\Omega_{2}^{2}=d\theta^{2}+Sin^{2}\theta d\phi^{2}$ (2.2)

Here, $d\Omega_{2}^{2}$ denotes the line element of unit two dimensional
sphere. For this static and the spherically symmetric spacetime, the value of
lapse function f(r) is defined as, [21, 22]

$f(r)=1-\frac{2Mr^{2}}{r^{3}+\tilde{\omega}(r+\psi M)}$ (2.3)

after simplifying the function (2.3) takes the form as [6]

$f(r)=1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}=\frac{\Delta}{r^{5}}$
(2.4)

here M is the geometric mass of black hole , r is the radial coordinate,
$\psi$ denotes free parameter and $\tilde{\omega}$ is fixed parameter. Fixed
parameter $\tilde{\omega}$ and free parameter $\psi$ can take the value [20]

$\tilde{\omega}=\frac{118}{15\pi}$

and

$\psi\geq 0$

For a given M, $\tilde{\omega}$ and $\psi$, the horizon function $\Delta$
depends only on the radial coordinate r. Pluggingg the value of f(r) in
equation (2.1), as a result we get the following

$ds^{2}=-\frac{\Delta}{r^{5}}dt^{2}+\frac{r^{5}}{\Delta}dr^{2}+r^{2}d\theta^{2}+r^{2}Sin^{2}\theta
d\phi^{2}$ (2.5)

The trajectories of the particle for this static solution are more complicated
than the Schwarzschild case. In order to make the analysis simpler, the
geodesics of this metrics are generally studied in terms of an effective
Newtonian potential. The behavior of the real singularity at r = 0, depends on
the parameter $\tilde{\omega}$ and $\psi$, which we assume in such a manner
that there won’t be any space-like naked singularity in the spacetime. A
comparison study of the effective potentials for both the Newtonian and the
Schwarzschild case has been discussed in [23]. We can analyze the positions of
the horizons of this black hole by making $g^{rr}=0$, i.e

$f(r)=\frac{\Delta}{r^{5}}=0$ (2.6)

where $\Delta$ =
$r^{5}-2Mr^{4}+2M\tilde{\omega}r^{2}-\tilde{\omega}^{2}r-2\tilde{\omega}^{2}\psi
M$

The positive real zeros of horizon function $\Delta$ gives the position of the
horizon of this black hole. Figure 1 illustrates the profile of an improved
Schwarzschild horizon function $\Delta$ with respect to the radial coordinate
r for three different values of free parameter $\psi$ with a fixed parameter
$\tilde{\omega}=\frac{118}{15\pi}$. Left figure is for free parameter
$\psi=0.1,0.3,0.5$ while the right figure is for varying free parameter
$\psi=1.2,1.35,1.5$. Both plot shows the existence of one event horizon of the
improved Schwarzschild black hole, indicating one positive real root of the
lapse function. To get the horizon we can consider only one constant value of
fixed parameter whereas we can take different values of $\psi$ as mentioned
above.

$\begin{array}[]{cccc}\begin{@subfigure}\includegraphics[width=216.81pt]{e01.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=216.81pt]{e02.eps}\end{@subfigure}\end{array}$

Figure 1:  In 1, the behavior of horizon function $\Delta$ with respect to r
by changing $\psi$=0.1,0.3 and 0.5. In 1, the behavior of f(r) with respect to
r by changing $\psi$=1.2,1.35 and 1.5. Here $\tilde{\omega}=\frac{118}{15\pi}$
and M=1.

##  3 Four-dimensional geodesics

Due to the spherical symmetry of an improved Schwarzschild black hole, we can
study the motion of the test particles on the equatorial plane, i.e at
$\theta$ = $\frac{\pi}{2}$. Therefoe, we have three geodesic equations:

$\frac{d^{2}t}{d\lambda^{2}}+\frac{B(r)}{A(r)}\frac{dt}{d\lambda}\frac{dr}{d\lambda}=0$
(3.1)

$\frac{d^{2}r}{d\lambda^{2}}+A(r)B(r)\Big{(}\frac{dt}{d\lambda}\Big{)}^{2}-\frac{B(r)}{A(r)}\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}+rA(r)\Big{(}\frac{d\phi}{d\lambda}\Big{)}^{2}=0$
(3.2)

$\frac{d^{2}\phi}{d\lambda^{2}}+\frac{2}{r}\frac{dr}{d\lambda}\frac{d\psi}{d\lambda}=0$
(3.3)

where $A(r)=-f(r)$ and
$B(r)$=$\frac{1}{r^{2}}\Big{[}-\frac{2\tilde{\omega}^{2}}{r^{3}}-\frac{5\tilde{\omega}^{2}\psi
M}{r^{4}}-M+\frac{3M\tilde{\omega}}{r^{2}}\Big{]}$
To solve this geodesic equations we will consider two approach in the next
section. One is effective Newtonian potential approach as studied in [24] and
another one is dynamical system analysis approach [25, 26].

##  4 Effective Newtonian potential approach

The lagrangian for particle motion for improved Schwarzschild black hole is
given by

$2\mathcal{L}=-\Big{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Big{)}\dot{t}^{2}+\frac{\dot{r}^{2}}{1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}}+r^{2}(\dot{\theta}^{2}+Sin^{2}\theta\dot{\phi}^{2})$
(4.1)

here overdot denotes the differentiation with respect to $\lambda$ (affine
parameter). Here lagrangian $\mathcal{L}$ is not depend explicitly on the
coordinates t and $\phi$. Hence, due to this two cyclic coordinates we get two
conserved quantities, exemplary the energy (E) and the momentum (h) conjugate
to $\phi$,
Energy,

$E=g_{tt}\frac{dt}{d\lambda}=-f(r)\frac{dt}{d\lambda}=A(r)\frac{dt}{d\lambda}$
(4.2)

Momentum,

$2h=\frac{\partial{\mathcal{L}}}{\partial{\dot{\phi}}}=2r^{2}\dot{\phi}=constant$
(4.3)

here h indicates the total angular momentum of the particles. Recalling the
normalization condition

$g_{\mu\nu}\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda}=-\epsilon$ (4.4)

for timelike geodesics $\epsilon$ = 1 and for null geodesics $\epsilon$ = 0
[27]. On the equatorial plane we have,

$\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}=E^{2}-\frac{\Delta}{r^{5}}\Big{(}\frac{h^{2}}{r^{2}}+\epsilon\Big{)}=E^{2}+A(r)\Big{(}\frac{h^{2}}{r^{2}}+\epsilon\Big{)}$
(4.5)

Thus equations (4.2), (4.3) and (4.5) are required to describe the dynamics of
particle trajectories at the equatorial plane of an improved Schwarzschild
black hole.
Rewriting (4.5),

$\frac{1}{2}\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}=E_{eff}-V_{eff}$ (4.6)

where

$E_{eff}=\frac{E^{2}}{2}$ (4.7)

and

$V_{eff}=\frac{\Delta}{2r^{5}}\Big{(}\frac{h^{2}}{r^{2}}+\epsilon\Big{)}=\frac{1}{2}\Big{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Big{)}\Big{(}\frac{h^{2}}{r^{2}}+\epsilon\Big{)}$
(4.8)

Therefore, Eq. (4.6) represents the equation of motion of a particle having
unit mass and an effective energy $E_{eff}$ moving in a one-dimensional(1D)
effective potential $V_{eff}(r)$. Though E denotes the conserved energy of the
particle per unit mass, $V_{eff}(r)$ for the radial coordinate r responds to
$E_{eff}$ . This effective potential will be vanishes at the zeros of the
lapse function. The physically acceptable regions are given by the values of
those r for which $E_{eff}$ $>$ $V_{eff}(r)$ . Thus Eq. (4.6) represents the
energy equation for the radial coordinate r. This equation is required to
study the radial free fall and the stability of particle trajectories.
By plugging A(r), equation (4.5) takes the form:

$\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}=E^{2}-\epsilon+\frac{2M\epsilon}{r}-\frac{h^{2}}{r^{2}}+\frac{2Mh^{2}}{r^{3}}-\frac{2M\tilde{\omega}\epsilon}{r^{3}}+\frac{\tilde{\omega}^{2}\epsilon}{r^{4}}+\frac{2\tilde{\omega}^{2}\psi
M\epsilon}{r^{5}}-\frac{2M\tilde{\omega}h^{2}}{r^{5}}+\frac{\tilde{\omega}^{2}h^{2}}{r^{6}}+\frac{2\tilde{\omega}^{2}\psi
Mh^{2}}{r^{7}}$ (4.9)

The geometry of the geodesics in the equatorial plane $\theta$ =
$\frac{\pi}{2}$ can be determined by above Eq. (4.9). To determine the shape
of the trajectories we use (4.3) to express $\frac{dr}{d\lambda}$ as

$\frac{dr}{d\lambda}=\frac{dr}{d\phi}\frac{d\phi}{d\lambda}=\frac{h}{r^{2}}\frac{dr}{d\phi}$
(4.10)

Now, introducing $u=\frac{1}{r}$, we arrive at

$\Big{(}\frac{du}{d\lambda}\Big{)}^{2}=u^{4}\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}$
(4.11)

From Eq. (4.9)

$\displaystyle\Big{(}\frac{du}{d\lambda}\Big{)}^{2}$ $\displaystyle=$
$\displaystyle 2\tilde{\omega}^{2}\psi
Mh^{2}u^{11}+\tilde{\omega}^{2}h^{2}u^{10}+(2\tilde{\omega}^{2}\psi
M\epsilon-2M\tilde{\omega}h^{2})u^{9}+\tilde{\omega}^{2}\epsilon
u^{8}+(2Mh^{2}-2M\tilde{\omega}\epsilon)u^{7}$ (4.12) $\displaystyle-$
$\displaystyle h^{2}u^{6}+2M\epsilon u^{5}+(E^{2}-\epsilon)u^{4}$

Now, using equations (4.10), (4.11) and (4.12), we obtained the expression as
following :

$\displaystyle\Big{(}\frac{du}{d\phi}\Big{)}^{2}$ $\displaystyle=$
$\displaystyle 2\tilde{\omega}^{2}\psi
Mu^{7}+\tilde{\omega}^{2}u^{6}+\Big{(}\frac{2\tilde{\omega}^{2}\psi
M\epsilon}{h^{2}}-2M\tilde{\omega}\Big{)}u^{5}+\frac{\tilde{\omega}^{2}\epsilon}{h^{2}}u^{4}+\Big{(}2M-\frac{2M\tilde{\omega}\epsilon}{h^{2}}\Big{)}u^{3}$
(4.13) $\displaystyle-$ $\displaystyle
u^{2}+\frac{2M\epsilon}{h^{2}}u+\frac{E^{2}-\epsilon}{h^{2}}=S(u)$

This equation describes the trajectories of the test particles near an
improved Schwarzschild black hole. In the physical point of view, the radial
motion and the circular motion of the particles are two important aspects for
studying the particles trajectories. Therefore, in next section we will
discuss the particle’s trajectories in these two cases.

###  4.1 Radial motion

We need to study the radial geodesics for illustrating the essential features
of the space-time, for which angular momentum becomes zero. Also, for radial
motion, $\phi$ = constant. Equation (4.13) fails to provide us the information
about the radial trajectories as most of the terms blow up for h = 0.
Moreover, to study the radial trajectories of the particles we can consider
equation (4.5), which becomes

$\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}=E^{2}+A(r)\epsilon$ (4.14)

####  4.1.1 Motion of massive particles

For massive particles $\epsilon$ = 1. Eq. (4.14) gives

$\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}=E^{2}+A(r)=E^{2}-1+\frac{\tilde{\omega}^{2}}{r^{4}}+\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}+\frac{2M}{r}-\frac{2M\tilde{\omega}}{r^{3}}$ (4.15)

By differentiating (4.15) with respect to $\lambda$ and then dividing it by
2$\frac{dr}{d\lambda}$, we obtain

$\frac{d^{2}r}{d\lambda^{2}}=-\Big{(}\frac{2\tilde{\omega}^{2}}{r^{5}}+\frac{5\tilde{\omega}^{2}\psi
M}{r^{6}}+\frac{M}{r^{2}}\Big{)}+\frac{3M\tilde{\omega}}{r^{4}}$ (4.16)

Since the trajectories of massive particles are characterized by timelike
geodesics, so here we can assume the affine parameter proper time $\tau$
instead of $\lambda$ along the path. Hence the motion of massive particles can
studied by the above equation. The condition for attractive force per unit
mass is

$\frac{2\tilde{\omega}^{2}}{r^{5}}+\frac{5\tilde{\omega}^{2}\psi
M}{r^{6}}+\frac{M}{r^{2}}>\frac{3M\tilde{\omega}}{r^{4}}$ (4.17)

this condition is important to get the bound states of massive particles. The
particle can acquire kinetic energy during the gravitational interaction when
it goes through the gravitational field of any black hole. Then we can
estimate the change in the potential energy in that gravitational field of the
particle by taking into account particle’s rest position such that the radial
coordinate r changes to R. We can write Eq. (4.15) in terms of the affine
parameter $\tau$, that gives

$\Big{(}\frac{dr}{d\tau}\Big{)}^{2}=(E^{2}-1)+\frac{\tilde{\omega}^{2}}{r^{4}}+\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}+\frac{2M}{r}-\frac{2M\tilde{\omega}}{r^{3}}$ (4.18)

For r = R, $\frac{dr}{d\tau}$ = 0 and we can arrive at the equation for the
kinetic energy per unit mass of the particle in terms of the change in its
gravitational potential energy as follows:

$\frac{1}{2}\Big{(}\frac{dr}{d\tau}\Big{)}^{2}=\frac{\tilde{\omega}^{2}}{2}\Big{(}\frac{1}{r^{4}}-\frac{1}{R^{4}}\Big{)}+\tilde{\omega}^{2}\psi
M\Big{(}\frac{1}{r^{5}}-\frac{1}{R^{5}}\Big{)}+M\Big{(}\frac{1}{r}-\frac{1}{R}\Big{)}-M\tilde{\omega}\Big{(}\frac{1}{r^{3}}-\frac{1}{R^{3}}\Big{)}$
(4.19)

Now, we are interested to represent the trajectories of the particle in the
(t,r)-plane. With the help of Eqs. (4.2) and (4.15),

$\Big{(}\frac{dr}{dt}\Big{)}^{2}=\Big{(}E^{2}+A(r)\Big{)}\frac{A^{2}(r)}{E^{2}}$
(4.20)

It is evident that $\frac{dr}{dt}$ = 0 for all points when $E^{2}$ = f(r) but
f(r) $\neq$ 0 and E $\neq$ 0\. From Eq. (4.15) we get the effective potential
$V_{eff}$ as follows :

$V_{eff}=\frac{1}{2}\Bigg{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Bigg{)}$ (4.21)

$\begin{array}[]{cccc}\begin{@subfigure}\includegraphics[width=151.76964pt]{e1.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e1a.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e1b.eps}\end{@subfigure}\end{array}$

Figure 2:  In 2 and 2, The behavior of $V_{eff}$ with respect to r by changing
$\tilde{\omega}$ for a fixed $\psi$ = 1. In 2, The behavior of $V_{eff}$ with
respect to r by changing $\psi$ for a fixed $\tilde{\omega}$ = 1. Here, $M=1$.

Figure 2 depicts the behavior of $V_{eff}$ with respect to the radius r of the
orbits for five different values of $\tilde{\omega}$ and $\psi$ of the black
hole. From figure it is clear to us that, $V_{eff}$ is an exponentially
increasing function of r but remains negative valued.
Now, we can write from eq. (4.20)

$\frac{dr}{dt}=\frac{1}{E}\Bigg{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Bigg{)}\Bigg{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Bigg{)}^{\frac{1}{2}}$
(4.22)

The above equation won’t be blow up for $r\rightarrow\infty$. Hence, this
equation will only be a function of total energy E and therefore if one can
observe the behavior of the particle at infinite distance from the source then
obviously it will behaves like a free particle. However, $\frac{dr}{dt}$ and E
will be finite for a measurable distance from souce. So on (t,r)-plane, we can
utilize (4.22) in principle to get the correct trajectories.

####  4.1.2 Motion of photons

For photons $\epsilon$ = 0. Eq. (4.14) becomes

$\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}=E^{2}$ (4.23)

where E is defined in eq. (4.2). Also defining
$\frac{dr}{dt}=\frac{dr}{d\lambda}\frac{d\lambda}{dt}$ we get,

$\frac{dr}{dt}=\pm
A(r)=\mp\Bigg{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Bigg{)}$ (4.24)

The solution of this equation will be complicated in nature and therefore it
is difficult to describe the trajectories of the photons using this method. If
we can consider a dynamical system by taking the variables
$\frac{dr}{d\lambda}$ and $\frac{dt}{d\lambda}$, then Eq. (4.23) will be able
to gives us the information about the particle trajectories. We will show such
an analysis later.

###  4.2 Circular motion

It is clear from eq. (4.13) that at equilibrium position of circular orbits
$u=constant$, so $S(u)=0$ and $S^{\prime}(u)$ = 0. From Eq. (4.13) we have,

$\displaystyle S(u)$ $\displaystyle=$ $\displaystyle 2\tilde{\omega}^{2}\psi
Mu^{7}+\tilde{\omega}^{2}u^{6}+\Big{(}\frac{2\tilde{\omega}^{2}\psi
M\epsilon}{h^{2}}-2M\tilde{\omega}\Big{)}u^{5}+\frac{\tilde{\omega}^{2}\epsilon}{h^{2}}u^{4}+\Big{(}2M-\frac{2M\tilde{\omega}\epsilon}{h^{2}}\Big{)}u^{3}$
(4.25) $\displaystyle-$ $\displaystyle
u^{2}+\frac{2M\epsilon}{h^{2}}u+\frac{E^{2}-\epsilon}{h^{2}}$

and

$\displaystyle S^{\prime}(u)$ $\displaystyle=$ $\displaystyle
14\tilde{\omega}^{2}\psi
Mu^{6}+6\tilde{\omega}^{2}u^{5}+5\Big{(}\frac{2\tilde{\omega}^{2}\psi
M\epsilon}{h^{2}}-2M\tilde{\omega}\Big{)}u^{4}+\frac{4\tilde{\omega}^{2}\epsilon}{h^{2}}u^{3}+3\Big{(}2M-\frac{2M\tilde{\omega}\epsilon}{h^{2}}\Big{)}u^{2}$
(4.26) $\displaystyle-$ $\displaystyle 2u+\frac{2M\epsilon}{h^{2}}$

After plugging the condition $S^{\prime}(u)=0$ and $S(u)=0$ in (4.26), we get
an expression of angular momentum (h) and energy (E) of the particle obeying
circular motion, respectively, as follows:

$\displaystyle h^{2}=\frac{\epsilon(5\tilde{\omega}^{2}\psi
Mu^{4}+2\tilde{\omega}^{2}u^{3}-3M\tilde{\omega}u^{2}+M)}{u(1-7\tilde{\omega}^{2}\psi
Mu^{5}-3\tilde{\omega}^{2}u^{4}+5M\tilde{\omega}u^{3}-3Mu)}$ (4.27)

and

$\displaystyle E^{2}=\frac{2\epsilon(1-2\tilde{\omega}^{2}\psi
Mu^{5}-\tilde{\omega}^{2}u^{4}+2M\tilde{\omega}u^{3}-2Mu)}{(1-7\tilde{\omega}^{2}\psi
Mu^{5}-3\tilde{\omega}^{2}u^{4}+5M\tilde{\omega}u^{3}-3Mu)}$ (4.28)

####  4.2.1 Motion of massive particles

For massive particles Eqs. (4.27) and (4.28) takes the form as,

$\displaystyle h^{2}=\frac{(5\tilde{\omega}^{2}\psi
Mu^{4}+2\tilde{\omega}^{2}u^{3}-3M\tilde{\omega}u^{2}+M)}{u(1-7\tilde{\omega}^{2}\psi
Mu^{5}-3\tilde{\omega}^{2}u^{4}+5M\tilde{\omega}u^{3}-3Mu)}$ (4.29)

and

$\displaystyle E^{2}=\frac{2(1-2\tilde{\omega}^{2}\psi
Mu^{5}-\tilde{\omega}^{2}u^{4}+2M\tilde{\omega}u^{3}-2Mu)}{(1-7\tilde{\omega}^{2}\psi
Mu^{5}-3\tilde{\omega}^{2}u^{4}+5M\tilde{\omega}u^{3}-3Mu)}$ (4.30)

$\begin{array}[]{cccc}\begin{@subfigure}\includegraphics[width=151.76964pt]{e2.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e2a.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e2b.eps}\end{@subfigure}\end{array}$

Figure 3:  In 3 and 3, The behavior of $E^{2}$ with respect to $\frac{r}{M}$
by changing $\tilde{\omega}$ for a fixed $\psi$ = 1. In 3, The behavior of
$E^{2}$ with respect to $\frac{r}{M}$ by changing $\psi$ for a fixed
$\tilde{\omega}$ = 1. Here, $M=1$.

The radius of the circular (bound) orbits can be calculated from eqs. (4.29)
and (4.30). Moreover, the energy (E) of the particle will be a constant of
motion for radius r (=$\frac{1}{u}$) of a circular orbit. However, when the
radius r varies E also vary. The behavior of $E^{2}$ vs specific radius
$\frac{r}{M}$ for massive particles executing circular motion around an
improved Schwarzschild black hole, for five different values of
$\frac{\tilde{\omega}}{M}$ and $\frac{\psi}{M}$, is plotted in figure 3. Also,
the variation of angular momentum $h^{2}$ with respect to specific radius
$\frac{r}{M}$ for five different values of $\frac{\tilde{\omega}}{M}$ and
$\frac{\psi}{M}$, is depicted in figure 4.

$\begin{array}[]{cccc}\begin{@subfigure}\includegraphics[width=151.76964pt]{e3.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e3a.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e3b.eps}\end{@subfigure}\end{array}$

Figure 4:  In 4 and 4, The behavior of $h^{2}$ with respect to $\frac{r}{M}$
by changing $\tilde{\omega}$ for a fixed $\psi$ = 1. In 4, The behavior of
$h^{2}$ with respect to $\frac{r}{M}$ by changing $\psi$ for a fixed
$\tilde{\omega}$ = 1. Here, $M=1$.

In figures 3, 3, 4 and 4 the energy E and angular momentum h is a decreasing
function of $\frac{r}{M}$ for different values of $\frac{\tilde{\omega}}{M}$
but constant free parameter $\psi$ while in figure 3 and 4, both E and h is a
decreasing function of $\frac{r}{M}$ for different values of $\frac{\psi}{M}$
but constant fixed parameter $\tilde{\omega}$. In figures 3 and 4, energy and
angular momentum increases for large values of $\frac{\tilde{\omega}}{M}$,
indicating that for test particles to maintain in a stable circular orbit, a
more energy E and angular momentum h for a larger values of
$\frac{\tilde{\omega}}{M}$ is required. In figures 3 and 4, energy and angular
momentum increases for large values of $\frac{\psi}{M}$ respectively,
indicating that for test particles to maintain in a stable circular orbit, a
more energy E and angular momentum h for a larger values of $\frac{\psi}{M}$
is needed. However, in 3 and 4, for particles to maintain in a stable circular
orbit, a more energy and angular momentum for a smaller values of
$\frac{\tilde{\omega}}{M}$ is requisite. For increasing values of
$\frac{r}{M}$, both E and h of all the plots are a increasing function.

####  4.2.2 Motion of photons

It is obvious that we can’t calculate E and h of photons from eqs. (4.27) and
(4.28). From (4.13),

$\Big{(}\frac{du}{d\phi}\Big{)}=2\tilde{\omega}^{2}\psi
Mu^{7}+\tilde{\omega}^{2}u^{6}-2M\tilde{\omega}u^{5}+2Mu^{3}-u^{2}+\frac{E}{h^{2}}=T(u)$
(4.31)

Hence,

$T^{\prime}(u)=14\tilde{\omega}^{2}\psi
Mu^{6}+6\tilde{\omega}^{2}u^{5}-10M\tilde{\omega}u^{4}+6Mu^{2}-2u$ (4.32)

For $T^{\prime}(u)=0$, we get

$7\tilde{\omega}^{2}\psi
Mu^{5}+3\tilde{\omega}^{2}u^{4}-5M\tilde{\omega}u^{3}+3Mu-1=0$ (4.33)

It is a ploynomial equation of degree five, so it has five roots.
For $T(u)=0$, we obtain

$\frac{E}{h^{2}}=-2\tilde{\omega}^{2}\psi
Mu^{7}-\tilde{\omega}^{2}u^{6}+2M\tilde{\omega}u^{5}-2Mu^{3}+u^{2}$ (4.34)

After putting one of the solution of (4.33) in (4.34), We can arrive at the
relation between energy E and angular momentum h of photons in terms of the
geometric mass M of the improved Schwarzschild black hole, the free parameter
$\psi$ and the fixed parameter $\tilde{\omega}$.

###  4.3 Stability of orbits for massive particles

We get freom (4.8),

$V_{eff}=\frac{1}{2}\Bigg{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Bigg{)}\Big{(}1+\frac{h^{2}}{r^{2}}\Big{)}$
(4.35)

The graph for effective potential vs $\frac{r}{M}$ for different values of
$\frac{\tilde{\omega}}{M}$ and $\frac{\psi}{M}$, is shown in figure 5. The
behavior of $V_{eff}$ is not similar to Newtonian potentials [23]. Here, we
observe that for all values of r, potential remains always negative which
describes a bound system. From figure we see that effective potential graph
has no maxima and minima. This clearly tell us that for each curve there exist
only point of inflextion which indicates the existence of only marginally
stable circular orbit ( known as innermost stable circular orbit, abbreviated
as ISCO). Also the nature of this effective potential is different from the
Schwarzschild black hole case [28]. If we can make $\tilde{\omega}=0$ and
$\psi=0$ in (4.35), it returns the well known Schwarzschild effective
potential. The detailed study on effective potential and the probable orbits
in Schwarzschild space-time has been discussed on [29, 30].

$\begin{array}[]{cccc}\begin{@subfigure}\includegraphics[width=151.76964pt]{e4.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e4a.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=151.76964pt]{e4b.eps}\end{@subfigure}\end{array}$

Figure 5:  In 5 and 5, The behavior of $V_{eff}$ with respect to r by changing
$\tilde{\omega}$ for a fixed $\psi$ = 1. In 5, The behavior of $V_{eff}$ with
respect to r by changing $\psi$ for a fixed $\tilde{\omega}$ = 1. Here, $M=1$.

Stable circular orbits and unstable circular orbits gives the value of radius
which represents the position of local minima and local maxima of the
potential respectively. For stable circular orbits, we have $S(u)=0$,
$S^{\prime}=0$ and $S^{\prime\prime}=0$. We get from eq. (4.26)

$S^{\prime\prime}(u)=84\tilde{\omega}^{2}\psi
Mu^{5}+30\tilde{\omega}^{2}u^{4}+40M\Big{(}\frac{\tilde{\omega}^{2}\psi}{h^{2}}-\tilde{\omega}\Big{)}u^{3}+\frac{12\tilde{\omega}^{2}}{h^{2}}u^{2}+12M\Big{(}1-\frac{\tilde{\omega}}{h^{2}}\Big{)}u-2$
(4.36)

For $S^{\prime\prime}(u)=0$,

$84\tilde{\omega}^{2}\psi
Mh^{2}u^{5}+30\tilde{\omega}^{2}h^{2}u^{4}+40M\Big{(}\tilde{\omega}^{2}\psi-\tilde{\omega}h^{2}\Big{)}u^{3}+12\tilde{\omega}^{2}u^{2}+12M\Big{(}h^{2}-\tilde{\omega}\Big{)}u-2h^{2}=0$
(4.37)

Therefore, for massive particles, to evaluate the innermost stable circular
orbit (ISCO) near an improved Schwarzschild black hole we need to solve the
polynomial of degree 5 in u. We can calculate the radius of ISCO of massive
particles for improved Schwarzschild black hole from eq. (4.37), provided
$\tilde{\omega}\psi\neq h^{2}$ and $h^{2}\neq\tilde{\omega}$. An observer can
see this track if this minimum radius will be bigger than the radius of event
horizon.
If $\psi=0$, $\tilde{\omega}\psi=h^{2}$ and $h^{2}=\tilde{\omega}$ in eq.
(4.37), then ISCO of improved Schwarzschild black hole is given by,

$r^{2}_{ISCO}=\frac{60\tilde{\omega}^{2}h^{2}}{-12\tilde{\omega}\pm\sqrt{144\tilde{\omega}^{2}+240\tilde{\omega}^{2}h^{4}}}$
(4.38)

###  4.4 Stability of orbits for photons

To discuss the stability of photon orbits, we get from (4.5)

$\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}-A(r)\frac{h^{2}}{r^{2}}=E^{2}$ (4.39)

it takes the following form :

$\frac{1}{h^{2}}\Big{(}\frac{dr}{d\lambda}\Big{)}^{2}+\bar{V}_{eff}=\bar{E}_{eff}$
(4.40)

where
$V_{eff}=\frac{1}{r^{2}}\Bigg{(}1-\frac{\tilde{\omega}^{2}}{r^{4}}-\frac{2\tilde{\omega}^{2}\psi
M}{r^{5}}-\frac{2M}{r}+\frac{2M\tilde{\omega}}{r^{3}}\Bigg{)}$ and
$E_{eff}=\frac{E^{2}}{L^{2}}$
If the photon follows circular paths, then we estimate the minimum radius
($r_{mc}$) for the stable circular orbits of photon by considering the
condition $T^{\prime\prime}(u)=0$, which gives

$84\tilde{\omega}^{2}\psi
Mu^{5}+30\tilde{\omega}^{2}u^{4}-40M\tilde{\omega}u^{3}+12Mu-2=0$ (4.41)

It is a ploynomial equation of degree five. Radius $r_{mc}$ of stable circular
orbits can be estimated from (4.41) if $\psi\neq 0$ and $\tilde{\omega}\neq
o$. When free parameter $\psi$ of improved Schwarzschild black hole vanishes
then (4.41) reduces to a polynomial of degree four. Therefore $r_{mc}$ can be
calculated as,

$\displaystyle r_{mc}$ $\displaystyle>$
$\displaystyle\frac{3M}{2}+\frac{1}{2}\sqrt{9M^{2}+\frac{1202^{\frac{1}{3}}\tilde{\omega}(2M^{2}-\tilde{\omega})}{\tilde{\omega}^{2}\mathcal{P}}+\frac{\mathcal{P}}{62^{\frac{1}{3}}\tilde{\omega}^{2}}}$
(4.42) $\displaystyle+$
$\displaystyle\frac{1}{2}\sqrt{18M^{2}-\frac{1202^{\frac{1}{3}}\tilde{\omega}(2M^{2}-\tilde{\omega})}{\tilde{\omega}^{2}\mathcal{P}}-\frac{\mathcal{P}}{62^{\frac{1}{3}}\tilde{\omega}^{2}}+\frac{216M^{3}-160M\tilde{\omega}}{\mathcal{D}}}$
(4.43)

where,
$\mathcal{P}=\Big{(}-30240M^{2}\tilde{\omega}^{2}-\sqrt{914457600M^{4}\tilde{\omega}^{4}-4(1440M^{2}\tilde{\omega}-720\tilde{\omega}^{2})^{3}}\Big{)}^{\frac{1}{3}}$
and
$\mathcal{D}=4\sqrt{9M^{2}+\frac{1202^{\frac{1}{3}}\tilde{\omega}(2M^{2}-\tilde{\omega})}{\tilde{\omega}^{2}\mathcal{P}}+\frac{\mathcal{P}}{62^{\frac{1}{3}}\tilde{\omega}^{2}}}$

##  5  Dynamical systems analysis

In this section we will define geodesics equations in terms of dynamical
variables and thereafter study its phase space. A detailed study on particle
trajectories using dynamical system approach can be found in [31, 32]. Let us
define three dynamical variable as follows :

$\frac{dt}{d\lambda}=X$ (5.1) $\frac{dr}{d\lambda}=Y$ (5.2)
$\frac{d\phi}{d\lambda}=Z$ (5.3)

Now, geodesic eqs. (3.1) - (3.3) in terms of these three new variables,

$\frac{dX}{d\lambda}+\frac{B(r)}{A(r)}XY=0$ (5.4)
$\frac{dY}{d\lambda}+A(r)B(r)X^{2}-\frac{B(r)}{A(r)}Y^{2}+rA(r)W^{2}=0$ (5.5)

and

$\frac{dZ}{d\lambda}+\frac{1}{r}YZ=0$ (5.6)

X, Y and Z are ralated by the following equation

$A(r)X^{2}-\frac{1}{A(r)}Y^{2}+r^{2}Z^{2}=-\epsilon$ (5.7)

As W can be calculated from (5.7), so studying of eqs. (5.4), (5.5) and (5.6)
is not necessary. Using Eqs. (5.7) and (5.5), we get

$\frac{dY}{d\lambda}+\frac{A(r)}{r}(rB(r)-A(r))X^{2}+\frac{1}{rA(r)}(A(r)-rB(r))Y^{2}-\frac{\epsilon
A(r)}{r}=0$ (5.8)

###  5.1 Real non-linear dynamical system

$\frac{dX}{d\lambda}=M(X,Y,r)$ (5.9) $\frac{dY}{d\lambda}=N(X,Y,r)$ (5.10)
$\frac{dr}{d\lambda}=Y$ (5.11)

where $M(X,Y,r)=-\frac{B(r)}{A(r)XY}$
$N(X,Y,r)=\frac{\epsilon
A(r)}{r}-\frac{A(r)}{r}(rB(r)-A(r))X^{2}-\frac{1}{rA(r)}(A(r)-rB(r))Y^{2}$

###  5.2 Fixed point determination

To get the fixed point, we have to solve eqs. (5.9), (5.10) and (5.11), for
which Y, M and N becomes zero. Hence, equilibrium positions of particle
trajectories can be obtained from these solutions. For phase space
trajectories we consider the fixed point is $(X_{0},Y_{0}=0,r_{0})$. For a
black hole with a given $\psi$, $\tilde{\omega}$ and mass M, the fixed point
in the (X-Y) plane for null geodesics is (0, 0) whereas timelike geodesics
produces definite fixed points. For timelike geodesics,

$X_{0}=\pm\sqrt{\frac{1}{r_{0}B(r_{0}-A(r_{0}))}}$ (5.12)

###  5.3 Phase space trajectories analysis

The phase evolution of the dynamical system in the (X, Y ) phase plane is
obtained from the following equation :

$\frac{dY}{dX}=\frac{N}{M}$ (5.13)

here $M\neq 0$ and in case of null geodesics it gives the information of phase
space evolution except at the point (0, 0).
Plugging the values of M and N in (5.13) and after integrating we obtain

$g_{1}X^{2}-g_{2}Y^{2}-g_{3}=0$ (5.14)

where, $g_{1}=\frac{A(r)C(r)}{3r}$,
$g_{2}=\frac{C(r)}{rA(r)}+\frac{B(r)}{2A(r)}$,
$g_{3}=\frac{\epsilon}{A}(r){r}$ and $C(r)=rB(r)-A(r)$. For null geodesics
$g_{3}=0$.

$\begin{array}[]{cccc}\begin{@subfigure}\includegraphics[width=195.12767pt]{e5.eps}\end{@subfigure}\begin{@subfigure}\includegraphics[width=195.12767pt]{e6.eps}\end{@subfigure}\end{array}$

Figure 6:  In 6, the plot for timelike geodesics for three sample equations
$10X^{2}+20Y^{2}-700=0$, $20X^{2}+30Y^{2}-700=0$, $25X^{2}+35Y^{2}-500=0$. In
6, the plot for null geodesics for three sample equations
$-10X^{2}+20Y^{2}=0$, $-20X^{2}+30Y^{2}=0$, $-25X^{2}+35Y^{2}=0$.

Plot 6 depicts the possible trajectories of timelike particles for different
values of $g_{1}$, $g_{2}$ and $g_{3}$ on (X,Y)-plane. We see that the orbits
are elliptical (periodic bpund) in nature. In plot 6, the orbits of null like
particles for different values of $g_{1}$ and $g_{2}$ has been shown. The
phase trajectories of null particles are necessarily straight lines passing
through the origin (0, 0), with their slopes get changes with the varying
$g_{1}$ and $g_{2}$.

##  6 Conclusions and remarks

In this paper, We have investigated the timelike and null geodesics near an
improved Schwarzschild black hole. We have determined the position of event
horizon by plotting the lapse function of this black hole and the plot tells
us that only one event horizon is physically acceptable. To investigate the
circular motion and stability of circular orbits, we took inverse radial
coordinate r and energy E and angular momentum h of massive particles has been
calculated. The nature of the particle trajectories depend upon the energy,
angular momentum, free parameter $\psi$ and fixed parameter $\tilde{\omega}$.
We have noticed that the enargy and angular momentum of massive particles is a
decreasing function of specific radius for different values of
$\frac{\tilde{\omega}}{M}$ and $\frac{\psi}{M}$ but for increasing values of
specific radius both energy and angular momentum are increasing function. We
have seen that the radius of the innermost stable circular orbit (ISCO) of
massive particles is completely defined in terms of their angular momentum h
and fixed parameter $\tilde{\omega}$, free parameter $\psi$ and mass M of the
black hole. Since it is impossible to obtain the possible trajectories using
Newtonian potential method, so simultaneously we have studied dynamical
systems. Also, We have noticed that the geodesics of timelike (massive)
particles near an improved Schwaezschild black hole possesses only periodic
bound orbits and for these orbits there are definite fixed point. Moreover,
the null geodesics produces a unique fixed point and these trajectories are
terminating orbits in nature.

Data availibility
Data sharing is not applicable to this article, as no data sets were generated
or analyzed during the current study.

Acknowledgements
The author is thankful to Dr. Sudhaker Upadhyay, department of physics, K. L.
S. College, Nawada, Bihar 805110, India, for various suggestions which
developed the presentation of the paper.

## References

  * [1] K. Schwarzschild, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.) 1916, 189 (1916).
  * [2] C.Rovelli, Contribution to the Special Issue 2000 of the Journal of Mathematical Physics, e-Print Archive: hep-th/9910131.
  * [3] K.Wilson, J.Kogut, Phys. Rep. 12C, 75, (1974).
  * [4] R. Moti and A. Shojai, International Journal of Modern Physics AVol. 35, No. 04, 2050016 (2020).
  * [5] R. Moti and A. Shojai, Phys. Lett. B 793, 313 (2019).
  * [6] Wajiha Javed, Muhammad Aqib and Ali Övgün, New Astronomy Volume 96, October 2022, 101827.
  * [7] Page D N, Kubiznak D, Vasudevan M and Krtous P 2007 Phys. Rev. Lett. 98 061102.
  * [8] Eva Hackmann, Claus Lämmerzahl, Phys. Rev. D 78, 024035 (2008).
  * [9] Cruz N, Olivares M and Villanueva J R 2005 Class. Quant. Grav. 22 1167.
  * [10] Abdujabbarov A and Ahmedov B 2010 Phys. Rev. D 81 044022.
  * [11] D. V. Galtsov and E. Masar, Classical Quantum Gravity 6, 1313 (1989).
  * [12] S. Chakraborty and L. Biswas, Classical Quantum Gravity 13, 2153 (1996).
  * [13] N. Ozdemir, Classical Quantum Gravity 20, 4409 (2003).
  * [14] F. Ozdemir, N. Ozdemir, and B. T. Kaynak, Int. J. Mod. Phys. A 19, 1549 (2004).
  * [15] A. N. Aliev and D. V. Galtsov, Sov. Astron. Lett. 14, 48 (1988).
  * [16] Stuchlik Z and Calvani M 1991 Gen. Relativ. Grav. 23 507.
  * [17] Kraniotis G 2004 Class. Quant. Grav. 21 4743.
  * [18] Fujita R and Hikada W 2009 Class. Quant. Grav. 26 135002.
  * [19] Hackmann E, Kagramanova V, Kunz J and L¨ammerzahl C 2008 Phys. Rev. D 78 124018.
  * [20] Alfio Bonanno, Martin Reuter, Renormalization group improved black hole space–times, Phys. Rev. D 62 (2000) 043008.
  * [21] A. Bonanno and M. Reuter, Phys. Rev. D 62, 043008 (2000).
  * [22] A. Rincon and G. Panotopoulos, Phys. Dark Univ.30, 100639 (2020).
  * [23] Carroll S 2004 Spacetime and Geometry (Addison Wesley)
  * [24] Paper mentioned before in [9, 10, 16, 17, 18, 19].
  * [25] Uzan J-P and Lehoucq R 2001 Eur. J. Phys. 22 371.
  * [26] Dahia F, Romero C, da Silva L F P and Tavakol R 2007 J. Math. Phys. 48 072501.
  * [27] See the sign convention part of Carroll [23] and Hartle (Hartle J B 2003 Gravity, An Introduction to Einstein’s General Relativity (Pearson Education).
  * [28] Mariam Bouhmadi-López et al JCAP07(2020)066
  * [29] arXiv:1911.04305v1.
  * [30] Stuart L. Shapiro and Saul A. Teukolsky 1983 Black Holes, White Dwarfs, and Neutron Stars (WILEY-VCH Verlag GmbH $\&$ Co. KGaA).
  * [31] Guha S and Chakraborty S 2010 Gen. Relativ. Grav. 42 1739.
  * [32] Sarbari Guha and Pinaki Bhattacharya 2012 J. Phys.: Conf. Ser. 405 012017.


# Determining Atmospheric Electric Fields using MGMR3D

T. N. G. Trinh<EMAIL_ADDRESS>Physics Education Department, School of
Education, Can Tho University, Campus II, 3/2 Street, Ninh Kieu District, Can
Tho City, Viet Nam  O. Scholten<EMAIL_ADDRESS>Kapteyn Astronomical
Institute, University of Groningen, P.O. Box 72, 9700 AB Groningen,
Netherlands  Interuniversity Institute for High-Energy, Vrije Universiteit
Brussel, Pleinlaan 2, 1050 Brussels, Belgium  S. Buitink  Astrophysical
Institute, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium
Department of Astrophysics/IMAPP, Radboud University Nijmegen, Nijmegen, The
Netherlands  K. D. de Vries  Astrophysical Institute, Vrije Universiteit
Brussel, Pleinlaan 2, 1050 Brussels, Belgium  P. Mitra  Astrophysical
Institute, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium
Faculty of Physics, University of Warsaw - Warsaw, Poland  T. Phong Nguyen
Department of Physics, Can Tho University, 3/2 Street, Ninh Kieu, Can Tho City
94000, Vietnam  D. T. Si  Department of Education and Training, 3/2 Street,
Ninh Kieu District, Can Tho City, Viet Nam

###### Abstract

Cosmic-ray particles impinging on the atmosphere induce high-energy particle
cascades in air, an Extensive Air Shower (EAS), emitting coherent radio
emission. This emission is affected by the presence of strong electric fields
during thunderstorm conditions. To reconstruct the atmospheric electric field
from the measured radio footprint of the EAS we use an analytic model for the
calculation of the radio emission, MGMR3D. In this work we make an extensive
comparison between the results of a microscopic model for radio emission,
CoREAS, to obtain an improved parametrization for MGMR3D in the presence of
atmospheric electric fields, as well as confidence intervals. The approach to
extract the electric field structure is applied successfully to an event with
a complicated radio footprint measured by LOFAR during thunderstorm
conditions. This shows that, with the improved parametrization, MGMR3D can be
used to extract the structure of the atmospheric electric field.

cosmic rays; shower parameters; atmospheric electric fields; radio emission;
extensive air showers

##  I Introduction

Lightning is a common phenomena but the detailed understanding of its
generation and development is still unknown. One of the reasons why this topic
is still under investigation is that atmospheric electric fields inside
thunderclouds are difficult to measure. Existing measurements, from balloons
and airplanes, are limited because the measurements depend on the paths of the
balloons or the aircrafts will affect the nature of the thunderclouds. A new
method to determine atmospheric electric fields is using their influence on
radio emission emitted from extensive air showers [1]. Unlike the balloon and
airplane measurements, this unique tool is not limited by the wind conditions
and is sensitive to a large part of the cloud.

When a cosmic ray enters the atmosphere of the earth, it interacts with air
molecules and generates a particle cascade, called an extensive air shower
(EAS). The particles in the EAS move with velocities near the speed of light
and are concentrated in the thin shower front. In an EAS during fair-weather
conditions, called fair-weather showers, the electrons and positrons are
deflected in opposite directions by the Lorentz force caused by the
geomagnetic field. This induces an electric current in the shower front that
is transverse to the shower axis. Since the electric current changes as a
function of height, due to the change of the number of particles in the EAS,
radio waves are emitted [2, 3, 4]. A secondary, yet important, contribution to
the radio emission is caused by an excess of electrons in the shower front due
to knock-out from atmospheric molecules by shower particles. This process
creates a radio pulse that is linearly polarized but oriented radially to the
shower axis [5, 6]. This charge excess emission interferes with the
geomagnetic emission that is linearly polarized along the direction of the
Lorentz force.

An EAS that occurs during thunderstorm conditions when there are strong
atmospheric electric fields, is called a thunderstorm shower, and produces
radio emission that differs considerably from that of fair weather showers [1,
7, 8]. In thunderstorm showers, the atmospheric electric field exerts an
electric force on the electrons and positron that is usually much stronger
than the Lorentz force and this affects the radio emission. The component of
the electric field perpendicular to the shower axis will change the direction
and magnitude of the transverse current changing the shape of the radio
footprint. The electric field parallel to the shower axis accelerates the
electrons or positrons, depending on its direction, and thus they gain more
energy. As a result, these particles will generate additional low-energy
particles. However, these additional particles travel with a velocity smaller
than the speed of light and thus they trail far behind the shower front. Thus,
their radiation is not added coherently in the frequency range from 30 MHz to
80 MHz at which typical cosmic-ray airshower radio detectors operate [7].

There are several models to describe radio emission from air showers.
Microscopic models such as ZHAires [9] and CoREAS [10] are based on full
Monte-Carlo simulation codes. Macroscopic models such as MGMR [2], EVA [3] and
Radio Morphing [11] calculate the emission of the bulk of electrons and
positrons described as currents. SELFAS2 [12] is a mix of macroscopic and
microscopic approaches. Recently, MGMR3D [13] has been introduced that uses a
multi-dimensional parametrization of the current density in the air showers to
calculate the radiation field using Maxwell equations.

Because atmospheric electric fields influence the radio emission from air
showers, the radio footprint, as is measured on the ground, can be used to
determine the strength, direction, and the altitude dependence of the
atmospheric electric field along the path of the air showers [14]. Since we
have only a forward model, i.e. calculate the footprint by assuming an
electric field configuration, we have to apply a search algorithm to extract
the field structure. This only works with a fast model calculation that is
deterministic, i.e. non Monte-Carlo based. MGMR3D is potentially such a model.
However, MGMR3D relies on a parameterization of the emitting currents. To
gauge these parameters we perform an extensive comparison between MGMR3D and
the full Monte-Carlo code CoREAS.

In Ref. [13, 8] a parametrization was proposed based on CONEX-MC [15]
simulations. In particular it had been assumed that the induced drift
velocity, driving the currents, is inversely proportional to the air density,
$v_{d}\propto 1/\rho$. This assumes that the drift velocity is inversely
proportional to the friction force, the collision frequency. In Ref. [16] it
was shown that an inverse dependence of the drift velocity on $v_{d}\propto
1/\sqrt{\rho}$ yields a much better agreement with the results of microscopic
calculations. The $1/\sqrt{\rho}$ dependence is characteristic of an
equilibrium situation, terminal velocity, where the acceleration of the fast
electron by the Lorentz force is compensated by the average friction due to
collisions with air molecules [16]. This has led to a considerable improved
parametrization for the description of fair weather showers in MGMR3D [16].
This different proportionality for the drift velocity has made us to revisit
the parametrizations that enter in MGMR3D in the presence of atmospheric
electric fields, by making an extensive comparison with a large number of
CoREAS calculations for random shower directions and various structures of the
atmospheric electric fields. The improved parametrization allows us to extract
the fields for more complicated configurations that were inaccessible before.

In Section II we present a very short review of the most important aspects of
CoREAS that describes radio emission on a microscopic level and which we
consider as ‘the truth’ to which we gauge the parameters that enter in the
modeling of the currents entering in MGMR3D, as discussed in Section III. We
first review the essentials of the parametrization for fair-weather showers
[16] after which we present the discussion of the parametrization in the
presence of electric fields. In Section IV we apply the improved model to data
for the event which was not possible to analyze in Ref. [8] to extract the
atmospheric electric field.

##  II CoREAS simulations

To gauge the parameters in the MGMR3D model, discussed in Section III, we
perform extensive comparisons with the radio emission footprints from EAS as
calculated from CoREAS, which we regard as the Monte-Carlo truth. CoREAS is a
plug-in for the particle simulation code CORSIKA [17] and is based on a
microscopic description of the radiation mechanism. In CoREAS the radio
emission from each electron and positron in the shower is calculated, without
making any assumptions on the type of radiation. In other words, CoREAS does
not simulate the transverse-current and charge-excess components separately,
but it produces the complete radiation field from the emissions of the
individual particles. Therefore, CoREAS simulations are compute intensive with
running times of the order of days. Using the electric field EFIELD option
[18] in CORSIKA it is possible to simulate air showers passing through strong
electric-field regions. For this work the radiation profile is calculated in
the shower plane (the plane perpendicular to the shower axis passing through
the point where the shower axis touches ground), with the $\hat{x}$-axis in
the direction of ${\bf{v}\times\bf{B}}$, along the direction of the Lorentz
force, and $\hat{y}$-axis along
${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$. Here $\mathbb{v}$ is the
velocity of the shower and $\mathbb{B}$ is the geomagnetic field. We use the
shower plane instead of the ground plane since it allows for an easier
interpretation. To construct a two-dimensional radio footprint, we run
simulations for 160 antenna positions in a star-shaped pattern in the shower
plane with eight arms, where each arm contains 20 antennas with a spacing of
25 m [19].

To compare with the results from MGMR3D, we extract the values for the Stokes
parameters at the 160 antenna positions and thus capture both intensity and
polarization of the radio emission. Because the aim of this work is to develop
a realistic approach to extract electric fields from LOFAR [20] data, we
construct the Stokes parameters from the CoREAS simulations in the frequency
band from 30 MHz to 80 MHz. For each position in each simulated shower, the
Stokes parameters are calculated as

$I=\frac{1}{n}\sum_{i=0}^{n-1}\left(\left|\varepsilon_{i,{\bf{v}\times\bf{B}}}\right|^{2}+\left|\varepsilon_{i,\bf{v}\times\left({\bf{v}\times\bf{B}}\right)}\right|^{2}\right),$
(1)
$Q=\frac{1}{n}\sum_{i=0}^{n-1}\left(\left|\varepsilon_{i,{\bf{v}\times\bf{B}}}\right|^{2}-\left|\varepsilon_{i,\bf{v}\times\left({\bf{v}\times\bf{B}}\right)}\right|^{2}\right),$
(2)
$U+iV=\frac{2}{n}\sum_{i=0}^{n-1}\left(\varepsilon_{i,{\bf{v}\times\bf{B}}}\varepsilon^{*}_{i,\bf{v}\times\left({\bf{v}\times\bf{B}}\right)}\right).$
(3)

$\varepsilon_{i,x}$ is the complex-valued signal radiation field along
direction $x$, where $i$ denotes the sample number and
$x={\bf{v}\times\bf{B}}$ or $x=\bf{v}\times\left({\bf{v}\times\bf{B}}\right)$.
$n$ is the number of time samples centered around the pulse-peak position.
Stokes $I$ is the intensity of the radio emission. Stokes $Q$ and $U$ are used
to derive the linear polarization angle

$\psi=\frac{1}{2}\mathrm{tan}^{-1}\left(\frac{U}{Q}\right),$ (4)

and Stokes $V$ is the intensity of the circularly polarized fraction of the
radiation.

In order to perform chi-square fitting with realistic sensitivities, error
bars (confidence levels) were assigned to the values of the Stokes parameters
for the CoREAS simulated showers, given by

$\sigma_{k,\left(I,Q,U,\mathrm{or}\,V\right)}=2\sqrt{\frac{\xi_{I}}{n}(I_{k}+\xi_{I})},$
(5)

where $n$ is the number of time samples as enters in the equations for the
Stokes parameters Eq. 3, $I_{k}$ is the value of Stokes $I$ for a test antenna
at position $k$ and $\xi_{I}$ is taken at the level of the instrumental noise
measured in LOFAR antennas [19]. This definition of the error bars of the
Stokes parameter guarantees that the relative error bars at positions near the
core, where the intensity is usually large, is small while at large distances
where the intensity is small, the relative error bars become large.

##  III Parameterization in MGMR3D

In MGMR3D [13], the Maxwell’s equations are solved in the far-field
approximation to obtain the radio footprint, using a parametrized charge-
current density that depends on atmospheric depth, radial distance from the
shower axis, and distance from the shower front. In earlier work, presented in
Ref. [8], the parametrization of the charge-current density was obtained
through a comparison of selected cases with the results of CONEX-MC and CoREAS
calculations. The resulting parametrization was adequate for extracting the
structure of relatively simple atmospheric electric fields, however for a
thunderstorm EAS with a complicated radio footprints there were are large
discrepancies between the MGMR3D and CoREAS results [8].

Following the approach presented in Ref. [16] we tune the parametrization in
MGMR3D by making a statistical comparison of the calculated radio footprint
with the results of microscopic calculations for a large numbers of showers
passing through different electric-field structures. The advantage of such an
approach is that a comparison is made at the level of the important
observables for later applications. In addition this approach allows to asses
the intrinsic accuracy of obtaining the structure of the atmospheric electric
field.

Since an extensive discussion of the parametrization for fair-weather showers
is presented in Ref. [16], we will in this paper only summarize the essentials
of the parametrization for fair-weather showers in Section III.1. The results
[16] show that a good agreement between MGMR3D and CoREAS is obtained for
fair-weather showers. The parametrization of the charge-current cloud in the
presence of electric fields is discussed in Section III.2.

###  III.1 Parameterization for fair-weather showers

Under fair-weather conditions only the Lorentz force,
$\vec{F}_{\perp}=e\vec{v}\times\vec{B}$, is acting on the charged particles in
the EAS, causing the electrons and positrons to drift in opposite directions
and thus induce an electric current. As presented in [16, 21] the induced
drift velocity depends on the square-root of the air density and can be
parametrized as

$\vec{v}_{\perp}(h)=\vec{F}_{\perp}/F_{0}\frac{\sqrt{\rho(X_{\rm
max})/\rho(h)}}{((X_{\rm max}-X_{v})/(X_{z}(h)-X_{v})+3.)}\;,$ (6)

with $X_{v}=50$ g/cm2 and $F_{0}=62.5$ keV/m. This -on first sight- surprising
dependence on $\sqrt{\rho}$ is due to the fact that at the equilibrium
sideways ‘friction’ force, proportional to $\rho v_{\perp}^{2}$, equals the
transverse component of the Lorentz force. Since only the Lorentz force is
acting, the drift velocity $\vec{u}_{\perp}$ is still non-relativistic and it
is parameterized as

$\vec{u}_{\perp}(h)=c\vec{v}_{\perp}(h)/\sqrt{1+v_{\perp}^{2}/u_{0}^{2}}\;,$
(7)

where the parameter $u_{0}$ is taken much larger than typical fair-weather
drift velocities, resulting in a linear scaling of the drift velocity with the
Lorentz force.

In MGMR3D the number of charged particles at a radial position $r_{s}$ from
the shower axis and a longitudinal position $d_{h}$ from the shower front is
written as $Nf(d_{h})w/r_{s}$, with $w(r_{s})$ the normalized radial particle
distribution, and $f(d_{h})$ the normalized longitudinal particle
distribution. The current is the product of the charged particle number and
the drift velocity, where both depend on the penetration depth in the
atmosphere, $X_{z}(h)$. The current is thus expressed as

$\vec{J}_{\perp}(t_{s},x_{s},y_{s},h)={w(r_{s})\over
r_{s}}\,f(d_{h},r_{s})\,N_{c}(X_{z})\vec{u}_{\perp}(h)\;,$ (8)

at a radial distance $r_{s}=\sqrt{x_{s}^{2}+y_{s}^{2}}$ from the shower axis.
The parametrization of the longitudinal shower profile for the current in the
shower front is based on the Gaisser-Hillas formula [22] for the dependence of
the number of charged particles on $X_{z}$,

$N_{c}(X_{\rm z})=\left({X_{\rm z}-X_{0}\over X_{\rm
max}-X_{0}}\right)^{X_{\rm max}-X_{0}\over\gamma}e^{X_{\rm max}-X_{\rm
z}\over\gamma}$ (9)

where $\gamma$ is a parameter controlling the width of the distribution and
$X_{0}$ the first-interaction point. $X_{\rm max}$ is the atmospheric depth
where the number of secondary particles reaches a maximum in the air shower.

The spatial distribution of the charged particle density is parametrized as

$w(r_{s})=N_{w}\,\xi(\xi+1)^{-2.5}\;;f(d_{h},r_{s})=N_{f}\,{\eta\over
e^{\sqrt{\eta}}+1}\;,$ (10)

with $\eta={d_{h}/\lambda}$ where $d_{h}$ is the distance from the shower
front, and $\xi=r_{s}/R_{0}$ is a scaled radial distance where we introduced
the radiation radius $R_{0}$. The scaling factors $N_{w}$ and $N_{f}$ are
chosen such that the distributions are normalized, $\int w(r)\,dr=1$ and
$\int_{0}^{\infty}f(d_{h},r_{s})dd_{h}=1$.

It was observed that the optimum value for $R_{0}$ depends on the distance
from $X_{\rm max}$ to the shower impact point. For small distances there is an
almost linear dependence until a saturation value of $R_{0}=50$ m is reached
at a distance of about 5 km, independent of zenith angle. A full account of
the parametrization for fair weather showers in MGMR3D is presented in [16].

###  III.2 Parameterization for thunderstorm showers

Due to convection flows in clouds there is a build up of electric charge that
has a layered structure when considering large clouds. In some cases this can
spontaneously cause a lightning discharge. The charge layer structure depends
on detailed cloud dynamics and temperatures, see for a recent reference [23].
To extract the atmospheric field configuration we thus parameterize a semi-
realistic atmospheric field configuration as built in three different layers
where in each layer the electric field may have a different orientation and
strength. Including the boundary heights we thus arrive at nine parameters for
the field in addition to two parameters for the location of the shower core as
well as one for the energy of the cosmic ray and one for $X_{\rm max}$. We
assume that the arrival direction of the cosmic ray is known from the arrival
time of the radio signals in the different antennas. Since this adds up to a
large number of parameters it is imperative to limit the introduction of
additional parameters that are not likely to greatly influence the structure
of the extracted electric fields. The values for the radiation radius $R_{0}$,
the amount of charge excess, parameters defining the structure of the shower
such as the first interaction depth $X_{0}$, and the width of the shower
profile $\lambda$ are therefore kept fixed to the generic values used for
fair-weather showers. The values of these parameters are given in Table. 5 in
the Appendix (except for $R_{0}$ that was discussed in the previous section).

There are two parameters that are more specific for the cases in which there
is an electric field present. One is related to how quickly the induced
transverse currents adjust to an electric field that changes with height,
discussed in Section III.2.1. Another parameter is $u_{0}$ as introduced in
Eq. (7). In Section III.2.2 the optimal values for these additional parameters
are obtained by fitting an ensemble of showers simulated in CoREAS, following
the approach introduced in [16].

To verify the implementation of atmospheric electric fields in MGMR3D and to
test the superposition principle we present in Section III.2.3 a comparison
with CoREAS calculations is performed for more intricate field configurations.

####  III.2.1 Height-dependent electric fields

In the presence of atmospheric electric fields the net force acting on the
particles changes. The component of the electric field parallel to the shower
axis is assumed to be well below the runaway breakdown [24]. This component
can increase the number of secondary non-ultra-relativistic (with energies
well below a GeV) electrons while decreasing the number positrons or the other
way around, depending on the its orientation. However, these particles are
generally non relativistic and thus they trail far behind the shower front
[7]. For this reason, they do not contribute much to the radio radiation in
the frequency range from 30 MHz to 80 MHz [7] and thus this field component
can be ignored.

We have noticed that electric fields in excess of 50 kV/m, that are transverse
to the shower, do affect the structure of the shower profile. We see an
increase in the number of charged particles starting at the height where the
electric field perpendicular to the shower axis, $\vec{E}_{\perp}$, is
present. As a result, the value of $X_{\rm max}$ is modified in the presence
of a strong electric field. The reason for this effect is still unknown. To
model this we modify the parametrization of the shower profile as

$\displaystyle N_{c,E}(X_{\rm z})=N_{c}(X_{\rm
z})\left(1+0.0015\,\frac{\rho(0)}{\rho(h)}\frac{F_{\perp}^{2}}{|v\times
B|^{2}}\right)\;,$ (11)

where $F_{\perp}$ is given in Eq. (12).

As has already been argued, $\vec{E}_{\perp}$, has a large effect on the radio
emission. In the presence of atmospheric electric fields the transverse force
acting on the shower particles is

$\vec{F}=q(\vec{E}_{\perp}+\bf{v}\times\bf{B})\;.$ (12)

This changes both the magnitude and the direction of the drift velocity which
is parallel to $\vec{F}$. In strong atmospheric electric fields, the net
forces become large and thus the drift velocity increases. The simulations
show that the radio emission does not increase proportionally and shows
saturation effect for strong fields [7]. This can be understood from
relativistic effects where the energy is no longer proportional to the square
of the velocity. To account for this we follow Ref. [7] and introduce the
saturation effects in the drift velocity, Eq. (7), by means of the parameter
$u_{0}$. In Section III.2.2, Table 2, the value for this parameter is
discussed.

Atmospheric electric fields may vary strongly with height. Thus another
important parameter is the time, or distance, it takes for the drift velocity
to adapt to the changes in the electric field. To parametrize this we replace
$\vec{F}_{\perp}$ in Eq. (6) by an effective force, $\vec{\tilde{F}}_{\perp}$,
defined as

$\displaystyle\vec{\tilde{F}}(h)$ $\displaystyle=$ $\displaystyle\vec{F}(h)+$
(13)
$\displaystyle\sum_{j}\frac{1}{1+e^{(a+(h-H_{s}-h_{j})\,b)}}e\Delta\vec{E}_{j}\;,$

where

$a=D_{Es}\,H_{s}\,\rho_{\rm ground}/X_{EM}\;,$ (14)

and

$b=D_{Es}\,\rho(h_{j})/X_{EM}\;.$ (15)

Here $\Delta\vec{E}_{j}$ is the change in the electric field at height $h_{j}$
and $X_{EM}$ is the interaction length for an electromagnetic particle. The
parameter $D_{Es}$ (for E-smooth) governs the distance over which the drift
velocity (through the effective force) adapts to the change in the electric
field, if $D_{Es}$ is large the second form will reduce to a step function
$\theta(h+H_{s}-h_{j})$. The drift velocity will only adapt to the applied
force after a certain relaxation time and the parameter $H_{s}=700$ is chosen
such that the average shift in height, as shown later in Fig. 2 is minimal.

Figure 1: The ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$-component of the
current profile for $D_{Es}=4,6,8$ and $10$ as a function of atmospheric
depth.

Fig. 1 shows how the current along the
${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ direction changes when the EAS
with $X_{\rm max}=500$ g/cm2 passes through an electric field boundary for
different values of $D_{Es}$. An electric field of 50 kV/m is taken between a
height of 4.5 km ($X=590$ g/cm2) and the ground, oriented along the
${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$-direction. It shows that the
current changes faster when the value of $D_{Es}$ is larger. The drop-off of
the current at large values of $X$, i.e. close to the ground, is due to the
decrease in the number of shower particles.

To determine the optimal values for the parameters $D_{Es}$ and $u_{0}$ we
simulate a set of 94 vertical showers passing through a one-layer electric
field using the microscopic code CoREAS. For each shower there is an electric
field present from a certain height $h$ to the ground with certain strength
$E$ and angle $\alpha$ where $\alpha$ is the angle between the electric field
and ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$-direction. The height $h$,
the strength $E$, the angle $\alpha$, as well as $X_{\rm max}$ are chosen
randomly.

To select the optimal value for $D_{Es}$, we choose four different values for
$D_{Es}$, 4, 6, 8, and 10. As shown in Table. 1, the mean value of $\chi^{2}$
does not change much when changing $D_{Es}$ and thus we do not need to
determine $D_{Es}$ to higher precision. While keeping this parameter fixed and
using the height $h$, the strength $E$ and the angle $\alpha$ obtained from
CoREAS, we perform a chi-square fit of $X_{\rm max}$. We also keep $H_{s}$
fixed at 700 and $u_{0}=10$. By using a Levenberg-Marquardt minimization
procedure based on a steepest descent method, we optimize the value of $X_{\rm
max}$ by minimizing

$\chi^{2}=\sum_{\text{antenna}}\sum_{S=I}^{Q,U,V}\left(\frac{S_{k,C}-f_{r}S_{k,3D}}{\sigma_{k}^{2}}\right)\,,$
(16)

where $S_{k,C}$ are the Stokes parameters from CoREAS calculations for antenna
at position $k$ and $S_{k,3D}$ are the ones from MGMR3D. $\sigma_{k}$ is the
uncertainty as defined in Eq. (5) and $f_{r}$ is the normalization factor for
the radio intensity. The results indicate that the value of $D_{Es}$ has only
a minor effect on the fitted values of $X_{\rm max}$. However, as shown in
Table. 1, the $\chi^{2}$ has a shallow minimum for $D_{Es}=8$ . Therefore, we
keep $D_{Es}$ fixed to this value for all subsequent calculations. The results
of fitting also show that the value of $\chi^{2}$ does not depend on $X_{\rm
max}$ and the electric field configuration.

$D_{Es}$ | 4 | 6 | 8 | 10
---|---|---|---|---
Mean $\chi^{2}$ | 0.183 | 0.178 | 0.175 | 0.178
Standard deviation | 0.156 | 0.150 | 0.146 | 0.145
Table 1: Mean $\chi 2$ and its standard deviation for different values of
$DEs$

The other parameter which is kept fixed for all thunderstorm showers is
$u_{0}$ as introduced in Eq. (7). Following the same procedure as for
$D_{Es}$, fitting $X_{\rm max}$ while keeping $u_{0}$, $h,E$ and $\alpha$
fixed, the mean $\chi^{2}$ and the standard deviation are calculated for
several values of $u_{0}$ in the range from 0.3 to 10. We also see that the
parameter $u_{0}$ hardly influences the extracted values of $X_{\rm max}$. As
shown in Table 2, the mean value of $\chi^{2}$ decreases when $u_{0}$
increases. However, when $u_{0}$ is equal or larger than 1, it does not affect
the quality of the fit anymore, so we keep it fixed at the value of 10.

$u_{0}$ | 0.3 | 0.6 | 1 | 3 | 7 | 10
---|---|---|---|---|---|---
Mean $\chi^{2}$ | 0.188 | 0.177 | 0.176 | 0.176 | 0.175 | 0.176
Standard deviation | 0.164 | 0.149 | 0.147 | 0.146 | 0.147 | 0.146
Table 2: Mean $\chi 2$ and its standard deviation for different values of
$u_{0}$

####  III.2.2 Accuracy of determined electric fields

We first investigate the accuracy of extracting one-layer electric fields. For
all 94 vertical showers passing through a one-layer electric field simulated
by CoREAS as discussed in the previous section, we follow a fitting procedure
using MGMR3D to reconstruct the structure of the field and compare this with
the true values used in the CoREAS calculation. In this section the true
electric field is taken homogeneous from a top height $h_{\rm true}$ to the
ground with a strength $E_{\rm true}$, making an angle $\alpha_{\rm true}$
with the ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ axis. The
reconstructed values, obtained by fitting the radio footprint using MGMR3D,
carries a subscript ${\rm reco}$.

Since we have observed that fitting the field structure and $X_{\rm max}$ at
the same time sometimes results in ill-converging fit, we have taken the
following approach. We perform 20 reconstructions for fixed values for $X_{\rm
max}$ that vary in steps of 20 g/cm2 between 500 g/cm2 and 900 g/cm2. For each
construction, we fit the three parameters of the electric-field structure
$h,E$ and $\alpha$ while keeping $X_{\rm max}$ fixed. As discussed in Ref.
[8], from the Stokes parameters, we know the power value and the orientation
of the polarization vector, however not the sign of the electric field vector.
As a result, there are two solutions that give almost the same value of
$\chi^{2}$ but differ in the sign of the induced current and thus that of the
electric-field direction. For this reason, we change the obtained angle
$\alpha$ by 180∘ and fit again for $h,E$ and $\alpha$ while keeping $X_{\rm
max}$ fixed. We select the electric field structure and the value of $X_{\rm
max}$ that give the smallest value for $\chi^{2}$.

Figure 2: $\Delta h=h_{\rm reco}-h_{\rm true}$ as a function of the true
height, $h_{\rm true}$. The black squares give the average value per kilometer
bin while the error-bars denote the standard deviation.

In Fig. 2 the difference $\Delta h=h_{\rm reco}-h_{\rm true}$ is shown as
function of the true height where the electric field starts, $h_{\rm true}$,
chosen randomly between 3 km and 10 km. Lower values were not considered as
the agreement for these cases is close to perfect. Per height-bin of 1 km the
mean values and the standard deviations are given. For heights lower than 9 km
a rather good agreement is obtained between the reconstructed height, $h_{\rm
reco}$, and the true height where the electric field starts, $h_{\rm true}$
with a standard deviation of less than 1 km. For heights between 9 km and 10
km, the spread in $\Delta h$ increases to reach a standard deviation of about
1.3 km although the mean difference is vanishingly small. The explanation is
that at large altitudes there are not many particles, as the shower is still
very young, and we thus loose sensitivity. Even though the spread in heights
is increasing from 0.2 km at the height of 3.5 km to 1.3 km at 9.5 km it
should be noted that the relative errors, $\Delta h/h_{\rm true}$, stay more
constant, ranging from 5% to 13%, in a similar range as for the single layer
case.

Figure 3: $\Delta\alpha=\alpha_{\rm reco}-\alpha_{\rm true}$ as a function of
the work of the net force $W$. The black squares give the average value per
100-MeV bin while the error-bars denote the standard deviation.

The true value of the electric field $E_{\rm true}$ is also chosen randomly
between 5 kV/m and 70 kV/m while the angle $\alpha$ is randomly selected from
0 to 360∘. In Fig. 3 we plot $\Delta\alpha=\alpha_{\rm reco}-\alpha_{\rm
true}$ as a function of the work of the net force $W=F.H$. Here $F$ is the net
force of the Lorentz force and the true electric force and $H$ is the
thickness of the layer. For these one-layer cases, $H=h_{\rm true}$. When the
work of the net force of a layer is large, the amount of radio emission
emitted from this layer is also large and thus we are more sensitive to this
layer. As a result, the orientation of the electric field in this layer is
determined more accurately (or $\Delta\alpha$ is small) for large values of
$W$. Fig. 3 shows that the angle of the field is very well reconstructed, with
$\Delta\alpha=\alpha_{\rm reco}-\alpha_{\rm true}<4^{\circ}$. For work less
than 100 MeV, the deviation of $\Delta\alpha$ is about 1∘, while for larger
work of the net force, it almost vanishes.

Figure 4: $\Delta E=E_{\rm reco}-E_{\rm true}$ as a function of the true
strength of the electric field, $E_{\rm true}$. The black squares give the
average value per 10-kV/m bin while the error-bars denote the standard
deviation.

Fig. 4 shows the difference between the reconstructed strength of the electric
field and the true value, $\Delta E=E_{\rm reco}-E_{\rm true}$. When the
electric field is less than 40 kV/m, the reconstructed values from the MGMR3D
calculation agrees well with the true value, with differences less than 5
kV/m. When the field strength is large, the differences become large since the
current saturates at about 50 kV/m and thus we loose the sensitivity to the
field strength. However, the relative errors are less than 5% when the field
strenght is weaker than 50 kV/m and vary from 8% to 11% when the field
strength is large.

Figure 5: $\Delta X_{\rm max}=X_{\rm max}$reco$-X_{\rm max}$true as a function
of the true height $h_{\rm true}$. The black squares give the average value
per kilometer bin while the error-bars denote the standard deviation.

Fig. 5 shows the difference between the reconstructed values of $X_{\rm max}$
and the true values, $\Delta X_{\rm max}=X_{\rm max}$reco$-X_{\rm max}$true.
For heights lower than 7 km, there seems to be no sizable systematic
differences. The mean discrepancies are less than 10 g/cm2, comparable to the
fair-weather results for $X_{\rm max}$. For heights from 7 km to 10 km, the
mean difference is a bit larger, from 30 g/cm2 to 40 g/cm2, and increases with
the true height, but the standard deviation is smaller.

One other important parameter is the normalization factor $f_{r}$ for the
Stokes parameters between MGMR3D and CoREAS which is given in Eq. (16). For
fair-weather showers, the mean normalization factor is equal to 1. For
thunderstorm showers, the mean value of the normalization factor of these 94
showers is 0.63 and the standard deviation is 0.24. This means that the
radiation in MGMR3D is overestimated as compared to the CoREAS result. This
could be a reflection of the fact that the reconstructed electric fields are a
bit over estimated, as shown in Fig. 4.

As was shown in this section, there are no sizable systematic differences
between the reconstructed parameters for the electric field structure and
$X_{\rm max}$ using a MGMR3D fit and the Monte-Carlo truth as expressed using
CoREAS. The systematic differences are mostly smaller than the reconstruction
accuracy.

####  III.2.3 Multiple layers

To obtain test results for semi-realistic cases we have simulated showers
passing through a two-layer electric field. Layer 1 extents from height
$h_{1}$ to the height $h_{2}<h_{1}$. Layer 2 extents from height $h_{2}$ to
the ground. The height for layer 1, $h_{1}$, is selected in the range from 5
km to 10 km. The height for layer 2, $h_{2}$, is chosen in the range from 2 km
to 4.5 km. In reality the field strength in layer 1 is usually larger than the
one in layer 2. Thus, the field strength in layer 1, $E_{1}$, and layer 2,
$E_{2}$, are chosen randomly from 2 kV/m to 60 kV/m and from 2 kV/m to 30
kV/m, respectively, each with arbitrary angles $\alpha_{1}$ and $\alpha_{2}$
with respect to the ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$-axis. The
field strength $E_{2}$ is chosen smaller than $E_{1}$ since in reality the
electric field in the cloud is usually stronger than the one between the cloud
to ground. Within these constraints the parameters are chosen randomly to
simulate 91 vertical showers in CoREAS. In order to find the reconstructed
values using MGMR3D, following a similar procedure as was taken for the one-
layer case, we perform 20 separate parameter searches with different values of
$X_{\rm max}$. We fit the six parameters of the two-layer electric field while
keeping $X_{\rm max}$ fixed. We then repeat the procedure changing the values
of $\alpha$ by 180∘.

Figure 6: $\Delta h=h_{\rm reco}-h_{\rm true}$ as a function of the true
height, $h_{\rm true}$, for layer 1 and layer 2. The black squares give the
average value per kilometer bin while the error-bars denote the standard
deviation.

Fig. 6 shows the difference between the reconstructed and the true height,
$\Delta h=h_{\rm reco}-h_{\rm true}$, for both layers. We observe very similar
differences, of the same order of magnitude as shown in Fig. 2 for the single-
layer case. For small heights in layer 2 the standard deviation is of the
order of 100 m, but is larger, of the order of 1000 m, for larger heights in
layer 1. On average the top height is reconstructed at a too low altitude.
Similar to what is seen in the single-layer case, although the spread in
heights for layer 1 is much larger than that for layer 2, the relative errors
vary less, from 5% to 13%.

Figure 7: $\Delta E=E_{\rm reco}-E_{\rm true}$ as a function of the true
height, $h_{\rm true}$, for layer 1 and layer 2. The black squares give the
average value per $10-kV/m$ bin while the error-bars denote the standard
deviation.

Since we have found that the difference between reconstructed and true
electric field strength correlates most strongly with the upper height of the
layer, we show in Fig. 7 $\Delta E=E_{\rm reco}-E_{\rm true}$ as function of
the height. We observe that the standard deviation of $\Delta E$ differs for
the two layers. While the mean $\Delta E$ is constant and slightly negative
for the lower layer, the mean value for $\Delta E$ is positive and decreases
with height for the upper layer. The spreads in layer 1 are larger than that
in layer 2 because the field strength in layer 1 is stronger than the one in
layer 2 and for large field strength we loose sensitivity as shown in the one-
layer case.

Figure 8: $\Delta\alpha=\alpha_{\rm reco}-\alpha_{\rm true}$ as a function of
the true height, $h_{\rm true}$, for layer 1 and layer 2. The black squares
give the average value per 50-MeV bin while the error-bars denote the standard
deviation.

Fig. 8 shows that there is a correlation between $\Delta\alpha$ and the work
of the net force $W$. When the work of the net force of a layer is larger than
50 MeV, which is usually seen in layer 1, the amount of radio emission emitted
from this layer is large and thus we are quite sensitive to this layer.
Therefore, the orientation of the electric field in this layer is determined
rather accurately or $\Delta\alpha$ is small. In contrast, when the work of
the net force is smaller than 50 MeV, which is usually seen in layer 2, the
amount of radiation from this layer is small and thus we loose sensitivity to
the polarization or the spread in $\Delta\alpha$ is large.

Figure 9: $\Delta X_{\rm max}=X_{\rm max}$reco$-X_{\rm max}$true as a function
of the true height of the top layer $h_{1\rm true}$. The black squares give
the average value per kilometer bin while the error-bars denote the standard
deviation.

Fig. 9 shows that the differences for $X_{\rm max}$ are on average positive
with a mean difference of about 50 g/cm2, about equal to the standard
deviation. The fact that $X_{\rm max}$ is deeper in the atmosphere is probably
correlated with the, on average, negative values for $\Delta h$ shown in Fig.
6.

Similar to the one-layer showers, the mean norm factor and its standard
deviation of the two-layer showers are $0.62$ and $0.19$, respectively.

To conclude, in general the reconstructed values are scattered around the true
values and the strongest bias is seen in the obtained values for $X_{\rm max}$
that are, however, of little interest in the study of atmospheric electric
fields.

##  IV  Extracted electric field configuration from LOFAR data using MGMR3D

To test the proposed approach for reconstructing the atmospheric electric
field from data, we have applied the reconstruction procedure to LOFAR data
for a cosmic ray event that has been recorded under thunderstorm conditions
[8]. We selected an event for which the radio footprint is particularly
complicated, to the extent we could not reach satisfactory results in an
earlier application of the procedure, event 1 from Ref. [8]. The shower axis
for this event is inclined with a zenith angle of $\theta=39.2^{\circ}$. As
shown in Ref. [8], there were large differences between CoREAS and the
reconstruction for this event using the older version of MGMR3D.

We have applied the present reconstruction method to recover the electric
field for this event. Similar to what has been done in Ref. [8], we model the
electric field with three layers. The core position is kept fixed at the
position estimated from the particle-detector array LORA [25]. We perform 20
reconstructions for fixed values for $X_{\rm max}$ that vary in steps of 20
g/cm2 between 500 g/cm2 and 900 g/cm2. For each construction, we fit nine
parameters of the electric-field structure while keeping $X_{\rm max}$ fixed.
At the end we select the electric field structure and the value of $X_{\rm
max}$ that give the best fit result. This approach is similar to the
reconstruction procedure used in Section III.2.2 and Section III.2.3. The
results of the MGMR3D calculation that gives the best fit is shown in Fig. 10.
The electric field reconstructed from this calculation is plugged in CoREAS to
confirm the procedure. Since in CoREAS, $X_{\rm max}$ is an output which
cannot be chosen before running simulations, we simulate 20 showers with the
reconstructed electric-field structure and select the simulation which gives
the smallest value of $\chi^{2}$, shown in Fig. 11. The values of the
extracted electric field parameters, the value of $X_{\rm max}$, the
reconstructed energy and the values of $\chi^{2}$ are given in Table. 3.

With the new version of MGMR3D, the event can be reconstructed rather well by
a three-layer electric field. As shown in Fig. 10 and Fig. 11, the differences
between MGMR3D and CoREAS are relatively small. It can be seen from the first
panels of these two figures that Stokes I, i.e. the radiation intensity, given
by both MGMR3D and CoREAS fits the LOFAR data well, except for distances near
the core of the shower. At these close distances, the intensity given by
CoREAS is a bit lower than the data, with the peak shifted to a slightly
larger distance. The intensity at small distances is rather small because of a
destructive interference between the radiation from layer 2 and layer 3. The
peak in the Stokes I at a distance of 175 m from the shower axis is reproduced
well in the MGMR3D calculation. Stokes $Q$ and Stokes $U$ represent the linear
polarization of the radiation. The second and the third panels of Fig. 10 and
Fig. 11 show the ratio $Q/I$ and $U/I$, respectively. The polarization of the
radiation is complicated, changing with the distances from the shower axis.
Stoke $V$, presenting the circular polarization, is plotted in the last panels
of Fig. 10 and Fig. 11 in terms of the ratio $V/I$. There is a large amount of
circular polarization which is reflected in the change in the direction of the
electric field between three layers. The linear polarization and the circular
polarization are also reproduced quite well except at small distances, within
100 m from the core of the shower where the radio intensity is rather small.
As shown in Table. 3 the MGMR3D and CoREAS results correspond to an energy of
the cosmic ray of $1\cdot 10^{8}$ GeV and $1.07\cdot 10^{8}$ GeV. This energy
differs from what is deduced from the particle detectors, $2.67\cdot 10^{7}$
GeV. The reason for this could be that the response of the LORA detectors is
affected by the thunderstorm.

Since we can determine the components of the electric fields perpendicular to
the direction of the shower, we can derive the purely horizontal components of
the field. We decompose $\textbf{E}_{\perp}$ into two components
$\textbf{E}_{\textbf{vxz}}$ and $\textbf{E}_{\textbf{vx(vxz)}}$ along vxz and
vxz and shown in Table. 4. The purely horizontal components
$\textbf{E}_{\textbf{vxz}}$ determined for this event are small because this
shower is inclined, with the zenith angle of $\theta=39.2^{\circ}$, since one
expects the atmospheric electric field to have a strong vertical component.
This component in layer 1 is larger than that in other lower layers, which is
one would expect because the charge layers are not purely horizontal or the
event occurred at the edge of the charged layer. The heights where the
electric fields change could be the positions of different charged layers.
There could be a negative-charge layer at 3.38 km and positive-charge layers
at 8 km and 1.85 km altitude. This is similar to what has been determined from
Lightning Mapping Array (LMA) observations of a flash that occurred near the
LOFAR core in June 2019. Details about the interpretation of this flash can be
found in Ref. [8].

The result of performing a reconstruction for this complicated-footprint shows
that we can get a stable result from the new version of MGMR3D that captures
the main structures seen in intensity and polarization. Using the
reconstructed fields in the microscopic calculation gives results for the
complicated radio footprint that show the same, very non-trivial, structures.

Layer | h [km] | E [kV/m] |  $\alpha$ [∘]
---|---|---|---
1 | 8.012 | 34.900 | 157.3
2 | 3.381 | 99.377 | -64.6
3 | 1.851 | 71.756 | 102.2
$X_{\rm max}$3D[g/cm2] | 860
$X_{\rm max}$C[g/cm2] | 848
Energy3D[GeV] | $1.0\times 10^{8}$
EnergyC[GeV] | $1.07\times 10^{8}$
$\chi^{2}_{\rm 3D}$ | $2.39$
$\chi^{2}_{\rm C}$ | $3.76$
Table 3: Nine extracted electric field parameters, $X_{\rm max}$ and the reconstructed energy of event 1 Layer | h [km] |  $\textbf{E}_{\textbf{vxz}}$ [kV/m] |  $\textbf{E}_{\textbf{vx(vxz)}}$ [kV/m]
---|---|---|---
1 | 8.012 | 24.721 | 24.635
2 | 3.381 | -5.547 | -99.222
3 | 1.851 | -12.460 | 70.665
Table 4: The components $\textbf{E}_{\textbf{vxz}}$ and
$\textbf{E}_{\textbf{vx(vxz)}}$ of the electric fields determined from event 1
Figure 10: The results of MGMR3D calculation for Stokes parameters (filled
blue dots) are compared to LOFAR data (filled red circles) for event 1.
$\sigma$ denotes one standard deviation error. Figure 11: Same as Fig. 10 but
for the comparison between CoREAS simulation and LOFAR data for event 1

##  V Summary and Conclusions

In this work, we have made an comprehensive comparison between the radio
emission during thunderstorm conditions from MGMR3D and from CoREAS. As a
result, we improve the parameterization for MGMR3D in the presence of
atmospheric electric fields.

With the parameters determined, we have generated a large number of radio
footprints for showers passing through a layered electric field using the
microscopic code CoREAS. From these footprints we have reconstructed the
electric field structure using MGMR3D and compared these with the true
structure that was used for generating the footprint. This shows that the
layer heights can be reconstructed with an accuracy of better that $\pm 1$ km,
the field strengths with $\pm 10$ kV/m, and orientation angles within a few
degrees. Exceptions are the orientation of weak fields for cases in which
there are also strong fields involved as well as cases where the top height of
the electric field layer lies at a height where the shower is still young with
relatively few charged particles. In the region where we are sensitive to the
parameters of the electric fields, the differences between MGMR3D and CoREAS
are rather small for all vertical showers passing through a one-layer electric
field or a two-layer one. Large discrepancies are observed for the cases where
we loose sensitivity to the parameters of the electric fields seen in both
MGMR3D and CoREAS.

We have applied the approach to extract electric field structure to an event
having a complicated radio footprint as measured by LOFAR and we are able to
reconstruct the main features of this event. Therefore, it can be concluded
that MGMR3D can be used to reconstruct the structure of the electric fields by
using the radio emission emitted from air showers passing through
thunderclouds. This method of determining atmospheric electric fields will
help to study about the process of lightning initiation and propagation.

## Acknowledgement

This research is funded by Vietnam National Foundation for Science and
Technology Development (NAFOSTED) under grant number 103.01-2019.378. K. D. de
Vries was supported by the European Research Council under the EU-ropean
Unions Horizon 2020 research and innovation programme (KdV, grant agreement No
805486).

##  Appendix A Parameter values

Parameter | $J_{Q}$ | $X_{0}$ | $\lambda$ |
---|---|---|---|---
Value | 0.21 | 100 | 100 |
Table 5: Fixed values of parameters

##  Appendix B Programming details

The program can be downloaded from MGMR3D-v3 as a zip file. Make sure you run
version 3.

## References

  * Schellart  _et al._ [2015] P. Schellart _et al._ , Probing atmospheric electric fields in thunderstorms through radio emission from cosmic-ray-induced air showers, Phys. Rev. Lett. 114, 165001 (2015).
  * Scholten  _et al._ [2008] O. Scholten, K. Werner, and F. Rusydi, A macroscopic description of coherent geo-magnetic radiation from cosmic-ray air showers, Astroparticle Physics 29, 94 (2008).
  * Werner  _et al._ [2012] K. Werner, K. D. de Vries, and O. Scholten, A realistic treatment of geomagnetic cherenkov radiation from cosmic ray air showers, Astroparticle Physics 37, 5 (2012).
  * Scholten and Werner [2009] O. Scholten and K. Werner, Macroscopic model of geomagnetic-radiation from air showers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 604, S24 (2009), {ARENA} 2008.
  * Askaryan [1962] G. Askaryan, Excess Negative Charge of the Electron-Photon Shower and Coherent Radiation Originating from It. Radio Recording of Showers under the Ground and on the Moon, J. Phys. Soc. Japan Vol: 17, Suppl. A-III (1962).
  * de Vries  _et al._ [2010] K. D. de Vries _et al._ , The lateral distribution function of coherent radio emission from extensive air showers: Determining the chemical composition of cosmic rays, Astroparticle Physics 34, 267 (2010).
  * Trinh  _et al._ [2016] T. N. G. Trinh _et al._ , Influence of atmospheric electric fields on the radio emission from extensive air showers, Phys. Rev. D 93, 023003 (2016).
  * Trinh  _et al._ [2020] T. N. G. Trinh, O. Scholten, S. Buitink, U. Ebert, B. M. Hare, P. R. Krehbiel, H. Leijnse, A. Bonardi, A. Corstanje, H. Falcke, T. Huege, J. R. Hörandel, G. K. Krampah, P. Mitra, K. Mulrey, A. Nelles, H. Pandya, J. P. Rachen, L. Rossetto, C. Rutjes, S. ter Veen, and T. Winchen, Determining electric fields in thunderclouds with the radiotelescope lofar, Journal of Geophysical Research: Atmospheres 125, e2019JD031433 (2020), e2019JD031433 10.1029/2019JD031433.
  * Alvarez-Muñiz  _et al._ [2012] J. Alvarez-Muñiz, W. R. C. Jr., and E. Zas, Monte carlo simulations of radio pulses in atmospheric showers using {ZHAireS}, Astroparticle Physics 35, 325 (2012).
  * Huege  _et al._ [2013] T. Huege, M. Ludwig, and C. W. James, Simulating radio emission from air showers with coreas, AIP Conference Proceedings 1535, 128 (2013).
  * Zilles  _et al._ [2020] A. Zilles, O. Martineau-Huynh, K. Kotera, M. Tueros, K. de Vries, W. Carvalho, V. Niess, N. Renault-Tinacci, and V. Decoene, Radio Morphing: towards a fast computation of the radio signal from air showers, Astropart.Phys. 114, 10 (2020).
  * Marin and Revenu [2012] V. Marin and B. Revenu, Simulation of radio emission from cosmic ray air shower with {SELFAS2}, Astroparticle Physics 35, 733 (2012).
  * Scholten  _et al._ [2018] O. Scholten, T. N. G. Trinh, K. D. de Vries, and B. M. Hare, Analytic calculation of radio emission from parametrized extensive air showers: A tool to extract shower parameters, Phys. Rev. D 97, 023005 (2018).
  * Schellart  _et al._ [2014] P. Schellart _et al._ , Polarized radio emission from extensive air showers measured with lofar, Journal of Cosmology and Astroparticle Physics 2014 (10), 014.
  * Bossard  _et al._ [2001] G. Bossard, H. J. Drescher, N. N. Kalmykov, S. Ostapchenko, A. I. Pavlov, T. Pierog, E. A. Vishnevskaya, and K. Werner, Cosmic ray air shower characteristics in the framework of the parton-based gribov-regge model nexus, Phys. Rev. D 63, 054030 (2001).
  * Mitra [2021] P. Mitra,  _High Precision Reconstruction of Air Shower Properties With Dense Radio Arrays_ , Ph.D. thesis, Vrije Universiteit Brussel (2021).
  * Heck  _et al._ [1998] D. Heck _et al._ ,  _CORSIKA: a Monte Carlo code to simulate extensive air showers_ (1998).
  * S. Buitink  _et al._ [2007] S. Buitink _et al._ , Amplified radio emission from cosmic ray air showers in thunderstorms, A&A 467, 385 (2007).
  * Buitink  _et al._ [2014] S. Buitink _et al._ , Method for high precision reconstruction of air shower ${X}_{\mathrm{max}}$ using two-dimensional radio intensity profiles, Phys. Rev. D 90, 082003 (2014).
  * van Haarlem  _et al._ [2013] M. P. van Haarlem _et al._ , Lofar: The low-frequency array, A&A 556, A2 (2013).
  * [21] Terminal velocity, https://www.grc.nasa.gov/WWW/K-12/airplane/termv.html.
  * Gaisser and Hillas [1977] T. K. Gaisser and A. M. Hillas, Reliability of the Method of Constant Intensity Cuts for Reconstructing the Average Development of Vertical Showers, in  _International Cosmic Ray Conference_ , International Cosmic Ray Conference, Vol. 8 (1977) p. 353.
  * Salvador  _et al._ [2021] A. Salvador, N. Pineda, J. Montanyà, J. A. López, and G. Solà, Thunderstorm charge structures favouring cloud-to-ground lightning, Atmospheric Research 257, 105577 (2021).
  * Dwyer and Uman [2014] J. R. Dwyer and M. A. Uman, The physics of lightning, Physics Reports 534, 147 (2014), the Physics of Lightning.
  * Thoudam  _et al._ [2014] S. Thoudam _et al._ , Lora: A scintillator array for {LOFAR} to measure extensive air showers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 767, 339 (2014).