# Aggregation of self-propelled particles with sensitivity to local order

Kunal Bhattacharya  Department of Industrial Engineering and Management, Aalto
University School of Science, 00076, Finland  Department of Computer Science,
Aalto University School of Science, 00076, Finland  Abhijit Chakraborty
<EMAIL_ADDRESS>Complexity Science Hub Vienna,
Josefstaedter Strasse 39, 1080 Vienna, Austria  Graduate School of Advanced
Integrated Studies in Human Survivability, Kyoto University, 1 Nakaadachi-cho,
Yoshida, Sakyo-ku, Kyoto 606-8306, Japan

###### Abstract

We study a system of self-propelled particles (SPPs) in which individual
particles are allowed to switch between a fast aligning and a slow nonaligning
state depending upon the degree of the alignment in the neighborhood. The
switching is modeled using a threshold for the local order parameter. This
additional attribute gives rise to a mixed phase, in contrast to the ordered
phases found in clean SPP systems. As the threshold is increased from zero, we
find the sudden appearance of clusters of nonaligners. Clusters of nonaligners
coexist with moving clusters of aligners with continual coalescence and
fragmentation. The behavior of the system with respect to the clustering of
nonaligners appears to be very different for values of low and high global
densities. In the low density regime, for an optimal value of the threshold,
the largest cluster of nonaligners grows in size up to a maximum that varies
logarithmically with the total number of particles. However, on further
increasing the threshold the size decreases. In contrast, for the high density
regime, an initial abrupt rise is followed by the appearance of a giant
cluster of nonaligners. The latter growth can be characterized as a continuous
percolation transition. In addition, we find that the speed differences
between aligners and nonaligners is necessary for the segregation of aligners
and nonaligners.

##  I Introduction

Collective motion observed in diverse natural and artificial systems has been
the subject of numerous experimental and theoretical investigations. Systems
that have been studied include fish schools Huth and Wissel (1994), birds
flocks Cavagna  _et al._ (2010), bacterial colonies Czirók  _et al._ (1996),
human crowds Silverberg  _et al._ (2013), as well as, synthetic microswimmer
assemblies Elgeti  _et al._ (2015) and robotic swarms Rubenstein  _et al._
(2014). Local interactions in such systems are understood to lead to the
emergence of global order or flocking states. This has been demonstrated in
self-propelled particle (SPP) models where particles are attributed the
tendency to align their direction of motion with their immediate spatial
neighbours in the presence of noise Vicsek  _et al._ (1995). Recent studies
have also focussed on the possible effects of environmental and individual-
level inhomogeneity on the flocking dynamics Chepizhko  _et al._ (2013);
Copenhagen  _et al._ (2016); Yllanes  _et al._ (2017). For example, disorder
is introduced in SPP models in the form of spatially distributed obstacles
Chepizhko  _et al._ (2013); Chepizhko and Peruani (2013), or, a finite
fraction of the particles are made non-aligners Copenhagen  _et al._ (2016);
Yllanes  _et al._ (2017). The dynamics in these systems shows the development
of phases with complex features like quasi-long-range order Chepizhko  _et
al._ (2013) and self-sorting Copenhagen  _et al._ (2016). In natural flocks,
the latter type of inhomogeneity could result from differences in signalling
and receptive behaviour or conflict in intentions. In an otherwise homogeneous
flock, behavioural shifts at individual levels could imply a certain fraction
of flock members spontaneously modifying their nature of motion.

Flocks of living organisms are known to arrange themselves into cohesive and
sometimes segregated units while performing activities like foraging and
migration Conradt and Roper (2005); Petit and Bon (2010); Bhattacharya and
Vicsek (2014). This is achieved through consensus decision making by flock
members while performing the activities that, in turn is a consequence of
mechanisms at the level of individuals. Such behavioural transitions between
different states have been documented in various species Buhl  _et al._
(2006); Daruka (2009); Miller and Gerlai (2012); Ginelli  _et al._ (2015).
Suitable modifications to simple SPP models have proven to be useful in
reproducing the spatio-temporal features of flocks with behavioural shifts.
Models have considered additional attributes to SPPs, like adaptive speed Buhl
_et al._ (2006); Li and Wang (2007); Farrell  _et al._ (2012), random fields
Bhattacharya and Vicsek (2010), and transition rates Ginelli  _et al._ (2015).

Experiments on fish schooling Katz  _et al._ (2011); Tunstrøm  _et al._ (2013)
and bacterial suspensions Cisneros  _et al._ (2011) have shown that the
individual speeds can vary depending on the local order parameter
(polarization). A model motivated by these experiments considered SPPs with
alignment interactions and a power-law dependence of the speed on the local
polarization Mishra  _et al._ (2012); Singh and Mishra (2020). This showed the
nucleation of static clusters and an inverse correlation between the speed and
the local density. Notably, for SPP systems without alignment such a speed-
density relationship leads to a motility-induced phase separation (MIPS) Cates
and Tailleur (2015) whereby two phases with distinct densities coexist in the
system. This is known to arise from the feedback between the slowing down and
crowding of the particles. In general, the variability in the speed both in
the absence Tailleur and Cates (2008); Fily and Marchetti (2012); Martin  _et
al._ (2021) and in the presence Farrell  _et al._ (2012); McCandlish  _et al._
(2012); Barré  _et al._ (2015); Sese-Sansa  _et al._ (2018); Van Der Linden
_et al._ (2019) of alignment has been shown to result in novel complex
phenomena in models and experiments.

In this paper we study an SPP model where particles can switch between a fast
aligning state and a slow non-aligning state depending on the local
orientational order parameter. An aligner becomes a non-aligner once the local
polarization falls below a threshold $\phi_{th}$, and conversely, a non-
aligner becomes an aligner if the local polarization rises above $\phi_{th}$.
Using the model we illustrate a mechanism where processing of local
information allows an SPP system to simultaneously organize into polarized
moving clusters as well as aggregations. The collisions between clusters play
a crucial role in such phase separation as we explain later. With our
numerical analyses we primarily focus on characterizing the clustering
behaviour.

In the absence of a threshold, or equivalently with $\phi_{th}\to 0$, the
model expectedly shows a order to disorder transition with the increase in
noise Vicsek  _et al._ (1995); Grégoire and Chaté (2004). We find that the
introduction of a finite threshold has a complex interplay with this
transition. In the steady state, the model with $\phi_{th}>0$ permits clusters
of aligners to coexist with those of non-aligners. The dynamics are found to
crucially depend on the level of the noise, the value of $\phi_{th}$ and the
overall density. At low noise, the aggregation behaviour of non-aligners can
be broadly categorized into two different regimes. For low enough densities,
an optimal value of $\phi_{th}$ is found to limit the growth of the largest
cluster of non-aligners; and at higher densities the latter is able to grow
macroscopically large when $\phi_{th}$ is increased.

Recent studies have considered SPP models relevant to the understanding of
epidemic and information spreading in populations of motile agents Paoluzzi
_et al._ (2018, 2020). The particles could either irreversibly or reversibly
switch between motile and non-motile states, and collectively exhibited
fractal aggregation and MIPS. These models considered switching rules based on
logic gates involving the states of colliding particles. In contrast, the
particles in our model change their states depending upon the orientations of
their neighbours.

The outline of the paper is as follows. We explain the details of the model in
section II. In section III we discuss the results of our numerical
investigation, and in section IV we conclude with a summary and final
observations.

##  II The Model

Figure 1: Snapshots of the dynamics are shown. The top row corresponds to
different parameter sets. The first three (from left) snapshots illustrate the
effect of increasing the density ($\rho$), such that, in the third we find an
emerging giant cluster of non-aligners. The aligners are denoted in blue (dark
gray) and the non-aligners in orange (light gray). The direction of a vector
points to the instantaneous direction of motion. Particles belonging to the
largest cluster of non-aligners are marked with magenta (gray). In the fourth,
the effect of increasing the noise ($\eta$) is shown. In the bottom row, A1
and A2 illustrate the mechanism of the emergence of cluster of non-aligners
from the collision of two moving clusters of aligners. The snapshots B1 and B2
show the growth of a cluster of non-aligners when a moving cluster of aligners
gets impacted. For both the pairs (A1,A2) and (B1,B2), the separation in time
is $50$ steps which is equivalent to $2.5$ units measured in the time scale
$r_{0}/v_{1}$. The number of particles in all the snapshots is $N=1024$. For
the bottom row, $\rho=1.2$, $\eta=0.15$ and $\phi_{th}=0.6$. Movies
corresponding to the first three snapshots (from left) are provided in the
Supplemental Material sm .

We consider $N$ self-propelled particles to be moving on a two dimensional
square area of linear size $L$ under periodic boundary conditions. The global
density of the system is given by $\rho=N/L^{2}$. At discrete times $t$, the
state of the $i$th particle is given by its position $\mathbf{r}_{i}^{t}$,
angle of the direction of motion $\theta_{i}^{t}$, and $s_{i}^{t}$, denoting
aligner ($s_{i}^{t}=1$) or non-aligner ($s_{i}^{t}=0$). The variables are
updated in the following way. First, the set of neighbours (${\cal N}_{i}$) of
the $i$th particle is enumerated, which comprises of all the particles that
are within a distance of $r_{0}$ from $i$. Then the local order parameter
($\phi_{i}$), that is the average normalized velocity within the
neighbourhood, is calculated as

$\phi_{i}^{t}=\frac{1}{1+k_{i}}\lvert\mathbf{n}_{i}^{t}+\sum_{j\in\mathcal{N}_{i}}\mathbf{n}_{j}^{t}\rvert.$
(1)

Here, $\mathbf{n}_{i}^{t}=(\cos\theta_{i}^{t},\sin\theta_{i}^{t})$ is a unit
vector pointing in the direction of motion of $i$, and $k_{i}$ is the number
of neighbours of $i$. Whether the particle $i$ would have the tendency to
align its direction of motion with its neighbours is decided depending on
$\phi_{i}^{t}$:

$s_{i}^{t+1}=\begin{cases}0,&\text{if }\phi_{i}^{t}\leq\phi_{th}\\\
1,&\text{if }\phi_{i}^{t}>\phi_{th},\end{cases}$ (2)

where $\phi_{th}$ is a parameter in the model. Lastly, the angle of heading
and the position are updated according to:

$\displaystyle\theta_{i}^{t+1}$
$\displaystyle=\arg\bigg{[}\mathbf{n}_{i}^{t}+s_{i}^{t+1}\sum_{j\in\mathcal{N}_{i}}\mathbf{n}_{j}^{t}+\alpha\sum_{j\in{\cal
N}_{i}}f_{ij}^{t}\hat{\mathbf{r}}_{ij}^{t}\bigg{]}+\eta\xi_{i}^{t},$ (3)
$\displaystyle\mathbf{r}_{i}^{t+1}$
$\displaystyle=\mathbf{r}_{i}^{t}+v(s_{i}^{t+1})\mathbf{n}_{i}^{t+1},$ (4)

where $v(s_{i}^{t+1})$ is the magnitude of the velocity depending on whether
the particle $i$ is an aligner or a non-aligner, $f_{ij}^{t}$ is the
interaction force between $i$ and its neighbour $j$,
$\hat{\mathbf{r}}_{ij}^{t}$ is a unit vector from $i$ towards $j$, $\alpha$ is
the strength of the interaction, $\xi_{i}^{t}\in[-\pi,\pi]$ is a delta-
correlated angular noise, and $\eta$ is the amplitude of the noise. For an
aligner, the dynamics represented in Eqs. 3 and 4 is similar to the Vicsek
model Vicsek  _et al._ (1995) with an additional short-range interaction
between particles Couzin  _et al._ (2002); Grégoire  _et al._ (2003); Grégoire
and Chaté (2004). For a non-aligner the second term inside the brackets in Eq.
3 is rather absent.

The force $f_{ij}$ depends on the distance of separation
$r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$, and comprises of finite repulsive
and attractive terms:

$f_{ij}=\begin{cases}-f_{r}\frac{r_{e}-r_{ij}}{r_{e}},&\text{if
}r_{ij}<r_{e}\\\ \frac{r_{ij}-r_{e}}{r_{a}-r_{e}},&\text{if
}r_{e}<r_{ij}<r_{a}\\\ \frac{r_{0}-r_{ij}}{r_{0}-r_{a}},&\text{if
}r_{a}<r_{ij}<r_{0},\end{cases}$ (5)

where $r_{e}$ is the equilibrium distance, $r_{a}$ is the distance at which
attraction is maximum, and, $f_{r}$ is the relative magnitude of the repulsive
force when $r_{ij}=0$. The two-body interaction helps to maintain a finite
packing density of particles similar to some of the systems McCandlish  _et
al._ (2012); Van Der Linden  _et al._ (2019) where MIPS-like phenomena is
evidenced.

In the simulations, we fix the following values for the parameters: $r_{0}=1$,
$r_{a}=0.625$, $r_{e}=0.25$, $f_{r}=1000$, and $\alpha=0.1$. For the non-
aligners, we will consider only the repulsive interaction to be present. Most
of the results, unless otherwise specified, are obtained with aligner speed,
$v_{1}=0.05$, non-aligner speed, $v_{0}=0.005$, and noise, $\eta=0.20$. Note,
in Eq. 1, $k_{i}=0$ implies $\phi_{i}\equiv 1$, and hence a non-aligner can
switch to the aligning state if it becomes isolated. To ensure that switching
arises only as a consequence of interactions in the neighbourhood, we prevent
isolated non-aligners from switching states. However, for the parameter ranges
that we investigate, this additional rule does not influence the macroscopic
behaviour of the system.

Figure 2: Sizes of the largest clusters of non-aligners ($M_{0}$) and aligners
($M_{1}$) are characterized at density $\rho=0.5$. (a) The dependence of
$M_{0}$ on noise amplitude ($\eta$) and threshold ($\phi_{th}$) is shown as a
heat map for $N=2^{10}$. The dashed line represents the equation
$\phi_{th}=\phi^{*}(\eta)$ where $\phi^{*}(\eta)$ is given by Eq. 6. (b) The
rate of switching per unit time per particle is plotted against $\phi_{th}$
for $N=2^{10}$ at $\eta=0.2$. Here, the aligners are switching to non-
aligners, and vice-versa. (Inset b) The variation of $M_{0}$ with $\phi_{th}$
for three noise amplitudes. The maximum of $M_{0}$ for $\eta=0.2$ occurs at
$\phi_{th}\sim 0.6$ which corresponds to the crossing of the switching rates.
(c) The dependence of $M_{0}$ on $N$ is shown, where $M_{0}$ is measured at
$\phi_{th}=0.6$ and $\eta=0.2$. In addition to the non-aligner speed
$v_{0}=0.005$, the dependence is also shown for $v_{0}=0.0$. The dashed lines
represent ordinary least squares fits having the form, $M_{0}=c_{0}+c_{1}\log
N$. For $v_{0}=0.005$, $c_{0}=-225(24)$ and $c_{1}=60(3)$; for $v_{0}=0.0$,
$c_{0}=-333(22)$ and $c_{1}=70(3)$. (d) The dependence of size of the largest
cluster of aligners ($M_{1}$) on $\phi_{th}$ at $\eta=0.2$. The different
symbols correspond to different values of $N$, as indicated in the legend.
Data collapse is obtained by scaling $M_{1}$ by the corresponding
$N^{\zeta_{1}}$ with $\zeta_{1}=0.78$.

##  III Results

At low noise and in the absence of a threshold ($\phi_{th}=0$) the system is
in a globally ordered state with a single macroscopically large cluster of
aligners ($s_{i}=1$). With the introduction of the switching behaviour
($\phi_{th}>0$), state of the particles become sensitive to fluctuations
occurring locally, and as a result the non-aligners ($s_{i}=0$) start
appearing in the system that are eventually separated from aligners due to the
difference in speeds. In the steady state, we find the system to be phase
separated into moving clusters of aligners and diffusing clusters of non-
aligners. If we observe the system in the very dilute limit, $\rho<0.1$ and
with $\phi_{th}=0$, we observe a phase with very small-sized clusters of
aligners due to the short-range two-body force Grégoire  _et al._ (2003). This
is different from the gaseous phase predicted for the original Vicsek model
Solon  _et al._ (2015); Chaté (2020). This also implies that for the higher
densities and for finite $\phi_{th}$ the large clusters of aligners have a
higher chance to coexist alongside clusters of non-aligners.

Our definition of a cluster is based on connecting neighbouring particles that
are in similar states ($r_{ij}<r_{0}$ and $s_{i}=s_{j}$). In Fig. 1 (top row)
we show snapshots of steady state configurations resulting from different
parameter values. Large fluctuations in $\phi_{i}^{t}$ primarily occur as a
result of collisions between clusters. Clusters of non-aligners form and grow
when moving clusters of aligners collide between themselves or with clusters
of non-aligners. This process is illustrated in the bottom row of Fig. 1.
Similarly, when a moving cluster of aligners grazes a cluster of non-aligners,
particles at the boundary of the latter switch their states in a short time
span to become aligners. Switching of particles at the boundary of a cluster
of non-aligners also happen as random events. To describe the generic
properties of the system we measure the sizes of the largest clusters of
aligners and non-aligners as functions of overall density and speed of the
particles.

Figure 3: (a) The dependence of size of the largest cluster of non-aligners
($M_{0}$) on the threshold ($\phi_{th}$) at noise, $\eta=0.2$, and density,
$\rho=1.6$. The different symbols correspond to different system sizes ($N$),
as indicated in the legend. The curves reveal a system size dependent
crossover occurring at $N\sim 2^{10}$. For larger systems $M_{0}$ has an
initial steep rise and then a gradual increase. (b) A susceptibility function
$\chi$ corresponding to $M_{0}$ is shown, where
$\chi=N\sqrt{\langle(M_{0}/N)^{2}\rangle-\langle(M_{0}/N)\rangle^{2}}$. From
around the crossover system size a second peak in $\chi$ starts becoming the
dominant maximum. The first peak corresponds to $\phi_{th}=\phi^{*}$ and
second peak corresponds to the second increase of $M_{0}$ at
$\phi_{th}=\phi_{th,c}$. (c) Curves from (a) with system sizes from $N\sim
2^{10}$ and above, are rescaled. By plotting $(M_{0}/N)N^{\beta/\nu}$ versus
$(\phi_{th}-\phi_{th,c})N^{1/\nu}$ the validity of Eq. 7 is illustrated. Here,
$\phi_{th,c}=0.74$. (Inset) Shows $(M_{0}/N)N^{0.5}$ versus $\phi_{th}$. The
collapse shows that the initial rise of $M_{0}$ at $\phi_{th}=\phi^{*}$ occurs
according to $M_{0}\sim N^{0.5}$. (d) Collapse of the second peak of
susceptibility (from (b)) is obtained following a procedure similar to that in
(c). This demonstrates the scaling ansatz in Eq. 8.

###  III.1 Low density regime – dependence on threshold and noise

In Fig 2 we show the behaviour of system at a density $\rho=0.5$. The
aggregation of non-aligners is only possible when the corresponding clean SPP
system is in the ordered phase. A small but finite $\eta$ ensures the presence
of clusters of aligners moving in different directions which can collide and
allow clusters of non-aligners to nucleate. The latter can not happen when
$\eta$ is large and the system is in a gas-like phase where large clusters of
aligners are absent. This is evidenced in Fig. 2(a) where we show the
dependence of the size of the largest cluster of non-aligners, $M_{0}$ on
$\eta$ and $\phi_{th}$. The plot also shows that $M_{0}$ attains its maximum
around $\phi_{th}=0.6$ and $\eta=0.2$.

For the individual particles in a cluster of aligners $\phi_{i}^{t}$ is high
in the ordered phase. However, during the collisions $\phi_{i}^{t}$ for
particles at the border of the colliding peripheries decreases momentarily. If
the drop in $\phi_{i}^{t}$ is less than $\phi_{th}$, then aligners switch to
become non-aligners. Therefore, an increase of $\phi_{th}$ leads to an
increase of switching events. For aligners we measure the rate of switching as
the number of switches to non-aligning states per unit time per aligner.
Similarly, we measure the switching rate for non-aligners. As Fig. 2(b) shows
this rate for aligners increases with $\phi_{th}$ and decreases for non-
aligners. Initially, at low $\phi_{th}$ the rate is much higher for non-
aligners implying non-aligners do not persist and proliferate, but when formed
almost instantly switch back to become aligners. But at higher values of
$\phi_{th}$ the switching rate for aligners overtakes that for non-aligners.
This aspect is also reflected when $M_{0}$ versus $\phi_{th}$ is examined
detail. As shown in Fig. 2 (Inset b) for different $\eta$’s, $M_{0}$ has a
sharp rise at a $\phi_{th}$ and then reaches a maximum. For $\eta=0.2$ the
maximum occurs at around $\phi_{th}=0.6$ and coincides with the point where
the rates cross each other. The sharp rise in $M_{0}$ is also found to be
noise dependent. Note that the switching rate for the non-aligners becomes
relatively a constant when $\phi_{th}$ is large. This is because inside the
bulk of a cluster the average separation is $r_{e}$ which also implies the
number of neighbours for a particle is $k\sim r_{0}^{2}/(r_{e}/2)^{2}$. With
$k$ randomly oriented neighbours inside a cluster of non-aligners the local
order parameter takes the typical value of $\phi\sim 1/\sqrt{k}\sim
r_{e}/(2r_{0})$ in the steady state Ginelli (2016).

An approximation for the value of the $\phi_{th}$ for which $M_{0}$ sharply
rises can be found in the following way. We consider a non-aligner on a
colliding boundary as illustrated in Fig. 1 (bottom row). We assume that on an
average half of the neighbours are non-aligners, and the rest half are
aligners. Therefore, the value of the local order parameter can be
approximated as

$\phi^{*}(\eta)=\frac{1}{2}\frac{\sin\eta\pi}{\eta\pi},$ (6)

where the expression to the right is half the polar order in an SPP system at
low noise and at sufficiently high densities Peruani  _et al._ (2010);
Dossetti  _et al._ (2009). The equation $\phi_{th}=\phi^{*}(\eta)$ is shown as
a dashed line in Fig. 2(a).

In Fig. 2(c) we find that the maximum values of $M_{0}$ increases as $\log N$.
In the same plot we show the case of $v_{0}=0$ where the non-aligners can only
rotate but not move. The dependence on $N$ is qualitatively similar in the
both the cases. While for $v_{0}=0$ the ejection of aligners from a cluster of
non-aligners occurs due to random switching events, with $v_{0}>0$ there is an
additional diffusion of the non-aligner particles before the switchings
happen. In Sec. III.4 we show the dependence of $M_{0}$ on $v_{0}$.

The abrupt increase in $M_{0}$ as $\phi_{th}$ increases in the low noise
regime coincides with a decrease in the size of the largest cluster of
aligners ($M_{1}$). This is visible in Fig 2(d) where $M_{1}$ is plotted as a
function of $\phi_{th}$ for different $N$ at $\eta=0.2$. By tuning the
exponent in the relation $M_{1}\sim N^{\zeta_{1}}$ we obtain the best collapse
for different $N$ with $\zeta_{1}=0.78$.

Figure 4:  (a) Normalized distributions of cluster sizes of non-aligners for
different system sizes, indicated in the legend. The distributions have been
obtained at $\eta=0.2$, $\rho=1.6$ and $\phi_{th}=0.74=\phi_{th,c}$. (b) A
collapse of the distributions from (a), obtained by plotting
$n(m_{0})N^{\epsilon 2}$ versus $m_{0}/N^{\epsilon_{1}}$ and tuning
$\epsilon_{1}$ and $\epsilon_{2}$. In (b) $\epsilon_{1}=0.96(2)$ and
$\epsilon_{2}=1.98(1)$. The collapse indicates a power-law
$n(m_{0})\sim{m_{0}}^{-\tau}$ in the large $N$ limit with
$\tau=\epsilon_{2}/\epsilon_{1}=2.06(3).$

###  III.2 High density regime

In the high density and low noise regime, we observe the largest cluster of
non-aligners grows with $\phi_{th}$ whose size is of the order of $N$ when
$\phi_{th}$ is unity. An incipient cluster in this regime is shown in Fig. 1
(top row - third from left). In addition, we find the behaviour of the system
to be strongly dependent on $N$. In Fig. 3(a) we plot the fraction $M_{0}/N$
for different system sizes. From the different curves we observe that the
generic dependence of $M_{0}/N$ on $\phi_{th}$ is different in smaller and
larger values of $N$ with a crossover occurring at $N\sim 2^{10}$. For the
smaller systems the growth in the largest cluster primarily occurs when
$\phi_{th}$ is between 0.6 and 0.7. For the larger systems there is an initial
rapid increase in $M_{0}$ that is similar to that observed at low densities,
and is followed by a gradual further growth.

Considering the growth of the larger cluster of non-aligners to be a
percolation phenomenon that occurs with respect to the tuning of $\phi_{th}$,
we calculate a susceptibility corresponding to the order parameter $M_{0}/N$,
given by $\chi=N\sigma$, with
$\sigma^{2}=\langle(M_{0}/N)^{2}\rangle-\langle(M_{0}/N)\rangle^{2}$ Radicchi
and Fortunato (2009), where angular brackets denote averaging in the steady
state. The plot of $\chi$ as a function of $\phi_{th}$ in Fig. 3(b)
demonstrates a crossover in the finite-size effect. For $N<2^{10}$ there is
only a single maximum that shifts to the left of the $\phi_{th}$ axis on
increasing $N$. With $N\geq 2^{10}$ we find the emergence of a second peak in
$\chi$ which is the dominant maximum as $N$ increases further. Ideally, for a
given $N$, the position of the (second) maximum of $\chi$ is expected to
provide the critical thresholds (pseudocritical point) $\phi_{th,c}(N)$.
Observing that accurately locating the second maximum can be difficult for the
smaller system sizes we circumvent the problem in the following way. For the
finite-size effects in percolation we assume the relations,

$\displaystyle M_{0}/N$ $\displaystyle=$ $\displaystyle
N^{-\beta/\nu}\mathcal{F}[(\phi_{th}-\phi_{th,c})N^{1/\nu}],$ (7)
$\displaystyle\chi$ $\displaystyle=$ $\displaystyle
N^{\gamma/\nu}\mathcal{G}[(\phi_{th}-\phi_{th,c})N^{1/\nu}],$ (8)

where $\phi_{th,c}$ is the critical threshold in the infinite size limit
($N\to\infty$), and $\nu$, $\beta$, and $\gamma$ are the critical exponents
characterizing a second order percolation transition. Using the above two
relations and the fact that $\chi=N\sigma$ we get a hyperscaling relation,

$\gamma/\nu=1-\beta/\nu.$ (9)

At different values of $\phi_{th}$ we fit power laws to the data corresponding
to $M_{0}/N$ versus $N$. This gives us a set of trial values for the exponent
$\beta/\nu$. Similarly, we obtain a set of trial values of $\gamma/\nu$ from
$\chi$ versus $N$ at different $\phi_{th}$. Then we obtain the critical point
$\phi_{th,c}$ by locating the $\phi_{th}$ at which $\beta/\nu$ and
$\gamma/\nu$ satisfies Eq. 9. This method yields $\beta/\nu=0.035(9)$,
$\gamma/\nu=0.965(4)$, and $\phi_{th,c}=0.740(5)$. The error estimates in the
exponents correspond to power law fits at $\pm 0.005$ from $\phi_{th,c}$. (In
the Appendix we provide expressions for $M_{0}$ and $\phi_{th,c}$ from a
reaction-limited description.)

To determine $\nu$, we first scale the y-axis of Fig. 3(a) by multiplying with
$N^{\beta/\nu}$. Then upon fixing a value of $(M_{0}/N)N^{\beta/\nu}$ around
$1.0$ we obtain the corresponding values of $\phi_{th}$ for different $N$. We
estimate the value of $1/\nu$ from the slope of the line fitted with
$\log|\phi_{th}\left(N\right)-\phi_{th,c}|$ versus $\log N$. Repeating the
process for different values of $(M_{0}/N)N^{\beta/\nu}$, we get
$\nu=2.16(3)$. Using the above values for the scaling exponents and
$\phi_{th,c}$ we obtain data collapse for $M_{0}/N$ and $\chi$ shown in Fig.
3(c) and Fig. 3(d), respectively. The collapse of the curves for different $N$
when $\phi_{th}$ is close to $\phi_{th,c}$ shows that the scaling forms in Eq.
7 and 8 hold true.

Similar to the case of low density, there is an initial rapid increase in
$M_{0}$ at $\phi_{th}=\phi^{*}$. This is characterized by solely scaling the
$M_{0}/N$ axis, and collapsing the curves for different $N$ as shown in the
inset of Fig. 3(c). The scaling shows that as $\phi_{th}$ crosses $\phi^{*}$,
$M_{0}$ abruptly increases from $O(1)$ to $O(\sqrt{N})$. The latter increase
in $M_{0}$ also coincides with a fall in $M_{1}$ (not shown). For
$\phi_{th}<\phi^{*}$ we find $M_{1}\sim N^{\zeta_{2}}$ with $\zeta_{2}=0.80$
that is quite close to $\zeta_{1}$.

Additionally, we also studied the cluster size distribution of the non-
aligners. We obtained the statistics by observing the systems at
$\phi_{th}=\phi_{th,c}$. The normalized distributions $n(m_{0})$ for the
cluster sizes ($m_{0}$) of non-aligners are plotted in Fig. 4(a) for different
$N$. We assume power-law distributed cluster sizes and finite-size effects are
present, such that $n(m_{0})\sim m_{0}^{-\tau}f\left(m_{0}/M_{0}\right)$,
where $\tau$ is the Fisher exponent. The system being at $\phi_{th,c}$ we
expect $M_{0}\sim N^{\epsilon_{1}}$, with $\epsilon_{1}=1-\beta/\nu$, and
therefore, the above power-law can be recast into the following form:

$n(m_{0})=N^{-\epsilon_{2}}D\left(m_{0}/N^{\epsilon_{1}}\right),$ (10)

where, the scaling function $D(x)\sim x^{-\tau}$ for $x\to 0$, and $D(x)$
decreases faster than a power law for $x\gg 1$. This implies
$\tau=\epsilon_{2}/\epsilon_{1}$. In Fig. 4(b) we plotted
$n(m_{0})N^{\epsilon_{2}}$ versus $m_{0}/N^{\epsilon_{1}}$ and tuned the
values of $\epsilon_{1}$ and $\epsilon_{2}$ to get a collapse of the
distribution for different $N$. The latter allows us to validate Eq. 10. We
get the best collapse for $\epsilon_{1}=0.96(2)$ and $\epsilon_{2}=1.98(1)$,
which implies $\tau=2.06(3)$.

Figure 5: The following quantities of interest are plotted as functions of
density ($\rho$) at noise, $\eta=0.2$, and threshold, $\phi_{th}=0.6$:
$M_{0}/M_{0,\textrm{max}}$ – size of the largest cluster of non-aligners
normalized by its maximum value within the range of investigation,
$M_{1}/M_{1,\textrm{max}}$ – normalized size of the largest cluster of
aligners, and, $\Phi/\Phi_{\textrm{max}}$ – normalized polarization for the
system. The dashed vertical line indicates the density $\rho_{CP}=1.42$ that
corresponds to the critical filling factor for continuum percolation.  Figure
6: (a) Dependence of size of the largest cluster of aligners ($M_{1}$) on non-
aligners speed ($v_{0}$), for three different thresholds ($\phi_{th}$). The
aligner speed, $v_{1}=0.05$, noise, $\eta=0.2$, number of particles, $N=1024$,
and density $\rho=0.5$. (b) Dependence of size of the largest cluster of non-
aligners ($M_{0}$).

###  III.3 Generic dependence on density

After observing that switching of particle states in ordered flocks produces
two distinct type of mixed phases – non-percolating and percolating depending
on the density $\rho$, we investigate how the relevant quantities continuously
vary as functions of $\rho$. In Fig. 5 we show the variation of $M_{0}$ and
$M_{1}$ at low noise ($\eta=0.2$) and substantial presence of switching
activity ($\phi_{th}=0.6$). In addition to the sizes of largest clusters we
also show the polarization $\Phi$ Cavagna  _et al._ (2010); Copenhagen  _et
al._ (2016) which is a measure of the degree of global order in the system,
and is defined as
$\Phi=\langle\lvert(1/N)\sum_{i=1,N}\mathbf{n}_{i}^{t}\rvert\rangle$.

We compare $M_{0}$, $M_{1}$ and $\Phi$ scaled by their respective maximum
values with $\rho$ being in the range $0.01$ to $10.0$. It is known that in a
pure SPP system for a fixed $\eta$, long-range order vanishes when density is
lowered Czirók and Vicsek (2000). This phenomena partially underlies the
dynamics observed in our case where switching is allowed. The relatively low
degree of order as reflected in the values of $M_{1}/M_{1,max}$ and
$\Phi/\Phi_{max}$ at very low densities is result of the lack of order in the
pure system. The maximum in the values of $M_{1}$ and $\Phi$ occur at around
$\rho\gtrsim 0.1$. In this regime, the absence of switching implies that the
system has high degree of order and a macroscopically large cluster of
aligners always exists in the steady state. When switching is present we
observe the non-percolating mixed phase where $M_{0}\sim\log N$. A further
increase in the density would result in the enhancement of order and aligner
cluster size in the pure system. However, for the mixed phase the aligner
cluster size and order decreases with further increase in density. The cluster
size for non-aligners increases monotonically. The possibility of the
formation of a giant cluster of non-aligners may be considered similar to the
situation of continuum percolation (CP) of overlapping disks with radii
$R=r_{0}/2$. Noting that the filling fraction is defined as $\pi R^{2}\rho$
and that the estimate for the critical value is around $1.12$ Mertens and
Moore (2012), the corresponding critical density is $\rho_{\textrm{CP}}=1.42$.
This density, therefore, would signify the bordering between non-percolating
and percolating mixed phases. The latter is evident in the case of
$\rho=1.6>\rho_{\textrm{CP}}$.

###  III.4 Dependence on speed differences

Lastly, we study how the difference between the speeds of the aligners and
non-aligners governs the evolution of system. We fix the aligner speeds to
$v_{1}=0.05$ and vary the speed of the non-aligners, $v_{0}$. We plot $M_{1}$
and $M_{0}$ as functions of $v_{0}$ for three different values of $\phi_{th}$
in Fig. 6. It is apparent that the mixed phase ($\phi_{th}$=0.6) where
macroscopically large aligner and non-aligner clusters coexist is delicately
dependent on the value $v_{0}$. As $v_{0}$ approaches $v_{1}$, $M_{0}$ is
found to decrease and $M_{1}$ is found to increase. The diffusion of non-
aligners occurs with a diffusion constant that is proportional to $v_{0}^{2}$
Peruani  _et al._ (2010); Peruani and Baer (2013). Therefore, the rate of
ejection of aligners from the boundary of a cluster of non-aligners also
increases as the non-aligner speed is increased. In addition, as the relative
difference of speeds vanishes, the non-aligners formed after a collision
effectively fail to segregate and to eventually proliferate. As a result of
the above the effect of a finite $\phi_{th}$ diminishes and large clusters of
non-aligners are rarely observed. These results show that the formation of
mixed phases that is controlled by $\phi_{th}$ is also dependent on the
difference in speeds. We have shown the dependence on $v_{0}$ in the non-
percolating regime, but the indications are similar for the percolating case
as well.

##  IV Conclusions

We studied a system where self-propelled particles were allowed to switch
states between fast aligners and slow non-aligners based on the degree of
alignment in their neighbourhood. In the steady state, the system segregated
into separate clusters of aligners and non-aligners. In the mixed phase, the
largest cluster of aligners was found to vary algebraically with the system
size. However, depending on the density of the system, the aggregation of the
non-aligners appeared to be very different. For the low densities, the largest
cluster of non-aligners reached a maximum size for an optimal noise and an
optimal threshold. For high densities, after an initial abrupt increase, a
giant percolating cluster could emerge with the increase in the threshold.
Also, the behaviour for small system sizes appeared to be very different. The
boundary between the density regimes roughly coincided with the density
corresponding to the critical filling factor for a continuum percolation
transition. Irrespective of the density, the separation of speeds seemed to be
a necessary condition for the model to display the segregation of non-
aligners. When the speeds become comparable, large clusters of non-aligners
are predominantly absent. Although the appearance of the giant cluster
conforms to a set of finite-size scaling hypotheses, the transition could be
non-universal Liu and Pandey (1997) with dependence on the noise amplitude and
density, through the functions $\mathcal{F}$ and $\mathcal{G}$ in Eqs. 7 and
8.

The percolation of clusters was recently studied Kyriakopoulos  _et al._
(2019) in the classical Vicsek model. Unlike our model, the SPPs in the Vicsek
model always remain aligners (without switching) and attraction-repulsion
forces are absent. The authors investigated the global connectivity of
clusters with increase in the global density ($\rho$) along both the
longitudinal and the transverse directions with respect to the direction of
global order. They estimated a critical density,
$\rho_{c}=1.96~{}(>\rho_{CP})$. Similar to the current model, if we denote the
size of the largest cluster (of aligners) in the Vicsek model as $M_{1}$, then
near to $\rho_{c}$, the dependence on $N$ may be characterized by using
$M_{1}\sim N^{\zeta}$. Using the reported Kyriakopoulos  _et al._ (2019)
values of the different critical exponents, $\zeta$ is found be in the range
$0.95-1.00$. In our model, for the density regimes investigated, and when the
clusters of aligners are macroscopically large ($\phi_{th}<\phi^{*}$) we find
$\zeta$ to be in the range $0.78-0.80$. Taken together, we expect that our
model was investigated at densities which are still lower than the critical
density that would be needed for the percolation of clusters of aligners if
switching is absent. Also note, that in our model, clusters of non-aligners
are formed mainly due to the collisions between clusters and the speed
difference between aligners and non-aligners.

Recent advances in living active matter have found that modifications in
individual behaviour through the sensing of local densities lead to the
formation of regions of orientational disorder and aggregations Ling  _et al._
(2019); Leggett  _et al._ (2019); Klamser  _et al._ (2021); Bäuerle  _et al._
(2018). Similar observations are made in experiments with active colloidal
systems employing different methods to program the particle motion, like
optical feedback loops and field modulations Bäuerle  _et al._ (2018);
Lavergne  _et al._ (2019); Soma  _et al._ (2020). Therefore, the observed
macroscopic behaviour in our model could be of relevance, for example, to
active colloids with setups allowing the particles to sense and respond to the
average orientation of neighbours Bäuerle  _et al._ (2020), to the design and
control of robot swarms Werfel  _et al._ (2014); Rubenstein  _et al._ (2014),
and in general, to systems exhibiting both polar order and MIPS-related
behaviour Van Der Linden  _et al._ (2019); Crosato  _et al._ (2019); Sesé-
Sansa  _et al._ (2021). Also, owing to the additional state variable in our
model, the latter can be contrasted with the study of clustering and
percolation in the classical Vicsek model Kyriakopoulos  _et al._ (2019), and
similarities with models on information spreading in motile collectives can be
further explored Paoluzzi  _et al._ (2018, 2020).

Also, there appears to be scope for additional complexity in the current
model. We have assumed that switching behaviour is symmetric in terms of
having a single threshold for both aligners as well as non-aligners. In a more
general scheme, there can be two different thresholds for particles of either
types. Currently, the additive noise causes the moving clusters to collide
which generates the non-aligners. It would be interesting to test the model at
zero-noise Chakraborty and Bhattacharya (2016) but with other forms of
disorder like boundaries, obstacles Chepizhko  _et al._ (2013) and quenched
non-aligners Yllanes  _et al._ (2017). The phase separation between aligners
and non-aligners could be also studied in systems without self-propulsion, for
example, in Brownian walkers with a velocity alignment interaction Dossetti
and Sevilla (2015) and in the Vicsek model on the lattice Bhattacherjee  _et
al._ (2014).

##  V Appendix

###  V.1 Mean field description

Here we provide a reaction-limited description of the system in the high
density regime based on our simulations and neglecting the spatial
correlations. We consider the system to be consisting of the following types
of particles: particles that are part of the giant cluster of non-aligners
($M_{0}$), aligners ($N_{1}$), and non-aligners that are not part of the giant
cluster ($N_{0}^{\prime}$), such that

$M_{0}+N_{0}^{\prime}+N_{1}=N.$ (11)

As seen above with regards to Fig. 3, as $\phi_{th}$ crosses $\phi^{*}$,
$M_{0}$ becomes $O(\sqrt{N})$. On further increasing $\phi_{th}$, when the
latter reaches $\phi_{th,c}$, the incipient giant cluster is observed. We
model the growth of this cluster using the following equation:

$\frac{\mathrm{d}M_{0}}{\mathrm{d}t}=(\phi_{th}-\phi^{*})AM_{0}N_{1}-(1-\phi_{th})BM_{0}.$
(12)

The first term to the right accounts for the collision between the aligner
particles and the giant cluster of non-aligners by which particles are added
to the perimeter of the latter. Given that it is a percolating cluster which
is growing, and is far from being circular in shape the perimeter is assumed
to vary as $M_{0}$. A typical such cluster is shown in the Fig. 1 (top row -
third from the left). In general, for a percolating cluster the relation
between the perimeter (or hull) and the mass is given by
$H_{0}\sim{M_{0}}^{x}$, with $x=d_{h}/d_{f}$. Here $d_{f}$ and $d_{h}$ are the
fractal and hull dimensions, respectively, and can be computed using the
relations $d_{f}=d-\beta/\nu_{L}$ and $d_{h}=1+1/\nu_{L}$ Saleur and
Duplantier (1987); Stauffer and Aharony (2003). In the current context, the
exponents from our numerical calculations in Sec. III.2 would indicate
$x\approx 1$. The exponent $\nu_{L}=\nu/d$, where the dimensionality of space,
$d=2$. For classical percolation in $d=2$, $d_{f}=91/48$ and $d_{h}=7/4$,
which gives $x\simeq 0.92$ Stauffer and Aharony (2003). On the other, clusters
of aligners are mostly small in size in this regime, therefore, the dependence
is taken to be proportional to the total number of aligners, $N_{1}$. The
colliding aligners are expected to have their local order parameter
distributed around $\phi^{*}$. Keeping other factors unchanged, as $\phi_{th}$
approaches $\phi^{*}$ from below and eventually crosses $\phi^{*}$, more and
more aligners are expected to switch their states. We approximate this
dependence on $\phi_{th}$ as, $\phi_{th}-\phi^{*}$. The coefficient, $A$
contains factors in the collision rate including the speed of the group of
aligners Peruani  _et al._ (2010); Peruani and Baer (2013) and is independent
of $\phi_{th}$. The second term accounts for the loss of particles by which
particles near the perimeter of the giant cluster switch to aligning states
and detach from the latter. While within the bulk, the local order parameter
would typically be around $r_{e}/(2r_{0})$ (discussion on Fig. 2(b)), for the
non-aligners at the boundary we expect larger fluctuations and hence a flatter
distribution extending up to unity. We assume that the rate of switching is
proportional to $1-\phi_{th}$. The coefficient $B$ takes into account other
factors independent of $\phi_{th}$.

We observe in our simulations that apart from the ones in the giant cluster,
the non-aligners are formed in the process of collision between small clusters
of non-aligners. The resulting clusters of non-aligners are also small and are
not stable. We model this process using the following equation:

$\frac{\mathrm{d}N_{0}^{\prime}}{\mathrm{d}t}=(\phi_{th}-\phi^{*})kA{N_{1}}^{2}-(1-\phi_{th})BN_{0}^{\prime}.$
(13)

The first term accounts for formation of the small clusters of non-aligners.
We assume that the coefficient in the collision rate accounting for factors
independent of time and $\phi_{th}$ is only different by a multiplicative
constant from the coefficient in Eq. 12. The second term is similar to the
loss term in $\mathrm{d}{M_{0}}/\mathrm{d}t$.

Figure 7: Results from the simulation in the high density regime ($\rho=1.6$,
$\eta=0.2$, and $N=2^{13}$) is compared with the steady state solutions from
the model. The dashed line denotes the function
$f(\phi_{th})=(1-\phi_{th})/(\phi_{th}-\phi^{*})$. The circles show the
dependence on $\phi_{th}$ for the number of aligners ($N_{1}$) scaled by a
constant ($b$). Choosing $b/N\approx 1/12$ shows the validity of Eq. 14.
Similarly, the triangles and the squares show the ranges of validity for the
relation $N_{0}^{\prime}=kN_{1}$ and Eq. 15, respectively where $k=1$. The
plot of $N_{1}/b$ in the low density regime is shown using the crosses. Here,
we take $\phi^{*}=0.5$ which provides the best fit instead of $\phi^{*}=0.47$
which is obtained from Eq. 6.

We are interested in the steady state dependence of $M_{0}$ on $\phi_{th}$.
Therefore, we set the right hand side of Eq. 12 to zero. Assuming a non-zero
finite solution for $M_{0}$ we get the following steady solution for $N_{1}$:

$N_{1}^{s}=bf(\phi_{th}),$ (14)

where $f(\phi_{th})=(1-\phi_{th})/(\phi_{th}-\phi^{*})$ and $b=B/A$.
Similarly, by equating the right hand side of Eq. 13 to zero, we get
$N_{0}^{\prime s}=kN_{1}^{s}$. Next, the steady state solution for $M_{0}$ is
found by using $N_{1}^{s}$ and $N_{0}^{\prime s}$ in Eq. 11:

$M_{0}^{s}=N-b(1+k)f(\phi_{th}).$ (15)

In Fig. 7 we compare the steady state solutions for $N_{1}$, $N^{\prime}_{0}$
and $M_{0}$ with the numerical results. As mentioned above, the derived
expression appears to be valid for $\phi_{th}>\phi_{th,c}$ when the incipient
giant cluster is already present. Below $\phi_{th,c}$, clusters of non-
aligners keep continually forming and fragmenting, clusters of aligners are
relatively larger in size, and spatial correlations cannot be neglected in the
description of the dynamics.

Considering that the collision coefficient $A$ is inversely proportional to
the area ($L^{2}$) Peruani  _et al._ (2010); Peruani and Baer (2013) we can
also rewrite Eq. 15 as

$M_{0}^{s}=N\left[1-(c/\rho)f(\phi_{th})\right],$ (16)

where $c/\rho=b(1+k)/N$. The above expressions for $M_{0}^{s}$ imply that when
$N$ increases at a fixed density or when $\rho$ increases, the transition to a
global connectivity becomes faster. This expression, however, is not valid in
the limit $\rho\to 0$. For the latter limit in Eq. 12, the second term would
be the dominant term for all values of $\phi_{th}$ and $M_{0}$ would be zero
in the steady state. Alternately, the dynamics at low density cannot be
described by Eqs. 12 and 13, and as a result, Eq. 14 does not hold. The
deviation of $N_{1}^{s}/b$ at low density from the function $f(\phi_{th})$ is
shown in Fig. 7.

We obtain an estimate of $\phi_{th,c}$ by assuming that once the critical
threshold is exceeded, the number of non-aligners inside the largest cluster
becomes greater than the number of non-aligners present outside:
$M_{0}^{s}>N_{0}^{\prime s}$. Here using the expressions $M_{0}^{s}$ and
$N_{0}^{\prime s}$, we get a lower bound on $\phi_{th}$,

$\phi_{th,c}=\phi^{*}+\frac{3c}{3c+2\rho}(1-\phi^{*}),$ (17)

where we have set $k=1$ as shown in Fig. 7. Therefore, in this reaction
dominated description, the limits $N\to\infty$ at a fixed density or
$\rho\to\infty$ imply that $\phi_{th,c}\to\phi^{*}$ from above. Finally, given
that the function $f(\phi_{th})$ is analytic in the range
$\phi^{*}<\phi_{th}\leq 1$, we can express as a first order approximation:
$M_{0}^{s}(\phi_{th})-M_{0}^{s}(\phi_{th,c})\sim(\phi_{th}-\phi_{th,c})$ for
$\phi_{th}\gtrsim\phi_{th,c}$, giving $\beta{\textrm{(mean-field)}}=1.$

###### Acknowledgements.

KB and AC acknowledge helpful discussions with Shakti N. Menon, Sasidevan
Vijayakumar, and Sumanta Kundu.

## References

  * Huth and Wissel (1994) A. Huth and C. Wissel, Ecological modelling 75, 135 (1994).
  * Cavagna  _et al._ (2010) A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale, Proceedings of the National Academy of Sciences 107, 11865 (2010).
  * Czirók  _et al._ (1996) A. Czirók, E. Ben-Jacob, I. Cohen, and T. Vicsek, Physical Review E 54, 1791 (1996).
  * Silverberg  _et al._ (2013) J. L. Silverberg, M. Bierbaum, J. P. Sethna, and I. Cohen, Physical Review Letters 110, 228701 (2013).
  * Elgeti  _et al._ (2015) J. Elgeti, R. G. Winkler, and G. Gompper, Reports on progress in physics 78, 056601 (2015).
  * Rubenstein  _et al._ (2014) M. Rubenstein, A. Cornejo, and R. Nagpal, Science 345, 795 (2014).
  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Physical review letters 75, 1226 (1995).
  * Chepizhko  _et al._ (2013) O. Chepizhko, E. G. Altmann, and F. Peruani, Physical review letters 110, 238101 (2013).
  * Copenhagen  _et al._ (2016) K. Copenhagen, D. A. Quint, and A. Gopinathan, Scientific reports 6, 31808 (2016).
  * Yllanes  _et al._ (2017) D. Yllanes, M. Leoni, and M. Marchetti, New Journal of Physics 19, 103026 (2017).
  * Chepizhko and Peruani (2013) O. Chepizhko and F. Peruani, Physical review letters 111, 160604 (2013).
  * Conradt and Roper (2005) L. Conradt and T. J. Roper, Trends in ecology & evolution 20, 449 (2005).
  * Petit and Bon (2010) O. Petit and R. Bon, Behavioural Processes 84, 635 (2010).
  * Bhattacharya and Vicsek (2014) K. Bhattacharya and T. Vicsek, Journal of the Royal Society Interface 11, 20140674 (2014).
  * Buhl  _et al._ (2006) J. Buhl, D. J. Sumpter, I. D. Couzin, J. J. Hale, E. Despland, E. R. Miller, and S. J. Simpson, Science 312, 1402 (2006).
  * Daruka (2009) I. Daruka, Proceedings of the Royal Society of London B: Biological Sciences 276, 911 (2009).
  * Miller and Gerlai (2012) N. Miller and R. Gerlai, PLoS One 7, e48865 (2012).
  * Ginelli  _et al._ (2015) F. Ginelli, F. Peruani, M.-H. Pillot, H. Chaté, G. Theraulaz, and R. Bon, Proceedings of the National Academy of Sciences 112, 12729 (2015).
  * Li and Wang (2007) W. Li and X. Wang, Physical Review E 75, 021917 (2007).
  * Farrell  _et al._ (2012) F. D. C. Farrell, M. C. Marchetti, D. Marenduzzo, and J. Tailleur, Phys. Rev. Lett. 108, 248101 (2012).
  * Bhattacharya and Vicsek (2010) K. Bhattacharya and T. Vicsek, New Journal of Physics 12, 093019 (2010).
  * Katz  _et al._ (2011) Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin, Proceedings of the National Academy of Sciences 108, 18720 (2011).
  * Tunstrøm  _et al._ (2013) K. Tunstrøm, Y. Katz, C. C. Ioannou, C. Huepe, M. J. Lutz, and I. D. Couzin, PLoS Comput Biol 9, e1002915 (2013).
  * Cisneros  _et al._ (2011) L. H. Cisneros, J. O. Kessler, S. Ganguly, and R. E. Goldstein, Physical Review E 83, 061907 (2011).
  * Mishra  _et al._ (2012) S. Mishra, K. Tunstrøm, I. D. Couzin, and C. Huepe, Physical Review E 86, 011901 (2012).
  * Singh and Mishra (2020) J. P. Singh and S. Mishra, Physica A: Statistical Mechanics and its Applications 544, 123530 (2020).
  * Cates and Tailleur (2015) M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Tailleur and Cates (2008) J. Tailleur and M. Cates, Physical review letters 100, 218103 (2008).
  * Fily and Marchetti (2012) Y. Fily and M. C. Marchetti, Physical review letters 108, 235702 (2012).
  * Martin  _et al._ (2021) D. Martin, J. O’Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur, and F. van Wijland, Physical Review E 103, 032607 (2021).
  * McCandlish  _et al._ (2012) S. R. McCandlish, A. Baskaran, and M. F. Hagan, Soft Matter 8, 2527 (2012).
  * Barré  _et al._ (2015) J. Barré, R. Chétrite, M. Muratori, and F. Peruani, Journal of Statistical Physics 158, 589 (2015).
  * Sese-Sansa  _et al._ (2018) E. Sese-Sansa, I. Pagonabarraga, and D. Levis, EPL (Europhysics Letters) 124, 30004 (2018).
  * Van Der Linden  _et al._ (2019) M. N. Van Der Linden, L. C. Alexander, D. G. Aarts, and O. Dauchot, Physical review letters 123, 098001 (2019).
  * Grégoire and Chaté (2004) G. Grégoire and H. Chaté, Physical review letters 92, 025702 (2004).
  * Paoluzzi  _et al._ (2018) M. Paoluzzi, M. Leoni, and M. C. Marchetti, Physical Review E 98, 052603 (2018).
  * Paoluzzi  _et al._ (2020) M. Paoluzzi, M. Leoni, and M. C. Marchetti, Soft Matter 16, 6317 (2020).
  * Couzin  _et al._ (2002) I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks, Journal of theoretical biology 218, 1 (2002).
  * Grégoire  _et al._ (2003) G. Grégoire, H. Chaté, and Y. Tu, Physica D: Nonlinear Phenomena 181, 157 (2003).
  * (40) See supplemental material at http://link.aps.org/supplemental/10.1103/PhysRevE.105.044124 for numerical simulations, relevant length scales and timescales, and movies of the dynamics, which includes refs. [67,68] .
  * Solon  _et al._ (2015) A. P. Solon, J.-B. Caussin, D. Bartolo, H. Chaté, and J. Tailleur, Physical Review E 92, 062111 (2015).
  * Chaté (2020) H. Chaté, Annual Review of Condensed Matter Physics 11, 189 (2020).
  * Ginelli (2016) F. Ginelli, The European Physical Journal Special Topics 225, 2099 (2016).
  * Peruani  _et al._ (2010) F. Peruani, L. Schimansky-Geier, and M. Baer, The European Physical Journal Special Topics 191, 173 (2010).
  * Dossetti  _et al._ (2009) V. Dossetti, F. Sevilla, and V. Kenkre, Physical Review E 79, 051115 (2009).
  * Radicchi and Fortunato (2009) F. Radicchi and S. Fortunato, Physical review letters 103, 168701 (2009).
  * Czirók and Vicsek (2000) A. Czirók and T. Vicsek, Physica A: Statistical Mechanics and its Applications 281, 17 (2000).
  * Mertens and Moore (2012) S. Mertens and C. Moore, Physical Review E 86, 061109 (2012).
  * Peruani and Baer (2013) F. Peruani and M. Baer, New Journal of Physics 15, 065009 (2013).
  * Liu and Pandey (1997) Y. Liu and R. B. Pandey, Physical Review B 55, 8257 (1997).
  * Kyriakopoulos  _et al._ (2019) N. Kyriakopoulos, H. Chaté, and F. Ginelli, Physical Review E 100, 022606 (2019).
  * Ling  _et al._ (2019) H. Ling, G. E. Mclvor, J. Westley, K. van der Vaart, R. T. Vaughan, A. Thornton, and N. T. Ouellette, Nature communications 10, 1 (2019).
  * Leggett  _et al._ (2019) S. E. Leggett, Z. J. Neronha, D. Bhaskar, J. Y. Sim, T. M. Perdikari, and I. Y. Wong, Proceedings of the National Academy of Sciences 116, 17298 (2019).
  * Klamser  _et al._ (2021) P. P. Klamser, L. Gómez-Nava, T. Landgraf, J. W. Jolles, D. Bierbach, and P. Romanczuk, Frontiers in Physics 9 (2021).
  * Bäuerle  _et al._ (2018) T. Bäuerle, A. Fischer, T. Speck, and C. Bechinger, Nature communications 9, 3232 (2018).
  * Lavergne  _et al._ (2019) F. A. Lavergne, H. Wendehenne, T. Bäuerle, and C. Bechinger, Science 364, 70 (2019).
  * Soma  _et al._ (2020) R. Soma, B. Nakayama, M. Kuwahara, E. Yamamoto, and T. Saiki, Applied Physics Letters 117, 221601 (2020).
  * Bäuerle  _et al._ (2020) T. Bäuerle, R. C. Löffler, and C. Bechinger, Nature communications 11, 1 (2020).
  * Werfel  _et al._ (2014) J. Werfel, K. Petersen, and R. Nagpal, Science 343, 754 (2014).
  * Crosato  _et al._ (2019) E. Crosato, M. Prokopenko, and R. E. Spinney, Physical Review E 100, 042613 (2019).
  * Sesé-Sansa  _et al._ (2021) E. Sesé-Sansa, D. Levis, and I. Pagonabarraga, Physical Review E 104, 054611 (2021).
  * Chakraborty and Bhattacharya (2016) A. Chakraborty and K. Bhattacharya, EPL (Europhysics Letters) 116, 48001 (2016).
  * Dossetti and Sevilla (2015) V. Dossetti and F. J. Sevilla, Physical review letters 115, 058301 (2015).
  * Bhattacherjee  _et al._ (2014) B. Bhattacherjee, K. Bhattacharya, and S. S. Manna, Frontiers in Physics 1, 35 (2014).
  * Saleur and Duplantier (1987) H. Saleur and B. Duplantier, Physical review letters 58, 2325 (1987).
  * Stauffer and Aharony (2003) D. Stauffer and A. Aharony,  _Introduction to percolation theory: revised second edition_ (Taylor & Francis, Philadelphia, 2003).
  * Allen and Tildesley (1987) M. Allen and D. Tildesley,  _Computer Simulation of Liquids, Oxford: Clarendon Press_ (1987).
  * Hoshen and Kopelman (1976) J. Hoshen and R. Kopelman, Physical Review B 14, 3438 (1976).


# Charged Wormhole Solutions in 4D Einstein-Gauss-Bonnet Gravity

Piyachat Panyasiripan<EMAIL_ADDRESS>School of Science, Walailak
University, Nakhon Si Thammarat, 80160, Thailand  Fereshteh Felegary
<EMAIL_ADDRESS>School of Science, Walailak University, Nakhon Si
Thammarat, 80160, Thailand  Faculty of Physics, Shahrood University of
Technology, Shahrood, Iran  Phongpichit Channuie<EMAIL_ADDRESS>School of Science, Walailak University, Nakhon Si Thammarat, 80160, Thailand
College of Graduate Studies, Walailak University, Nakhon Si Thammarat, 80160,
Thailand

###### Abstract

In the present work, we construct models of static wormholes in 4-dimensional
Einstein-Gauss-Bonnet (4D EGB) gravity with an (an)isotropic energy momentum
tensor (EMT) and a Maxwell field as supporting matters for the wormhole
geometry assuming a constant redshift function. We obtain exact spherically
symmetric wormhole solutions in 4D EGB gravity for an isotropic and
anisotropic matter sources under the effect of electric charge. The solutions
are more generalized in the sense of incorporating the charge contributions.
Furthermore, we examine the null, weak and strong energy conditions at the
wormhole throat of radius $r=r_{0}$. We demonstrate that at the wormhole
throat, the classical energy conditions are violated by arbitrary small
amount. Additionally, we analyze the wormhole geometry incorporating
semiclassical corrections through embedding diagrams. We discover that the
wormhole solutions are in equilibrium as the anisotropic and hydrostatic
forces cancel each other.

##  I Introduction

The concept of wormholes has long fascinated scientists and the public alike,
offering tantalizing possibilities of shortcuts through spacetime, allowing
for rapid travel between distant regions of the universe. In the realm of
theoretical physics, wormholes emerge as solutions to Einstein’s field
equations of general relativity, representing bridges that connect separate
points in spacetime. The first theoretical exploration of wormholes attracted
considerable attention soon after the discovery of General Relativity, with
initial work by Flamm Flamm:1916 . This investigation was further advanced by
Einstein and Rosen Einstein:1935tc , who introduced the Einstein-Rosen bridge
to describe wormholes at that time. However, the early investigations into
wormholes were primarily theoretical exercises with no practical connection to
physical reality. This changed when Morris and Thorne conducted their study,
demonstrating the potential existence of wormholes and how they could be
utilized for interstellar travel Morris:1988cz .

Subsequently, inspired by the ideas of Morris and Thorne Morris:1988cz , there
has been intense activity in investigating wormholes within modified theories
of gravity and higher-dimensional gravity theories Harko:2013yb ;
Zangeneh:2014noa ; Banerjee:2020uyi ; Shaikh:2016dpl ; Bronnikov:2015pha .
Moreover, there has been research on wormholes in extended theories of gravity
DeFalco:2023twb ; Rosa:2018jwp ; Rosa:2022osy ; Rosa:2023guo ; Rosa:2021yym ;
Rosa:2023olc , focusing on developing astrophysical techniques to detect them
DeFalco:2020afv ; DeFalco:2021klh ; DeFalco:2021btn ; DeFalco:2023kqy , the
reconstruction of their solutions DeFalco:2021ksd , and their existence in
higher dimensions Deshpande:2022zfm . Among many other interesting scenarios,
wormhole solutions in Einstein-Gauss-Bonnet theory represent a fascinating
extension of the traditional wormhole concept by incorporating higher-
dimensional and higher-order curvature terms, see Refs.Bhawal:1992sz ;
Mehdizadeh:2020nrw ; Ilyas:2023rde ; Zubair:2023wjh ; Rani:2016xpa ;
Mehdizadeh:2015jra ; Dotti:2007az ; Maeda:2008nz . However, in 4D, the Gauss-
Bonnet (GB) term is a topological invariant and does not affect gravitational
dynamics. Recently, Glavan and Lin Glavan:2019inb  proposed a method to
regularize the Gauss-Bonnet equations dimensionally, resulting in a 4D metric
theory that circumvents Lovelock’s theorem and avoids Ostrogradsky
instability. Their approach involves formulating the equations in D
dimensions, rescaling the coupling constant $\alpha\to\alpha/(D-4)$, and then
taking the limit as $D\to 4$. This allows the GB term to contribute
nontrivially to gravitational dynamics, leading to what is known as 4D
Einstein-Gauss-Bonnet (EGB) gravity. This regularization process was initially
considered by Tomozawa Tomozawa:2011gp , who included finite one-loop quantum
corrections to Einstein gravity. Wormhole solutions were also investigated in
4D EGB gravity. For instance, Ref.Jusufi:2020yus  discussed wormhole solutions
in 4D EGB gravity for isotropic and anisotropic matter sources. In
Ref.Kokubu:2020lxs ; Zhang:2020kxz , the possibility for eternal existence of
thin-shell wormholes has been investigated. Moreover, the stability of thin-
shell wormholes Godani:2022jwz  and Yukawa-Casimir wormholes Shweta:2023hvm
has been investigated. Moreover, wormholes within galactic dark matter halos
has been recently explored in Ref.Hassan:2024xyx .

Among these fascinating structures, charged traversable wormholes stand out
due to their additional complexity and potential implications. A charged
wormhole, influenced by electromagnetic fields, must satisfy the Einstein-
Maxwell equations, which integrate the principles of general relativity with
electromagnetism. Therefore, the work of Ref.Jusufi:2020yus  can be extended
to include also the electromagnetic field as an additional source. Charged
wormholes appear in many theories of gravity, often providing intriguing
solutions that combine both electromagnetic and gravitational fields to
explore the stability, structure, and physical properties of these exotic
spacetime configurations. Regarding charged solutions, the authors of
Ref.Mehdizadeh:2018smu  study wormhole solutions in Einstein-Cartan theory
(ECT) in the presence of electric charge; while, the exact solutions of
wormhole with extra fields such as scalar field and electric charge have been
addressed in Ref.Kim:2001ri .

This work is organized as follows: We give a short review on charged EGB
gravity & wormhole solutions in Sec.II. Subsequently, the embedding diagrams
of the wormholes are illustrated in Sec.III. In Sec.IV, we study the energy
conditions of the charge wormhole solutions. We also use the stability
analysis to examine whether the wormhole solutions are in equilibrium in
Sec.V. Our conclusions are drawn in the last section.

##  II A short review on charged EGB gravity & Wormhole solutions

Recently, 4D Einstein–Gauss–Bonnet (4DEGB) gravity has attracted considerable
attention. Initially, we introduce Lovelock’s theorem as a foundational
concept, focusing particularly on Gauss–Bonnet terms within the gravity
action, which are extensively reviewed in Ref.Fernandes:2022zrq . However, our
discussion closely follows the approach outlined in Ref.Glavan:2019inb .
Starting with the standard Einstein-Gauss-Bonnet action without a cosmological
constant, we incorporate the Maxwell field, modifying the Gauss-Bonnet
coupling constant $\alpha\rightarrow\alpha/(D-4)$, as described in
Ref.Fernandes:2020rpa .

$S=\int
d^{D}x\sqrt{-g}\left[\frac{1}{16\pi}R+\frac{\alpha}{D-4}\mathcal{L}_{\text{GB}}-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}\right]+\mathcal{S}_{\text{matter}},$
(1)

where $g$ denotes the determinant of the metric $g_{\mu\nu}$ and $\alpha$ is
the Gauss-Bonnet coupling coefficient with dimension $[\rm length]^{2}$ and we
assume that it is positive here,
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the usual Maxwell
tensor. The term $\mathcal{L}_{\text{GB}}$ is the Lagrangian defined by

$\mathcal{L}_{\text{GB}}=R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}-4R^{\mu\nu}R_{\mu\nu}+R^{2}.$
(2)

In Eq.(1), $S_{\rm matter}$ represents the matter fields within the theory. By
varying this action with respect to the metric $g_{\mu\nu}$ and the field
strength tensor $F_{\mu\nu}$, the following field equations can be derived
EslamPanah:2020hoj :

$\displaystyle G_{\mu\nu}+\frac{\alpha}{D-4}H_{\mu\nu}$ $\displaystyle=$
$\displaystyle 8\pi T_{\mu\nu},$ (3) $\displaystyle\nabla_{\mu}F_{\mu\nu}$
$\displaystyle=$ $\displaystyle 0\,.$ (4)

where $G_{\mu\nu}$ is the usual Einstein tensor, $T_{\mu\nu}=T^{\rm
m}_{\mu\nu}+T^{\rm em}_{\mu\nu}$, $T^{\rm
m}_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}\mathcal{S}_{m}\right)}{\delta
g^{\mu\nu}}$ is the energy momentum tensor of matter, while $H_{\mu\nu}$ and
$T^{\rm em}_{\mu\nu}$ read

$\displaystyle H_{\mu\nu}$ $\displaystyle=$ $\displaystyle
2\Bigr{(}RR_{\mu\nu}-2R_{\mu\sigma}{R}{{}^{\sigma}}_{\nu}-2R_{\mu\sigma\nu\rho}{R}^{\sigma\rho}-R_{\mu\sigma\rho\delta}{R}^{\sigma\rho\delta}{{}_{\nu}}\Bigl{)}-\frac{1}{2}g_{\mu\nu}\mathcal{L}_{\text{GB}},$
(5) $\displaystyle T^{\rm em}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\Big{(}F^{\lambda}_{\mu}F_{\nu\lambda}-\frac{1}{4}g_{\mu\nu}F^{\alpha\beta}F_{\alpha\beta}\Big{)}\,.$
(6)

Here, $R$ denotes the Ricci scalar, $R_{\mu\nu}$ the Ricci tensor,
$H_{\mu\nu}$ is the Lanczos tensor, and $R_{\mu\sigma\nu\rho}$ represents the
Riemann tensor. It is straightforward to demonstrate that the electromagnetic
stress–energy tensor, $T^{\rm em}_{\mu\nu}$, is traceless. To construct a
wormhole within the framework of charged 4D Einstein-Gauss-Bonnet (EGB)
gravity, we consider the general static, spherically symmetric $D$-dimensional
metric discussed in Mehdizadeh:2015jra ; Zubair:2023wjh , given by

$ds^{2}=-e^{2\Phi(r)}dt^{2}+\frac{dr^{2}}{1-\frac{b(r)}{r}}+r^{2}d\Omega_{D-2}^{2}.$
(7)

where

$d\Omega^{2}_{D-2}=d\theta^{2}_{1}+\sum^{D-2}_{i=2}\prod^{i-1}_{j=1}\sin^{2}\theta_{j}\;d\theta^{2}_{i}\;,$

where the function $\Phi(r)$ represents the redshift function of an infalling
body, which must be finite everywhere to prevent the presence of an event
horizon. Conversely, $b(r)$ denotes the spatial shape function of the wormhole
geometry, determining its shape in the embedding diagram Mehdizadeh:2015jra .
It is crucial that $b(r)$ satisfies the boundary condition $b(r=r_{0})=r_{0}$
at the throat $r_{0}$, where $r_{0}\leq r\leq\infty$. Moreover, to ensure the
traversability of the wormhole, the function $b(r)$ must adhere to the
flaring-out condition, which is derived from the embedding calculation and is
expressed as follows:

$\frac{b(r)-rb^{\prime}(r)}{b^{2}(r)}>0.$ (8)

We consider an anisotropic fluid for the matter source defined by the stress
energy tensor

$\displaystyle T^{\nu}_{i}$ $\displaystyle=$
$\displaystyle(\rho+p_{t})u^{\nu}u_{i}+p_{t}g^{\nu}_{i}+(p_{r}-p_{t})\chi_{i}\chi^{\nu},$
(9)

where, $u_{\nu}$ represents the 4-velocity and $\chi_{\nu}$ denotes the unit
spacelike vector in the radial direction, where $\rho(r)$ stands for the
energy density, and $p_{r}(r)$ and $p_{t}(r)$ denote the radial and transverse
pressures, respectively. The field equations in Einstein-Gauss-Bonnet (EGB)
gravity can be derived using the metric given in Eq. (7), as detailed in
Zubair:2023wjh . The corresponding electromagnetic potential is expressed as:

$\displaystyle A^{\mu}=-\frac{q}{r}dt\delta^{\mu}_{t}\,.$ (10)

The only non-vanishing components of the Faraday tensor take the form:

$\displaystyle F_{01}=-F_{10}=-\frac{q}{r^{2}},$ (11)

where $q$ denotes an integration constant associated with the electric charge.
By employing the metric given in Eq. (7) along with the stress tensor from Eq.
(9), in the limit as $D\to 4$, the components of the field equations (3) can
be expressed as:

$\displaystyle 8\pi\Big{(}\rho(r)+\frac{q^{2}}{8\pi
r^{4}}\Big{)}=\frac{b(r)}{r^{3}}\Big{(}1-\frac{\alpha
b(r)}{r^{3}}\Big{)}+\frac{1}{r^{2}}\left(\frac{2\alpha
b(r)}{r^{3}}+1\right)\left(b^{\prime}(r)-\frac{b(r)}{r}\right)\,,$ (12)
$\displaystyle 8\pi\Big{(}p_{r}(r)-\frac{q^{2}}{8\pi
r^{4}}\Big{)}=-\frac{b(r)\left(1-\frac{\alpha
b(r)}{r^{3}}\right)}{r^{3}}+\frac{2}{r}\left(1-\frac{b(r)}{r}\right)\left(\frac{2\alpha
b(r)}{r^{3}}+1\right)\Phi^{\prime}(r)\,,$ (13)

and

$\displaystyle 8\pi\Big{(}p_{t}(r)+\frac{q^{2}}{8\pi r^{4}}\Big{)}$
$\displaystyle=$ $\displaystyle-\frac{\alpha
b(r)^{2}}{r^{6}}+\frac{2\alpha}{r^{4}}\left(1-\frac{b(r)}{r}\right)\left(b(r)-rb^{\prime}(r)\right)\Phi^{\prime}(r)$
(14) $\displaystyle+$
$\displaystyle\left(1-\frac{b(r)}{r}\right)\left(\frac{2\alpha
b(r)}{r^{3}}+1\right)\left(\frac{\left(b(r)-rb^{\prime}(r)\right)\Phi^{\prime}(r)}{2r(r-b(r))}+\Phi^{\prime\prime}(r)-\Phi^{\prime}(r)^{2}\right)$
$\displaystyle+$
$\displaystyle\left(1-\frac{b(r)}{r}\right)\left(1-\frac{2\alpha
b(r)}{r^{3}}\right)\left(\frac{b(r)-rb^{\prime}(r)}{2r^{2}(r-b(r))}-\frac{\Phi^{\prime}(r)}{r}\right)\,,$

where the prime denotes a derivative with respect to the radial coordinate
$r$. In this context, we have five unknown functions of $r$, i.e., $\rho(r)$,
$p_{r}(r)$, $p_{t}(r)$, $b(r)$ and $\Phi(r)$. In the following, we examine
wormhole solutions in 4D Einstein-Gauss-Bonnet Gravity.

###  II.1 Isotropic solution

To simplify the problem, we consider a constant redshift function,
$\Phi(r)=\text{const.}$, implying a wormhole solution with zero tidal force.
This choice simplifies calculations and yields intriguing exact wormhole
solutions. Recent studies Fernandes:2020nbq ; Hennigar:2020lsl  have
demonstrated that by taking the trace of the field equations (3), one obtains
the following simple form:

$R+\frac{\alpha}{2}\mathcal{L}_{\text{GB}}=-8\pi T,$ (15)

where the trace $T=T^{\mu}_{\nu}$. Using the metric form (7) for isotropic
fluid matter we can use the relation $p_{r}(r)=p_{t}(r)=\omega\rho(r)$ obtain
the following condition

$\frac{2b^{\prime}(r)}{r^{2}}+8\pi\left[-\rho(r)+3\omega\rho(r)\right]=0.$
(16)

Utilizing Eq.(16) and (12), we obtain a differential equation for the shape
function as

$\displaystyle\frac{r^{4}(3\omega+1)b^{\prime}(r)}{3\omega-1}+\alpha
b(r)\left(2rb^{\prime}(r)-3b(r)\right)-q^{2}r^{2}=0.$ (17)

Solving the above equation for $b(r)$ which has the following form

$b(r)=-\frac{r^{3}(3\omega+1)}{2\alpha(3\omega-1)}\Bigg{(}1\pm\sqrt{1+\frac{4\alpha{\cal
A}}{r^{3}r_{0}}}\Bigg{)},$ (18)

where we have defined a new parameter ${\cal A}$ as

${\cal
A}=\Big{(}\frac{q^{2}(3\omega-1)^{2}(r-r_{0})}{r(3\omega+1)^{2}}+\frac{(3\omega-1)\left(-\alpha+3\omega\left(\alpha+r_{0}^{2}\right)+r_{0}^{2}\right)}{(3\omega+1)^{2}}\Big{)}/,,$
(19)

where we have applied the condition $b(r=r_{0})=r_{0}$, noting that when
$q=0$, our results reduce to those presented in Ref.Jusufi:2020yus . It’s
important to emphasize that the solution given by Eq.(18) is valid only for
$\omega\neq\pm 1/3$. The $\pm$ sign in Eq.(18) denotes two distinct branches
of solutions. In studies such as Boulware:1985wk ; Wheeler:1985nh , it has
been demonstrated that Einstein-Gauss-Bonnet (EGB) black holes with the $+$
branch are unstable, exhibiting a ghost graviton degree of freedom, whereas
those with the $-$ branch are stable and free of ghosts. In our scenario, in
the limit $\alpha\to 0$, the positive branch leads to

$\frac{b(r)}{r}=-\frac{r^{2}(3\omega+1)}{(3\omega-1)\alpha}-\frac{r_{0}}{r}-\frac{q^{2}(3\omega-1)(r-r_{0})}{r^{2}r_{0}(3\omega+1)}\ldots$
(20)

which represents a wormhole solution within de-Sitter/anti-de Sitter
spacetimes, depending on the sign of $\alpha$. Conversely, in the limit as
$\alpha\to 0$, the negative branch transitions to

$\frac{b(r)}{r}=\frac{r_{0}}{r}+\frac{q^{2}(3\omega-1)(r-r_{0})}{r^{2}r_{0}(3\omega+1)}+\ldots,$
(21)

and the standard Morris-Thorne wormhole is obtained when the charge vanishes.
Henceforth we restrict to the negative branch in which case the 4D EGB
wormhole metric reads

$ds^{2}=-dt^{2}+\frac{dr^{2}}{1+\frac{r^{2}(3\omega+1)}{2\alpha(3\omega-1)}\left(1\pm\sqrt{1+\frac{4\alpha\mathcal{A}}{r^{3}r_{0}}}\right)}+r^{2}d\Omega^{2}_{2},$
(22)

where the unit sphere line element is given by
$d\Omega_{2}^{2}=d\theta^{2}+\sin^{2}\theta,d\phi^{2}$. It’s important to note
that the above analysis holds true for the interval $\omega\geq-1/3$ and
$\omega\neq 1/3$. For values $\omega<-1/3$, the sign of the solutions
reverses.

Figure 1: Figures show the behaviors of $b(r)$ and $b^{\prime}(r)$ as a
function of $r$ for our exact isotropic wormhole. We have used the constants
$\alpha=0.001$ and $r_{0}=1$ along with $\omega=0.5$.

To obtain a TW we need to compute the redshift function. We have to impose the
EoS $p_{r}(r)=\omega\rho(r)$. From Eq.(12) & (13), one finds

$\displaystyle\Phi^{\prime}(r)=\frac{r^{4}\omega b^{\prime}(r)+2\alpha r\omega
b(r)b^{\prime}(r)+r^{3}b(r)-3\alpha\omega b(r)^{2}-\alpha
b(r)^{2}-q^{2}r^{2}}{2r(r-b(r))\left(2\alpha b(r)+r^{3}\right)}\,.$ (23)

Substituting $b(r)$ from Eq.(21), we come up with

$\displaystyle\Phi^{\prime}(r)$ $\displaystyle=$
$\displaystyle\frac{3q^{2}r\omega-8q^{2}r_{0}\omega+q^{2}(-r)+3rr_{0}^{2}\omega+rr_{0}^{2}}{2r\left(r-r_{0}\right)\left(-3q^{2}\omega+q^{2}+3rr_{0}\omega+rr_{0}\right)}$
(24)
$\displaystyle-\frac{\alpha}{2\left(r^{5}\left(r-r_{0}\right)r_{0}(3\omega+1)\left(-3q^{2}\omega+q^{2}+3rr_{0}\omega+rr_{0}\right)\right)}\times$
$\displaystyle\times\Big{(}\big{(}3q^{2}r\omega-3q^{2}r_{0}\omega+q^{2}(-r)+q^{2}r_{0}+3rr_{0}^{2}\omega+rr_{0}^{2}\big{)}$
$\displaystyle\quad\big{(}9q^{2}r\omega^{2}-15q^{2}r_{0}\omega^{2}+6q^{2}r\omega-14q^{2}r_{0}\omega-3q^{2}r+q^{2}r_{0}+9rr_{0}^{2}\omega^{2}+12rr_{0}^{2}\omega+3rr_{0}^{2}\big{)}\Big{)}$
$\displaystyle+O\left(\alpha^{2}\right)\,.$

Let us consider solely the first term, we have

$\displaystyle\Phi^{\prime}(r)\simeq\frac{3q^{2}r\omega-8q^{2}r_{0}\omega+q^{2}(-r)+3rr_{0}^{2}\omega+rr_{0}^{2}}{2r\left(r-r_{0}\right)\left(-3q^{2}\omega+q^{2}+3rr_{0}\omega+rr_{0}\right)}\,.$
(25)

We choose $\omega$ such that $\Phi^{\prime}(r)=0$ to avoid the appearance of a
horizon. This can be seen if

$\displaystyle\omega\simeq\frac{r\left(q^{2}-r_{0}^{2}\right)}{3q^{2}r-8q^{2}r_{0}+3rr_{0}^{2}}\,.$
(26)

We next assume the small charge approximation and expand (26) to obtain

$\displaystyle\omega\simeq-\frac{1}{3}+\frac{2q^{2}(3r-4r_{0})}{9rr_{0}^{2}}+{\cal
O}\left(q^{3}\right)\,.$ (27)

Interestingly, it is worth mentioning here that due to the present of a charge
contribution, $\omega$ naturally deviates from $-\tfrac{1}{3}$.

###  II.2 Anisotropic fluid scenario

We shall argue that there are two ways to find an anisotropic solution. Taking
into account Eqs. (13-14) along with an interesting EoS Moraes:2017dbs ,

$p_{t}(r)=\omega_{t}p_{r}(r),$ (28)

leads to the following expression

$\displaystyle r\left(r^{3}-2\alpha
b(r)\right)b^{\prime}(r)+b(r)\left(2\alpha(\omega_{t}+2)b(r)-r^{3}(2\omega_{t}+1)\right)+2q^{2}r^{2}(\omega_{t}+1)=0,$

where the EoS parameter $\omega_{t}$ and the redshift function $\Phi(r)$ are
constants. Solving the last differential equation we get the following shape
function

$b(r)=\frac{r^{3}}{2\alpha}\left[1\pm\sqrt{\frac{r^{4}-4\alpha
q^{2}}{r^{4}}+4\alpha^{2}c_{1}r^{2\omega_{t}-2}}\right],$ (29)

where $c_{1}$ is a constant of integration, which can be fixed using the
condition $b(r=r_{0})=r_{0}$, Then the shape function modifies to

$b(r)=\frac{r^{3}}{2\alpha}\left[1\pm\sqrt{1-\frac{4\alpha
q^{2}}{r^{4}}+\frac{4\alpha
r^{2\left(\omega_{t}-1\right)}}{r_{0}^{2(\omega_{t}+1)}}\left(\alpha+q^{2}-r_{0}^{2}\right)}\right].$
(30)

Similarly for the $+$ve branch sign in the limit $\alpha\to 0$, we obtain

$\frac{b(r)}{r}=\frac{r^{2}}{\alpha}+\frac{r^{2\omega_{t}}}{r_{0}^{2(\omega_{t}+1)}}\left(q^{2}-r_{0}^{2}\right)-\frac{q^{2}}{r^{2}}+...$
(31)

which is a wormhole solution in a de-Sitter/ anti-de Sitter spacetimes
depending on the sign of $\alpha$. On the other hand, in the limit $\alpha\to
0$, the $-$ve goes over

$\frac{b(r)}{r}=-\frac{r^{2\omega_{t}}}{r_{0}^{2(\omega_{t}+1)}}\left(q^{2}-r_{0}^{2}\right)+\frac{q^{2}}{r^{2}}+\ldots,$
(32)

which reduces to the standard Morris-Thorne wormhole. Henceforth, we restrict
to the $-$ve branch for anisotropic fluid in the 4D EGB takes the form

$ds^{2}=-dt^{2}+\frac{dr^{2}}{1-\frac{r^{2}}{2\alpha}\left[1-\sqrt{1-\frac{4\alpha
q^{2}}{r^{4}}+\frac{4\alpha
r^{2\left(\omega_{t}-1\right)}}{r_{0}^{2(\omega_{t}+1)}}\left(\alpha+q^{2}-r_{0}^{2}\right)}\right]}+r^{2}d\Omega_{2}^{2}.$
(33)

Figure 2: Figures show the behaviors of $b(r)$ and $b^{\prime}(r)$ as a
function of $r$ for an isotropic wormhole. We have used the constants
$\alpha=0.001$ and $r_{0}=1$ along with $\omega_{t}=-1/3$.

###  II.3 Model with a specific energy density

To determine the shape function, one can consider a specific energy density
profile. Here, we will adopt the following energy density Lobo:2012qq :

$\rho(r)=\rho_{0}\left(\frac{r_{0}}{r}\right)^{\beta},$ (34)

where $\beta$ and $\rho_{0}$ are constants. Solving the Eq. (12), we obtain
the following shape function

$b(r)=-\frac{r^{3}}{\alpha}\left(1\pm\sqrt{1+4\alpha\left(-\frac{8\pi\rho_{0}\left(\frac{r_{0}}{r}\right)^{\beta}}{\beta-3}+\frac{-q^{2}+\alpha
c_{1}r}{r^{4}}\right)}\right),$ (35)

in which $c_{1}$ is a constant of integration. Using the condition
$b(r=r_{0})=r_{0}$, we obtain

$c_{1}=\frac{\alpha+q^{2}+r_{0}^{2}}{\alpha
r_{0}}+\frac{8\pi\rho_{0}r_{0}^{4}}{\alpha(\beta-3)r_{0}}.$ (36)

Substituting the above expression into the Eq. (35), one obtains the two
different branches of the solution, which are

$b(r)=-\frac{r^{3}}{\alpha}\left(1\pm\sqrt{1-\frac{4\alpha
q^{2}}{r^{4}}+\frac{4\alpha}{r^{3}r_{0}}\left(\alpha+q^{2}+\frac{8\pi\rho_{0}r_{0}^{4}}{\beta-3}+r_{0}^{2}\right)-\frac{32\pi\alpha\rho_{0}}{\beta-3}\left(\frac{r_{0}}{r}\right)^{\beta}}\right).$
(37)

The branch with the negative sign is stable and we select it as the physical
solution, under the condition $\beta\geq 4$. However, when $\beta=3$, there
appears to be a singularity. It’s important to note that there is a third case
where the sign of the solution flips, specifically when $\beta=1$ or $2$. The
negative branch solutions for $\beta\geq 4$ are asymptotically flat; for
instance, choosing $\beta=4$ illustrates this behavior in the limit $\alpha\to
0$, where the negative branch gives

$\frac{b(r)}{r}=-\frac{q^{2}}{r^{2}}\Big{(}1-\frac{r}{r_{0}}\Big{)}-\frac{8\pi\rho_{0}r^{4}\left(\frac{r_{0}}{r}\right)^{\beta}+rr_{0}\left(\beta+8\pi\rho_{0}r_{0}^{2}-3\right)}{(\beta-3)r^{2}}+\ldots$
(38)

On the other hand, in the limit $\alpha\to 0$, the positive branch goes over

$\frac{b(r)}{r}=-\frac{r^{2}}{\alpha}+\frac{q^{2}}{r^{2}}\Big{(}1-\frac{r}{r_{0}}\Big{)}+\frac{8\pi\rho_{0}r^{4}\left(\frac{r_{0}}{r}\right)^{\beta}-rr_{0}\left(\beta+8\pi\rho_{0}r_{0}^{2}-3\right)}{(\beta-3)r^{2}}+\ldots,$
(39)

corresponds to the wormhole solution in de-Sitter/anti-de Sitter spacetimes.
Finally, it can be verified that for $\beta=1,2$, our solutions (37) are not
asymptotically flat. Using the shape function (37), the wormhole metric is
defined by the expression:

$ds^{2}=-dt^{2}+\frac{dr^{2}}{1+\frac{r^{2}}{\alpha}\left(1-\sqrt{1+\zeta}\right)}+r^{2}d\Omega^{2}_{2},$
(40)

where

$\zeta=-\frac{4\alpha
q^{2}}{r^{4}}+\frac{4\alpha}{r^{3}r_{0}}\left(\alpha+q^{2}+\frac{8\pi\rho_{0}r_{0}^{4}}{\beta-3}+r_{0}^{2}\right)-\frac{32\pi\alpha\rho_{0}}{\beta-3}\left(\frac{r_{0}}{r}\right)^{\beta},$

under the condition $\beta\geq 4$. However, as previously noted, caution is
required since for $\beta=1,2$, the sign of the solution flips. As a special
case, one can consider the limit $\beta\to 0$, where the energy density
becomes constant, i.e., $\rho=\rho_{0}$.

Figure 3: Figures show the behaviors of $b(r)$ and $b^{\prime}(r)$ as a
function of $r$ for a model with a specific energy density. We have used the
constants $\alpha=0.001,\,\beta=4,\,\rho_{0}=0.01,\,q=0.5,\,r_{0}=1$.

##  III Embedding diagram

In this subsection, we examine the embedding diagrams to depict the charged
wormhole solutions in 4D Einstein-Gauss-Bonnet Gravity, focusing on an
equatorial slice $\theta=\pi/2$ at a fixed moment in time $t={\rm constant}$.
To achieve this, we consider the metric expressed as

$\displaystyle ds^{2}=\frac{dr^{2}}{1-\frac{b(r)}{r}}+r^{2}d\phi^{2}\,.$ (41)

We embed the metric (7) into three-dimensional Euclidean space to visualize
this slice. The spacetime can then be described using cylindrical coordinates
as

$\displaystyle ds^{2}=dz^{2}+dr^{2}+r^{2}d\phi^{2}\,.$ (42)

From the above equations, we can write

$\displaystyle\frac{dz}{dr}=\pm\sqrt{\frac{r}{r-b(r)}-1}\,,$ (43)

where $b(r)$ is given by Eq.(18), Eq.(30) and Eq.(LABEL:Anla). However, in
this section, we only focus on exploring the geometrical properties of these
matrices Eq.(18) and Eq.(30) via the embedding diagram. Invoking numerical
techniques allows us to illustrate the wormhole shape given in Fig.4.

Figure 4: The embedding diagram of wormholes geometry along the equatorial
plane ($t={\rm const.},\theta=\pi/2$). The specific case of a constant
redshift $\Phi^{\prime}(r)=0$ with $\alpha=0.001,q=0.3$ and $r_{0}=1$, we have
drawn embedding diagram for isotropic case and an anisotropic wormhole
illustrated on the left and right panel, respectively. We have considered the
numerical values of $\omega=0.5$ and $\omega_{t}=-1/3$ and $\beta=4$.

##  IV Energy conditions

In this section, we provide a detailed description of the energy conditions,
which are sets of inequalities based on the energy-momentum tensor.
Specifically, we begin by finding wormhole solutions for the weak energy
condition (WEC), defined as $T_{\mu\nu}U^{\mu}U^{\nu}$, where $U^{\mu}$ is a
timelike vector. For the given diagonal energy-momentum tensor, the WEC
implies:

$\displaystyle\rho(r)\geq 0\quad{\rm and}\quad\rho(r)+P_{i}(r)\geq 0\,.$ (44)

Next, the null energy condition (NEC) is defined by
$T_{\mu\nu}k^{\mu}k^{\nu}$, where $k^{\mu}$ is a null vector. For the given
diagonal energy-momentum tensor, the NEC implies:

$\displaystyle\rho(r)\geq 0\quad{\rm and}\quad\rho(r)+P_{i}(r)\geq 0\,.$ (45)

The strong energy condition (SEC) asserts that
$T_{\mu\nu}U^{\mu}U^{\nu}\geq\tfrac{1}{2}Tg_{\mu\nu}U^{\mu}U^{\nu}$ for any
timelike vector $U^{\mu}$, ensuring that gravity is attractive. The SEC
implies:

$\displaystyle\rho(r)+\sum_{i}P_{i}(r)\geq 0\quad{\rm
and}\quad\rho(r)+P_{i}(r)\geq 0\,.$ (46)

Note that while the WEC or SEC imply the NEC, any violation of the NEC also
indicates a violation of the SEC and WEC. In the following, we present the
energy conditions for all specific cases.

###  IV.1 Isotropic solution

We begin with the isotropic solution. Using the key Eqs.(12-14), we can
compute $\rho$ and $p_{r}$ to obtain

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle
B^{-1}\Big{(}3(A-1)r_{0}r^{4}(3\omega+1)^{2}-6\alpha
r(3\omega-1)\big{(}-\alpha+3\omega\big{(}\alpha+q^{2}+r_{0}^{2}\big{)}$ (47)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-q^{2}+r_{0}^{2}\big{)}+4\alpha
q^{2}r_{0}(1-3\omega)^{2}\Big{)}\,,$ $\displaystyle p_{r}$ $\displaystyle=$
$\displaystyle B^{-1}\Big{(}3(A-1)r_{0}r^{4}(3\omega+1)^{2}-6\alpha
r(3\omega-1)\big{(}-\alpha+3\omega\big{(}\alpha+q^{2}+r_{0}^{2}\big{)}$ (48)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-q^{2}+r_{0}^{2}\big{)}+4\alpha
q^{2}r_{0}(1-3\omega)^{2}\Big{)}\,,$

where

$\displaystyle A$ $\displaystyle\equiv$
$\displaystyle\sqrt{1+\frac{4\alpha(3\omega-1)\left(q^{2}\left(r-r_{0}\right)(3\omega-1)+r\left(-\alpha+3\omega\left(\alpha+r_{0}^{2}\right)+r_{0}^{2}\right)\right)}{r^{4}r_{0}(3\omega+1)^{2}}},$
(49) $\displaystyle B$ $\displaystyle\equiv$ $\displaystyle 8\pi\alpha
Ar^{4}r_{0}(1-3\omega)^{2}(3\omega+1)\,.$ (50)

In particular, it is evident from Fig.(6) that for $\alpha>0$, the NEC, and
consequently the WEC, are violated at the throat due to the flaring-out
condition.

Figure 5: Figures show that the behavior of NEC and SEC diagrams have been
plotted for the isotropic wormhole against $r$. Here also NEC is violated,
this implies WEC is also violated. For plotting the constants are
$\alpha=0.001$ and $r_{0}=1$ along with $\omega=0.5$.

###  IV.2 Anisotropic solution

Using the shape function provided in Eq.(30), we can compute the energy-
momentum components. Specifically, we determine the energy density and the
radial component as follows:

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\frac{1}{8\pi\alpha
Br^{4}}\Big{(}r_{0}^{-2(\omega_{t}+1)}\big{(}\alpha(2(B-2)\omega_{t}+B-8)r^{2\omega_{t}+2}\big{(}\alpha+q^{2}-r_{0}^{2}\big{)}$
(51) $\displaystyle\quad\quad\quad\quad+3(B-1)r^{4}r_{0}^{2\omega+2}+4\alpha
q^{2}r_{0}^{2\omega_{t}+2}\big{)}\Big{)}\,,$

and for the tangential component

$\displaystyle
p_{t}=\frac{\omega_{t}}{8\pi}\Big{(}r^{2(\omega_{t}-1)}r_{0}^{-2(\omega_{t}+1)}\left(\alpha+q^{2}-r_{0}^{2}\right)\Big{)}\,,$
(52)

where

$\displaystyle B\equiv\sqrt{1+4\alpha
r^{2(\omega_{t}-1)}r_{0}^{-2(\omega_{t}+1)}\left(\alpha+q^{2}-r_{0}^{2}\right)-\frac{4\alpha
q^{2}}{r^{4}}}.$ (53)

In the plots shown in Fig.(6), we observe that at the wormhole throat,
$r=r_{0}$, the energy condition $(\rho+p_{r})|_{r=r_{0}}<0$ and the condition
$(\rho+p_{r}+2p_{t})|_{r=r_{0}}<0$ are satisfied by arbitrarily small values.

Figure 6: Figures show that the behavior of NEC and SEC diagrams have been
plotted for the anisotropic wormhole against $r$. Here also NEC is violated,
this implies WEC is also violated. For plotting the constants are
$\alpha=0.001$ and $r_{0}=1$ along with $\omega_{t}=-1/3$.

###  IV.3 Model with a specific energy density

Using the shape function provided in Eq.(LABEL:Anla), we can compute the
energy-momentum components. Specifically, we determine the energy density and
the radial component as follows:

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\rho_{0}\left(\frac{r_{0}}{r}\right)^{\beta}\,,$ (54)
$\displaystyle p_{r}$ $\displaystyle=$
$\displaystyle\frac{1}{8\pi}\Bigg{(}\frac{1}{\alpha}+\frac{\alpha+q^{2}+r_{0}^{2}}{r^{3}r_{0}}-\frac{1}{\alpha}{\cal
B}(q,\alpha,\rho_{0})+\frac{8\pi\rho_{0}}{(\beta-3)r^{3}}\left(\rho_{0}^{3}-r^{3}\left(\frac{\rho_{0}}{r}\right)^{\beta}\right)\Bigg{)}\,,$
(55)

and for the tangential component

$\displaystyle p_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi\alpha(\beta-3)Br^{3}r_{0}}\Bigg{(}2r^{3}r_{0}\Bigg{(}-8\pi\alpha\rho_{0}\left(\frac{r_{0}}{r}\right)^{\beta}\Bigg{(}\sqrt{\frac{{\cal
C}(q,\alpha,\rho_{0})}{(\beta-3)r^{4}r_{0}}}-4\Bigg{)}$ (56)
$\displaystyle+\beta\Bigg{(}4\pi\alpha\rho_{0}\left(\frac{r_{0}}{r}\right)^{\beta}\Bigg{(}\sqrt{\frac{{\cal
C}(q,\alpha,\rho_{0})}{(\beta-3)r^{4}r_{0}}}-2\Bigg{)}+\sqrt{\frac{{\cal
C}(q,\alpha,\rho_{0})}{(\beta-3)r^{4}r_{0}}}-1\Bigg{)}-3\sqrt{\frac{{\cal
C}(q,\alpha,\rho_{0})}{(\beta-3)r^{4}r_{0}}}+3\Bigg{)}$
$\displaystyle-\alpha\left((\beta-3)\left(\alpha+q^{2}+r_{0}^{2}\right)+8\pi\rho_{0}r_{0}^{4}\right)\left(\sqrt{\frac{{\cal
C}(q,\alpha,\rho_{0})}{(\beta-3)r^{4}r_{0}}}+2\right)\Bigg{)}\,,$

where we have defined new functions

$\displaystyle{\cal B}(q,\alpha,\rho_{0})$ $\displaystyle\equiv$
$\displaystyle\sqrt{-\frac{4\alpha
q^{2}}{r^{4}}+\frac{4\alpha\left(\alpha+q^{2}+\frac{8\pi\rho_{0}r_{0}^{4}}{\beta-3}+r_{0}^{2}\right)}{r^{3}r_{0}}-\frac{32\pi\alpha\rho_{0}\left(\frac{r_{0}}{r}\right)^{\beta}}{\beta-3}+1}\,,$
(57) $\displaystyle{\cal C}(q,\alpha,\rho_{0})$ $\displaystyle\equiv$
$\displaystyle 4\alpha
r\left((\beta-3)\left(\alpha+q^{2}+r_{0}^{2}\right)+8\pi\rho_{0}r_{0}^{4}\right)-4\alpha(\beta-3)q^{2}r_{0}$
(58)
$\displaystyle+r^{4}r_{0}\left(\beta-32\pi\alpha\rho_{0}\left(\frac{r_{0}}{r}\right)^{\beta}-3\right).$

Figure 7: Figures show that the behavior of NEC and SEC diagrams have been
plotted for the anisotropic wormhole against $r$. Here also NEC is violated,
this implies WEC is also violated. For plotting the constants are
$\alpha=0.001$ and $r_{0}=1$ along with $\beta=4$.

In the plots shown in Fig.(7), we also find that at the wormhole throat,
$r=r_{0}$, the energy condition $(\rho+p_{r})|_{r=r_{0}}<0$ and the condition
$(\rho+p_{r}+2p_{t})|_{r=r_{0}}<0$ are satisfied by arbitrarily small values.

##  V Stability analysis

The generalized Tolman–Oppenheimer–Volkov (TOV) equation extends the
traditional TOV equation to include anisotropic pressures and exotic forms of
matter. This can help examine aspects of wormhole stability, although it might
not cover all unique features. By applying the generalized TOV equation to a
wormhole, one can analyze equilibrium conditions under gravitational forces
and pressure gradients from exotic matter. The equilibrium condition is given
by the generalized TOV equation, expressed as follows:

$\displaystyle\frac{dp_{r}}{dr}+\frac{1}{2}\big{(}\rho+p_{r}\big{)}\Phi^{\prime}+\frac{2}{r}\big{(}p_{r}-p_{t}\big{)}=0\,.$
(59)

This helps in understanding the balance required to maintain a stable
wormhole. For detailed study, we refer to Refs.Oppenheimer:1939ne ;
Tolman:1939jz ; Sokoliuk:2021rtv ; Dutta:2023wfg . For our purpose, we recast
the above equation as

$\displaystyle F_{a}+F_{g}+F_{h}=0\,,$ (60)

which provides the equilibrium condition for the wormholes. Here we have
defined new quantities:

$\displaystyle F_{h}=-\frac{dp_{r}}{dr}\,,\quad\quad
F_{g}=-\frac{1}{2}\big{(}\rho+p_{r}\big{)}\Phi^{\prime}\,,\quad\quad
F_{a}=\frac{2}{r}\big{(}p_{t}-p_{r}\big{)}\,.$ (61)

Here, $F_{a}$ represents the force resulting from the anisotropic matter of
the wormhole, $F_{g}$ denotes the gravitational force, and $F_{h}$ is the
hydrostatic force. $F_{a}$ arises due to modifications in the gravitational
Lagrangian of the Einstein-Hilbert action. However, as indicated by Eq.(60),
it is evident that for the system to remain in equilibrium, the sum of these
three forces must equate to zero.

In the anisotropic solution, the behaviors of $F_{a},\,F_{g}$, and $F_{h}$ are
illustrated in Fig.(8) for specific parameter values where the energy
conditions are met. Notably, $F_{g}$ is zero because $\Phi$ is assumed to be
constant, meaning the gravitational force has no impact on our model. The
figure shows that the other two forces are equal in magnitude but opposite to
each other. This indicates that the equilibrium of forces is achieved through
the combined effect of the three force terms, supporting the system’s
stability. The same behavior is observed in a model with a specific energy
density.

Figure 8: The three forces for the equilibrium conditions are plotted against
$r$; Left panel: for $\omega_{t}=-1/3,\,\alpha=0.001,\,r_{0}=1,\,q=0.3$, Right
panel: for $\alpha=0.001,\,\beta=4,\,\rho_{0}=0.01,\,q=0.3,\,r_{0}=1$. Here
$F_{g}=0$ for a constant redshift function. Note that the case $q=0$ is also
shown for reference.

##  VI Conclusion

In the present work, we constructed models of static wormholes within the
framework of in 4-dimensional Einstein-Gauss-Bonnet (4D EGB) gravity an
(an)isotropic energy momentum tensor (EMT) and a Maxwell field as supporting
matters for the wormhole geometry assuming a constant redshift function. We
obtained exact spherically symmetric wormhole solutions in 4D EGB gravity for
an isotropic and anisotropic matter sources under the effect of electric
charge. Furthermore, we examined the null, weak and strong energy conditions
at the wormhole throat of radius $r=r_{0}$. We demonstrated that at the
wormhole throat, the classical energy conditions are violated by arbitrary
small amount. Additionally, we analyze the wormhole geometry incorporating
semiclassical corrections through embedding diagrams.

From the viewpoint of the equilibrium condition for wormhole configuration, we
considered the modified Tolman-Oppenheimer-Volkov equation. The equation can
be devided into three parts: the anisotropic force $F_{a}=2(p_{t}-p_{r})/r$,
the hydrostatic force $F_{h}=-dp_{r}/dr$ and the force due to gravitational
$F_{g}=-\big{(}\rho+p_{r}\big{)}\Phi^{\prime}/2$, where the last one vanishes
in our model since the redshift function is assumed to be constant. Taking the
derivative of radial pressure, it is easy to check that $F_{h}=-F_{a}$ and
hence we can deduce that the wormhole solutions are in equilibrium as the
anisotropic and hydrostatic forces cancel each other

Furthermore, the results in this study generalize previous discussions in
Ref.Jusufi:2020yus  by extending them beyond the limit $q\rightarrow 0$. The
results also provide a broader generalization to Morris-Thorne wormholes of
General Relativity, which are encompassed as a special case in the limits
$\alpha\rightarrow 0$ and $q\rightarrow 0$. Exploring the possibility of
extending wormholes to include rotation and applying them to more general
Lovelock gravity theories [93] are current and intriguing research challenges
actively under consideration. The observational aspects of wormhole solutions
are also worthy of study, such as through the utilization of gravitational
lensing effects, see e.g., Ref.Mehdizadeh:2018smu , and shadows and photon
spheres in static and rotating traversable wormholes Tangphati:2023mpk  as
well as a charge one Saleem:2024kld .

###### Acknowledgements.

This work is financially supported by Thailand NSRF via PMU-B under grant
number PCB37G6600138.

## References

  * (1) L. Flamm, “Beitrage zur Einsteinschen Gravitationstheorie,” Phys. Z. 17 (1916) 448
  * (2) A. Einstein and N. Rosen, Phys. Rev. 48 (1935), 73-77
  * (3) M. S. Morris and K. S. Thorne, Am. J. Phys. 56 (1988), 395-412
  * (4) T. Harko, F. S. N. Lobo, M. K. Mak and S. V. Sushkov, Phys. Rev. D 87, no.6, 067504 (2013)
  * (5) M. K. Zangeneh, F. S. N. Lobo and N. Riazi, Phys. Rev. D 90, no.2, 024072 (2014)
  * (6) A. Banerjee, M. K. Jasim and S. G. Ghosh, Annals Phys. 433, 168575 (2021)
  * (7) R. Shaikh and S. Kar, Phys. Rev. D 94, no.2, 024011 (2016)
  * (8) K. A. Bronnikov and A. M. Galiakhmetov, Grav. Cosmol. 21, no.4, 283-288 (2015)
  * (9) V. De Falco and S. Capozziello, Phys. Rev. D 108, no.10, 104030 (2023)
  * (10) J. L. Rosa, J. P. S. Lemos and F. S. N. Lobo, Phys. Rev. D 98, no.6, 064054 (2018)
  * (11) J. L. Rosa and P. M. Kull, Eur. Phys. J. C 82, no.12, 1154 (2022)
  * (12) J. L. Rosa, N. Ganiyeva and F. S. N. Lobo, Eur. Phys. J. C 83, no.11, 1040 (2023)
  * (13) J. L. Rosa, Phys. Rev. D 104, no.6, 064002 (2021)
  * (14) J. L. Rosa, R. André and J. P. S. Lemos, Gen. Rel. Grav. 55, no.5, 65 (2023)
  * (15) V. De Falco, E. Battista, S. Capozziello and M. De Laurentis, Phys. Rev. D 101, no.10, 104037 (2020)
  * (16) V. De Falco, E. Battista, S. Capozziello and M. De Laurentis, Phys. Rev. D 103, no.4, 044007 (2021)
  * (17) V. De Falco, M. De Laurentis and S. Capozziello, Phys. Rev. D 104, no.2, 024053 (2021)
  * (18) V. De Falco, Phys. Rev. D 108, no.2, 024051 (2023)
  * (19) V. De Falco, E. Battista, S. Capozziello and M. De Laurentis, Eur. Phys. J. C 81, no.2, 157 (2021)
  * (20) R. Deshpande and O. Lunin, Nucl. Phys. B 996, 116355 (2023)
  * (21) B. Bhawal and S. Kar, Phys. Rev. D 46, 2464-2468 (1992)
  * (22) M. R. Mehdizadeh, Eur. Phys. J. C 80, no.4, 310 (2020)
  * (23) M. Ilyas and K. Bamba, JCAP 10, 038 (2023)
  * (24) M. Zubair and M. Farooq, Eur. Phys. J. C 83, no.6, 507 (2023)
  * (25) S. Rani and A. Jawad, Adv. High Energy Phys. 2016, 7815242 (2016)
  * (26) M. R. Mehdizadeh, M. Kord Zangeneh and F. S. N. Lobo, Phys. Rev. D 91, no.8, 084004 (2015)
  * (27) G. Dotti, J. Oliva and R. Troncoso, Phys. Rev. D 76, 064038 (2007)
  * (28) H. Maeda and M. Nozawa, Phys. Rev. D 78, 024005 (2008)
  * (29) P. G. S. Fernandes, P. Carrilho, T. Clifton and D. J. Mulryne, Class. Quant. Grav. 39 (2022) no.6, 063001
  * (30) D. Glavan and C. Lin, Phys. Rev. Lett. 124 (2020) no.8, 081301
  * (31) Y. Tomozawa, [arXiv:1107.1424 [gr-qc]].
  * (32) P. G. S. Fernandes, Phys. Lett. B 805 (2020), 135468
  * (33) B. Eslam Panah, K. Jafarzade and S. H. Hendi, Nucl. Phys. B 961 (2020), 115269
  * (34) K. Jusufi, A. Banerjee and S. G. Ghosh, Eur. Phys. J. C 80 (2020) no.8, 698
  * (35) T. Kokubu and T. Harada, Universe 6, no.11, 197 (2020)
  * (36) C. Y. Zhang, C. Niu, W. L. Qian, X. Wang and P. Liu, Chin. J. Phys. 83, 527-538 (2023)
  * (37) N. Godani, D. V. Singh and G. C. Samanta, Phys. Dark Univ. 35, 100952 (2022)
  * (38) Shweta, U. K. Sharma and A. K. Mishra, Int. J. Geom. Meth. Mod. Phys. 20, no.08, 2350140 (2023)
  * (39) Z. Hassan and P. K. Sahoo, Annalen Phys. 536, 2400114 (2024)
  * (40) M. R. Mehdizadeh and A. H. Ziaie, Phys. Rev. D 99, no.6, 064033 (2019)
  * (41) S. W. Kim and H. Lee, Phys. Rev. D 63, 064014 (2001)
  * (42) P. H. R. S. Moraes and P. K. Sahoo, Phys. Rev. D 97, no.2, 024007 (2018)
  * (43) S. Capozziello, T. Harko, T. S. Koivisto, F. S. N. Lobo and G. J. Olmo, Phys. Rev. D 86, 127504 (2012)
  * (44) J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374-381 (1939)
  * (45) R. C. Tolman, Phys. Rev. 55, 364-373 (1939)
  * (46) O. Sokoliuk and A. Baransky, Eur. Phys. J. C 81, no.8, 781 (2021)
  * (47) A. Dutta, D. Roy, N. J. Pullisseri and S. Chakraborty, Eur. Phys. J. C 83, no.6, 500 (2023)
  * (48) D. G. Boulware and S. Deser, Phys. Rev. Lett. 55, 2656 (1985)
  * (49) J. T. Wheeler, Nucl. Phys. B 268, 737-746 (1986)
  * (50) P. G. S. Fernandes, P. Carrilho, T. Clifton and D. J. Mulryne, Phys. Rev. D 102, no.2, 024025 (2020)
  * (51) R. A. Hennigar, D. Kubizňák, R. B. Mann and C. Pollack, JHEP 07, 027 (2020)
  * (52) F. S. N. Lobo, F. Parsaei and N. Riazi, Phys. Rev. D 87, no.8, 084030 (2013)
  * (53) T. Tangphati, P. Channuie, K. Bamba and D. Momeni, [arXiv:2310.16916 [gr-qc]].
  * (54) R. Saleem, M. I. Aslam and S. Shahid, Eur. Phys. J. C 84, no.5, 480 (2024)