# The centrifugal acceleration and the Y-point of the Pulsar Magnetosphere

Shinpei Shibata Faculty of Science, Yamagata University
Kojirakawa 1-4-12, Yamagata 990-8560, JAPAN  Shota Kisaka Physics Program,
Graduate School of Advanced Science and Engineering, Hiroshima University,
Higashi-Hiroshima 739-8526, Japan

###### Abstract

We investigate the centrifugal acceleration in an axisymmetric pulsar
magnetosphere under the ideal-MHD approximation. We solved the field-aligned
equations of motion for flows inside the current sheet with finite thickness.
We find that flows coming into the vicinity of a Y-point become super fast.
The centrifugal acceleration takes place efficiently, and most of the Poynting
energy is converted into kinetic energy. However, the super fast flow does not
provide enough centrifugal drift current to open the magnetic field. Opening
of the magnetic field is possible by the plasmas that are accelerated in the
azimuthal direction with a large Lorentz factor in the closed field region. We
find that this acceleration takes place if the field strength increases toward
the Y-point from inside. The accelerated plasma is transferred from the closed
field region to the open field region by magnetic reconnection with plasmoid
emission. We also estimate the Lorentz factor to be reached in the centrifugal
wind.

Rotation powered pulsars (1408) — High energy astrophysics (739) —
Magnetohydrodynamics (1964) — Plasma astrophysics (1261)

##  1 Introduction

Particle acceleration by pulsars is a long-standing issue in astrophysics. In
an early stage of investigation, the idea of the relativistic centrifugal wind
is intensively studied (Michel, 1969; Goldreich & Julian, 1970; Li & Melrose,
1994). Since the centrifugal acceleration is based on a corotational motion of
the magnetospheric plasmas, plasma density is assumed to be high, and the
ideal-MHD condition holds everywhere. In reality, however, acceleration
mechanism depends on plasma density, and more specifically on where and how
much electron-positron pairs are created. Then particle acceleration by field-
aligned electric field also interests us to explain the pulsed emission from
radio to gamma-ray (Jackson, 1976; Cheng et al., 1976; Takata et al., 2011;
Ruderman & Sutherland, 1975; Shibata, 1997; Timokhin, 2010). Recently,
Particle-in-Cell (PIC) simulations are applied to understand the global
structure of the pulsar magnetosphere (Chen & Beloborodov, 2014; Belyaev,
2015; Hakobyan et al., 2023; Bransgrove et al., 2023; Hu & Beloborodov, 2022).
They are able to treat both the field-aligned particle acceleration and the
centrifugal acceleration. However, the mechanism of the centrifugal
acceleration still remains unresolved. In this paper, we revisit this problem,
and compare with recent PIC simulations.

The model we consider is an axisymmetric steady magnetosphere with a magnetic
dipole moment at the center parallel to the angular velocity. The plasma
number density is assumed to be much higher than the Goldreich-Julian density
everywhere, so that the ideal-MHD approximation holds. We do not deal with
pair creation processes.

The magnetospheric plasma tends to corotate with the star. Strict corotation
would lead us to a divergence of the Lorentz factor at the light cylinder,
i.e.,

$v_{\varphi}=\varpi\Omega_{*}\rightarrow c\ \ \ \mbox{as }\varpi\rightarrow
c/\Omega_{*}\equiv R_{\rm L},$ (1)

where $\Omega_{*}$ is the angular velocity of the star, $\varpi$ is the axial
distance, and $R_{L}$ is the radius of the light cylinder. In reality, just
within the light cylinder, the inertia would be so large that the centrifugal
drift current would open the magnetic field lines. The plasma accelerated in
azimuthal direction with a large Lorentz factor would finally be thrown away
along the open field lines. This is the idea of the centrifugal acceleration,
just like a trebuchet.

However, after some works (eg., Michel, 1969; Begelman & Li, 1994; Takahashi &
Shibata, 1998), it is concluded that the centrifugal acceleration is
inefficient, and flows on the open field lines are Poynting energy dominant.
There is toroidal magnetic field $B_{\varphi}$ due to the global current loop
as shown in Figure 1. The azimuthal velocity does not follow the corotation
but does the iso-rotation law, $v_{\varphi}=\Omega_{*}\varpi+\kappa
B_{\varphi}$, where $\kappa$ is a scalar function to be determined. The open
magnetic field lines trail backward due to rotation, resulting in negative
$\kappa B_{\varphi}$, meaning departure from the corotation. Thus the
slingshot acceleration fails. A typical Lorentz factor to be reached is
$\sigma_{*}^{1/3}$ (Michel, 1969, and elsewhere) where $\sigma_{*}=B^{2}/4\pi
nmc^{2}$, is the magnetization parameter, and $B$ and $n$ are the field
strength and the number density, respectively. An expression
$\sigma_{*}=\gamma_{\rm max}/2{\cal M}$ may be useful, where $\gamma_{\rm
max}=eB_{\rm L}/mc\Omega_{*}$ is the maximum reachable Lorentz factor, ${\cal
M}=n/(\Omega_{*}B_{\rm L}/2\pi ec)$ is the multiplicity of the plasma, and
$B_{\rm L}$ is the light cylinder field. The value of $\sigma_{*}$ gives a
typical Lorentz factor when most of the Poynting energy is converted into
plasma. A flow with the Lorentz factor of $\sigma_{*}^{1/3}$ is obviously
Poynting energy dominant.

Figure 1: Schematic picture of an axisymmetric magnetosphere with close and
open magnetic field lines. The small arrows indicate poloidal currents, which
form a loop in each hemisphere starting from the star in lower latitudes and
go back to the star in higher latitudes. The top of the close field lines is a
junction of the current sheets in the force-free limit, called a Y-point.

The poloidal current loop has an essential role to carry Poynting energy away
from the neutron star, since it generates $B_{\varphi}$ and the Poynting flux,
$c{\boldsymbol{E}}_{\perp}\times{\boldsymbol{B}}_{\varphi}/4\pi$, where
${\boldsymbol{E}}_{\perp}$ is the electric field across the magnetic field
produced by the electromotive force of the star. Due to symmetry, however, the
equatorial plane is the place where the toroidal field vanishes. In a region
where the toroidal field becomes very weak, the plasma motion may become very
close to the pure corotation to get large Lorentz factors. Such a place is
where the equatorial current sheet touches the closed field region, a vicinity
of a so-called Y-point (Figure 1).

Meanwhile, highly accurate PIC simulations are available to provide detailed
analysis for the centrifugal acceleration (Chen & Beloborodov, 2014; Belyaev,
2015; Hakobyan et al., 2023; Bransgrove et al., 2023; Hu & Beloborodov, 2022).
It has been shown that most of the space outside the equatorial current sheet
is Poynting energy dominant. However, Poynting energy is converted to kinetic
energy in the vicinity of the Y-point and the equatorial current sheet. It is
interesting that Hu & Beloborodov (2022) show an acceleration in azimuthal
direction such as expected in the centrifugal acceleration. Another
interesting feature of their simulations is ejection of plasmoids.

In this paper, we study the centrifugal acceleration that is expected in a
vicinity of the Y-point. It will be shown that super-fast flows appear inside
the thin current layer. The centrifugal drift current opens the poloidal
magnetic fields in the vicinity of the Y-point. We will propose a self-
consistent structure of the centrifugal wind.

##  2 Magnetic field opening by the centrifugal drift current

Let us estimate the Lorentz factor to open the magnetic field lines by
centrifugal drift currents. The region we consider is a vicinity of the
Y-point. In the following sections, we use the cylindrical coordinate
$(z,\varpi,\varphi)$ with the unit vectors denoted by
$({\boldsymbol{e}}_{z},{\boldsymbol{e}}_{\varpi},{\boldsymbol{e}}_{\varphi})$.
The centrifugal drift velocity, which has opposite directions depending on the
charge sign, would be

${\boldsymbol{v}}_{\rm d}=c\frac{{\boldsymbol{F}}\times{\boldsymbol{B}}}{\pm
eB^{2}}\approx\pm\frac{\gamma mc^{2}}{eB_{\rm Y}R_{\rm
Y}}c{\boldsymbol{e}}_{\varphi},$ (2)

where the centrifugal force on an electron of mass $m$ and the Lorentz factor
$\gamma$ is evaluated as ${\boldsymbol{F}}=(\gamma mc^{2}/R_{\rm
Y}){\boldsymbol{e}}_{\varpi}$, $R_{\rm Y}\sim R_{\rm L}$ is the axial distance
of this region, and $B_{\rm Y}$ is the magnetic field that is weakened
significantly as compared with $B_{\rm L}$. The drift current density may be
$j_{\varphi}=2en|{\boldsymbol{v}}_{\rm d}|\approx 2nc(\gamma mc^{2}/B_{\rm
Y}R_{\rm Y})$. Applying Ampère’s law for a closed loop in the $\varpi-z$ plane
around the Y-point, we have $(4\pi/c)j_{\varphi}\Delta=2B_{\rm p}^{\rm(out)}$
for a condition of opening the magnetic field, where $B_{\rm
p}^{\rm(out)}\approx B_{\rm L}$ is the poloidal field strength outside the
current layer, and $\Delta$ is a typical thickness. With these conditions, we
have the Lorentz factor to open the magnetic field,

$\gamma_{c}\approx\frac{(B_{\rm p}^{\rm(out)})^{2}}{4\pi nmc^{2}}\frac{R_{\rm
Y}}{\Delta}\frac{B_{\rm Y}}{B_{\rm p}^{\rm(out)}}=\sigma_{*}\frac{R_{\rm
Y}}{\Delta}\frac{B_{\rm Y}}{B_{\rm p}^{\rm(out)}}.$ (3)

We expect a thin layer, $\Delta\ll R_{\rm Y}$, and weak magnetic field,
$B_{\rm Y}\ll B_{\rm L}$. If $(R_{\rm Y}/\Delta)(B_{\rm Y}/B_{\rm L})\approx
1$, $\gamma_{\rm c}\approx\sigma_{*}$. It is one of our objectives to find
$\Delta$ and $B_{\rm Y}$ in a self-consistent manner, considering the
magnetohydrodynamics inside the current layer.

##  3 Current distribution in the current sheet

In the following calculations of field-aligned flows, we assume an
axisymmetric magnetosphere with the ideal-MHD condition.

Let us consider two flows: one runs just outside the equatorial current layer,
the other runs inside. We call them Flow A and Flow B, respectively, as shown
in Figure 2. What is distinctive of Flow B from Flow A is significant decrease
in the poloidal field strength as is flows near the Y-point. We expect that
the centrifugal acceleration is inefficient in Flow A, but efficient in Flow
B.

Figure 2: Two types of field-aligned flows. Flow A and Flow B, respectively,
runs outside and inside the current layer with the thickness $\Delta$. The
dashed curves are their stream lines. The solid curves are the magnetic field
lines. The thin red curves are the poloidal current lines. The $\varpi$-axis
indicates the axial distance. The thin dashed line $\varpi=x_{\rm A}$ is the
Alvén surface. The thin dotted line $\varpi=R_{\rm L}$ is the light cylinder.

Flows such as Flow A and Flow B can be solved by using the field-aligned
equations (Okamoto, 1978; Camenzind, 1986; Takahashi, 1991), provided that a
poloidal field is given. It must be noticed here that the solution of field-
aligned flows determines values of $B_{\varphi}$ and in turn the current
function $I=\varpi B_{\varphi}$, where stream lines of the poloidal current is
given by contours of $I$. The field-aligned equations are solved with a given
set of boundary conditions and the condition that the solution passes through
the fast critical point. As a result, $I$ is uniquely determined. On the other
hand, as was done in obtaining the force-free solution (Contopoulos et al.,
1999; Timokhin, 2006), the trans-field equation that gives the poloidal field
structure determines the current function $I$ so that the field lines are
regular on the Alfvén surface, which coincides with the light cylinder in the
force-free limit. Thus, the both systems of equations can determine $I$. The
pulsar problem is only solved simultaneously the field-aligned equations and
the trans-field equation, i.e., the two $I$’s obtained in the two ways must be
identical. This has never been done so far.

The force-free solution may be a good approximation for the field structure in
high and middle latitudes, because the field-aligned solutions give weak
centrifugal acceleration there. This is supported by the PIC simulations and
the Relativistic Magnetohydrodynamic simulations, which give the results very
similar to the force-free solution.

The problem is how the current runs inside the current layer and in the
vicinity of the Y-point.

For flows such as Flow A, the poloidal current lines almost coincide with the
flow stream lines, since the centrifugal acceleration is weak. If a strong
acceleration takes place on Flow B, the current stream lines crosses the flow
stream line causing the acceleration. In this case, some current lines
converge to the Y-point as indicated in Figure 2, forming a very thin current
sheet. This seems a realization of the infinitely thin current sheet that is
required by the regularity on the light cylinder in the force-free
approximation. Even if inertia is taken into account in the trans-field
equation, the inertial effect may not be significant on the Alfvén surface;
the acceleration become significant beyond the Alfvén surface. Therefore, the
thin current sheet on the Alfvén surface would be a genuine property of the
centrifugal wind, ensuring the regularity of the poloidal magnetic field.

##  4 Field-aligned flows in the current layer

In this section, we solve the field-aligned equations for Flow A and B. The
field-aligned equations are composed of the three conservation laws and the
iso-rotation law. Given a stream function $\psi$ of the poloidal magnetic
field, the particle flux $g(\psi)$, the angular momentum $\ell(\psi)$, and the
energy $\epsilon(\psi)$ are functions of $\psi$ only. In the iso-rotation law,
the angular velocity $\Omega(\psi)$ is also a function of $\psi$. These
quantities are all determined by a set of boundary conditions and the
condition that the flow passes through a fast-critical point.

The field-aligned equations are algebraic equations, and are joined into a
single expression called the Bernoulli function $\epsilon(\xi,v)$. Here, we
normalize the axial distance $\varpi$ and the poloidal 4-velocity $u_{\rm p}$
by the values at the Alfvén point, i.e., $\xi=\varpi/\varpi_{\rm A}$ and
$v=u_{\rm p}/u_{\rm A}$.

The field-aligned equations are algebraic equations, and are joined into a
single expression called the Bernoulli function $\epsilon(\xi,v)$. Here, we
normalize the axial distance $\varpi$ and the poloidal 4-velocity $u_{\rm p}$
by the values at the Alfvén point, i.e., $\xi=\varpi/\varpi_{\rm A}$ and
$v=u_{\rm p}/u_{\rm A}$. A non-dimensional form of the Bernoulli function
$\epsilon(\xi,v)$ with the definition of the Alfvén point (Okamoto, 1978;
Camenzind, 1986; Takahashi, 1991) is given below, and it is notable that a
solution is obtained as a curve $\epsilon(\xi,v)=$constant on the $\xi$-$v$
plane. The Bernoulli function $\epsilon(\xi,v)$ has two non-dimensional
parameters,

$\lambda=\frac{\Omega(\psi)\ell(\psi)}{\epsilon(\psi)},$ (4)

and

$\sigma=\frac{\Omega(\psi)^{2}B_{\rm p}\varpi^{2}}{4\pi mc^{3}g(\psi)},$ (5)

and it has an singularity at the Alfvén point, $\varpi_{\rm
A}=\sqrt{\lambda}c/\Omega(\psi)$. The 4-velocity at the Alfvén point becomes
$u_{\rm A}=(1-\lambda)\sigma_{\rm A}/\lambda$, where $\sigma_{\rm A}$ is the
value of $\sigma$ at the Alfvén point. The Bernoulli function has a form

$\left(\frac{\epsilon(\xi,v)}{mc^{2}}\right)^{2}={\frac{(1+u_{\rm
p}^{2})}{(1-\lambda)^{2}}}\ \
{\frac{[1-\lambda\xi^{2}-(1-\lambda)v\xi^{2}]^{2}}{(1-v\xi^{2})^{2}-\lambda\xi^{2}(1-v)^{2}}},$
(6)

where

$u_{\rm p}=\frac{1-\lambda}{\lambda}\sigma v.$ (7)

A family of solutions is easily drawn as contours of $\epsilon(\xi,v)$ for a
given set of $\lambda$ and $\sigma$.

We chose one solution among them that passes through the fast-critical point;
hereafter we denote the fast point as $(\xi_{\rm f},v_{\rm f})$. At the same
time, the flow energy $\epsilon_{\rm f}=\epsilon(\xi_{\rm f},v_{\rm f})$ is
obtained as if it is an eigenvalue.

The Bernoulli function $\epsilon(\xi,v)$ has the two parameters, $\lambda$ and
$\sigma$. However, the value of $\lambda$ is determined implicitly by given
inner boundary conditions. Firstly we assume $\lambda$, and obtain a critical
solution. Then we examine whether the solution satisfies the set of given
boundary conditions. If not, we change $\lambda$ until the boundary conditions
are satisfied. The inner boundary conditions are those for a plasma source
region. More specifically, for Flow A and Flow B, they may be such that the
plasma is injected well within the light cylinder $\xi\ll 1$ with the Lorentz
factor $\gamma\sim 1$.

Thus, $\sigma$ is the only parameter that determines the flow. It is
decomposed into two factors:

$\displaystyle\sigma$ $\displaystyle=$ $\displaystyle\sigma_{0}\hat{B},$ (8)
$\displaystyle\sigma_{0}$ $\displaystyle=$ $\displaystyle\frac{B_{\rm L}}{4\pi
mcg(\psi)}=\gamma_{\rm max}\left(\frac{v_{\rm p}/c}{{\cal M}}\right)_{\rm
in}\left(\frac{\Omega_{*}}{\Omega(\psi)}\right)\approx\frac{\gamma_{\rm
max}}{{\cal M}_{\rm in}},$ (9) $\displaystyle\hat{B}$ $\displaystyle=$
$\displaystyle\frac{B_{\rm p}\varpi^{2}}{B_{\rm L}R_{\rm L}^{2}},$ (10)

where $g(\psi)=n\kappa=nv_{\rm p}/B_{\rm p}$, and ${\cal M}_{\rm in}$ is the
parameter of the plasma number density at the inner boundary, i.e., ${\cal
M}_{\rm in}=n/[\Omega(\psi)B_{\rm p}/2\pi ce]$, evaluated at the inner
boundary. The configuration of the poloidal magnetic field only comes into the
parameter $\hat{B}\propto B_{p}\varpi^{2}$.

The key parameter that makes difference between Flow A and Flow B is
$\hat{B}$. The poloidal magnetic field would be something like dipole near the
star, and changes to something like radial near and beyond the light cylinder.
Mimicking this, a simple function for $\hat{B}$ is introduced as follows. The
magnetic stream function for a dipole field is $\psi=B_{\rm L}R_{\rm
L}^{3}\varpi/r^{2}$, where $r=(z^{2}+\varpi^{2})^{1/2}$. For the dipole field,
$\hat{B}$ is given by

$\hat{B}_{\rm dip}(x)=({\psi/B_{\rm L}R_{\rm L}^{2}})\left[4-3({\psi x/B_{\rm
L}R_{\rm L}^{2}})^{2/3}\right]^{1/2},$ (11)

where $x=\varpi/R_{\rm L}$. Since $\hat{B}$ is constant for a radial field, we
assume that it gradually changes to a constant, $\hat{B}_{\rm
radial}=\hat{B}_{\rm dip}(x_{\rm c})$ beyond a certain point $x_{\rm c}$. Then
we assume that $\hat{B}$ for Flow A is

$\hat{B}_{\rm a}(x)=S\left(\frac{x_{c}-x}{\delta_{c}}\right)\hat{B}_{\rm
dip}(x)+S\left(\frac{x-x_{c}}{\delta_{c}}\right)\hat{B}_{\rm radial},$ (12)

where $S(x)=1/[\exp(-x)+1]$ is the sigmoid function (see Figure 3).

Figure 3: Models of $\hat{B}(x)\propto B_{\rm p}\varpi^{2}$ for Flow A, Flow B
and Flow C.

As an example of Flow A, we take $\gamma_{\rm max}=10^{7}$ and ${\cal M}_{\rm
in}=10^{3}$, and therefore $\sigma_{0}=10^{4}$. The solutions in the $\xi$-$v$
plane are shown in Figure 4, where the critical solution is colored in red. It
is known that all the solutions pass through the Alfvén point. The denominator
of (6) must be positive, so that there is a forbidden region, which is bounded
by the two curves,

$v_{+}=\frac{1+\sqrt{\lambda}\xi}{\xi(\xi+\sqrt{\lambda})},\ \ \ \ \
v_{-}=\frac{1-\sqrt{\lambda}\xi}{\xi(\xi-\sqrt{\lambda})},$ (13)

indicated by the dashed curves in the figure. The critical points are the
points where the two curves, ${\partial\epsilon}/{\partial v}=0$ and
${\partial\epsilon}/{\partial\xi}=0$ meet. These curves are also shown with
the dots in Figure 4. Since $\hat{B}$ is essentially constant beyond the
Alfvén point, the X-type critical point, namely the fast point, goes to
infinity as in the case of the radial flow. After an iteration, the value of
$\lambda$ is determined so that the injection point locates at about
$\xi=0.1$. The obtained parameters of the critical solution is given in Table
1.

Figure 4: A family of the solutions and the critical solution (red) of Flow A. The dashed curves indicate the boundaries of the forbidden region. The green and orange dots indicate where $\partial\epsilon/\partial v=0$ and $\partial\epsilon/\partial\xi=0$, respectively.  Table 1:  Some parameters for the critical solution. Type | $\sigma_{0}$ | $\log\xi_{f}$ | $\log v_{f}$ | $\epsilon_{f}$ | $1-\lambda$ | $1-x_{\rm A}$ | $\gamma$
---|---|---|---|---|---|---|---
Flow A | $10^{4}$ | $\infty$ | $1.3275$ | $0.9905\times 10^{4}$ | $1.02\times 10^{-4}$ | $5.15\times 10^{-5}$ | $21.4{\rm(a)}$
Flow B | $10^{4}$ | $0.0525$ | $0.9075$ | $3.039\times 10^{3}$ | $4.0\times 10^{-4}$ | $2.0\times 10^{-4}$ | ${2.95\times 10^{3}}{\rm(b)}$
Flow C (c) | $10^{4}$ | $1$ | $1$ | $10^{2}$ | $4.0\times 10^{-4}$ | $2.0\times 10^{-4}$ | $10^{2}$
(a): evaluated at $\xi=100$, (b): evaluated at $\xi=3.14$, (c):$\lambda$ is
assumed.

The critical solution gives $\gamma=21.4$ at $\xi=100$, which is consistent
with the Michel’s solution $\gamma\rightarrow\epsilon_{\rm f}^{1/3}=21.5$ as
$\xi\rightarrow\infty$ for the radial flow. The total energy $\epsilon_{\rm
f}$ is $\approx\sigma_{0}$. In the left panel of Figure 5, we plot $\gamma$ as
a function of $\log\xi$ with $\epsilon_{\rm f}/mc^{2}$, which is constant.
Flow A is a Poynting dominant flow. The azimuthal velocity deviates from the
corotation as seen in the right panel of Figure 5.

Figure 5:  Properties of Flow A. Left panel: The particle energy $\gamma$
(solid curve), and the total energy $\epsilon_{\rm f}$ (dashed line) as
functions of $\log\xi$. Right panel: The azimuthal velocity $v_{\varphi}$ as a
function of $x=\varpi/R_{\rm L}$. The dashed curve indicates the corotation
for $x<1$ and the free emission tangent to the light cylinder for $x>1$.

Next, we consider Flow B that comes into a weak field region near the Y-point.
We expect a significant decrease of $\hat{B}$ beyond the Alfvén point. We
model $\hat{B}$ by

$\hat{B}_{\rm b}(x)=\left[(1-D)\ S\left(\frac{x_{\rm d}-x}{\delta_{\rm
d}}\right)+D\right]\hat{B}_{\rm a},$ (14)

where $D$ and $x_{\rm d}$, respectively, determine how much and where the
magnetic field decreases. If $D=1$, there is no decrease, but if $D<1$,
$\hat{B}$ drops down to $D\hat{B}_{\rm a}$ beyond $x_{\rm d}$. In the
following example, we use $D=0.01$, $x_{\rm d}=1.2$, and $\delta_{d}=0.1$, as
shown in Figure 3. The value of $\sigma_{0}$ is the same as Flow A, i.e.,
$\sigma_{0}=10^{4}$.

Figure 6: The same as Figure 4, but for the flow B.

There is a topological change in the solution curves as shown in Figure 6. An
X-type critical point, the fast point, appears just beyond the Alfvén point;
$\xi_{\rm f}=1.13$. The Lorentz factor glows up around the fast point from the
value at the injection ($\gamma=1.1$) to $2.95\times 10^{3}$, which
corresponds to 98% of the total energy $\epsilon_{\rm f}$ (Figure 7). In
contrast with Flow A, deviation of the azimuthal velocity from the corotation
is weaker, and $v_{\varphi}$ peaks at $\sim 0.85c$. More importantly the curve
eventually almost follows the analytic curve of the tangential ejection with
the light speed at the Alfvén radius. The flow is effectively corotation up to
the Alfvén point. The total energy $\epsilon$ is about 30% of Flow A. This is
reasonable because the energy flux of Flow B is on its way of decreasing down
to zero on the equator.

Figure 7: The same as Figure 5, but for Flow B.

It is known that a decrease of $\hat{B}$ causes efficient conversion from the
Poynting energy to kinetic energy in the context of the AGN jets that have a
collimated geometry (Camenzind, 1989; Begelman & Li, 1994).

##  5 Centrifugal drift current of Flow B

Our next question is whether the super fast flow such as Flow B can provide
enough centrifugal drift current to open magnetic field lies. We take the same
way as we derive (3), but we use $v_{\rm d}=\pm(c/e)(\gamma
mv_{\varphi}^{2}/R_{\rm Y}B_{\rm Y})$ and $nv_{\rm p}/B_{\rm p}=g(\psi)$ to
evaluate the density. Then we have

$\frac{\Delta}{R_{\rm Y}}=\frac{(v_{\rm p}/c)}{(v_{\varphi}/c)^{2}}\
\frac{\sigma_{0}}{\gamma}\ \frac{B_{\rm Y}}{B_{\rm p}}\ \frac{B_{\rm
p}^{\rm(out)}}{B_{\rm L}}.$ (15)

If the drift current is enough large, then $\Delta/R_{Y}$ must be much smaller
than unity. However, this expression indicates that $\Delta/R_{Y}$ will never
be small. The last factor $B_{\rm p}^{\rm(out)}/{B_{\rm L}}$ is about unity.
The factor ${B_{\rm Y}}/{B_{\rm p}}$ will not be much smaller than unity. The
effective field to produce the drift current is $z$-component of the poloidal
field. The region with a large azimuthal velocity would be where the flow
comes into the Y-point, and the ratio of the $z$-component to the poloidal
component of the flow will not be so small. Moreover, if the Y-point is
actually T-point as suggested by Contopoulos et al. (2024), it becomes unity.
The ratio $\sigma_{0}/\gamma$ can be small, if $\gamma\gg\sigma_{0}$. However,
a very efficient acceleration gives at most $\gamma\sim\sigma_{0}$. Since
$v_{\varphi}$ is something around $\sim 0.85$, the factor $(v_{\rm
p}/c)/(v_{\varphi}/c)^{2}$ will not be very small. In conclusion, accelerated
flows such as Flow B will not produce enough current to open the magnetic
field. This may sound strange when one compares with (3). However, as is seen
in the expression $n=g(\psi)B_{\rm p}/v_{\rm p}$, the plasma density decreases
due to the decrease of the magnetic field or increase of the cross section of
the flow. As a result, even if the centrifugal acceleration becomes efficient,
the drift current density does not increase.

##  6 Plasma injection at the Y-point

The reason why sufficient drift current to open magnetic field lines cannot be
obtained is due to decrease in flow density as the flow comes into the weak
field region around the Y-point. The magnetic field should be opened in the
centrifugal wind. One possible and very likely way to avoid the density
decrease is an injection of plasma at the Y-point. Let us consider a flow
injected from a top of the closed field region into the current layer, such as
Flow C in Figure 8.

The parameters characterize Flow C would be such that $\hat{B}$ stays small
throughout from injection (Figure 3), and that $\sigma_{0}$ is also small. In
the definition (9) of $\sigma_{0}$, the density parameter ${\cal M}_{\rm in}$
is calculated at the injection point where the magnetic field is weak around
the Y-point. Therefore, ${\cal M}_{\rm in}$ is much larger than a multiplicity
based on a typical Goldrech-Julian density, $\Omega_{*}B_{\rm L}/2\pi ce$.

For an example of Flow C, we suppose that the magnetic field strength at the
injection point is one hundredth of the field outside the current layer. Then,
$\sigma_{0}=10^{4}$ of Flow A and B is reduced to $\sigma_{0}=100$ for Flow C,
and $\hat{B}=\hat{B}_{\rm c}=\mbox{constant}=0.01$ for Flow C.

Figure 8: Flow C injected just within the Alfvén point, and Flow D is a flow
before the injection.

A family of the solutions is shown in Figure 9. We find another topology for
Flow C. Although $\hat{B}$ is constant, the topology is different from the
Michel’s radial-field solution. The two curves
$\partial\epsilon/\partial\xi=0$ and $\partial\epsilon/\partial v=0$ meet at
large distance, where an O-type critical point appears. On the other hand, the
X-type critical point seems to go to the Alfvén point as shown in a closeup
view (the right panel of Figure 9). The numerical solutions suggest that both
${\partial\epsilon}/{\partial v}=0$ and ${\partial\epsilon}/{\partial\xi}=0$
curves approach the Alfvén point tangentially to the $v$-axis. Thus, the
critical solution becomes the vertical line coinciding with the $v$-axis as
far as $v$ is not so large. Other solutions are two-valued functions of $\xi$,
and are inappropriate. When $\xi=1$ with $v\not=1$, the Bernoulli function
becomes

$\left(\frac{\epsilon(1,v)}{mc^{2}}\right)^{2}=\frac{1}{1-\lambda}\left[1-\left(\frac{1-\lambda}{\lambda}\right)^{2}\sigma^{2}v^{2}\right].$
(16)

It follows from this expression that in the limit $1-\lambda\rightarrow 0$, in
other words, the centrifugal-driven limit, the fast energy is $\epsilon_{\rm
f}=(1-\lambda)^{-1/2}$, regardless of $\sigma$. Figure 10 is drawn for
$1-\lambda=10^{-4}$, and therefore $\epsilon_{\rm f}=100$. In this case, the
injection point is always very closed to the Alfvén point, as we intend. The
topology of this type is not new but has been studied by Takahashi (1991), and
applied in the context of a disc wind around compact stars.

Figure 9: Left panel: The critical solution for Flow C injected just inside of
the Alfvén point. Right panel: A close-up view around the Alfvén point.

The Lorentz factor of the flow and the azimuthal velocity are shown in Figure
10. It can be seen that the flow energy is injected totally as kinetic energy.
Since $\lambda=\Omega\ell/\epsilon\approx 1$, angular momentum is brought in
as well.

Figure 10:  Properties of Flow C. Left panel: The Lorenz factor (solid line)
and the energy of the flow (dashed line) as functions of $\log\xi$. Right
panel: The azimuthal velocity as a function of $x$.

In contrast with Flow A and Flow B, we cannot determine the value of $\lambda$
and therefore $\gamma=(1-\lambda)^{-1/2}$ since the critical solutions are
degenerated. Determination of $\lambda$ is postponed until a self-consistent
treatment of plasma injection and determination of the poloidal field.

Pair creation in the vicinity of the Y-point is one of the injection process,
but a most efficient way would be magnetic reconnection. At the top of the
closed field region, a high density corotating plasma with a large Lorentz
factor accumulates, and is likely to be transferred to the open field region
via magnetic reconnection. In this case the value of $\sigma_{0}/\gamma$ in
(15) would be very small, meaning that the toroidal current is enough strong
to open the magnetic field, where note again that the face value of
$\sigma_{0}$ is reduced in Flow C. The reconnection process will be associated
with plasmoid emission (the dashed small loop in Figure 8). The injection may
be intermittent. We suppose that the plasmoid formation seen in PIC
simulations corresponds to this injection process.

Since the MHD equation of motion in the cold limit is equivalent to the
equation of motion of individual particles, it is likely that the motion of
plasmoids follows the solution of Flow C.

##  7 Pre-acceleration in the closed field region

Before the injection, the plasma that resides in the closed field region must
be accelerated in the azimuthal direction. In a steady state, such plasmas
migrate from the inner magnetosphere to the top of the closed region. This
type of flow can also be described by the field-aligned equation. Let us call
this Flow D. We again assume $\sigma_{0}=10^{4}$, and look for a solution with
the expected properties. The only parameter we can change is $\hat{B}$. We
find that the expected flow appears if $\hat{B}$ increases as $\xi\rightarrow
1$. It is know in the force-free model that the poloidal field strength
diverges as the Y-point is approached to the light cylinder (Lyubarskii, 1990;
Uzdensky, 2003; Timokhin, 2006). This increase of magnetic field plays a key
role in Flow D.

Figure 11: $\hat{B}$ (left) and the family of the solutions (right) of Flow D.

To give an illustrative solution for Flow D, we use the same model (14) for
$\hat{B}$, but we take a set of parameters, $D=10^{5}$, $x_{d}=0.9$, and
$\delta_{\rm d}=0.1$, mimicking a quick increase of $\hat{B}$ just before the
Alfvén point. We denote this model by $\hat{B}_{\rm d}(x)$ as shown in Figure
11. The right panel shows the family of the solutions for $\sigma_{0}=10^{4}$
and $1-\lambda=10^{-5}$. The value of $\lambda$ is chosen so that $\gamma=1$
at injection. As shown in Figure 12, the Lorentz factor grows up to
$\approx\hat{B}_{\rm d}(x_{\rm A})\approx 10^{5}$ as $\xi\rightarrow 1$, and
the azimuthal velocity follows the corotation. Since $\sigma_{0}/\gamma\ll 1$,
the toroidal current is large enough to open the magnetic field.

After some numerical works, we find that $\epsilon\approx\hat{B}_{\rm
d}(x_{\rm A})$, $1-\lambda\approx\epsilon^{-1}$, and
$\gamma\rightarrow\hat{B}_{\rm d}(x_{\rm A})$, provided that the injection
point is well within the light cylinder, the Lorentz factor is $\approx 1$
there, and the azimuthal velocity follows the corotation throughout the flow.

Although the increase of the poloidal field strength in the force-free model
is due to the imposed condition of the closed-to-open-field structure via an
electromagnetic force balance, we find that this increase is related to
increase of inertia and in turn the centrifugal drift current to open the
magnetic field. Therefore, a self-consistent treatment of the field-aligned
equation and the trans-field equation is expected to give the terminal Lorentz
factor and the saturated field strength at the top of the closed field region.
Furthermore, the transition from Flow D and Flow C would take place, and
therefore magnetic reconnection must be included in the self-consistent
treatment.

Figure 12: The Lorenz factor (solid line) and the energy of the flow (dashed
line) as functions of $\log\xi$, and the azimuthal velocity as a function of
$x$ for Flow D.

##  8 Discussion

Our results seem aligned well with the recent results of PIC simulations (Chen
& Beloborodov, 2014; Belyaev, 2015; Hakobyan et al., 2023; Bransgrove et al.,
2023; Hu & Beloborodov, 2022), in which kinetic energy dominant flows appear
around the equatorial current sheet, and plasmoids are ejected from the top of
the last closed field lines. More importantly, Hu & Beloborodov (2022) showed
that plasmas are accelerated in the azimuthal direction with high Lorentz
factors near the top of the closed field region, and that the plasmas are
released away freely by reconnection and plasmoid emission. These two steps
correspond to Flow D and Flow C, respectively. Comparison between their
results and our results would be a future issue.

A problem found in the previous sections is that our model does not predict
the thickness of the current sheet or the Lorentz factor to be reached. This
is because the fast critical point degenerates into the Alfvén point under the
circumstances of equipartition between the magnetic field and the plasma
energy densities ($\sigma\sim 1$). We need to undertake reconnection
processes. However, before that, we can make an estimate by assuming that the
thickness is of order of Larmor radius (Asano et al., 2004; Hoshino, 2020). We
employ (3), where we set $B_{\rm p}^{\rm(out)}\approx B_{\rm L}$, and the
density is the value in the Y-point, denoted by $n_{\rm Y}$, which is
distinguished from the density outside the current sheet, $n$. With the
condition $\Delta\sim(eB/\gamma mc^{2})^{-1}$, where we assume $B\approx
B_{\rm L}$ as a mean value of a strong field in the closed region and a weak
field in the open region, (3) yields

$1\approx\frac{2{\cal M}\gamma^{2}(mc^{2})^{2}}{(eB_{\rm L}R_{\rm
L})^{2}}\frac{B_{\rm L}}{B_{\rm Y}}\frac{n_{\rm Y}}{n}.$ (17)

The continuity implies $(B_{\rm L}/B_{\rm Y})(n_{\rm Y}/n)\approx 1$ for the
both sides of the Y-point. Then we have

$\gamma\approx\frac{\gamma_{\rm max}}{(2{\cal M})^{1/2}},\ \mbox{and
}\frac{\Delta}{R_{\rm L}}\approx\frac{1}{(2{\cal M})^{1/2}}.$ (18)

We also find that the drift velocity is well below $c$ by a factor of
$1/(2{\cal M})^{1/2}$. We suggest that the highest energy can be much higher
than a simple estimate of $\gamma_{\rm max}/{\cal M}$. We may hopefully verify
the efficiency factor ${\cal M}^{-1/2}$ by the future PIC simulations.

A significant difference between the present model and the previous PIC
simulations is the location of the Y-point. In our previous particle
simulations (see Figure 5 of Yuki & Shibata (2012)), we have implied that
field-aligned potential drop inside the light cylinder makes a superrotation,
and the Y-point move inward. The PIC simulation by Hu & Beloborodov (2022)
also shows supercorotation. In their simulation, there is field-aligned
acceleration around the null surface. The field-aligned acceleration,
superrotation, and subrotation are interact to each other to establish the
proper net loss of the angular momentum as shown in Shibata & Kisaka (2021).

Another possible reason for the Y-point within the light cylinder is that
$\lambda$ becomes much smaller than unity. This case can happen if the plasma
injected into the magnetosphere has high density and bring a large amount of
energy. The introduced kinetic and thermal energy contribute partly to drive
wind.

It is certain that the acceleration in the azimuthal direction takes place in
the closed field region, so that if ${\cal M}$ is not so large, where Larmor
radius becomes comparable with the light radius, then reconnection would start
within the light cylinder. This also can be a reason that the Y-point seems to
locate inside the light cylinder.

Contopoulos et al. (2024) discussed an increase of the poloidal magnetic field
just inside the Y-point. The increase would be caused by the close-open
structure imposed by hand in the force-free model. We have found that the
increase is essential for the acceleration in the azimuthal acceleration in
the sub-Alfvénic flow. Although the electromagnetic force balance with the
close-open structure requires the increase of the magnetic field, the increase
results in the particle acceleration, and the force balance will be achieved
with plasma inertia. In reality, the closed-open structure is maintained by
the intermittent plasmoid emission. The system is then sub-steady.

##  9 Conclusions

We investigate the centrifugal acceleration in an axisymmetric pulsar
magnetosphere under the ideal-MHD approximation. Anticipating the poloidal
field structure we have solved the field-aligned equations for different field
lines to obtain a hit toward a self-consistent view of the centrifugal wind.
The poloidal field structure comes into the flow model through $\hat{B}\propto
B_{\rm p}\varpi^{2}$.

Flows on the open field lines running apart form the equatorial current sheet,
called Flow A, are Poyinting energy dominant, and there is no efficient
centrifugal acceleration. Flow B that comes into the vicinity of a Y-point and
goes out in a current layer with finite thickness is characterized by decrease
of $\hat{B}$ just beyond the Alfvén point. Flow B becomes a super fast flow,
and the centrifugal acceleration takes place efficiently. However, the
centrifugal drift current is not enough to open the magnetic field lines. Flow
D is a flow coming from an inner part of the magnetosphere and reaches the top
of the closed field region. The key feature of Flow D is an increase of the
poloidal field strength toward the Alfvén point, the loop top. The azimuthal
velocity follows the corotaion, and the plasma gains a very large Lorentz
factor, $\gamma\sim\hat{B}(x_{\rm A})$. The corotating plasma with a large
Lorentz factor produced in Flow D is expected to be injected as Flow C via
magnetic reconnection. Flow C is a super fast flow, and goes out in the
equatorial current layer. The region around the junction of the two flows has
a centrifugal drift current large enough to open field lines. A jump from Flow
D to Flow C require magnetic reconnection and plasmoid emission.

Although we have not investigate the reconnection process, an simple estimate
suggests that the Lorentz factor of the centrifugal driven outflow goes up to
$\approx\gamma_{\rm max}{\cal M}^{-1/2}$. In a global point of view, the
centrifugal acceleration takes place in a small fraction of $\Delta/R_{\rm
L}\approx{\cal M}^{-1/2}$, and the most of the spin-down power is carried by
the Poynting energy in the open magnetic flux. The role of the centrifugal
force is to open the magnetic field and to establish the poloidal current
loop.

The centrifugal acceleration takes place in Flow B and Flow D just outside and
inside of the Alfvén point, respectively. The poloidal current lines cross the
stream lines there as indicated in in Figure 13. As a result, the current
sheets shrinks at the Y-point. This may be a favored poloidal current
distribution if the regularity of the poloidal field on the Alfvén surface
requires an infinitely thin current as seen in the force-free solution
(Contopoulos et al., 1999; Timokhin, 2006).

Figure 13: A schematic picture of the vicinity of the Y-point. The locations
of Flow A, B, C and D are indicated. The dashed curves are the streamlines,
and the thin curves are the poloidal current lines.

This work is supported by KAKENHI 22K03681 (SS, SK), 21H01078, and 22H01267
(SK). SS would like to thank Professor Masahiro Hoshino for fruitful
discussions on magnetic reconnection.

## References

  * Asano et al. (2004) Asano, Y., Mukai, T., Hoshino, M., et al. 2004, Journal of Geophysical Research (Space Physics), 109, A02212. doi:10.1029/2003JA010114
  * Begelman & Li (1994) Begelman, M. C. & Li, Z.-Y. 1994, ApJ, 426, 269. doi:10.1086/174061
  * Belyaev (2015) Belyaev, M. A. 2015, MNRAS, 449, 2759
  * Bransgrove et al. (2023) Bransgrove, A., Beloborodov, A. M., & Levin, Y. 2023, ApJ, 958, L9. doi:10.3847/2041-8213/ad0556
  * Camenzind (1986) Camenzind, M. 1986, A&A, 156, 137
  * Camenzind (1989) Camenzind, M. 1989, Accretion Disks and Magnetic Fields in Astrophysics, 156, 129. doi:10.1007/978-94-009-2401-7_14
  * Cerutti et al. (2015) Cerutti, B., Philippov, A., Parfrey, K., et al. 2015, MNRAS, 448, 606. doi:10.1093/mnras/stv042
  * Chen & Beloborodov (2014) Chen, A. Y., & Beloborodov, A. M. 2014, ApJ, 795, LL22. doi:10.1088/2041-8205/795/1/L22
  * Cheng et al. (1976) Cheng, A., Ruderman, M., & Sutherland, P. 1976, ApJ, 203, 209. doi:10.1086/154068
  * Contopoulos et al. (1999) Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351
  * Contopoulos et al. (2024) Contopoulos, I., Ntotsikas, D., & Gourgouliatos, K. N. 2024, MNRAS, 527, L127. doi:10.1093/mnrasl/slad153
  * Goldreich & Julian (1970) Goldreich, P. & Julian, W. H. 1970, ApJ, 160, 971. doi:10.1086/150486
  * Hakobyan et al. (2023) Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023, ApJ, 943, 105. doi:10.3847/1538-4357/acab05
  * Hoshino (2020) Hoshino, M. 2020, ApJ, 900, 66. doi:10.3847/1538-4357/aba59d
  * Hu & Beloborodov (2022) Hu, R. & Beloborodov, A. M. 2022, ApJ, 939, 42. doi:10.3847/1538-4357/ac961d
  * Jackson (1976) Jackson, E. A. 1976, ApJ, 206, 831. doi:10.1086/154446
  * Li & Melrose (1994) Li, J. & Melrose, D. B. 1994, MNRAS, 270, 687. doi:10.1093/mnras/270.3.687
  * Lyubarskii (1990) Lyubarskii, Y. E. 1990, Soviet Astronomy Letters, 16, 16
  * Michel (1969) Michel, F. C. 1969, ApJ, 158, 727. doi:10.1086/150233
  * Okamoto (1978) Okamoto, I. 1978, MNRAS, 185, 69. doi:10.1093/mnras/185.1.69
  * Ruderman & Sutherland (1975) Ruderman, M. A. & Sutherland, P. G. 1975, ApJ, 196, 51. doi:10.1086/153393
  * Shibata & Kisaka (2021) Shibata, S. & Kisaka, S. 2021, MNRAS, 507, 1055. doi:10.1093/mnras/stab2206
  * Shibata (1997) Shibata, S. 1997, MNRAS, 287, 262. doi:10.1093/mnras/287.2.262
  * Takahashi & Shibata (1998) Takahashi, M. & Shibata, S. 1998, PASJ, 50, 271. doi:10.1093/pasj/50.2.271
  * Takahashi (1991) Takahashi, M. 1991, PASJ, 43, 569
  * Takata et al. (2011) Takata, J., Wang, Y.,
  * Timokhin (2006) Timokhin, A. N. 2006, MNRAS, 368, 1055
  * Timokhin (2010) Timokhin, A. N. 2010, MNRAS, 408, 2092
  * Uzdensky (2003) Uzdensky, D. A. 2003, ApJ, 598, 446. doi:10.1086/378849
  * Yuki & Shibata (2012) Yuki, S., & Shibata, S. 2012, PASJ, 64, 43


42 2024 233

# Guiding drones by information gain

Alouette van Hove Corresponding Author. Department of Geosciences, University
of Oslo, Norway  Kristoffer Aalstad  Department of Geosciences, University of
Oslo, Norway  Norbert Pirk  Department of Geosciences, University of Oslo,
Norway

###### Abstract

The accurate estimation of locations and emission rates of gas sources is
crucial across various domains, including environmental monitoring and
greenhouse gas emission analysis. This study investigates two drone sampling
strategies for inferring source term parameters of gas plumes from atmospheric
measurements. Both strategies are guided by the goal of maximizing information
gain attained from observations at sequential locations. Our research compares
the myopic approach of infotaxis to a far-sighted navigation strategy trained
through deep reinforcement learning. We demonstrate the superior performance
of deep reinforcement learning over infotaxis in environments with non-
isotropic gas plumes.

##  1 Introduction

As a result of global warming, permafrost soils are thawing and will
potentially release large amounts of greenhouse gases such as methane and
carbon dioxide into the atmosphere [1]. Understanding the scale and
distribution of these greenhouse gas sources helps improve the reliability of
climate projections.

Determining the location and flux (i.e. emission rate) of greenhouse gas
sources remains challenging as they cannot be observed directly. Instead,
these source term parameters must be inferred through the analysis of
observational data, such as gas concentrations and other atmospheric
variables. This inverse problem is known as source term estimation (STE) [2].
The Bayesian inversion approach combines observational data with an
atmospheric plume model, accounting for uncertainties in both, to yield a
probability distribution over possible source term parameters. Sequentially
applying Bayes’ theorem to each observation updates the parameters’
probability distribution [3] so as to constrain uncertainty in the source term
estimates.

Observations can be acquired using static sensors (e.g. measurement towers),
mobile sensors (e.g. drones), or a combination of the two. Although mobile
sensors offer the advantage of flexibility and adaptability, the challenges of
devising the optimal sampling strategy become even more apparent. Sample
trajectories can be pre-planned in an offline manner, for example, by
selecting a sweeping “lawnmower” pattern, or they can be adapted to the
environment in an online manner [4]. The latter is known as informative path
planning (IPP).

In this paper, we explore two approaches to IPP for addressing the STE
problem. The first approach involves a modified version of infotaxis, an
online search strategy initially introduced by [5] [5] for source
localisation. Our work extends this algorithm to estimate multiple source term
parameters. The second approach uses deep reinforcement learning (DRL) [6] to
develop a sampling strategy. We train a mobile drone-based sensor to obtain
the most certain estimates of the source term parameters within a finite
number of observations, simulating the constraints of a drone’s limited
battery life. This study aims to determine whether a DRL strategy with an
information-based reward can outperform the infotaxis strategy in the context
of STE.

Recent studies have explored the application of DRL for guiding mobile sensors
in source seeking tasks [7, 9, 11, 8, 10]. However, these studies focused
solely on source localisation, neglecting source strength estimation.
Moreover, task completion in [7, 9, 11, 10] was determined based on the
proximity of the mobile sensor to the source. This is impractical for the
real-world identification of unknown greenhouse gas sources whose location is
uncertain. This study seeks to jointly address both the source localisation
and strength problem. While [12] [12] did consider the STE problem, they aimed
to minimize the number of observations needed for accurate STE, rather than
minimizing uncertainty within a single sampling flight. In line with [13]
[13], we apply an information-based reward function. Whereas this earlier
study estimated the source strength at a known location using tabular RL, we
apply deep RL herein to account for the larger parameter space.

The study’s primary contributions are as follows: (a) Application of DRL to
address the STE problem under pertinent constraints for greenhouse gas mapping
in the field - considering unknown source locations, fixed flight time due to
battery limitations, and noisy observations by the sensor. (b) Comparing the
performance of DRL with an information-based reward to infotaxis. (c)
Comparison between two neural network architectures: a feed forward structure
with fully connected layers (as utilized by [7, 9, 11, 8, 10]) and a
convolutional neural network.

##  2 Problem description

The navigation of mobile sensors for STE is a sequential decision-making
problem characterized by uncertainty and partial observability. The mobile
sensor cannot directly perceive the true state of the emissions but instead
receives incomplete and uncertain observations related to gas concentration.
In a mathematical framework, the problem can be represented as a Partially
Observable Markov Decision Process (POMDP) [14], that can be framed as a
belief Markov Decision Process (belief-MDP) where states are replaced by
belief states. In this context, the mobile sensor - _the agent_ \- maintains a
probability distribution over the unknown source term parameters - _its
belief_ \- while navigating a two-dimensional grid. At each time step, the
agent can choose from a set of five actions: {$\rightarrow$, $\leftarrow$,
$\downarrow$, $\uparrow$, stay}. After executing an action, the agent receives
a noisy observation which it uses to update its beliefs about the source term
parameters through Bayesian inference. Solving the navigation problem for STE
with mobile sensors means finding the optimal strategy - _the policy_ \- for
mapping belief states to actions. The agent’s objective is to maximize a
reward function that incentivizes it to select actions leading to more certain
source term estimations.

###  2.1 Observation model

A greenhouse gas source emits detectable particles that are dispersed through
domain $\mathcal{A}$ under turbulent transport conditions. We use the
advection-diffusion plume model presented in [5] to compute the mean number of
detected particles, or “hits”, by an agent at position
$\mathbf{x}=\left[x,y\right]\in\mathcal{A}$ during a sampling time interval
$\Delta t$. The source term parameter vector $\boldsymbol{\theta}$ includes
the source location
$\mathbf{x}_{\mathrm{s}}=\left[x_{\mathrm{s}},y_{\mathrm{s}}\right]\in\mathcal{A}$
and its source strength (i.e. flux) $\phi\in\mathbb{R}^{+}$. In two
dimensions, the mean number of hits detected by an agent with body radius $r$
is given by Eq. (1)

$\begin{split}\mu(\boldsymbol{\theta},\mathbf{x})&=\phi\Delta
t\frac{1}{\ln(\frac{\lambda}{r})}\exp\left(\frac{-(y-y_{\mathrm{s}})V}{2D}\right)\\\
&K_{0}\left(\frac{\lVert\mathbf{x}-\mathbf{x}_{\mathrm{s}}\rVert}{\lambda}\right)\,,\end{split}$
(1)

where the wind has been taken to blow in the negative $y$-direction. $K_{0}$
is the modified Bessel function of the second kind of order 0. In this model,
the gas particles have a finite lifetime $\tau$. Mean wind speed $V$ and
effective diffusivity $D$ (i.e. turbulent and molecular diffusivity) affect
the dispersion length scale $\lambda$ of the particles by Eq. (2)

$\lambda=\sqrt{\frac{D\tau}{1+\frac{V^{2}\tau}{4D}}}\,.$ (2)

An example of a mean hit map is shown on the left half of Figure 1. The
maximum number of hits is fixed to an upper bound, $h_{\textrm{max}}$. The
mean hit map becomes isotropic in the absence of wind.

The stochastic data generating process for the number of hits detected by the
sensor during sampling interval $\Delta t$ is modelled by a Poisson
distribution. The probability that the sensor detects $h$ hits is given by Eq.
(3)

$\Pr(h|\mu)=\frac{\mu^{h}}{h!}\exp(-\mu)\,.$ (3)

The right half of Figure 1 shows a possible measurement map. The difference
between the left half and right half of the map demonstrates the challenge of
STE in turbulent conditions. Observations are sparse and noisy, and gradient-
based sampling strategies for STE are therefore often not optimal and can fail
in turbulent flows.

Figure 1: Example of a hit map: (left half) map of the mean number of hits
$\mu$ calculated by the analytic plume model Eq. (1), and (right half) a
possible map of noisy sensor measurements $h$ from Eq. (3).

###  2.2 Bayesian inference

Each hit observation provides information about the source term parameters.
The agent applies Bayesian inference after each observation to maintain an up-
to-date belief about the source term parameters $\Pr(\boldsymbol{\theta})$,
where $\boldsymbol{\theta}$ is a discrete random variable of all possible
source locations and fluxes. In particular, Bayes’ rule (4) maps the prior
belief $\Pr(\boldsymbol{\theta}_{n})$ to the posterior belief
$\Pr(\boldsymbol{\theta}_{n}|h_{n})$ after a hit observation $h_{n}$ at time
$t_{n}$

$\Pr(\boldsymbol{\theta}_{n}|h_{n})=\frac{\Pr(h_{n}|\boldsymbol{\theta}_{n})\Pr(\boldsymbol{\theta}_{n})}{\Pr(h_{n})}\,,$
(4)

where the likelihood $\Pr(h_{n}|\boldsymbol{\theta}_{n})$ is computed by Eq.
(3) with $\mu$ given by Eq. (1). All probabilities are implicitly conditioned
on the plume model and background information, including agent’s location
$\mathbf{x}$ and wind speed $V$, among others. The evidence
$\Pr(h_{n})=\sum_{i}\Pr(h_{n}|\boldsymbol{\theta}_{n}^{(i)})\Pr(\boldsymbol{\theta}_{n}^{(i)})$
acts as a normalizing constant such that the posterior belief over all
possible $\boldsymbol{\theta}_{n}^{(i)}$ sums to one. Our agent applies Eq.
(4) sequentially in time by letting the posterior at time $t_{n}$ become the
prior at time $t_{n+1}$,
$\Pr(\boldsymbol{\theta}_{n+1})=\Pr(\boldsymbol{\theta}_{n}|h_{n})$, resulting
in dynamic belief updating.

##  3 Sampling strategies

Within the belief-MDP framework, the agent’s belief state is denoted as
$s=(\mathbf{x},\Pr(\boldsymbol{\theta}))$. In this study, both DRL and
infotaxis policies utilize a reward function based on information entropy.
Information entropy $H$ [15] is non-negative and measures the uncertainty of a
probability distribution. The information entropy (measured in nats) of a
belief state $s$ is given by Eq. (5)

$H(s)=-\sum_{i}\Pr(\boldsymbol{\theta}^{(i)})\ln\Pr(\boldsymbol{\theta}^{(i)})\,.$
(5)

Note that information entropy of a belief state is not directly affected by
the agent’s location $\mathbf{x}$, and relies only on
$\Pr(\boldsymbol{\theta})$, which in turn depends on the agent’s sampling
strategy and observations. The initial prior belief of the agent is a uniform
distribution over all possible source locations and fluxes, resulting in
maximum information entropy. Conversely, the information entropy of a
Kronecker delta distribution, where the probability is one for a specific
$\boldsymbol{\theta}^{(i)}$ and zero for all others, attains the minimum
information entropy of zero.

###  3.1 Infotaxis

Infotaxis [5] is a greedy policy which sequentially maximizes its reward
function. At each step, the agent selects an action $a$ that maximizes the
expected information gain $G$ [9], such that Eq. (6) follows

$\pi_{\textrm{infotaxis}}(s)=\arg\max_{a}G(s,a)\,.$ (6)

Expected information gain is the expected difference in information entropy
between the agent’s current belief state $s$ and all possible successor belief
states $s^{\prime}$, given by Eq. (7)

$G(s,a)=H(s)-\sum_{s^{\prime}}\Pr(s^{\prime}|s,a)H(s^{\prime})\,.$ (7)

The agent can reach different successor belief states depending on the hits it
receives after performing action $a$. Combining Eq. (6) and (7) gives Eq. (8)

$\pi_{\textrm{infotaxis}}(s)=\arg\min_{a}\sum_{s^{\prime}}\Pr(s^{\prime}|s,a)H(s^{\prime})\,.$
(8)

such that the policy can also be formulated as selecting the action that
transitions the agent to the successor belief state with the minimum expected
information entropy (i.e. minimum uncertainty).

###  3.2 Deep reinforcement learning

Infotaxis is a myopic approach, aiming to maximize the expected reward of the
next action. In contrast, RL trains a policy to maximize the cumulative reward
obtained throughout an entire _episode_ , considering a sequence of actions.
The essence of RL lies in identifying the optimal policy associated with the
optimal value function $v_{*}(s)$. This represents the maximum expected
_cumulative_ reward that can be achieved from a given belief state. The large
number of possible beliefs $\Pr(\boldsymbol{\theta})$ in a belief state $s$
prevents the use of exact methods for obtaining $v_{*}(s)$. Consequently, we
approximate $v_{*}(s)$ through a parameterized function
$\hat{v}(s,\mathbf{w})$, in the form of a deep neural network [16] with weight
vector $\mathbf{w}$. The Bellman optimality equation for the approximate value
function [6] is given by Eq. (9)

$\hat{v}(s,\mathbf{w}_{*})=\max_{a}\sum_{s^{\prime},r}\Pr(s^{\prime},r|s,a)\left[r+\hat{v}(s^{\prime},\mathbf{w}_{*})\right]\,,$
(9)

where $r$ is the reward. We use a reward function based on information gain
(7): $r=-H(s^{\prime})$. This reward function depends only on successor state
$s^{\prime}$, such that Eq. (10) follows as the Bellman optimality equation
for this study

$\hat{v}(s,\mathbf{w}_{*})=\min_{a}\sum_{s^{\prime}}\Pr(s^{\prime}|s,a)\left[H(s^{\prime})+\hat{v}(s^{\prime},\mathbf{w}_{*})\right]\,.$
(10)

Solving for the optimal weights $\mathbf{w_{*}}$ involves minimizing the
residual error, known as the Bellman optimality error (11) [9]

$\begin{split}L(\mathbf{w})=\mathbb{E}_{s}&\Biggl{[}\min_{a}\sum_{s^{\prime}}\Pr(s^{\prime}|s,a)\\\
&\left[H(s^{\prime})+\hat{v}(s^{\prime},\mathbf{w})\right]-\hat{v}(s,\mathbf{w})\Biggr{]}^{2}\,.\end{split}$
(11)

The expectation is taken over belief states $s$ visited when following policy
$\hat{\pi}$ derived from $\hat{v}$ and defined by Eq. (12)

$\hat{\pi}(s,\mathbf{w})=\arg\min_{a}\sum_{s^{\prime}}\Pr(s^{\prime}|s,a)\left[H(s^{\prime})+\hat{v}(s^{\prime},\mathbf{w})\right]\,.$
(12)

The functional $L(\mathbf{w})$ is the _loss function_ , and _training_ the
neural network entails optimizing the weights
$\mathbf{w}_{\star}=\arg\min_{\mathbf{w}}L(\mathbf{w})$ using stochastic
gradient descent, aiming to minimize the loss function. The subsequent DRL
policy is given by Eq. (12) evaluated with the optimized weights
$\pi_{\textrm{DRL}}(s,\mathbf{w}_{*})=\hat{\pi}(s,\mathbf{w}_{\star})$. Note
that by initializing the weights of the value function in proximity to zero,
the untrained DRL policy (12) closely resembles the infotaxis policy (8).

Our DRL algorithm is based on Deep Q-Networks (DQN) [17], an extension of
Q-learning to deep neural networks. While most RL algorithms are model-free
and use state-action values $q(s,a)$, our study is inspired by the approach of
[9] [9] and [11] [11] (with accompanying code [18]), and adopts a model-based
approach due to our exact knowledge of transition probabilities
$\Pr(s^{\prime},r|s,a)$. This allows us to work with $\hat{v}(s,\mathbf{w})$
instead of $\hat{q}(s,a,\mathbf{w})$.

Figure 2: Agent navigation in an environment with $\tilde{V}=2$ and
$\tilde{D}=2$, and source term parameters
$\tilde{\mathbf{x}}_{\mathrm{s}}=(9,1)$ and $\tilde{\phi}=2$: (left) DRL and
infotaxis sampling paths, and (right) (cumulative) entropy. Lower entropy
indicates a reduced uncertainty in the estimation of the source term
parameters.

We utilize a feedforward neural network and compare two different network
architectures: (a) an architecture with 3 fully connected hidden layers each
containing 32 units, and (b) an architecture with 4 hidden convolutional
layers with a filter size of 32 and a kernel size of 3. Average pooling with
kernel size 2 and stride 2 is applied. Both network architectures use
rectified linear units (ReLU) for activations, and a linear output layer. The
network sizes are tuned such that adding more weights does not increase their
performance any further. The weight parameters are trained over 20,000
episodes using stochastic gradient descent with mini-batches of size 128.

During preprocessing, we transform belief state
$s=(\mathbf{x},\Pr(\boldsymbol{\theta}))$ containing the agent’s position and
source term probability distribution ($N_{x}\times N_{y}\times N_{\phi}$
tensor), into a source term probability distribution centered on the agent
($(2N_{x}-1)\times(2N_{y}-1)\times N_{\phi}$ tensor) [18]. Tensor elements
outside the physical domain are set to zero. When used as an input for
architecture (b), the tensor can be conceived as an image representing spatial
data with dimensions $(2N_{x}-1)\times(2N_{y}-1)$ pixels and $N_{\phi}$
channels. The tensor is flattened when used as an input for architecture (a).

##  4 Results

We adopt a non-dimensional representation of observation model (1) whereby hit
counts are governed by dimensionless parameters: source flux $\tilde{\phi}$,
location $\tilde{\mathbf{x}}$, agent radius $\tilde{r}$, gas lifetime
$\tilde{\tau}$, mean wind speed $\tilde{V}$, and effective diffusivity
$\tilde{D}$. We assume that the source is located within a domain of 11 x 11
grid cells and has one of the following fluxes:
$\tilde{\phi}={\\{1,2,3,4,5\\}}$. As a result, we consider 605 possible source
term scenarios. The constants $\tilde{r}=0.5$ and $\tilde{\tau}=10,000,000$
are fixed. In sections 4.1 and 4.2, we maintain $\tilde{V}=2$ and
$\tilde{D}=2$, while other values for $\tilde{V}$ and $\tilde{D}$ are explored
in section 4.3. The agent starts in the center of the domain and performs 20
actions. The evaluation statistics presented herein are derived from 5,000
episodes.

###  4.1 Deep neural network

We assess the performance of a DRL policy approximated through a fully
connected (FC) architecture and a convolutional neural network (CNN) as
outlined in section 3.2. The FC architecture attains an average cumulative
entropy of $66~{}\textrm{nats}$, while the CNN performs slightly better with
an average of $65~{}\textrm{nats}$. In this study, we define _success_ as the
agent’s ability to estimate the correct source location and flux with a final
posterior probability of $\Pr\geq 0.5$. The FC architecture achieves success
in $70\%$ of episodes, while the CNN achieves this in $71\%$ of episodes. The
CNN architecture is used in subsequent sections of this study.

Figure 3: Relative Discrete Ranked Probability Score (DRPS) for source
location $\tilde{x}_{s}$ (top), $\tilde{y}_{s}$ (middle) and flux
$\tilde{\phi}$ (bottom). Circle radii correspond to the median (dark shade)
and 75${}^{\textrm{th}}$ percentile (light shade) relative DRPS. DRL results
are shown as (purple) left half circles and infotaxis results are shown as
(green) right half circles.

###  4.2 DRL versus infotaxis

We observe that the infotaxis policy achieves success in $58\%$ of episodes
(compared to 71% for DRL). Figure 2 shows a scenario with two agents
navigating an environment with a gas source with $\tilde{\phi}=2$ at
$\tilde{\mathbf{x}}_{\mathrm{s}}=(9,1)$: one agents follows the infotaxis
policy, while the other adheres to the DRL policy. In the absence of receiving
any hits, both agents move downwind (negative $y$-direction). The DRL agent
follows a sweeping pattern downwind, whereas the infotaxis agent traverses
along the edge of the domain. Upon detecting hits, the entropy of the agent’s
belief state reduces more significantly compared to zero hit events. Over the
entire episode, the DRL agent receives more valuable information than the
infotaxis agent, leading to a greater reduction in uncertainty of its belief
state, and subsequently a superior (i.e. smaller) cumulative entropy.

We investigate whether DRL outperforms infotaxis across all source term
scenarios considered in this study. As demonstrated in Table 1, the success
rate increases with source strength, consistently favoring the DRL approach
over infotaxis. By evaluating the performance across different source
locations (not shown), we find that the DRL policy surpasses the success rate
of infotaxis in $76\%$ of source locations, primarily in the upper and right
segments of the domain. Conversely, infotaxis outperforms DRL in $24\%$ of
source locations, mainly in the lower-left segment of the domain. As
demonstrated in Figure 2, this coincides with the region that the infotaxis
agent navigates towards in the absence of any detected hits.

Table 1: Success rate of DRL and infotaxis evaluated across different flux values $\tilde{\phi}$ (training encompassed all source term scenarios). | $\tilde{\phi}=1$ | $\tilde{\phi}=2$ | $\tilde{\phi}=3$ | $\tilde{\phi}=4$ | $\tilde{\phi}=5$
---|---|---|---|---|---
DRL | $40\%$ | $65\%$ | $80\%$ | $84\%$ | $85\%$
infotaxis | $31\%$ | $49\%$ | $65\%$ | $72\%$ | $76\%$

###  4.3 Sensitivity to $\tilde{V}$ and $\tilde{D}$

Figure 3 illustrates a relative Discrete Ranked Probability Score [19]
(DRPS${}_{\textrm{posterior}}$ / DRPS${}_{\textrm{prior}}$) for DRL and
infotaxis policies across environments with varying known diffusivities
$\tilde{D}$ and wind speeds $\tilde{V}$. The agents are trained under these
specific environmental conditions. DRPS is a measure for belief state
accuracy, favoring lower scores. Note that DRPS${}_{\textrm{prior}}$ is the
initial uniform belief state and remains consistent across all test cases. We
observe improved estimations of source location in conditions with higher
winds. In cases where $\tilde{V}>0$, DRL outperforms infotaxis. Additionally,
we observe that the estimation of source location on the crosswind axis
($\tilde{x}_{\mathrm{s}}$) surpasses that on the along windward axis
($\tilde{y}_{\mathrm{s}}$), and that lower diffusivities are beneficial for
STE in such cases. Under isotropic conditions ($\tilde{V}=0$), DRL and
infotaxis perform comparably given the specified training parameters.

##  5 Discussion

The findings presented in this study are a first effort in applying DRL with
an information-based reward to solve the STE problem with mobile gas sensors.
While infotaxis adopts a myoptic perspective, maximizing information gain
solely for the upcoming action, DRL has the capability to develop a far-
sighted strategy for action selection, considering information gain across the
entire sampling trajectory. This study demonstrates that adopting a DRL
approach has the potential to improve the accuracy of STE.

Our findings show that DRL outperforms infotaxis in non-isotropic plume
environments. The policy allows the agent to select actions that might not be
immediately optimal, but are chosen in anticipation of achieving better
rewards in the future. Comparative analysis shows that the DRL agent opts to
stay (downwind) centralized when its belief over the source term parameters
remains uncertain. In contrast, the infotaxis agent tends to move towards the
edge of the domain under similar circumstances. This trajectory could bring
the infotaxis agent closer to the source and valuable information, but also
poses a risk of leading the agent in the opposite direction, further away from
the source and information about the source term parameters.

DRL’s superiority over infotaxis does not extend to isotropic plumes in this
study, possibly due to the increased symmetry of the hit map and hence limited
directional information. In contrast, recent work by [9] [9] showed DRL’s
efficiency in source localisation within isotropic plumes. Key distinctions
between our studies are that [9] exclusively focuses on localisation, omitting
flux estimation. In their case, the flux strength was known which assisted the
agent in source seeking, whereas in our study the search for informative
observations is further complicated by the unknown source strength. Their
policy employed a $-1$ reward, terminating episodes upon the agent’s arrival
at the source - information not available to our agent. Moreover, episodes
initiated with the agent receiving a hit, while episodes can start with zero
hit observations in our study.

[9] [9] find that the performance of space-aware infotaxis is close to the
optimal policy obtained by DRL for source seeking. This variation on
traditional infotaxis incorporates an expected distance measure between agent
and source, and is designed to efficiently select the shortest path to the
source. Our approach differs: we aim to strategically select a 20-step path
that optimally constrains the source term parameters through observations. In
the environment studied currently, the optimal policy indeed leads the agent
towards the source. Nonetheless, unlike [9] [9], our agent entertains multiple
possible flux strengths, such that the agent may place equal probability to
the source being close by with a low flux strength as it being further away
with a higher flux strength. Furthermore, we deliberately refrain in this
study from rewarding actions that select the shortest path upfront, as our
long-term aim is to develop a framework that can perform in various
environments. In particular, future work will include domains with multiple
sources. In this case, we anticipate that the shortest path to one of the
sources is not the optimal policy to constrain the uncertainty over the
parameters characterising all sources.

##  6 Outlook

As we look ahead to the practical implementation of our DRL framework, we
identify several challenges and possibilities: (a) The framework presented in
this study is trained for known wind conditions and diffusion. Real-world
application requires training across a range of atmospheric conditions, which
can be inferred from observations by mobile and ground-based sensors, or a
combination of the two [20]. (b) The current observation model (1) assumes
steady-state atmospheric conditions, a simplification that may not hold in
real-world scenarios. More complex, time-dependent models, as employed by [20]
[20], may enhance the accuracy of the inference results. However, these models
are computationally expensive, posing challenges for online usage. (c) We
currently assume the presence of a single source within the domain. However,
in reality, each “grid cell” represents a possible greenhouse gas source or
sink. Future research aims to extend the framework to map heterogeneous
greenhouse gas fluxes across a domain. (d) Applying the framework in field
studies requires upscaling to larger domains and parameter spaces. Future work
should investigate whether the FC and CNN architectures in this study can
sustain their comparable performance within a realistic computational budget.

Despite the anticipated challenges in implementing the framework in real-world
scenarios, our study represents a significant advancement in the methodology
for STE of gas sources. It demonstrates encouraging outcomes achieved through
the guidance of drones by information gain strategies.

## Acknowledgements

This work was supported by the Research Council of Norway (project 301552
“Upscaling hotspots – understanding the variability of critical land-
atmosphere fluxes to strengthen climate models (Spot-On)” and project 333232
“Strategies for Circular Agriculture to reduce GHG emissions within and
between farming systems across an agro-ecological gradient (CircAgric-GHG)”).
This work is a contribution to the strategic research initiative LATICE
(Faculty of Mathematics and Natural Sciences, University of Oslo, project
UiO/GEO103920) as well as the Centre for Computational and Data Science
(dScience, University of Oslo).

## References

  * [1] E… Schuur, A.. McGuire, C. Schädel, G. Grosse, J.. Harden, D.. Hayes, G. Hugelius, C.. Koven, P. Kuhry, D.. Lawrence, S.. Natali, D. Olefeldt, V.. Romanovsky, K. Schaefer, M.. Turetsky, C.. Treat and J.. Vonk  “Climate change and the permafrost carbon feedback”  In _Nature_ , 2015  DOI: 10.1038/nature14338
  * [2] Michael Hutchinson, Hyondong Oh and Wen-Hua Chen  “A review of source term estimation methods for atmospheric dispersion events using static or mobile sensors”  In _Information Fusion_ , 2017  DOI: 10.1016/j.inffus.2016.11.010
  * [3] Simo Särkkä and Lennart Svensson  “Bayesian Filtering and Smoothing”  Cambridge University Press, 2023  DOI: 10.1017/9781108917407
  * [4] Adam Francis, Shuai Li, Christian Griffiths and Johann Sienz  “Gas source localization and mapping with mobile robots: A review”  In _Journal of Field Robotics_ , 2022  DOI: 10.1002/rob.22109
  * [5] Massimo Vergassola, Emmanuel Villermaux and Boris I. Shraiman  “‘Infotaxis’ as a strategy for searching without gradients”  In _Nature_ , 2007  DOI: 10.1038/nature05464
  * [6] R.. Sutton and A.. Barto  “Reinforcement Learning: An Introduction”  The MIT Press, 2018
  * [7] Bardienus P. Duisterhof, Srivatsan Krishnan, Jonathan J. Cruz, Colby R. Banbury, William Fu, Aleksandra Faust, Guido C… Croon and Vijay Janapa Reddi  “Learning to Seek: Autonomous Source Seeking with Deep Reinforcement Learning Onboard a Nano Drone Microcontroller”  In _ArXiv_ , 2019
  * [8] Lingxiao Wang, Shuo Pang and Jinlong Li  “Olfactory-Based Navigation via Model-Based Reinforcement Learning and Fuzzy Inference Methods”  In _IEEE Transactions on Fuzzy Systems_ , 2021  DOI: 10.1109/TFUZZ.2020.3011741
  * [9] Aurore Loisy and Christophe Eloy  “Searching for a source without gradients: how good is infotaxis and how to beat it”  In _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 2022  DOI: 10.1098/rspa.2022.0118
  * [10] Yong Zhao, Bin Chen, XiangHan Wang, Zhengqiu Zhu, Yiduo Wang, Guangquan Cheng, Rui Wang, Rongxiao Wang, Ming He and Yu Liu  “A deep reinforcement learning based searching method for source localization”  In _Information Sciences_ , 2022  DOI: 10.1016/j.ins.2021.12.041
  * [11] Aurore Loisy and Robin A. Heinonen  “Deep reinforcement learning for the olfactory search POMDP: a quantitative benchmark”  In _The European Physical Journal E_ , 2023  DOI: 10.1140/epje/s10189-023-00277-8
  * [12] Minkyu Park, Pawel Ladosz and Hyondong Oh  “Source Term Estimation Using Deep Reinforcement Learning With Gaussian Mixture Model Feature Extraction for Mobile Sensors”  In _IEEE Robotics and Automation Letters_ , 2022  DOI: 10.1109/LRA.2022.3184787
  * [13] A. van Hove, K. Aalstad and N. Pirk  “Using reinforcement learning to improve drone-based inference of greenhouse gas fluxes”  In _Nordic Machine Intelligence_ , 2023  DOI: 10.5617/nmi.9897
  * [14] Leslie Pack Kaelbling, Michael L. Littman and Anthony R. Cassandra  “Planning and acting in partially observable stochastic domains”  In _Artificial Intelligence_ , 1998  DOI: 10.1016/S0004-3702(98)00023-X
  * [15] David J C MacKay  “Information Theory, Inference, and Learning Algorithms”  Cambridge University Press, 2003
  * [16] Mykel J. Kochenderfer, Tim A. Wheeler and Kyle H. Wray  “Algorithms for decision making”  The MIT Press, 2022
  * [17] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A. Rusu, Joel Veness, Marc G. Bellemare, Alex Graves, Martin Riedmiller, Andreas K. Fidjeland, Georg Ostrovski, Stig Petersen, Charles Beattie, Amir Sadik, Ioannis Antonoglou, Helen King, Dharshan Kumaran, Daan Wierstra, Shane Legg and Demis Hassabis  “Human-level control through deep reinforcement learning”  In _Nature_ , 2015  DOI: 10.1038/nature14236
  * [18] A. Loisy and C. Eloy  “OTTO: A Python package to simulate, solve and visualize the source-tracking POMDP”  In _Journal of Open Source Software_ , 2022  DOI: 10.21105/joss.04266
  * [19] G. Candille and O. Talagrand  “Evaluation of probabilistic prediction systems for a scalar variable”  In _Quarterly Journal of the Royal Meteorological Society_ , 2005  DOI: 10.1256/qj.04.71
  * [20] N. Pirk, K. Aalstad, S. Westermann, A. Vatne, A. van Hove, L.. Tallaksen, M. Cassiani and G. Katul  “Inferring surface energy fluxes using drone data assimilation in large eddy simulations”  In _Atmospheric Measurement Techniques_ , 2022  DOI: 10.5194/amt-15-7293-2022


# Open data for Moroccan license plates for OCR applications : data
collection, labeling, and model construction.

Abdelkrim Alahyane 1 Mohamed El Fakir 1 Saad Benjelloun 1 Ikram Chairi 1

( 1 MSDA group, Mohammed VI Polytechnic University (UM6P), Lot 660, Hay Moulay
Rachid, Ben Guerir 43150 , Morocco)

###### Abstract

Significant number of researches have been developed recently around
intelligent system for traffic management, especially, OCR based license plate
recognition, as it is considered as a main step for any automatic traffic
management system. Good quality data sets are increasingly needed and produced
by the research community to improve the performance of those algorithms.
Furthermore, a special need of data is noted for countries having special
characters on their licence plates, like Morocco, where Arabic Alphabet is
used. In this work, we present a labeled open data set of circulation plates
taken in Morocco, for different type of vehicles, namely cars, trucks and
motorcycles. This data was collected manually and consists of 705 unique and
different images. Furthermore this data was labeled for plate segmentation and
for matriculation number OCR. Also, As we show in this paper, the data can be
enriched using data augmentation techniques to create training sets with few
thousands of images for different machine leaning and AI applications. We
present and compare a set of models built on this data. Also, we publish this
data as an open access data to encourage innovation and applications in the
field of OCR and image processing for traffic control and other applications
for transportation and heterogeneous vehicle management.

##  1 introduction

The world is knowing a continuously increase of road traffic since early
2000s. Traffic management is one of the main solutions to avoid congestion and
road accident. With the recent technological evolution concerning IoT and
smart sensors, many researches have addressed road traffic issues using data
generated by transportation and mobility systems[1].

One of the most developed sub-field of smart mobility, is Automatic Number
Plate Recognition (ANPR), which is an integration of Artificial Intelligence
along with computer vision and pattern recognition [2, 3]. ANPR consists of
three main stages: Number Plate Localization, character segmentation and
Optical character recognition (OCR).

OCR plays a main role in ANPR as it transforms character into encoded text
information. Different methods for OCR were proposed in the literature going
from Thresholding [4] to Machine Learning methods [3]. To improve performance
of recognition systems, especially for Machine Learning methods, algorithms
need to be trained on a big amount of data representing real situation of cars
in road traffic. For this reason, different platforms provide access to data
set for car plate images like Kaggle, Github or different open data web site.
However, as each country has its own way to characterize plate numbers, those
data bases can not be used for all countries, especially the ones using Arabic
alphabet like Morocco. Up to today, there is no data set containing Moroccan
vehicle plates in open data platforms.

In this work, we present a first Moroccan vehicle plate data set containing
705 unique labeled images for cars, motocycles and trucks. The data set was
collected by taking randomly pictures for vehicles in different context.We
present in this article also different way to manipulate this data set and to
possibility use data augmentation techniques to improve learning of Machine
Learning recognition algorithms. We present also baseline models for plate
segmentation and recognition using the built data set. The rest of the paper
is organized as follows. In section 2, we present the data set, in section 3
we focus on describing the labeling step and in section 4 we present baseline
models for the plate segmentation and OCR tasks. We conclude the article in
section 5.

##  2 Presentation of the data set

The natural first step in building plate segmentation and recognition models
is to collect an adequate data-set bringing together vehicle’s images with
relatively clear licence plates taken in different contexts and environments.
The goal is to capture all possible variations related to the surroundings of
real world plates and thus render the built AI models more robust and
efficient.

The set of images collected must at least respond to two major practical
issues. On the one hand, it must allow the construction of models ensuring the
segmentation of license plates in images of vehicles. On the other hand, it is
supposed to grant the detection and automatic recognition of the characters
constituting the plates by training suitable OCR models.

###  2.1 Data for license plate segmentation

The performance of models ensuring the plate segmentation task depends deeply
on the quantity and quality of the collected images used for training. As
mentioned above, multiple environments (views, backgrounds, scenes…) must be
taken into consideration when collecting the data to ensure that models can
properly generalize their learning. Therefore, we were interested in
collecting and building a set of good quality images with different contexts.
To achieve this, the data was gathered by randomly capturing pictures of
vehicles in Moroccan streets, the resulting data-set contains 705 unique
images of different types of vehicles (motorcycles, cars and trucks) with two
main categories of Moroccan license plates: square and rectangular plates.
Some examples of those different vehicle types are shown in figure 1, while
the number of images corresponding to each type of vehicle are specified in
table 1. Another important point to mention is that some of the images contain
multiple licence plates.

Figure 1: Sample of 8 images taken from the original data-set Table 1: Number of images aggregated by vehicle type. Vehicles types | Cars | motorcycles | Trucks
---|---|---|---
Number of images | 633 | 35 | 37

The shooting angles of the images as well as the zoom degrees are also
multiple and very diversified. This always aims to capture as much as
possible, the variability of conditions in which the trained model could
operate and predict.

As it will be detailed later, this data set was labeled in order to be used to
train machine learning algorithm for plates segmentation.

Finally we note that the data can be enriched using data augmentation
techniques such as adding noise, changing contrast, flipping and rotating, to
create training sets with four to five thousands (4000 to 5000) images, to
train robust machine leaning and AI applications.

###  2.2 Data for plates OCR

After collecting this data-set, rich enough in images of vehicles that can be
used to train a plate segmentation model, we cropped each image to extract
only the license plates. This aims to build a data-set of license plate images
which will be subsequently annotated to train OCR models allowing the
automatic segmentation and recognition of characters.

The plates obtained contains a set of 17 characters, 10 digits (from 0 to 9),
and 7 Arabic letters corresponding to the letters currently used in Moroccan
license plates. The distribution of these different labels is shown in the
histogram of figure 3. It should also be noted, as mentioned previously, that
the plates are in two main categories: the rectangular plates represented by
the plate (a) of figure 3, then the square plates, generally used for
motorcycles and some trucks, represented by the plate (b) of the same figure.

Figure 2: Labels distribution in the plate data-set.

Figure 3: License plate types.

###  2.3 Data preparation

After gathering the data-sets, the next step is to prepare the images to be
used to solve the two problems mentioned previously. Indeed, the construction
of high-performance models requires a significant amount of data, that’s why
data augmentation can be a good solution. This is a technique allowing a
significant increase in the quantity of images available to us. It is based on
the principle of generating new observations or images starting from the
existing ones. In our case, we performed the following data augmentation
operations:

  * •

Flipping the image vertically and horizontally.

  * •

Shearing, which means shifts part of the image like a parallelogram.

  * •

Adding noise, we choose to use a noise sampled from Gaussian distributions.

  * •

Adjust images brightness.

  * •

Apply affine transformations and scaling.

Data augmentation has demonstrated its effectiveness especially on the
training of neural networks which requires a significant amount of images.
After completing this step we obtained 4935 images.

####

We limit ourselves to these preparation operations and potential future users
of these data-sets can make all the appropriate pre-processing steps according
to the models that will be built and their respective architectures. Other
suggestions for effective preparation steps are as follows: Re-perform the
data augmentation operations to further increase not only the quantity of
images used for learning but also to add noise in order to improve
generalization, especially if the models used are neural networks and if the
computing power is encouraging. Some other approaches and algorithms can be
more precise by filtering the input images using noise reduction and smoothing
filters, for others it may be necessary to remove the background color from
the images to reduce the noise. If one wants to reduce the number of
parameters to estimate, thereby reducing the complexity of the model, it may
be more interesting to perform gray scaling operations or binarization.

In conclusion, the choice of data preparation operations depends on many
factors, including the nature of the models, their characteristics and the
available computing power.

##  3 Data Labeling

Labeling data is a major step in a supervised machine learning tasks, as the
model’s output depends on the labels we introduce during the learning phase.
The quality of the training data determines the quality of the model, the same
is true for the quality of annotation and its correctness.

There are several types of annotations for image segmentation. In our case we
use Bounding Boxes. Using an annotation tool we draw rectangular boxes around
the plates to segment them in the vehicle images. The traced boxes are then
represented by specific files containing the coordinates of the x and y axes
in the upper left corner and the coordinates of the x and y axes in the lower
right corner of the rectangle.

The data set we present is labeled in three formats, to give researchers and
users the freedom to use multiple models that require different annotation
formats. These formats are as following:

###  3.1 XML format

For each vehicle image, we generate an XML file to annotate the plates it
contains. The same process will be performed for the images of the plates used
for the OCR in order to segment the characters. This format allows us to
transfer a lot of useful information about the bounding boxes in the images
namely: Xmin,Ymin (upper left corner coordinates) and Xmax, Ymax (lower right
corner coordinates). The XML format also contains the name of the class
corresponding to each traced box. Furthermore, we can evaluate on a scale of 0
to 1 whether the object is difficult or truncated. For reasons of simplicity
we have chosen not to specify these last two parameters. An example of XML
content describing the annotation information for a single bounding box in the
image is shown in Figure 4.

Figure 4: The object tag of the XML file, it contains the labeling information
specific to a Bounding Box, including coordinates and label name.

###  3.2 Text file format

The second annotation format is a Text file adapted to YOLO [7] requirements
(normalized coordinates). This format is very simple, again we generate a text
file for each image, each line contains an object labeled in the image by
carrying the coordinates of the center of its boundary box, the height and the
width of the rectangular box and the detected class. An example of a line from
the text file representing an annotated object in the image is shown in figure
5.

Figure 5: Line of annotation text file with clarifications concerning the
various information it contains.

###  3.3 CSV format

The third format that we make available to researchers is the CSV format,
which carries like that of XML: the coordinates, the height and width of the
delimiting box, the class of the object and finally the name of the image. The
only exception is that here is a single CSV file concerning all the images,
with each line representing information concerning a single bounding box. An
example of the contents of the lines of the csv file is shown in figure 6.

Figure 6: Three (3) lines of the CSV file used to annotate the data-set.

##  4 Baseline models

The objective of this section is to present some models showing how one can
use the data set to build powerful AI applications. These models can also
serve as baseline to compare their performance with future models. Therefore,
we will build some machine and deep learning applications trying to solve the
two problems mentioned above, namely Plate segmentation and plate recognition
using OCR models.

###  4.1 Plate segmentation

Given its performance, speed and precision, the YOLOv3 [7] architecture is the
model chosen to be developed and trained on plate segmentation, we will use a
training-set of 650 images and a testing-test of 55 images. This set
corresponds to original images without data augmentation. Another baseline
model to consider is a pre-trained model on plate segmentation, which is the
Wpod-Net model [6]. It is a convolutional neural network model built from a
data-set containing images of foreign (i.e. non Moroccan) vehicle licence
plates. These vehicle images correspond to 10 different countries: Germany,
Vietnam, Japan, Thailand, Saudi Arabia, Russia, North Korea, USA, India and
China. The goal is to evaluate this model on our testing-set containing
vehicles with Moroccan licence plates, and use it as a baseline model to
compare these performances with those of YOLO v3.

####

We will use the accuracy metric on the testing-set as a performance measure to
identify the proportion of plates detected or segmented by each of the two
models. The results, as presented in table 2, show that the YOLO v3
architecture performs better than the Wpod-Net model despite the small number
of images that we used to train the YOLO neural network. Besides that, even
though the Wpod-Net is trained using a gigantic amount of vehicle images from
several countries and it has performed well in many experiments, it could not
outperform the YOLO model. This is explained by the fact that the Wpod-Net
model has never seen vehicle images of Moroccan license plates. So although
the plates have a lot of common characteristics, it was not able to detect all
of them. Unlike the YOLO model which, despite the small training set used,
succeeded in segmenting the overwhelming majority of the plates present in the
test-set.

Table 2: Results of used models. Models | Accuracy Training set | Accuracy Testing set
---|---|---
Wpod-Net | (pre-trained) | 0,7077
YOLOv3 | 0,9861 | 0,8506

####

This illustrates the importance of the data-set that we have grouped in the
Moroccan plate segmentation applications. It can be used, for example, to
perform transfer learning and improve already existing models that allow
segmentation of foreign plates.

###  4.2 Plate OCR

Concerning the recognition of the plates, we will use three different
approaches namely:

  * •

Histogram projection with k-NN classifier.

  * •

Histogram projection with CNN classifier.

  * •

YOLO v3, to detect characters in the plate and classify them simultaneously.

The thresholds of the histogram projection as well as the hyper-parameter $k$
are fixed by using the cross validation approach. On the other hand, the CNN
architecture used is MobileNetV2 by Google [5]. The training-set is composed
of characters from 650 plate images while the models are tested using
characters from 55 plate images. The results concerning the accuracy are shown
in table 3:

Table 3: Results of OCR used models. Models | Accuracy Training set | Accuracy Testing set
---|---|---
Histogram projection + k-NN classifier | 0,5128 | 0,2187
Histogram projection + CNN classifier | 0,6894 | 0,3775
YoloV3 | 0,9142 | 0,8238

We can therefore observe that the YOLO model again has the best performance.
It will then be a possible baseline architecture to compare with other models.

##  5 Conclusion

In this work we presented a Labeled data set of multiple-type vehicle images,
taken in Morocco. We presented the process of labeling this data for plate
segmentation and for plate’s number OCR using different formats. We put this
labeled data openly available for the research, startup and industrial
community to build models for applications related to traffic management or
car parks management in Morocco. We believe that initiatives such as ours can
catalyse the innovation and the technological development of AI solutions in
countries like Morocco by removing the burden linked to the non-availability
of open data and to the time-consuming tasks of collecting and labeling data.
We presented a set of Machine learning models that can be used as baseline
models and to which future users of this data set can compare and aim to
outperform by innovating in term of computational methods and algorithms. The
labeled data set, with data augmentation or without, can be downloaded at [8],
and all the implemented algorithms are available at [9].

## References

  * [1] I. Lana and J. Del Ser and M. Velez and E. I. Vlahogianni: Road Traffic Forecasting: Recent Advances and New Challenges. IEEE Intelligent Transportation Systems Magazine,10.1109/MITS.2018.2806634
  * [2] B. Tejas and D. Omkar and D. Rutuja and K. Prajakta and P. Bhakti: Number plate recognition and document verification using feature extraction OCR algorithm. At 2017 International Conference on Intelligent Computing and Control Systems (ICICCS), 10.1109/ICCONS.2017.8250683
  * [3] B. V. Kakani, D. Gandhi and S. Jani : Improved OCR based automatic vehicle number plate recognition using features trained neural network. At 8th International Conference on Computing, Communication and Networking Technologies (ICCCNT), Delhi, 2017, pp. 1-6, 10.1109/ICCCNT.2017.8203916
  * [4] Sushruth Shastry, Gunasheela G, Thejus Dutt, Vinay D S and Sudhir Rao Rupanagudi : i - A novel algorithm for Optical Character Recognition (OCR), IEEE 2013.10.1109/iMac4s.2013.6526442
  * [5] M. Sandler, A. Howard, M. Zhu, A. Zhmoginov and L. Chen : MobileNetV2: Inverted Residuals and Linear Bottlenecks, 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, Salt Lake City, UT, 2018, pp. 4510-4520, 10.1109/CVPR.2018.00474
  * [6] Montazzolli, Sergio and Jung, Claudio. (2018). License Plate Detection and Recognition in Unconstrained Scenarios.
  * [7] Redmon, Joseph and Farhadi, Ali. (2018). YOLOv3: An Incremental Improvement.
  * [8] MSDA open datasets webpage, https://msda.um6p.ma/msda_datasets. University Mohammed VI Polytechnic. Last accessed Dec 2020.
  * [9] Abdelkrim Alahyane, Mohamed Elfakir, Python Code for plate segmentation and plate numbers OCR, drive.google.com/drive/folders/1xbRTp9QXsKWlbYoqu2zFKQK3TjRXUHfs.