\documentclass[12pt,letterpaper]{article}
%graphics & geometry
%%Check if we are running under latex or pdflatex
\ifx\pdftexversion\undefined
\usepackage[dvips]{graphicx}
\usepackage[dvips,usenames]{color}
\usepackage[letterpaper,margin=1in,nohead,dvips]{geometry}
\else
\usepackage[pdftex]{graphicx}
\usepackage[pdftex,dvipsnames,usenames]{color}
\usepackage[letterpaper,margin=1in,nohead,pdftex]{geometry}
\fi
%refs
\usepackage[round]{natbib}
\bibliographystyle{plainnat}
%language
\usepackage[british]{babel}
\usepackage{t1enc}
%font
\usepackage{mathptm}
\usepackage{xspace}
\renewcommand{\sfdefault}{phv}
%math
\usepackage{amsmath,amssymb}
\usepackage{array}
\renewcommand{\boldsymbol}[1]{\pmb{#1}}
\newcommand{\diff}[1][x]%
{\mathchoice{\frac{d#1}{dt}}%
{d#1/dt}{d#1/dt}{d#1/dt}}
%figures
%\newenvironment{oldfigureref}[2]%
%{\begin{figure}[#2]%
%\renewcommand{\thefigure}{\ref{#1}\addtocounter{figure}{-1}}}%
%{\makeatletter%
%\renewcommand \thefigure %
%{\@arabic\c@figure}%
%\makeatother%
%\end{figure}}
%shortcuts
\newcommand{\ie}{i.e.\@\xspace}
\newcommand{\eg}{e.g.\@\xspace}
\newcommand{\etc}{etc.\@\xspace}
\newcommand{\eqn}{Eqn.\@\xspace}
\newcommand{\eqns}{Eqns.\@\xspace}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\e}[1]{\ensuremath{{\cdot\!10^{#1}}}\xspace}
\begin{document}
\title{
\textsc{Linear stability analysis of locust model}}
\author{Authors\footnote{}
}

\maketitle

\section{The Model Equations}

\begin{subequations}\label{LocustModelEqs}
\begin{eqnarray}
s_t+ \nabla\cdot (v_s s) &=& -f_2(\rho) s +f_1(\rho)g\\
g_t+ \nabla\cdot (v_g g) &=& f_2(\rho) s -f_1(\rho)g
\end{eqnarray}
\end{subequations}
where
$s(x,t), g(x,t)$ are solitarious and gregarious locusts, and $\rho(x,t) = s(x,t) +g(x,t)$
is the total local density.
The motion of the locusts is governed by Morse Potential forces, with corresponding potentials $Q_s,Q_g$ such that 

\begin{equation}
v_s=-\nabla (Q_s*\rho),\quad v_g=-\nabla (Q_g*\rho)
\end{equation}
and
\begin{equation}\label{Eq:potentials}
Q_s= R_s e^{-|x|}, \quad Q_g=R_g e^{-|x|/r_g} -A_g e^{-|x|/r_a}
\end{equation}
This represents the fact that solitarious locusts repel from other locusts, whereas
gregarious locusts are attracted to others (with possibly some hard-core repulsion
or repulsion at very close distance).

We also model the exchange between the locusts using the functions
\begin{subequations}
\begin{eqnarray}
f_1(\rho)&=& \frac{1}{1+\rho^2}\\
f_2(\rho)&=& \delta_2\frac{(\rho/k_2)^2}{1+(\rho/k_2)^2}
\end{eqnarray}
\end{subequations}

We will consider a 1D version and assume that there exist(s) (a) Homogeneous 
steady state solution, such that
$\rho_0, s_0, g_0$ are constant, $\rho_0=s_0+g_0$.

\section{A few facts}

Here we concentrate a few facts that make the subsequent calculations clearer. Most of these can be easily established by the appropriate manipulations.

\begin{itemize}

\item It does not seem easy to solve for the homogenous steady states explicitly. However, we can say that $f_1$ is a strictly decreasing function of $\rho$ and $f_2$ is a strictly increasing function of $\rho$. Thus we know that the derivatives $f_1'(\rho), f_2'(\rho)$ are, respectively, negative and positive real values at any steady state.

\item There is a family of homogeneous steady states parameterized by the total density $\rho$. In particular, one can choose any value $\rho$ and have the steady state
$$s = \frac{\delta_2 \rho^3 (1+\rho^2)}{k_2^2 + \rho^2 + \delta_2 \rho^2 +\delta_2 \rho^4}, \qquad g = \frac{\rho (k_2^2+\rho^2)}{k_2^2 + \rho^2 + \delta_2 \rho^2 +\delta_2 \rho^4}$$
In the small $\rho$ limit,
$$ s \approx \frac{\delta_2}{k_2^2} \rho^3, \qquad g \approx \rho - \frac{\delta_2}{k_2^2} \rho^3$$
so that $g \gg s$.
In the large $\rho$ limit,
$$ s \approx \rho - \frac{1}{\delta_2} \frac{1}{\rho}, \qquad g \approx \frac{1}{\delta_2} \frac{1}{\rho}$$
so that $s \gg g$.
Chad thinks that the small $\rho$ HSS should be unstable for reasonable parameter choices (since a primarily gregarious population at low density shouldn't persist as such) and also that the large $\rho$ HSS should be unstable for reasonable parameter choices (since a primarily solitarious population at high density shouldn't persist as such).

\item Assuming either an infinite domain with no locusts ``far away'' or possibly no-flux boundaries on a finite domain leads to the conservation condition
$$M=\int_{-\infty}^{\infty} [s(x,t)+g(x,t)] dx=\int_{-\infty}^{\infty} \rho(x,t) dx.$$
(We expect that the total number of locusts is fixed, as they are merely redistributing and switching phase, but not entering, leaving nor born/die.)

\item The velocities can be expressed either via the force or the potential kernels, in two equivalent expressions i.e. 
\begin{equation}
v_s=F_s*\rho=-\nabla (Q_s*\rho),\quad v_g=F_g*\rho=-\nabla (Q_g*\rho)
\end{equation}
where $F_i=-\nabla Q_i$. 
This follows from the fact (e.g. in 1D) that 
$$F_s*\rho=\int_{-\infty}^{\infty} F_s(x-z)\rho(z) dz= \int_{-\infty}^{\infty}(-\nabla Q_s(x-z))\rho(z) dz=-\nabla\int_{-\infty}^{\infty}Q_s(x-z)\rho(z) dz = -\nabla (Q_s*\rho)$$
where we have taken the derivative w.r.t. $x$ outside of the integral w.r.t. $z$. 


\item For nonlocal attraction-repulsion forces, the kernel $F_i$ has to be an odd function to depict the fact that there is antisymmetry in the forces induced due to a neighbor to the left and to the right. This means that when $\rho(x,t)=\rho_0$=constant, we have that

$$v_s=F_s*\rho_0 =-\nabla (Q_s*\rho_0)=0$$
This says that when the density is constant, there is no net flow of any locusts, and it follows from the fact that

$$F_s*\rho_0=\int_{-\infty}^{\infty} F_s(x-z)\rho_0 dz =\rho_0 \int_{-\infty}^{\infty} F_s(y) dy=0,$$
the latter being an integral of an odd function over the infinite domain.

\item
The following convolution is simply related to the Fourier transform of the kernel $Q$:
$$Q*e^{ikx}= {\hat Q}(k) e^{ikx}.$$
This follows from writing out the convolution integral and some simplifications.

\item
Based on the above, we also compute that
$$\frac{\partial^2}{\partial x^2} [Q*e^{ikx}]=-k^2 {\hat Q}(k) e^{ikx}.$$
This follows by elementary differentiation of the previous formula.

\item The potentials  in Eqn.~\ref{Eq:potentials} are made up of even functions of the form $Q(x)=A e^{|x/a|}$ (as well as sums or differences of such exponentials). The Fourier transform of such a function is found to be
$${\hat Q(k)} = 2\frac{Aa}{1+(ak)^2}. $$
This result follows by direct calculation of 
$$\int_{-\infty}^{\infty} Q(x)e^{-ikx}dx=2A\int_{0}^{\infty} e^{x/a}\cos(kx)dx$$
with the usual two integration by parts. See also Leverentz et al (2009) or any table of Fourier transforms for this standard result.

\item Based on the above, the Fourier transforms of the potentials  in Eqn.~\ref{Eq:potentials}
are:
\begin{equation}\label{Eq:FTofpotentials}
{\hat Q_s(k)}= 2\frac{R_s}{1+k^2}, \quad {\hat Q_g(k)}=2\left( \frac{R_gr_g}{1+(r_gk)^2}-
\frac{A_gr_a}{1+(r_ak)^2}\right)
\end{equation}
In particular, we observe that ${\hat Q_s(k)}>0$ for all $k$ is a monotonic decreasing function (since solitarious locusts are assumed to have pure repulsion) whereas there exist values of $k$ for which ${\hat Q_g(k)}<0$ because gregarious locusts are in a regime where, left on their own, they would tend to collapse to a cluster. 

\end{itemize}

\section{Analysis of growth of small perturbations of the model}

We now consider the behaviour of small perturbations of Eqs.~\ref{LocustModelEqs}.
Specifically, we assume that
\begin{subequations}\label{Perturbations}
\begin{eqnarray}
s(x,t)&=&s_0+s'(x,t)=s_0+\mathcal{S}e^{ikx}e^{\sigma t}\\
g(x,t)&=&g_0+g'(x,t)=g_0+\mathcal{G}e^{ikx}e^{\sigma t}
\end{eqnarray}
\end{subequations}
where $\mathcal{S}, \mathcal{G}$ are (small) amplitudes, and $k$ is the wavenumber. 
We ask whether there are conditions under which $\sigma >0$ signifying growth of these spatially heterogeneous perturbations.
We also define
$$\rho(x,t)=\rho_0+\rho'(x,t) =s_0+g_0+ s'(x,t)+g'(x,t)=\rho_0+\mathcal{P}e^{ikx}e^{\sigma t}$$
Then by conservation, it follows that
$$\int_{-\infty}^{\infty} [s'(x,t)+g'(x,t)] dx=0.$$
[This seems to indicate that we are not free to impose arbitrary perturbations, as we have to ensure mass conservation.]

We substitute the perturbations \ref{Perturbations} into Eqs.~\ref{LocustModelEqs}. 
Let us observe first that when doing so, the only surviving linearized convection term in the equation for $s$ has the form
$$\frac{\partial}{\partial x}\left(s_0\frac{\partial}{\partial x}[Q_s*\rho'] \right)=s_0
\frac{\partial^2}{\partial x^2}[Q_s*\rho'].$$
This follows from the fact that $Q_s*\rho_0=0$ by facts in a previous section. (We also dropped the small nonlinear term that has both $s'$ and $\rho'$.)

For the linearization step, we observe that 
\begin{eqnarray}
f_1(\rho)g-f_2(\rho)s&=&f_1(\rho_0+\rho')(g_0+g')-f_2(\rho_0+\rho')(s_0+s') \nonumber\\
&\approx& [f_1(\rho_0)g_0 -f_2(\rho_0)s_0] + f_1(\rho_0)g'-f_2(\rho_0)s' +
f_1'(\rho_0)\rho'g_0 - f_2'(\rho_0)\rho's_0 \nonumber\\
&=& -\alpha s' +\beta g' \nonumber\\
\mbox{where  }~\alpha&=&f_2(\rho_0)+f_2'(\rho_0)s_0-f_1'(\rho_0)g_0 \\
\beta&=& f_1(\rho_0)-f_2'(\rho_0)s_0+f_1'(\rho_0)g_0
\end{eqnarray}
In the above, we have used the fact that the terms in square brackets are zero as they are steady state terms. We have also collected all other terms as coefficients of the $s'$ and $g'$ terms, and neglected any higher order terms such as $s'g', etc$.

The equations (in 1D) are then 
\begin{subequations}\label{LinearizedEqs}
\begin{eqnarray}
s'_t- s_0\frac{\partial^2}{\partial x^2}[Q_s*\rho'] &=& -\alpha s' +\beta g'\\
g'_t- g_0\frac{\partial^2}{\partial x^2}[Q_g*\rho'] &=& \alpha s' -\beta g'
\end{eqnarray}
\end{subequations}
%where $f'_i(\rho_0)=(d f_i(\rho)/d\rho)|_{\rho=\rho_0}$ are derivatives of the kinetic functions w.r.t. the density, evaluated at the HSS. 
%(From Chad: I could be wrong about this, but in the equations above, should the first plus signs -- after the time derivative terms -- be minus signs because the velocity is minus the gradient of Q? If so, then this sign change would propagate below, I imagine...)

Using the explicit forms of the perturbations and the facts previously established, we obtain

\begin{subequations}\label{PertEqs}
\begin{eqnarray}
\sigma\mathcal{S}+ s_0(\mathcal{S+G})k^2 {\hat Q_s} &=& -\alpha \mathcal{S} +\beta \mathcal{G}\\
\sigma\mathcal{G}+ g_0(\mathcal{S+G})k^2 {\hat Q_g} &=& \alpha  \mathcal{S} -\beta \mathcal{G}
\end{eqnarray}
\end{subequations}

These can be written in matrix form
\[ \left( \begin{array}{cc}
\sigma +s_0 k^2 {\hat Q_s} + \alpha &  s_0 k^2 {\hat Q_s}-\beta \\
g_0 k^2 {\hat Q_g}-\alpha & \sigma +g_0 k^2 {\hat Q_g} + \beta \end{array} \right)
\left( \begin{array}{c}
\mathcal{S}   \\
\mathcal{G}  \end{array} \right)
=
\left( \begin{array}{c}
0   \\
0  \end{array} \right)
\] 
For nontrivial solutions, the determinant of the above matrix must be zero. This gives rise to the characteristic equation

\[
\left(\sigma +s_0 k^2 {\hat Q_s} + \alpha \right)\left( \sigma +g_0 k^2 {\hat Q_g} + \beta \right)- (s_0 k^2 {\hat Q_s}-\beta )( g_0 k^2 {\hat Q_g}-\alpha)=0
\]
Simplification leads to the form
\[
\sigma^2 +B\sigma +C=0.
\]
where
$$B=B(k)=k^2[s_0{\hat Q_s}+g_0{\hat Q_g}]+(\alpha+\beta),$$
$$C=C(k)=(s_0 k^2 {\hat Q_s} + \alpha)(g_0 k^2 {\hat Q_g} + \beta)-
(s_0 k^2 {\hat Q_s}-\beta )( g_0 k^2 {\hat Q_g}-\alpha)$$
The expression in $C$ is of the form
$$(a+\alpha)(b+\beta)-(a-\beta)(b-\alpha) $$
Where $a=s_0 k^2 {\hat Q_s}, b=g_0 k^2 {\hat Q_g}$. Expanding the products and simplifying leads to
$$C= (ab+a\beta+b\alpha+\alpha\beta)-(ab-\beta b-\alpha a+\alpha\beta)=a \beta+b \beta+a \alpha+b \alpha=(a+b)(\alpha+\beta)$$
i.e. we have found (collecting both coefficients here) that 
$$B=B(k)=k^2[s_0{\hat Q_s}+g_0{\hat Q_g}]+(\alpha+\beta),$$
$$C=C(k)=k^2(s_0 {\hat Q_s}+g_0{\hat Q_g})(\alpha+\beta)$$

We now observe that $B(k)$ is the sum of two terms and $C(k)$ is the product of the same two terms, which implies that the characteristic equation factors as follows:
$$
\sigma^2 +B\sigma +C=(\sigma +k^2(s_0 {\hat Q_s}+g_0{\hat Q_g}))(\sigma + \alpha+\beta)=0$$
We have thus found the two eigenvalues explicitly, namely
$$\sigma_1=-k^2(s_0 {\hat Q_s}+g_0{\hat Q_g}),\quad \sigma_2=-(\alpha+\beta)$$
Using the formulas for the coefficients $\alpha,\beta$ leads to
\begin{eqnarray}
\sigma_2&=&-(f_2(\rho_0)+f_2'(\rho_0)s_0-f_1'(\rho_0)g_0+ f_1(\rho_0)-f_2'(\rho_0)s_0+f_1'(\rho_0)g_0)\\
&=&-(f_2(\rho_0)+f_1(\rho_0))
\end{eqnarray}



We now observe that in the well-mixed case, when there is no spatial redistribution, we want the homogeneous steady state to be stable. This is equivalent to saying that $\sigma_2<0$. Note that in that same case, $k=0$, so a second zero eigenvalue is obtained, as expected from the mass conservation condition. 

In order for instability to spatially heterogeneous perturbations, we require that $\sigma_1>0$. Hence the condition for aggregation is:

$$k^2(s_0 {\hat Q_s}+g_0{\hat Q_g})<0$$
Since $k$ is real, this condition is equivalent to
$$(s_0 {\hat Q_s}+g_0{\hat Q_g})<0.$$

We point out that this condition is remarkably analogous to the finding in Section 2.2 of Leverentz et al's (2009) calculation for the stability of a single species swarm.
We also note that this condition does not depend on the nature of the functions $f_1, f_2$ beyond the steady state swarm densities that they determine. It thus depends only on the attraction-repulsion potentials and the relative amounts of solitary and gregarious locusts.

\section{How aggregation could occur}

In view of the above, let us initially consider only the redistribution dynamics, i.e., take $f_1=f_2=0$ and assume that $s_0, g_0, \rho_0=s_0+g_0$ are set by some arbitrary mechanism, e.g. by introducing a fixed number of each type of locust into an arena where they interact with each other by attraction/repulsion.
Let us define the parameter $\phi=g_0/\rho_0, 1-\phi = s_0/\rho_0$, such that $0 \le \phi \le 1$. We are interested in an initial population consisting mainly of solitarious locusts, where a few gregarious locusts are introduced (small $\phi$). We ask how many gregarious locusts are needed in order to lead to an aggregation instability.
Dividing by $N$, restate the instability condition as
$$(1-\phi){\hat Q_s}(k)+\phi{\hat Q_g}(k)<0.$$
then in the case $\phi=0$ it is clear that no instability is possible, since we have established that ${\hat Q_s}(k)>0$ for all $k$. We also know that for $\phi=1$ there exist values of $k$
such that the inequality is satisfied. Let $k^*$ be the wavenumber corresponding to the minimum value of ${\hat Q_g}(k)$ (which is negative, as argued above). Then we can adjust $\phi$ such that at that $k^*$ and $\phi_{crit}$ the above expression is an exact equality, and at all other $k$ values, the inequality is violated. Now consider a slightly larger value of $\phi$ (slightly more gregarious locusts). This ($\phi>\phi_{crit}$) will result in a range of wavenumbers satisfying the instability condition, and aggregation will ensue.

Now let us return to the full dynamics, and consider some initial homogeneous steady state $s_0,g_0$ that is determined by the dynamics of switching with $f_1, f_2$ nonzero. Then according to the above, aggregation will occur only once $g_0$ is sufficiently large relative to $s_0$. This could result either from some noise that leads by chance to a high enough local perturbation of $g_0$, or by a bistable mechanism that drives the population into a second steady state with relatively large $g_0$ value (relative to $s_0$).

Base on these ideas it should be relatively easy to construct situations that do or do not lead to an aggregation of locusts, either by adjusting the kinetic functions or by directly manipulating the relative levels of each type of locust.
\end{document}