# Can a Fractional Order Delay Differential Equation be Chaotic Whose Integer-
Order Counterpart is Stable?

Sachin Bhalekar111Corresponding Author Email<EMAIL_ADDRESS>Deepa Gupta  School of Mathematics and Statistics, University of Hyderabad,
Hyderabad, 500046 India

###### Abstract

For the fractional order systems

$D^{\alpha}x(t)=f(x),\quad 0<\alpha\leq 1,$

one can have a critical value of $\alpha$ viz $\alpha_{*}$ such that the
system is stable for $0<\alpha<\alpha_{*}$ and unstable for
$\alpha_{*}<\alpha\leq 1$. In general, if such system is stable for some
$\alpha_{0}\in(0,1)$ then it remains stable for all $\alpha<\alpha_{0}.$
In this paper, we show that there are some delay differential equations

$D^{\alpha}x(t)=f(x(t),x(t-\tau))$

of the fractional order which behave in an exactly opposite way. These systems
are unstable for higher values of fractional order and stable for the lower
values.
The striking observation is the example which is chaotic for $\alpha=0.27$ but
stable for $\alpha=1$. This cannot be observed in the fractional differential
equations (FDEs) without delay.
We provide the complete bifurcation scenarios in the scalar FDEs.

##  1 Introduction

Various equations arising in the theory of Applied Mathematics played a vital
role in modeling natural systems [1, 2, 3, 4, 5, 6]. The time variable is
continuous in most of these systems and one can use the ordinary or partial
differential equations to model the dynamics. The rate of change of the
physical quantity under consideration is described by the ``derivatives'' in
these models. The remarkable results in the history of Science viz. Newtons's
laws [7, 8], Einstein's theory [9, 10], heat conduction [11], chemical
kinetics [12] etc made use of the derivative operators.
In spite of a wide range of applicability, the derivative operator fails in
modeling the ``memory and hereditary properties'' in the natural systems
because it is a ``local operator''. Vito Volterra in 1959 introduced the
theory of non-local operators [13] which proved advantageous in modeling
memory and hereditary properties.
The fractional derivative(FD) [14, 15] and the delay [16, 17] are the
popularly used mathematical concepts which are very useful in achieving this
goal. The fractional differential equation can be viewed as an integral
equation with the kernel $(t-\tau)^{\alpha-1}/\Gamma(\alpha)$ (or similar).
The flexible order $\alpha$ of the derivative can be adjusted according to the
physical system to achieve more accuracy in the modeling. On the other hand,
the delay differential equation(DDE) contains the value of state variable at
the past time to include the memory of the system [18, 19].
Fractional differential equations have a wide applications in the studies of
viscoelastic materials [20, 21] as well as in many fields of science and
engineering including fluid flow, rheology, electrical networks and
electromagnetic theory [22, 15, 23]. Even a non-linear oscillation of an
earthquake can be modeled with the fractional derivatives [24]. These
operators are also proved useful in modeling various phenomena in
bioenginering and related area [25, 26].
The delay appears in the natural system because there is always a time
duration of certain hidden processes (e.g. duration of infectious period)
[27]. The delay differential equations have applications in various areas of
life sciences such as in population dynamics, epidemiology, immunology,
physiology and neural networks [18, 28, 29]. So, the fractional order delay
differential equations (FDDE) containing FD as well as the delay seem to be
more realistic.
Although the fractional derivative and delay are more practical, the problems
become infinite dimensional [17, 30]. The mathematical analysis is more
involved and one has to be careful while applying the classical results.
The stability and bifurcation analyses are the principal branches of Dynamical
Systems [31, 32]. One has to provide these results before applying these
useful tools (viz. FD and delay) in the real-life applications. In the seminal
work, Bhalekar [33] provided the stability and bifurcation analyses of the
FDDE system

$D^{\alpha}x(t)=ax(t)+bx(t-\tau).$

He provided three regions viz. unstable(U), stable(S) and the delay dependent
stable(DS) regions (cf. Figure 1).
In the delay dependent stable region, the system is stable for some values of
delay $\tau$ and become unstable for other values.
In this work, we extend the analysis in [33] and provide more details on the
bifurcation regions in the delay dependent stable region.
The manuscript is organized as below: Section 2 provides the basic definitions
and results useful in understanding this paper. Section 3 analyses the
stability and bifurcations. Section 4 contains illustrative examples. In the
Section 5, we provide some new observations. The conclusions are given in
Section 6.

##  2 Preliminaries

In this section, we provide some basic definitions described in the literature
[14, 16, 15, 34, 35].

######  Definition 1 (Fractional Integral).

For any $f\in\mathcal{L}^{1}(0,b)$ the Riemann-Liouville fractional integral
of order $\upmu>0$, is given by

$\textit{I}^{\upmu}f(t)=\dfrac{1}{\Gamma(\upmu)}\int_{0}^{t}(t-\tau)^{\upmu-1}f(\tau)d\tau,\quad
0<t<b.$

######  Definition 2 (Caputo Fractional Derivative).

For $f\in\mathcal{L}^{1}(0,b)$, $0<t<b$ and $m-1<\upmu\leq m$,
$m\in\mathbb{N}$, the Caputo fractional derivative of function $f$ of order
$\upmu$ is defined by,

$\textit{D}^{\upmu}f(t)=\begin{cases}\frac{d^{m}}{dt^{m}}f(t),\textit{ if
}\quad\upmu=m\\\ \textit{I}^{m-\upmu}\dfrac{d^{m}f(t)}{dt^{m}},\textit{ if
}\quad m-1<\upmu<m.\end{cases}$

Note that for $m-1<\upmu\leq m$, $m\in\mathbb{N},$

$\textit{I}^{\upmu}\textit{D}^{\upmu}f(t)=f(t)-\sum_{k=0}^{m-1}\dfrac{d^{k}f(0)}{dt^{k}}\dfrac{t^{k}}{k!}.$

######  Definition 3 (Equilibrium Point).

Consider the generalized delay differential equation

$D^{\alpha}{x}(t)=f(x(t),x(t-\tau)),\quad 0<\alpha\leq 1,$ (1)

where $\tau\textgreater 0$, $f:E\rightarrow\mathbb{R}$,
$E\subseteq\mathbb{R}^{2}$ is open and $f\in C^{1}(E)$.

A steady state solution of equation (1) is called an equilibrium point.

Note that $x^{*}$ is an equilibrium point if and only if

$f(x^{*},x^{*})=0.$ (2)

The linearization of equation (1) near $x^{*}$ is given by [33].

$D^{\alpha}x(t)=ax(t)+bx(t-\tau),$ (3)

where $a=\partial_{1}f(x^{*},x^{*})$ and $b=\partial_{2}f(x^{*},x^{*})$ are
partial derivatives of $f$ with respect to the first and second variables
respectively evaluated at $x(t)=x(t-\tau)=x^{*}.$ Near the equilibrium point,
the qualitative behaviour of the trajectories of equations (3) and (1) are
same.

##  3 Bifurcation Analysis

######  Theorem 1.

[33] Consider the fractional order linear delay differential equation (3).
Case 1 If b $\in(-\infty,-|a|)$ then the stability region of $x^{*}$ in
$(\tau,a,b)$ parameter space is located between the plane $\tau=0$ and

$\tau(a,b,\alpha)=\dfrac{\arccos\Bigg{(}\dfrac{\Bigg{(}acos\Big{(}\dfrac{\alpha\pi}{2}\Big{)}+\sqrt{b^{2}-a^{2}\sin^{2}\Big{(}\dfrac{\alpha\pi}{2}}\Big{)}\Bigg{)}\cos\dfrac{\alpha\pi}{2}-a}{b}\Bigg{)}}{\Bigg{(}a\cos\Big{(}\dfrac{\alpha\pi}{2}\Big{)}+\sqrt{b^{2}-a^{2}\sin^{2}\Big{(}\dfrac{\alpha\pi}{2}\Big{)}}\Bigg{)}^{1/\alpha}}.$
(4)

The equation undergoes Hopf bifurcation at this value.
Case 2 If $b\in(-a,\infty)$ then $x^{*}$ is unstable for any $\tau\geq 0.$
Case 3 If $b\in(a,-a)$ and $a<0$ then $x^{*}$ is stable for any $\tau\geq 0.$
The curve $\tau(a,b,\alpha)$ given by equation (4) is called the boundary
curve.

This theorem provides three regions in $ab-$plane (cf. Figure 1 ).
viz. S- stable region-Case 3.
U- unstable region-Case 2 and
DS- delay dependent stable region-Case 1.

Figure 1:  Stability regions of equation (3)

In this work, we provide more details of the region DS in Figure 1.
The bifurcation regions are shown in Figure 2 and can be described as below:
1. U: If $a+b>0$ then system (3) is unstable $\forall\tau\geq 0.$
2. S: If $a<0$ and $a<b<-a$ then the system (3) is asymptotically stable $\forall\tau\geq 0.$
3. DS1: In this region, we get the stable region in $\alpha\tau-$ plane as shown in Figure 3.
4. DS2: The stability properties in this region are as shown in Figure 4.
5. DS3: The stable region in $\alpha\tau$-plane for this case is shown in Figure 5.

Figure 2: The complete bifurcation diagram of equation (3)

###  3.1 Bifurcation Scenarios

In this section, we provide the details of various stable regions described in
previous section.
Region DS1: It is observed that, for given $a\in\mathds{R}$, $\exists$ a
function $a_{0}(a)$ such that $-\infty<b<a_{0}(a)\Rightarrow$ the boundary
curve $\tau(a,b,\alpha)$ in the $\alpha\tau-$ plane is monotonically
increasing as given in Figure 3. Let $\tau_{*}=\tau(a,b,1)$. For fixed
$\tau_{0}\in(0,\tau_{*}]$, $\exists$ critical value of $\alpha$, say
$\alpha_{0}$ such that equilibrium point is asymptotically stable for
$\alpha_{0}<\alpha\leqslant 1$ and unstable for $0<\alpha<\alpha_{0}$. If
$\tau>\tau_{*}$ then the system is unstable for all $\alpha\in(0,1]$.
Note This case is exactly opposite of non-delayed case.

####  3.1.1 Estimation for the function $a_{0}(a)$:

We observed that for some values of $a\in\mathds{R}$, there exists functions
$a_{0}(a)$ and $a_{1}(a)$ such that the curve $\tau(a,b,\alpha)$ is
monotonically increasing for $b\in(-\infty,a_{0}(a))$ and in the interval
$a_{0}(a)<b<a_{1}(a)$, the curve has local maxima. This local maxima occurs at
$\alpha=1$ if $b\shortrightarrow a_{0}$.
Therefore, $a_{0}$ is the function such that-

$\lim_{b\shortrightarrow
a_{0}}\frac{\partial\tau(a,b,\alpha)}{\partial\alpha}\mid_{\alpha=1}=0.$

Solving this equation, we get,

$-\arccos(\frac{-a}{a_{0}})=\frac{\sqrt{a_{0}^{2}-a^{2}}(\frac{\pi}{2})}{a\frac{\pi}{2}+\sqrt{a_{0}^{2}-a^{2}}\log(\sqrt{a_{0}^{2}-a^{2}})}$

or,

$\frac{-a}{a_{0}}=\cos\Big{[}\frac{-\sqrt{a_{0}^{2}-a^{2}}(\frac{\pi}{2})}{a\frac{\pi}{2}+\sqrt{a_{0}^{2}-a^{2}}log(\sqrt{a_{0}^{2}-a^{2}})}\Big{]}.$
(5)

Equation (5) is the implicit expression for the function $a_{0}(a)$.

Figure 3:  Boundary curve (4) in the region DS1

Region DS2:
Case 1 We assume that $a<0$.
If $a_{0}<b<a_{1}$ then the stable curve $\tau(a,b,\alpha)$ is not monotonic
(cf. Figure(4)) and $\exists$ a local maxima at some $\alpha_{**}<1$. In this
case there exists $\tau_{*}=\tau(a,b,1)$ and $\tau_{**}=\tau(a,b,\alpha_{**})$
such that-
$\bullet$ If $0<\tau<\tau_{*}$ then-
(i) there exists critical value of $\alpha$, say $\alpha_{0}$ such that
$0<\alpha<\alpha_{0}$, we get unstable solution and for $\alpha_{0}<\alpha\leq
1$ we get stable solution.
$\bullet$ If $\tau_{*}<\tau<\tau_{**}$ then-
(i) there exist $\alpha_{1}$ and $\alpha_{2}$ such that for
$0<\alpha<\alpha_{1}$ and $\alpha_{2}<\alpha\leq 1$ we get unstable solution
and for $\alpha_{1}<\alpha<\alpha_{2}$ we get stable solution (stable window).
$\bullet$ If $\tau>\tau_{**}$ then equilibrium is unstable for all
$\alpha\in(0,1].$

####  3.1.2 Estimation of function $a_{1}(a):$

Let $a<0$. As we take the values of $b>a_{0}$, the left end of critical curve
i.e. $\tau(a,b,0)$ is zero upto some value $a_{1}(a)$. For $b>a_{1}(a)$,
$\lim_{\alpha\shortrightarrow 0}\tau(a,b,\alpha)=\infty$. It is observed that,
at this particular value $b=a_{1}(a),$

$\tau(a,a_{1}(a),0)=\pi.$

Solving this equation, we get ${a_{1}(a)=a-1}$.
Case 2 We assume that $a>0.$
It is observed that the stable region in $\alpha\tau$-plane is shown in Figure
4 for $a_{0}(a)<b<-a.$
Note that, $\nexists$ $a_{1}(a)$ if $a>0.$

Figure 4:  Boundary curve (4) in the region DS2

Region DS3: This region exists only for $a<0$.
If $a_{1}(a)<b<a<0$ then stable curve $\tau(a,b,\alpha)$ is monotonically
decreasing in $\alpha\tau-$plane as shown in Figure 5.
Note that $\lim_{\alpha\rightarrow 0}\tau(a,b,\alpha)=\infty,$ in this case.
For any $\tau_{0}>0$, $\exists$ $\alpha_{0}\in(0,1)$ satisfying
$\tau(a,b,\alpha_{0})=\tau_{0}$ such that $0<\alpha<\alpha_{0}\Rightarrow$
stable and $\alpha_{0}<\alpha\leq 1\Rightarrow$ unstable solutions.
This case is similar to the non-delayed FDEs.

Figure 5:  Boundary curve (4) in the region DS3

##  4 Illustrative Examples

Example 1 Consider the fractional order linear DDE

$D^{\alpha}x(t)=ax(t)+bx(t-\tau).$ (6)

In this case, a=-1.5. Therefore, $a_{0}=-4.31582$ and $b\in(-\infty,a_{0}).$
This parameter set is in the region DS1 of Figure 2. Therefore, the stable
region for this case is given by Figure 3.
The stable solution for $(\alpha,\tau)=(0.8,0.15)$ is shown in Figure 6. If we
take $(\alpha,\tau)=(0.4,0.3)$ then it belongs to the unstable region and this
is verified in Figure 7.
In the similar way we have also verified this by taking $a$ to be positive
value. In this case we take
$f(x(t),x(t-\tau))=2x(t)-x(t)^{2}-x(t-\tau)^{3}-4x(t-\tau)$.
At the equilibrium point zero we get $a=2$, $b=-4$ and $a_{0}=-2.33343$. By
fixing $\tau=0.12$ we get critical value $\alpha=0.583977$. So, we get stable
solution when $\alpha=0.7$ and for $\alpha=0.4$ we get unstable solution.

Figure 6:  Stable solution of (6) for $\alpha=0.8$ and $\tau=0.15$ Figure 7:
Unbounded solution of (6) for $\alpha=0.4$ and $\tau=0.3$

Example 2 Let us consider the non-linear system

$D^{\alpha}x(t)=-10.9x(t-\tau)-x(t-\tau)^{3}-x(t)^{2}-9x(t),$ (7)

where $0<\alpha\leq 1.$
Here, $f(x,y)=-9x-x^{2}-10.9y-y^{3}$. At equilibrium $x^{*}=0$, we have,
$a=\partial_{1}f(0,0)=-9$ and $b=\partial_{2}f(0,0)=-10.9$. For this $a$, we
have $a_{0}=-11.5972$ and $a_{1}=-10$. Therefore, $b\in(a_{0},a_{1})$.
Therefore, the local behaviour of $x^{*}$ is according to the bifurcation
region DS2. The stable region will be as shown in Figure 4 with
$\tau_{*}=0.413437$, $\tau_{**}=0.524926$ and $\alpha_{**}=0.653835.$
(1) If we fix $\tau=0.38<\tau_{*}$, then we get the critical value
$\alpha_{0}=0.384137$. So, we get stable solution for $\alpha=0.45$ (cf.
Figure 8 ) and unstable solution for $\alpha=0.3$ (cf. Figure 9).

Figure 8: Stable solution of equation (7) for $\alpha=0.45$ Figure 9: The
divergent solution of equation (7) for $\alpha=0.3$ and $\tau=0.38$

(2) When $\tau=0.45$ which is in between $\tau_{*}$ and $\tau_{**}$, we get
$\alpha_{1}=0.454118$ and $\alpha_{2}=0.919559$. Unbounded solutions are shown
in Figures 10 and 11 when $\alpha=0.41$ and $\alpha=0.95$ respectively and
bounded solution shown in Figure 12 for $\alpha=0.7.$
(3) When $\tau=0.54>\tau_{**}$ and $\alpha=0.7$ we get unstable solution (cf.
Figure 13).

Figure 10: Unstable solution of (7) for $\alpha=0.41<\alpha_{1}$ and
$\tau=0.45$ Figure 11: Unstable solution of (7) for $\alpha=0.95>\alpha_{2}$
and $\tau=0.45$ Figure 12: Stable solution of (7) for $\alpha=0.7$ which is in
between $\alpha_{1}$ and $\alpha_{2}$ for $\tau=0.45$ Figure 13: Unstable
solution of (7) for $\tau=0.54>\tau_{**}$ and $\alpha=0.7$

We also verify the results in the region DS2 with $a>0$. In this case
$f(x,y)=9x-9.03y$. So, near the equilibrium point $0$, $a=9$, $b=-9.03$,
$a_{0}=-9.06263$, $\tau_{*}=0.110865$ and $\tau_{**}=0.111041$. Therefore,
when $\tau=0.1$ we get critical value $\alpha=0.906205$. So, we get stable
solution for $\alpha=0.94$ and unstable for $\alpha=0.8$. Moreover, when
$\tau=0.111$ the corresponding $\alpha_{1}=0.986534$ and $\alpha_{2}=0.99529$.
So, we get divergent solutions for $\alpha=0.8$, $1$ and convergent solution
for $\alpha=0.9$. Also, when we go beyond $\tau_{**}$ say $0.12$ we get
unbounded solution for $\alpha=0.9.$

Example 3 For the region DS3, we take a non linear FDDE (1) with
$f(x(t),x(t-\tau))=-x(t)-x(t)^{2}-1.5x(t-\tau)-x(t-\tau)^{3}.$
So, near the equilibrium point $0$ we have $a=-1$ and $b=-1.5$ and $a_{1}=-2$.
Therefore, $b\in(a_{1},a)$. For the fixed $\tau_{0}=2.5$, we get the critical
value $\alpha_{0}=0.904463$. Therefore, for the smaller value $\alpha=0.4,$ we
obtain stable curve (cf. Figure 14) and for the higher values of $\alpha$,
e.g. $0.91$ we get divergent solution (cf. Figure (15)).

Figure 14: For the Example 3 when $\alpha=0.4$ and $\tau=2.5$ we get
convergent solution near equilibrium $0$ Figure 15:  When $\alpha=0.91$ and
$\tau=2.5$ we get divergent solution

##  5 New Observations

In this section, we provide the answer to the question in the title of this
paper. For the FDEs without delay, the stability of the system at higher value
of fractional order $\alpha$ does not get disturbed for the smaller values of
the order. As we have discussed in section 3, in the region DS1 of Figure 2,
the fractional order system may be stable for higher orders and unstable for
the lower order values. In the following example, the parameters $a$ and $b$
of the linearized equation near the equilibrium $x^{*}=2$ belongs to the
region DS1. We show that, the system is stable for $\alpha=1$ and show chaos
for $\alpha=0.27$. To the best of authors' knowledge, this is first time
reported in the literature. It can be observed that the chaotic attractor in
this example is very similar to that in the discrete dynamical system viz.
Henon map [36].
Example Consider the nonlinear FDDE

$D^{\alpha}x(t)=x(t)-x^{2}(t)+5x(t-\tau)-x(t-\tau)^{3},\quad 0<\alpha\leq 1.$
(8)

It is clear that $x^{*}=2$ is an equilibrium point.
The linearized equation at $x^{*}$ is

$D^{\alpha}x(t)=-3x(t)-7x(t-\tau)$

.
Therefore, $a=-3$ and $b=-7$.
Here, $a_{0}=-4.10117.$
Therefore, $b<a_{0}$ and the stability region for $x^{*}$ is DS1 (as shown in
Figure 2).
If we fix $\tau=0.31$ then the corresponding critical value of fractional
order is $\alpha_{*}=0.93777.$
Therefore, For $\alpha=1$, we get the stable solution as shown in Figure 16.
Further, for $\alpha=0.27<\alpha_{*}$, we get a chaotic solution as shown in
Figure 17. The chaotic attractor in this case is as shown in Figure 18. The
maximum Lyapunov exponent [37] in this case is $9.885276.$

Figure 16: The stable solution of (8) with $\alpha=1$ (integer order) Figure
17: Time series of (8) for $\alpha=0.27$ Figure 18: The chaotic attractor of
(8) at $\alpha=0.27$

##  6 Conclusions

In this paper, we extended the analysis in [33] and revealed the various
bifurcation regions in the delay dependent stable region.
We showed that, there are three subregions viz. DS1, DS2 and DS3. The region
DS3 has analogous bifurcation properties as the non-delayed FDEs whereas the
region DS1 shows the exactly contrast properties. The region DS2 has richer
dynamics and show various bifurcation scenarios including the stable window.
The remarkable acheivement is the example which shows the stable solutions for
integer order $\alpha=1$ and chaos for smaller value $\alpha=0.27$.
We hope that this manuscript will open new avenues in the research in
fractional order delay differential equations.

##  7 Acknowledgments

S. Bhalekar acknowledges the University of Hyderabad for Institute of
Eminence-Professional Development Fund (IoE-PDF) by MHRD (F11/9/2019-U3(A)).
D. Gupta thanks University Grants Commission for financial support
(No.F.82-44/2020(SA-III)).

## References

  * [1] M. Rafikov, J. M. Balthazar, H. Von Bremen, Mathematical modeling and control of population systems: applications in biological pest control, Applied Mathematics and Computation 200 (2) (2008) 557–573.
  * [2] R. Haberman, Mathematical models: mechanical vibrations, population dynamics, and traffic flow, SIAM, 1998.
  * [3] J. D. Logan, Applied mathematics, John Wiley & Sons, 2013.
  * [4] P. M. Morse, H. Feshbach, Methods of theoretical physics, American Journal of Physics 22 (6) (1954) 410–413.
  * [5] T. Manninen, M.-L. Linne, K. Ruohonen, Developing itô stochastic differential equation models for neuronal signal transduction pathways, Computational Biology and Chemistry 30 (4) (2006) 280–291.
  * [6] M. Braun, M. Golubitsky, Differential equations and their applications, Vol. 1, Springer, 1983.
  * [7] C. T. O’Sullivan, Newton’s law of cooling—a critical assessment, American Journal of Physics 58 (10) (1990) 956–960.
  * [8] I. Newton, Philosophiae naturalis principia mathematica, Vol. 1, G. Brookman, 1833\.
  * [9] K. Eriksson, K. Eriksson, D. Estep, P. Hansbo, C. Johnson, Computational differential equations, Cambridge University Press, 1996.
  * [10] A. Einstein, The photoelectric effect, Ann. Phys 17 (132) (1905) 4.
  * [11] J. Douglas, H. H. Rachford, On the numerical solution of heat conduction problems in two and three space variables, Transactions of the American mathematical Society 82 (2) (1956) 421–439.
  * [12] J. I. Steinfeld, J. S. Francisco, W. L. Hase, Chemical kinetics and dynamics, Prentice Hall Upper Saddle River, NJ, 1999.
  * [13] V. Volterra, E. T. Whittaker, Theory of functionals and of integral and integro-differential equations, Dover publications, 1959.
  * [14] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Elsevier, 1998.
  * [15] A. Kilbas, Theory and applications of fractional differential equations.
  * [16] M. Lakshmanan, D. V. Senthilkumar, Dynamics of nonlinear time-delay systems, Springer Science & Business Media, 2011.
  * [17] H. L. Smith, An introduction to delay differential equations with applications to the life sciences, Vol. 57, Springer New York, 2011.
  * [18] G. A. Bocharov, F. A. Rihan, Numerical modelling in biosciences using delay differential equations, Journal of Computational and Applied Mathematics 125 (1-2) (2000) 183–199.
  * [19] K. Gopalsamy, Stability and oscillations in delay differential equations of population dynamics, Vol. 74, Springer Science & Business Media, 2013.
  * [20] F. Mainardi, Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models, World Scientific, 2010.
  * [21] P. J. Torvik, R. L. Bagley, On the appearance of the fractional derivative in the behavior of real materials, Journal of Applied Mechanics 51 (2) (1984) 298–298.
  * [22] L. Debnath, Recent applications of fractional calculus to science and engineering, International Journal of Mathematics and Mathematical Sciences 2003 (54) (2003) 3413–3442.
  * [23] M. Saqib, I. Khan, S. Shafie, Application of fractional differential equations to heat transfer in hybrid nanofluid: modeling and solution via integral transforms, Advances in Difference Equations 2019 (1) (2019) 1–18.
  * [24] C. G. Koh, J. M. Kelly, Application of fractional derivatives to seismic analysis of base-isolated models, Earthquake engineering & structural dynamics 19 (2) (1990) 229–241.
  * [25] R. Magin, Fractional calculus in bioengineering, part 1, Critical Reviews in Biomedical Engineering 32 (1) (2004).
  * [26] R. Magin, M. Ovadia, Modeling the cardiac tissue electrode interface using fractional calculus, Journal of Vibration and Control 14 (9-10) (2008) 1431–1442.
  * [27] F. A. Rihan, D. A. Rahman, S. Lakshmanan, A. Alkhajeh, A time delay model of tumour–immune system interactions: Global dynamics, parameter estimation, sensitivity analysis, Applied Mathematics and Computation 232 (2014) 606–623.
  * [28] F. A. Rihan, D. Abdelrahman, F. Al-Maskari, F. Ibrahim, M. A. Abdeen, Delay differential model for tumour-immune response with chemoimmunotherapy and optimal control, Computational and mathematical methods in medicine 2014 (2014).
  * [29] C. Baker, G. Bocharov, C. Paul, F. Rihan, Modelling and analysis of time-lags in some basic patterns of cell proliferation, Journal of mathematical biology 37 (4) (1998) 341–371.
  * [30] R. Garrappa, E. Kaslik, On initial conditions for fractional delay differential equations, Communications in Nonlinear Science and Numerical Simulation 90 (2020) 105359.
  * [31] J. D. Meiss, Differential dynamical systems, SIAM, 2007.
  * [32] J. K. Hale, H. Koçak, Dynamics and bifurcations, Vol. 3, Springer Science & Business Media, 2012.
  * [33] S. Bhalekar, Stability and bifurcation analysis of a generalized scalar delay differential equation, Chaos: An Interdisciplinary Journal of Nonlinear Science 26 (8) (2016) 084306.
  * [34] S. B. Bhalekar, Stability analysis of a class of fractional delay differential equations, Pramana 81 (2) (2013) 215–224.
  * [35] K. Diethelm, The analysis of fractional differential equations, Springer Berlin, Heidelberg, 2010.
  * [36] J. A. Gallas, Structure of the parameter space of the hénon map, Physical Review Letters 70 (18) (1993) 2714.
  * [37] S. Kodba, M. Perc, M. Marhl, Detecting chaos from a time series, European Journal of Physics 26 (1) (2004) 205.


By [37, Theorem 5.1.2], $\mathcal{S}(v)\in\mathcal{X}$ (solutions actually
enjoy further Hölder regularity). We have the following standard
representation formula: let $\Phi$ be the fundamental solution to the
diffusion equation
$\partial_{t}u-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}u)=0$
on $[0,T]\times({\mathbb{R}}^{d})^{N}$; that is,

$\Phi(t,\boldsymbol{x})\vcentcolon=\frac{1}{\sqrt{\det(4\pi(\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})t)}\,}\,e^{-\tfrac{(\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})^{-1}\boldsymbol{x}\cdot\boldsymbol{x}}{4t}}=\frac{1}{\sqrt{1+Nd\beta/\sigma}\,}\frac{1}{(4\pi\sigma
t)^{\frac{Nd}{2}}}e^{-\tfrac{\left(\boldsymbol{\mathsf{I}}-\frac{\beta}{\sigma+Nd\beta}\boldsymbol{\mathsf{J}}\right)\boldsymbol{x}\cdot\boldsymbol{x}}{4\sigma
t}}.$

Then for $v\in\mathcal{X}$ and $i\in\\{1,\dots,N\\}$, we have

$\mathcal{S}(v)^{i}(t,\boldsymbol{x})\vcentcolon=(\Phi(T-t,\mspace{1.8mu}\cdot\mspace{1.8mu})\star
g^{i})(\boldsymbol{x})+\int_{t}^{T}(\Phi(T-s,\mspace{1.8mu}\cdot\mspace{1.8mu})\star(f^{i}-w^{i}(s-t,\mspace{1.8mu}\cdot\mspace{1.8mu}))(\boldsymbol{x})\,\mathrm{d}s,$

where $\star$ denotes the convolution operator over $({\mathbb{R}}^{d})^{N}$
and $w^{i}\vcentcolon=\frac{1}{2}|D_{i}v^{i}|^{2}+\sum_{j\neq
i}D_{j}v^{j}\cdot D_{j}v^{i}$. Then
$\lVert\mathcal{S}(v)^{i}\rVert_{\infty}\lesssim\lVert
g\rVert_{\mathcal{X}}+(T-\tau)(\lVert f\rVert_{\mathcal{X}}+\lVert
v\rVert_{\mathcal{X}}^{2})$, where the implied constant depends only on $N$
and $d$. Also, for $j\in\\{1,\dots,N\\}$ and $k\in\\{1,\dots,d\\}$

$\begin{multlined}D_{x^{jk}}\mathcal{S}(v)^{i}(t,\boldsymbol{x})=(\Phi(T-t,\mspace{1.8mu}\cdot\mspace{1.8mu})\star
D_{x^{jk}}g^{i})(\boldsymbol{x})+\int_{t}^{T}(\Phi(T-s,\mspace{1.8mu}\cdot\mspace{1.8mu})\star
D_{x^{jk}}f^{i})(\boldsymbol{x})\,\mathrm{d}s\vspace{-7pt}\\\
-\int_{t}^{T}(D_{x^{jk}}\Phi(T-s,\mspace{1.8mu}\cdot\mspace{1.8mu})\star
w^{i}(s-t,\mspace{1.8mu}\cdot\mspace{1.8mu}))(\boldsymbol{x})\,\mathrm{d}s;\end{multlined}D_{x^{jk}}\mathcal{S}(v)^{i}(t,\boldsymbol{x})=(\Phi(T-t,\mspace{1.8mu}\cdot\mspace{1.8mu})\star
D_{x^{jk}}g^{i})(\boldsymbol{x})+\int_{t}^{T}(\Phi(T-s,\mspace{1.8mu}\cdot\mspace{1.8mu})\star
D_{x^{jk}}f^{i})(\boldsymbol{x})\,\mathrm{d}s\vspace{-7pt}\\\
-\int_{t}^{T}(D_{x^{jk}}\Phi(T-s,\mspace{1.8mu}\cdot\mspace{1.8mu})\star
w^{i}(s-t,\mspace{1.8mu}\cdot\mspace{1.8mu}))(\boldsymbol{x})\,\mathrm{d}s;$

then, using that $\int_{({\mathbb{R}}^{d})^{N}}\lvert
D_{x^{jk}}\Phi(t,\boldsymbol{y})|\,\mathrm{d}\boldsymbol{y}\leqslant
Nd/{\sqrt{\pi\sigma t}\,}$, one has $\lVert
D\mathcal{S}(v)^{i}\rVert_{\infty}\lesssim\lVert
g\rVert_{\mathcal{X}}+(T-\tau)\lVert
f\rVert_{\mathcal{X}}+\sqrt{T-\tau}\,\lVert v\rVert^{2}_{\mathcal{X}}$, the
implied constant depending only on $\sigma$, $N$ and $d$. Furthermore, given
$v,\overline{v}\in\mathcal{X}$ similar computations show that
$\lVert\mathcal{S}(v)-\mathcal{S}(\overline{v})\rVert_{\mathcal{X}}\lesssim(1+\sqrt{T-\tau}\,)\sqrt{T-\tau}\max\\{\lVert
v\rVert_{\mathcal{X}},\lVert\overline{v}\rVert_{\mathcal{X}}\\}\lVert
v-\overline{v}\rVert_{\mathcal{X}}$, with implied constant depending only on
$\sigma$, $N$ and $d$. Letting
$\mathcal{B}_{c}\vcentcolon=\\{v\in\mathcal{X}:\ \lVert
v\rVert_{\mathcal{X}}\leqslant c\\}$, we deduce that
$\mathcal{S}\colon\mathcal{B}_{c}\to\mathcal{B}_{c}$ is a contraction for any
$c=c(\lVert g\rVert_{\mathcal{X}},\sigma,N,d)$ sufficiently large and
$\tau=\tau(c,\sigma,N,d,\beta,\lVert f\rVert_{\mathcal{X}})$ such that
$T-\tau$ is sufficiently small. It follows by the Banach–Caccioppoli theorem
that $\mathcal{S}$ has a unique fixed point $u\in\mathcal{B}_{c}$, which by
the definition of $\mathcal{S}$ solves system (1.1) in the mild sense. Then
note that by standard parabolic Schauder estimates we can bootstrap up to the
desired regularity.

Step 2. In order to see that such a solution extends to the whole time horizon
$[0,T]$, we consider the following a priori estimates, that hold for classical
solutions that are $C^{3}$ in space and $C^{1}$ in time (therefore they will
apply to the short-time solutions obtained in the previous step, since
$\gamma\geqslant 1$). Recall the representation formula for the solution of
the Nash system: for any
$(s,\boldsymbol{z})\in[\tau,T]\times({\mathbb{R}}^{d})^{N}$,

(A.1)
$u^{i}(s,\boldsymbol{z})=\int_{s}^{T}\int_{({\mathbb{R}}^{d})^{N}}\Bigl{(}\,\frac{1}{2}|D_{i}u^{i}(t,\boldsymbol{x})|^{2}+f^{i}(\boldsymbol{x})\Bigr{)}\rho(t,\mathrm{d}\boldsymbol{x})\mathrm{d}t+\int_{({\mathbb{R}}^{d})^{N}}g^{i}(\boldsymbol{x})\,\rho(T,\mathrm{d}\boldsymbol{x}),$

where $\rho=\rho^{s,\boldsymbol{z}}$ solves
$\partial_{t}\rho-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}\rho)+\sum_{1\leqslant
j\leqslant N}\mathrm{div}_{x^{j}}(D_{j}u^{j}\rho)=0$ on
$[s,T]\times({\mathbb{R}}^{d})^{N}$ with
$\rho|_{t=s}=\delta_{\boldsymbol{z}}$. This is easily obtained by testing the
equation of $\rho$ by $u^{i}$. In particular,

$u^{i}(s,\boldsymbol{z})\geqslant\int_{s}^{T}\int_{({\mathbb{R}}^{d})^{N}}f^{i}(\boldsymbol{y})\,\rho(t,\mathrm{d}\boldsymbol{x})\mathrm{d}t+\int_{({\mathbb{R}}^{d})^{N}}g^{i}(\boldsymbol{x})\,\rho(T,\mathrm{d}\boldsymbol{x}).$

On the other hand, since
$-\partial_{t}u^{i}-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}u^{i})+\sum_{1\leqslant
j\leqslant N,\,j\neq i}D_{j}u^{j}\cdot D_{j}u^{i}\leqslant f^{i}$, if
$\hat{\rho}^{i}$ solves instead
$\partial_{t}\hat{\rho}^{i}-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}\hat{\rho}^{i})+\sum_{1\leqslant
j\leqslant N,\,j\neq i}\mathrm{div}_{x^{j}}(D_{j}u^{j}\hat{\rho}^{i})=0$ then

$u^{i}(s,\boldsymbol{z})\leqslant\int_{s}^{T}\int_{({\mathbb{R}}^{d})^{N}}f^{i}(\boldsymbol{x})\,\hat{\rho}^{i}(t,\mathrm{d}\boldsymbol{x})\mathrm{d}t+\int_{({\mathbb{R}}^{d})^{N}}g^{i}(\boldsymbol{x})\,\hat{\rho}^{i}(T,\mathrm{d}\boldsymbol{x}).$

Then $\sup_{t\in[\tau,T]}\lVert
u^{i}(t,\cdot)\rVert_{\infty}\leqslant(T+1)\max\\{\lVert
f^{i}\rVert_{\infty},\lVert g^{i}\rVert_{\infty}\\}$.

We now want to control the gradient of $u^{i}$. As a preliminary step, we show
that the quantity

$\mathfrak{c}=\mathfrak{c}(s,\boldsymbol{z})\vcentcolon=\int_{s}^{T}\int_{({\mathbb{R}}^{d})^{N}}\sum_{1\leqslant
i,j\leqslant
N}|D_{j}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}\,\mathrm{d}\rho(t)\mathrm{d}t,$

can be bounded uniformly in $(s,\boldsymbol{z})$. Note indeed that
$v^{i}\vcentcolon=\frac{1}{2}|u^{i}|^{2}$ satisfies

$\begin{multlined}-\partial_{t}v^{i}-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}v^{i})+\sigma\\!\sum_{1\leqslant
j\leqslant N}\\!|D_{j}u^{i}|^{2}+\beta\biggl{|}\sum_{1\leqslant j\leqslant
N}D_{j}u^{i}\biggr{|}^{2}\\\
+\frac{1}{2}D_{i}u^{i}D_{i}v^{i}+\\!\sum_{\begin{subarray}{c}1\leqslant
i,j\leqslant N\\\ j\neq
i\end{subarray}}\\!D_{j}u^{j}D_{j}v^{i}=f^{i}u^{i},\end{multlined}-\partial_{t}v^{i}-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}v^{i})+\sigma\\!\sum_{1\leqslant
j\leqslant N}\\!|D_{j}u^{i}|^{2}+\beta\biggl{|}\sum_{1\leqslant j\leqslant
N}D_{j}u^{i}\biggr{|}^{2}\\\
+\frac{1}{2}D_{i}u^{i}D_{i}v^{i}+\\!\sum_{\begin{subarray}{c}1\leqslant
i,j\leqslant N\\\ j\neq i\end{subarray}}\\!D_{j}u^{j}D_{j}v^{i}=f^{i}u^{i},$

so that by testing the equation of $\rho$ by $v^{i}$ one gets

$\begin{multlined}\sigma\int_{\tau}^{T}\int_{({\mathbb{R}}^{d})^{N}}\sum_{1\leqslant
j\leqslant
N}|D_{j}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}\,\mathrm{d}\rho(t)\mathrm{d}t\vspace{-5pt}\\\
\leqslant\int_{({\mathbb{R}}^{d})^{N}}\frac{1}{2}|g^{i}|^{2}\,\mathrm{d}\rho(T)+\int_{\tau}^{T}\int_{({\mathbb{R}}^{d})^{N}}\frac{1}{2}|D_{i}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})\,\mathrm{d}\rho(t)\mathrm{d}t+\int_{\tau}^{T}\int_{({\mathbb{R}}^{d})^{N}}f^{i}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})\,\mathrm{d}\rho(t)\mathrm{d}t.\end{multlined}\sigma\int_{\tau}^{T}\int_{({\mathbb{R}}^{d})^{N}}\sum_{1\leqslant
j\leqslant
N}|D_{j}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}\,\mathrm{d}\rho(t)\mathrm{d}t\vspace{-5pt}\\\
\leqslant\int_{({\mathbb{R}}^{d})^{N}}\frac{1}{2}|g^{i}|^{2}\,\mathrm{d}\rho(T)+\int_{\tau}^{T}\int_{({\mathbb{R}}^{d})^{N}}\frac{1}{2}|D_{i}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})\,\mathrm{d}\rho(t)\mathrm{d}t+\int_{\tau}^{T}\int_{({\mathbb{R}}^{d})^{N}}f^{i}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})\,\mathrm{d}\rho(t)\mathrm{d}t.$

Then, letting $A\vcentcolon=\max\\{\lVert f^{i}\rVert_{\infty},\lVert
g^{i}\rVert_{\infty}\\}$ and noting that (A.1) yields

$\int_{\tau}^{T}\int_{({\mathbb{R}}^{d})^{N}}\frac{1}{2}|D_{i}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}\,\mathrm{d}\rho(t)\mathrm{d}t\leqslant
2(T+1)A,$

we obtain

(A.2)
$\mathfrak{c}\leqslant\overline{\mathfrak{c}}:=\Bigl{(}\,\frac{1}{2}+T(T+1)(2T+3)\Bigr{)}\frac{NA^{2}}{\sigma}.$

Consider now that $\tilde{w}\vcentcolon=\sum_{1\leqslant i,j\leqslant
N}|D_{j}u^{i}|^{2}$ solves

$\begin{multlined}-\partial_{t}\tilde{w}-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}\tilde{w})+\sigma\\!\sum_{1\leqslant
i,j,k\leqslant N}|D^{2}_{jk}u^{i}|^{2}+\beta\biggl{|}\sum_{1\leqslant
i,j,k\leqslant N}D^{2}_{jk}u^{i}\biggr{|}^{2}+\sum_{1\leqslant k\leqslant
N}D_{k}u^{k}D_{k}\tilde{w}\\\ =\sum_{1\leqslant i,j\leqslant
N}D_{j}f^{i}D_{j}u^{i}-\sum_{\begin{subarray}{c}1\leqslant i,j,k\leqslant N\\\
k\neq
i\end{subarray}}D^{2}_{kj}u^{k}D_{k}u^{i}D_{j}u^{i}.\end{multlined}-\partial_{t}\tilde{w}-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}\tilde{w})+\sigma\\!\sum_{1\leqslant
i,j,k\leqslant N}|D^{2}_{jk}u^{i}|^{2}+\beta\biggl{|}\sum_{1\leqslant
i,j,k\leqslant N}D^{2}_{jk}u^{i}\biggr{|}^{2}+\sum_{1\leqslant k\leqslant
N}D_{k}u^{k}D_{k}\tilde{w}\\\ =\sum_{1\leqslant i,j\leqslant
N}D_{j}f^{i}D_{j}u^{i}-\sum_{\begin{subarray}{c}1\leqslant i,j,k\leqslant N\\\
k\neq i\end{subarray}}D^{2}_{kj}u^{k}D_{k}u^{i}D_{j}u^{i}.$

By testing the equation of $\rho$ by $w$, using Young and Hölder inequality,
we deduce that

$\begin{multlined}\lVert\tilde{w}(s,\cdot)\rVert_{\infty}\leqslant\frac{1}{2}\sum_{1\leqslant
i,j\leqslant
N}\int_{({\mathbb{R}}^{d})^{N}}|D_{j}g^{i}|^{2}\,\mathrm{d}\rho(T)+\frac{1}{2}\sum_{1\leqslant
i,j\leqslant
N}\int_{s}^{T}\int_{({\mathbb{R}}^{d})^{N}}|D_{j}f^{i}|^{2}\,\mathrm{d}\rho(t)\mathrm{d}t\\\
+\int_{s}^{T}\Biggl{(}1+\frac{1}{\sigma}\\!\int_{({\mathbb{R}}^{d})^{N}}\sum_{\begin{subarray}{c}1\leqslant
i,j\leqslant N\\\ j\neq
i\end{subarray}}|D_{j}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}\,\mathrm{d}\rho(t)\Biggr{)}\lVert\tilde{w}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})\rVert_{\infty}\,\mathrm{d}t.\end{multlined}\lVert\tilde{w}(s,\cdot)\rVert_{\infty}\leqslant\frac{1}{2}\sum_{1\leqslant
i,j\leqslant
N}\int_{({\mathbb{R}}^{d})^{N}}|D_{j}g^{i}|^{2}\,\mathrm{d}\rho(T)+\frac{1}{2}\sum_{1\leqslant
i,j\leqslant
N}\int_{s}^{T}\int_{({\mathbb{R}}^{d})^{N}}|D_{j}f^{i}|^{2}\,\mathrm{d}\rho(t)\mathrm{d}t\\\
+\int_{s}^{T}\Biggl{(}1+\frac{1}{\sigma}\\!\int_{({\mathbb{R}}^{d})^{N}}\sum_{\begin{subarray}{c}1\leqslant
i,j\leqslant N\\\ j\neq
i\end{subarray}}|D_{j}u^{i}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})|^{2}\,\mathrm{d}\rho(t)\Biggr{)}\lVert\tilde{w}(t,\mspace{1.8mu}\cdot\mspace{1.8mu})\rVert_{\infty}\,\mathrm{d}t.$

(For the sake of accuracy, we point out that one should consider
$\rho^{s,\boldsymbol{z}}$ such that
$\tilde{w}(s,\boldsymbol{z})\geqslant\lVert\tilde{w}(s,\cdot)\rVert_{\infty}-\varepsilon$
for $\varepsilon>0$, and let $\varepsilon\to 0$ in the following
computations.) Whence, for

(A.3) $s\geqslant
s^{*}\vcentcolon=T-(2(T+\sigma^{-1}\overline{\mathfrak{c}}))^{-1}$

one has $\sup_{t\in[s^{*}\vee
0,T]}\lVert\tilde{w}(t,\cdot)\rVert_{\infty}\leqslant C$, where $C$ depends
only on $\sigma$, $N$, $\max_{ij}\\{\lVert D_{j}f^{i}\rVert_{\infty},\lVert
D_{j}g^{i}\rVert_{\infty}\\}$ and $T$. Then by iterating the estimate a finite
number of times $\kappa$, one obtains
$\sup_{t\in[\tau,T]}\lVert\tilde{w}(t,\cdot)\rVert_{\infty}\leqslant
2^{\kappa}C$.

This is sufficient to conclude that $\lVert
u\rVert_{\mathcal{X}_{\tau}}\leqslant C^{\prime}$ for some constant
$C^{\prime}$ independent of $\tau$. Define at this point
$\tau_{*}\vcentcolon=\inf\\{\tau\in[0,T]:\ \text{\eqref{ENS_NS} has a unique
solution belonging to $\mathcal{X}_{\tau}$}\\}$. We proved that $\tau_{*}<T$;
seeking for a contradiction, suppose that $\tau_{*}>0$. Let
$\varepsilon\in(0,T-\tau_{*})$ to be announced and consider a solution
$u\in\mathcal{X}_{\tau_{*}+\varepsilon}$. We have proved that $\lVert
u\rVert_{\mathcal{X}_{\tau_{*}+\varepsilon}}\leqslant C^{\prime}$. Then one
can redo the initial fixed point argument with $u(\tau_{*}+\varepsilon,\cdot)$
in lieu of $g$ and conclude that $u$ extends to the horizon
$[\tau^{\prime},T]$, provided that
$\tau_{*}+\varepsilon-\tau^{\prime}\leqslant\delta$ for some $\delta>0$ which
is independent of $\tau_{*}$ and $\varepsilon$. So one can fix
$\delta_{0}=\frac{\tau_{*}}{2}\wedge\delta$,
$\varepsilon=\frac{1}{2}({\delta_{0}}\wedge(T-\tau_{*}))$ and
$\tau^{\prime}=\tau_{*}-\frac{\delta_{0}}{2}\in(0,\tau_{*})$, thus
contradicting the definition of $\tau_{*}$ ∎

##  Appendix B An estimate on the space-time regularity for the heat equation

The Lemma below shows how to transfer weighted Lipschitz regularity in space
of solution to heat equations to Hölder regularity in time, following a
standard idea proposed in [30]. We state it for forward equations for
simplicity. Note again that Hölder continuity is independent of $N$.

######  Lemma B.1.

Let $f:[0,T]\times({\mathbb{R}}^{d})^{N}\to{\mathbb{R}}$, $\sigma>0$,
$\beta\geqslant 0$, and $u$ be a classical solution of

$\partial_{t}u-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}u)=f\qquad\text{in
$(0,T)\times({\mathbb{R}}^{d})^{N}$ }.$

Assume also that

$|u(t,\boldsymbol{x})-u(t,\boldsymbol{y})|\leqslant
c_{1}\lVert\boldsymbol{x}-\boldsymbol{y}\rVert^{i}\qquad\forall\,t\in[0,T],\
\boldsymbol{x},\boldsymbol{y}\in({\mathbb{R}}^{d})^{N},$

and $\|f\|_{\infty},\|u\|_{\infty}\leqslant c_{2}$ for some $c_{1},c_{2}>0$
and $i=1,\ldots,N$. Then,

$|u(\tau,\boldsymbol{z})-u(s,\boldsymbol{z})|\leqslant(c_{1}+4(\sigma+\beta)dc_{2})|\tau-s|^{\frac{1}{3}}+c_{2}|\tau-s|\qquad\forall\,\tau,s\in[0,T],\
\boldsymbol{z}\in({\mathbb{R}}^{d})^{N}.$

###### Proof.

Assume first that $s=0$, $\boldsymbol{z}=0$ and $u(0,0)=0$. Let $\rho>0$ and

$v(t,\boldsymbol{x})=c_{1}\rho+\frac{4(\sigma+\beta)dc_{2}}{\rho^{2}}\,t+c_{2}t+\frac{c_{2}}{\rho^{2}}(\|x\|^{i})^{2}.$

Since

$\partial_{t}v-\operatorname{tr}((\sigma\boldsymbol{\mathsf{I}}+\beta\boldsymbol{\mathsf{J}})D^{2}v)-f=\frac{4(\sigma+\beta)dc_{2}}{\rho^{2}}+c_{2}-\frac{(\sigma+\beta)dc_{2}}{\rho^{2}}\Bigl{(}4-\frac{2}{N}\Bigr{)}-f\geqslant
0$

everywhere, and

$\begin{split}&v(0,\boldsymbol{x})\geqslant c_{1}\rho\geqslant
c_{1}\|\boldsymbol{x}\|^{i}\geqslant u(0,\boldsymbol{x})\qquad\text{on
$\\{\|\boldsymbol{x}\|^{i}\leqslant\rho\\}$},\\\
&v(t,\boldsymbol{x})\geqslant\frac{c_{2}}{\rho^{2}}(\|x\|^{i})^{2}\geqslant
u(t,\boldsymbol{x})\qquad\text{on
$[0,T]\times\\{\|\boldsymbol{x}\|^{i}=\rho\\}$},\end{split}$

by the comparison principle we get $u(t,\boldsymbol{x})\leqslant
v(t,\boldsymbol{x})$ on
$[0,T]\times\\{\|\boldsymbol{x}\|^{i}\leqslant\rho\\}$. For $\tau>0$, pick
then $\rho=\tau^{\frac{1}{3}}$, so that the previous inequality yields

$u(\tau,0)\leqslant
c_{1}\tau^{\frac{1}{3}}+4(\sigma+\beta)dc_{2}\tau^{\frac{1}{3}}+c_{2}\tau.$

Using $-v$ as a subsolution yields the analogous bound from below. It is now
straightforward to conclude the estimate for general $s,\boldsymbol{z}$. ∎

## References

  * [1] S. Ahuja.  Wellposedness of Mean Field Games with common noise under a weak monotonicity condition.  SIAM Journal on Control and Optimization, 54(1):30–48, Jan 2016\.
  * [2] B. Andrews and J. Clutterbuck.  Time-interior gradient estimates for quasilinear parabolic equations.  Indiana Univ. Math. J., 58(1):351–380, 2009.
  * [3] A. Aurell, R. Carmona, and M. Laurière.  Stochastic graphon games: II. The linear-quadratic case.  Appl. Math. Optim., 85(3):Paper No. 26, 33, 2022.
  * [4] P. Baldi.  Stochastic calculus.  Universitext. Springer, Cham, 2017.  An introduction through theory and exercises.
  * [5] M. Bansil and A. R. Mészáros.  Hidden monotonicity and canonical transformations for mean field games and master equations, 2024.  arXiv:2403.05426.
  * [6] A. Bensoussan and J. Frehse.  $C^{\alpha}$-regularity results for quasilinear parabolic systems.  Comment. Math. Univ. Carolin., 31(3):453–474, 1990.
  * [7] A. Bensoussan and J. Frehse.  Smooth solutions of systems of quasilinear parabolic equations.  volume 8, pages 169–193. 2002.  A tribute to J. L. Lions.
  * [8] A. Bensoussan, T. Huang, and M. Laurière.  Mean field control and mean field game models with several populations.  Minimax Theory Appl., 3(2):173–209, 2018.
  * [9] C. Bertucci, J.-M. Lasry, and P.-L. Lions.  On Lipschitz solutions of mean field games master equations, 2023.  arXiv:2302.05218.
  * [10] V. I. Bogachev, N. V. Krylov, M. Röckner, and S. V. Shaposhnikov.  Fokker-Planck-Kolmogorov equations, volume 207 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 2015.
  * [11] P. E. Caines and M. Huang.  Graphon mean field games and their equations.  SIAM J. Control Optim., 59(6):4373–4399, 2021.
  * [12] P. Cardaliaguet.  Notes on mean field games.  Lecture notes, 2012.  From P.-L. Lions’ lectures at Collège de France.
  * [13] P. Cardaliaguet, M. Cirant, and A. Porretta.  Splitting methods and short time existence for the master equations in mean field games.  J. Eur. Math. Soc. (JEMS), 25(5):1823–1918, 2023.
  * [14] P. Cardaliaguet, F. Delarue, J.-M. Lasry, and P.-L. Lions.  The master equation and the convergence problem in mean field games, volume 201 of Annals of Mathematics Studies.  Princeton University Press, Princeton, NJ, 2019.
  * [15] P. Cardaliaguet and A. Porretta.  An introduction to mean field game theory.  In Mean field games, volume 2281 of Lecture Notes in Math., pages 1–158. Springer, Cham, [2020] ©2020.
  * [16] R. Carmona and F. Delarue.  Probabilistic theory of mean field games with applications. I, volume 83 of Probability Theory and Stochastic Modelling.  Springer, Cham, 2018.  Mean field FBSDEs, control, and games.
  * [17] R. Carmona and F. Delarue.  Probabilistic theory of mean field games with applications. II, volume 84 of Probability Theory and Stochastic Modelling.  Springer, Cham, 2018.  Mean field games with common noise and master equations.
  * [18] M. Cirant.  Multi-population mean field games systems with Neumann boundary conditions.  J. Math. Pures Appl. (9), 103(5):1294–1315, 2015.
  * [19] M. Cirant and A. Porretta.  Long time behavior and turnpike solutions in mildly non-monotone mean field games.  ESAIM Control Optim. Calc. Var., 27:Paper No. 86, 40, 2021.
  * [20] M. Cirant and D. F. Redaelli.  Some remarks on linear-quadratic closed-loop games with many players, 2024\.  arXiv:2401.06534, to appear in Dyn. Games Appl.
  * [21] M. F. Djete.  Large population games with interactions through controls and common noise: convergence results and equivalence between open-loop and closed-loop controls.  ESAIM Control Optim. Calc. Var., 29:Paper No. 39, 42, 2023.
  * [22] J. Frehse.  Remarks on diagonal elliptic systems.  In Partial differential equations and calculus of variations, volume 1357 of Lecture Notes in Math., pages 198–210. Springer, Berlin, 1988.
  * [23] A. Friedman.  Stochastic differential games.  J. Differential Equations, 11:79–108, 1972.
  * [24] W. Gangbo, A. R. Mészáros, C. Mou, and J. Zhang.  Mean field games master equations with nonseparable Hamiltonians and displacement monotonicity.  Ann. Probab., 50(6):2178–2217, 2022.
  * [25] P. J. Graber and A. R. Mészáros.  On monotonicity conditions for mean field games.  J. Funct. Anal., 285(9):Paper No. 110095, 45, 2023.
  * [26] Z. Hedrlín.  On integration in compact metric spaces.  Commentationes Mathematicae Universitatis Carolinae, 002(4):17–19, 1961.
  * [27] M. Huang, R. P. Malhamé, and P. E. Caines.  Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle.  Commun. Inf. Syst., 6(3):221–251, 2006.
  * [28] J. Jackson.  On quasilinear parabolic systems and FBSDEs of quadratic growth.  Ann. Appl. Probab., 34(1A):357–387, 2024.
  * [29] J. Jackson and L. Tangpi.  Quantitative convergence for displacement monotone mean field games with controlled volatility, 2023.  arXiv:2304.04543.
  * [30] S. N. Kružkov.  Results on the nature of the continuity of solutions of parabolic equations, and certain applications thereof.  Mat. Zametki, 6:97–108, 1969.
  * [31] D. Lacker.  On the convergence of closed-loop Nash equilibria to the mean field game limit.  Ann. Appl. Probab., 30(4):1693–1761, 2020.
  * [32] D. Lacker and L. Le Flem.  Closed-loop convergence for mean field games with common noise.  Ann. Appl. Probab., 33(4):2681–2733, 2023.
  * [33] D. Lacker and A. Soret.  A label-state formulation of stochastic graphon games and approximate equilibria on large networks.  Math. Oper. Res., 48(4):1987–2018, 2023.
  * [34] O. A. Ladyzenskaja, V. A. Solonnikov, and N. N. Ural’ceva.  Linear and quasilinear equations of parabolic type.  Translated from the Russian by S. Smith. Translations of Mathematical Monographs, Vol. 23. American Mathematical Society, Providence, R.I., 1968.
  * [35] J.-M. Lasry and P.-L. Lions.  Mean field games.  Jpn. J. Math., 2(1):229–260, 2007.
  * [36] P.-L. Lions.  Cours du Collège de France, available at https://www.college-de-france.fr/fr/agenda/cours/jeux-champ-moyen-suite-2, 2008\.
  * [37] A. Lunardi.  Analytic semigroups and optimal regularity in parabolic problems, volume 16 of Progress in Nonlinear Differential Equations and their Applications.  Birkhäuser Verlag, Basel, 1995.
  * [38] A. R. Mészáros and C. Mou.  Mean field games systems under displacement monotonicity.  SIAM J. Math. Anal., 56(1):529–553, 2024.
  * [39] C. Mou and J. Zhang.  Mean field game master equations with anti-monotonicity conditions, 2022\.  arXiv:2201.10762.
  * [40] G. Peyré and M. Cuturi.  Computational optimal transport.  Foundations and Trends in Machine Learning, 11(5-6):355–607, 2019\.
  * [41] A. Porretta and E. Priola.  Global Lipschitz regularizing effects for linear and nonlinear parabolic equations.  J. Math. Pures Appl. (9), 100(5):633–686, 2013.
  * [42] D. Possamaï and L. Tangpi.  Non-asymptotic convergence rates for mean-field games: weak formulation and McKean-Vlasov BSDEs, 2021.  arXiv:2105.00484.
  * [43] C. Villani.  Optimal transport, volume 338 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, 2009.  Old and new.


# State estimation of a moving frequency source from observations at multiple
receivers

Michela Mancini
Guggenheim School of Aerospace Engineering
Georgia Institute of Technology
Atlanta, GA, 30332
<EMAIL_ADDRESS>
&Anton Leykin
School of Mathematics
Georgia Institute of Technology
Atlanta, GA, 30332
&John A. Christian
Guggenheim School of Aerospace Engineering
Georgia Institute of Technology
Atlanta, GA, 30332

###### Abstract

The task of position and velocity estimation of a moving transmitter (with
either a known or unknown frequency) is a common problem arising in many
different application domains. Based on the Doppler effect, this work presents
a direct solution using only the frequency measured by a multitude of
receivers with a known state. A natural rewriting of the problem as a system
of polynomial equations allows for the use of homotopy continuation to find
the global solution without any _a priori_ information about the frequency
source. We show that the data from six or seven receivers is sufficient in
case of known or unknown frequency, respectively. After a brief development of
the mathematics, two simple examples are provided: (1) position and velocity
estimation of a vocalizing dolphin emitting an acoustic signal and (2) initial
orbit determination of a satellite emitting an electromagnetic signal.

##  1 INTRODUCTION

Despite the criticism that C. A. Doppler encountered after publishing his
famous principle in 1842 ([1, 2]), the Doppler effect has subsequently proved
to be an important phenomenon in many different application domains. Indeed,
we can find examples of its use in many different disciplines, including
medicine ([3]), astronomy ([4]), satellite orbit determination ([5]), and many
others. In this work, we consider how the Doppler effect can be exploited to
simultaneously determine the position and velocity of an object emitting a
signal—acoustic or electromagnetic—with either a known or unknown transmit
frequency.

Consider a particle with unknown state (position and velocity) that emits a
signal. Suppose this signal is subsequently detected by a receiver with a
known state (position and velocity). If there is relative motion between the
two objects, then the Doppler effect tells us that the frequency measured by
the receiver may be different from the frequency at the transmitter. A
visualization of this phenomenon is given in Fig. 1. We will show that the
classical equation for the Doppler effect may be rearranged to yield a
polynomial in the unknown transmitter state. If the frequency at the source is
known, six receivers suffice to determine the position and velocity of the
transmitter up to finitely many possibilities. If it is unknown, seven
receivers are needed, and the unknown frequency can be estimated as well. The
polynomial system obtained can be solved directly via homotopy continuation,
without the need of any initial guess.

Figure 1: Illustration of the Doppler effect. If the transmitter (black dot)
and a receiver have a positive range rate (i.e., are moving away from one
another) then the measured frequency will decrease (e.g., blue dot).
Conversely, if the transmitter and receiver have a negative range rate (i.e.
are approaching one another) then the measured frequecy will increase (e.g.,
green or red dots).

The literature presents many circumstances where localizing an unknown
transmitter is encountered. Among the marine applications, for example, it can
be used to determine the position and velocity of a vocalizing underwater
animal ([6]). It may also be used to localize nodes within an underwater
sensor network, which are widely used to collect data for underwater
exploration, environmental monitoring, natural disaster prevention ([7]).
Doppler-based state estimation may also be used to track or navigate moving
vehicles on the ground, in air, or in space.

We have examples of many previous solutions in all these fields, including for
underwater device detection ([8, 9]), drone localization and tracking ([10,
11]), and initial orbit determination (IOD) ([12, 13, 14]). However, to the
best of our knowledge, most of the past approaches process Doppler
measurements within the framework of a sequential filter (e.g., some variant
of a Kalman filter) that requires an initial guess of the transmitter
state—with it often being unclear how such an _a priori_ state estimate might
be obtained. An attempt to solve the initialization problem was recently
studied for the problem of spacecraft IOD ([14]), but the runtimes for the
general problem were prohibitive.

In this work we introduce a fast and simple solution to the generic Doppler-
only state estimation problem based on homotopy continuation. This allows the
analyst to find the global solution in the absence of _a priori_ state
information, thus providing a means of obtaining the initial guess required by
earlier work. We provide Macaulay2 source code and illustrate this method on
simple underwater and space examples.

##  2 Polynomial formulation

Consider a transmitter at the unknown location
$\textbf{\emph{r}}\in\mathbb{R}^{3}$ moving with the unknown velocity
$\textbf{\emph{v}}\in\mathbb{R}^{3}$ and emitting a signal with a frequency
$f\in\mathbb{R}_{>0}$. Assume that the signal is detected by a receiver,
placed at the position $\textbf{\emph{r}}_{i}$, moving with velocity
$\textbf{\emph{v}}_{i}$. This situation is shown in Fig. 1. Due to Doppler
effect, the frequency seen by the receiver $f_{i}$ will be related to the
frequency at the source through the relation

$f_{i}=\left(1-\frac{\dot{\rho}_{i}}{c}\right)f$ (1)

where $c$ is the signal’s propagation speed in the given medium, and the range
rate $\dot{\rho}_{i}$ can be obtained taking the derivative of the range
$\rho_{i}$ between the transmitter and the receiver

$\rho_{i}=\sqrt{\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)^{T}\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)}$
(2)

that yields

$\dot{\rho}_{i}=\frac{1}{\rho_{i}}\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)^{T}\left(\textbf{\emph{v}}_{i}-\textbf{\emph{v}}\right)$
(3)

A few simple manipulations lead to the polynomial expression we seek.
Rearrangement of Eqs. (1)–(3) allows us to write the product
$\rho_{i}\dot{\rho}_{i}$ as

$\rho_{i}\dot{\rho}_{i}=\sqrt{\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)^{T}\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)}\left(\frac{c\left(f-f_{i}\right)}{f}\right)=\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)^{T}\left(\textbf{\emph{v}}_{i}-\textbf{\emph{v}}\right)$
(4)

Squaring the middle and right-hand term of Eq. (4), we obtain a polynomial
expression relating the known parameters
$(\textbf{\emph{r}}_{i},\textbf{\emph{v}}_{i},f_{i})$ to the unknown
parameters $(\textbf{\emph{r}},\textbf{\emph{v}},f)$:

$c^{2}\left(f-f_{i}\right)^{2}\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)^{T}\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)-f^{2}\left[\left(\textbf{\emph{r}}_{i}-\textbf{\emph{r}}\right)^{T}\left(\textbf{\emph{v}}_{i}-\textbf{\emph{v}}\right)\right]^{2}=0$
(5)

The repetition of Eq. (5) for seven receivers produces a polynomial system of
seven equations in seven unknowns. If the transmit frequency $f$ is known, six
receivers may be used to solve for the remaining six unknowns.

Finally, it is important to note the $c^{2}$ coefficient in the left-most term
of Eq. (5). In applications with high propagation speeds (e.g., speed of
light), scaling can have a significant impact on the ability to solve such
polynomial systems on a finite precision computer. In such cases, rescaling of
the polynomials may provide considerable performance improvements.

##  3 Numerical Solution with Homotopy Continuation

The polynomial system formed by a multitude of frequency observations—each
represented by its own instantiation of Eq. (5)—can be solved in many
different ways. One family of techniques is homotopy continuation ([15, 16]).
In this work, we use an implementation of parameter homotopy provided by the
software package NAG4M2 ([17]) in the computer algebra system Macaulay2
([18]).

Consider a square system, $N$ equations in $N$ unknowns, of polynomials that
we aim to solve via homotopy continuation. The procedure begins by selecting a
_start system_ $\mathcal{S}$ which belongs to the same family as the _target
system_ $\mathcal{F}$ that we desire to solve—which, in this particular case,
is a system formed by multiple copies of Eq. (5). We intelligently choose the
start system to be one that is easily solvable and then proceed to compute the
corresponding solution (e.g., [19]). Then, beginning from the known solution
to the start system $\mathcal{S}$, homotopy continuation follows a path (in
the space of the systems) from start to target while deforming the
corresponding solutions—typically using an adaptive predictor-corrector
method. Specifically, we track the solutions of the system

$\mathcal{H}=(1-t)\,\mathcal{S}+t\,\mathcal{F}$ (6)

as the continuation parameter $t$ is varied from 0 to 1.

Construction of the family of polynomial systems is a key factor in this
approach. This is critical because the number of solutions of a generic system
within the polynomial family corresponds to the number of continuation paths
that must be tracked. For instance, many blackbox solvers implement the idea
of _total degree homotopy_. Total degree homotopy, an approach which works for
any square system, is attractive because it results in a simple start system
whose solutions are easy to write down. However, the number of continuation
paths equals the total degree in Bézout’s theorem (see chapter 18 of [20]),
which is often much larger than the number of target solutions we seek. The
consequence is that total degree homotopy often results in more paths than are
necessary. Moreover, tracking these excess paths (which are divergent)
significantly increases the computational cost. Reliance on total degree
homotopy likely led to exceptionally long runtimes in a similar approach
attempted recently ([14]).

A more nuanced approach is to use the natural parametric family arising in our
specific problem. Once the solutions to a generic system are known we use
_parameter homotopy_ to deform them into the target solutions. We can choose a
random start system in the family—this way the start and target systems are
structurally the same differing only in the values of the parameters.
Depending on the application, the speedup can be considerable. We created and
solved the start system using an approach based on monodromy ([19]). Once the
solutions to the start system are known, they can be used solve any other
instance of the problem in the family. A short example code that shows the
setup of the polynomial system and the calls necessary to perform the tracking
in Macaulay2 is provided in the Appendix.

We now apply these ideas to two different scenarios. The first is the generic
case, where both the source and receivers are permitted to move. The second is
a simplified special case where the source is moving, but the receivers are
stationary.

###  3.1 General case: moving receivers

Consider the polynomial system given by Eq. (5) specialized for seven
different receivers with an unknown transmit frequency $f$. According to
Bézout’s theorem, we can expect up to $6^{7}=279,936$ solutions in the form of
complex-valued 7-tuples
$\left(\textbf{\emph{r}},\,\textbf{\emph{v}},\,f\right)$. However, the actual
number of solutions that we found with techniques based on Gröbner basis
([21]) is only $672$. Using the total degree homotopy, almost all the paths
will diverge and many hours (on contemporary laptop computers) are required to
converge to the 672 solutions. Instead, parameter homotopy allows us to track
only $672$ solutions and solve the desired polynomial system more quickly. We
can repeat a similar reasoning for the case of known frequency $f$, where the
total degree homotopy requires the tracking of $4^{6}=4,096$ solutions as
6-tuples $\left(\textbf{\emph{r}},\,\textbf{\emph{v}}\right)$ , while the
parameter homotopy allows us to track the actual number of 128 solutions.

We have experimentally demonstrated that most of the 672 (for unknown $f$) or
128 (for known $f$) complex solutions can be eliminated without using
observations from any additional receivers. In fact, even if they all satisfy
the polynomial system obtained with Eq. (5), not all satisfy the original
(i.e., non-squared) relation from Eq. (4). However, it is not guaranteed that
this will eliminate all but one candidates. If an additional observation is
available, we can discard the remaining “wrong” solutions by checking the
residuals of Eq. (4) written for the additional receiver. The screening
procedure described here was sufficient to uniquely identify the desired
solution in all of the experiments we performed.

###  3.2 Special case: stationary receivers

Although the general case encompasses all classes of observer motion (both
moving and stationary), the special case of stationary observers
($\textbf{\emph{v}}_{i}=\textbf{0}$) simplifies the polynomial systems and
permits a more efficient solution. In particular, setting
$\textbf{\emph{v}}_{i}=\textbf{0}$ causes some of the terms in the original
polynomial vanish. Furthermore, we observe that if
$(\textbf{\emph{r}},\,\textbf{\emph{v}},\,f)$ is a solution to the polynomial
system, then the same must be true for
$(\textbf{\emph{r}},\,-\textbf{\emph{v}},\,f)$. Counting the number of
solutions for this special case, we found that the number of solutions of the
polynomial system with stationary observers is 296 for an unknown transmit
frequency $f$. However, thanks to the symmetry of the problem, we can discard
half of them and track only 148 paths. Even better, if the frequency at the
source is known, then we only need to track 24 paths. This reduction in the
number of paths increases the speed of the solver.

##  4 Experiments

The efficacy of the proposed solution is demonstrated through two illustrative
examples.

###  4.1 Example 1: Underwater localization with stationary receivers

The first illustrative experiment is the simultaneous position and velocity
estimation of a vocalizing _tursiops truncatus_ (the common bottlenose
dolphin) swimming within a volume equipped with a network of acoustic sensors
(e.g., hydrophones). We assume that the receivers are stationary, with well-
known positions from prior surveys. Recognizing that these dolphins have
typical swim speeds of 1-4 m/s ([22]) and their whistles are around 7.4-17.2
kHz ([23]), we simulate a dolphin with the following state that is unknown _a
priori_ to the sensor network:

$\textbf{\emph{r}}=\begin{bmatrix}-5.23\\\ 5.28\\\
-15.00\end{bmatrix}m,\qquad\textbf{\emph{v}}=\begin{bmatrix}1.38\\\ 1.53\\\
0.22\end{bmatrix}m/s,\qquad f=15\,kHz$

When provided perfect measurements (i.e., no noise), we recover the dynamical
state (position and velocity) of the source to machine precision in less than
0.5 seconds for unknown $f$ ($N=7$ observations) and in less than 0.1 seconds
for known $f$ ($N=6$ observations). Due to the simplicity of the problem’s
mathematical structure, an approach similar to that of [24] may lead to an
exceptionally fast solver with execution time less than a millisecond per
problem instance. Implementation of such a speed-up was not pursued further
here as the current sub-second runtimes were fast enough for our proof-of-
concept purposes.

We also conducted a Monte Carlo analysis to study the affect of imperfect
receiver state knowledge and noisy frequency measurements. Figure 3 shows the
results from 1,000 numerical experiments for both a known ($N=6$) and unknown
($N=7$) whistle frequency. In each experiment, we apply additive, zero-mean
Gaussian noise to each of the observers’ position ($\sigma_{r}=1.5\,cm$) and
to each of the observed frequencies ($\sigma_{f}=0.1\,Hz$).

###  4.2 Example 2: Orbit determination with moving receivers

The second experiment consists of the initial orbit determination (IOD) of a
satellite emitting an electromagnetic signal. One receiver is placed aboard a
satellite in a circular mid-Earth orbit (MEO) and all of the remaining
receivers are placed on the Earth’s surface (and rotate with it). We require
one orbiting receiver since the state observability is poor when all the
receivers are on the Earth’s surface. The orbital elements ([25]) for the
transmitting spacecraft are assumed to be

$a=12,000\,km,\qquad e=0.1,\qquad
i=20^{\circ},\qquad\Omega=200^{\circ},\qquad\omega=20^{\circ},\qquad\nu=0^{\circ}$

where $a$ is the semi-major axis, $e$ is the eccentricity, $\Omega$ is the
right ascension of the ascending node, $i$ is the inclination, $\omega$ is the
argument of the periapsis and $\nu$ is the true anomaly at epoch. These
orbital elements may be converted to three-dimensional position and velocity
using standard equations found in most modern astrodynamics textbooks ([25,
26, 27]). With perfect (noise free) receiver state knowledge and frequency
measurements, we are able to recover the 6-Degree-of-Freedom dynamical state
(and, hence, the six orbital elements) to nearly machine precision.

As with the underwater case, we may study how performance is affected by
imperfect knowledge of the receiver states and errors in the measured
frequency. Results from 1,000 numerical experiments are shown in Fig. 3 for
both a known ($N=6$) and unknown ($N=7$) transmit frequency. In each
experiment, we apply additive, zero-mean Gaussian noise to each of the
observers’ position and velocity and to the observed frequencies
($\sigma_{f}=0.5\,Hz$). Specifically, we used $\sigma_{r}=5\,cm$ and
$\sigma_{v}=1\,mm/s$ for the receivers on the Earth’s surface, and
$\sigma_{r}=1\,m$ and $\sigma_{v}=2\,cm/s$ for the receiver on orbit.

Figure 2: State errors in position (left column) and velocity (right column)
for example of a vocalizing dolphin (underwater application). Frequency
estimate errors are shown at the bottom for the case of $N=7$ receivers.
Results are from a 1,000 trial Monte Carlo analysis.

Figure 3: Orbital element errors for example of initial orbit determination of
a transmitting satellite (space application). Frequency estimate errors are
shown at the bottom for the case of $N=7$ receivers. Results are from a 1,000
trial Monte Carlo analysis.

##  5 Conclusions

The position and velocity of an object emitting an acoustic or electromagnetic
signal may be inferred from only the apparent frequency as measured by a
multitude of known receivers. This state estimation problem may be written as
a system of polynomials, which permits an efficient solution via homotopy
continuation. The frequency at the source may be either known or unknown. We
highlight the importance of using parameter homotopy, instead of the total
degree homotopy, and provide a strategy to filter out multiple solutions.

Two example applications show the feasibility of this approach, both in the
perfect (noise-free) and imperfect (noisy) scenarios. For the imperfect case,
we model errors in the receiver states (both position and velocity) and noise
is applied to the measured frequency. The proposed method provided a noisy
estimate of the transmitter state consistent with the amount of noise in the
measurements, both in the case of known and unknown frequency.

##  Appendix A Example Code

This appendix provides example code using Macaulay2, which is a popular
software system for solving polynomial systems and other problems in algebraic
geometry. For more information on Macaulay2, see [18].

The following code example shows how to build a polynomial system based on
multiple instantiations of Eq. (5).

⬇

\-------

\-- import the necessary packages

\-------

needsPackage "NumericalAlgebraicGeometry"

needsPackage "MonodromySolver"

\-------

\-- declare unknowns (corresponding to attributes of freq. source)

\-------

variableMatrix = gateMatrix{toList vars(x_1..x_3,u_1..u_3,f)}

\-- unknown source position (3D vector)

rr = matrix{{x_1},{x_2},{x_3}}

\-- unknown source velocity (3D vector)

vv = matrix{{u_1},{u_2},{u_3}}

\-- unknown source transmit frequency

ff = matrix{{f}}

\-------

\-- declare the parameters of the system

\-------

\-- fix the number of receivers for the unknown frequency case

nObs=7;

\-- receivers’ positions (3 x nObs matrix)

posObservers = gateMatrix apply(3, i -> apply(nObs, j -> posi_(i,j)))

\-- receivers’ velocities (3 x nObs matrix)

velObservers = gateMatrix apply(3, i -> apply(nObs,j -> veli_(i,j)))

\-- measured frquencies (nObs-D vector)

freqObservers = gateMatrix apply(nObs, i -> {freqi_(i)})

\-- speed of the signal

ci = gateMatrix{{signalSpeed_(i)}}

\-- variable of all the parameters of the system (1 x 50 vector)

parameterMatrix = gateMatrix{flatten entries(posObservers^{0,1,2})} | gateMatrix{flatten entries(velObservers^{0,1,2}) | flatten entries freqObservers | flatten entries ci}

\-------

\-- build the system of equations

\-------

fsquared = gateMatrix(ff*ff);

csquared = gateMatrix(ci*ci);

DPconstraints = flatten apply(nObs, i-> (

rri := gateMatrix posObservers_{i};

vvi := gateMatrix velObservers_{i};

ffi := gateMatrix{freqObservers#i};

dri := rri - rr;

dvi := vvi - vv;

dfi := ff-ffi;

drdv := transpose dri * dvi;

\-- Eq. (5) in the text for the i-th receiver

csquared * (transpose dri * dri) * (dfi * dfi)- (drdv * drdv) * fsquared

));

\-- create the system described by the given parameters, unknowns and
constraints

G = gateSystem(parameterMatrix, variableMatrix, gateMatrix DPconstraints)

Given the polynomial system G produced by the above code, we may write another
script to track the solutions using parameter homotopy continuation. It is
assumed that the parameters describing the start and target system have
already been loaded in the $1\times(7\,\texttt{nObs}+1)$ vectors
startParameters and targetParameters. The same holds for the solutions of the
start system, contained in the matrix startSolutions, of dimension (number of
unknowns $\times$ number of solutions). In practice, we obtain startSolutions
using MonodromySolver, a package of Macaulay2 based on [19].

⬇

\-------

\-- track the solutions

\-------

\-- declare the continuation parameter t, where t=0 for the start system and
t=1 for the target system

declareVariable t;

\-- define the homotopy given by the chosen start system and target system:
see Eq. (6) for reference

H = sub(gateMatrix G, parameters G,

(1-t) * startParameters + t * targetParameters);

Hgate = gateHomotopy(H, vars G, t);

\-- track the paths from t=0 to t=1

targetSolutions = trackHomotopy(Hgate, startSolutions);

## Acknowledgment

The authors thank Rickey Huang for valuable discussions related to the number
of polynomial solutions and the implementation of Macaulay2 code. We also
thank Tim Duff for assistance with developing the monodromy code used to
acquire the solutions to our start system.

## Competing insterests declaration

Competing interests: the authors declare none.

## References

  * [1] C. A. Doppler and F. J. Studnica, “Ueber das farbige licht der doppelsterne und einiger anderer gestirne des himmels.,” In Commission bei Borrosch & André, 1842.
  * [2] D. D. Nolte, “The fall and rise of the Doppler effect,” Physics Today, vol. 73, pp. 30–35, 03 2020.
  * [3] J. D. Kaunitz, “The doppler effect: A century from red shift to red spot,” Digestive Diseases and Sciences, vol. 61, pp. 340–341, 2016.
  * [4] R. Lambourne, “The doppler effect in astronomy,” Physics Education, vol. 32, p. 34, jan 1997.
  * [5] W. H. Guier and G. C. Weiffenbach, “Theoretical analysis of doppler radio signals from earth satellites,” Nature, vol. 181, pp. 1525–1526, 1958\.
  * [6] D. Gillespie, L. Palmer, J. Macaulay, C. Sparling, and G. Hastie, “Passive acoustic methods for tracking the 3d movements of small cetaceans around marine structures,” PLoS ONE, 2020.
  * [7] E. Felemban, F. K. Shaikh, U. M. Qureshi, A. A. Sheikh, and S. B. Qaisar, “Underwater sensor network applications: A comprehensive survey,” International Journal of Distributed Sensor Networks, vol. 11, no. 11, p. 896832, 2015.
  * [8] D. Amrita and D. Mou, “On accurate localization of sensor nodes in underwater sensor networks: a doppler shift and modified genetic algorithm based localization technique,” Evolutionary Intelligence, vol. 14, no. 11, pp. 119–131, 2021.
  * [9] Z. Gong, C. Li, F. Jiang, and J. Zheng, “Auv-aided localization of underwater acoustic devices based on doppler shift measurements,” IEEE Transactions on Wireless Communications, vol. 19, no. 4, pp. 2226–2239, 2020\.
  * [10] Y.-T. Chan and F. Jardine, “Target localization and tracking from doppler-shift measurements,” IEEE Journal of Oceanic Engineering, vol. 15, no. 3, pp. 251–257, 1990.
  * [11] A. Famili and J.-M. J. Park, “Rolatin: Robust localization and tracking for indoor navigation of drones,” in 2020 IEEE Wireless Communications and Networking Conference (WCNC), pp. 1–6, 2020.
  * [12] R. B. Patton, “Orbit determination from single pass doppler observations,” IRE Transactions on Military Electronics, vol. MIL-4, no. 2/3, pp. 336–344, 1960.
  * [13] J. Christian, C. Ertl, K. Horneman, and A. Lovell, “Doppler-only initial orbit determination for an orbiting transmitter,” 01 2022.
  * [14] C. A. Ertl, Solutions and Analysis of Multivariate Polynomial Systems for Geolocation and Initial Orbit Determination.  Phd dissertation, Department of Mechanical, Aerospace and Nuclear Engineering. Rensselaer Polytechnic Institute, Troy, New York, April 2023.
  * [15] A. Morgan, Solving polynomial systems using continuation for engineering and scientific problems.  Englewood Cliffs, NJ: Prentice Hall Inc., 1987.
  * [16] A. J. Sommese and C. W. Wampler, II, The numerical solution of systems of polynomials.  World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2005.
  * [17] A. Leykin, “Numerical algebraic geometry,” J. Softw. Algebra Geom., vol. 3, pp. 5–10, 2011.
  * [18] D. R. Grayson and M. E. Stillman, “Macaulay2, a software system for research in algebraic geometry.” Available at http://www.math.uiuc.edu/Macaulay2/, 2023.
  * [19] T. Duff, C. Hill, A. Jensen, K. Lee, A. Leykin, and J. Sommars, “Solving polynomial systems via homotopy continuation and monodromy,” IMA Journal of Numerical Analysis, vol. 39, 09 2016.
  * [20] J. Harris, Algebraic Geometry. A First Course, vol. 133 of Graduate Texts in Mathematics.  New York, NY: Springer, 1992.
  * [21] D. A. Cox, J. Little, and D. O’Shea, Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra.  Undergraduate Texts in Mathematics, Cham: Springer, 4 ed., 2015.
  * [22] B. Wursig and M. Wursig, “Behavior and ecology of the bottlenose dolphin (tursiops truncatus) in the south atlantic,” Fishery Bulletin, vol. 77, no. 2, pp. 399–412, 1979.
  * [23] J. N. Oswald, J. Barlow, and T. F. Norris, “Acoustic identification of nine delphinid species in the eastern tropical pacific ocean,” Marine Mammal Science, vol. 19, no. 1, pp. 20–037, 2003.
  * [24] P. Hruby, T. Duff, A. Leykin, and T. Pajdla, “Learning to solve hard minimal problems,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 5532–5542, 2022.
  * [25] R. Bate, D. Mueller, J. White, and W. Saylor, Fundamentals of Astrodynamics, 2nd Ed. Mineola, NY: Dover, 2020.
  * [26] D. Vallado, Fundamentals of Astrodynamics and Applications, 3rd Ed. Hawthorne, CA: Microcosm Press, 2007.
  * [27] H. Curtis, Orbital Mechanics for Engineering Students, 4th Ed. Cambridge, MA: Butterworth-Heinemann, 2020.