# Mathematical modeling the wildfire propagation in a Randers space

Hengameh R. Dehkordi
Center of Mathematics, Computation and Cognition,
Federal University of ABC,
Santo André, Brazil ORCID ID: 0000-0002-1738-3373,
<EMAIL_ADDRESS><EMAIL_ADDRESS>Phone: +55 (11)
4996-8332

###### Abstract

The devastating effects of wildfires on the wildlife and their impact on human
lives and properties are undeniable. This shows the importance of studying the
spread of wildfire, predicting its behavior and presenting more reliable
models for its propagation. Here, by using the validity of the Huygens’
envelope principle for wavefronts in Randers spaces, we present some models
for the propagation of wildfire in an $n$-dimensional smooth manifold under
the presence of wind. In the models, trajectories of fire particles are
tracked and the equations that give the wavefront at each time are provided.
Furthermore, we determine the paths and points of great importance in the
process of wildfire management, called strategic paths and points. Finally, we
consider two examples of spreading the wildfire in some agricultural land or
woodland, for the sake of illustration.

MSC 2010: 53B40; 53B50; 83C57; 83C80

Keywords: Randers metric; Huygens’ Principle; wavefronts; wave rays; causal
structure; analogue gravity.

##  1 Introduction

Every year, wildfires cause significant damage to the wildlife, jungles,
grasslands, agricultural lands, and natural resources; and threaten
infrastructures, properties, and human lives [10]. Sometimes, the renovation
of destroyed regions and the recovery of damaged faunas are impossible.
Indeed, wildfires have strongly negative ecological effects and every year
engulf millions hectares of rainforests [11]. Global warming and carbon
dioxide released into the atmosphere due to wildfires are other issues that
can not be ignored [27].

Providing a more accurate and reliable model for spreading the wildfire in
time plays an important role in the wildfire management strategies and without
any doubt reduces both the financial and life losses due to wildfire. Finsler
geometry, which is a classical branch of differential geometry started by P.
Finsler in 1918, is a strong tool to model some real phenomena in anisotropic
or inhomogeneous media [12, 17, 28]. On the propagation of waves and tracking
the wave rays, several authors have already applied the Finsler metric [13,
19]. As a particular case, the propagation of wildfire waves in dimension $2$
was studied in [19], in which the author showed that, for a wildfire
spreading, the wavefronts and wave rays are respectively the geometric spheres
and geodesics of the corresponding Finsler metric. Very recently, the validity
of Huygens’ principle was verified for wavefronts in Finsler spaces of any
dimension $n$ [9].

The Huygens’ principle is frequently used for modeling the growth of wildfire
in dimension $2$. A literature review shows that several authors used some
fixed elliptical template fields, such as double ellipse, lemniskata, oval
shape, and tear shape, as spherical wavefronts in the process of applying the
Huygens’ principle [1, 7, 14, 19, 22]. However, the spherical wavefronts used
in this process have some deviations from the spherical wavefront in reality
[19]. Consequently, sometimes, the model presented based on these spherical
wavefronts is neither accurate enough nor reliable. The reason why these
template fields are not sometimes so suitable is that the curvature of space
is assumed to be zero while, in reality, it is not zero. Therefore, the more
the curvature differs from zero, the less the model approximates the
propagation [15]. By the way, in the models presented by using the Finsler
metric, the curvature of space is taken into account. In fact, the Finsler
metric, by taking the curvature into account, provides geometric spheres with
negligible deviation from the spherical wavefronts created by the wildfire.

The use of simulators, such as Phoenix, IGNITE, Bushfire, FireMaster, FARSITE
and Prometheus is another technique widely used by the researchers to predict
the behavior of the fire [16, 21, 26]. To the best of our knowledge, FARSITE
and Prometheus are considered to be the best ones among these simulators [25].
The challenge we confront dealing with simulators is that they produce errors
during the process. Therefore, regarding simulators, new problems appear which
are reducing the errors and time of computations [2].

In [9], it was shown that for a given propagation of waves, if one finds a
Finsler metric $F$ in such a way that one of the preimages of Finsler distance
function coincides with the wavefront at some time $t$, then the Finsler
distance function provides a model for that propagation. In fact, the Finsler
geodesics and distance function are applied to modeling the propagation of
waves. In our work, by using the results of [9], we suggest two different
strategies to predict the progress of fire waves in the case of Randers
metrics. In other words, in order to provide the model of a given propagation,
one may use the Randers geometric spheres, and then apply the Huygens’
principle, or solve the system of Randers geodesic equations, which are indeed
paths of fire particles. By the way, by the spread of wildfire in a Randers
space, it means a fire spreading across a smooth manifold $M$ under the
presence of some wind $W$. The wind here is a smooth vector field $W$ such
that $|W|<1$.

In this work, through three main results, that is Theorems 3.2, 3.5 and 3.6,
we provide the equations of wave rays and wavefronts at each time $\tau$.
Moreover, the equations of strategic paths and points are presented. Let us
explain that by strategic paths we mean the paths along which the fire engulfs
more regions or it reaches to some special zone that has some priority
(protecting fauna, houses, etc.) and should be protected from the fire. In
other words, firefighters and equipment should be located along the strategic
paths in order to control the fire or prevent it from progressing toward some
special direction or area. In fact, having strategic paths increases the
chance of success in wildfire management. After finding the model of
propagation and strategic paths, by using some information on the behavior of
wildfires, one determines some points where the firefighters and equipment
should be located to attack the fire. Such points, that are along strategic
paths, are called strategic points. It should be pointed out that time is so
important when it comes to responding to an emergency incident and also our
sources, forces, and equipment are limited against natural phenomena such as
wildfires. Therefore, finding the strategic points, which leads to saving time
and expense, is vital to wildfire management, especially in wildfires of big
scale and magnitude that are happening every year around the world, such as
the United States, Australia, and Brazil. Hence, it is fair to claim that the
study of strategic points and paths is of great importance in the firefighting
process and really demands more attention. Whereas, to the best of our
knowledge, there is no study related to such paths or points.

###  1.1 Hypotheses, methodology, and Outline of the paper

Throughout this work, it is assumed that a wildfire is sweeping across some
space $M$ which is a smooth manifold of dimension $n$ and some fuel has been
distributed uniformly and smoothly through $M$ Also the temperature and
moisture are constant at all points of $M$. A mild wind $W$, that is a smooth
vector field, is blowing in $M$. Although the wind might or not be space-
dependent, it must be time-independent at intervals of time. Also, it is
assumed that the fire is stopped before it creates singularities or cut loci.
To be closer to reality, the focus of our work is on the dimension $3$ (for
instance $M$ can be any open subset of $\mathbb{R}^{3}$), however our results
are valid for any dimension $n$.

On the methodology, we show that the model of propagation is merely determined
by considering some translated ellipsoids. This way, it is enough to determine
the equation of such ellipsoids from the experimental or laboratory data.
Afterwards, equations of wavefronts and wave rays associated with the
propagation are determined, and finally the model is presented. Following this
strategy, we investigate the spread of wildfire across some space under the
presence of wind which is a constant, Killing or smooth vector field.

The rest of this paper is organized as follows. In Section 2, we recall some
preliminaries on Randers spaces and Hyugens’ principle. In Section 3, we
proceed with the main results and certain discussions on the spread of fire
waves. In Section 4, we provide some examples in which certain wildfire is
spreading in some agricultural land, woodland or forest under the influence of
wind. In Section 5, concluding remarks and some ideas for future works are
provided.

##  2 preliminaries

Here, for the sake of self-contained paper, we first provide a brief review of
the Randers geometry, and then recall the Huygens’ envelope principle.

Let $M$ be a smooth manifold, $p=(x_{1},...,x_{n})\in M$ a point and $T_{p}M$
the space tangent to $M$ at $p$. Assume that $V=(v_{1},...,v_{n})\in T_{p}M$
is a vector according to the canonical basis $\\{\frac{\partial}{\partial
x_{i}}\\}_{i=1}^{n}$ for $T_{p}M$ and $TM$ the tangent bundle, i.e. the
collection of all vectors tangent to $M$ or in other words

$TM=\underset{p\in M}{\dot{\cup}}\\{(p,V):V\in T_{p}M\\}.$

A Riemannian metric on $M$ is a smooth function $h$ such that to each point
$p\in M$, $h$ assigns a positive-definite inner product $h_{p}:T_{p}M\times
T_{p}M\to\mathbb{R}$. The smoothness condition on $h$ refers to the fact that
the function $p\in M\to h_{p}(\frac{\partial}{\partial
x_{i}},\frac{\partial}{\partial x_{j}})\in\mathbb{R}$ must be smooth. Now,
consider a Riemannian metric $\mathfrak{r}$ on $M$, $\mathfrak{r}:TM\times
TM\to[0,\infty)$, and a $1$-form $\beta:TM\to[0,\infty)$ such that
$\mathfrak{r}(\beta^{*},\beta^{*})<1$, where $\beta^{*}$ stands for the dual
vector of $\beta$. Considering $\alpha(.)=\sqrt{\mathfrak{r}(.,.)}$, then
$F=\alpha+\beta$ is called the Randers metric. It is well known that given any
Randers metric $F$ on $M$, there is a Riemannian metric $h$ on $M$ and some
vector field $W$, which satisfies $h(W,W)<1$, related to this Randers metric
$F$. We call the pair $(W;h)$ the Zermelo data associated with $F$. Due to
Zermelo, the problem which associates the Randers metric $F$ with the
Riemannian metric $h$ is called the Zermelo’s navigation problem. For a
literature review on the Zermelo’s problem of navigation see [20]. The Randers
metric $F$ and the associated Zermelo data are satisfied in the following
equation:

$F(V)=\alpha(V)+\beta(V)=\frac{\sqrt{h^{2}(W,V)+\lambda
h(V,V)}}{\lambda}-\frac{h(W,V)}{\lambda},$ (2.1)

where $\lambda=1-h(W,W)$ [5]. As a result of Eq. (2.1), one has

$F(V)=1\ \text{if and only if}\ h(V-W,V-W)=1$ (2.2)

(see the details in section $1.1.2$ of [5]).

Given a Randers space $(M,F)$ and some point $p\in M$, the indicatrix of
radius $\tau$ at $p$ is the set

$\mathcal{I}_{{}_{F}}^{\tau}=\\{V\in T_{p}M:F(V)=\tau\\}.$

Indeed, the indicatrix $\mathcal{I}_{{}_{F}}^{\tau}$ is the collection of all
end points of vectors tangent to $M$ at $p$ and length $\tau$, and therefore
it is a hypersurface in $T_{p}M$. We show the indicatrix of radius $1$ with
$\mathcal{I}_{{}_{F}}$. By the way, given the equation of Randers indicatrix,
there are some techniques to determine the equation of its Randers metric (for
instance, see [4], p. 13). Given any Riemannian manifold $(M,h)$, one defines
the Riemannian indicatrix of radius $\tau$, $\mathcal{I}_{h}^{\tau}$ similar
to the Randers indicatrix $\mathcal{I}_{{}_{F}}^{\tau}$.

Given some Randers space $(M,F)$ and a piece-wise smooth curve
${\gamma}:[0,1]\longrightarrow M$, the length of ${\gamma}$ is defined as
$L[{\gamma}]:=\int_{0}^{1}F({\gamma^{\prime}}(t))dt$. Similar to the
Riemannian space, the distance from a point $p\in M$ to another point $q\in M$
in the Randers space $(M,F)$ is

$d(p,q):=\inf_{{\gamma}}\int_{0}^{1}F({\gamma^{\prime}}(t))dt,$ (2.3)

where the infimum is taken over all piece-wise smooth curves
${\gamma}:[0,1]\longrightarrow M$ joining $p$ to $q$. A smooth curve in a
Randers manifold is called a geodesic (shortly $F$-geodesic) if it is locally
the shortest time path connecting any two nearby points on this curve. Given a
compact subset $A\subset M$, we define the Randers distance function
$\rho:M\longrightarrow\mathbb{R}$ with $\rho(p)=d(A,p)$. One proves that
$\rho$ is locally Lipschitz continuous and therefore it is differentiable
almost everywhere [24].

Given some smooth vector field $W$ on $M$, the flow of $W$ is the smooth map
$\varphi:(-\epsilon,\epsilon)\times M\longrightarrow M$ such that for all
$p\in M$, $\varphi^{p}(t):=\varphi(t,p)$ is an integral curve of $W$, that is
$\frac{d\varphi^{p}}{dt}(t)|_{t=0}=W(p)$ [18]. A vector field $W$ on a
Riemannian space is called Killing if and only if its flow is an isometry of
$(M,h)$. One can also say that $W$ is Killing if and only if
$\mathcal{L}_{W}h=0$, where $\mathcal{L}$ is the Lie derivative.

######  Lemma 2.1.

[23] Assume that $(M,h)$ is a Riemannian manifold. Given a unitary Riemannian
geodesic ($h$-geodesic) ${\gamma_{{}_{h}}}:(-\epsilon,\epsilon)\longrightarrow
M$ and a Killing vector field $W$, the unitary $F$-geodesics are
${\gamma_{{}_{F}}}(t)=\varphi(t,{\gamma_{{}_{h}}}(t))$, where
$\varphi:(-\epsilon,\epsilon)\times U\longrightarrow M$ is the flow of $W$
through ${\gamma_{{}_{h}}}(t)$.

Given some Randers space $(M,F)$ and a submanifold $A\subset M$, we say that a
vector $V$ is  $F$-orthogonal to $A$, and write $V\underset{F}{\perp}A$, if
for every vector $U$ tangent to $A$ one has $g_{V}(V,U)=0$, where

$g_{V}(V,U):=\frac{1}{2}\bigg{(}\frac{\partial^{2}}{\partial t\partial
s}F^{2}(V+tV+sU)\bigg{)}_{s=t=0}.$

Similarly, given the Riemannian manifold $(M,h)$, a vector $V$ is
$h$-orthogonal to $A$, i.e. $V\underset{h}{\perp}A$, if for every vector $U$
tangent to $A$ one has $h(V,U)=0$. The following corollary states the relation
between $F$-orthogonality and $h$-orthogonality.

######  Corollary 2.2.

[8] Given a Randers manifold $(M,F)$, assume that $(h;W)$ is the Zermelo data
associated with it. Then, for any two non-zero vectors $U$ and $V$ tangent to
$M$, $U\underset{F}{\perp}V$ if and only if $h(U,\frac{V}{F(V)}-W)=0$.

###  2.1 Huygens’ principle

Assume that $P$ is a source which emits waves. Given any time $t>0$, we
consider the collection of all points to which the wave reaches at time $t$.
This collection is called the wavefront at time $t$ [3]. In the case that $P$
is a single point, the wavefront is called the spherical wavefront. Given a
wavefront $B$, assume that each point on $B$ acts as a new source that emits
spherical wavefronts. At any time later, a surface tangent to all of the
spherical wavefronts is called the envelope of $B$. By the wave ray it means
the shortest time path that connects any point of $B$ to the wavefront at any
time later. We recall the Huygens’ Theorem as follows:

######  Theorem 2.3.

[3] Let $\phi_{p}(t)$ be the wavefront of the point $p$ after time $t$. For
every point $q$ of this wavefront, consider the wavefront after time $s$, i.e.
$\phi_{q}(s)$. Then, the wavefront of point $p$ after time $s+t$,
$\phi_{p}(s+t)$, will be the envelope of wavefronts $\phi_{q}(s)$, for
$q\in\phi_{p}(t)$.

We recall the following result from [9] which is used throughout this work
several times. By the way, since a Randers metric is a particular case of a
Finsler metric, any result on Finsler is valid for Randers, as well. It is
assumed that, in a Finsler space $(M,F)$, a wildfire is spreading and sweeping
some area $U\subset M$ in the interval of time from $t=s>0$ to $t=r$. It is
also assumed that $U$ is a smooth manifold and $d$ is the Finsler distance
function.

######  Proposition 2.4.

[9] Let $\rho:M\longrightarrow\mathbb{R}$ with $\rho(p)=d(A,p)$ where $A$ is a
compact subset of $M$ and $\rho(U)=[s,r]$, where $0<s<r$. Suppose that
$\rho^{-1}(s)$ is the wavefront at time $0$ and there is no cut loci in
$\rho^{-1}(s,r)$. Then, for each $t\in[s,r]$, $\rho^{-1}(t)$ is the wavefront
at time $t-s$ and the Huygens’ principle is satisfied by the wavefronts

$\\{\rho^{-1}(t)\\}_{t\in[s,r]}.$

Furthermore, the track of each fire particle is a geodesic of $F$ and also it
is orthogonal to each wavefront $\rho^{-1}(t)$ at time $t-s$.

This theorem says that once one finds the Finsler metric associated with a
wildfire spreading in a smooth manifold $M$, by using the distance function
$\rho$, the model of propagation is provided. As a result from Proposition
2.4, one has the following corollary.

######  Corollary 2.5.

Let $\rho:M\longrightarrow\mathbb{R}$ with $\rho(p)=d(A,p)$ where $A$ is a
compact subset of $M$ and $\rho(U)=[0,r]$. Suppose that there is no cut loci
in $\rho^{-1}(0,r)$ and $\rho^{-1}(0)$ is the wavefront at time $0$. Then, for
each $t\in[0,r]$, $\rho^{-1}(t)$ is the wavefront at time $t$ and the Huygens’
principle is satisfied by the wavefronts

$\\{\rho^{-1}(t)\\}_{t\in[0,r]}.$

Furthermore, the track of each fire particle is a geodesic of $F$ and also it
is orthogonal to each wavefront $\rho^{-1}(t)$ at time $t$.

###### Proof.

In fact, since $\rho$ is differentiable almost everywhere It is enough to take
the limit of the function $\rho$ to extend the results of Proposition 2.4 to
the whole interval $[0,r]$. ∎

##  3 Modeling the propagation by using the Randers geometry

Throughout this section, we consider the hypothesis stated in the Introduction
and establish some results. We consider three different states for the wind:
constant, Killing, and smooth vector filed.

Here, it is shown that for modeling a propagation one first needs to find the
equation of some ellipsoid - as several authors have used elliptic fields in
the case of dimension $2$ [1]. From this ellipsoid, that is actually our
indicatrix, we calculate the Riemannian metric, then the wave rays, and
finally we provide the model. Before presenting the main results, we state and
prove Lemma 3.1 which says that the Randers indicatrix of radius $\tau$ is
translation of the Riemannian indicatrix of the same radius by vector $\tau
W$.

######  Lemma 3.1.

Given a Randers space $(M,F)$, let $(h;W)$ be the Zermelo data associated with
$F$ and $\mathcal{I}_{{}_{F}}^{\tau}$ and $\mathcal{I}_{{}_{h}}^{\tau}$ the
Randers and Riemannian indicatrices of radius $\tau$, respectively. Then,
$\mathcal{I}_{{}_{F}}^{\tau}=\mathcal{I}_{{}_{h}}^{\tau}+\tau W.$

###### Proof.

One has $V\in\mathcal{I}_{{}_{h}}^{\tau}$ if and only if $|V|_{{}_{h}}=\tau$
if and only if $|\frac{V}{\tau}|_{{}_{h}}=1$ if and only if
$\frac{V}{\tau}\in\mathcal{I}_{{}_{h}}$, where
$|\,.\,|_{{}_{h}}=\sqrt{h(\,.,.\,)}$. In other words,
$\mathcal{I}_{{}_{h}}^{\tau}=\tau\mathcal{I}_{{}_{h}}$. Similarly, one shows
$\mathcal{I}_{{}_{F}}^{\tau}=\tau\mathcal{I}_{{}_{F}}$. Furthermore, according
to the relation 2.2, $F(V)=1$ if and only if $|V-W|_{{}_{h}}=1$, or in other
words $V\in\mathcal{I}_{{}_{F}}$ if and only if $V-W\in\mathcal{I}_{{}_{h}}$
which leads to $\mathcal{I}_{{}_{F}}=\mathcal{I}_{{}_{h}}+W$. Finally, from
these relations we have

$\displaystyle\mathcal{I}_{{}_{F}}^{\tau}$
$\displaystyle=\tau\mathcal{I}_{{}_{F}}=\tau(\mathcal{I}_{{}_{h}}+W)=\tau\mathcal{I}_{{}_{h}}+\tau
W$ $\displaystyle=\mathcal{I}_{{}_{h}}^{\tau}+\tau W,$

which is the desired relation. ∎

From this lemma, given a Randers indicatrix $\mathcal{I}_{F}$, one can find
the Riemannian indicatrix $\mathcal{I}_{h}$, and vice versa. Consequently,
assuming that $Q_{h}$ is the quadratic equation of $\mathcal{I}_{h}$, the
Riemannian metric is $h=\frac{1}{2}\operatorname{Hess}Q_{h}$ and, by using Eq.
(2.1), one finds the Randers metric $F$.

###  3.1 Constant wind

In this section, for a constant wind $W$, we first find the wavefronts and
wave rays of the propagation, and then the strategic paths and points.

####  3.1.1 Wavefronts and Wave rays

######  Theorem 3.2.

Assume that a wildfire is spreading in some space $M$ while the wind
$W=(0,W_{2},W_{3})$ is blowing across $M$ and $A$ is the wavefront at time
$0$. Then:

  * (i)

Given any point $p$ in $M$, the spherical wavefront of radius $\tau$ and
center $p$ is

$\displaystyle Q_{{}_{h}}(\frac{u}{\tau},\frac{v}{\tau},\frac{w}{\tau})+\tau
W+p,$ (3.1)

where $Q_{h}(u,v,w)$ is given by the Eq. (3.4).

  * (ii)

Given any point $p$ in $A$, the wave ray emanating from $p$ is the straight
line given by ${\gamma_{{}_{F}}}(t)=p+tV$, where $V$ is a vector such that
$|V-W|_{{}_{h}}=1$ and $V-W\underset{h}{\perp}A$, where $h=\verb*|Hess|\,Q$.

  * (iii)

The wavefront at time $\tau$ is the set

$\\{p+\tau V\ :\ p\in A,|V-W|_{{}_{h}}=1,V-W\underset{h}{\perp}A\\}.$

###### Proof.

$(i):$ First of all, observe that we have the same conditions at all points of
$M$. Hence, the spherical wavefronts of the radius $\tau$ and centers at any
points of $M$ are the same geometric objects, up to a translation.
Consequently, the metric $F$ associated with the spherical spheres does not
depend on the point, however it depends on the direction. In other words, $F$
is a Minkowski-Randers type. As a consequence, the Randers geodesics are
straight lines and the space is of zero curvature. Therefore, the spherical
wavefront of radius $\tau$ and center $p$ coincides with
$\mathcal{I}^{\tau}_{{}_{F}}$, that is the indicatrix of radius $\tau$ and
center $0_{p}\in T_{p}M$. Let us find the equation of
$\mathcal{I}^{\tau}_{{}_{F}}$. Assume that $h$ is the metric in Zermelo data
$(h;W)$ whose Zermelo’s solution is $F$. By Lemma 3.1, $\mathcal{I}_{h}$ is a
translation of $\mathcal{I}_{F}$ by $-W$. Therefore, let’s find
$\mathcal{I}_{h}$, in order to find $\mathcal{I}_{{}_{F}}$. Assume that
$Q_{{}_{h}}$ is the quadratic equation of $\mathcal{I}_{{}_{h}}$, that is

$Q_{{}_{h}}(u,v,w)=[u,v,w]\hbar[u,v,w]^{T}=1,$

where $\hbar$ is the matrix of the metric $h$ and $B^{T}$ is the transpose of
the matrix $B$. Observe that $\hbar$ is a symmetric and positive-definite
matrix, since it is the matrix of the metric $h$. Also, the fact that
$Q_{{}_{h}}$ at each point has the same equation implies that $\hbar$ has
constant components. We show through the following lemma that $Q_{{}_{h}}$ is
a rotated ellipsoid.

######  Lemma 3.3.

Let $Q_{{}_{h}}(u,v,w)$ be a quadratic equation such that

$Q_{{}_{h}}(u,v,w)=[u,v,w]\hbar[u,v,w]^{t}=1,$

where $\hbar$ is a symmetric and positive-definite matrix with constant
components. Then, $Q_{{}_{h}}$ is the equation of a rotated ellipsoid.

###### Proof.

Since $\hbar$ is a symmetric and positive-definite matrix, by Spectral
theorem, there exists some orthogonal matrix $P$ (i.e. $PP^{T}=P^{T}P=I$) such
that $\hbar=P^{T}\mathcal{D}P$, where
$\mathcal{D}=diag(\lambda_{1},\lambda_{2},\lambda_{3})$ is a diagonal matrix.
Therefore, one can write

$\begin{array}[]{ll}1=\begin{bmatrix}u\\\ v\\\
w\end{bmatrix}^{T}P^{T}\mathcal{D}P\begin{bmatrix}u\\\ v\\\
w\end{bmatrix}&=(P\begin{bmatrix}u\\\ v\\\ w\\\
\end{bmatrix})^{T}\mathcal{D}P\begin{bmatrix}u\\\ v\\\ w\\\ \end{bmatrix}\\\
&\\\ &=\begin{bmatrix}u^{\prime}\\\ v^{\prime}\\\ w^{\prime}\\\
\end{bmatrix}^{T}\mathcal{D}\begin{bmatrix}u^{\prime}\\\ v^{\prime}\\\
w^{\prime}\\\
\end{bmatrix}=\lambda_{1}{u^{\prime}}^{2}+\lambda_{2}{v^{\prime}}^{2}+\lambda_{3}w{{}^{\prime}}^{2},\end{array}$

where we put $P\begin{bmatrix}u\\\ v\\\ w\\\
\end{bmatrix}=\begin{bmatrix}u^{\prime}\\\ v^{\prime}\\\ w^{\prime}\\\
\end{bmatrix}$ and the last equality is because $\mathcal{D}$ is a diagonal
matrix. Furthermore, since $\hbar$ is positive definite, the elements on the
principle diagonal of $\mathcal{D}$ are positive real numbers and we can write

$1=\lambda_{1}{u^{\prime}}^{2}+\lambda_{2}{v^{\prime}}^{2}+\lambda_{3}w{{}^{\prime}}^{2}=(\frac{u^{\prime}}{a})^{2}+(\frac{v^{\prime}}{b})^{2}+(\frac{w^{\prime}}{c})^{2},$
(3.2)

where $a=\frac{1}{\sqrt{\lambda_{1}}}$, $b=\frac{1}{\sqrt{\lambda_{2}}}$,
$c=\frac{1}{\sqrt{\lambda_{3}}}$. Eq. (3.2) is the equation of an ellipsoid in
the vector space $T_{p}M$ with basis $\mathbb{B}=\\{P\frac{\partial}{\partial
x},P\frac{\partial}{\partial y},P\frac{\partial}{\partial z}\\}$, which is the
rotation of the canonical basis. Therefore, it is deduced that $Q_{h}$ is some
rotated ellipsoid in the system with basis of $\\{\frac{\partial}{\partial
x},\frac{\partial}{\partial y},\frac{\partial}{\partial z}\\}$. ∎

From Lemma 3.3, the quadratic equation of spherical wavefronts with respect to
the Riemannian metric $h$, $Q_{h}(u,v,w)$, is some rotated ellipsoid given by
Eq. (3.2). Let us figure out what is the angle of rotation. Observe that,
according to the following facts: the heat goes toward above, i.e. $z$-axis,
the fuel has been distributed uniformly through $M$, and $W$ is in the
$yz$-plane; it is deduced that $Q_{h}$ is rotated around $x$-axis. That is the
matrix of rotation is

$\displaystyle P=R_{x}(\alpha)=\left(\begin{matrix}1&0&0\\\
0&\cos\alpha&-\sin\alpha\\\ 0&\sin\alpha&\cos\alpha\end{matrix}\right).$ (3.3)

Here, we used a right-handed coordinate system and a right-handed rotation
through an angle $\alpha$ around $x$-axis. To summarize we have:

$\displaystyle
Q_{{}_{h}}(u,v,w)=(\frac{u}{a})^{2}+(\frac{v\cos\alpha-w\sin\alpha}{b})^{2}+(\frac{v\sin\alpha+w\cos\alpha}{c})^{2}=1,$
(3.4)

were $a,b$, $c$ and $\alpha$ are constant real numbers and will be determined
from the experimental data. Finally, the Riemannian indicatrix is

$\mathcal{I}_{{}_{h}}=\\{(u,v,w)\ :Q_{{}_{h}}(u,v,w)=1\\}$

and therefore, by Lemma 3.1, one obtains

$\mathcal{I}_{{}_{F}}^{\tau}=\\{(u,v,w)\in\mathbb{R}^{3}\
:Q_{{}_{h}}(\frac{u}{\tau},\frac{v}{\tau},\frac{w}{\tau})=1\\}+\tau W.$ (3.5)

Eq. (3.5), that is the equation of spherical wavefront at time $\tau$ and
center $p$, is translation of the rotated ellipsoid given by Eq. (3.4), in
which the constant numbers are determined from the experimental data, and so
the proof of item $(i)$ is complete.

$(ii)$: Given some point $p\in A$, by Proposition 2.4 and since geodesics of
Minkowski spaces are straight lines, the wave rays emanating from the point
$p$ are unit speed straight lines ${\gamma_{{}_{F}}}(t)=p+tV$, provided that
vectors $V$ are $F$-orthogonal to $A$. Therefore, by Corollary 2.2 and
relation 2.2, we must have $V-W\underset{h}{\perp}A$ and $|V-W|_{{}_{h}}=1$.

$(iii)$: By Proposition 2.4, the wavefront at time $\tau$ is
$\rho^{-1}(\tau)$, where $\rho:M\to\mathbb{R}$, $\rho(p)=d(A,p)$. Given
$q\in\rho^{-1}(\tau)$, assume that ${\gamma_{{}_{F}}}$ is a unit speed
geodesic (wave ray) that minimizes the distance from $A$ to $q$. That is it
emanates from some point $p\in A$ and reaches to $q$ at time $\tau$, or in
other words,

$\tau=d(A,q)=d(p,q)=d({\gamma_{{}_{F}}}(0),{\gamma_{{}_{F}}}(\tau)).$

Note that, by item $(ii)$, ${\gamma_{{}_{F}}}$ is the straight line
${\gamma_{{}_{F}}}(t)=p+tV$ with $|V-W|_{{}_{h}}=1$ and
$V-W\underset{h}{\perp}A$. It concludes that

$\rho^{-1}(\tau)=\\{q={\gamma_{{}_{F}}}(\tau)=p+\tau V\ :\ p\in
A,|V-W|_{{}_{h}}=1,V-W\underset{h}{\perp}A\\}.$

∎

######  Corollary 3.4.

Assume that a wildfire is spreading in some space $M$ while the wind
$W=(0,W_{2},W_{3})$ is blowing across $M$ and $\mathcal{VF}_{\tau_{1}}$ is the
wavefront at time $\tau_{1}\geq 0$. Then the wavefront at time $\tau_{2}$ is
the envelope of ellipsoids
$Q_{{}_{h}}(\frac{u}{\tau},\frac{v}{\tau},\frac{w}{\tau})+\tau W+p,$ where
$\tau=\tau_{2}-\tau_{1}$, $p\in\mathcal{VF}_{\tau_{1}}$, and
$Q_{{}_{h}}(u,v,w)$ is given by Eq. 3.4.

###### Proof.

We consider the Randers distance function $\rho:M\to\mathbb{R},$
$\rho(p)=d(\mathcal{VF}_{0},p)$, where $\mathcal{VF}_{0}$ is the wavefront at
time $0$. Therefore, the proof is a direct consequence of Proposition 2.4 and
item $(i)$ of Theorem 3.2. ∎

It is notable that there may not be any relations between the angle $\alpha$
and angles that the wind makes with the axes. We can only say that the
stronger the wind, the closer the angle $\alpha$ to the angle that the wind
makes with the $z$-axis. In fact, since the heat goes towards above, the
indicatrix is an ellipsoid with the major axis along the $z$-axis, before the
influence of the wind. Then, the wind makes it rotate and translate toward the
direction of $W$. Consequently, for the stronger wind we have the smaller
deviation between the angle of rotation and the direction of $W$.

####  3.1.2 Strategic paths and points in the case of constant wind

Here, we determine the equations of strategic paths for two different
situations. The situations are as follows.

  * 1-

All the points of $M$ have the same priority from the fire fighting point of
view.

Therefore, the strategic path is the path along which the fire engulfs more
area.

  * 2-

Some special area or point has higher priority and the objective is protecting
it against the wildfire.

Hence, for a given point, the strategic path is the path through which the
fire particles reach to this point at some time $\tau$. Or, for a given area
$B$, the strategic path is the path through which the fire meets $\partial B$
for the first time at time $\tau$, where $\partial B$ is the frontier of $B$.
By the way, depending on the case, we may have more than one strategic path.

######  Proposition 1.

Assume that a wildfire is spreading in the space $M$, the constant wind
$W=(0,W_{2},W_{3})$ is blowing across $M$, and $A$ is the wavefront at time
$0$. Then:

  * (i)

In the case that all points of $M$ have the same priority, the strategic path
is $\gamma_{{}_{F}}(t)=p+tV$ provided that $p\in A$ and the Euclidean norm of
$V$, i.e. $v_{1}^{2}+v_{2}^{2}+v_{3}^{2}$, is the maximum among all the
vectors $V$ satisfying $|V-W|_{{}_{h}}=1$ and $V-W\underset{h}{\perp}A$.

  * (ii)

Given any point $q$, the strategic path that reaches to $q$ is
$\gamma_{{}_{F}}(t)=q+(t-\tau)V$ where $|V-W|_{{}_{h}}=1,$
$V-W\underset{h}{\perp}A$, and $\tau$ is the time of the wavefront to which
$q$ belongs.

  * (iii)

Given any area $B$, the strategic path which reaches to $B$ is
$\gamma_{{}_{F}}(t)=q+(t-\tau)V$ such that
$|V-W|_{{}_{h}}=1,V-W\underset{h}{\perp}A$, and $\tau$ is the time of the
wavefront that intersects $B$ for the first time and $q$ is the point of
intersection.

###### Proof.

$(i)$: First, since the strategic path is some wave ray of the fire that
emanates from $A$, by the item $(ii)$ of Theorem 3.2, it must be
$\gamma_{{}_{F}}(t)=p+tV$ such that $p\in A,|V-W|_{{}_{h}}=1,$ and
$V-W\underset{h}{\perp}A$. Next, the fact that the strategic path is the wave
ray through which the fire engulfs more region implies that the Euclidean
velocity of such a wave ray is the maximum among all these paths
$\gamma_{{}_{F}}$. In other words, $v_{1}^{2}+v_{2}^{2}+v_{3}^{2}$ is the
maximum among all such vectors $V$ satisfying in above conditions.

$(ii)$:  For a given point $q$, assuming that $\tau$ is the time when the
wavefront meets $q$ for the first time, there exists a unique wave ray that
emanates from some point belonging to $A$ and reaches to $q$ at time $\tau$
(up to a reparameterization). Because the wave rays are integral curves of the
gradient of distance function [9]. Hence, from item $(ii)$ of Theorem 3.2, the
path of this wave ray must be a unit speed straight line that emanates form
$A$, $F$-orthogonally, and reaches to the wavefront at time $\tau$. Therefore,
it is not difficult to show that $\gamma_{{}_{F}}(t)=q+(t-\tau)V$ where
$|V-W|_{{}_{h}}=1,$ and $V-W\underset{h}{\perp}A$.

$(iii)$:  First we find the wavefront that meets $B$ for the first time and
then the intersection of this wavefront and $B$. Assuming that $q$ is the
point of intersection, the rest of the proof is similar to that of item
$(ii)$.

∎

###  3.2 Wind as Some Smooth Vector Field

Here, through Theorem 3.5, we provide the model for the case that the wind is
a  Killing vector field. The positive aspect of this kind of vector field is
that there exists a direct relation between the wave rays of propagation and
geodesics of the Riemannian metric $h$. Therefore, in order to find the wave
rays, one just needs to solve the system of $h$-geodesic equations. To see the
system of geodesic equations of a Riemannian metric, see Chapter $6$ of [18].
In Theorem 3.6, we provide the propagation model for the case that the wind is
some smooth vector field, not necessarily Killing.

######  Theorem 3.5.

Assume that a wildfire is spreading in the space $M$, the wind $W$, which is a
Killing vector filed, is blowing across $M$ and $A$ is the wavefront at time
$0$. Then:

  * (i)

Given any point $p$ in $A$, the wave rays emanating from $p$ are
${\gamma_{{}_{F}}}(t):=\varphi(t,{\gamma_{{}_{h}}}(t))$, where $\varphi$ is
the flow of $W$ and $\gamma_{{}_{h}}$ is the $h$-geodesic such that
${\gamma_{{}_{h}}}(0)=p$, $|{\gamma_{{}_{h}}}^{\prime}(t)|=1$, and
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}(0)\underset{h}{\perp}A$.

  * (ii)

The spherical wavefront at time $\tau$ and center of some point $p\in M$ is
the set

$\\{\varphi(\tau,{\gamma_{{}_{h}}}(\tau)):{\gamma}_{{}_{h}}\text{ is the unit
speed $h$-geodesic that }{\gamma}_{{}_{h}}(0)=p\\}$.

  * (iii)

The wavefront at time $\tau$ is the set

$\\{\varphi(\tau,{\gamma_{{}_{h}}}(\tau)):\gamma_{{}_{h}}\text{ is an
$h$-geodesic that
}d\varphi_{{}_{p}}{\gamma_{{}_{h}}}^{\prime}\underset{h}{\perp}A\,\text{and}\,|{\gamma_{{}_{h}}}^{\prime}|=1\\}.$

###### Proof.

Before proceeding with the proofs of items $(i),(ii)$, and $(iii)$, let’s see
what the indicatrix $Q_{h}$ and Riemannian metric $h$ are. Given $p\in M$, one
is faced with Theorem 3.2 on $T_{p}M$ with constant wind $W(p)$ on it.
Although, since the direction of $W$ depends on $p$, one cannot assume that
$W$ is always in $yz$-plane. That is the indicatrix (as a subset of $T_{p}M$)
is some translation of an ellipsoid rotated along $x$-, $y$-, and $z$\- axes.
Therefore, the indicatrix is the translation of the rotated ellipsoid

$\displaystyle
Q_{{}_{h}}(p;(u,v,w))=(\frac{u^{\prime}}{a(p)})^{2}+(\frac{v^{\prime}}{b(p)})^{2}+(\frac{w^{\prime}}{c(p)})^{2}=1,$
(3.6)

by the vector $W(p)$, where $\begin{bmatrix}u^{\prime}\\\ v^{\prime}\\\
w^{\prime}\\\ \end{bmatrix}=P(p)\begin{bmatrix}u\\\ v\\\ w\\\ \end{bmatrix}$,
$P(p)=R_{z}({\theta(p)})R_{y}({\beta(p)})R_{x}({\alpha(p)})$, and

$\displaystyle R_{x}(\alpha(p))=\left(\begin{matrix}1&0&0\\\
0&\cos\alpha(p)&-\sin\alpha(p)\\\
0&\sin\alpha(p)&\cos\alpha(p)\end{matrix}\right),R_{y}(\beta(p))=\left(\begin{matrix}\cos\beta(p)&0&\sin\beta(p)\\\
0&1&0\\\ -\sin\beta(p)&0&\cos\beta(p)\end{matrix}\right),$ $\displaystyle
R_{z}(\theta(p))=\left(\begin{matrix}\cos\theta(p)&-\sin\theta(p)&0\\\
\sin\theta(p)&\cos\theta(p)&0\\\ 0&0&1\end{matrix}\right).$

Here $a,b,c,\alpha,\beta,$ and $\theta$ are smooth functions and at each point
$p$ they are determined from the experimental data. By the way, since $P(p)\in
SO(3)$, $\mathbb{B}=\\{P\frac{\partial}{\partial x},P\frac{\partial}{\partial
y},P\frac{\partial}{\partial z}\\}$ is a rotation of the canonical basis
$\\{\frac{\partial}{\partial x},\frac{\partial}{\partial
y},\frac{\partial}{\partial z}\\}$ and therefore it is an orthogonal basis for
the space. It is not difficult to show that the matrix of Riemannian metric is
$\hbar(p)=(P^{T}\mathcal{D}P)(p)$, where
$\mathcal{D}(p)=diag(\lambda_{1}(p),\lambda_{2}(p),\lambda_{3}(p))$.

$(i)$: First, by Proposition 2.4, the wave rays are unitary $F$-geodesics that
are $F$-orthogonal to $A$. To find unitary $F$-geodesics, by Lemma 2.1 it is
enough to find unitary $h$-geodesics. Therefore, for any $F$-geodesic
$\gamma_{{}_{F}}(t)$ and its corresponding $h$-geodesic
${\gamma_{{}_{h}}}(t)$, from this property of flow that $\varphi(0,p)=p$, one
has
${\gamma_{{}_{F}}}(0)=\varphi(0,{\gamma_{{}_{h}}}(0))={\gamma_{{}_{h}}}(0)$;
that means both ${\gamma_{{}_{F}}}$ and $\gamma_{{}_{h}}$ have the same
initial point. Hence, given $p\in A$, to find the wave rays emanating from
$p$, it suffices to find the unitary $h$-geodesics $\gamma_{{}_{h}}(t)$
emanating from $p$ provided that
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}(0)\underset{h}{\perp}A$. The last
relation satisfies the condition that $\gamma_{{}_{F}}(t)$ must be
$F$-orthogonal to $A$. Because by Corollary 2.2,
${\gamma_{{}_{F}}}^{\prime}(0)\underset{F}{\perp}A$ if and only if
${\gamma_{{}_{F}}}^{\prime}(0)-W\underset{h}{\perp}A$. Furthermore, from the
chain rule we have

${\gamma_{{}_{F}}}^{\prime}(0)=d\varphi(0,{\gamma_{{}_{h}}}(0))+d\varphi(0,{\gamma_{{}_{h}}}(0)){\gamma_{{}_{h}}}^{\prime}(0)=W(p)+d\varphi_{{}_{p}}{\gamma_{{}_{h}}}^{\prime}(0).$

$(ii)$: By Proposition 2.4, the spherical wavefront of radius $\tau$ and
center $p$ is the set $\rho^{-1}(\tau)$, where
$\rho(.):=d(p,.):M\to\mathbb{R}$. Given $q\in\rho^{-1}(\tau)$, assume that
${\gamma_{{}_{F}}}(t)$ is the equation of a wave ray from $p$ to $q$. Hence,
by Proposition 2.4, ${\gamma_{{}_{F}}}(t)$ is a unitary $F$-geodesic that
minimizes the distance from $p$ to $q$, or in other words

$d(p,q)=d(\gamma_{{}_{F}}(0),\gamma_{{}_{F}}(\tau))=\tau.$

Therefore,

$\rho^{-1}(\tau)=\\{{\gamma_{{}_{F}}}(\tau):{\gamma_{{}_{F}}}(t)\text{ is an
$F$-geodesic that
}d(p,{\gamma_{{}_{F}}}(\tau))=\tau,F({\gamma_{{}_{F}}}^{\prime})=1,{\gamma_{{}_{F}}}(0)=p\\}.$

Now, by item $(i)$, ${\gamma_{{}_{F}}}(t)=\varphi(t,{\gamma_{{}_{h}}}(t))$
where ${\gamma_{{}_{h}}}(t)$ is a unit speed $h$-geodesic such that
${\gamma_{{}_{h}}}(0)=p$, concluding the proof of item $(ii)$.

$(iii)$: The wavefront at time $\tau$, by Proposition 2.4, is the set
$\rho^{-1}(\tau)$, where $\rho(.)=d(A,.):M\to\mathbb{R}$. In other words,

$\rho^{-1}(\tau)=\\{q:\ d(A,q)=\tau\\}.$

Suppose ${\gamma_{{}_{F}}}(t)$ is the equation of some unit speed $F$-geodesic
that minimizes the distance from $A$ to $q\in\rho^{-1}(\tau)$. That is
${\gamma_{{}_{F}}}(0)\in A$, $q={\gamma_{{}_{F}}}(\tau)$,
$\gamma_{{}_{F}}^{\prime}(0)\underset{F}{\perp}A$ and

$d(A,q)=d({\gamma_{{}_{F}}}(0),{\gamma_{{}_{F}}}(\tau))=\tau.$

Moreover, by item $(i)$,
${\gamma_{{}_{F}}}(t)=\varphi(t,{\gamma_{{}_{h}}}(t))$ where
${\gamma_{{}_{h}}}(t)$ is a unitary $h$-geodesic such that
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}(0)\underset{h}{\perp}A$.
Therefore, we have

$\rho^{-1}(\tau)=\\{\varphi(\tau,{\gamma_{{}_{h}}}(\tau)):\,\gamma_{{}_{h}}(t)\text{
is an
}\text{h-geodesic},\,|{\gamma_{{}_{h}}}^{\prime}|_{{}_{h}}=1,d\varphi_{{}_{p}}{\gamma_{{}_{h}}}^{\prime}\underset{h}{\perp}A\,\\},$

which concludes the proof.

∎

######  Remark 1.

We have a corollary similar to Corollary 3.4. That is to model a wildfire
propagation in some space $M$ under the presence of the wind $W$ which is
Killing, one can find the spherical wavefronts, centers at different points of
some wavefront, by item $(ii)$ of Theorem 3.5 and then the envelope.

####  3.2.1 Strategic paths and points in the case of Killing vector field

Through Proposition 2 below, we verify the strategic paths and points in the
case of wind being a Killing vector field.

######  Proposition 2.

Assume that some wildfire is spreading in the space $M$ and the wind $W$, that
is a Killing vector filed, is blowing throughout $M$ and $A$ is the wavefront
at time $0$. Then:

  * (i)

In the case that all points of $M$ have the same priority, given any time
$\tau$, the strategic path until time $\tau$ is
$\varphi(t,{\gamma_{{}_{h}}}(t))$, where $\gamma_{{}_{h}}$ is the $h$-geodesic
such that $|{\gamma_{{}_{h}}}(\tau)-{\gamma_{{}_{h}}}(0)|$ is maximum among
all $h$-geodesics that $|{\gamma_{{}_{h}}}^{\prime}(t)|_{{}_{h}}=1$ and
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}\underset{h}{\perp}A$.

  * (ii)

Given any point $q$ belonging to the wavefront at time $\tau$, the strategic
path which reaches to $q$ is $\varphi(t,{\gamma_{{}_{h}}}(t))$, where
$\gamma_{{}_{h}}$ is the $h$-geodesic such that
$|{\gamma_{{}_{h}}}^{\prime}(t)|=1$,
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}\underset{h}{\perp}A$ and
$\varphi(\tau,{\gamma_{{}_{h}}}(\tau))=q$.

  * (iii)

Given any area $B$, the strategic path that reaches to it is
$\varphi(t,{\gamma_{{}_{h}}}(t))$ where $\gamma_{{}_{h}}$ is the unitary
$h$-geodesic such that
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}(0)\underset{h}{\perp}A$ and
$\varphi(\tau,{\gamma_{{}_{h}}}(\tau))=q$. Here, $\tau$ is the time of the
wavefront that intersects $B$ for the first time and $q$ is the point of
intersection.

###### Proof.

To prove item $(i)$, each strategic path must be an $F$-geodesic. Therefore,
by item $(i)$ of Theorem 3.5, the strategic path must be
$\varphi(t,{\gamma_{{}_{h}}}(t))$, where $\varphi$ is the flow of $W$ and
$\gamma_{{}_{h}}$ is the unitary $h$-geodesic such that
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}(0)\underset{h}{\perp}A$. However,
as the wave rays might be some curves, one can not claim that there exists a
unique strategic path that remains valid for any time $t$. Indeed, a wave ray
might be a strategic path just before some time $\tau$ and after this time one
has to choose another wave ray as the strategic path. Since the strategic path
is the path through which the fire engulfs more area, the Euclidean length of
the strategic path is the maximum. That is
$|{\gamma_{{}_{h}}}(\tau)-{\gamma_{{}_{h}}}(0)|$ must be maximum among all
wave rays joining $A$ to the wavefront at time $\tau$, closing the proof of
item $(i)$.

To prove item $(ii)$, similar to item $(i)$, the strategic path is
$\varphi(t,{\gamma_{{}_{h}}}(t))$, where $\gamma_{{}_{h}}$ is the unitary
$h$-geodesic such that
$d\varphi_{{}_{p}}{\gamma^{\prime}_{{}_{h}}}(0)\underset{h}{\perp}A$. However,
among all such wave rays, we have to find the one that passes through $q$ at
time $\tau$. In other words, the wave ray that satisfies the condition
$\varphi(\tau,{\gamma_{{}_{h}}}(\tau))=q$ is what we need.

The proof of item $(iii)$, is similar to that of item $(iii)$ of Proposition
1. ∎

In the next theorem, we present the model for wildfire propagation when the
wind $W$ is a smooth vector field, not necessarily Killing.

######  Theorem 3.6.

Assume that some wildfire is spreading in a space $M$ across which the wind
$W$, that is a smooth vector field, is blowing. Suppose that $A$ is the
wavefront at time $0$. Then:

  * (i)

The wave rays are unit speed $F$-geodesics that are $F$-orthogonal to $A$,
where $F$ is the Randers metric whose indicatrix at each point $p\in M$ is
translation of the rotated ellipsoid, given by Eq. (3.6), by $W(p)$.

  * (ii)

Given some point $p\in M$, the spherical wavefront of radius $\tau$ and center
$p$ is:

$\\{{\gamma_{{}_{F}}}(\tau):{\gamma_{{}_{F}}}(\tau)\text{ is the unit speed
$F$-geodesic that }{\gamma_{{}_{F}}}(0)=p\\}$.

  * (iii)

The wavefront at time $\tau$ is the set:

$\\{{\gamma_{{}_{F}}}(\tau):{\gamma_{{}_{F}}}(\tau)\text{ is the unit speed
$F$-geodesic that }{\gamma_{{}_{F}}}^{\prime}(0)\underset{F}{\perp}A\\}.$

###### Proof.

To prove item $(i)$, it should be noted that once one has the Randers metric
$F$ on $M$, by Proposition 2.4, the wave rays are unitary $F$-geodesics that
are $F$-orthogonal to $A$. The metric $F$ is given by Eq. (2.1) which demands
finding the Riemannian metric $h$. To find $h$, one has to find the equation
of Riemannian/Randers indicatrix, as it was explained after Lemma 3.1. By
following the same argument as that of Theorem 3.5, it can be shown that at
each point $p\in M$ the Randers indicatrix is translation of the rotated
ellipsoid, given by Eq. (3.6), by the vector $W(p)$. It completes the proof of
item $(i)$.

The proofs of items $(ii)$ and $(iii)$ are similar to those of items $(ii)$
and $(iii)$ of Theorem 3.5, respectively, in which one does not involve
$\gamma_{{}_{h}}(t)$ and $\varphi(.,.)$ in the arguments. ∎

####  3.2.2 Strategic paths and points in the case of wind being a smooth
vector field

Here, we verify the strategic paths and points in the case that wind is not
necessarily a Killing vector field.

######  Corollary 3.7.

Assume a wildfire is spreading in some space $M$, the wind $W$, that is a
smooth vector filed, is blowing throughout $M$ and $A$ is the wavefront at
time $0$. Then:

  * (i)

If all the points of $M$ have the same priorities, given any time $\tau$, the
strategic path till time $\tau$ is the unitary $F$-geodesic
$\gamma_{{}_{F}}(t)$ that is $F$-orthogonal to $A$, provided that
$|{\gamma_{{}_{F}}}(\tau)-{\gamma_{{}_{F}}}(0)|$ is maximum among all such
$F$-geodesics.

  * (ii)

Given any point $q$ belonging to the wavefront at time $\tau$, the strategic
path which reaches to $q$ is the unitary $F$-geodesic $\gamma_{{}_{F}}(t)$
that is $F$-orthogonal to $A$ and ${\gamma_{{}_{F}}}(\tau)=q$.

  * (iii)

Given any area $B$, the strategic path that meets $B$ for the first time is
the unitary $F$-geodesic $\gamma_{{}_{F}}(t)$ which is $F$-orthogonal to $A$
and ${\gamma_{{}_{F}}}(\tau)=q$. Here, $\tau$ is the time of the wavefront
that intersects $B$ for the first time and $q$ is the point of intersection.

###### Proof.

One applies Theorem 3.6 and follows a proof similar to that of Proposition 2,
in which the flow of $W$ and Riemannian geodesic $h$ do not get involved. ∎

##  4 Examples

Here, we provide two examples in which we assume that a wildfire is spreading
across some agricultural land or woodland $M\subset\mathbb{R}^{3}$. In fact we
assume that $M$ is some open subset of $\mathbb{R}^{3}$. It is assumed that
the fuel has been distributed smoothly across $M$ and some wind $W$ is blowing
throughout $M$. In Example 4.1, the wind is constant and in Example 4.2, it is
a Killing vector field. In both examples, the fire starts from ${A}\subset M$,
which is a certain point, smooth curve or cylinder. The objective in these
examples is finding the wavefronts, wave rays and strategic paths. It is
assumed that $M$ contains some flat field of fuel $D$, that is $D$ is a
$2$-dimensional subspace of $M$ such that its slope is negligible. In fact, if
one walks around in $D$, he would feel some going up and down but negligible.
For instance, $D$ might be the bed of agricultural land, a forest, and so
forth.

######  Example 4.1.

Assume that a wildfire is spreading in some agricultural land that $M$ and the
wind $W=(0,1/3,1/6)$ is blowing across $M$. We consider three different cases
for $A$, that is the wavefront at time $0$, as follows:

  * Case 1.

$A$ is some point in $D$.

  * Case 2.

$A$ is the path of closed curve

$C(s)=\big{(}\frac{1}{4}\cos s(\cos s+6),\frac{4}{13}\sin s(-\sin
s+3),0\big{)},s\in[0,2\pi].$ (4.1)

  * Case 3.

$A$ is the image of surface

$S(s_{1},s_{2})=\big{(}\frac{1}{4}\cos s(\cos s+6),\frac{4}{13}\sin s(-\sin
s+3),s_{2}\big{)},s_{1}\in[0,2\pi],\ s_{2}\in[0,2].$ (4.2)

Before verifying the cases, let’s see what the Riemannian metric and
indicatrix are. Since the wind is a constant vector field, by Theorem 3.2, the
indicatrix is translation of the rotated ellipsoid given by Eq. (3.4), by $W$.
Assume that, from experimental data, we are given the constant numbers in Eq.
(3.4) as follows:

$\displaystyle a=1/2,\ b=1\ ,c=2,\hskip 5.69054pt$
$\displaystyle\alpha=\pi/6.$

Therefore, the equation of the spherical wavefront is Eq. (4.3) below:

$64u^{2}+13v^{2}+7w^{2}-(26/3+\sqrt{3})v+(2\sqrt{3}-7/3)w-6\sqrt{3}vw=\frac{\sqrt{3}}{3}+\frac{635}{36}.$
(4.3)

Hence, the matrix of the Riemannian metric is

$\displaystyle\hbar=\left(\begin{matrix}4&0&0\\\
0&\frac{13}{16}&-\frac{3\sqrt{3}}{16}\\\
0&-\frac{3\sqrt{3}}{16}&\frac{7}{16}\end{matrix}\right).$

In the sequel, we find the wavefronts, wave rays and strategic paths for each
of the cases listed above, separately.

Case $1$

We consider the point $A$ as the origin of the coordinate system. By Theorem
3.2, the wavefront at time $\tau$ is translation of

$4\big{(}\frac{u}{\tau}\big{)}^{2}+\frac{13}{16}\big{(}\frac{v}{\tau}\big{)}^{2}+\frac{7}{16}\big{(}\frac{w}{\tau}\big{)}^{2}-\frac{3\sqrt{3}}{8}\big{(}\frac{v}{\tau}\big{)}\big{(}\frac{w}{\tau}\big{)}=1,$

by $\tau W$.

The set of wave rays, by item $(ii)$ of Theorem 3.2, is

$\\{\gamma_{{}_{F}}(t)=tV:\ |V-W|_{{}_{h}}=1\\},$

where $V=(v_{1},v_{2},v_{3})$ and the equality $|V-W|_{{}_{h}}=1$ is
equivalent to

$64{v_{1}}^{2}+13v_{2}^{2}+7v_{3}^{2}-(\frac{26}{3}+\sqrt{3})v_{2}+(2\sqrt{3}-\frac{7}{3})v_{3}-6\sqrt{3}v_{2}v_{3}=\frac{\sqrt{3}}{3}+\frac{635}{36}.$
(4.4)

In order to track the strategic path which reaches to the wavefront at time
$\tau$, it suffices to find the point $q=(u,v,w)$ for which
$\sqrt{u^{2}+v^{2}+w^{2}}$ is maximum among all other points belonging to this
wavefront. Afterwards, the strategic path is
$\gamma_{{}_{F}}(t)=t\frac{q}{\tau}$. The Fig. 1 depicts, in $D$, wavefronts
at $t=1\,,2\,,\cdots,\,10$ and the strategic path from $t=0$ to $t=10$.

Figure 1:  some of wavefronts and strategic path in $D$ from $1$ to $10$ for
Case 1

Case $2$

To find the wavefront at time $\tau$, one may use the spherical wavefronts
given by Eq. (4.3) and centers at different points belonging to $C([0,\pi])$,
and then apply the Huygens’ principle.

The wave ray that emanates from some point $p=C(s_{p})$ is
$\gamma_{{}_{F}}(t)=p+tV$, where $V=(v_{1},v_{2},v_{3})$ is the solution of
the following equations system at point $p$.

$\left\\{\begin{array}[]{l}\langle V-W,C^{\prime}(s_{p})\rangle=0\\\
|V-W|_{{}_{h}}=1.\\\ \end{array}\right.$ (4.5)

Here, in the system (4.5), $|V-W|_{{}_{h}}=1$ is equivalent to Eq. (4.4) and
the first equation is equivalent to

$-v_{1}\big{(}\sin 2s_{p}+6\sin
s_{p}\big{)}+\big{(}\frac{13}{16}(v_{2}-\frac{1}{3})-\frac{3\sqrt{3}}{16}(v_{3}-\frac{1}{6})\big{)}\frac{1}{13}\big{(}-\sin
2s_{p}+3\cos s_{p}\big{)}=0.$

In the case that all points of wavefront have the same priority, to find the
strategic path that reaches to the wavefront at time $10$ one has to find the
wave ray ${\gamma}_{{}_{F}}(t)=p+tV$ for which $v_{1}^{2}+v_{2}^{2}+v_{3}^{2}$
is maximum. Fig. 2 depicts some of the wavefronts from time $1$ to $10$, the
strategic path (in purple color) from till time $10$, and the path (in black
color) through which the fire is progressing slower.

Figure 2: some of wavefronts from $1$ to $10$ and strategic path (purple) for
Case 2 in $D$

Case $3$

In order to find the wavefront at time $\tau$, one way is finding the
spherical wavefronts with centers at points $p\in S([0,2\pi]\times[0,2])$, and
then the hypersurface which is tangent to all of these spherical wavefronts.

To discover the wave ray that emanates from some point $p=S(s_{1,p},s_{2,p})$,
i.e. $\gamma_{{}_{F}}(t)=p+tV$, one has to find vectors $V$ that satisfies the
following system

$\left\\{\begin{array}[]{l}\langle V-W,\frac{\partial S}{\partial
s_{1}}\rangle=0,\\\ \langle V-W,\frac{\partial S}{\partial s_{2}}\rangle=0,\\\
|V-W|_{{}_{h}}=1,\\\ \end{array}\right.$ (4.6)

provided that $V-W$ is in the direction of the propagation. The first equation
in the system 4.6 is equivalent to

$-\frac{2}{5}v_{1}\sin s_{1}(2\cos s_{1}+1)+\frac{1}{128}\cos s_{1}(2\sin
s_{1}+1)\big{(}13(v_{2}-\frac{1}{3})+3\sqrt{3}(v_{3}-\frac{1}{6})\big{)}=0,$

and the second equation is equivalent to

$\frac{3\sqrt{3}}{7}(v_{2}-\frac{1}{3})=v_{3}-\frac{1}{6}.$

Once one has the wave rays, the set $\\{\gamma_{{}_{F}}(\tau)=p+\tau V:\ p\in
S([0,\pi]\times[0,2])\\}$ is the wavefront at time $\tau$.

If all the points of $M$ have the same priority, the strategic path that
reaches to the wavefront at time $\tau$ is the wave ray
${\gamma}_{{}_{F}}(t)=p+tV$ for which $v_{1}^{2}+v_{2}^{2}+v_{3}^{2}$ is
maximum.

######  Example 4.2.

Assume that a wildfire is spreading in some agricultural land and the wind $W$
is blowing across it. we have $W(p)=k(y,0,0)$, where $p=(x,y,z)\in M$ and $k$
is some small enough constant number so that the origin of tangent space
remains inside the indicatrix. Similar to the Example 4.1, we consider three
different cases for the wavefront at time $0$, i.e. $A$. In other words, $A$
is a point or some curve given by Eq. (4.1) or surface given by Eq. (4.2). As
the first step, one has to find the indicatrix. Since the wind is some smooth
vector field, by Theorems 3.5 and 3.5, the indicatrix must be translation of
the rotated ellipsoid, given by Eq. (3.6), by the vector $W$. Keeping in mind
Eq. (3.6), assume that from experimental data we are given $a=1,\
b=\frac{1}{2},\ c=2$, $\alpha=\theta=0$, and $\beta=y$. Therefore, the
indicatrix is translation of the following ellipsoid by $W$:

$\displaystyle Q_{{}_{h}}(u,v,w)=({u\cos y+w\sin y})^{2}+4v^{2}+(\frac{-u\sin
y+w\cos y}{2})^{2}=1.$ (4.7)

From Eq. (4.7), the Riemannian metric is $h=P^{T}\mathcal{D}P$, where
$\mathcal{D}=diag(1,\frac{1}{2},2)$ and $P(x,y,z)=R_{y}(y)$. From the fact
that $\frac{\partial}{\partial x}(P_{ki}P_{kj})=0$, it is deduced that $W$ is
a Killing vector field. It is easy to see that the flow of $W$ is as
$\varphi(t,p)=kt(y,0,0)+p$, where $p=(x,y,z)$. We apply Theorem 3.5 to deal
with the propagation and obtain wavefronts, wave rays, and strategic paths for
each case.

Case $1$

The wave rays emanating from some point $p$ and the unitary velocity vector
$V$ are

$\gamma_{{}_{F}}(t)=kt(y(t),0,0)+(x(t),y(t),z(t)),$

where $\gamma_{{}_{h}}(t)=(x(t),y(t),z(t))$ is solution of the following
system

$\displaystyle\left\\{\begin{array}[]{l}x^{\prime\prime}+y^{\prime}z^{\prime}=0,\\\
y^{\prime\prime}=0,\\\ z^{\prime\prime}-x^{\prime}y^{\prime}=0,\\\
\end{array}\right.$ (4.11)

with the initial conditions $\gamma_{{}_{h}}(0)=p=(x,y,z)$,
$\gamma_{{}_{h}}^{\prime}(0)=V-W$ and $|V-W|_{{}_{h}}=1$. After some
calculations, if $v_{2}\neq 0$, we have

$\displaystyle\left\\{\begin{array}[]{l}x(t)=x-\frac{v_{3}}{v_{2}}+t\frac{v_{3}}{v_{2}}\cos
v_{2}+t\frac{v_{1}-ky}{v_{2}}\sin v_{2},\\\ y(t)=tv_{2}+y,\\\
z(t)=z-\frac{v_{1}-ky}{v_{2}}-t\frac{v_{1}-ky}{v_{2}}\cos
v_{2}+t\frac{v_{3}}{v_{2}}\sin v_{2}.\\\ \end{array}\right.$ (4.15)

and if $v_{2}=0$,

$\displaystyle\left\\{\begin{array}[]{l}x(t)=(v_{1}-ky)t+x,\\\ y(t)=y,\\\
z(t)=v_{3}t+z,\\\ \end{array}\right.$ (4.19)

are component of $\gamma_{{}_{h}}(t)$, provided that $|V-W|_{{}_{h}}=1$. One
can easily verify that $|V-W|_{{}_{h}}=1$ is equivalent to

$4(v_{1}-ky)^{2}+\frac{13}{16}v_{2}^{2}-\frac{3\sqrt{3}}{8}v_{2}v_{3}+\frac{7}{16}v_{3}^{2}=1.$
(4.20)

Once one has the curves $\gamma_{{}_{h}}(t)=(x(t),y(t),z(t))$, which are
satisfied in the systems 4.15 or 4.19, the wavefront at time $\tau$ is the set

$\\{k\tau(y(\tau),0,0)+(x(\tau),y(\tau),z(\tau))\\}.$

The strategic path for the situation that all points belonging to the space
have the same priority is the wave ray
$\gamma_{{}_{F}}(t)=kt(y,0,0)+\gamma_{{}_{h}}(t)$ provided that
$|k\tau(y,0,0)+\gamma_{{}_{h}}(\tau)-p|$ is maximum among all the curves
$\gamma_{{}_{h}}(t)$. If one wants to find the wave ray that reaches to some
point $q$ belonging to the wavefront at time $\tau$, it is
$\gamma_{{}_{F}}(t)=kt(y,0,0)+\gamma_{{}_{h}}(t)$ such that
$\gamma_{{}_{h}}(\tau)=q-k\tau(y,0,0)$.

Case $2$

For Case $2$, the wave ray which emanates from some point $p=C(s_{p})$ is
$\gamma_{{}_{F}}(t)=kt(y(t),0,0)+(x(t),y(t),z(t))$ such that
$\gamma_{{}_{h}}(t)=(x(t),y(t),z(t))$ is a solution of the systems 4.15 or
4.19 and it also satisfies $|\gamma^{\prime}_{{}_{h}}(0)|=1$,
$\gamma_{{}_{h}}(0)=p$, and
$\gamma_{{}_{h}}^{\prime}(0)\underset{h}{\perp}C^{\prime}(s)$ in the direction
of the propagation.

The wavefront at time $\tau$ is the set

$\\{k\tau(y(\tau),0,0)+(x(\tau),y(\tau),z(\tau))\\}.$

The calculations related to strategic paths, and wave rays and wavefront in
the Case $3$ are similar to the previous case and the Case $3$ of Example 4.1.

##  5 Conclusion and Final Remarks

In this work, we presented a model for wildfire propagation in some smooth
manifold under the presence of wind. The results and proofs can be generalized
to any dimension $n$, however to be closer to reality, we focused on the
dimension $3$. Three different cases for the wind were considered; a constant,
Killing, and smooth vector field. These three cases were separately studied
through Theorems 3.2, 3.5, and 3.6. This way, equations of wavefronts and wave
rays were provided. To summarize the results of these theorems, for a wildfire
spreading in some space, we provided two ways to find the model of propagation
as follows.

  * (i)

Using the Huygens’ principle. It means, given the wavefront at some time
$t_{0}\geq 0$, we find the spherical wavefronts of radius $r>0$ and center of
each at some point belonging to this wavefront. Then we find the envelope of
these spherical wavefronts; that is the wavefront at time $t_{0}+r$. This
method is more applicable when we want to use some computer programming and
software for modeling the propagation.

  * (ii)

Using geodesics. The unitary Randers geodesics that start orthogonally from
some wavefront are wave rays and all of them reach to the next wavefront at
the same time. It means given some wavefront $A$, after time $\tau$, the
wavefront will be the set

$\\{{\gamma_{{}_{F}}}(\tau)\ \\},$

where $\gamma_{{}_{F}}(t)$ are unitary Randers geodesics that emanates from
$A$ orthogonally.

One can find Randers geodesics without even calculating the metric. In
Theorems 3.2 and 3.5, the formulas of Randers geodesics, for the cases of wind
being a constant or Killing vector field, are provided. For the case of wind
being a smooth vector field, by Theorem 3.6 and by formulas and relations
provided in Section 2.3 of [6], one obtains the system of Randers geodesic
equations from the system of Riemannian geodesic equations. Therefore, given a
propagation, it suffices to find the Riemannian metric, and then the system of
geodesic equations associated with it, and finally the wavefronts.

By the way, concepts of strategic paths and points were introduced and their
equations were determined for different situations. These paths and points are
those which are of great importance in wildfire management strategies; by
providing the places where the firefighters and equipment should be located.
Two examples, in which some wildfire is spreading in some agricultural land,
were given to illustrate the results.

In this work, it was assumed that the wind $W$ is time-independent. However,
one can apply our results for the case that the wind is some time-dependent
vector field $W(t,p),$ $t\in[0,s]$, assuming that $\\{t_{1},....,t_{n}\\}$ is
a partition of $[0,s]$ such that at each interval
$[t_{i},t_{i+1}]\subset[0,s]$, the wind is the time-independent vector field
$W^{i}(p)$. Therefore, at each interval $[t_{i},t_{i+1}]$, we are faced with a
new propagation and have to apply one of the Theorems 3.2, 3.5 or 3.6 to
modeling it. Generally, there is no relations between the metric related to
one interval with that of another one.

By the way, in the case that we also need the model of wildfire propagation in
dimension $2$, that is on the land, it suffices to consider the intersection
of wavefronts/propagation with the land.

For further works in the future, one may work on the cases that the wildfire
creates singularities or cut loci. Also , it is interesting to see how the
wildfire spreads when the fuel has not been distributed smoothly through
space. Moreover, studying the propagation of wildfire on a mountain or some
slope is welcome as it is related to certain real situations. In this regard,
one problem can be finding the relations between propagation of wildfires on
two domains under the same conditions but different slopes. By the way, as
another problem, one may study the model associated with a strong wind, that
is when the origin of tangent space does not locate inside the indicatrix.

## References

  * [1] Anderson, D. H., Catchpole, E. A., De Mestre, N. J., and Parkes, T. Modelling the spread of grass fires.  Anziam J. 23, 4 (1982), 451–466.
  * [2] Arca, B., Duce, P., Laconi, M., Pellizzaro, G., Salis, M., and Spano, D. Evaluation of farsite simulator in mediterranean maquis.  IJWF 16, 5 (2007), 563–572.
  * [3] Arnold, V. I. Mathematical methods of classical mechanics, vol. 60.  New York, Springer, 2013.
  * [4] Bao, D., Chern, S. S., and Shen, Z. An introduction to Riemann-Finsler geometry.  Springer Science & Business Media, New York, Springer, 2012.
  * [5] Bao, D., Robles, C., and Shen, Z. Zermelo navigation on Riemannian manifolds.  J. Differential Geom. 66, 3 (2004), 377–435.
  * [6] Cheng, X., and Shen, Z. Finsler geometry.  An approach via Randers spaces (New York, Springer, 2012).
  * [7] Cunbin, L., Jing, Z., Baoguo, T., and Ye, Z. Analysis of forest fire spread trend surrounding transmission line based on rothermel model and huygens principle.  Int. J. Multimedia Ubiquitous Eng. 9, 9 (2014), 51–60.
  * [8] Dehkordi, H. R. Finsler Transnormal functions and singular foliations of codimension 1.  PhD thesis, IME, University of Sao paulo, Brazil, 2018.
  * [9] Dehkordi, H. R., and Saa, A. Huygens’ envelope principle in Finsler spaces and analogue gravity.  Classical Quant. Grav. 36, 8 (2019), 085008.
  * [10] Dore, S., Kolb, T. E., Montes-Helu, M., Eckert, S. E., Sullivan, B. W., Hungate, B. A., Kaye, J. P., Hart, S. C., Koch, G. W., and Finkral, A. Carbon and water fluxes from ponderosa pine forests disturbed by wildfire and thinning.  Ecol. Appl. 20 (2010), 663–683.
  * [11] Flannigan, M. D., Amiro, B. D., Logan, K. A., Stocks, B. J., and Wotton, B. M. Forest fires and climate change in the 21 st century.  Mitig. Adapt. Strateg. Glob. Chang. 11, 4 (2006), 847–859.
  * [12] Giambò, R., Giannoni, F., and Piccione, P. Genericity of nondegeneracy for light rays in stationary spacetimes.  Comm. Math. Phys. 287, 3 (2009), 903–923.
  * [13] Gibbons, G. W., and Warnick, C. M. The geometry of sound rays in a wind.  Contemp. Phys. 52, 3 (2011), 197–209.
  * [14] Glasa, J., and Halada, L. On mathematical foundations of elliptical forest fire spread model.  Chapter 12 (2009), 315–333.
  * [15] Hilton, J. E., Miller, C., Sharples, J. J., and Sullivan, A. L. Curvature effects in the dynamic propagation of wildfires.  IJWF 25, 12 (2017), 1238–1251.
  * [16] Johnston, P., Milne, G., and Klemitz, D. Overview of bushfire spread simulation systems.  Bushfire CRC Project B 6 (2005).
  * [17] Kopacz, P. Application of planar randers geodesics with river-type perturbation in search models.  Applied Mathematical Modelling 49 (2017), 531–553.
  * [18] Lee, J. M. Introduction to Riemannian manifolds.  New York, Springer, 2018.
  * [19] Markvorsen, S. A Finsler geodesic spray paradigm for wildfire spread modelling.  Nonlinear Anal. Real World Appl. 28 (2016), 208–228.
  * [20] Muresan, M. On Zermelo’s navigation problem with mathematica.  J. Funct. Anal. 9, 3-4 (2014), 349–355.
  * [21] Papadopoulos, G. D., and Pavlidou, F. N. A comparative review on wildfire simulators.  IEEE Syst J. 5, 2 (2011), 233–243.
  * [22] Richards, G. D. A general mathematical framework for modeling two-dimensional wildland fire spread.  IJWF 5, 2 (1995), 63–72.
  * [23] Robles, C. Geodesics in randers spaces of constant curvature.  Trans. Amer. Math. Soc. 359 (2005).
  * [24] Shen, Z. Lectures on Finsler geometry.  Singapura, World Scientific, 2001.
  * [25] Sullivan, A. L. Wildland surface fire spread modelling, 1990–2007. 3: Simulation and mathematical analogue models.  IJWF 18, 4 (2009), 387–403.
  * [26] Tymstra, C., Bryce, R. W., Wotton, B. M., Taylor, S. W., Armitage, O., et al. Development and structure of prometheus: the canadian wildland fire growth simulation model.  Natural Resources Canada, Canadian Forest Service, Northern Forestry Centre, Information Report NOR-X-417 (Edmonton, AB, 2010).
  * [27] Wiedinmyer, C., and Neff, J. C. Estimates of ${CO}^{2}$ from fires in the united states: implications for carbon management.  Carbon Balance Manag. 2, 1 (2007), 10.
  * [28] Yajima, T., and Nagahama, H. Finsler geometry for nonlinear path of fluids flow through inhomogeneous media.  Nonlinear Anal. Real World Appl. 25 (2015), 1–8.


# Fluctuations of a swarm of Brownian bees

Maor Siboni<EMAIL_ADDRESS>Racah Institute of Physics, Hebrew
University of Jerusalem, Jerusalem 91904, Israel  Pavel Sasorov
<EMAIL_ADDRESS>Institute of Physics CAS, ELI Beamlines, 182 21
Prague, Czech Republic  Keldysh Institute of Applied Mathematics, Moscow
125047, Russia  Baruch Meerson<EMAIL_ADDRESS>Racah Institute of
Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel

###### Abstract

The “Brownian bees” model describes an ensemble of $N$ independent branching
Brownian particles. When a particle branches into two particles, the particle
farthest from the origin is eliminated so as to keep the number of particles
constant. In the limit of $N\to\infty$, the spatial density of the particles
is governed by the solution of a free boundary problem for a reaction-
diffusion equation. At long times the particle density approaches a
spherically symmetric steady state solution with a compact support. Here we
study fluctuations of the “swarm of bees” due to the random character of the
branching Brownian motion in the limit of large but finite $N$. We consider a
one-dimensional setting and focus on two fluctuating quantities: the swarm
center of mass $X(t)$ and the swarm radius $\ell(t)$. Linearizing a pertinent
Langevin equation around the deterministic steady state solution, we calculate
the two-time covariances of $X(t)$ and $\ell(t)$. The variance of $X(t)$
directly follows from the covariance of $X(t)$, and it scales as $1/N$ as to
be expected from the law of large numbers. The variance of $\ell(t)$ behaves
differently: it exhibits an anomalous scaling $(1/N)\ln N$. This anomaly
appears because all spatial scales, including a narrow region near the edges
of the swarm where only a few particles are present, give a significant
contribution to the variance. We argue that the variance of $\ell(t)$ can be
obtained from the covariance of $\ell(t)$ by introducing a cutoff at the
microscopic time $1/N$ where the continuum Langevin description breaks down.
Our theoretical predictions are in good agreement with Monte-Carlo simulations
of the microscopic model. Generalizations to higher dimensions are briefly
discussed.

##  I Introduction

This work touches two distinct basic concepts of nonequilibrium statistical
physics: branching Brownian motion (BBM) and nonequilibrium steady states
(NESSs). BBM – a combination of Brownian motion of particles and their random
branching – was extensively studied in the past McKean ; Bramson , and it
continues to attract a great deal of attention from physicists BD2009 ;
Mueller ; Ramola1 ; DMS  and mathematicians.

In its turn, NESSs of ensembles of reacting and diffusing particles occupy an
important niche of nonequilibrium statistical mechanics JL1993 ; JL2004 ;
Bodineau2010 ; Hurtado2013 ; M2015 . The NESSs capture in a simple way
different aspects of physics of more complicated dissipative systems, both
living and non-living.

Here we study a simple but nontrivial example of a NESS in a particle-
conserving version of branching Brownian motion (BBM), which is known under
the name of “Brownian bees” bees1 ; bees2 ; MS2021 . The model involves $N$
independent particles (bees) located in a $d$-dimensional space. On a small
time interval $\Delta t$ each particle branches into two particles with
probability $\Delta t$. It also performs, with the complementary probability
$1-\Delta t$, continuous-time Brownian motion with diffusion constant $1$.
When a branching event occurs, the particle which is farthest from the origin
is immediately removed, so as to keep the number of particles constant at all
times. The NESS of this relatively simple system has some interesting and non-
trivial properties that will be the focus of this paper.

An additional motivation to study the Brownian bees model is that it is a
close relative of a whole family of Brunet-Derrida $N$-particle models:
branching Brownian models with selection, the studies of which were initiated
in Refs. BDMM2006 ; BDMM2007 . A Brunet-Derrida $N$-particle model consists of
a population of $N$ particles which undergo BBM. When a branching event
occurs, the particle of the lowest fitness is removed. The models differ among
themselves by their _fitness function_ – thus mimicking different aspects of
biological selection – and by the dimension of space. The original papers
BDMM2006 ; BDMM2007  considered discrete-time processes, but much of the
subsequent work BG2010 ; DR2011 ; BBS2013 ; M2016 ; BZ2018 ; DMFPSL ,
including the studies of the Brownian bees model bees1 ; bees2 ; MS2021  and
its barocentric variant ABLT ; barocentric , has focused on the continuous-
time branching Brownian motion.

Here we consider the Brownian bees model on the line, $d=1$. It has been
rigorously shown in Ref. bees1  that, in the limit of $N\to\infty$, the
coarse-grained spatial density $u(x,t)\geq 0$ of the particles, normalized by
$N$, is governed by a deterministic free boundary problem:

$\displaystyle\partial_{t}u(x,t)=\partial_{x}^{2}u(x,t)+u(x,t)\,,\quad|x|\leq\ell(t)\,,$
(1) $\displaystyle u(x,t)=0\,,\quad|x|>\ell(t)\,,$ (2)
$\displaystyle\int_{-\ell(t)}^{\ell(t)}u(x,t)\,dx=1\,.$ (3)

$u(x,t)$ is continuous at $|x|=\ell(t)$, and an initial condition should be
specified. As one can see, the compact support of $u(x,t)$, at all $t>0$, is
centered at the origin. Effectively, there are two absorbing walls, at
$x=\pm\ell(t)$, which move in synchrony so as to keep the number of particles
constant at all times.

Furthermore, it has been proven in Ref. bees2  that, at long times, the
solution of the deterministic problem (1)-(3) approaches a unique steady state

$U(x)=\left\\{\begin{array}[]{lr}{\displaystyle\frac{1}{2}}\cos
x\,,&~{}~{}~{}|x|\leq\ell_{0},\\\ \\\
0\,,&~{}~{}~{}|x|>\ell_{0},\end{array}\right.$ (4)

where $\ell_{0}=\pi/2$. Correspondingly, $\ell(t)$ approaches $\ell_{0}$.

A natural next step is to study _fluctuations_ in the Brownian bees model at
large but finite $N$, caused by the random character of the elemental
processes of the branching Brownian motion. _Persistent_ fluctuations,
including large deviations, of the swarm size $\ell(t)$ have been recently
considered in Ref. MS2021 . Here we study typical (that is small) Gaussian
fluctuations of the swarm center of mass $X(t)$ and the swarm radius
$\ell(t)$. Our goals in this work are to calculate the steady-state two-time
covariances of $X(t)$ and $\ell(t)$, which we call $g_{X}(\tau)$ and
$g_{\ell}(\tau)$, respectively. Here $\tau$ is the time between two
measurements of the corresponding quantities in the steady state. To this end
we introduce a pertinent Langevin equation for the Brownian bees model and
linearize it around the deterministic steady state solution (4). Solving the
linearized equation, we obtain explicit expressions for small random
deviations of $X(t)$ and $\ell(t)$ from their deterministic stationary values
$0$ and $\ell_{0}$, respectively. Then we use these solutions to compute the
functions $g_{X}(\tau)$ and $g_{\ell}(\tau)$.

We find that the covariance $g_{X}(\tau)$ remains finite at $\tau\to 0$ and
yields the variance of $X(t)$ which scales as $1/N$, as to be expected from
the law of large numbers.

The covariance $g_{\ell}(\tau)$ behaves markedly differently. At $\tau\ll 1$
it exhibits a logarithmic scaling with $\tau$, giving an unexpected example of
a system with $1/f$ noise. Furthermore, $g_{\ell}(\tau)$ formally diverges
logarithmically as $\tau\to 0$. As we argue, the continuum Langevin
description of the Brownian bees breaks down when $\tau$ becomes comparable
with the microscopic time of the model, $1/N$. Using the microscopic time as a
cutoff, we evaluate, with a logarithmic in $N$ accuracy, the variance which
scales as $(1/N)\ln N$. This and other theoretical predictions are in good
agreement with Monte Carlo simulations of the microscopic model that we
performed.

Here is a plan of the remainder of this paper. In Sec. II we introduce a
coarse-grained stochastic formulation of the problem, based on the Langevin
equation for the BBM. In Sec. III we formulate a linearized version of the
problem by expanding the solution around the (deterministic) steady-state
solution (4). In the same Sec. III we obtain the solution to the linearized
problem at long times. This solution is then used in Secs. IV and V to
calculate the two-time covariances, and then the variances, of $X(t)$ and
$\ell(t)$, respectively. The theoretical predictions are compared with results
of Monte-Carlo simulations. We summarize and briefly discuss our results in
Sec. VI. The Appendix includes a brief description of the Monte-Carlo
simulation algorithm.

##  II Langevin equation

To account for fluctuations resulting from the intrinsic randomness of the
elemental processes of the BBM, one must go beyond the deterministic equation
(1). At large $N$ typical fluctuations are small, and they are captured by the
Langevin equation: a stochastic partial differential equation vK ; Gardiner .
In our case the Langevin equation can be written as

$\partial_{t}u(x,t)=\partial_{x}^{2}u(x,t)+u(x,t)+\tilde{R}(u,x,t),\;\;|x|\leq\ell(t).$
(5)

This equation replaces the mean-field equation (1), whereas Eqs. (2) and (3)
remain unchanged. Now the solution $u(x,t)$ and the size of its support
$\ell(t)$ are random functions of their variables. The multiplicative noise
term $\tilde{R}(u,x,t)$ comes from two independent multiplicative noises: the
branching noise $\tilde{R}_{b}(u,x,t)$ and the noise of Brownian motion
$\tilde{R}_{d}(u,x,t)$:

$\displaystyle\tilde{R}(u,x,t)=\tilde{R}_{b}(u,x,t)+\tilde{R}_{d}(u,x,t)\,,$
$\displaystyle\tilde{R}_{b}(u,x,t)=\frac{\sqrt{u}}{\sqrt{N}}\eta(x,t)\,,$ (6)
$\displaystyle\tilde{R}_{d}(u,x,t)=\frac{1}{\sqrt{N}}\partial_{x}\left[\chi(x,t)\sqrt{2u}\right]\,,$
(7)

where $\eta(x,t)$ and $\chi(x,t)$ are two independent Gaussian white noises
with zero mean:

$\displaystyle\left\langle\eta(x_{1},t_{1})\eta(x_{2},t_{2})\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle\chi(x_{1},t_{1})\chi(x_{2},t_{2})\right\rangle$ (8)
$\displaystyle=$ $\displaystyle\delta(x_{1}-x_{2})\delta(t_{1}-t_{2})\,.$

A few words are in order about the origin of $R_{b}$ and $R_{d}$. The
branching noise $R_{b}$ can be derived from the exact master equation for the
branching process. For typical fluctuations in a system of many particles one
can use the van Kampen system size expansion to approximate the master
equation by a Fokker-Planck equation, see e.g. Ref. vK ; Gardiner . In the
equivalent language of the Langevin equation one obtains Eq. (6). The Brownian
motion noise term $R_{d}$ is best known in the context of a large-scale and
long-time description of the dynamics of a lattice gas of independent random
walkers, see e.g. Ref. Spohn .

Differentiating the conservation law (3) with respect to time and using Eq.
(5), we obtain

$\partial_{x}u[-\ell(t),t]-\partial_{x}u[\ell(t),t]=1+\int\limits_{-\ell(t)}^{\ell(t)}\tilde{R}\,dx\,.$
(9)

We will use this equation instead of Eq. (3) at $t>t_{0}$, where $t_{0}$ is
the initial time.

One should bear in mind that a Langevin description provides a _macroscopic_
description of fluctuations in microscopic models like the one we are dealing
with here. In this particular case it is accurate only at times much longer
than the typical time between two consecutive branching events: $\Delta t\gg
1/N$, and at distances much larger than the typical interparticle distance.
The latter is of order $1/N$ in the bulk of the swarm, but much larger – of
order $1/\sqrt{N}$ – at the edges of the swarm, where the macroscopic density,
see Eq. (4), approaches zero.

##  III Linearization and solution

###  III.1 Linearization

The small parameter $1/\sqrt{N}\ll 1$ in the noise term in Eq. (5) calls for a
perturbation expansion. In the leading order one can simply linearize Eq. (5)
around the deterministic steady state (4). We write

$\displaystyle u(x,t)=U(x)+v(x,t)\,,\quad|v|\ll 1\,;\quad
v\propto\frac{1}{\sqrt{N}}\,,$ (10)
$\displaystyle\ell(t)=\frac{\pi}{2}+\delta\ell(t)\,,\quad|\delta\ell|\ll
1\,;\quad\delta\ell\propto\frac{1}{\sqrt{N}}.$ (11)

The linearized versions of Eqs. (5), (2) and (9) are

$\displaystyle\partial_{t}v(x,t)-\partial_{x}^{2}v(x,t)-v(x,t)=\tilde{R}(x,t)\,,$
(12)
$\displaystyle-\frac{1}{2}\delta\ell(t)+v\left(\pm\frac{\pi}{2},t\right)=0\,,$
(13)
$\displaystyle\partial_{x}v\left(-\frac{\pi}{2},t\right)-\partial_{x}v\left(\frac{\pi}{2},t\right)=\int\limits_{-\pi/2}^{\pi/2}\tilde{R}(x,t)\,dx\,,$
(14)

where we have denoted $\tilde{R}(x,t)\equiv\tilde{R}(U(x),x,t)$. At the
initial time $t=t_{0}$ we have

$\int_{-\pi/2}^{\pi/2}v(x,t_{0})\,dx=0\,.$ (15)

In the steady-state regime, that we are interested in, the results will depend
neither on the initial condition $v(x,t_{0})$, nor on $t_{0}$ itself.
Therefore, we will set $v(x,t_{0})=0$ and ultimately send $t_{0}$ to
$-\infty$.

The two conditions in Eq. (13) can be rewritten as

$\displaystyle
v\left(-\frac{\pi}{2},t\right)=v\left(\frac{\pi}{2},t\right)\,,$ (16)
$\displaystyle\delta\ell(t)=2v\left(\frac{\pi}{2},t\right)\,.$ (17)

Equation (16) allows us to formally extend $v$ and $\tilde{R}$ as functions of
$x$ periodically to the whole $x$-axis. As follows from Eq. (14), the
derivative $\partial_{x}v$ experiences jumps at the points $x=\pi/2+m\pi$,
$m=0,\pm 1\,,\dots$. Each jump is equal to

$\partial_{x}v\left(\frac{\pi}{2}+0,t\right)-\partial_{x}v\left(\frac{\pi}{2}-0,t\right)=\int\limits_{-\pi/2}^{\pi/2}\tilde{R}(x,t)\,dx\,.$
(18)

We can explicitly account for these jumps by modifying the source term in Eq.
(12):

$\partial_{t}v(x,t)-\partial_{x}^{2}v(x,t)-v(x,t)=R(x,t)\,,$ (19)

where the new source term is $R=R_{b}+R_{d}$,

$\displaystyle R_{i}(x,t)$ $\displaystyle=$
$\displaystyle\\!\tilde{R}_{i}(x,t)$ (20) $\displaystyle-$
$\displaystyle\\!\\!\left(\int\limits_{-\pi/2}^{\pi/2}\tilde{R}_{i}(x,t)\,dx\\!\right)\\!\\!\sum_{m\in\mathbb{Z}}\delta\left(x-\frac{\pi}{2}-m\pi\right)$

and $i=b,d$. For this $R$ (and also separately for $R_{b}$ and $R_{d}$) the
condition

$\int\limits_{-\pi/2+\Delta}^{\pi/2+\Delta}R(x,t)\,dx=0$ (21)

is obeyed automatically for all $\Delta$. To avoid uncertainty with the
positions of the delta-function, we can shift the interval of interest
$\left[-\pi/2,\pi/2\right]$ by an infinitely small $\Delta$. Equation (14) now
reads:

$\partial_{x}v\left(-\frac{\pi}{2}+0,t\right)-\partial_{x}v\left(\frac{\pi}{2}+0,t\right)=0\,.$
(22)

The linearized problem is now completely defined by Eq. (19), the conditions
(16), (17) and (22), and a zero initial condition $v(x,t_{0})=0$.

###  III.2 Expanding over eigenfunctions

We can solve the linear equation (19) by expanding the solution over the
eigenfunctions $V_{n}(x)e^{-\lambda_{n}t}$ of the linear operator on the left
hand side of Eq. (19), subject to periodic boundary conditions (16) and (22).
The complete set of eigenfunctions and eigenvalues is the following:

$\left\\{1,\cos(2nx),\,\sin(2nx)\right\\}\,,\quad\quad\lambda_{n}=4n^{2}-1\,,$
(23)

where $n=1,2,3,\dots$. Let us introduce the Green’s function
$G(x,x^{\prime};t,t^{\prime})$ of the problem, that is the solution of Eq.
(19) with the source term $R=\delta(x-x^{\prime})\delta(t-t^{\prime})$,
subject to the boundary conditions (16) and (22) and a zero initial condition
at $t=t^{\prime}$. The Green’s function has the form

$\displaystyle G(x,x^{\prime};t,t^{\prime})$ $\displaystyle=$
$\displaystyle\left[\frac{e^{t-t^{\prime}}}{\pi}+\frac{2}{\pi}\sum\limits_{n=1}^{\infty}\cos
2n(x-x^{\prime})\,e^{-\lambda_{n}(t-t^{\prime})}\right]$ (24)
$\displaystyle\times$ $\displaystyle\theta(t-t^{\prime})\,,$

where $\theta(t-t^{\prime})$ is the step function. Using the Green’s function,
we can write an explicit expression for $v(x,t)$ and, by virtue of Eq. (17),
for $\delta\ell(t)$. Sending $t_{0}$ to $-\infty$, we arrive at the following
results:

$v(x,t)\\!=\\!\frac{2}{\pi}\sum\limits_{n=1}^{\infty}\int\limits_{0}^{\infty}e^{-\lambda_{n}t^{\prime}}\\!\\!\\!\\!\\!\\!\\!\\!\int\limits_{-\pi/2+0}^{\pi/2+0}\\!\\!\\!\\!\\!\\!R(x^{\prime},t-t^{\prime})\,\cos
2n(x-x^{\prime})\,dx^{\prime}dt^{\prime},$ (25)
$\delta\ell(t)\\!=\\!\frac{4}{\pi}\sum\limits_{n=1}^{\infty}(-1)^{n}\int\limits_{0}^{\infty}e^{-\lambda_{n}t^{\prime}}\\!\\!\\!\\!\\!\\!\\!\\!\int\limits_{-\pi/2+0}^{\pi/2+0}\\!\\!\\!\\!\\!\\!R(x,t-t^{\prime})\cos
2nx\,dx\,dt^{\prime}.$ (26)

We will use these expressions in the subsequent sections IV and V to analyze
the fluctuations of the swarm center of mass (COM) and of the swarm radius,
respectively.

Figure 1: Fluctuations of the COM $X(t)$. Top panel: the variance of $X(t)$
versus $N$ from Monte-Carlo simulations (symbols) and from theoretical
prediction (37) (solid line). Bottom panel: the covariance $g_{X}(\tau)$ vs.
$\tau$ at $N=2500$ from simulations (symbols) and from Eq. (32) (solid line).
The dashed line is the small-$\tau$ asymptotic
$g_{X}(\tau)=0.6346...-2.4674...\,\tau$.

##  IV Fluctuations of the center of mass

In the Langevin description, the coordinate of the COM of the system is

$X(t)=\int_{-\ell(t)}^{\ell(t)}xu(x,t)dx=\int_{-\pi/2}^{\pi/2}xv(x,t)\,dx\,,$
(27)

where, to the leading order, we can use our result (25) for $v(x,t)$. Let us
evaluate the COM’s covariance in the steady state,

$g_{X}(\tau)\equiv\textrm{cov(}X(t),X(t+\tau)\textrm{)}=\left<X(0)X(\tau)\right>\,.$
(28)

Combining Eqs. (25) and (27), we obtain

$\displaystyle X(t)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}f_{n}(t)\,,$
$\displaystyle f_{n}(t)\\!$ $\displaystyle=$
$\displaystyle\\!\int\limits_{0}^{\infty}e^{-\lambda_{n}t^{\prime}}\\!\\!\\!\\!\\!\\!\int\limits_{-\pi/2+0}^{\pi/2+0}\\!\\!\\!\\!R(x,t-t^{\prime})\,\sin(2nx)\,dx\,dt^{\prime}.$
(29)

Plugging $X(t)$ into Eq. (28) and performing the averaging over the noises
$\eta$ and $\chi$, we obtain after some algebra

$g_{X}(\tau)=\frac{1}{N}\sum_{n,\,m\geq 1}s^{(X)}_{nm}\,,$ (30)

where

$s^{(X)}_{nm}\\!=\\!\frac{-16(n^{2}+m^{2})\,\exp\left[-\left(4n^{2}-1\right)|\tau|\right]}{(4(n+m)^{2}-1)(4(n-m)^{2}-1)(2n^{2}+2m^{2}-1)}\,.$
(31)

Fortunately, the infinite sum over $m$ can be evaluated exactly using
“Mathematica”. This leaves us with a single infinite sum:

$g_{X}(\tau)=\frac{1}{N}\sum_{n=1}^{\infty}r_{n}^{(X)}\,\exp\left[-\left(4n^{2}-1\right)|\tau|\right]\,,$
(32)

where

$r_{n}^{(X)}=\frac{8\left(64n^{6}-32n^{4}+n^{2}\right)-2\pi\left(4n^{2}-1\right)^{2}\sqrt{4n^{2}-2}\,\coth\left(\pi\sqrt{n^{2}-\frac{1}{2}}\right)}{\left(4n^{2}-1\right)^{2}\left(128n^{6}-128n^{4}+34n^{2}-1\right)}\,.$
(33)

As $n\to\infty$, the coefficient $r_{n}^{(X)}$ behaves as $1/n^{4}$.
Therefore, the infinite sum in Eq. (32) converges for any $\tau$ including
$\tau=0$. At long times, $|\tau|\gg 1$, $g_{X}(\tau)$ decays exponentially
with the decay rate determined by the ground state eigenvalue $\lambda_{1}=3$:

$g_{X}(|\tau|\gg 1)\simeq\frac{1}{N}\,r_{1}^{(X)}e^{-3|\tau|}\,,$ (34)

where

$r_{1}^{(X)}=\frac{264-18\sqrt{2}\pi\coth\left(\frac{\pi}{\sqrt{2}}\right)}{297}=0.61321\dots\,.$
(35)

In its turn, the quantity $g_{X}(\tau=0)$ yields the variance of the COM
fluctuations:

$\left<X^{2}\right>=g_{X}(\tau=0)=\frac{1}{N}\sum_{n=1}^{\infty}r_{n}^{(X)}\,.$
(36)

Evaluating the infinite sum over $n$ numerically, we obtain

$\text{var}\,X(t)=\left<X^{2}\right>=\frac{0.6346\dots}{N}\,.$ (37)

Figure 1 (top panel) compares this theoretical prediction with the variance of
the COM, measured in Monte-Carlo simulations at different $N$ (see the
Appendix for a brief description of the simulation algorithm). Figure 1
(bottom panel) compares the theoretically predicted covariance $g_{X}(\tau)$,
Eq. (30), with the covariance measured in the simulations at $N=2500$ and
different $\tau$. An excellent agreement is observed in both cases.

The $1/N$ scaling of the variance of the COM is to be expected from the law of
large numbers. It also reflects the fact that large-scale hydrodynamic modes
(with wavelengths comparable to the swarm size) dominate the contribution to
the variance. We will encounter quite a different situation in the next
section, which deals with fluctuations of the swarm radius $\ell(t)$.

##  V Fluctuations of the swarm radius

Averaging Eq. (26) over the noise, we obtain $\left<\delta\ell(t)\right>=0$.
Therefore, at least to the accuracy of the linear theory in $1/\sqrt{N}$, we
have $\left<\ell(t)\right>=\pi/2$. The covariance of $\ell(t)$ is given by

$g_{\ell}(\tau)=\left<\ell(0)\ell(\tau)\right>-\left<\ell\right>^{2}=\left<\delta\ell(0)\delta\ell(\tau)\right>\,.$
(38)

It is convenient to calculate the contributions to $g_{\ell}(\tau)$ from the
branching noise and from the noise of Brownian motion separately. Therefore we
define $\delta\ell=\delta\ell_{b}+\delta\ell_{d}$, where

$\delta\ell_{i}(t)\\!=\\!\frac{4}{\pi}\sum\limits_{n=1}^{\infty}(-1)^{n}\\!\\!\int\limits_{0}^{\infty}e^{-\lambda_{n}t^{\prime}}\\!\\!\\!\\!\\!\\!\\!\\!\int\limits_{-\pi/2+0}^{\pi/2+0}\\!\\!\\!\\!\\!\\!R_{i}(x,t-t^{\prime})\cos
2nx\,dx\,dt^{\prime}$ (39)

and $i=b,d$. We can write

$\displaystyle g_{\ell}(\tau)$ $\displaystyle=$ $\displaystyle
g_{\ell}^{(b)}(\tau)+g_{\ell}^{(d)}(\tau)$ (40) $\displaystyle=$
$\displaystyle\left<\delta\ell_{b}(0)\,\delta\ell_{b}(\tau)\right>+\left<\delta\ell_{d}(0)\,\delta\ell_{d}(\tau)\right>\,.$

Then, using Eq. (39) and performing the averaging, we obtain

$g_{\ell}^{(i)}(\tau)=\frac{8}{\pi^{2}N}\sum_{n,\,m\geq
1}\frac{\exp\left[-\left(4n^{2}-1\right)|\tau|\right]}{4(n^{2}+m^{2})-2}s^{(i)}_{nm}\,,$
(41)

where

$\displaystyle s^{(b)}_{nm}$ $\displaystyle=$
$\displaystyle-\frac{1}{4(n+m)^{2}-1}-\frac{1}{4(n-m)^{2}-1}$ (42)
$\displaystyle+$ $\displaystyle\frac{2}{4n^{2}-1}+\frac{2}{4m^{2}-1}+2$

and

$s^{(d)}_{nm}=8nm\left[\frac{1}{4(n+m)^{2}-1}-\frac{1}{4(n-m)^{2}-1}\right]\,.$
(43)

Again, the infinite sums over $m$ can be evaluated exactly with “Mathematica”,
leading to the expression

$g_{\ell}^{(i)}(\tau)=\frac{64}{\pi
N}\sum_{n=1}^{\infty}r_{n}^{(\ell)}\exp\left[-\left(4n^{2}-1\right)|\tau|\right]\,,$
(44)

where

$r_{n}^{(\ell)}=\frac{n^{2}\sqrt{4n^{2}-2}\coth(\pi\sqrt{n^{2}-1/2})}{64n^{4}-32n^{2}+1}$
(45)

for both $i=b$ and $i=d$. Somewhat surprisingly, the branching noise and the
noise of the Brownian motion give identical contributions to the covariance of
$\ell(t)$. Therefore,

$g_{\ell}(\tau)=\frac{128}{\pi
N}\sum_{n=1}^{\infty}r_{n}^{(\ell)}\exp\left[-\left(4n^{2}-1\right)|\tau|\right]\,.$
(46)

For $\tau\neq 0$ the infinite sum over $n$ in Eq. (46) converges. In
particular, the long-time decay of $g_{\ell}(\tau)$ with $\tau$ is again
purely exponential, $\sim e^{-3|\tau|}$, with the same ground state eigenvalue
$\lambda_{1}=3$ as for the COM.

But what happens at $\tau=0$? The large-$n$ asymptotic of $r_{n}^{(\ell)}$
from Eq. (45) is

$r_{n}^{(\ell)}\simeq\frac{1}{32n}\,.$ (47)

Therefore, at $\tau=0$ the sum over $n$ in Eq. (46) diverges logarithmically
at $n\to\infty$, implying an infinite variance. On the other hand, our Monte
Carlo simulations (see below) clearly show a finite variance of $\ell$. How to
resolve this paradox?

Figure 2: Fluctuations of the swarm radius $\ell(t)$. Top panel: the variance
of $\ell(t)$ versus $N$ as measured in Monte Carlo simulations. Bottom panel:
$N$ times the variance of $\ell(t)$ vs. $\ln N$. Symbols indicate simulation
results. The dashed line shows the function $(2/\pi)\ln N+a$, where $a\simeq
1.29\pm 0.03$, see the text after Eq. (50).

Let us return to the covariance $g_{\ell}(\tau)$ and recall that the Langevin
description, which we adopted here as an approximation to the microscopic
model, is applicable only at “macroscopic” times $\tau\gg 1/N$. Since the
divergence of $g_{\ell}(\tau)$ as $\tau\to 0$ is only logarithmical, we can
introduce a cutoff at $\tau\sim 1/N$ in the hope to obtain correct result for
the variance of $\ell(t)$ with logarithmic accuracy. Therefore, we proceed in
the following way. For $1/N\ll|\tau|\ll 1$ we use the large-$n$ asymptotic
(47) of $r_{n}^{(\ell)}$ and replace the summation over $n$ in Eq. (46) by
integration. This gives

$g_{\ell}(\tau)\simeq\frac{2}{N}\frac{e^{|\tau|}}{\pi}\Gamma(0,4|\tau|)\,,$
(48)

where $\Gamma(0,4|\tau|)$ is the incomplete Gamma function. In fact, we must
use the $|\tau|\ll 1$ asymptotic of Eq. (48), which leads us to

$g_{\ell}(\tau)\simeq\frac{2}{\pi
N}\ln\frac{1}{|\tau|}\,,\quad\frac{1}{N}\ll|\tau|\ll 1\,.$ (49)

The logarithmic growth of correlations of $\ell(t)$ with a decrease of the
time difference $\tau$ is a noticeable phenomenon. In the frequency space it
corresponds to a $1/f$ noise. The $1/f$ noise has been observed in a plethora
of physical and biological systems 1fnoise , and its appearance in the
fluctuations of $\ell(t)$ is unexpected.

Figure 3: Fluctuations of the swarm radius $\ell(t)$. Symbols: the covariance
$g_{\ell}(\tau)$ measured in the Monte Carlo simulations for $N\\!=\\!2500$.
Solid line: theoretical prediction (46). Dashed line: the asymptotic (49).

Now we can evaluate the variance of $\ell(t)$ by introducing a cutoff in Eq.
(49) at $\tau=1/N$. We obtain, with logarithmic accuracy,

$\textrm{var }\ell(t)\simeq\frac{2}{\pi}\frac{\ln N}{N}\,.$ (50)

This result is markedly different from Eq. (37) by the presence of the large
logarithmic factor $\ln N$, and this result is fully supported by our Monte-
Carlo simulations of the microscopic model. Figures 2 and 3 show the
simulation results for the fluctuations of $\ell(t)$, and compare them with
theoretical predictions (46), (49) and (50). As one can see from the bottom
panel of Fig. 2, the slope of the straight line, predicted by Eq. (50), agrees
very well with the simulations. To eliminate the offset, we introduced a
single adjustable offset parameter $a$. It corresponds to a numerical factor
under the logarithm which is beyond the logarithmic accuracy of Eq. (50).

##  VI Summary and Discussion

We studied analytically and numerically fluctuations of a swarm of Brownian
bees in the limit of large but finite $N$. We focused on two fluctuating
quantities: the COM of the swarm $X(t)$ and the swarm radius $\ell(t)$. We
employed a first-principles Langevin equation for this system, linearized it
around the deterministic steady state solution, and calculated the two-time
covariances of $X(t)$ and $\ell(t)$. As we found, the variance of $X(t)$
behaves “normally”: it scales as $1/N$ as to be expected from the law of large
numbers. But the variance of the swarm radius $\ell(t)$ exhibits an anomalous
scaling $(1/N)\ln N$. This anomaly is caused by significant contributions to
the variance from all spatial scales: from scales comparable to the swarm size
down to the scale $\sim 1/\sqrt{N}$ of a narrow region near the edges of the
swarm where only a few bees are present. We have also shown that, at high but
not too high frequencies, the power spectrum of the fluctuations of $\ell(t)$
exhibits the $1/f$ noise. Our Monte-Carlo simulations of the microscopic
system are in good agreement with theory.

A generalization to higher dimensions is especially interesting for the
fluctuations of the swarm radius $\ell(t)$. One should expect that a small
fraction of particles near the (circular or spherical) swarm boundary give a
significant contribution to the variance of $\ell(t)$ in higher dimensions as
well. Furthermore, in all dimensions, the average steady-state particle flux
to the boundary is the same (and equal to $N$). Therefore, the deterministic
steady-state particle density $u=U_{d}(\mathbf{x})$ behaves in a universal way
near the swarm boundary. Does this universality lead to a similar behavior of
the variance of $\ell(t)$ vs. $N$ in different dimensions? Figure 4 depicts
the results of our Monte Carlo simulations for $d=2$ and $3$, where we
measured the variance of $\ell(t)$ at different $N$. For comparison, also
shown is the variance of $\ell(t)$ for $d=1$. Amazingly, as one can see from
the figure, the slope $2/\pi$ of the straight line
$N\times\text{var}\,\ell(t)$ vs. $\ln N$, predicted by Eq. (50) for $d=1$,
holds for $d=2$ and $3$ as well. This remarkable coincidence certainly
deserves a closer look at via an extension of our analytical calculations to
higher dimensions.

Figure 4: The variance of the swarm radius $\ell(t)$ for $d=1,2$ and $3$ as
measured in the Monte Carlo simulations. Shown is $N$ times the variance of
$\ell(t)$ vs. $\ln N$. The symbols indicate simulation results for $d=1$, $2$
and $3$ (from bottom to top, respectively). The three dashed lines show the
function $(2/\pi)\ln N+a_{d}$ with the same slope $2/\pi$ and $d$-dependent
offsets: $a_{1}=1.29\pm 0.03$, $a_{2}=2.42\pm 0.01$, and $a_{3}=3.35\pm 0.05$.

It is instructive to compare the Brownian bees model with stochastic models,
describing “pulled” reaction fronts propagating into an unstable empty state.
The best known example of such models is the stochastic Fisher-Kolmogorov-
Petrovsky-Piskunov (F-KPP) equation. The velocity of a stochastic F-KPP front
fluctuates around its mean value, and the variance of these fluctuations
scales as $1/\ln^{3}N$, where $N\gg 1$ is the typical number of particles in
the front region Derridaetal . This F-KPP anomaly is much stronger than the
one in the Brownian bees model. This is because the F-KPP front velocity
fluctuations are _solely_ determined by a small fraction of particles located
at the leading edge of the front Derridaetal .

In conclusion, we should remind the reader that the Langevin equation provides
an accurate description only of typical, small fluctuations of reacting and
diffusing systems of particles. Therefore, our present work is limited to
typical fluctuations of the Brownian bees. It would be very interesting, and
challenging, to also study _large deviations_ of $X(t)$ and $\ell(t)$, which
correspond to the distribution tails of these quantities. Large deviations of
_persistent_ fluctuations of $\ell(t)$ have been recently addressed in Ref.
MS2021  by employing the optimal fluctuation method which bypasses the
Langevin description.

## Acknowledgments

We are grateful to N. R. Smith for useful discussions. The research of M.S.
and B.M. was supported by the Israel Science Foundation (grant No. 1499/20).
The research of P.S. was supported by the project “High Field Initiative”
(CZ.02.1.01/0.0/0.0/15_003/0000449) of the European Regional Development Fund.

##  Appendix A Monte-Carlo simulations

We performed continuous-time Monte-Carlo simulations (see e.g. Ref. MCbook )
of the microscopic model. The small time intervals $\Delta t$ between two
consecutive binary branching events were drawn from the exponential
distribution with rate parameter $N$. During these time intervals $\Delta t$
the particle positions were advanced according to the Brownian motion, given
by the normal distribution $N(x_{i},2\Delta t)$, where $x_{i}$ is the
particle’s initial position, and $i=1,2,\dots,N$.

Next, one particle is chosen randomly, from the uniform distribution, to be
the particle which is about to branch. The particle farthest from the origin
is then relocated to the position of the chosen particle. This is equivalent
to the removal of the farthest particle because of the birth of a “new”
particle. This process can be thought of as relabeling the particles.

The measurements were performed, after the system reached a steady state, at
fixed time steps $\delta t=1/N$. The simulation time was chosen to be
sufficiently large, so that the precision of all measurements is at least
within $1\%$ with $95\%$ confidence.

## References

  * (1) H.P. McKean, Comm. Pure Appl. Math. 28, 323 (1975).
  * (2) M.D. Bramson, Mem. Am. Math. Soc. vol. 44 (285) (1983).
  * (3) É. Brunet and B. Derrida, EPL 87, 60010 (2009); J. Stat. Phys. 143, 420 (2011).
  * (4) A. H. Mueller and S. Munier, Phys. Rev. E 90, 042143 (2014).
  * (5) K. Ramola, S. N. Majumdar, and G. Schehr, Phys. Rev. Lett. 112, 210602 (2014); Chaos, Solitons and Fractals 74, 79 (2015); Phys. Rev. E 91, 042131 (2015).
  * (6) B. Derrida, B. Meerson and P. V. Sasorov, Phys. Rev. E 93, 042139 (2016).
  * (7) G. Jona-Lasinio, C. Landim, and M. E. Vares, Probab. Theory Relat. Fields 97, 339 (1993).
  * (8) G. Basile and G. Jona-Lasinio, Int. J. Mod. Phys. B 18, 479 (2004).
  * (9) T. Bodineau and M. Lagouge, J. Stat. Phys. 139, 201 (2010).
  * (10) P. I. Hurtado, A. Lasanta, and A. Prados, Phys. Rev. E 88, 022110 (2013).
  * (11) B. Meerson, J. Stat. Mech. 2015, 05004.
  * (12) J. Berestycki, É. Brunet, J. Nolen, and S. Penington, arXiv:2006.06486.
  * (13) J. Berestycki, É. Brunet, J. Nolen, and S. Penington, arXiv:2005.09384.
  * (14) B. Meerson and P. Sasorov, Phys. Rev. E 103, 032140 (2021).
  * (15) É. Brunet, B. Derrida, A. H. Mueller and S. Munier, Europhys. Lett., 76, 1 (2006).
  * (16) É. Brunet, B. Derrida, A. H. Mueller and S. Munier, Phys. Rev. E 76, 041104 (2007).
  * (17) J. Bérard and J.-B. Gouéré, Comm. Math. Phys. 298, 323 (2010).
  * (18) R. Durrett and D. Remenik, Ann. Probab. 39, 2043 (2011).
  * (19) J. Berestycki, N. Berestycki, and J. Schweinsberg, Ann. Probab. 41, 527 (2013).
  * (20) P. Maillard, Probab. Theory Related Fields 166, 1061 (2016).
  * (21) N. Berestycki and L. Z. Zhao, Ann. Appl. Probab. 28, 651 (2018).
  * (22) A. De Masi, P. A. Ferrari, E. Presutti and N. Soprano-Loto, in “Stochastic Dynamics out of Equilibrium”, Springer Proc. Math. Stat., vol. 282, p. 523 (Springer, Cham, 2019).
  * (23) L. Addario-Berry, J. Lin, and T. Tendron, arXiv:2006.04743.
  * (24) In the variant of the Brownian bees model, considered in Ref. ABLT , one removes, upon branching, the particle farthest from the _center of mass_ of the system.
  * (25) N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 2001).
  * (26) C. W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 2004).
  * (27) H. Spohn, Large Scale Dynamics of Interacting Particles (Springer-Verlag, New York, 1991).
  * (28) P. H. Handel and A. L. Chung, Noise in Physical Systems and $1/f$ Fluctuations (AIP, New York, 1993).
  * (29) É. Brunet, B. Derrida, A. H. Mueller, and S. Munier, Phys. Rev. E 73, 056126 (2006).
  * (30) D. P. Kroese, T. Taimre and Z. I. Botev, Handbook of Monte Carlo Methods, Wiley Series in Probability and Statistics (John Wiley and Sons, New York, 2011).