# Modelling and Search-Based Testing of Robot Controllers Using Enzymatic
Numerical P Systems

Radu Traian Bobe Florentin Ipate Ionuţ Mihai Niculescu Department of Computer
Science, Faculty of Mathematics and Computer Science, University of Bucharest,
Str Academiei 14, Bucharest, 010014, Romania<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The safety of the systems controlled by software is a very important area in a
digitalized society, as the number of automated processes is increasing. In
this paper, we present the results of testing the accuracy of different lane
keeping controllers for an educational robot. In our approach, the robot is
controlled using numerical P systems and enzymatic numerical P systems. For
tests generation, we used an open-source tool implementing a search-based
software testing approach.

###### keywords:

tests generation, numerical P systems, enzymatic numerical P systems, search-
based software testing, cyber-physical systems, membrane computing

##  1 Introduction

Due to the remarkable technological progress of late years, software
applications tend to have a considerable role in solving most problems of
everyday life. The medical, financial or automotive fields are just three of
the main areas in which software products are intensively used. Given the
importance of these areas in every individual’s life, ensuring product quality
and functionality is an essential step in the development process. The safety
of software systems for large-scale use is ensured by testing. Software
testing aims to validate the fulfillment of the requirements defined for the
developed product, as well as to identify possible unwanted behaviors
triggered by simulating certain operational contexts.

In this paper, we propose an approach for testing two different lane keeping
controllers designed to move an educational robot called E-puck [9]. Both
controllers are based on numerical P systems, introduced by G. Păun and R.
Păun in [12]. We also provide an equivalent version for the models using
enzymatic numerical P systems, an extension of numerical P systems, defined by
A. Pavel et al. in [10]. For this experiment we used some reliable tools which
will be introduced in the following sections.

The paper is structured as described: Section 2 presents the P system variants
to be used in the paper. Section 3 introduces the working environment,
including the tools used. Section 4 describes the models and the main
differences between them, while Section 5 illustrates the testing approach
along with the results. In the end, Section 6 presents the future work and
conclusions.

##  2 Preliminaries

Membrane computing is a field of research introduced by Gh. Păun in [11, 13].
The computational paradigm was originally inspired by the structure and
functionality of the living cells. Several classes of membrane systems (P
systems) have been later defined and investigated, being classified according
to the structure of the membranes as cell-like, tissue-like and neural-like P
systems. Membrane computing has made significant breakthroughs in the last
decades in fields like computer science, economics or biology. Depending on
the requirements, extensions of the main concept were introduced and our
experiment involves two of these types: numerical P systems and enzymatic
numerical P systems.

###  2.1 (Enzymatic) Numerical P System

The (enzymatic) numerical P systems [12, 10] are computational models that
only inherit the membrane structure from the membrane systems, more exactly a
cell-like membrane structure. The membranes contain variables and their values
are processed by the programs on every time unit. The whole system is
synchronized by a global clock in discrete time units.

The (enzymatic) numerical P system (EN P system) is defined by the tuple:

$\varPi\
\mathit{=(m,H,\mu,(Var_{1},Pr_{1},Var_{1}(0)),\dots,(Var_{m},Pr_{m},Var_{m}(0)))}$
(1)

where:

  * •

$\mathit{m}\geq 1$ is degree of the system $\Pi$ (the number of membranes);

  * •

$\mathit{H}$ is an alphabet of labels;

  * •

$\mu$ is membrane structure;

  * •

$\mathit{Var_{i}}$ is a set of variables from membrane $\mathit{i,1\leq i\leq
m}$;

  * •

$\mathit{Var_{i}(0)}$ is the initial values of the variables from region
$\mathit{i,1\leq i\leq m}$;

  * •

$\mathit{Pr_{i}}$ is the set of programs from membrane $\mathit{i,1\leq i\leq
m}$.

The program $\mathit{Pr_{l_{i},i}}$, $1\leq\mathit{l_{i}}\leq\mathit{m_{i}}$
has one of the two following forms:

    1. i)

non-enzymatic

$\mathit{F_{l_{i},i}(x_{1,i},\dots,x_{k,i})\rightarrow
c_{1,i}|v_{1}+c_{2,i}|v_{2}+\dots+c_{m_{i},i}|v_{m_{i}}}$

where $\mathit{F_{l_{i},i}(x_{1,i},\dots,x_{k,i})}$ is the production
function, $\mathit{c_{1,i}|v_{1}+c_{2,i}|v_{2}+\dots+c_{m_{i},i}|v_{m_{i}}}$
is the repartition protocol, and $\mathit{x_{1,i},\dots,x_{k,i}}$ are
variables from $Var_{i}$. Variables $v_{1},v_{2}\dots v_{m_{i}}$ can be from
the region where the programs are located, and to its upper and inner
compartments, for a particular region $i$. If a compartment contains more than
one program, only one will be chosen in non-deterministically manner.

    2. ii)

enzymatic

$\mathit{F_{l_{i},i}(x_{1,i},\dots,x_{k,i})|_{e_{i}}\rightarrow
c_{1,i}|v_{1}+c_{2,i}|v_{2}+\dots+c_{m_{i},i}|v_{m_{i}}}$

where $\mathit{e_{i}}$ is an enzymatic variable from $\mathit{Var_{i}}$,
$\mathit{e_{i}\notin\\{x_{1,i},\dots,x_{k,i},v_{1},\dots,v_{m_{i}}\\}}$. The
program can be applied at time $\mathit{t}$ only if
$\mathit{e_{i}>min(x_{1,i}(t),\dots,x_{k,i}(t))}$. The programs that meet this
condition in a region will be applied in parallel.

When the program is applied by the system at time $t\geq 0$, the computed
value

$\mathit{q_{l_{i},i}(t)=\frac{F_{l_{i},i}(x_{1,i}(t),\dots,x_{k,i}(t))}{\displaystyle\sum_{j=1}^{n_{i}}c_{j,i}}}$

represents the unitary portion that will be distributed to the variables
$\mathit{v_{1},\dots,v_{n}}$, proportional to coefficients
$\mathit{c_{1,i},\dots,c_{m_{i},i}}$, where
$\mathit{c_{j,i}\in\mathbf{N}^{+}}$ and the received values will be
$\mathit{q_{l_{i},i}(t)\cdot c_{1,i},\dots,q_{l_{i},i}(t)\cdot c_{m_{i},i}}$.

The values of variables, from $\mathit{t-1}$, present in the production
functions are consumed, reset to zero, and their new value is the sum of the
proportions distributed to variable through the repartition protocols, if they
appear in them or remain at the value zero.

##  3 Experimental environment

In this section we will provide brief descriptions of the tools we integrated
in our experiment. Firstly, we used an open-source software which allows the
simulation of numerical P systems and enzymatic numerical P systems. The
simulator is called PeP and will be introduced later in this section. Since we
don’t have the physical education robot involved in this study, we also used a
dedicated platform for robot simulations, called Webots. For tests generation,
we used a tool which won the SBST Tool Competition 2022 [5]. We will discuss
later in this section the arguments for using a search-based testing tool.

###  3.1 PeP simulator

PeP simulator [15] is an open-source product developed by A.Florea and C.Buiu,
used for simulations based on numerical P systems and enzymatic numerical P
systems. The program is written in Python and receives numerical P systems as
an input file. The input file includes the membrane structure and the contents
of each membrane, being stored in memory and executed.
PeP can be used as a stand-alone tool for simple simulations and run from the
command line with some options, like the number of simulations steps or a csv
document generation containing the values at each step of the simulation. As
observed in [15], the tool comes with a set of basic input files examples,
both numerical P systems and enzymatic numerical P systems.

Besides the simplicity of running this tool, another advantage which can be
taken into account when using PeP is that it can be used as an integrated
module in complex projects. We used this approach in our experiment in order
to make a controller accepted by the robot simulation platform and able to
receive information from the platform. In our lane keeping experiments, the
simulation ends when the robot drives off the generated lane or when the lane
is kept until the end.

###  3.2 Webots and E-puck

Webots is a robotics simulation software which allows the user to construct a
complex environment for programming, modelling and simulating mobile robots.
The environment can include multiple scene objects with different properties
which can be set from the graphic interface or from the generation files [8].
In addition, the robots can be equiped with a large number of objects called
nodes, like sensors, camera, GPS, LED, light sensor etc.
In our approach, additionally to the original equipments of E-puck, we used a
GPS attached to the turret slot in order to examine the coordinates at each
step of the simulation. The scenes, called ”worlds” in Webots, are containing
the road that the robot will try to follow. The roads are generated with
Ambiegen, a tool that will be described later in this section. Each world is
defined by a _.wbt_ file. The objects can be edited in this file and we used
this option in order to place the road object in the scene with coordinates
exported from Ambiegen. Additional functionalities, like sensors or GPS can
also be added to the robot by editing the world file which will be imported in
Webots, obtaining the visual representation of the scene.

As mentioned before, for this experiment we used a robot widely known from
educational and research purposes, called E-puck. At the moment, the robot has
some capabilities that are not implemented in Webots, but considering the fact
that both hardware and software components of E-puck are open source, this
remains a challenging opportunity [2].

E-puck has eight infrared proximity sensors placed around the body [9]. For
lane keeping simulation, we used just six of them: the two sensors placed in
front and the four placed two on each side. This aspect can be easily adapted
by changing the membrane structure and creating new membranes if more sensors
are needed or deleting a few of them if required. Each sensor has a
corresponding membrane in the numerical P system model and the association was
made in the controller. The robot has two motors attached to the body along
with two wheels, and the speed value is also changeable from the controller.

###  3.3 Ambiegen

Ambiegen is an open-source tool that utilizes evolutionary search for the
generation of test scenarios for autonomous systems. It can be used in
experiments involving lane keeping assist systems and robots navigating a room
with obstacles [7]. The software is developed in Python and uses evolutionary
search [17] for tests generation. The main goal of Ambiegen in this approach
is to generate roads as test cases in order to challenge E-puck to keep the
lane. The tool exports the roads in separate text files as a sequence of
points, representing the road spine. From this points, we can build the road
with a proportional size to E-puck.

Challenging different LKAS (Lane Keeping Assist Systems) involves a large
diversity of road topologies in order to detect the behavior in limit
situations, such as narrow curves. Ambiegen figures out the solution for
diversity by using a multi-objective genetic algorithm for search-based test
generation, called Non-dominated Sorting Genetic Algorithm-II (NSGA-II) [3].
In Ambiegen implementation, NSGA-II has two-objectives: to increase the fault
revealing power of test cases and to preserve their diversity [6]. This multi-
objective approach, which combines roads generation with a high attention to
diversity along with remarkable results at the competition mentioned above,
attracted our curiosity to integrate Ambiegen with Webots and testing E-puck
on the roads resulted.

###  3.4 Experimental Procedure

Considering the above information, we will detail the way we worked with the
presented tools. PeP and Ambiegen are developed in Python and so is the robot
controller.
First of all, we could easily integrate PeP with E-puck controller using the
PeP module which allowed us to parse the numerical P system model as an input
file for controller. Achieving this, the model membranes were associated with
controller variables. Names and constant values (e.g., robot cruise speed)
were taken using a text file containing membrane’s values of the variables.
The values were chosen empirically.

Next is a pseudocode version of the main loop in our controller, which
performs the simulation steps.

Algorithm 1 Simulation steps performing algorithm

1:repeat

2: for i=1 to number_of_sensors do

3: $\textit{sensor\\_membrane}(i)\leftarrow\textit{value}(i)$

4: run one simulation step

5: read lw, rw  from P system

6: $\textit{leftMotor}\leftarrow\textit{lw}$

7: $\textit{rightMotor}\leftarrow\textit{rw}$

8:until the end of the road or E-puck goes out of the road

Another challenge for us was to move the robot on the roads exported from
Ambiegen with the above presented approach. Ambiegen exports the roads as
.json files along with informations like test outcome, maximum curvature
coefficient etc. We took the road points from the file and wrote them in the
world file.

Webots provides the possibility to extend the set of scene nodes by adding
custom nodes created by users. The mechanism is called PROTO and is described
in [14]. After a node is extending with the PROTO interface, it can be
instantiated from the Webots graphic interface.

We used this technique to retrieve the points forming the spines of the roads
generated by Ambiegen and putting them into the wayPoints of the Road node.
Using javascript, used as scripting language by the PROTO, we constructed new
nodes illustrating the roads from Ambiegen. Then, in the graphic interface of
Webots, the road is represented in accordance with the road from Ambiegen.
With minimal Python code additions we plotted each generated road with the
corresponding spine to confirm that the shape illustrated in Webots respects
the original one.

##  4 Models

In this section we will present two models used to control the robot, the core
of the controller. The controller receives data from proximity sensors, that
measure distances to obstacles from the environment, to determine the
direction of movement of a differential two wheeled robot, E-puck, in our
case.

The proximity sensor has a range of 4 cm; if the obstacles are further than
this limit the sensor returns the value of 0. The proximity sensors are placed
on the left and right side of the robot in the direction of its movement at
different angles.

The first model was taken from [4] and adapted. The equations that calculate
the linear and angular velocity are shown below:

$\displaystyle\mathit{leftSpeed}=$
$\displaystyle\mathit{cruiseSpeed}+\displaystyle\sum_{i=1}^{n}\mathit{weightLeft_{i}}\cdot\mathit{prox_{i}}$
$\displaystyle\mathit{rightSpeed}=$
$\displaystyle\mathit{cruiseSpeed}+\displaystyle\sum_{i=1}^{n}\mathit{weightRight_{i}}\cdot\mathit{prox_{i}}$

The $\mathit{leftSpeed}$ and $\mathit{rightSpeed}$ are the speeds of the two
wheels of the robot. The enzymatic numerical P system described below
encapsulates this behavior.
The first model is defined as follows:

$\varPi\mathit{{}_{M_{1}}=(m,H,\mu,(Var_{1},Pr_{1},Var_{1}(0)),\dots,(Var_{m},Pr_{m},Var_{m}(0)))}$

where:

  * •

$\mathit{m=k\cdot 3+3,k=6}$, where $\mathit{k}$ is the number of proximity
sensors;

  * •

$\mathit{H=\\{s,s_{c}\\}\cup\displaystyle\bigcup_{i=1}^{k}\\{c_{i},s_{i},w_{i}\\}}$;

  * •

$\mathit{\mu=[[[]_{s_{1}}[]_{w_{1}}]_{c_{1}}\dots[[]_{s_{k}}[]_{w_{k}}]_{c_{k}}[]_{s_{c}}]_{s};}$

  * •

$\mathit{Var_{s}=\\{x_{s_{l}},x_{s_{r}}\\},Var_{s_{c}}=\\{x_{s_{c}}\\},}$
$\mathit{Var_{c_{i}}=\\{x_{c_{i},s_{l}},x_{c_{i},s_{r}},x_{c_{i},w_{l}},x_{c_{i},w_{r}},e_{c_{i}}\\},1\leq
i\leq k,}$
$\mathit{Var_{s_{i}}=\\{x_{s_{i},i}\\},1\leq i\leq k,}$
$\mathit{Var_{w_{i}}=\\{x_{w_{i},w_{l}},x_{w_{i},w_{r}},e_{w_{i}}\\},1\leq
i\leq}k$;

  * •

$\mathit{Var_{i}(0)=0,1\leq i\leq k}$;

  * •

$\mathit{Pr_{s}=\\{0\cdot x_{s_{l}}\cdot x_{s_{r}}\rightarrow
1|x_{s_{l}}+1|x_{s_{r}}\\}}$;
$\mathit{Pr_{s_{c}}=\\{3x_{s_{c}}\rightarrow
1|x_{s_{c}}+1|x_{s_{l}}+1|x_{s_{r}}\\}}$;
$\mathit{Pr_{c_{i}}=\\{x_{c_{i},s_{l}}\cdot
x_{c_{i},w_{l}}|_{e_{c_{i}}}\rightarrow 1|x_{s_{l}},}$
$\mathit{x_{c_{i},s_{r}}\cdot x_{c_{i},w_{r}}|_{e_{c_{i}}}\rightarrow
1|x_{s_{r}}\\},1\leq i\leq k}$;
$\mathit{Pr_{s_{i}}=\\{3x_{s_{i},i}\rightarrow
1|x_{s_{i},i}+1|x_{c_{i},s_{l}}+1|x_{c_{i},s_{r}}\\},1\leq i\leq k}$;
$\mathit{Pr_{w_{i}}=\\{2x_{w_{i},w_{l}}|_{e_{w_{i}}}\rightarrow
1|x_{w_{i},w_{l}}+1|x_{c_{i},w_{l}},}$
$\mathit{2x_{w_{i},w_{r}}|_{e_{w_{i}}}\rightarrow
1|x_{w_{i},w_{r}}+1|x_{c_{i},w_{r}}\\},1\leq i\leq k}$;

The meaning of the variables from the model is the following:

  * $\circ$

$\mathit{x_{s_{l}}}$ and $\mathit{x_{s_{r}}}$ from the region $\mathit{s}$
represent $\mathit{leftSpeed}$ and $\mathit{rightSpeed}$, the sum of the
products are accumulated in $\mathit{s}$ ;

  * $\circ$

$\mathit{x_{s_{c}}}$ from the compartment $\mathit{s_{c}}$ is
$\mathit{cruiseSpeed}$;

  * $\circ$

each pair of weights, $\mathit{weightLeft_{i}}$ and
$\mathit{weightRight_{i}}$, resides in the regions $\mathit{w_{i}}$,
$\mathit{1\leq i\leq k}$;

  * $\circ$

for each proximity sensor, $\mathit{prox_{i}}$, a compartment is defined,
namely $s_{i}$, containing a single variable, $\mathit{x_{s_{i},i}}$,
$\mathit{1\leq i\leq k}$;

  * $\circ$

the products are calculated by two distinct programs,
$\mathit{weightLeft_{i}\cdot prox_{i}}$, and $\mathit{weightRight_{i}\cdot
prox_{i}}$, $\mathit{1\leq i\leq k}$, in the compartments $\mathit{c_{i}}$.

The second model is an improvement on the first one. First we define the
function

$f(x)=\begin{cases}1,&\text{if }x=0\\\ 0,&\text{otherwise}\end{cases}$

This function will be used in the equations describing the behavior of the
model and in the production functions from the programs.

The equations describing the behavior are:

$\displaystyle\mathit{weightLeft}=$
$\displaystyle\displaystyle\sum_{i=1}^{n}\mathit{weightLeft_{i}}\cdot\mathit{prox_{i}}$
$\displaystyle\mathit{weightRight}=$
$\displaystyle\displaystyle\sum_{i=1}^{n}\mathit{weightRight_{i}}\cdot\mathit{prox_{i}}$
$\displaystyle\mathit{leftSpeed}=$
$\displaystyle\mathit{cruiseSpeed}\cdot\mathit{weightLeft+f(weightLeft)}\cdot\mathit{cruiseSpeed}$
$\displaystyle\mathit{rightSpeed}=$
$\displaystyle\mathit{cruiseSpeed}\cdot\mathit{weightRight+f(weightRight)}\cdot\mathit{cruiseSpeed}$

The model is defined as follows:

$\varPi\mathit{{}_{M_{2}}=(m,H,\mu,(Var_{1},Pr_{1},Var_{1}(0)),\dots,(Var_{m},Pr_{m},Var_{m}(0)))}$

where:

  * •

$\mathit{m=3k+3,k=6;}$

  * •

$\mathit{H=\\{s,w,s_{c}\\}\cup\displaystyle\bigcup_{i=1}^{k}\\{c_{i},s_{i},w_{i}\\}}$;

  * •

$\mathit{\mu=[[[[]_{s_{1}}[]_{w_{1}}]_{c_{1}}\dots[[]_{s_{k}}[]_{w_{k}}]_{c_{k}}[]_{s_{c}}]_{w}]_{s};}$

  * •

$\mathit{Var_{s}=\\{x_{s_{l}},x_{s_{r}}\\},Var_{w}=\\{x_{w_{l}},x_{w_{r}},e_{w}\\},Var_{s_{c}}=\\{x_{s_{c}}\\},}$
$\mathit{Var_{c_{i}}=\\{x_{c_{i},s_{l}},x_{c_{i},s_{r}},x_{c_{i},w_{l}},x_{c_{i},w_{r}},e_{c_{i}}\\},1\leq
i\leq k,}$
$\mathit{Var_{s_{i}}=\\{x_{s_{i},i}\\},1\leq i\leq k,}$
$\mathit{Var_{w_{i}}=\\{x_{w_{i},w_{l}},x_{w_{i},w_{r}},e_{w_{i}}\\},1\leq
i\leq k}$;

  * •

$\mathit{Var_{i}(0)=0,1\leq i\leq k}$;

  * •

$\mathit{Pr_{s}=\\{0\cdot x_{s_{l}}\cdot x_{s_{r}}\rightarrow
1|x_{s_{l}}+1|x_{s_{r}}\\}}$;
$\mathit{Pr_{w}=\\{x_{s_{c}}\cdot x_{w_{l}}+f(x_{w_{l}})\cdot
x_{s_{c}}|_{e_{w}}\rightarrow 1|x_{s_{l}},}$
$\mathit{x_{s_{c}}\cdot x_{w_{r}}+f(x_{w_{r}})\cdot
x_{s_{c}}|_{e_{w}}\rightarrow 1|x_{s_{r}}\\}}$;
$\mathit{Pr_{s_{c}}=\\{x_{s_{c}}\rightarrow 1|x_{s_{c}}\\}}$;
$\mathit{Pr_{c_{i}}=\\{x_{c_{i},s_{l}}\cdot
x_{c_{i},w_{l}}|_{e_{c_{i}}}\rightarrow 1|x_{s_{l}},}$
$\mathit{x_{c_{i},s_{r}}\cdot x_{c_{i},w_{r}}|_{e_{c_{i}}}\rightarrow
1|x_{s_{r}}\\},1\leq i\leq k}$;
$\mathit{Pr_{s_{i}}=\\{3x_{s_{i},i}\rightarrow
1|x_{s_{i},i}+1|x_{c_{i},s_{l}}+1|x_{c_{i},s_{r}}\\},1\leq i\leq k}$;
$\mathit{Pr_{w_{i}}=\\{2x_{w_{i},w_{l}}|_{e_{w_{i}}}\rightarrow
1|x_{w_{i},w_{l}}+1|x_{c_{i},w_{l}},}$
$\mathit{2x_{w_{i},w_{r}}|_{e_{w_{i}}}\rightarrow
1|x_{w_{i},w_{r}}+1|x_{c_{i},w_{r}}\\},1\leq i\leq k}$;

As we see from the definition of the $\varPi\mathit{{}_{M_{1}}}$ and its
equations, in the proximity of an obstacle, roadside, the speed of the wheel
on the side with the obstacle increases according to the weights. In certain
circumstances, particularly related to the geometry of a tight road curve, the
sum between travel speed and the weights results in a rotation in the opposite
direction of the obstacle, but a slight forward motion still continues,
causing the robot to leave the road.

The second model, $\varPi\mathit{{}_{M_{2}}}$, overcomes this problem by
calculating the product of the travel speed and the sum of the weights. In a
similar situation, the second model performs an angular rotation, in the
opposite direction to the obstacle, and when the proximity sensors no longer
detect obstacles, it continues moving forward with a constant velocity. This
behavior is modeled by compartment $\mathit{w}$.

In the next section, it can be seen that the $\varPi\mathit{{}_{M_{2}}}$ model
has better behavior in similar situations compared to the first model.

##  5 Simulation Results

Having presented the way we integrate the tools along with the models we will
detail in this section the testing stage. Ambiegen offers the possibility to
configure the setup for its genetic algorithms, choosing the values for
parameters like population size, number of generations, mutation rate or
crossover rate. Also, aspects like the amount of time allocated for tests
generation, map size, out of bound percent from which the test is considered
as failed can be set easily from the command line, when executing the main
Python file used in [5].

Test 1

Test 2

Test 3

Test 4

Figure 1: Illustration of different road tests

After series of trials we kept the default values, so we used the population
size 100 with 75 generations. Mutation rate is 0.4 and crossover rate is 1.
These values can be changed from the internal configuration file of Ambiegen.
From the command line we chose the time budget allocated to generation and
execution to be 1800 seconds and the map size to 200x200 meters which is the
default value. For the out of bound percentage we also kept the default value
(95%). After each simulation, Ambiegen exported roads spines coordinates as
text files and we also plotted each generated road, as mentioned before. Then
we could easily chose different roads based on the number of curves and their
angle as the main criteria to diversify the tests that were supposed to be
given to E-puck controller in Webots.

Next we report some roads along with the simulation results. The road is
marked with grey, whilst the trajectory using the first model is represented
with red. The trajectory of the second model (the improved one) is colored
with green.

In Figure 1 we can see the results of different roads tests. We observe that
Test 1, being the simplest test of the above presented, is passed by both
models. In the next scenarios, the complexity of the road increases and only
the improved model can pass.

From all the experiments, we noticed that the first model (the one marked with
red in results) cannot pass a huge majority of tests with curves, whilst the
improved model performs a rotation movement when the road is curved and for
this reason it manages to advance until the end of the road. Additional tests
were added to [16].

##  6 Conclusions and Future Work

In this paper we presented an approach to test different enzymatic numerical P
systems models using modern tools and search-based generated tests. We
evaluated our approach on a research and educational robot called E-puck,
virtually represented in Webots simulator. We set up our working environment
incorporating the tools with different scripts created to ensure a smooth
integration between them and also a better data processing. We formally
described each model involved and then showed the differences between the
lane-keeping controllers resulted when using each of them. As future work, we
will investigate the possibility to dynamically calculate the values for
weights, which are at the moment empirically assigned. Based on the controller
behavior during the previous tests, the weights values will be automatically
adapted. Also, we will try to develop a method to generate more complex roads
in order to better challenge the controllers.

## References

  * [1]
  * [2] Micael S Couceiro, Patricia A Vargas & Rui P Rocha (2014): _Bridging the reality gap between the Webots simulator and e-puck robots_.  Robotics and Autonomous Systems 62(10), pp. 1549–1567, 10.1016/j.robot.2014.05.007.
  * [3] Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal & TAMT Meyarivan (2002): _A fast and elitist multiobjective genetic algorithm: NSGA-II_.  IEEE transactions on evolutionary computation 6(2), pp. 182–197, 10.1109/4235.996017.
  * [4] Andrei George Florea & Cătălin Buiu (2017): _Modelling multi-robot interactions using a generic controller based on numerical P systems and ROS_.  In: 2017 9th International Conference on Electronics, Computers and Artificial Intelligence (ECAI), pp. 1–6, 10.1109/ECAI.2017.8166411.
  * [5] Alessio Gambi, Gunel Jahangirova, Vincenzo Riccio & Fiorella Zampetti (2022): _SBST tool competition 2022_.  In: Proceedings of the 15th Workshop on Search-Based Software Testing, pp. 25–32, 10.1145/3526072.3527538.
  * [6] Dmytro Humeniuk, Giuliano Antoniol & Foutse Khomh (2022): _AmbieGen tool at the SBST 2022 Tool Competition_.  In: Proceedings of the 15th Workshop on Search-Based Software Testing, pp. 43–46, 10.1145/3526072.3527531.
  * [7] Dmytro Humeniuk, Foutse Khomh & Giuliano Antoniol (2023): _AmbieGen: A Search-based Framework for Autonomous Systems Testing_.  arXiv preprint arXiv:2301.01234, 10.48550/arXiv.2301.01234.
  * [8] Olivier Michel (2004): _Cyberbotics ltd. webots™: professional mobile robot simulation_.  International Journal of Advanced Robotic Systems 1(1), p. 5, 10.5772/5618.
  * [9] Francesco Mondada, Michael Bonani, Xavier Raemy, James Pugh, Christopher Cianci, Adam Klaptocz, Stephane Magnenat, Jean-Christophe Zufferey, Dario Floreano & Alcherio Martinoli (2009): _The e-puck, a robot designed for education in engineering_.  In: Proceedings of the 9th conference on autonomous robot systems and competitions, IPCB: Instituto Politécnico de Castelo Branco, pp. 59–65.
  * [10] Ana Pavel, Octavian Arsene & Catalin Buiu (2010): _Enzymatic numerical P systems-a new class of membrane computing systems_.  In: 2010 IEEE Fifth International Conference on Bio-Inspired Computing: Theories and Applications (BIC-TA), IEEE, pp. 1331–1336, 10.1109/BICTA.2010.5645071.
  * [11] Gheorghe Păun (2002): _Membrane computing: an introduction_.  Springer Science & Business Media, 10.1007/978-3-642-56196-2.
  * [12] Gheorghe Păun & Radu Păun (2006): _Membrane computing and economics: Numerical P systems_.  Fundamenta Informaticae 73(1-2), pp. 213–227.  Available at http://content.iospress.com/articles/fundamenta-informaticae/fi73-1-2-20.
  * [13] Gheorghe Păun, Grzegorz Rozenberg & Arto Salomaa, editors (2010): _The Oxford Handbook of Membrane Computing_.  Oxford University Press.
  * [14] Webots reference manual https://cyberbotics.com/doc/reference/proto.
  * [15] A. G. Florea and C. Buiu, “PeP - an open-source software simulator of Numerical P systems and Numerical P systems with Enzymes,” 2017\. [Online]. Available: https://github.com/andrei91ro/pep.
  * [16] Github simulation results folder https://github.com/radubobe/Research/tree/main/Modelling%20and%20testing%20robot%20controllers%20using%20ENPS/Simulation%20results.
  * [17] Darrell Whitley, Soraya Rana, John Dzubera & Keith E Mathias (1996): _Evaluating evolutionary algorithms_.  Artificial intelligence 85(1-2), pp. 245–276, 10.1016/0004-3702(95)00124-7.


# Interpretation of the puzzling gamma-ray spectrum of the Geminga halo

Kun Fanga Xiao-Jun Bia,b aKey Laboratory of Particle Astrophysics, Institute
of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
bUniversity of Chinese Academy of Sciences, Beijing 100049, China

###### Abstract

The gamma-ray halo around Geminga is formed owing to the slow diffusion of the
electrons released by the Geminga PWN. The latest HAWC and HESS observations
exhibit complex features in the TeV gamma-ray spectrum of the Geminga halo. We
first show that the new results cannot be interpreted by the commonly used
simple model, where a single power-law injection spectrum and an energy index
of $\delta=1/3$ for the diffusion coefficient are assumed. We then propose a
two-population electron injection model based on the x-ray observations of the
Geminga PWN, which consists of a population of freshly accelerated electrons
escaping from the PWN through rapid outflows and a population trapped longer
inside the PWN before escaping. The two-population model interprets the HAWC
and HESS data well, and the goodness of fit improves significantly compared
with the single power-law injection model. It also predicts a different energy
dependency of the gamma-ray profile from the single power-law model, which
could be tested by LHAASO in the coming future. We note that a $\delta$
slightly larger than 1 is needed to fit the HAWC and HESS data consistently.
We also discuss the possible improvements by adopting the two-zone diffusion
model.

##  I Introduction

After electrons and positrons111Electrons will denote both electrons and
positrons hereafter if not specified. escape from some middle-aged pulsars (or
pulsar wind nebulae, PWNe hereafter), they diffuse very inefficiently in the
surrounding interstellar medium (ISM). The diffusion coefficient is several
hundred times smaller than the typical value in the Galaxy, and the
accumulated electrons generate observable gamma-ray halos through the inverse
Compton scattering (ICS) of the background photons, which are known as the
pulsar halos Sudoh  _et al._ (2019); Giacinti  _et al._ (2020). As the gamma-
ray morphology traces the spatial distribution of the parent electrons, pulsar
halos can be good indicators of electron propagation in localized regions of
the Galaxy.

The TeV halo around Geminga is the first discovered and so far best-studied
pulsar halo Abeysekara  _et al._ (2017). The unexpected slow-diffusion
environment indicated by the pulsar halo revives the discussion on the
possibility of Geminga as the source of the cosmic positron excess Hooper  _et
al._ (2017); Fang  _et al._ (2018); Profumo  _et al._ (2018); Tang and Piran
(2019); Xi  _et al._ (2019); Di Mauro  _et al._ (2019); Wang  _et al._ (2021).
However, there is no clear conclusion due to the lack of knowledge of the
electron injection spectrum and features of the slow-diffusion zone. In the
original paper of HAWC Abeysekara  _et al._ (2017), a simple model with a
single power-law electron injection spectrum and Kolmogorov’s energy
dependency of the diffusion coefficient is enough to interpret the
observation. We refer to this model as the “simple model” hereafter.

The recent gamma-ray spectrum measurements of HAWC and HESS provide more
detailed information about the Geminga halo Zhou (2021); Torres Escobedo  _et
al._ (2021); Mitchell  _et al._ (2021). The latest HAWC spectrum indicates a
possible bump feature around 10 TeV Torres Escobedo  _et al._ (2021).
Meanwhile, the spectrum unexpectedly climbs again below $\approx 3$ TeV. This
low-energy feature is also confirmed by the observation of HESS, which
measures the gamma-ray spectrum within $None$ around Geminga Mitchell  _et
al._ (2021). Given the high angular resolution of HESS, the low-energy gamma-
ray component is very likely associated with Geminga. We will show that these
new features can no longer be interpreted by the “simple model.”

In Sec. II, we briefly introduce the calculation of the gamma-ray spectrum of
the Geminga halo and present the difficulties of interpreting the new spectral
features with the “simple model.” Considering the electron injection process
from the Geminga PWN implied by the x-ray observations, we propose a two
population injection model in Sec. III to better fit the new measurements of
the gamma-ray spectrum. In Sec. IV, we further discuss the effect of two-zone
diffusion. Finally, we conclude in Sec. V.

##  II Difficulties of the “simple model”

Firstly, we introduce the routine calculation of the gamma-ray spectrum of
pulsar halos. We obtain the electron number density around the pulsar by
solving the propagation equation and then do the line-of-sight integration to
get the electron surface density. The gamma-ray spectrum is derived from the
electron surface density through the ICS calculation.

The electron propagation equation can be expressed by

$\frac{\partial N(E_{e},\textbf{\em r},t)}{\partial
t}=\nabla\cdot[D(E_{e})\nabla N(E_{e},\textbf{\em
r},t)]+\frac{\partial[b(E_{e})N(E_{e},\textbf{\em r},t)]}{\partial
E_{e}}+Q(E_{e},\textbf{\em r},t)\,,$ (1)

where $N$ is the electron number density and $E_{e}$ is the electron energy.
The diffusion coefficient takes the form of $D(E_{e})=D_{0}(E_{e}/{\rm
100~{}TeV})^{\delta}$, where we assume $D_{0}=3.2\times 10^{27}$ cm2 s-1 in
the present work as measured by HAWC Abeysekara  _et al._ (2017). The second
and third terms on the right-hand side are the energy-loss and source terms,
respectively. Synchrotron radiation and ICS dominate the energy losses of
high-energy electrons. We take the local magnetic field strength (3 $\mu$G,
Minter and Spangler (1996)) for the synchrotron component. We adopt the method
introduced in Ref. Fang  _et al._ (2021a) and the seed photon field given in
Ref. Abeysekara  _et al._ (2017) to get the ICS component.

The source function takes the form of

$Q(E_{e},\textbf{\em r},t)=\left\\{\begin{aligned}
&q(E_{e})\,\delta(\textbf{\em r}-\textbf{\em r}_{s})\,[(t_{s}+t_{\rm
sd})/(t+t_{\rm sd})]^{2}\,,&t\geq 0\\\ &0\,,&t<0\end{aligned}\right.\,,$ (2)

where $q(E_{e})$ is the current electron injection spectrum, $\textbf{\em
r}_{s}$ and $t_{s}$ are the position and age of Geminga, respectively, and
$t_{\rm sd}$ is the pulsar spin-down time scale, which is set to be 10 kyr.
The time profile of the source function is assumed to follow the pulsar spin-
down luminosity, and $t=0$ corresponds to the birth time of Geminga.

The solution of Eq. (1) can be expressed by

$N(E_{e},\textbf{\em r},t)=\int_{R^{3}}d^{3}\textbf{\em r}_{0}\int_{t_{\rm
ini}}^{t}dt_{0}\,\frac{b(E_{e}^{\star})}{b(E_{e})}\frac{1}{(\pi\lambda^{2})^{3/2}}{\rm
exp}\left[-\frac{(\textbf{\em r}-\textbf{\em
r}_{0})^{2}}{\lambda^{2}}\right]\,Q(E_{e}^{\star},\textbf{\em
r}_{0},t_{0})\,,$ (3)

where

$E_{e}^{\star}\approx\frac{E_{e}}{[1-b_{0}E_{e}(t-t_{0})]}\,,\quad\lambda^{2}=4\int_{E_{e}}^{E_{e}^{\star}}\frac{D(E^{\prime}_{e})}{b(E^{\prime}_{e})}dE^{\prime}_{e}\,,$
(4)

and $t_{\rm ini}={\rm max}\\{t-1/(b_{0}E_{e}),0\\}$. We integrate $N$ over the
line of sight from Earth to the vicinity of the pulsar and get the electron
surface density $S_{e}(\theta)=\int_{0}^{\infty}N(l_{\theta})dl_{\theta}$,
where $\theta$ is the angle observed away from the pulsar, $l_{\theta}$ is the
length in that direction, and $N(l_{\theta})$ is the electron number density
at a distance of $\sqrt{d^{2}+l_{\theta}^{2}-2dl_{\theta}\cos\theta}$ from the
pulsar, where $d$ is the distance between the pulsar and Earth. The gamma-ray
surface brightness $S_{\gamma}(\theta,E_{\gamma})$ is derived from the
electron number density and the standard calculation of ICS Blumenthal and
Gould (1970). Finally, we can get the gamma-ray spectrum
$F_{\theta_{0}}(E_{\gamma})$ within an arbitrary angular radius $\theta_{0}$
around the pulsar by
$F_{\theta_{0}}(E_{\gamma})=\int_{$$}^{\theta_{0}}S_{\gamma}(\theta,E_{\gamma})2\pi\theta
d\theta$.

Figure 1: Left: the gamma-ray spectra of the Geminga halo given by the ”simple
model” (a single power-law injection spectrum and Kolmogorov’s energy index of
the diffusion coefficient for electrons), compared with the HAWC Torres
Escobedo  _et al._ (2021) and HESS Mitchell  _et al._ (2021) measurements.
Right: the same as the left, while the energy index of the diffusion
coefficient, $\delta$, is set as a free parameter to fit the HAWC and HESS
data. The best-fit $\delta$ is 1.15.

For the “simple model,” the injection spectrum takes the form of
$q(E_{e})\propto E_{e}^{-p}$. The normalization of the injection spectrum can
be obtained by the relation of $\int_{\rm
1GeV}^{\infty}q(E_{e})E_{e}dE_{e}=\eta L$, where $L$ is the current pulsar
spin-down luminosity, and $\eta$ is the conversion efficiency from the spin-
down energy to the electron energy. The pulsar age, distance, and spin-down
luminosity are 342 kyr, 250 pc, and $3.25\times 10^{34}$ erg s-1, which can be
found in the Australia Telescope National Facility catalog Manchester  _et
al._ (2005). We set $p$ and $\eta$ as the free parameters in the following
fit. The energy index of the diffusion coefficient, $\delta$, is set to be
$1/3$ as predicted by Kolmogorov’s theory.

We fit the “simple model” to the HAWC gamma-ray spectrum by the minimum chi-
square method and then compare the best-fit model with the HESS spectrum. The
best-fit parameters are $p=2.39$ and $\eta=5.46$222The required conversion
efficiency is larger than 1, while it is not a serious problem for a power-law
injection spectrum with $p>2.0$. In this case, the energy of the electron
spectrum is concentrated at GeV energy range, where the injection spectrum is
not constrained. A low-energy spectral break or cutoff can avoid a too large
$\eta$., respectively. The result is presented in the left panel of Fig. 1.
There are two serious problems in the result. Firstly, the best-fit model to
the HAWC data fits poorly to the HESS spectrum, especially for the low-energy
data. It means that the HAWC and HESS results cannot be consistently
interpreted with the “simple model.” The gamma-ray flux ratio $F_{$$}/F_{$$}$
predicted by the “simple model” increases with energy, while the measurements
show a contrary tendency. Secondly, the fit to the HAWC spectrum alone is also
poor—the reduced $\chi^{2}$ is $\approx 3.5$. A single power-law injection
spectrum cannot reproduce the possible spectral bump and high-energy cutoff
features.

The small $F_{$$}/F_{$$}$ below $\approx 3$ TeV indicates that the spatial
distribution of low-energy electrons could be more concentrated than expected.
A $\delta$ larger than the Kolmogorov value is the most straightforward
solution to this problem. We add $\delta$ as a free parameter and
simultaneously fit the HAWC and HESS spectrum. The best-fit parameters are
$p=2.39$, $\eta=5.21$, and $\delta=1.15$. The required $\delta$ is
significantly larger than $1/3$ ($1/2$) as predicted by Kolmogorov’s
(Kraichnan’s) theory and close to the value of Bohm diffusion. As shown in the
right panel of Fig. 1, a consistent interpretation of the two spectra seems to
be achieved. However, the reduced $\chi^{2}$ of the fit is 1.85, and the model
is rejected by the goodness-of-fit test at a confidence level of $97\%$. The
main reason is still the poor fit to the HAWC spectrum.

##  III A two-population injection model

###  III.1 Description of the model

Electron-positron pairs produced by pulsars are accelerated to very high
energies in the PWNe and then released into the ISM to generate pulsar halos.
Thus, observations of the non-thermal radiation of PWNe are helpful for
inferring the electron injection spectrum. The Geminga PWN has been observed
in x-rays by XMM-Newton and Chandra telescopes with high spatial resolution
Caraveo  _et al._ (2003); Pavlov  _et al._ (2006, 2010); Posselt  _et al._
(2017). The most distinct feature of the PWN is the three-tail structure,
consisting of two lateral tails and an axial tail. The latest Chandra
observation indicates that the lateral tails, which directly connect to the
pulsar, are more likely to be interpreted by outflows induced by polar jets of
the pulsar than the limb-brightened shell of the PWN Posselt  _et al._ (2017);
Hui  _et al._ (2017). The axial tail could be interpreted as a crushed torus
Posselt  _et al._ (2017) or the bow-shock nebula Hui  _et al._ (2017).

According to the image and spectral measurements of the Geminga PWN, we
propose a two-population electron injection model and illustrate it with the
sketch in the left of Fig. 2. Obviously, electrons can be injected through the
lateral outflows, and we name this population Pop A. As there is no evidence
of synchrotron cooling in the lateral tails, the velocity of the outflows must
be larger than $\sim 1000$ km s-1 Posselt  _et al._ (2017). Moreover, Ref. Hui
_et al._ (2017) suggests that the outflows are mild relativistic. The observed
brightness difference between the northern and southern tails could be
interpreted by Doppler boosting of the high-speed flows. All these indicate
that Pop A is fresh electrons efficiently escaping from the acceleration
sites, which is consistent with the very hard x-ray spectrum of the lateral
tails. Besides, the energy of the parent electrons of the gamma-ray halo can
reach $\sim 100$ TeV. The synchrotron lifetime of 100 TeV electrons is
$\approx 200$ yr considering the $20$ $\mu$G magnetic field in the Geminga PWN
Posselt  _et al._ (2017), which means that part of the high-energy electrons
must escape from the PWN within 200 yr after being accelerated. Therefore, Pop
A is expected to explain the high-energy part of the halo spectrum.

We assume a super-exponentially-cutoff power law for the injection spectrum of
Pop A:

$q_{A}(E_{e})\propto E_{e}^{-p}\,{\rm
exp}\left[-\left(\frac{E_{e}}{E_{c}}\right)^{2}\right]\,.$ (5)

The electron spectral indices of the northern and southern tails of the x-ray
PWN are 0.67 and 1.04, respectively Posselt  _et al._ (2017), and we take the
average value of 0.85 for $p$. The cutoff term in Eq. (5) describes the
acceleration limit of the PWN, the form of which is suggested by Ref.
Zirakashvili and Aharonian (2007). The cutoff energy $E_{c}$ is set as a free
parameter. Another free parameter of Pop A is the conversion efficiency
$\eta_{A}$, which corresponds to the normalization of the injection spectrum.

Meanwhile, electrons may also escape from the axial tail in a less efficient
way, and we name this population Pop B. The x-ray observations show a tendency
of spectral softening along the axial tail Posselt  _et al._ (2017), which may
be due to the synchrotron cooling of electrons. This implies that Pop B may be
trapped longer in the PWN than Pop A before being injected into the ISM.
Suppose the trapping time is $\sim 1000$ yr, a spectral break at $\sim 10$ TeV
is expected for the electron injection spectrum. We expect that Pop B
dominates the low-energy part of the gamma-ray spectrum.

Compared with Eq. (5), a spectral break is added to the injection spectrum of
Pop B:

$q_{B}(E_{e})\propto
E_{e}^{-p}\,\left[1+\left(\frac{E_{e}}{E_{b}}\right)^{s}\right]^{\Delta
p/s}\,{\rm exp}\left[-\left(\frac{E_{e}}{E_{c}}\right)^{2}\right]\,,$ (6)

where the break energy $E_{b}$ is set as a free parameter. The smooth
parameter $s$ has little effect on the result and is set to be $5$. As the
injection spectrum of Pop A is the spectrum of freshly accelerated electrons,
we assume that Pop B shares the same $p$ and $E_{c}$ with Pop A. The energy
range of x-ray observations is 0.3-8 keV, corresponding to the electron energy
of $\approx 25-130$ TeV. As this energy range is between $E_{b}$ and $E_{c}$
(we will show below that $E_{c}$ is larger than 130 TeV as required by the
fit), we may estimate $-p+\Delta p$ by the x-ray spectral index of the axial
tail. We adopt the index of the region relatively far away from the pulsar
(the A2$+$A3 region333In Ref. Posselt  _et al._ (2017), the A4 region is the
farthest axial region to the pulsar, while the spectral measurement of the A4
region may be contaminated by a known star. shown in Ref. Posselt  _et al._
(2017)) and get $\Delta p=-2.85$. This electron spectrum is steeper than that
predicted by a constant-injection-cooling scenario and flatter than that of a
pure cooling scenario, which may be ascribed to the complex particle transport
in the bow-shock nebula. In addition, the conversion efficiency for Pop B,
$\eta_{B}$, is set to be a free parameter.

Figure 2: Left: sketch of the two-population injection model based on the
morphology of the Geminga x-ray PWN Posselt  _et al._ (2017). Right: best-fit
gamma-ray spectra to the HAWC and HESS data with the two-population injection
model. For each population, the spectrum within $None$ ($None$) field around
the pulsar is drawn with solid (dotted) line.

###  III.2 Fitting result and discussion

In summary, the free parameters of the two-population model are $\eta_{A}$,
$\eta_{B}$, $E_{c}$, $E_{b}$, and $\delta$. We fit the two-population model to
both the HAWC and HESS data and present the fitting result in the right panel
of Fig. 2. The reduced $\chi^{2}$ of the best-fit result is 1.10, which
improves significantly compared with the one-population models in Sec. II. The
high-energy cutoff term of Pop A can well reproduce the high-energy HAWC data,
and the hard power-law term of Pop A gives a better fit to the spectral
flattening just below $\approx 10$ TeV of the HAWC data. The gamma-ray
spectrum generated by Pop B dominates the low-energy range and account for the
spectral climb of the HAWC and HESS data below $\approx 5$ TeV.

The best-fit parameters are $\eta_{A}=0.0725$, $\eta_{B}=0.325$, $E_{c}=143$
TeV, $E_{b}=3.93$ TeV, and $\delta=1.25$. The total conversion efficiency from
the pulsar spin-down energy to the injected electron energy is $\approx 40\%$,
which is in a reasonable range. The observations of the x-ray PWN indicate
that the Geminga PWN can accelerate electrons up to $\sim 100$ TeV. The best-
fit $E_{c}$ is larger than 100 TeV, consistent with the x-ray observations.
The best-fit breaking energy of Pop B corresponds to a synchrotron cooling
time of $\approx 5000$ yr, indicating that Pop B may be trapped inside the PWN
for thousands of years before escaping. The best-fit $\delta$ is similar to
that obtained with the single power-law injection model in Sec. II.

As the electron injection process of Pop A is evident, we may also give a
rough estimate of the injection rate of Pop A based on the x-ray observations.
The brighter part of the southern tail measured by Ref. Pavlov  _et al._
(2010) is used for the estimate. The unabsorbed luminosity in the $0.3-8$ keV
band of this region is $\mathcal{L}=1.96\times 10^{29}$ erg s-1. The outflow
velocity is assumed as $v\sim c$. As the outflow is believed to be bent by the
ram pressure, only the length in the initial jet direction of this region is
useful for estimating the injection rate, which is $l\approx
75"\times\pi/$$\times 250~{}{\rm pc}=0.09$ pc. We take the approximation
provided by Ref. Longair (1994) for the calculation of synchrotron emission,
which can be expressed by

$\mathcal{L}=\int_{E_{e,1}}^{E_{e,2}}b(E_{e})\mathcal{Q}(E_{e})dE_{e}\,,$ (7)

where $E_{e,1}$ and $E_{e,2}$ are the electron energies corresponding to
synchrotron peak frequencies of $0.3$ and $8$ keV, respectively, and $b$ is
the synchrotron energy-loss rate. The electron energy spectrum $\mathcal{Q}$
has the same form as Eq. (5), while the normalization is determined by Eq.
(7). The total electron energy in this bright region is
$\mathcal{E}=\int\mathcal{Q}(E_{e})E_{e}dE_{e}$. Finally, the total injection
rate of the lateral outflows is estimated by $2\mathcal{E}v/l$, and the
conversion efficiency from the pulsar spin-down luminosity to Pop A is
$2\mathcal{E}v/l/L\approx 2\%$, where $L$ is the current spin-down luminosity
of Geminga. This efficiency is about three times smaller than the best-fit
$\eta_{A}$. However, the value of $\mathcal{L}$ is based on the assumption of
isotropic emission. Considering the large angle between the outflows and the
line of sight, the real $\mathcal{L}$ could be significantly larger, and the
conversion efficiency estimated by the x-ray observation could be in better
agreement with the fitting result.

Figure 3: Flux ratio of $F_{$$}/F_{$$}$ predicted by different models. The
blue, red, and green lines correspond to the one-population single power-law
injection model in the right of Fig. 1, the two-population injection model in
Fig. 2, and the two-population injection model with two-diffusion assumption
in Fig. 4, respectively.

Apart from the improvement of the goodness of fit, the two-population model
also predicts a different energy dependency of the gamma-ray flux ratio
$F_{$$}/F_{$$}$ from the single power-law model. As shown in Fig. 3, the
single power-law model in the right of Fig. 1 predicts a monotonically
decreasing flux ratio, while the flux ratio corresponds to the two-population
model in Fig. 2 decreases with energy in the low-energy range and increases
again above $\sim 10$ TeV. The difference is due to the different injection
spectra of these two models. In the case of continuous injection, the
injection spectrum may significantly affect the spatial distribution of
electrons. We give a qualitative explanation below.

Considering a constant injection case, the electron distribution can be
expressed by

$N(E_{e},r)=\int_{E_{e}}^{\infty}dE_{0}\,\frac{b(E_{0})}{b(E_{e})}\,\frac{1}{(\pi\lambda^{2})^{3/2}}\,{\rm
exp}\left(-\frac{r^{2}}{\lambda^{2}}\right)\,Q(E_{0})\,.$ (8)

It can be seen that $N(E_{e},r)$ is a superposition of electron distributions
with different initial energies $E_{0}$, and $Q(E_{0})$ is the weight of the
superposition. For $\delta>1$, electrons with larger energies have more
extended distributions, and the gamma-ray flux ratio $F_{$$}/F_{$$}$ tends to
decrease with energy. However, the two-population model assumes a cutoff term
for the injection spectrum. The number of electrons above the cutoff energy
decreases sharply, which means that the spatially extended component of the
electron distribution is seriously reduced. Thus, electron distributions with
higher energies are less extended, and $F_{$$}/F_{$$}$ increases with energy.

The energy dependency of the gamma-ray flux ratio can be a good criterion for
models. Due to the large uncertainties of the HESS spectrum above $\approx 5$
TeV, the energy dependency in high energies is not well constrained. The
LHAASO experiment Bai  _et al._ (2019) is expected to provide morphology
measurements for gamma-ray pulsar halos in a wide energy range in the coming
future, which may give clear judgments to different models.

##  IV Two-zone diffusion

In the above calculations, the electron diffusion is assumed to be
homogeneous. However, the slow-diffusion process around Geminga should not be
typical in the Galaxy. Considering the possible origins of the slow-diffusion
environment, the two-zone diffusion model may be a more reasonable assumption
(e.g., Ref. Fang  _et al._ (2018)). The slow-diffusion zone could be ascribed
to the streaming instability induced by the electrons escaping from Geminga
Evoli  _et al._ (2018); Mukhopadhyay and Linden (2021). It could also be
generated by an external source, such as the parent SNR of Geminga, which may
provide enough energy for the slow-diffusion environment Fang  _et al._
(2019). Either the interpretation suggests a slow-diffusion zone with a size
of $\sim 50$ pc. Furthermore, the multiwavelength gamma-ray spectrum of
another pulsar halo LHAASO J0622$+$3755 also indicates that the slow-diffusion
zone may not be larger than $\sim 50$ pc Fang  _et al._ (2021b). Below we
discuss the effects of two-zone diffusion, based on the two-population
injection model in Sec. III.

Firstly, the two-zone diffusion model can fit the observations with a $\delta$
smaller than the one-zone diffusion model required. When electrons escape from
the slow-diffusion zone, they rapidly spread in the ISM due to the much larger
diffusion coefficient outside. In this case, part of electrons is still
accumulated nearby the pulsar (e.g., within $None$ field around the pulsar),
while the electrons injected earlier may escape very far away from the pulsar
(e.g., out of $None$ field around the pulsar). This feature is clearly
illustrated by Fig. 1 of Ref. Fang  _et al._ (2018). Thus, the two-zone
diffusion model predicts a larger $F_{$$}/F_{$$}$ than the one-zone diffusion
case for the same $\delta$. For $\delta<1$, low-energy electrons have larger
propagation scales and can escape from the slow-diffusion zone more easily. It
means that the increase of $F_{$$}/F_{$$}$ is more significant for the low-
energy range of the gamma-ray spectrum, which may explain the large
$F_{$$}/F_{$$}$ indicated by the HAWC and HESS data $\lesssim 5$ TeV.

Figure 4: Left: spectra calculated under the two-zone diffusion assumption,
based on the two-population injection model in Fig. 2. Right: wide band gamma-
ray spectra extrapolated from the two-zone model in the left and the one-zone
model in Fig. 2, compared with the Fermi-LAT data obtained by Ref. Di Mauro
_et al._ (2019).

We assume $\delta=0.5$ that is predicted by Kraichnan’s theory and adopt the
two-zone diffusion model to explain the HAWC and HESS data. The diffusion
coefficient takes the form of

$D(E_{e},\textbf{\em r})=\left\\{\begin{aligned}
&D_{1}(E_{e})\,,\quad|\textbf{\em r}-\textbf{\em r}_{s}|<r_{\star}\\\
&D_{2}(E_{e})\,,\quad|\textbf{\em r}-\textbf{\em r}_{s}|\geq r_{\star}\\\
\end{aligned}\right.,$ (9)

where $D_{1}$ is the slow-diffusion coefficient used above, $D_{2}$ is the
typical diffusion coefficient of the Galaxy Yuan  _et al._ (2017), and
$r_{\star}$ is the size of the slow-diffusion zone. The gamma-ray profile
measured by HAWC indicates that $r_{\star}$ should not be smaller than
$\approx 25$ pc for high-energy electrons, while we find that a smaller
$r_{\star}$ is needed to interpret the large $F_{$$}/F_{$$}$ in the low-energy
range. Therefore, we assume $r_{\star}=30$ pc for $E_{e}\geq 30$ TeV, and
$r_{\star}=15$ pc for $E_{e}<30$ TeV. In fact, an energy-dependent $r_{\star}$
is also suggested by the self-excited scenario of slow diffusion Evoli  _et
al._ (2018). We show the comparison between the two-zone model and the
observations in the left panel of Fig. 4. The low-energy data are well
explained as expected. Most of the free parameters used here are similar to
the model in Fig. 2, while a larger conversion efficiency of 1.25 is required
for Pop B due to the smaller $r_{\star}$ in low energies. We also show the
energy dependency of the flux ratio of this model in Fig. 3, which has a
similar tendency with the model in Fig. 2.

Besides, the two-zone diffusion model may have a better agreement with the GeV
observation of the Geminga halo. The GeV Geminga halo is expected to be very
extended, and the spectral measurement is challenging due to the large
uncertainty of the gamma-ray background Xi  _et al._ (2019); Di Mauro  _et
al._ (2019). In the right panel of Fig. 4, we give a comparison between the
models in this work and the Fermi-LAT spectrum of the Geminga halo obtained by
Ref. Di Mauro  _et al._ (2019). The GeV spectrum is significantly affected by
the history of the diffusion pattern. For the self-generated scenario of slow
diffusion, the diffusion coefficient began to decrease to the current level
$\sim 10^{5}$ yr after the birth of Geminga Mukhopadhyay and Linden (2021).
For the alternative scenario proposed by Ref. Fang  _et al._ (2019), a time
scale of $\sim 10^{5}$ yr may also be needed for Geminga to approach the post-
shock region, where a slow-diffusion environment is expected. Thus, we assume
a $10^{5}$ yr delay for the emergence of the slow-diffusion zone. As shown in
Fig. 4, the GeV spectrum extrapolated from the one-zone diffusion model in
Fig. 2 is several times higher than the Fermi-LAT data, while the GeV fluxes
predicted by the two-zone diffusion model is significantly suppressed and
consistent with the Fermi-LAT spectrum. It should be emphasized that the above
comparison is qualitative. The calculation of the GeV spectrum is affected by
the low-energy electron injection spectrum, the time evolution and energy
dependency of $r_{\star}$, and the energy dependency of the diffusion
coefficient, all of which are currently not constrained.

##  V Conclusion

In this work, we first show that the latest HAWC and HESS gamma-ray spectra of
the Geminga halo can no longer be interpreted by the “simple model”, where a
single power-law injection spectrum and Kolmogorov’s energy index of the
diffusion coefficient ($\delta=1/3$) are assumed. Although a larger $\delta$
may account for the unexpected large flux ratio between the HESS and HAWC data
below $\approx 5$ TeV, the single power-law assumption for electron injection
is disfavored by the goodness-of-fit test. The reason is that the gamma-ray
spectrum derived from this simple assumption cannot reproduce the complex
spectral features of the measurements.

The Geminga PWN is the source of the electrons that light up the Geminga halo,
and the x-ray PWN and the TeV gamma-ray halo are generated by electrons with
almost the same energy range. Thus, it is meaningful to estimate the electron
injection spectrum by the x-ray observations. We propose a two-population
injection model based on the image and spectral measurements of the Geminga
x-ray PWN. One population (Pop A) is the freshly accelerated electrons that
escape from the PWN through rapid outflows, corresponding to the lateral tails
of the x-ray PWN. The spectrum of Pop A consists of a hard power-law term and
a high-energy cutoff term that describes the acceleration limit. The other
population (Pop B) is the electrons trapped longer in the PWN before escaping,
which corresponds to the axial tail of the x-ray PWN. A spectral break is
further assumed in the injection spectrum of Pop B due to the synchrotron
cooling inside the PWN.

The two-population injection model can consistently fit the HAWC and HESS
spectra with a reduced $\chi^{2}$ of 1.10, compared with $\chi^{2}=1.85$ for
the single power-law injection model. The high-energy HAWC spectrum is
reproduced by the cutoff term of the injection spectrum, and the HESS and low-
energy HAWC spectral features are well interpreted by the superposition of Pop
A and B. The required injection rate of Pop A is consistent with that roughly
derived from the x-ray observation, which further supports this model.
Intriguingly, a $\delta$ slightly larger than 1 is needed to fit the data. It
is larger than that predicted by Kolmogorov’s or Kraichnan’s theory and closer
to the value of Bohm diffusion. The two-population model also predicts a
different energy dependency of the gamma-ray profile from the single power-law
model, which will be tested by the energy-dependent morphological measurements
of LHAASO.

We also discuss the effect of two-zone diffusion on electron propagation,
which could be a more reasonable scenario considering the possible origins of
slow diffusion. Compared with the one-zone diffusion case, a smaller $\delta$
is needed to interpret the HAWC and HESS data. Besides, the GeV spectrum
predicted by the two-zone diffusion model may have better consistency with the
Fermi-LAT observation of the Geminga halo, although the uncertainties of the
data analysis and spectrum calculation in the GeV band are both large at
present.

###### Acknowledgements.

We thank P. Yin, X. Chen, and S. Xi for the helpful discussions. This work is
supported by the National Natural Science Foundation of China under Grants No.
12175248 and No. 12105292.

## References

  * Sudoh  _et al._ (2019) T. Sudoh, T. Linden, and J. F. Beacom, Phys. Rev. D 100, 043016 (2019), arXiv:1902.08203 [astro-ph.HE] .
  * Giacinti  _et al._ (2020) G. Giacinti, A. Mitchell, R. López-Coto, V. Joshi, R. Parsons, and J. Hinton, Astron. Astrophys. 636, A113 (2020), arXiv:1907.12121 [astro-ph.HE] .
  * Abeysekara  _et al._ (2017) A. Abeysekara _et al._ (HAWC), Science 358, 911 (2017), arXiv:1711.06223 [astro-ph.HE] .
  * Hooper  _et al._ (2017) D. Hooper, I. Cholis, T. Linden, and K. Fang, Phys. Rev. D 96, 103013 (2017), arXiv:1702.08436 [astro-ph.HE] .
  * Fang  _et al._ (2018) K. Fang, X.-J. Bi, P.-F. Yin, and Q. Yuan, Astrophys. J. 863, 30 (2018), arXiv:1803.02640 [astro-ph.HE] .
  * Profumo  _et al._ (2018) S. Profumo, J. Reynoso-Cordova, N. Kaaz, and M. Silverman, Phys. Rev. D 97, 123008 (2018), arXiv:1803.09731 [astro-ph.HE] .
  * Tang and Piran (2019) X. Tang and T. Piran, Mon. Not. Roy. Astron. Soc. 484, 3491 (2019), arXiv:1808.02445 [astro-ph.HE] .
  * Xi  _et al._ (2019) S.-Q. Xi, R.-Y. Liu, Z.-Q. Huang, K. Fang, and X.-Y. Wang, Astrophys. J. 878, 104 (2019), arXiv:1810.10928 [astro-ph.HE] .
  * Di Mauro  _et al._ (2019) M. Di Mauro, S. Manconi, and F. Donato, Phys. Rev. D 100, 123015 (2019), arXiv:1903.05647 [astro-ph.HE] .
  * Wang  _et al._ (2021) S.-H. Wang, K. Fang, X.-J. Bi, and P.-F. Yin, Phys. Rev. D 103, 063035 (2021), arXiv:2101.01438 [astro-ph.HE] .
  * Zhou (2021) H. Zhou (HAWC), PoS ICRC2019, 832 (2021).
  * Torres Escobedo  _et al._ (2021) R. Torres Escobedo _et al._ (HAWC), PoS ICRC2021, 842 (2021).
  * Mitchell  _et al._ (2021) A. M. W. Mitchell, S. Caroff, J. Hinton, and L. Mohrmann (H.E.S.S.), PoS ICRC2021, 780 (2021), arXiv:2108.02556 [astro-ph.HE] .
  * Minter and Spangler (1996) A. H. Minter and S. R. Spangler, Astrophys. J. 458, 194 (1996).
  * Fang  _et al._ (2021a) K. Fang, X.-J. Bi, S.-J. Lin, and Q. Yuan, Chin. Phys. Lett. 38, 039801 (2021a), arXiv:2007.15601 [astro-ph.HE] .
  * Blumenthal and Gould (1970) G. Blumenthal and R. Gould, Rev. Mod. Phys. 42, 237 (1970).
  * Manchester  _et al._ (2005) R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, Astron. J. 129, 1993 (2005), arXiv:astro-ph/0412641 .
  * Caraveo  _et al._ (2003) P. A. Caraveo, G. F. Bignami, A. De Luca, S. Mereghetti, A. Pellizzoni, R. Mignani, A. Tur, and W. Becker, Science 301, 1345 (2003).
  * Pavlov  _et al._ (2006) G. G. Pavlov, D. Sanwal, and V. E. Zavlin, Astrophys. J. 643, 1146 (2006), arXiv:astro-ph/0511364 [astro-ph] .
  * Pavlov  _et al._ (2010) G. G. Pavlov, S. Bhattacharyya, and V. E. Zavlin, Astrophys. J. 715, 66 (2010), arXiv:1002.1109 [astro-ph.HE] .
  * Posselt  _et al._ (2017) B. Posselt, G. Pavlov, P. Slane, R. Romani, N. Bucciantini, A. Bykov, O. Kargaltsev, M. Weisskopf, and C. Y. Ng, Astrophys. J. 835, 66 (2017), arXiv:1611.03496 [astro-ph.HE] .
  * Hui  _et al._ (2017) C. Y. Hui, J. Lee, A. K. H. Kong, P. H. T. Tam, J. Takata, K. S. Cheng, and D. Ryu, Astrophys. J. 846, 116 (2017), arXiv:1708.04179 [astro-ph.HE] .
  * Zirakashvili and Aharonian (2007) V. N. Zirakashvili and F. Aharonian, Astron. Astrophys. 465, 695 (2007), arXiv:astro-ph/0612717 .
  * Longair (1994) M. S. Longair,  _High-energy astrophysics. Vol. 2: Stars, the galaxy and the interstellar medium_ (Cambridge University Press, Cambridge, 1994).
  * Bai  _et al._ (2019) X. Bai _et al._ , (2019), arXiv:1905.02773 [astro-ph.HE] .
  * Evoli  _et al._ (2018) C. Evoli, T. Linden, and G. Morlino, Phys. Rev. D 98, 063017 (2018), arXiv:1807.09263 [astro-ph.HE] .
  * Mukhopadhyay and Linden (2021) P. Mukhopadhyay and T. Linden, (2021), arXiv:2111.01143 [astro-ph.HE] .
  * Fang  _et al._ (2019) K. Fang, X.-J. Bi, and P.-F. Yin, Astrophys. J. 884, 124 (2019), arXiv:1906.08542 [astro-ph.HE] .
  * Fang  _et al._ (2021b) K. Fang, S.-Q. Xi, and X.-J. Bi, Phys. Rev. D 104, 103024 (2021b), arXiv:2107.02140 [astro-ph.HE] .
  * Yuan  _et al._ (2017) Q. Yuan, S.-J. Lin, K. Fang, and X.-J. Bi, Phys. Rev. D 95, 083007 (2017), arXiv:1701.06149 [astro-ph.HE] .