# Dynamics of ultrarelativistic charged particles with strong radiation
reaction.
II. Entry into Aristotelian equilibrium

Yangyang Cai  Samuel E. Gralla  Department of Physics, University of Arizona,
1118 E 4th Street, U.S.A  Vasileios Paschalidis  Department of Physics,
University of Arizona, 1118 E 4th Street, U.S.A  Department of Astronomy,
University of Arizona, 933 N Cherry Ave, U.S.A

###### Abstract

As first proposed by Gruzinov, a charged particle moving in strong
electromagnetic fields can enter an equilibrium state where the power input
from the electric field is balanced by radiative losses. When this occurs, the
particle moves at nearly light speed along special directions called the
principal null directions (PNDs) of the electromagnetic field. This
equilibrium is “Aristotelian” in that the particle velocity, rather than
acceleration, is determined by the local electromagnetic field. In paper I of
this series, we analytically derived the complete formula for the particle
velocity at leading order in its deviation from the PND, starting from the
fundamental Landau-Lifshitz (LL) equation governing charged particle motion,
and demonstrated agreement with numerical solutions of the LL equation. We
also identified five necessary conditions on the field configuration for the
equilibrium to occur. In this paper we study the entry into equilibrium using
a similar combination of analytical and numerical techniques. We simplify the
necessary conditions and provide strong numerical evidence that they are also
sufficient for equilibrium to occur. Based on exact and approximate solutions
to the LL equation, we identify key timescales and properties of entry into
equilibrium and show quantitative agreement with numerical simulations. Part
of this analysis shows analytically that the equilibrium is linearly stable
and identifies the presence of oscillations during entry, which may have
distinctive radiative signatures. Our results provide a solid foundation for
using the Aristotelian approximation when modeling relativistic plasmas with
strong electromagnetic fields.

††preprint: APS/123-QED

##  I Introduction

Every electromagnetic field defines, algebraically at each point, a pair of
(possibly identical) light-speed velocities, its principal null directions
(PNDs) [1]. Recently, it has become clear that this mathematical notion has an
elegant physical manifestation: ultrarelativistic charged particles follow the
PNDs. The phenomenon appears to be quite universal in that it emerges in
different regimes and with different mechanisms regulating the ultimate
particle speed, such as classical radiation reaction in magnetically dominated
fields relevant to astrophysics [2, 3, 4, 5, 6, 7] or quantum radiation
reaction in nearly null fields relevant to laser-plasma physics [8, 9, 10, 11,
12]. This regime is Aristotelian in that the particle velocity, rather than
acceleration, is determined by the local electromagnetic field.

In paper I of this series [13] we initiated a detailed study of Aristotelian
motion for classical charged particles in strong external fields. We
considered the fundamental Landau-Lifshitz (LL) equation [14], which includes
both Lorentz force and self-force. We adopted the approximation that a
particle is nearly, but not exactly, moving on a PND, and derived equations
for its velocity at leading order in the deviation from the PND. We identified
precise conditions on the field configuration that are necessary for the
equilibrium to occur. Finally, we demonstrated numerical agreement of this
approximation with full solutions of the LL equation in the appropriate
regime, using a new numerical code.

In this paper we will use similar analytical and numerical techniques to study
the entry into Aristotelian equilibrium. Our main results are: (1) the
necessary conditions identified in paper I are in fact sufficient for
equilibrium to occur; (2) the equilibrium is linearly stable; (3) in some
parameter ranges there are oscillations during the approach to equilibrium,
whose properties we study analytically. Together with the findings of paper I,
these results provide a definite prescription for using the Aristotelian
approximation in practice.

This paper is organized as follows. In Sec. II we review the Aristotelian
equilibrium and relate the assumptions and results of paper I [13] to the
simple versions given originally by Gruzinov [15, 4]. In Sec. III we review
the LL equation and introduce notation. In Sec. IV we show an example of
Aristotelian equilibrium in a helical field configuration. In Sec. V we
analytically study the approach to equilibrium and demonstrate agreement with
numerical simulations. In Sec. VI we perform a large numerical parameter
survey that validates our conditions for entry into equilibrium. Finally, in
Sec. VII we summarize our results, focusing on a simple prescription for using
the Aristotelian approximation in astronomical modeling. Appendix A provides
details of our numerical scheme. We use Gaussian units with the speed of light
set equal to one.

##  II Aristotelian Equilibrium

In this section we review the properties of the PNDs and the Aristotelian
equilibrium. The PNDs [1] $\ell_{+}^{\mu}$ and $\ell_{-}^{\mu}$ of an
electromagnetic field $F_{\mu\nu}$ are the solutions to the pointwise
eigenvalue equation

$\displaystyle F^{\mu}{}_{\nu}\ell_{\pm}^{\nu}=\pm E_{0}\ell_{\pm}^{\nu}.$ (1)

The explicit solution in terms of electric and magnetic fields is

$\displaystyle\ell^{\mu}_{\pm}$ $\displaystyle=(1,\vec{v}_{\pm}),$ (2)
$\displaystyle\vec{v}_{\pm}$
$\displaystyle=\frac{\vec{E}\times\vec{B}\pm(B_{0}\vec{B}+E_{0}\vec{E})}{B^{2}+E_{0}^{2}},$
(3)

where $E_{0}$ and $B_{0}$ are given in terms of the invariants
$P=\vec{B}^{2}-\vec{E}^{2}$ and $Q=\vec{E}\cdot\vec{B}$ as

$\displaystyle E_{0}$ $\displaystyle=\sqrt{\sqrt{(P/2)^{2}+Q^{2}}-P/2}$ (4)
$\displaystyle B_{0}$
$\displaystyle=\textrm{sign}(Q)\sqrt{\sqrt{(P/2)^{2}+Q^{2}}+P/2}.$ (5)

When the PNDs are degenerate ($\ell^{\mu}_{+}=\ell^{\mu}_{-}$), the eigenvalue
vanishes and the field is null ($E_{0}=B_{0}=0$), in which case the electric
and magnetic fields are orthogonal and equal in magnitude in any Lorentz
frame. When the PNDs are distinct, the eigenvalues are $\pm E_{0}$, and
$E_{0}>0$ is the magnitude of the electric field in any frame where the
electric and magnetic fields are parallel. Similarly, $|B_{0}|$ is the
magnitude of the magnetic field in such a frame, with $B_{0}$
positive/negative when the fields are aligned/antialigned.

The PNDs define integral curves $\vec{x}_{\pm}(t)$ by the equation
$d\vec{x}_{\pm}/dt=v_{\pm}$. The parameter $t$ is the arc length of the space
curve, since $v_{\pm}^{2}=1$. On each curve we may erect a Frenet Serret frame
$\\{\vec{\ell},\vec{n},\vec{k}\\}$ with $\vec{\ell}=v_{\pm}$ [13]. Since these
curves fill space, for each choice of $\pm$ we have a full orthonormal basis
for vector fields. In particular, we may decompose the velocity vector of a
charged particle as

$\displaystyle\vec{v}=v_{\ell}\vec{\ell}+v_{n}\vec{n}+v_{k}\vec{k},$ (6)

where we choose $\vec{\ell}=\vec{v}_{+}$ for positively charged particles and
$\vec{\ell}=\vec{v}_{-}$ for negatively charged particles. As this is an
orthonormal frame, the Lorentz factor is reconstructed by

$\displaystyle\gamma=\frac{1}{\sqrt{1-v_{\ell}^{2}-v_{n}^{2}-v_{k}^{2}}}.$ (7)

The Frenet-Serret vectors $\\{\vec{\ell},\vec{n},\vec{k}\\}$ and their
associated curvature $\kappa$ and torsion $\iota$ are local functions of the
electric and magnetic fields for each choice of $\pm$. This may be seen from
the defining equations,

$\displaystyle\vec{\ell}$ $\displaystyle=\vec{v}_{\pm},$ (8)
$\displaystyle\kappa\vec{n}$
$\displaystyle=(\vec{\ell}\cdot\vec{\nabla})\vec{\ell},\qquad(\vec{n}\cdot\vec{n}=1)$
(9) $\displaystyle\vec{k}$ $\displaystyle=\vec{\ell}\times\vec{n},$ (10)
$\displaystyle\iota\vec{n}$
$\displaystyle=-(\vec{\ell}\cdot\vec{\nabla})\vec{k}.$ (11)

The first vector $\vec{\ell}=v_{\pm}$ is determined by the values of $\vec{E}$
and $\vec{B}$ via Eqs. (8) and (3). The second vector $\hat{n}$ and the
curvature $\kappa$ involve first derivatives as well [Eq. (9)]. The third
vector $\vec{k}=\vec{\ell}\times\vec{n}$ also depends on first derivatives.
Finally, the torsion $\iota$ depends on first and second derivatives [Eq.
(10)]. We will also define the radius of curvature $R$,

$\displaystyle R=\frac{1}{\kappa}.$ (12)

Note that while $E_{0}$, $B_{0}$ and $\ell^{\mu}$ are invariant notions, the
projection of the null vector $\ell^{\mu}$ to the spatial velocity
$v_{\pm}=\vec{\ell}$ depends on the choice of Lorentz frame. The corresponding
Frenet-Serret basis, together with its associated curvature and torsion, are
similarly non-invariant. We will be formulating assumptions and deriving
results in terms of these non-invariant quantities; the interpretation is that
our results hold in frames satisfying our assumptions. Note that we will use
the phrase “PND” for both the invariant null direction in a spacetime sense
and the non-invariant spatial direction $\vec{v}_{\pm}$ in a given frame.
Context will make clear which notion is meant.

A particle of charge $q$ and mass $m$ defines a length scale $\mathcal{R}$ and
a field scale $\mathcal{E}$ by

$\displaystyle\mathcal{R}\equiv\frac{q^{2}}{m},\qquad\mathcal{E}\equiv\frac{3}{2}\frac{m^{2}}{|q|^{3}},$
(13)

with a conventional factor of $3/2$. The “classical electron radius”
$\mathcal{R}$ is the distance where the electrostatic self-energy of a point
charge equals its rest mass, and the “classical critical field” $\mathcal{E}$
is (three-halves times) the strength of the electric field at this location.
These represent typical scales at which the classical description of the
particle will break down. These quantities also define time and magnetic field
scales after multiplying by suitable factors of the speed of light (set here
to unity). We therefore assume

$\displaystyle E_{0},B_{0}\ll\mathcal{E}\qquad L,T\gg\mathcal{R},$ (14)

where $L$ and $T$ are typical length and time scales for the field
configuration to change significantly.

###  II.1 Original Derivation

We now summarize Gruzinov’s original arguments for Aristotelian motion [4, 15]
using the notation of our paper. Gruzinov considered the case where a particle
moves primarily along a PND, i.e.,

$\displaystyle\gamma$ $\displaystyle\gg 1,$ (15)
$\displaystyle\sqrt{v_{n}^{2}+v_{k}^{2}}$ $\displaystyle\ll 1.$ (16)

If we approximate the particle motion as a circle of radius equal to the
radius of curvature $R=1/\kappa$ of the PND, the power radiated is
$(2/3)q^{2}R^{2}\gamma^{4}$. Gruzinov assumed that this “curvature radiation”
power is balanced by the Lorentz force power $|q|E_{0}$,

$\displaystyle\frac{2}{3}q^{2}R^{2}\gamma^{4}=|q|E_{0}.$ (17)

Solving for $\gamma$ gives the Gruzinov formula for the Lorentz factor,

$\displaystyle\gamma_{g}$
$\displaystyle=\left(\frac{3E_{0}R^{2}}{2|q|}\right)^{1/4}.$ (18)

Gruzinov obtained an additional condition for this equilibrium based on the
idea that it can only occur when the particle has enough space to be
accelerated to this terminal Lorentz factor before the PND curves
significantly. The curvature radius $R$ sets the scale over which the uniform
field approximation breaks down, so the particle must be able to gain
$m\gamma_{g}$ energy in a region much smaller than $R$. The typical energy
gain over a region of size $D$ is $qE_{0}D$, so we obtain the condition

$\displaystyle m\gamma_{g}\ll qE_{0}R,$ (19)

which may equivalently be written

$\displaystyle\gamma_{g}\ll\frac{E_{0}R}{\mathcal{E}\mathcal{R}}$ (20)

or

$\displaystyle\mathcal{E}^{3/2}\mathcal{R}\ll E_{0}^{3/2}R.$ (21)

###  II.2 LL Derivation

In paper I [13] we sought to better understand the Aristotelian equilibrium by
studying the fundamental LL equation of charged particle dynamics. We again
assumed motion along a PND [Eqs. (15) and (16)]. However, the energy balance
condition (17) is inappropriate in this context since (1) it requires the
particle motion to be treated as circular, an uncontrolled approximation whose
compatibility with the LL equation is not obvious; and (2) the local LL
dynamics contains more information than just energy conservation. Instead, to
express the idea of energy balance, we assumed that the local change in energy
is small compared to the typical value set by the Lorentz force,

$\displaystyle m\left|\frac{d\gamma}{dt}\right|$ $\displaystyle\ll|q|E_{0}.$
(22)

We also assumed that the timescale $T$ for changes in the field is long
compared to the lengthscale $L$ for spatial changes in the field,

$\displaystyle T\gg L.$ (23)

We found that Gruzinov’s equilibrium emerged only after a final additional
assumption,

$\displaystyle|\iota|$
$\displaystyle\ll\frac{|q|}{m\gamma}\textrm{Max}\\{E_{0},|B_{0}|\\},$ (24)

where $\iota$ is the torsion of the PND. These three conditions (22), (23) and
(24) replace the assumption (17) of global power balance in approximate
circular motion. Eq. (22) guarantees approximate power balance locally at the
level of the LL equation (as opposed to globally, at the level of total
radiated energy), while Eqs. (23) and (22) reflect the approximate circular
motion with radius $R$.

Under these five assumptions (15), (16), (22), (23), (24), we found that the
field determines the velocities pointwise as

$\displaystyle\gamma$
$\displaystyle=\left(\frac{9}{4}\frac{R^{2}}{\mathcal{R}^{2}}\frac{E_{0}}{\mathcal{E}}\right)^{1/4}=\gamma_{g}$
(25) $\displaystyle v_{n}$
$\displaystyle=-\frac{1+\delta}{\gamma}\sqrt{\frac{E_{0}}{\mathcal{E}}}$ (26)
$\displaystyle v_{k}$
$\displaystyle=\frac{\delta}{\gamma}\frac{B_{0}}{E_{0}}\sqrt{\frac{E_{0}}{\mathcal{E}}},$
(27)

where we used the definitions (13) and also introduced

$\displaystyle\delta=\frac{E_{0}\mathcal{E}}{E_{0}^{2}+B_{0}^{2}}.$ (28)

We thus reproduced Gruzinov’s Lorentz factor (18) and derived new expressions
(26) and (27) for the drift velocities.

###  II.3 Conditions on the field configuration

The five assumptions (15), (16), and (22)–(24) can be recast as conditions
purely on the field configuration by using the results (25)–(27). Presenting
the results as five conditions $C_{i}\ll 1$, the conditions are

$\displaystyle C_{1}$
$\displaystyle=\frac{\mathcal{R}}{R}\sqrt{\frac{\mathcal{E}}{E_{0}}}\ll 1$
(29) $\displaystyle C_{2}$
$\displaystyle=\delta\frac{\mathcal{R}}{R}\sqrt{\frac{\mathcal{E}}{E_{0}}}\ll
1$ (30) $\displaystyle C_{3}$
$\displaystyle=|\iota|\sqrt{R\mathcal{R}}\left(\frac{E_{0}}{\mathcal{E}}\right)^{1/4}\frac{\mathcal{E}}{\textrm{Max}\\{E_{0},|B_{0}|\\}}\ll
1$ (31) $\displaystyle C_{4}$
$\displaystyle=\eta\sqrt{\frac{\mathcal{R}}{R}}\left(\frac{\mathcal{E}}{E_{0}}\right)^{3/4}\ll
1,$ (32) $\displaystyle C_{5}$ $\displaystyle=\frac{L}{T}\ll 1$ (33)

with

$\displaystyle\eta=\left|\vec{v}\cdot\vec{\nabla}R+\frac{R}{2E_{0}}\vec{v}\cdot\vec{\nabla}E_{0}\right|.$
(34)

In this last expression, it is understood that Eqs. (25)–(27) [as well as (7)]
are to be used to express $\vec{v}$ in terms of the local field configuration.
In obtaining (30) we have used the fact that
$\sqrt{v_{n}^{2}+v_{k}^{2}}\approx\sqrt{\delta}/\gamma$ under the assumptions
(14).

The names $C_{i}$ are derived from paper I. In this paper, $C_{1}$–$C_{5}$
correspond (respectively) to Eqs. (15), (16), (24), (22), and (23). In
particular, $C_{1}\ll 1$ and $C_{2}\ll 1$ express the approximate PND motion
(15) and (16), while $C_{4}\ll 1$ expresses the approximate local equilibrium
(22).

Our expression for $C_{4}$ differs from that presented in paper I, where we
assumed that $\eta\sim R/L$. Generically we will have $\eta\sim 1$. In this
paper we also include cases where $\eta$ is very small, since they arise quite
naturally in our numerical studies, and help to illustrate the differences
between our approach and that originally given by Gruzinov (Sec. II.1).

Notice that $C_{1}$ and $C_{2}$ are related to $C_{4}$ by the equations

$\displaystyle\eta C_{1}^{3/2}$ $\displaystyle=\frac{\mathcal{R}}{R}C_{4}\ll
C_{4}$ (35) $\displaystyle\eta C_{2}$
$\displaystyle\leq\sqrt{\frac{\mathcal{R}}{R}}C_{4}\ll C_{4}.$ (36)

In the second line we used the bound $\delta\leq\mathcal{E}/E_{0}$, which
follows from the definition (28) of $\delta$. In the generic case that
$\eta\sim 1$, we see that both $C_{1}$ and $C_{2}$ are small compared with
$C_{4}$. This means that $C_{4}\ll 1$ in fact implies that $C_{1}\ll 1$ and
$C_{2}\ll 1$, and generically one can ignore the $C_{1}$ and $C_{2}$
conditions. While $\eta$ can be a small number in special cases (such as the
circular fields studied in Sec. VI.1 below), it would have to be less than
$\sim\sqrt{\mathcal{R}/R}$ for the $C_{1}$ and $C_{2}$ conditions to become
relevant. For macroscopic fields this factor is very small; for example,
$\sqrt{\mathcal{R}/R}\approx 5\times 10^{-8}$ if $R$ is one meter. A field
configuration with $\eta$ this small would be extraordinarily fine-tuned. We
therefore conclude that for the macroscopic fields of relevance in
astrophysics, we may always ignore the $C_{1}$ and $C_{2}$ conditions, since
they are implied by the $C_{4}$ condition.

Let us therefore focus on the $C_{4}$ condition (32). This condition can be
rewritten as

$\displaystyle\eta^{2}\mathcal{R}\mathcal{E}^{3/2}\ll RE_{0}^{3/2},$ (37)

showing that it is in fact equivalent to Gruzinov’s condition (21) in the
generic case $\eta\sim 1$. This agreement is interesting since Gruzinov’s
condition (21) arose from reasoning about when a particle can enter
equilibrium, whereas our condition (37) arose from the assumption (22) that
equilibrium had been achieved. We will see numerically in Sec. VI that Eq.
(37) is necessary and sufficient for equilibrium, provided the other
conditions are satisfied.

Finally we discuss $C_{3}$, which satisfies

$\displaystyle C_{3}$
$\displaystyle=\frac{|\iota|}{\kappa}\eta^{-1}C_{4}\frac{E_{0}}{\textrm{Max}\\{E_{0},|B_{0}|\\}}.$
(38)

The last factor $E_{0}/\textrm{Max}\\{E_{0},|B_{0}|\\}$ is always $\lesssim 1$
and will be small for pulsars, where $E_{0}\ll|B_{0}|$. Again assuming the
generic case $\eta\sim 1$, Eq. (38) shows that large values of the torsion-to-
curvature ratio $|\iota|/\kappa$ are required for the condition $C_{3}\ll 1$
to be violated in a regime where the condition $C_{4}\ll 1$ is satisfied. For
fields with $|\iota|/\kappa\lesssim 1$, the $C_{3}$ condition can be ignored,
since it is implied by the $C_{4}$ condition. We discuss cases with large
$|\iota|/\kappa$ in Sec. IV below.

The conditions $C_{i}\ll 1$ arose in paper I as a minimal set of assumptions
under which the LL equation reduced to the simple results (25)–(27). Although
we do not attempt any rigorous mathematical proof, it seems clear from the
steps of the derivation there that these assumptions were all necessary for
the simple results to emerge. We therefore regard the conditions (29)–(33) as
necessary conditions for the equilibrium described by (25)–(27). One of the
main goals of this paper is to argue that they are also sufficient for
equilibrium to occur.

##  III Landau-Lifshitz Equation

The LL equation is reviewed in paper I [13]; here we briefly recap the physics
and introduce our notation for numerical simulations. In terms of $\vec{E}$
and $\vec{B}$ fields, the LL equation is111We have dropped terms involving the
derivative of the field strength, which are always negligible in the regime of
validity discussed below (39).

$\displaystyle\frac{d(\gamma\vec{v})}{d\tau}=\frac{\gamma q}{m}\Bigg{\\{}$
$\displaystyle\vec{E}+\vec{v}\times\vec{B}$
$\displaystyle\pm\frac{1}{\mathcal{E}}\left[(\vec{E}\cdot\vec{v})\vec{E}+(\vec{E}+\vec{v}\times\vec{B})\times\vec{B}\right]$
$\displaystyle\mp\frac{\gamma^{2}}{\mathcal{E}}\left[(\vec{E}+\vec{v}\times\vec{B})^{2}-(\vec{E}\cdot\vec{v})^{2}\right]\Bigg{\\}},$
(39)

where $\pm$ is the sign of $q$, $\vec{v}$ is the velocity,
$\gamma=1/\sqrt{1-v^{2}}$ is the Lorentz factor, and $\tau$ is the proper
time. As discussed originally by LL [14], this equation is valid provided that
the corrections to Lorentz force motion (the final two terms) are small
compared to the Lorentz force in the rest frame of the particle. However, in
the lab frame they may be of comparable or greater magnitude, and this is the
regime of interest to us.

For numerical purposes, we normalize the equation using the natural scales of
the equation, together with a dimensionless number $\chi$ that can be chosen
for convenience. Defining

$\displaystyle\tilde{\tau}$
$\displaystyle=\frac{3}{2}\chi^{2}\frac{\tau}{\mathcal{R}}$ (40)
$\displaystyle\tilde{\vec{E}}$ $\displaystyle=\frac{\vec{E}}{\chi\mathcal{E}}$
(41) $\displaystyle\tilde{\vec{B}}$
$\displaystyle=\frac{\vec{B}}{\chi\mathcal{E}}$ (42)
$\displaystyle\tilde{\vec{p}}$
$\displaystyle=\frac{\vec{p}}{m}=\gamma\vec{v},$ (43)

Eq. (39) becomes

$\displaystyle\frac{d\tilde{\vec{p}}}{d\tilde{\tau}}$
$\displaystyle=\frac{1}{\chi}\vec{f_{L}}\pm\left[(\tilde{\vec{E}}\cdot\tilde{\vec{p}})\tilde{\vec{E}}+\vec{f_{L}}\times\tilde{\vec{B}}\right]\mp\left[f_{L}^{2}-(\tilde{\vec{E}}\cdot\tilde{\vec{p}})^{2}\right]\tilde{\vec{p}},$
(44)

where $\vec{f_{L}}=\gamma\tilde{\vec{E}}+\tilde{\vec{p}}\times\tilde{\vec{B}}$
is the Lorentz force term. The Lorentz factor is related to the rescaled
momentum by

$\displaystyle\gamma$ $\displaystyle=\sqrt{1+\tilde{p}^{2}}.$ (45)

Once the equation is solved for $\tilde{\vec{p}}(\tilde{\tau})$, the position
can be recovered by a subsequent integration. For convenience we define a
normalized position vector $\tilde{\vec{x}}$,

$\displaystyle\tilde{\vec{x}}$
$\displaystyle=\frac{3}{2}\chi^{2}\frac{\vec{x}}{\mathcal{R}}=\frac{3}{2}\chi^{2}\frac{1}{\mathcal{R}}\int\vec{v}dt=\int\tilde{\vec{p}}\
d\tilde{\tau}.$ (46)

If necessary, one can also determine the lab-frame time $t$ from $dt=\gamma
d\tau$.

We will choose the value of $\chi$ so that lengths are measured in meters when
the electron is considered,

$\displaystyle\chi=\sqrt{\frac{2}{3}\frac{\mathcal{R}}{1{\rm m}}}\approx
4.33\times 10^{-8}.$ (47)

In particular for electrons we have

$\displaystyle\tilde{\vec{x}}$ $\displaystyle=\textrm{(position in meters)}$
(48) $\displaystyle\tilde{\tau}$ $\displaystyle=\textrm{(proper time in
meters)}$ (49) $\displaystyle\tilde{\vec{E}}$ $\displaystyle=\textrm{(electric
field in units of $10^{13}$ V/m)}$ (50) $\displaystyle\tilde{\vec{B}}$
$\displaystyle=\textrm{(magnetic field in units of $10^{8}$ G)},$ (51)

These are somewhat convenient for pulsars, which have radii of $\sim 10$km and
magnetic fields of $10^{8}$–$10^{15}$G.

Our numerical method is described in Appendix. A.

##  IV An example with torsion

In paper I we presented a numerical example of Aristotelian equilibrium in
which the torsion was precisely zero. Here we complement that case with an
example where the torsion is significant. We consider a simple field
configuration consisting of parallel electric and magnetic fields tangent to
helical curves that fill space. In terms of some constant $h>0$, the field
configuration is

$\displaystyle\vec{E}$
$\displaystyle=\frac{E_{0}}{\sqrt{h^{2}+x^{2}+y^{2}}}\left\\{-y,x,h\right\\}$
(52) $\displaystyle\vec{B}$
$\displaystyle=\frac{B_{0}}{\sqrt{h^{2}+x^{2}+y^{2}}}\left\\{-y,x,h\right\\}.$
(53)

The PNDs are also helical, and their curvature and torsion are

$\displaystyle\kappa$
$\displaystyle=\frac{\sqrt{x^{2}+y^{2}}}{x^{2}+y^{2}+h^{2}}$ (54)
$\displaystyle\iota$ $\displaystyle=\frac{h}{x^{2}+y^{2}+h^{2}},$ (55)

with ratio

$\displaystyle\frac{\iota}{\kappa}=\frac{h}{\sqrt{x^{2}+y^{2}}}.$ (56)

We choose field strengths $\tilde{E}_{0}=\tilde{B}_{0}=1$ and height
$\tilde{h}=10$. For electrons, the physical units of $\tilde{E}_{0}$ and
$\tilde{B}_{0}$ are given in Eqs. (50) and (51), while the value of
$\tilde{h}$ represents meters. We start the particle at $(1,0,0)$ and it
reaches equilibrium within several timescales $m/(qE)$. Fig. 1 shows a portion
of the ensuing trajectory, which closely hugs a PND while slowly drifting
outwards. Although the torsion of this PND is times larger than the curvature,
we still have $C_{3}\approx.02\ll 1$ and the formulas for the Aristotelian
equilibrium should still be valid. Fig. 2 shows indeed that the numerical
results match the predicted analytical formulas.

By suitable choice of the parameters $E_{0}$, $B_{0}$, and $h$, one can
arrange for $C_{3}$ to be arbitrarily large over an arbitrarily large region
of space. However, since $\kappa$ (54) falls off more slowly than $\iota$
(55), the value of $C_{3}$ always becomes small at large distances from the
$z$ axis [see Eq. (38)]. In numerical experiments we find that equilibrium
does not occur in the region of large $C_{3}$.222Although $C_{3}\ll 1$ is
required for the formulas (25)–(27) to apply, in principle the particle could
reach an equilibrium described by different formulas. In paper I we in fact
derived a more general formula for the Lorentz factor as a quartic equation
for $\gamma$ (see Eq. (61) of paper I) that does not rely on $C_{3}\ll 1$.
However, this formula does require the equilibrium assumption (22), and we
found in this case that the particle does not enter equilibrium in this sense
while $C_{3}$ is still large. We have not identified a field configuration
where the alternative Lorentz factor formula both applies and is significantly
different from the Gruzinov value. However, the particle still moves primarily
along the principal null direction, now with Lorentz factor growing in time,
indicating yet another regime of PND motion worthy of future exploration. In
this particular numerical experiment, the motion also has a slower outward
drift that eventually takes the particle to a region of small $C_{3}$ and
$C_{4}$, where its Lorentz factor finally settles down to the Gruzinov value
(18). These properties are illustrated in Fig. 3

These simple examples demonstrate that the equilibrium works as expected when
torsion is non-zero, as long as $C_{3}\ll 1$. In the remainder of the paper we
will set the torsion precisely to zero in order to simplify the discussion.

Figure 1: Trajectory of a particle under the helical field configuration (52)
and (53) with equal electric and magnetic field strengths. The particle’s
trajectory is almost tangent to the PND shown in the plot while drifting
outwards slowly. Notice that the $\tilde{z}$ axis has a compressed scale; the
torsion-to-curvature ratio is $\iota/\kappa\approx 10$.
Figure 2: Agreement of the analytical predictions (25)–(27) with numerical
simulation in a case with significant torsion. The trajectory is shown in Fig.
1. The fractional differences $\epsilon_{i}$ are defined as the difference
between the numerical and analytical value, divided by the analytical value,
for $i=\\{\gamma,v_{n},v_{k}\\}$. The bottom two plots show the torsion and
curvature of the PND at the particle position, normalized as
$\tilde{\kappa}=(2/3)\chi^{-2}\mathcal{R}\kappa$ and
$\tilde{\iota}=(2/3)\chi^{-2}\mathcal{R}\kappa$ according to the conventions
of Sec. III. The $x$-axis uses a timescale $\tau_{E}=m/(qE_{0})$.  Figure 3:
An example of evolution through a region where the torsion condition $C_{3}\ll
1$ is violated. Since $C_{3}\propto\iota/\sqrt{\kappa}$ diverges on the
symmetry axis for the nested helical configuration we consider, the torsion
condition is violated near the axis. Here we show an evolution beginning in
this region; the particle starts at $\tilde{\vec{x}}=\\{0.01,0,0\\}$, where
$C_{3}\approx 5.6$. The initial velocity is along the PND with initial Lorentz
factor equal to half the Gruzinov value. The particle still approximately
follows the PND, but with time-variable Lorentz factor, such that it does not
reach equilibrium where $C_{3}$ is large. However, it eventually drifts to a
region of small $C_{3}$ and settles down to the Aristotelian equilibrium.

##  V Approach to Equilibrium

We now study the process by which a particle enters equilibrium, using a
combination of analytical and numerical approaches.

###  V.1 Uniform field solution

At least for some short time, the motion of a particle can be determined in
the approximation that the field is uniform. In this approximation, we can
always work in a frame where $\vec{E}$ and $\vec{B}$ are parallel. The full
analytic solution to the LL equation was found in this case by [16]. We denote
the initial velocity by $(v_{1},v_{2},v_{3})$ with initial Lorentz factor
$\gamma_{0}$, such that the initial four-velocity is given by

$\displaystyle u^{\alpha}(0)=\gamma_{0}(1,v_{1},v_{2},v_{3}).$ (57)

Taking the field to be in the $z$ direction, the analytic solution to the LL
equation is

$\displaystyle\gamma(\tau)$
$\displaystyle=\frac{1}{2}\frac{(1+v_{3})+(1-v_{3})e^{-2\tau/\tau_{E}}}{\sqrt{1-v_{3}^{2}-(v_{1}^{2}+v_{2}^{2})e^{-2\tau/(\delta\tau_{E})}}}e^{\tau/\tau_{E}}$
(58) $\displaystyle v_{x}(\tau)$
$\displaystyle=A(\tau)e^{-\tau/\tau_{\perp}}\sin(\tau/\tau_{B}+\phi)$ (59)
$\displaystyle v_{y}(\tau)$
$\displaystyle=-\text{sign}(qB_{0})A(\tau)e^{-\tau/\tau_{\perp}}\cos(\tau/\tau_{B}+\phi)$
(60) $\displaystyle v_{z}(\tau)$
$\displaystyle=\text{sign}(q)\frac{(1+v_{3})-(1-v_{3})e^{-2\tau/\tau_{E}}}{(1+v_{3})+(1-v_{3})e^{-2\tau/\tau_{E}}}.$
(61)

We are using Cartesian coordinates with $v_{x}(\tau=0)=v_{1}$,
$v_{y}(\tau=0)=v_{2}$, and $v_{z}(\tau=0)=v_{3}$, and $\phi$ is the initial
direction of motion in the $xy$ plane ($\tan\phi=v_{2}/v_{1}$). Here $A(\tau)$
is defined by

$\displaystyle A(\tau)$
$\displaystyle=\frac{2\sqrt{v_{1}^{2}+v_{2}^{2}}}{(1+v_{3})+(1-v_{3})e^{-2\tau/\tau_{E}}},$
(62)

and we introduced three time scales

$\displaystyle\tau_{B}$ $\displaystyle=\frac{m}{|qB_{0}|},$ (63)
$\displaystyle\tau_{E}$ $\displaystyle=\frac{m}{|q|E_{0}}$ (64)
$\displaystyle\tau_{\perp}$
$\displaystyle=\frac{m}{|q|E_{0}}\frac{\delta}{\delta+1}.$ (65)

The definition of $\delta$ was given above in Eq. (28).

We see that the particle asymptotically approaches light speed in the $z$
direction, i.e., it eventually moves along the PND. We interpret this to mean
that radiation damps only the perpendicular momentum. The particle executes a
damped circular motion in the $xy$ plane (around the field direction), in
accordance with the low-$E_{0}$ intuition of synchrotron motion and associated
radiation damping. The period of oscillation is just the classical synchrotron
period $\tau_{B}$ (expressed here in proper time), while the decay time
$\tau_{\perp}$ generalizes the classical synchrotron damping result. In
particular, we have

$\displaystyle\tau_{\perp}\approx\begin{cases}\tau_{B}\frac{\mathcal{E}}{|B_{0}|}&|B_{0}|\gg
E_{0}\\\ \tau_{E}&E_{0}\gg|B_{0}|\end{cases},$ (66)

taking into account $E_{0}\ll\mathcal{E}$. The former case agrees in order of
magnitude with Eq. (A2) of Ref. [17].

To understand the dynamics it is convenient to expand the Lorentz factor at
early and late times,

$\displaystyle\gamma(\tau)\approx\begin{cases}\gamma_{0}\left(1-\tau/\tau_{\rm
drop}-\frac{v_{3}}{\gamma_{0}^{2}}\tau/\tau_{E}\right),&\tau\to 0\\\
\frac{1}{2}\sqrt{\frac{1+v^{3}}{1-v_{3}^{2}}}e^{\tau/\tau_{E}},&\tau\to\infty.\end{cases}$
(67)

where we introduce yet another timescale

$\displaystyle\tau_{\rm
drop}\equiv\frac{\delta}{\gamma_{0}^{2}(v_{1}^{2}+v_{2}^{2})}\tau_{E}.$ (68)

For sufficiently large initial velocity perpendicular to the field,
$\gamma_{0}(v_{1}^{2}+v_{2}^{2})\gg\delta$, we have $\tau_{\rm
drop}\ll\tau_{\rm E}$. In this case there is a sudden drop in Lorentz factor
on timescale $\tau_{\rm drop}$ as the particle loses its perpendicular
momentum, followed by a subsequent rise on timescale $\tau_{E}$ as the
particle gains parallel momentum. The minimal Lorentz factor is determined by
the details of (58); in the special case $\delta\gg 1$, the minimum occurs at
approximately $\sqrt{\delta}$.

In a realistic field configuration, the exponential rise in Lorentz factor is
ultimately cut off by the effects of non-uniform fields. If the conditions are
right for Aristotelian equilibrium, there will be a transition around the
equilibrium Lorentz factor. In order to understand this transition, we now
study the LL equation near the Aristotelian equilibrium solution.

###  V.2 Near equilibrium solution

To study the dynamics near equilibrium, we adopt the same assumptions as of
paper I, assuming time derivatives are perturbatively small (instead of
dropping them entirely). Following the steps of Sec. IVB of that reference,
except retaining time derivatives,333We also drop all terms involving $\iota$,
a step that was justified in a later section of paper I. When time derivatives
are dropped, Eqs. (69), (70), and (71) reduce, respectively, to Eqs. (43),
(54) and (55) of paper I with $\iota=0$. we arrive at three coupled equations
for $\gamma$, $v_{n}$ and $v_{k}$

$\displaystyle\frac{d\gamma}{d\tau}$ $\displaystyle=\frac{|q|\gamma
E_{0}}{m}-\frac{|q|\gamma^{3}}{m\mathcal{E}}(E_{0}^{2}+B_{0}^{2})(v_{n}^{2}+v_{k}^{2})$
(69) $\displaystyle\frac{dv_{n}}{d\tau}$
$\displaystyle=\frac{|q|B_{0}}{m}v_{k}-\frac{|q|E_{0}}{m}\frac{1+\delta}{\delta}v_{n}-\gamma\kappa$
(70) $\displaystyle\frac{dv_{k}}{d\tau}$
$\displaystyle=-\frac{|q|B_{0}}{m}v_{n}-\frac{|q|E_{0}}{m}\frac{1+\delta}{\delta}v_{k}.$
(71)

We now express these quantities as small perturbations of their equilibrium
values,

$\displaystyle\gamma$ $\displaystyle=\Gamma+\bar{\gamma}$ (72) $\displaystyle
v_{n}$ $\displaystyle=V_{n}+\bar{v}_{n}$ (73) $\displaystyle v_{k}$
$\displaystyle=V_{k}+\bar{v}_{k},$ (74)

where $\Gamma,V_{n},V_{k}$ are taken to be Eqs. (25)–(27), respectively.
Regarding $E_{0},B_{0},\kappa$ as constants and linearizing in the barred
quantities, we find

$\displaystyle\tau_{E}\kappa\frac{d\bar{\gamma}}{d\tau}$
$\displaystyle=-2\kappa\bar{\gamma}+2\bar{v}_{n}/\tau_{\perp}-\epsilon
2\bar{v}_{k}/\tau_{B}$ (75) $\displaystyle\frac{d\bar{v}_{n}}{d\tau}$
$\displaystyle=-\bar{v}_{n}/\tau_{\perp}+\epsilon\bar{v}_{k}/\tau_{B}-\bar{\gamma}\kappa$
(76) $\displaystyle\frac{d\bar{v}_{k}}{d\tau}$
$\displaystyle=-\epsilon\bar{v}_{n}/\tau_{B}-\bar{v}_{k}/\tau_{\perp},$ (77)

where the timescales $\tau_{i}$ were introduced in Eqs. (63)–(65) and we have
also written $\epsilon=\textrm{sign}(B_{0})$. Defining

$\displaystyle T$ $\displaystyle=\frac{\tau}{\tau_{E}}$ (78) $\displaystyle b$
$\displaystyle=\epsilon\frac{\tau_{E}}{\tau_{B}}=\frac{B_{0}}{E_{0}}$ (79)
$\displaystyle c$
$\displaystyle=\frac{\tau_{E}}{\tau_{\perp}}=\frac{1+\delta}{\delta}>1,$ (80)

this set of equations may also be written

$\displaystyle\frac{d\bar{X}}{dT}=M\bar{X},$ (81)

with

$\displaystyle\bar{X}=\begin{pmatrix}\bar{\gamma}\kappa\tau_{E}\\\
\bar{v}_{n}\\\ \bar{v}_{k}\end{pmatrix},\quad M=\begin{pmatrix}-2&2c&-2b\\\
-1&-c&b\\\ 0&-b&-c\\\ \end{pmatrix}.$ (82)

Making the ansatz

$\displaystyle\bar{X}=\bar{X}_{0}e^{\lambda T},$ (83)

presents the eigenvalue problem $M\bar{X}_{0}=\lambda\bar{X}_{0}$ (here
$\bar{X}_{0}$ is a constant independent of $T$). If $\lambda^{(i)}$ are the
eigenvalues and $\bar{X}^{(i)}_{0}$ the associated eigenvectors (for
$i=1,2,3$), then the general solution is

$\displaystyle\bar{X}(T)=\sum_{i}C_{i}\textrm{Re}[\bar{X}_{0}^{(i)}e^{\lambda^{(i)}T}],$
(84)

for three real constants $C_{i}$. The eigenvalues and eigenvectors depend only
on $E_{0}$ and $B_{0}$ and can be found numerically for any given values of
$E_{0}$ and $B_{0}$.

We now show that the equilibrium is linearly stable, i.e.,
$\textrm{Re}[\lambda^{(i)}]<0$. The eigenvalues $\lambda^{(i)}$ are the roots
of the characteristic polynomial (defined with a minus sign for convenience),

$\displaystyle f(\lambda)=-\textrm{det}(M-\lambda I)$
$\displaystyle=\lambda^{3}+2(c+1)\lambda^{2}+(b^{2}+6c+c^{2})\lambda+4(b^{2}+c^{2}).$
(85)

Notice that all the coefficients are real and positive. Thus $f$ is strictly
positive for $\lambda\geq 0$, implying that any real root is strictly
negative. We have therefore proven stability in the case of real roots only.
For complex roots, note that such roots must appear in a complex conjugate
pair, since we consider a real polynomial. We may therefore write

$\displaystyle
f(\lambda)=(\lambda-\lambda_{r})(\lambda-\alpha-i\beta)(\lambda-\alpha+i\beta),$
(86)

with $\lambda_{r},\alpha,\beta$ all real. Comparing the coefficient of
$\lambda^{2}$ with the polynomial (V.2), we find that

$\displaystyle 2\alpha=-\lambda_{r}-2(c+1).$ (87)

However, we find that $f$ is negative at $\lambda=-2(c+1)$:

$\displaystyle f(-2(c+1))=-2\left[(c+3)(c+2)c+b^{2}(c-1)\right]<0,$ (88)

noting from Eq. (80) that $c>1$. Since $f$ is positive for $\lambda\geq 0$,
there must be a real root between $\lambda=-2(c+1)$ and $\lambda=0$. In the
case (86) we consider, $\lambda_{r}$ is the single real root, so we have

$\displaystyle-2(c+1)<\lambda_{r}<0.$ (89)

It then follows from (87) that $\alpha$ is strictly negative, completing the
proof that $\textrm{Re}[\lambda^{(i)}]<0$.

Perturbations away from equilibrium are thus exponentially damped. The
qualitative behavior of approach to equilibrium is determined by the mode (or
pair of complex conjugate modes) with smallest $|\textrm{Re}[\lambda^{(i)}]|$,
which decays the slowest. If this mode has a non-zero imaginary part (i.e., it
is part of a complex-conjugate pair), then oscillations will accompany the
decay, with frequency equal to $|\textrm{Im}[\lambda^{(i)}]|$. Note also that
the properties of the decay are insensitive to the sign of $b$, since the
characteristic polynomial (V.2) depends only on $b^{2}$.

It is interesting to ask for what parameter ranges such oscillations will be
important. The discriminant of the cubic (V.2) is equal to

$\displaystyle\Delta$ $\displaystyle=-\frac{4}{27}[(3b^{2}-c^{2}+10c-4)^{3}+$
(90) $\displaystyle(c^{3}-15c^{2}+30c+9b^{2}c-45b^{2}-8)^{2}]$ (91)

When $\Delta<0$ there are complex-conjugate roots; otherwise all roots are
real. Both signs occur over the parameter ranges $b\in\mathbb{R},c>1$, so both
behaviors are possible. If $\Delta\geq 0$, there will be no oscillations. If
$\Delta<0$ there will be oscillations, and these will be the dominant behavior
if $\alpha<\lambda_{r}$ in the notation of (86).

One tractable case is when $\delta\gg 1$ so that $c\approx 1$. It is easy to
see that the discriminant is negative for $c=1$, so that there will be
complex-conjugate roots. One can then check that these oscillatory modes decay
more slowly than the pure exponential mode provided $|b|\geq 2/\sqrt{3}$, and
that they decay with a similar order of magnitude even for $|b|\leq
2/\sqrt{3}$. The oscillations will thus always be important in the case
$\delta\gg 1$.

Some further details are helpful for interpreting this case. As a function of
$b$, the dominant decay time (real part of eigenvalue with real part closest
to zero) ranges between $-1$ at $b=0$ and $-4$ at $|b|\to\infty$, and the
frequency of oscillations (imaginary part of complex eigenvalue) ranges
between (approximately) $1.32$ at $b=0$ and infinity at large $|b|$,
approaching linearly like $|b|$ as $b\to\pm\infty$. Back in the physical time
coordinate $\tau=\tau_{E}T$, we see that decay will typically take several to
tens of $\tau_{E}$, and the oscillation angular frequency will be at least
$1.32\tau_{E}^{-1}$, and will be approximately
$|b|\tau_{E}^{-1}=\tau_{B}^{-1}$ at large $|b|$. It is notable that these
oscillations exist even in the pure-E case ($\delta\gg 1,b=0$) and are thus
physically distinct from synchrotron motion in that limit.

Figure 4: An example of entry into equilibrium. The field configuration is
circular with $\tilde{E}_{0}=1$ and $\tilde{B}_{0}=10$. The particle begins at
$\tilde{x}=(1,0,0)$ with initial momenta $\tilde{p}=(-2.32,8.28,3.97)\times
10^{4}$. The particle quickly loses its perpendicular momentum on timescale
$\tau_{\rm drop}\approx 10^{-4}\tau_{E}$ (68) before accelerating along the
PND to near the equilibrium Lorentz factor over several $\tau_{E}$ and finally
approaching the equilibrium in a weakly damped exponential fashion. Along with
the numerical trajectory, we show the corresponding analytical solutions in
the uniform-field and near-equilibrium approximations. The uniform field
solution (58)–(61) has 5 parameters $E_{0},B_{0},v_{i}$, all of which are
fixed by the initial data. The near-equilibrium solution (84) has six
parameters parameters $E_{0},B_{0},\kappa,C^{(i)}$. The first three are chosen
according to the local field at the initial position and the last three are
chosen by fitting with data from $\tau=15\tau_{E}$ to the end of the evolution
at $\tau=30\tau_{E}$. The oscillation period is approximately $.61\tau_{E}$,
and the damping time scale is approximately $10\tau_{E}$ (the exponential
envelope is $e^{-0.1\tau/\tau_{E}}$). Only the Lorentz factor was used for
these fits; however, using the fit parameters, the perpendicular velocities
$v_{n}$ and $v_{k}$ display a similar level of agreement.

###  V.3 Numerical examples

We now show two numerical examples illustrating the features explored in the
previous subsections. We consider the case of parallel circular electric and
magnetic fields, i.e., Eqs. (52) and (53) with $h=0$. Since the electric and
magnetic fields are parallel in the lab frame, no boost is required to relate
the parallel-frame uniform field solution (Sec. V.1) to the lab-frame near-
equilibrium solution (Sec. V.2). We always place the particle in a region of
the circular field configuration where $C_{4}\ll 1$ and hence equilibrium is
expected to occur.

First we consider the case where the particle begins with a large momentum
perpendicular to the PND. Based on the analysis of the previous subsections,
we expect the particle to sharply lose its perpendicular momentum on a
timescale $\tau_{\rm drop}$ (68), then gain parallel momentum exponentially on
a timescale $\tau_{E}$ (64) until it nears the equilibrium Lorentz factor, and
finally approach equilibrium exponentially, with timescale and possible
oscillations determined from $\tau_{E},\tau_{B},\tau_{\perp}$ via the
eigenvalue problem discussed in Sec. V.2. These expectations are confirmed by
our numerical experiments; an example is shown in Fig. 4.

We next consider the case where the particle begins with a large momentum
$\gamma\gg\gamma_{g}$ parallel to the PND. Here the uniform field
approximation is not useful since the particle quickly reaches the region
where the field line bends away from its motion. However, at this stage we can
regard the particle as traveling through a new uniform field configuration
with some perpendicular momentum, which will be removed by radiation reaction
according to the intuition discussed in the first new paragraph below Eq.
(65). In this way the particle will lose energy until it either enters
equilibrium directly (“from above”) or finds itself in a uniform field
configuration that will accelerate it back up to near-equilibrium conditions,
after which it enters equilibrium “from below”. An example somewhat
intermediate between these two cases is shown in Fig. 5.

Figure 5: Entry into equilibrium after starting with large momentum parallel
to the PND. The field setup is the same as Fig. 4, except that the particle
starts with momentum entirely in the $y$ direction with Lorentz factor equal
to 100 times the equilibrium value, $\gamma_{0}=100\gamma_{g}$. The particle
quickly loses energy until near the equilibrium value and then approaches
equilibrium with a sinusoidally modulated exponential.

##  VI Numerical survey

Up until now, we have explored the properties of the Aristotelian equilibrium
and the manner in which particles enter the equilibrium. We now turn to the
general question of when particles will indeed enter the equilibrium. As
reviewed in Sec. II.3, paper I identified five necessary conditions $C_{i}\ll
1$ for equilibrium to occur. We now provide numerical evidence that these
conditions are indeed sufficient for it to occur, and we use the results to
gain a more quantitative understanding of how well they must be satisfied.

Exploring all five conditions would be computationally intractable. However,
given that $C_{4}\ll 1$ implies $C_{1}\ll 1$ and $C_{2}\ll 1$ outside of
extremely finely tuned field configurations (Eqs. (35) and (36)), we can
safely ignore $C_{1}$ and $C_{2}$. The case where $C_{3}\ll 1$ is violated is
potentially interesting for helical fields (as occur, e.g., in relativistic
jets), and we found that equilibrium does not occur in this case (Sec. IV).
Violations of the quasi-static assumption $C_{5}\ll 1$ are certainly
interesting (e.g., for laser fields), but outside the scope of this paper.
Instead, we will pick field configurations where $C_{3}=0$ and $C_{5}=0$
(exactly torsion-free and static, respectively) and explore the role of
$C_{4}$ in determining whether equilibrium occurs.

###  VI.1 Circular fields and the role of $\eta$

Eq. (34) defined a quantity $\eta$ that appears in the conditions for
equilibrium. For generic fields, $\eta\approx 1$ and factors of $\eta$ can be
dropped. To illustrate a situation where $\eta$ cannot be dropped, consider a
purely azimuthal field configuration with constant parallel electric and
magnetic fields. The field lines and PNDs are just circles of radius
$\rho=\sqrt{x^{2}+y^{2}}$. Since $\vec{\nabla}\rho=\vec{n}$, where $\vec{n}$
is the Frenet-Serret normal direction to the curve, $\eta$ can be calculated
as

$\displaystyle\eta=|\vec{v}\cdot\nabla\rho|=|v_{n}|.$ (92)

This quantity is small in equilibrium since $|v_{n}|\ll 1$ by the assumption
of motion along a PND. The small size of $\eta$ arises because the field
strength and radius of curvature are both constant along the PND direction
$\vec{\ell}$, a finely tuned arrangement. From Eqs. (26) and (25), we have

$\displaystyle|v_{n}|^{2}=\frac{2}{3}\frac{\mathcal{R}}{R}\sqrt{\frac{E_{0}}{\mathcal{E}}}(1+\delta)^{2}.$
(93)

The condition (37) then becomes

$\displaystyle R^{2}E_{0}\gg(1+\delta)^{2}\mathcal{R}^{2}\mathcal{E},$ (94)

where we drop a factor of $2/3$.

Eq. (94) is a necessary condition for equilibrium to occur. To check whether
it is also sufficient, we simulated a large number of trajectories with
different field configuration parameters and particle initial conditions.
Specifically, we fixed $\tilde{B}_{0}=0.1$ and chose the particle to begin on
the $x$ axis, choosing the other parameters randomly from uniform
distributions in the ranges $\log_{10}\tilde{E}_{0}\in(-4,0)$ for the log of
the field strength, $\gamma\in(1,\gamma_{g})$ for the initial Lorentz factor,
$\theta\in(0,\pi)$ and $\phi\in(0,2\pi)$ for the initial direction of motion
(where $\theta,\phi$ are spherical coordinates on the space of velocities),
and $\log_{10}x\in(-3,0)$ for the log of the initial position. The initial
radius of curvature $\tilde{R}$ is just $x$, so we sample an initial range
$\log_{10}\tilde{R}\in(-3,0)$.

For each run of the code we evolve for a maximum time of $\tau_{\rm
max}=(6\log\gamma_{g})\tau_{E}$, which is sufficient for the particle to enter
equilibrium if it is destined to do so.444According to Eq. (67), the time for
a particle to be accelerated to $\gamma_{g}$ is of order
$\tau_{E}\log\gamma_{g}$. For the circular field $E_{0}$ and hence $\tau_{E}$
is constant; however, for other field configurations we update the maximum
allowed run time after each time step to use the local value of $E_{0}$. We
also examined a selection of the trajectories that did not enter equilibrium
in this time and evolved them for longer, finding that indeed they never enter
equilibrium. At any given time in the trajectory, the particle is considered
to have entered equilibrium if
$\langle|(\gamma-\gamma_{g})|/\gamma_{g}\rangle<3\%$, where the average
$\langle...\rangle$ is calculated over the last $20\%$ of the computed
trajectory, using the local value of $\gamma_{g}$. Once a particle enters
equilibrium, we record the local electric field $E_{0}$ and curvature radius
$R$ and terminate the run. (If the particle never enters equilibrium, we make
no corresponding record.) Fig. 6 shows the values of $E$ and $R$ at which
equilibrium occurred in our random sample, showing clearly that the condition
(94) indeed controls whether equilibrium is obtained. We find that the left-
hand-side must be approximately 15 times larger than the right-hand side (red
line in the figure) for equilibrium to occur.

When $E_{0}\gg B_{0}$, we have $\delta\approx\mathcal{E}/E_{0}\gg 1$, and this
condition reduces to the original Gruzinov condition (21). (This agreement is
accidental, due to the specific form of $\eta$ for this field configuration.)
By contrast, when $E_{0}\lesssim B_{0}$, the condition (94) remains distinct
from the Gruzinov condition (21). This behavior is seen clearly in Fig. 6. The
Gruzinov condition is necessary, but not sufficient, in this special case. The
correct condition for entry into equilibrium is the modified condition (37).

Figure 6: Values of the local electric field and PND curvature radius at entry
into equilibrium in a large parameter survey. In this circular field
configuration, equilibrium is expected to occur when
$N=R^{2}E_{0}/[(1+\delta)^{2}\mathcal{R}^{2}\mathcal{E}]\gg 1$ [Eq. (94)]. The
numerics validate this condition: the red line is the curve $N=15$.

###  VI.2 More general field configurations

The circular field is a special case that cleanly illustrates the role of
$\eta$ in the condition $C_{4}\ll 1$ (32) for entry into equilibrium. In order
to test the condition more generally, we consider two more kinds of field
configuration. The “ellipse” field configuration simply changes the shape of
the field lines to ellipses instead of circles, while keeping everything else
the same. The “separate center” field configuration consists of circular
electric and magnetic field lines with equal field strengths $E_{0}=B_{0}$ but
with the center of the circles shifted by a distance $d$. We randomly varied
the field strength and center-distance while also randomly choosing initial
conditions as before.

The results of these surveys are shown along with the circular case in Fig. 7,
demonstrating clearly that the condition $C_{4}\ll 1$ is necessary and
sufficient for equilibrium to occur. The precise value of $C_{4}$ required
depends on the field configuration, ranging from $\sim 0.1$ to $\sim 0.01$ for
the cases considered here.

Figure 7: Validity of the main condition for entry into equilibrium [$C_{4}\ll
1$, (32) or (37)] in three numerical parameter surveys. The colored dots
correspond to the local field properties where a particle in the corresponding
survey entered equilibrium. Lines of constant $C_{4}$ have slope $-1$ on this
plot and the $y$-intercept is the log of $1/C_{4}^{2}$. For each field
configuration, the edge of the region of equilibria clearly has the predicted
slope of $-1$. We have drawn approximate reference lines for this edge to
guide help guide the eye. Entry to equilibrium occurs when $C_{4}\lesssim
5\%$.

##  VII Summary of results

We have provided a detailed understanding of the entry into Aristotelian
equilibrium, identifying the relevant timescales, giving analytical
descriptions where possible, and showing that numerical simulations match the
predicted behavior. We showed analytically that the equilibrium is linearly
stable and identified the presence of oscillations at entry for a large region
of parameter space. We performed numerical parameter surveys exploring the
conditions for equilibrium to occur. Combined with the results of paper I,
this study provides strong evidence that the conditions (29)–(33) are
necessary and sufficient for Eqs. (25)–(27) to describe the motion of a
charged particle. This provides a solid foundation for using the Aristotelian
approximation in astrophysical modeling.

## Acknowledgements

This work was supported in part by NSF grant PHY–1752809, and NSF Grants
PHY-1912619 and PHY-2145421, to the University of Arizona.

##  Appendix A Numerical Method

Our numerical scheme is based on the implicit 4th order Runge-Kutta-Nyström
method. We first review this method before presenting an improved iteration
method and an adaptive timestep version that we also used in some cases.

For a second-order ordinary differential equation

$\displaystyle\ddot{x}^{\alpha}=f^{\alpha}(t,x^{\alpha},\dot{x}^{\alpha}),$
(95)

the implicit 4th order Runge-Kutta-Nyström method can be expressed as follows.
First fix a time step $h$. If $x_{n}^{\alpha}$ and $\dot{x}_{n}^{\alpha}$
represent the value and derivative at time $t_{n}$, then the values at the
next step $t_{n+1}=t_{n}+h$ are obtained by

$\displaystyle x_{n+1}^{\alpha}$
$\displaystyle=x_{n}^{\alpha}+h\dot{x}_{n}^{\alpha}+h^{2}b_{i}k^{\alpha}_{i}$
(96) $\displaystyle\dot{x}_{n+1}^{\alpha}$
$\displaystyle=\dot{x}_{n}^{\alpha}+ha_{i}k^{\alpha}_{i},$ (97)

where $a_{i}$ and $b_{i}$ are

$\displaystyle a_{i}$ $\displaystyle=\left(\frac{1}{2},\ \frac{1}{2}\right)$
(98) $\displaystyle b_{i}$
$\displaystyle=\left(\frac{1}{4}+\frac{\sqrt{3}}{12},\
\frac{1}{4}-\frac{\sqrt{3}}{12}\right),$ (99)

and $k^{\alpha}_{i}$ must satisfy

$\displaystyle k^{\alpha}_{i}=f^{\alpha}$
$\displaystyle\big{(}t_{n}+c_{i}h,x^{\alpha}_{n}+c_{i}h\dot{x}_{n}^{\alpha}+h^{2}B_{ij}k^{\alpha}_{j},\dot{x}_{n}^{\alpha}+hA_{ij}k^{\alpha}_{j}\big{)}.$
(100)

In these expressions, repeated indices are summed. Here $c_{i}$ is the vector

$\displaystyle c_{i}=\begin{pmatrix}\frac{1}{2}-\frac{\sqrt{3}}{6}\\\
\frac{1}{2}+\frac{\sqrt{3}}{6}\end{pmatrix}$ (101)

and $A_{ij}$ and $B_{ij}$ are the matrices

$\displaystyle A_{ij}$
$\displaystyle=\begin{pmatrix}\frac{1}{4}&\frac{1}{4}-\frac{\sqrt{3}}{6}\\\
\frac{1}{4}+\frac{\sqrt{3}}{6}&\frac{1}{4}\end{pmatrix}$ (102) $\displaystyle
B_{ij}$
$\displaystyle=\begin{pmatrix}\frac{1}{36}&\frac{5}{36}-\frac{\sqrt{3}}{12}\\\
\frac{5}{36}+\frac{\sqrt{3}}{12}&\frac{1}{36}\end{pmatrix}.$ (103)

### Iteration method

Eq. (100) can be solved by a fixed point iteration method as follows. The
initial values of $k^{\alpha}_{i}$ are taken to be

$\displaystyle(k^{\alpha}_{i})^{0}=f^{\alpha}(t_{n},x^{\alpha}_{n},\dot{x}_{n}^{\alpha})\
\ (i=1,2).$ (104)

We then iteratively improve the guess by calculating

$\displaystyle(k^{\alpha}_{i})^{N+1}=f^{\alpha}$
$\displaystyle\big{(}t_{n}+c_{i}h,$ $\displaystyle
x^{\alpha}_{n}+c_{i}h\dot{x}_{n}^{\alpha}+h^{2}B_{ij}(k^{\alpha}_{j})^{N}$
$\displaystyle\dot{x}_{n}^{\alpha}+A_{ij}(k^{\alpha}_{j})^{N}\big{)}.$ (105)

However, for large Lorentz factors we found that this iteration method
sometimes converged very slowly or got stuck in a limit cycle. In these cases
we instead used overrelaxation with the secant method,

$\displaystyle y_{n}$
$\displaystyle=w\left[y_{n-1}-F(y_{n-1})\frac{y_{n-1}-y_{n-2}}{F(y_{n-1})-F(y_{n-2})}\right]$
(106) $\displaystyle+(1-w)y_{n-1},$ (107)

where $y_{N}=(k_{i}^{\alpha})^{N}$,
$F^{\alpha}_{i}(k_{i}^{\alpha})=f^{\alpha}(k^{\alpha}_{i})-k_{i}^{\alpha}$ and
$w\in[0,1]$ is the weight. We chose $w=0.05$. At each time step we try the
simple iteration method first, jumping to the Secant method only if the simple
iteration fails to reach a given tolerance (we chose fractional error of
$10^{-10}$) within a fixed number of iterations (we chose $10^{3}$).

###  A.1 Adaptive time step

For most portions of a given evolution the main timescales are $\tau_{E}$ and
$\tau_{B}$, defined in Eqs. (64) and (63). However, for some initial
conditions (e.g., Fig. 4) there can be abrupt changes in Lorentz factor on the
much smaller timescale $\tau_{\rm drop}$ defined in Eq. (68). In most cases we
use an adaptive time step scheme, beginning with
$h=\rm{min}(\tau_{B},\tau_{E})\times 1\%$ and updating as follows.

We focus on the Lorentz factor as a simple scalar representative of the
solution. Suppose we are at time $t$ in the evolution and consider the value
at time $t+2h$, where $h$ is the current step size. We can compute this either
by taking a single step of size $2h$ or by taking two steps of size $h$. We
will denote the resulting values for $\gamma$ by $\gamma_{\rm num}^{(2h)}$ and
$\gamma_{\rm num}^{(h)}$, respectively. Since our method has fourth-order
accuracy, these are related to the true value $\gamma_{\rm true}$ by

$\displaystyle\gamma_{\rm num}^{(h)}$ $\displaystyle\approx\gamma_{\rm
true}(t+2h)+\phi h^{4}$ (108) $\displaystyle\gamma_{\rm num}^{(2h)}$
$\displaystyle\approx\gamma_{\rm true}(t+2h)+16\phi h^{4},$ (109)

where $\phi$ is some unknown constant. We can estimate $\phi$ by solving,

$\displaystyle\phi\approx\frac{\gamma_{\rm num}^{(2h)}-\gamma_{\rm
num}^{(h)}}{15h^{4}}.$ (110)

Suppose we instead consider a new time step $h_{0}$ and evolve forward one
step to $\gamma^{(h_{0})}_{\rm num}$. Now we have

$\displaystyle\gamma^{(h_{0})}_{\rm num}\approx\gamma_{\rm true}(t+2h)+\phi
h_{0}^{4},$ (111)

If we want to achieve an accuracy of $\epsilon=\left|(\gamma_{\rm
num}^{(h_{0})}-\gamma_{\rm true})/\gamma_{\rm true}\right|$, the time step we
should use satisfies

$\displaystyle h_{0}^{4}=\frac{\epsilon\gamma_{\rm num}^{(h)}}{|\phi|}.$ (112)

Using the formula (110) for $\phi$, we conclude that an appropriate time step
is

$\displaystyle h_{0}=\zeta h\left(\frac{15\epsilon\gamma_{\rm
num}^{(h)}}{|\gamma_{\rm num}^{(2h)}-\gamma_{\rm num}^{(h)}|}\right)^{1/4}.$
(113)

The factor of $\zeta<1$ is to guarantee that $h_{0}$ is smaller than that
which is expected to exactly achieve the desired error tolerance $\epsilon$,
which corresponds to $\zeta=1$. We choose $\zeta=(14/15)^{1/4}$ in our code,
and $\epsilon=10^{-6}$.

At each step in the evolution we calculate the candidate new time step $h_{0}$
according to (113). Before adopting this as the new time step we consider two
potential adjustments. First, we ensure that $h$ does not change by more than
a factor of two. That is, if $h_{0}$ is larger than $2h$, then we use $2h$
instead; similarly, if $h_{0}$ is smaller than $h/2$, we use $h/2$ instead.
Finally, we ensure that the candidate new time step satisfies the condition
for the convergence of fixed point iteration, namely that the spectral radius
of the (six-dimensional) Jacobian matrix of $f^{\alpha}_{i}(k^{\beta}_{j})$ is
strictly less than $1$, i.e., the maximum of the absolute values of the
eigenvalues of this matrix is less than 1. (Here $f^{\alpha}_{i}$ denotes the
RHS of (100).) If this test fails, then we divide the step size in half and
try again, iterating until a step size satisfying the condition has been
found.

Figure 8: Test demonstrating our code’s error convergence at 4th order. The
black line represents the linear function with slope of 4; the equation given
in the plot.

###  A.2 Convergence test

We performed a convergence test for the code using the uniform field setup,
because it has exact analytical solutions. We choose the field strengths as
$\tilde{E}_{0}=0.1,\tilde{B}_{0}=1$. The initial value we used is
$\tilde{\vec{p_{0}}}=\\{100,400,-300\\}$. The Lorentz factor Eq. (58) at
$\tau=2\tau_{E}$ was used as the reference value. The code was executed at
different time steps $h$. The fractional difference:

$\displaystyle\epsilon=\left|\frac{\gamma_{\rm numeric}-\gamma_{\rm
analytic}}{\gamma_{\rm analytic}}\right|$ (114)

was then calculated. The result of the convergence test is shown in Fig. 8. It
is clear that our code exhibits 4th-order convergence.

## References

  * [1] J. L. Synge, Relativity: the special theory.  North-Holland Publishing Company Amsterdam, 1956.
  * [2] L. Mestel, J. Robertson, Y.-M. Wang, and K. Westfold, “The axisymmetric pulsar magnetosphere,” Monthly Notices of the Royal Astronomical Society 217 (1985), no. 3, 443–484.
  * [3] B.Finkbeiner, H.Herold, T. Ertl, and H. Ruder Astronomy and astrophysics 225 (1989) 479.
  * [4] A. Gruzinov, “Pulsar Emission Spectrum,” arXiv e-prints (Sept., 2013) arXiv:1309.6974, 1309.6974.
  * [5] T. Jacobson, “Structure of Aristotelian Electrodynamics,” Phys. Rev. D 92 (2015), no. 2, 025029, 1504.07311.
  * [6] J. Pétri, “Pulsar gamma-ray emission in the radiation reaction regime,” Mon. Not. Roy. Astron. Soc. 484 (2019), no. 4, 5669–5691, 1901.11439.
  * [7] G. Cao and X. Yang, “Three-dimensional dissipative pulsar magnetospheres with Aristotelian electrodynamics,” 1912.00335.
  * [8] A. Gonoskov and M. Marklund, “Radiation-dominated particle and plasma dynamics,” Phys. Plasmas 25 (2018), no. 9, 093109, 1707.05749.
  * [9] A. S. Samsonov, E. N. Nerush, and I. Y. Kostyukov, “Asymptotic electron motion in the strongly-radiation-dominated regime,” Phys. Rev. A 98 (2018), no. 5, 053858, 1807.04071.
  * [10] R. Ekman, T. Heinzl, and A. Ilderton, “Exact solutions in radiation reaction and the radiation-free direction,” New J. Phys. 23 (2021), no. 5, 055001, 2102.11843.
  * [11] A. Gonoskov, T. G. Blackburn, M. Marklund, and S. S. Bulanov, “Charged particle motion and radiation in strong electromagnetic fields,” 2107.02161.
  * [12] E. A.S.Samsonov and I.Yu.Kostyukov, “High-order corrections to the radiationfree dynamics of an electron in the strongly radiation-dominated regime,” Matter and Radiation at Extremes 8 (2022), no. 1, 014402, 2208.00673.
  * [13] Y. Cai, S. E. Gralla, and V. Paschalidis, “Dynamics of ultrarelativistic charged particles with strong radiation reaction. I. Aristotelian equilibrium state,” 2209.07469.
  * [14] L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields: Volume 2 (Course of Theoretical Physics Series) 4th Edition.  Butterworth-Heinemann, 1980.
  * [15] A. Gruzinov, “Aristotelian Electrodynamics solves the Pulsar: Lower Efficiency of Strong Pulsars,” 1303.4094.
  * [16] H. Heintzmann and E. Schrüfer, “Exact solutions of the Lorentz-Dirac equations of motion for charged particles in constant electromagnetic fields,” Phys. Rev. A 43 (1973), no. 3, 287–288.
  * [17] M. C. Begelman, R. D. Blandford, and M. J. Rees, “Theory of extragalactic radio sources,” Reviews of Modern Physics 56 (Apr., 1984) 255–351.


# Stochastic Time-Optimal Control Studies for Additional Food provided Prey-
Predator Systems involving Holling Type-IV Functional Response

D Bhanu Prakash  First Author. Email<EMAIL_ADDRESS>D K K Vamsi
Corresponding author. Email<EMAIL_ADDRESS>
###### Abstract

We consider an additional food provided prey-predator model exhibiting Holling
type IV functional response with combined continuous white noise and
discontinuous Lévy noise. We prove the existence and uniqueness of global
positive solutions for the considered model. By considering the quality and
quantity of additional food as control parameters, we formulate a time-optimal
control problem. We obtain the condition for the existence of an optimal
control. Furthermore, making use of the arrow condition of the sufficient
stochastic maximum principle, we characterize the optimal quality of
additional food and optimal quantity of additional food. Numerical results are
given to illustrate the theoretical findings with applications in biological
conservation and pest management.

keywords:  Stochastic optimal control; Time-optimal control; Holling type IV
response; Biological conservation; Pest management; Brownian motion; Lévy
noise

MSC 2020 codes:  37A50; 60H10; 60J65; 60J70; 60J76; 49K45;

##  1 Introduction

The complex natural ecosystems present around us kindled great attention of
many ecologists and mathematicians to the mathematical modeling of ecological
systems in the last few decades. The interaction among species in these
ecosystems can be of several forms like competition, mutual interference,
prey-predator interactions and so on. The very first ecological models are
framed from the pioneering works of Alfred J. Lotka [1] and Vito Volterra [2]
in 1925. Various complex models are framed and studied ever since.

The basic component of prey-predator systems is a functional response, which
is defined as the rate at which each predator captures prey [3]. A few
examples of functional responses include Holling functional responses [4],
Beddington–DeAngelis functional responses [5], Arditi-Ginzburg functional
responses [6], Hassell-Varley functional responses [7] and Crowley & Martin
functional response [8]. In this work, we study the models exhibiting Holling
type IV functional response. Some of the organisms that display Holling type
IV functional response in nature are found in [9, 10]. In this functional
response, the predator’s per capita rate of predation decreases at
sufficiently high prey density, due to either prey interference or prey
toxicity.

In recent decades, many pioneering dynamical modeling works [11, 12, 13, 14,
15, 16] reveal that the provision of additional food to predators plays a
vital role in controlling the dynamics of the system. Recently, the authors in
[14] have studied additional food provided deterministic prey-predator systems
involving Holling type IV functional response. Also, the authors in [17, 16]
studied the optimal control problems of deterministic prey-predator systems
involving Holling type IV functional response with the quality of additional
food and the quantity of additional food as the control parameters
respectively.

Often it is observed that the parameters in an ecosystem are effected by the
environmental fluctuations [18]. Therefore in recent years, many researchers
have drawn their attention to stochastic models. Most stochastic prey-predator
models are driven by the Brownian motion, which captures the continuous noise.
Authors in [19, 20] studied the deterministic and stochastic dynamics of a
modified Leslie-Gower prey-predator system with simplified Holling type IV
functional response. However, the sudden changes in environment like toxic
pollutants, floods, earthquakes and so on, cannot be captured by the Brownian
motion as it is a continuous noise. Hence, addition of a discontinuous noise,
like Lévy noise, to the prey-predator system with Brownian motion makes the
models more realistic. [21] uses the stochastic averaging method to analyze
the modified stochastic Lotka-Volterra models under combined Gaussian and
Poisson noise. [22] studies the dynamics and dynamics of a Stochastic One-
Predator-Two-Prey time delay system with jumps.

To the best of our knowledge, there is no study of additional food provided
stochastic prey-predator system with jumps. Secondly, the optimal control
studies of Stochastic Differential Equations with Jumps (SDEJ) were not
performed on prey-predator systems. Lastly, very few works involved Holling
type-IV response which incorporates the most important group defence property.
Motivated by these observations, in this work, we study the optimal control
problems for additional food provided stochastic Holling type IV prey-predator
systems under combined Gaussian and Lévy noise.

The article is structured as follows: Section 2 introduces the stochastic
prey-predator model with Holling type-IV functional response and additional
food with intra-specific competition among predators. The existence of global
positive solution for this model is briefly discussed in Section 3. The time-
optimal control problem is formulated and the optimal quality and quantity of
additional food is characterized in Section 4. Section 5 illustrates the key
findings of the analysis through numerical simulations in the context of both
biological conservation and pest management. Finally, Section 6 presents the
discussions and conclusions.

##  2 The Stochastic Model

Let $N$ and $P$ denote the biomass of prey and predator population densities
respectively. In the absence of predator, the prey growth is modelled using
logistic equation. Further, we assume that the prey species exhibit Holling
type-IV functional response towards predators. We also assume that the
predators are supplemented with an additional food of biomass A, which is
uniformly distributed in the habitat. Incorporating these assumptions, the
prey-predator dynamics with Holling type-IV functional response along with
additional food for predators can be described as:

$\begin{split}\frac{\mathrm{d}N(t)}{\mathrm{d}t}&=rN(t)\left(1-\frac{N(t)}{K}\right)-\Bigg{(}\frac{cN(t)}{(A\eta\alpha+a)(bN^{2}(t)+1)+N(t)}\Bigg{)}P(t)\\\
\frac{\mathrm{d}P(t)}{\mathrm{d}t}&=e\Bigg{(}\frac{N(t)+\eta
A(bN^{2}(t)+1)}{(A\eta\alpha+a)(bN^{2}(t)+1)+N(t)}\Bigg{)}P(t)-m_{1}P(t)\end{split}$
(1)

For a complete analysis of model (1), the reader is referred to Vamsi et. al.,
[15].

In addition, we also assume that the predators exhibit intra-specific
competition. We capture this competition in similar lines with [23, 24].
Accordingly, the system (1) gets transformed to the following system.

$\begin{split}\frac{\mathrm{d}N(t)}{\mathrm{d}t}&=rN(t)\left(1-\frac{N(t)}{K}\right)-\Bigg{(}\frac{cN(t)}{(A\eta\alpha+a)(bN^{2}(t)+1)+N(t)}\Bigg{)}P(t)\\\
\frac{\mathrm{d}P(t)}{\mathrm{d}t}&=e\Bigg{(}\frac{N(t)+\eta
A(bN^{2}(t)+1)}{(A\eta\alpha+a)(bN^{2}(t)+1)+N(t)}\Bigg{)}P(t)-m_{1}P(t)-\delta
P(t)^{2}\end{split}$ (2)

The biological descriptions of the various parameters involved in the systems
(1) and (2) are described in Table 1.

Parameter | Definition | Dimension
---|---|---
T | Time | time
N | Prey density | biomass
P | Predator density | biomass
r | Prey intrinsic growth rate | time-1
K | Prey carrying capacity | biomass
c | Maximum rate of predation | time-1
e | Maximum growth rate of predator | time-1
m1 | Predator mortality rate | time-1
$\delta$ | Death rate of predators due to intra-specific competition | biomass-1 time-1
A | Quantity of additional food for predators | biomass
b | Group defence in prey | biomass-2
Table 1: Description of variables and parameters present in the systems (1),
(2)

In order to reduce the complexity of the model, we non-dimensionalize the
system (2) using the following non-dimensional parameters.

$N=ax,\ P=\frac{ay}{c},\ \gamma=\frac{K}{a},\ \xi=\frac{\eta A}{a},\
\omega=ba^{2},\ m_{2}=\frac{c}{a\delta}$

. Accordingly, system (2) gets reduced to the following system.

$\begin{split}\frac{\mathrm{d}x}{\mathrm{d}t}&=rx\Bigg{(}1-\frac{x}{\gamma}\Bigg{)}-\Bigg{(}\frac{xy}{(1+\alpha\xi)(\omega
x^{2}+1)+x}\Bigg{)}\\\
\frac{\mathrm{d}y}{\mathrm{d}t}&=e\Bigg{(}\frac{x+\xi(\omega
x^{2}+1)}{(\alpha\xi+1)(\omega
x^{2}+1)+x}\Bigg{)}y-m_{1}y-m_{2}y^{2}\end{split}$ (3)

In real world scenarios, environmental fluctuations affect the dynamics of the
system. In order to capture these fluctuations, we introduce the
multiplicative white noise terms into (3). As in [25, 26, 23], we now suppose
that the intrinsic growth rate of prey and the death rate of predator are
mainly affected by environmental noise such that

$r\rightarrow r+\sigma_{1}dW_{1}(t),\ m_{1}\rightarrow
m_{1}+\sigma_{2}dW_{2}(t)$

where $W_{i}(t)\ (i=1,2)$ are the mutually independent standard Brownian
motions with $W_{i}(0)=0$ and $\sigma_{1}$ and $\sigma_{2}$ are positive
constants and they represent the intensities of the white noise.

Also, the system can go through huge, occasionally catastrophic disturbances.
Since white noise is a continuous noise, it cannot capture sudden
environmental changes. To cater to these, we also apply a discontinuous
stochastic process as Lévy jumps to model these abrupt natural phenomenon as
in [27, 22].

We now perturb $r$ and $m_{1}$ with discontinuous Lévy noise in addition to
the continuous white noise. So, we have

$r\rightarrow
r+\sigma_{1}\mathrm{d}W_{1}(t)+\int_{\mathbb{Y}}\gamma_{1}(v)\,\widetilde{N}(\mathrm{d}t,\mathrm{d}v),\
-m_{1}\rightarrow-
m_{1}+\sigma_{2}\mathrm{d}W_{2}(t)+\int_{\mathbb{Y}}\gamma_{2}(v)\,\widetilde{N}(\mathrm{d}t,\mathrm{d}v)$
(4)

According to the Lévy decomposition theorem [28], we have
$\widetilde{N}(t,dv)=N(t,dv)-\lambda(dv)t$, where $\widetilde{N}(t,dv)$ is a
compensated Poisson process and N is a Poisson counting measure with
characteristic measure $\lambda$ on a measurable subset $\mathbb{Y}$ of
$(0,+\infty)$ with $\lambda(\mathbb{Y})<\infty$. The distribution of Lévy
jumps Li(t) can be completely parameterized by $(a_{i},\sigma_{i},\lambda)$
and satisfies the property of infinite divisibility.

Now, by incorporating noise induced parameters (4) into the reduced
deterministic system of equations (3), we get the following additional food
provided stochastic prey-predator system exhibiting Holling type-IV functional
response along with the environmental fluctuations captured using the white
noise and Lévy noise.

$\begin{split}\mathrm{d}x(t)&=x(t)\Bigg{[}r\Bigg{(}1-\frac{x(t)}{\gamma}\Bigg{)}-\Bigg{(}\frac{y(t)}{(1+\alpha\xi)(\omega
x^{2}(t)+1)+x(t)}\Bigg{)}\Bigg{]}\mathrm{d}t+\sigma_{1}x(t)\mathrm{d}W_{1}(t)+x(t)\int_{\mathbb{Y}}\gamma_{1}(v)\,\widetilde{N}(\mathrm{d}t,\mathrm{d}v)\\\
\mathrm{d}y(t)&=y(t)\Bigg{[}e\Bigg{(}\frac{x(t)+\xi(\omega
x^{2}(t)+1)}{(1+\alpha\xi)(\omega
x^{2}(t)+1)+x(t)}\Bigg{)}-m_{1}-m_{2}y(t)\Bigg{]}+\sigma_{2}y(t)\mathrm{d}W_{2}(t)+y(t)\int_{\mathbb{Y}}\gamma_{2}(v)\,\widetilde{N}(\mathrm{d}t,\mathrm{d}v)\end{split}$
(5)

##  3 Existence of Global Positive Solution

In order to study the stochastic dynamics of (5), we first prove that the
system (5) has a unique global positive solution.

######  Theorem 1.

For any given initial value $X(\theta)=(x(\theta),y(\theta))\in
C([-\tau_{0},0],\mathbb{R}^{+^{2}})$, there exists a unique positive global
solution $(x(t),y(t))$ of system (5) on $t\geq 0$.

The above theorem for existence of solutions of (5) can be proved in similar
lines to the proof in [22] using the Lyapunov method.

##  4 Stochastic Time-Optimal Control Problems

In this section, we formulate and study the stochastic time-optimal control
problems for prey-predator systems involving Holling type-IV functional
response where the predator is provided with additional food.

###  4.1 Quality of additional food as control parameter

In this subsection, we characterize the optimal quality of additional food for
driving the system (5) to a desired equilibrium state in minimum time using
the stochastic maximum principle. So, we fix the quantity of additional food
$\xi>0$ to be a constant and choose the objective functional to be minimized
for this stochastic time optimal control problem as follows.

$J(\alpha)=E\Bigg{[}\int_{0}^{T}1dt\Bigg{]}.$ (6)

From the Sufficient Stochastic Maximum Principle [29] for the optimal control
problems of jump diffusion, we characterize the optimal solution of the
stochastic time optimal control problem with state space as the solutions of
(5) and the objective functional (6).

Let $(p^{*},q^{*},r^{*})$ be a solution of the adjoint equation in the unknown
processes $p(t)\in\mathbb{R}^{2},\ q(t)\in\mathbb{R}^{2\times 2},\
r(t,z)\in\mathbb{R}^{2}$ satisfying the backward differential equations

$\begin{split}dp_{1}(t)&=\Bigg{[}\left(-r+\frac{2rx}{\gamma}-\frac{2\omega(1+\alpha\xi)x+1}{((1+\alpha\xi)(\omega
x^{2}+1)+x)^{2}}\right)p_{1}(t)-\frac{(1-\omega
x^{2})(1+(\alpha-1)\xi)}{((1+\alpha\xi)(\omega
x^{2}+1)+x)^{2}}eyp_{2}(t)-\sigma_{1}q_{1}\\\
&-\int\gamma_{1}(v)r_{1}v_{1}(dz_{1})\Bigg{]}dt+q_{1}(t)dW_{1}(t)+q_{2}(t)dW_{2}(t)+\int_{\mathbb{R}^{n}}r_{1}\widetilde{N}(dt,dz)\\\
dp_{2}(t)&=-\Bigg{[}\frac{-x}{(1+\alpha\xi)(\omega
x^{2}+1)+x}p_{1}(t)+\left(\frac{e(x+\xi(\omega x^{2}+1))}{(1+\alpha\xi)(\omega
x^{2}+1)+x}-m_{1}-2m_{2}y\right)p_{2}(t)+\sigma_{2}q_{4}\\\
&+\int\gamma_{2}(v)r_{2}v_{2}(dz_{2})\Bigg{]}dt+q_{3}(t)dW_{1}(t)+q_{4}(t)dW_{2}(t)+\int_{\mathbb{R}^{n}}r_{2}\widetilde{N}(dt,dz)\\\
p_{1}(T)&=0,\ p_{2}(T)=0\end{split}$ (7)

The Hamiltonian associated with this control problem is defined as follows.

$\begin{split}H(t,x,y,\alpha,p,q,r)&=1+\Bigg{[}r(1-\frac{x}{\gamma})-\frac{y}{(1+\alpha\xi)(1+\omega
x^{2})+x}\Bigg{]}xp_{1}+\Bigg{[}\frac{e(x+\xi(\omega
x^{2}+1))}{(1+\alpha\xi)(1+\omega x^{2})+x}\\\
&-m_{1}-m_{2}y\Bigg{]}yp_{2}+\sigma_{1}xq_{1}+\sigma_{2}yq_{4}+x\int\gamma_{1}r_{1}v_{1}(dz_{1})+y\int\gamma_{2}(v)r_{2}v_{2}(dz_{2})\end{split}$
(8)

Let $U=\\{\alpha(t)|0\leq\alpha(t)\leq\alpha_{max}\forall t\in(0,t_{f}]\\}$
where $\alpha_{max}\in\mathbb{R}^{+}.$ Let $\alpha^{*}\in U$ with the
corresponding solution $(x^{*},y^{*})=(x(u^{*}),y(u^{*}))$.

From the Arrow condition in the sufficient stochastic maximum Principle [29],
we have

$\begin{split}\frac{\partial H}{\partial\alpha}\Big{|}_{\alpha^{*}}=0&\\\
\implies&\Bigg{[}-yxp_{1}\frac{-\xi(\omega x^{2}+1}{((1+\alpha\xi)(\omega
x^{2}+1)+x)^{2}}-\frac{eyp_{2}(x+\xi(\omega x^{2}+1))(\xi(\omega
x^{2}+1))}{((1+\alpha\xi)(\omega x^{2}+1)+x)^{2}}\Bigg{]}\Bigg{|}_{*}=0\\\
\implies&\Bigg{[}yxp_{1}-eyp_{2}(x+\xi(\omega
x^{2}+1))\Bigg{]}\Bigg{|}_{*}=0\\\
\implies&x^{*}p_{1}^{*}=ep_{2}^{*}(x^{*}+\xi(\omega x^{*^{2}}+1))\\\
\end{split}$

Hence the optimal control $\alpha^{*}$ should satisfy the following condition.

$x^{*}p_{1}^{*}=ep_{2}^{*}(x^{*}+\xi(\omega x^{*^{2}}+1))$ (9)

Since the analytical solution of (7) is complex to solve, we numerically
simulate these results in section 5.

###  4.2 Quantity of additional food as control parameter

In this subsection, we characterize the optimal quantity of additional food
for driving the system (5) to a desired equilibrium state in minimum time
using the stochastic maximum principle. So, we fix the quality of additional
food $\alpha>0$ to be a constant and choose the objective functional to be
minimized for this stochastic time optimal control problem as follows.

$J(\xi)=E\Bigg{[}\int_{0}^{T}1dt\Bigg{]}.$ (10)

From the Sufficient Stochastic Maximum Principle [29] for the optimal control
problems of jump diffusion, we characterize the optimal solution of the
stochastic time optimal control problem with state space as the solutions of
(5) and the objective functional (10).

Let $(p^{*},q^{*},r^{*})$ be a solution of the adjoint equation in the unknown
processes $p(t)\in\mathbb{R}^{2},\ q(t)\in\mathbb{R}^{2\times 2},\
r(t,z)\in\mathbb{R}^{2}$ satisfying the backward differential equations

$\begin{split}dp_{1}(t)&=\Bigg{[}\Big{(}-r+\frac{2rx}{\gamma}-\frac{2\omega(1+\alpha\xi)x+1}{((1+\alpha\xi)(\omega
x^{2}+1)+x)^{2}}\Big{)}p_{1}(t)-\frac{(1-\omega
x^{2})(1+(\alpha-1)\xi)}{((1+\alpha\xi)(\omega
x^{2}+1)+x)^{2}}eyp_{2}(t)-\sigma_{1}q_{1}\\\
&-\int\gamma_{1}(v)r_{1}v_{1}(dz_{1})\Bigg{]}dt+q_{1}(t)dW_{1}(t)+q_{2}(t)dW_{2}(t)+\int_{\mathbb{R}^{n}}r_{1}\widetilde{N}(dt,dz)\\\
dp_{2}(t)&=-\Bigg{[}\frac{-x}{(1+\alpha\xi)(\omega
x^{2}+1)+x}p_{1}(t)+\Big{(}\frac{e(x+\xi(\omega
x^{2}+1))}{(1+\alpha\xi)(\omega
x^{2}+1)+x}-m_{1}-2m_{2}y\Big{)}p_{2}(t)+\sigma_{2}q_{4}\\\
&+\int\gamma_{2}(v)r_{2}v_{2}(dz_{2})\Bigg{]}dt+q_{3}(t)dW_{1}(t)+q_{4}(t)dW_{2}(t)+\int_{\mathbb{R}^{n}}r_{2}\widetilde{N}(dt,dz)\\\
p_{1}(T)&=0,\ p_{2}(T)=0\end{split}$ (11)

The Hamiltonian associated with this control problem is defined as follows.

$\begin{split}H(t,x,y,\xi,p,q,r)&=1+\Bigg{[}r(1-\frac{x}{\gamma})-\frac{y}{(1+\alpha\xi)(1+\omega
x^{2})+x}\Bigg{]}xp_{1}+\Bigg{[}\frac{e(x+\xi(\omega
x^{2}+1))}{(1+\alpha\xi)(1+\omega x^{2})+x}\\\
&-m_{1}-m_{2}y\Bigg{]}yp_{2}+\sigma_{1}xq_{1}+\sigma_{2}yq_{4}+x\int\gamma_{1}r_{1}v_{1}(dz_{1})+y\int\gamma_{2}(v)r_{2}v_{2}(dz_{2})\end{split}$
(12)

Let $U=\\{\xi(t)|0\leq\xi(t)\leq\xi_{max}\forall t\in(0,t_{f}]\\}$ where
$\xi_{max}\in\mathbb{R}^{+}.$ Let $\xi^{*}\in U$ with the corresponding
solution $(x^{*},y^{*})=(x(\xi^{*}),y(\xi^{*}))$.

From the Arrow condition in the sufficient stochastic maximum Principle [29],
we have

$\begin{split}\frac{\partial H}{\partial\xi}\Big{|}_{\xi^{*}}=0&\\\
\implies&\Bigg{[}-yxp_{1}\frac{-\alpha(\omega x^{2}+1}{((1+\alpha\xi)(\omega
x^{2}+1)+x)^{2}}+\\\ &eyp_{2}\left(\frac{((1+\alpha\xi)(\omega
x^{2}+1)+x)(\omega x^{2}+1)-\alpha(\omega x^{2}+1)(x+\xi(\omega
x^{2}+1))}{((1+\alpha\xi)(\omega
x^{2}+1)+x)^{2}}\right)\Bigg{]}\Bigg{|}_{*}=0\\\ \implies&\Bigg{[}\alpha
xy(\omega x^{2}+1)p_{1}+eyp_{2}(\omega x^{2}+1)(x(1-\alpha)(\omega
x^{2}+1))\Bigg{]}\Bigg{|}_{*}=0\\\ \implies&\Bigg{[}\alpha
xp_{1}+ep_{2}(1+\omega x^{2}+x(1-\alpha))\Bigg{]}\Bigg{|}_{*}=0\\\
\implies&\alpha x^{*}p_{1}^{*}+ep_{2}^{*}(1+\omega
x^{*^{2}}+x^{*}(1-\alpha))=0\\\ \end{split}$

Hence the optimal control $\xi^{*}$ should satisfy the following condition.

$\alpha x^{*}p_{1}^{*}+ep_{2}^{*}(1+\omega x^{*^{2}}+x^{*}(1-\alpha))=0$ (13)

Since the analytical solution of (7) is complex to solve, we numerically
simulate these results in section 5.

###  4.3 Existence and Uniqueness of Solutions for the Forward Backward
Stochastic Differential Equations with Jumps (FBSDEJ)

We so far obtained the adjoint equations (7),(11) for the state equations (5)
and the objective functional (6), (10) using the sufficient stochastic maximum
principle respectively. Upon simplifying the results obtained from the arrow
condition (9), (13) from earlier two subsections, we see that the optimal
controls are given by

$\alpha^{*}=\frac{ep_{2}^{*}(1+x^{*}+\omega
x^{*^{2}})}{ep_{2}^{*}x^{*}-p_{1}^{*}x^{*}},\
\xi^{*}=\frac{x^{*}p_{1}^{*}-ex^{*}p_{2}^{*}}{ep_{2}^{*}(1+\omega x^{*^{2}})}$
(14)

In this section, we now prove the existence of optimal controls by proving the
existence of the solutions for the FBSDEJ ((5),(7),(11)) which establishes the
existence of $(x^{*},y^{*},p_{1}^{*},p_{2}^{*})$ for all simulation purposes.
Using the theorem in [30], we now prove the existence of the optimal controls
(14) in the following theorem.

######  Theorem 2.

For any $(x_{0},y_{0})\in\mathbb{R}^{+^{2}}$, the FBSDEJ ((5),(7),(11)) admits
an optimal stochastic control.

###### Proof.

Let $(X_{t})_{t\geq 0}$ be the solution of the Stochastic Differential
Equation with Jumps(SDEJ)

$\mathrm{d}X_{t}=b(X_{t})\mathrm{d}t+\sigma(X_{t})\mathrm{d}W(t)+\int_{\mathbb{R}}\Gamma(v)\widetilde{N}(\mathrm{d}t,\mathrm{d}v)$

Here the term $b(X_{t})$ denotes the drift coefficient, the term
$\sigma(X_{t})$ denotes the diffusion coefficient and the term $\Gamma(v)$
denotes the poisson term coefficient.

The theorem 1 in section 3 guaranties the monotonicity and Lipschitz
continuity of the drift coefficient, the diffusion coefficient and the poisson
term coefficient of the state equations (5).

Following the the existence and uniqueness theorem of FBDSDEJ in [30], we are
only left to prove the monotonicity and Lipschitz continuity of the drift and
diffusion terms of the adjoint system of equations (7).

From (7), due to the positivity of state variables guaranteed by theorem 1,
the drift term and the diffusion terms are given as follows.

$b(X_{t})\leq\begin{pmatrix}\left(\frac{2rx}{\gamma}\right)p_{1}(t)+\left(\omega
ex^{2}y\right)p_{2}(t)\\\ (x)p_{1}(t)+(m_{1}+2m_{2}y)p_{2}(t)\end{pmatrix},\
\sigma(X_{t})=\begin{pmatrix}q_{1}(t)&q_{2}(t)\\\
q_{3}(t)&q_{4}(t)\end{pmatrix}$

Since the drift coefficient is a linear combination of adjoint terms ($p_{1},\
p_{2}$), the monotonicity and Lipschitz continuity are guaranteed.

In addition to this, the diffusion coefficient is independent of the adjoint
terms $(p_{1},\ p_{2})$. Therefore, the monotonicity and Lipschitz continuity
are guarenteed for the diffusion coefficients.

Hence the existence of unique stochastic optimal controls are proved for
FBSDEJ ((5),(7),(11)). ∎

##  5 Numerical Simulations

In this section, we perform the extensive numerical simulations using python
by choosing the following parameters [22] for the model (5). $r=1.5,\
\gamma=12,\omega=15,\ e=0.4,m_{1}=0.15,\ m_{2}=0.01,\
\sigma_{1}=\sigma_{2}=0.02,\ \gamma_{1}=1,\gamma_{2}=1$. In these simulations,
white noise is simulated using the Box-Normal transformations and the poisson
noise is simulated using the poisson point processes [31]. The state equations
(5) and the adjoint equations (7),(11) are simulated using the Forward
Backward Doubly Stochastic Differential Equations with Jumps (FBDSDEJ) method.
The subplots in figures 1, 2 depict the optimal state trajectories, optimal
co-state trajectories, phase diagram, optimal quality of additional food and
the optimal quantity of additional food respectively.

####  A. Applications to Biological Conservation:

The subplot (1(a)) depicts the optimal state trajectory of the system (5) from
the initial state $(2,8)$ that stabilizes over time around the state
$(16,90)$. The subplot (1(c)) gives the phase diagram which shows the
trajectories are stabilized over high values of prey and predator. The
subplots (1(d)) and (1(e)) depicts the optimal quality and quantity of
additional food respectively. These plots show that the high quality of
additional food is required to achieve biological conservation. Even if the
quantity of additional food is lower, still we will be able to achieve
biological conservation with higher quality of additional food.

(a)

(b)

(c)

(d)

(e)

Figure 1: This ﬁgure depicts the simulations of time-optimal control problem
with respect to the control parameters in the context of biological
conservation.

####  B. Applications to Pest Management:

The subplot (2(a)) depicts the optimal state trajectory of the system (5) from
the initial state $(2,8)$ that reaches the nearly prey-elimination stage
around the state $(5,90)$. The subplot (2(c)) depicts this property more
clearly through the phase diagram where it reaches the lowest prey value over
the time. The subplots (2(d)) and (2(e)) depicts that a lesser quality of
additional food and a lower quantity of additional food is good enough to
achieve pest management where pest is viewed as prey.

(a)

(b)

(c)

(d)

(e)

Figure 2: This ﬁgure depicts the simulations of time-optimal control problem
with respect to the control parameters in the context of pest management.

##  6 Conclusion

This paper is about a stochastic prey-predator system exhibiting Holling type
IV functional response along with the combined white noise and Lévy noise. We
do the time-optimal control studies for the system, where the quality of
additional food and the quantity of additional food are considered as control
parameters respectively. To begin with, we extend the model by incorporating
multiplicative noise to both prey and predator. In theorem 1, we proved the
existence of a unique positive global solution of (5). Further, we formulated
the time-optimal control problem with the objective to minimize the final time
in which the system reaches the pre-defined state. Using the sufficient
stochastic maximum principle, we characterized the optimal control values. In
theorem 2, we proved that the existence and uniqueness of Forward Backward
Doubly Stochastic Differential Equations with Jumps (FBDSDEJ). Lastly, we
numerically simulated the theoretical findings and applied them in the context
of biological conservation and pest management.

Some of the salient features of this work include the following. Unlike the
most traditional papers, here we considered a stochastic time-optimal control
problem. As Intra-specific competition among predators is ineluctable, we also
explicitly incorporated the intra-specific competition into our model. This
paper mainly deals with the novel study of the time-optimal control problems
where the state equations involve both the discrete and continuous noise which
is challenging. In future research, we hope to add and study more
specificalities such as mutual interference, stochastic bifurcations, Markov
chain and partial differential systems.

### Financial Support:

This research was supported by National Board of Higher Mathematics(NBHM),
Government of India(GoI) under project grant - Time Optimal Control and
Bifurcation Analysis of Coupled Nonlinear Dynamical Systems with Applications
to Pest Management,
Sanction number: (02011/11/2021NBHM(R.P)/R$\&$D II/10074).

### Conflict of Interests Statement:

The authors have no conflicts of interest to disclose.

### Ethics Statement:

This research did not required ethical approval.

### Acknowledgments

The authors dedicate this paper to the founder chancellor of SSSIHL, Bhagawan
Sri Sathya Sai Baba. The corresponding author also dedicates this paper to his
loving elder brother D. A. C. Prakash who still lives in his heart.

## References

  * [1] Alfred James Lotka  “Elements of physical biology”  Williams & Wilkins, 1925
  * [2] Vito Volterra  “Variazioni e fluttuazioni del numero d’individui in specie animali conviventi”  C. Ferrari Venezia, 1927
  * [3] Mark Kot  “Elements of mathematical ecology”  Cambridge University Press, 2001
  * [4] Johan A Metz and Odo Diekmann  “The dynamics of physiologically structured populations”  Springer, 2014
  * [5] Stefan Geritz and Mats Gyllenberg  “A mechanistic derivation of the DeAngelis–Beddington functional response”  In _Journal of theoretical biology_ 314 Elsevier, 2012, pp. 106–108
  * [6] Roger Arditi and Lev R Ginzburg  “Coupling in predator-prey dynamics: ratio-dependence”  In _Journal of theoretical biology_ 139.3 Elsevier, 1989, pp. 311–326
  * [7] MP Hassell and GC Varley  “New inductive population model for insect parasites and its bearing on biological control”  In _Nature_ 223.5211 Nature Publishing Group, 1969, pp. 1133–1137
  * [8] Philip H Crowley and Elizabeth K Martin  “Functional responses and interference within and between year classes of a dragonfly population”  In _Journal of the North American Benthological Society_ 8.3 North American Benthological Society, 1989, pp. 211–221
  * [9] Shuichi Yano  “Cooperative web sharing against predators promotes group living in spider mites”  In _Behavioral Ecology and Sociobiology_ 66.6 Springer, 2012, pp. 845–853
  * [10] Melanie McClure and Emma Despland  “Defensive responses by a social caterpillar are tailored to different predators and change with larval instar and group size”  In _Naturwissenschaften_ 98.5 Springer, 2011, pp. 425–434
  * [11] P D N Srinivasu, B S R V Prasad and M Venkatesulu  “Biological control through provision of additional food to predators: a theoretical study”  In _Theoretical Population Biology_ 72.1 Elsevier, 2007, pp. 111–120
  * [12] P D N Srinivasu and B S R V Prasad  “Time optimal control of an additional food provided predator–prey system with applications to pest management and biological conservation”  In _Journal of mathematical biology_ 60.4 Springer, 2010, pp. 591–613
  * [13] P D N Srinivasu, D K K Vamsi and I Aditya  “Biological conservation of living systems by providing additional food supplements in the presence of inhibitory effect: a theoretical study using predator–prey models”  In _Differential Equations and Dynamical Systems_ 26.1 Springer, 2018, pp. 213–246
  * [14] P D N Srinivasu, D K K Vamsi and V S Ananth  “Additional food supplements as a tool for biological conservation of predator-prey systems involving type III functional response: A qualitative and quantitative investigation”  In _Journal of theoretical biology_ 455 Elsevier, 2018, pp. 303–318
  * [15] D K K Vamsi, Deva Siva Sai Murari Kanumoori and Bishal Chhetri  “Additional food supplements as a tool for biological conservation of biosystems in the presence of inhibitory effect of the prey”  In _Acta biotheoretica_ 68.3 Springer, 2020, pp. 321–355
  * [16] V S Ananth and D K K Vamsi  “An Optimal Control Study with Quantity of Additional food as Control in Prey-Predator Systems involving Inhibitory Effect”  In _Computational and Mathematical Biophysics_ 9.1 De Gruyter Open Access, 2021, pp. 114–145
  * [17] V S Ananth and D K K Vamsi  “Achieving Minimum-Time Biological Conservation and Pest Management for Additional Food provided Predator–Prey Systems involving Inhibitory Effect: A Qualitative Investigation”  In _Acta Biotheoretica_ 70.1 Springer, 2022, pp. 1–51
  * [18] Robert M May  “Stability and complexity in model ecosystems”  In _Stability and Complexity in Model Ecosystems_ Princeton university press, 2019
  * [19] Lin Li and Wencai Zhao  “Deterministic and stochastic dynamics of a modified Leslie-Gower prey-predator system with simplified Holling-type IV scheme”  In _Math. Biosci. Eng_ 18, 2021, pp. 2813–2831
  * [20] Dongsheng Xu, Ming Liu and Xiaofeng Xu  “Analysis of a stochastic predator–prey system with modified Leslie–Gower and Holling-type IV schemes”  In _Physica A: Statistical Mechanics and its Applications_ 537 Elsevier, 2020, pp. 122761
  * [21] Wantao Jia, Yong Xu, Dongxi Li and Rongchun Hu  “Stochastic Analysis of Predator–Prey Models under Combined Gaussian and Poisson White Noise via Stochastic Averaging Method”  In _Entropy_ 23.9 MDPI, 2021, pp. 1208
  * [22] Tingting Ma, Xinzhu Meng and Zhengbo Chang  “Dynamics and optimal harvesting control for a stochastic one-predator-two-prey time delay system with jumps”  In _Complexity_ 2019 Hindawi, 2019
  * [23] Sampurna Sengupta, Pritha Das and Debasis Mukherjee  “Stochastic non-autonomous Holling type-III prey-predator model with predator’s intra-specific competition”  In _Discrete & Continuous Dynamical Systems-B_ 23.8 American Institute of Mathematical Sciences, 2018, pp. 3275
  * [24] Erin N Bodine and Anne E Yust  “Predator-prey dynamics with intraspecific competition and an allee effect in the predator population”  In _Letters in Biomathematics_ 4.1, 2017, pp. 23–38
  * [25] Shyam Pada Bera, Alakes Maiti and Guruprasad Samanta  “Stochastic analysis of a prey–predator model with herd behaviour of prey”  In _Nonlinear Analysis: Modelling and Control_ 21.3, 2016, pp. 345–361
  * [26] Mustapha Belabbas, Abdelghani Ouahab and Fethi Souna  “Rich dynamics in a stochastic predator-prey model with protection zone for the prey and multiplicative noise applied on both species”  In _Nonlinear Dynamics_ 106.3 Springer, 2021, pp. 2761–2780
  * [27] A La Cognata, D Valenti, AA Dubkov and B Spagnolo  “Dynamics of two competing species in the presence of Lévy noise sources”  In _Physical Review E_ 82.1 APS, 2010, pp. 011121
  * [28] Bernt Øksendal and Agnes Sulem  “Stochastic Control of jump diffusions”  Springer, 2005
  * [29] Nils C Framstad, Bernt Øksendal and Agnes Sulem  “Sufficient stochastic maximum principle for the optimal control of jump diffusions and applications to finance”  In _Journal of optimization theory and applications_ 121.1 Springer, 2004, pp. 77–98
  * [30] AbdulRahman Al-Hussein and Boulakhras Gherbal  “Existence and uniqueness of the solutions of forward-backward doubly stochastic differential equations with Poisson jumps”  In _Random Operators and Stochastic Equations_ 28.4 De Gruyter, 2020, pp. 253–268
  * [31] Xiaoling Zou and Ke Wang  “Numerical simulations and modeling for stochastic biological systems with jumps”  In _Communications in Nonlinear Science and Numerical Simulation_ 19.5 Elsevier, 2014, pp. 1557–1568