Serena Dipierro & Enrico ValdinociUniversity of Western Australia,
35 Stirling Highway,
Crawley WA 6009, Australia<EMAIL_ADDRESS>
Elliptic partial differential equations

from an elementary viewpoint

CHAPTER: PREFACE

more to life than mathematics”,
Joan said. “But not much more”.

Greg Egan, Glory.

These notes are the outcome of some courses taught to undergraduate and graduate students
from the University of Western Australia,
the Pontifícia Universidade Católica do Rio de Janeiro and
the Indian Institute of Technology Gandhinagar
in 2021.

Far from aiming at being all-encompassing, the following pages wish to
shed some light on a number of selected topics in the theory
of elliptic partial differential equations with a style that should be
accessible to third-year undergraduate students, possibly under an
inspired mentorship, but might
also provide some interest to more advanced students
and possibly professional researchers. While all the
topics presented are of classical flavor,
the exposition and the way the material are organized is
perhaps rather original, with the intention of addressing several quite
difficult points with a style that is as self-contained as possible,
rigorous as well as intuitive,
and approachable without major prerequisites (indeed,
we only assume prior knowledge of the “basic” analysis, giving also
reference to useful results and theorems whenever the less advanced readers
may need to consolidate their backgrounds).

We stress again that the list of topics covered here is far from being exhaustive,
since we mainly dealt with classical problems related to the Laplace operator
and their natural counterpart for equations in nondivergence form
(the divergence structure case is not really covered here,
though the Divergence Theorem is obviously utilized ubiquitously, the variational structure of the Laplace-Beltrami
operator is discussed quite in detail
and several pages are devoted to the theory of capacity which possesses a variational essence).
Also, we do not address here nonlinear, singular/degenerate elliptic operators, fully nonlinear equations
and fractional/nonlocal elliptic equations.
We do not even linger too much on explicit solutions
and cheap tricks to find them (e.g., we do not repeat over and over very specific methods such as the separation of variables),
since, after all, it is very unlikely that one can solve explicitly[Moreover, explicit solutions
may provide a handy resource to develop an initial theory, but they usually do not exhaust the complexity
of the problem. For example, according to <cit.>,
it is quite common to run into “one of the serious problems with such exact solutions [...]: namely, they often do not determine all possible solutions
and indeed, may not even give the most relevant one”.]
a partial differential equation,
instead it is often more useful to understand the qualitative and quantitative properties
of the solution without solving explicitly the corresponding equation.

The many parts of the theory of elliptic partial differential equations which are missing from these pages
are by no means less important than the topics covered: the issue is that we had to make
a rather harsh selection of topics just to be able to collect all the material into finitely many pages
(arguably, Jorge Luis Borges would remark that a comprehensive treatment of elliptic equations
can be found in the Library of Babel <cit.>).

Some additional fun for the interested reader could come from
the treatment of a few natural problems with geometric flavor, such as
the Soap Bubble Theorem and Serrin's overdetermined problem. To treat them,
but this is also a general approach used everywhere in the book,
we tried to provide different perspectives and different
points of view, also providing physical motivations
to trigger creativity and imagination, and favoring a dynamic interplay between geometric and analytic arguments.

These notes also collect a number of classical, albeit not always well known,
topics, which we believe can serve as an excellent training camp to develop
some familiarity with elliptic equations.

Here are
a couple of features of these notes that we hope can be appealing for many readers. First, all is built in an artisanal way, perhaps not by following the most elegant, concise or general approach, but aiming at a possibly slow, but conceptually clear, sequence of strategic steps. Second, this set of notes is planned in a way that one can read starting almost any place and freely jump from one topic to another just by following a personal stream of thoughts.

Though a list of classical and modern references is given at the end,
this is not a book focused on the history
of mathematics, hence we will not address topics such as priority
of discoveries or progress of subject over the course of time.

As usual, these notes may contain errors or inconsistencies: if you find
any, please let us know. In general, we will be happy to receive
comments and criticisms to possibly improve this work.
After all, scientific knowledge is based on a dynamic flow
of information and we
will certainly cherish and treasure readers' feedback and advice.

All right, enough chitchat, please fasten your PFD,
it's time to start our journey together.

Serena and Enrico

Courtesy of Burbuqe Shaqiri.

CHAPTER: MOTIVATIONS

§ THE LAPLACIAN COMING ON THE SCENE

Without aiming at being exhaustive, we present here some classical circumstances in which partial differential
equations, and especially those relying on the Laplace operator, naturally arise to model interesting natural phenomena.
The objective is not much to use mathematics  to find the Answer to the Ultimate Question of Life, the Universe, and Everything
(which is well known to be $42$ anyway) but rather to develop some familiarity with some of the basic features
of partial differential equations, to see them in action and to build an intuition and an instinctive feeling
of the problems,
which will turn out handy in our voyager through the mighty jungle of technical and seemingly abstract mathematics.

§.§ The heat equation

Among the many occurrences in which elliptic equations naturally arise in nature, one of the most widely popularized comes from the stationary state of the heat equation.
In its most basic formulation, the setting of this equation is due to
Jean-Baptiste Joseph Fourier who[Interestingly, while changing the course of mathematics and physics, Fourier was acting as a full time politician, since in 1801 Napoléon Bonaparte decided to appoint him as Prefect of the Department of Isère, in the Alps. And perhaps it was not a completely trivial task to upset the First Consul (Emperor of the French from 1804) and dedicate oneself to mathematical investigations; in any case, playing around
with the heat equation at that point was not really part of Fourier's workload model.

See Figure <ref>
for a watercolor caricature of Fourier (by Julien-Léopold Boilly).] presented for the first time a partial differential equation to describe conductive diffusion of heat (and, to study this equation, he also introduced a marvelous instrument that was going to change forever the History of the Universe, namely, the Fourier Series).

Caricature of Joseph Fourier (Public Domain image from

The ansatz[“Ansatz” (plural: “Ansätze”) is a masculine German name.
As quite customary in German, this name is obtained by a preposition (“an”,
meaning “at”, or “on”, or “to”, or “against”, or pretty much whatever one likes) and another name
(“Satz”, plural: “Sätze”, meaning “sentence”, or “statement”). And as quite customary in German,
a precise translation goes well beyond the expertise of modest mathematicians of our caliber.
The “ansatz” could be literally a “starting point”. In jargon, it is frequently used to denote
an “educated guess” about a problem,
typically an additional simplifying assumption set at the beginning of an argument,
possibly to be verified later on (if it produces any result worthy of further consideration).

Maybe mathematics is not quite a deductive discipline, after all. Perhaps it is mostly the art of making good guesses.
Of course, anybody can make guesses. But for a good guess one needs to become acquainted with the problem
and to have developed that sort of familiarity that permits to separate the essential aspects from
the minor details (and, in trying to do this, sometimes
a little bit of luck doesn't hurt).] made by Fourier is that, at some[The space variable $x$ here is supposed to be in $\R^n$. Though the models presented refer to physical spaces, in which one can focus on the case $n=3$ (or even $n=2$, when dealing with plates, or $n=1$, when dealing with ropes of negligible thickness), whenever possible we prefer to work in (Euclidean) spaces of arbitrary dimension. When one does that, the notation becomes usually more effective and the essential treats of the main ideas are very often more transparent. Also, when it will be needed to focus our attention on the special case of particular dimensions, one will be able to understand why and to what extent these dimensions are different from the general case. In dealing with the arbitrary dimension case, we can't help quoting the first “general working principle” in the Preface of <cit.>: “PDE theory is (mostly) not restricted to two independent variables. Many texts describe PDE as if functions of the two variables $(x, y)$ or $(x, t)$ were all that matter. This emphasis seems to me misleading, as modern discoveries concerning many types of equations, both linear and nonlinear, have allowed for the rigorous treatment of these in any number of dimensions. I also find it unsatisfactory to classify partial differential equations: this is possible in two variables, but creates the false impression that there is some kind of general and useful classification scheme available in general”.]
point $x\in\R^n$, at an instant of time $t\in\R$, the variation of temperature $u(x,t)$ of some body is determined by the heat flow around the point under consideration, plus possibly additional heat[For the sake of precision, we should probably say here “heating or cooling sources”,
e.g., fireplaces, bonfires, stoves, radiators, gas heaters, infrared heaters, but also refrigerators, air conditioners, swamp boxes, chilled beams, or whatever heating or cooling device Fourier might have possessed at that time.
With this understanding, positive values of $f$ would correspond to heating systems and negative values to cooling systems.]
sources. Let us denote by $B(x,t)$ the heat
flux vector and by $f(x,t,u(x,t))$ the scalar intensity of the external[To make things simpler, one can suppose that $f$ depends only on $x$, as a heat source which is permanently switched on or off, or on $x$ and $t$, in case the source is turned on and off or its intensity is modified over time. The additional dependence on $u$ itself models the interesting case of a thermostat controlling a heating or cooling system.] heat sources.

A precise measure of the temperature precisely at a point $x$ is certainly rather difficult from a practical point of view, so it is useful sometimes to consider instead an average[As usual, we use
here (and repeatedly over this set of notes) the integral notation for averagesaverage, that is

$$ \fint_\Omega u(x,t)\,dx:=\frac1{|\Omega|}\,\int_\Omega u(x,t)\,dx,$$

being $|\Omega|$ the Lebesgue measure of $\Omega$.]
temperature, measured in some region of the space $\Omega$,

$$ U(t;\Omega):=\fint_\Omega u(x,t)\,dx.$$

In this setting, up to physical constants that we omit, we can equate the variation of $U$
with the heat flux through $\partial\Omega$, possibly adding to it the effect of the heat sources (if any) in $\Omega$,
namely we write that

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj}
\partial_t U(t;\Omega)=-\frac1{|\Omega|}\int_{\partial\Omega} B(x,t)\cdot\nu(x)\,d{\mathcal{H}}^{n-1}_x+\fint_\Omega f(x,t,u(x,t))\,dx.\end{equation}

In terms of notation, we are denoting here (and we will mostly keep this notation throughout the all set of notes)
by $\nu(x)$ the unit normal vector at $x\in\partial\Omega$ pointing outwards from $\Omega$ and
by $d{\mathcal{H}}^{n-1}_x$ the surface element[In many textbooks, the surface
elementsurface element is denoted by $ds$, or $dS$, or $d\Sigma$. The
notation $d{\mathcal{H}}^{n-1}_x$ comes from the fact that we are denoting by ${\mathcal{H}}^k$ the $k$-dimensional
Hausdorff measure, and, for smooth $(n-1)$-dimensional surfaces, the surface measure itself is precisely
equal to ${\mathcal{H}}^{n-1}$. See e.g. <cit.> for a thorough presentation of the Hausdorff measure.

The notation $d{\mathcal{H}}^{n-1}_x$ used here is perhaps a bit heavier than $ds$, or $dS$, or $d\Sigma$, or similar ones, but it has the benefit of being clearer and more explicit. Also, it allows us to consider surface integrals along less regular objects without having to introduce a new notation on a case-by-case basis.] on $\partial\Omega$. The minus sign appearing on the right hand side of (<ref>) is due to the fact that the normal $\nu$ points towards the exterior, while we are computing there the flux coming into the region $\Omega$.

By differentiating under the integral sign[In this chapter about motivations, the arguments are developed at a formal level, we feel free to exchange derivatives and integrals, we do not keep track of lower order terms in the expansions, we do not discuss convergence issues and existence of limits, etc. This is a rather customary approach when one deals with providing convincing, but not necessarily circumstantial, motivations for a problem with the objective of developing some intuition about it. We promise to try to be more rigorous from the next chapter on.] in (<ref>), we find that

\begin{equation*}\fint_\Omega\partial_t u(x,t)\,dx=-\frac1{|\Omega|}\int_{\partial\Omega} B(x,t)\cdot\nu(x)\,d{\mathcal{H}}^{n-1}_x+\fint_\Omega f(x,t,u(x,t))\,dx.\end{equation*}

Hence, by the Divergence[Typically reserve the notation “$\nabla$” for the vectors of the derivatives with respect to the space variables and “$\div$” for the corresponding divergence. No confusion should arise with respect to derivatives with respect to the time variable, which is usually denoted here by “$\partial_t$”.]

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj2}\fint_\Omega\partial_t u(x,t)\,dx=-\fint_{\Omega} \div B(x,t)\,dx+\fint_\Omega f(x,t,u(x,t))\,dx.\end{equation}

From this, Fourier took a further step by realizing that to make this identity manageable one needs to
take a constitutive law about the flux vector $B$, relating it to the temperature. Fourier's new ansatz was to suppose that the
flux of heat between two adjacent (infinitesimal) regions is proportional to the (infinitesimal) difference of their temperatures, namely

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj3}
B(x,t)=-\kappa(x,t)\,\nabla u(x,t).
\end{equation}

The proportional coefficient $\kappa$ is sometimes called “heat conduction coefficient”heat conduction coefficient. We take $\kappa$ to be positive:
with respect to this, we stress that the minus sign in (<ref>) is motivated by the fact that the heat flows from hotter regions to colder ones (hence in the opposite direction of the growth of the temperature function $u$).
The relation in (<ref>) is also sometimes[Actually, the same presentation
here could have been done to describe the transport of mass through diffusive means.
In this setting, one argue that the flux goes from regions of high concentration to regions of low concentration.
The ansatz corresponding to (<ref>) would be that the
magnitude of the flux is proportional to the concentration gradient.
In the context of mass diffusion, $\kappa$ is sometimes called “diffusion coefficient”
and (<ref>) is sometimes refereed to with the name of
“Fick's Law”Fick's Law after the German physician and physiologist
Adolf Eugen Fick.

Alternative choices for the constitutive relation are possible. For instance, one can replace (<ref>) with a more general equation in which the flux depends possibly in a nonlinear way
from the gradient of the temperature (in the case of Fourier's model, or of the transported mass in the case of Fick's model),
that is one could suppose that $B=-\Phi(\nabla u)$, for some function $\Phi:\R^n\to\R^n$. A typical example is the case in which $\Phi(\nabla u)=|\nabla u|^{p-2}\nabla u$, that is $B=-\kappa|\nabla u|^{p-2}\nabla u$,
meaning that the flux is proportional to a power of the gradient of $u$ (more precisely, the vector $B$ has the same direction of the gradient of $u$ and its magnitude is proportional to a power of the magnitude of the gradient of $u$):
this setting would lead to the so-called “$p$-Laplace equation”.

Another setting of interest is the one in which the flux is related to a pressure drop
(this framework is usually related to the so-called Darcy's LawDarcy's Law), say $B=-\kappa\nabla P$, where $P$ has some physical meaning of pressure. In this sense, the case $P=u$ reduces to (<ref>),
but other situations are of interest. For instance, one could assume that $P$ and $u$ are related by
a state equation of the form $P(x)=p(u(x))$. A natural possibility is to take the function $p$ to be a
power of $u$, e.g. $p=u^m$. This choice would lead to the so-called porous medium equation.

These notes are of elementary nature, hence we will not address the cases of the $p$-Laplace equation,
or of the porous medium equation, or other types of “anomalous diffusions”, such as the type of diffusion that takes into
account mass transfer from remote regions due to long-range interactions.
The reader interested in these more advanced topics may look e.g. to <cit.> and the references therein.

For the readers interested in using Darcy's Law to have a proper cup of coffee, see e.g. <cit.>
and the references therein.] called “Fourier's Law”.

Portrait of Pierre-Simon Laplace by Johann Ernst Heinsius (image from
licensed under the Creative Commons Attribution-Share Alike 4.0 International license).

By substituting (<ref>) into (<ref>) we find that

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj4}
\fint_\Omega\partial_t u(x,t)\,dx=\fint_{\Omega} \div \big(
\kappa(x,t)\,\nabla u(x,t)
\big)\,d{\mathcal{H}}^{n-1}_x+\fint_\Omega f(x,t,u(x,t))\,dx.\end{equation}

Since this identity holds true for all regions $\Omega$, we find a pointwise counterpart of (<ref>) by writing

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj6}
{\partial_t u(x,t)}= \div \big(
\kappa(x,t)\,\nabla u(x,t)
\big)+f(x,t,u(x,t)).\end{equation}

The case of homogeneous media in which $\kappa$ is constant (say, equal to $1$ up to a renormalization of units of measure)
is of particular interest: in this case (<ref>) boils down[We exploit the
standard notation for the LaplacianLaplacian given by

$$ \Delta u:=\div(\nabla u)=\sum_{j=1}^n \frac{\partial^2 u}{\partial x_j} .$$

In some textbook, the use of “$\Delta$” is replaced by “$\nabla^2$”
or by “$|\nabla|^2$”. Of course all notations are good. Personally, we have a preference for
the $\Delta$ notation since we find it simpler to read. Also, in a sense, it highlights the fact that the Laplace operator
has a “dignity” which is “independent from the one of the gradient”.
This “philosophical” point will be perhaps clarified by the geometric interpretation of the Laplacian that we will discuss
in the forthcoming Theorem <ref>.

The name of the Laplace operator comes from Pierre-Simon, marquis de Laplace, who introduced it
while studying celestial mechanics, in connection with the gravitational potential (this approach will
be exploited here when dealing with the fundamental solution in Section <ref>.
See Figure <ref>
for a portrait of Laplace (by Johann Ernst Heinsius) when he was 35 years old, with an easel that quite resembles
a blackboard,  working instrument and sign of distinction of every passionate mathematician.

Let us remark that a slightly different approach to the derivation of the heat equation is possible,
by considering “heat” instead of “temperature” as the main building block of the equation.
The disadvantage of having heat, rather than temperature, in a pivotal role is that it is perhaps a more vague and less intuitive concept
than temperature (just because temperature seems a notion we are so familiar with, e.g. in view of weather forecasts).
However, the notion of heat relates more directly to thermal energy and constitutes an extensive (i.e., additive)
property. In this perspective,
the heat equation reflects an energy budget
in which the heat change in time is equal to heat produced (source) plus heat entering through the boundary (flux).
Namely, if one is willing to take heat, instead of temperature, as the building block for the heat equation,
the balance equation in (<ref>) reads at the level of
variation of heat, instead of averaged temperature
(and note that the constitutive relation in the Fourier's Law (<ref>)
is already a relation between the notions of heat and temperature, since it assumes the
heat flux to be proportional to the gradient of temperature).

The final heat equation in (<ref>) would remain essentially
unchanged since the changes of heat content $Q$ directly relates to the changes in temperature
via the relation $\partial_t Q= c \rho \partial_tu$, where $\rho$ is mass density and $c$ is the specific heat capacity (a constant that depends on the internal properties of the material, which can be tabulated through experiments).]

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj7}
\partial_t u(x,t) = \Delta u(x,t)+f(x,t,u(x,t)).\end{equation}

Equation (<ref>) (or variations of it) is typically[Often,
equations as in (<ref>) are referred to as “parabolic”.
The name comes from the following classification method.
A general linear differential operator of second order in the variables $(X_1,\dots,X_N)$ has the form

$$ \sum _{i,j=1}^{N} a_{i,j}{\frac {\partial ^{2}}{\partial X_{i}\partial X_{j}}}.$$

Then one classifies the operator, and the corresponding partial differential equation,
depending on the sign of the eigenvalues of the coefficient matrix $a_{i,j}$.
More specifically, when all the eigenvalues of $a_{i,j}$ have the same sign (either are all strictly
positive or all strictly negative), the equation is named “elliptic”elliptic.
When one eigenvalue is zero and all the other eigenvalues have the same sign
(either all strictly
positive or all strictly negative), the equation is referred to as “parabolic”parabolic.
When all the eigenvalues have the same sign except one which has the opposite sign
(i.e. if either there is only one negative eigenvalue and all the rest are positive,
or there is only one positive eigenvalue and all the rest are negative),
the equation is named “hyperbolic”hyperbolic.
This classification is not exhaustive: the remaining cases not classified here are typically quite
hard to study, require specific techniques and no general theory is available for those.

In the case of the heat equation in (<ref>), one has $N=n+1$, $X=(x,t)$ and

$$ a_{ij}=\begin{dcases}1 & {\mbox{ if }}i=j\in\{1,\dots,n\},\\
0&{\mbox{ otherwise }}.\end{dcases}$$

The corresponding eigenvalues are therefore $0$, with multiplicity $1$, and $1$, with multiplicity $n$,
hence the equation in (<ref>) is parabolic.]
called “the heat equation”heat equation.

Suppose now that $f$ in (<ref>) is independent of time, i.e. $f=
f(x,u(x,t))$. In this framework,
of particular interest are certainly stationary solutions of (<ref>),
i.e. solutions $u=u(x)$ which are also independent of time: these solutions are actually equilibria of (<ref>), since they provide solutions of (<ref>)
which remain the same at every instant of time. Of course, when $f$ and $u$ are independent of time,
equation (<ref>) boils down to[With respect to the classification
presented in footnote <ref>, equation (<ref>),
(<ref>) and (<ref>) are of elliptic type
(actually, they provide us with the “paradigmatic” form of elliptic equations).
To check their ellipticity, one exploits the setting of footnote <ref>
with $N=n$, $X=x$ and

1 & {\mbox{ if }}i=j,\\
0 &{\mbox{ if }}i\ne j.
\end{dcases}$$

In this scenario, all the eigenvalues of the coefficient matrix
are equal to $1$, hence strictly positive.

Further insight on the notion of ellipticity will be given on pages CLASSIFICATIONFOOTN2
and CLASSIFICATIONFOOTN3.

The above definition of $\delta_{ij}$ will be also used throughout all these notes (sometimes
this definition of $\delta_{ij}$ is referred to with the name of Kronecker notationKronecker notation).]

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj8}
\Delta u(x)+f(x,u(x))=0.\end{equation}

In particular, when $f$ is independent of $u$, equation (<ref>) reduces to

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj9}
\Delta u(x)+f(x)=0,\end{equation}

which is sometimes called Poisson's equation.

The case in which $f$ in (<ref>) vanishes identically produces

\begin{equation}\label{DAGB-ADkrVoiweLL4re2346ytmngrrUj10}
\Delta u(x)=0,\end{equation}

which is called Laplace's equationLaplace's equation.

Solutions of (<ref>) are named “harmonic functions”harmonic function.
The study of harmonic functions is essential, one way or another, in virtually any aspect of mathematics
(including mathematical analysis, mathematical physics, complex analysis, geometry, probability, finance, statistics, you name it). Also, the theory of harmonic functions is a clear example of adamantine beauty and elegance
in which human creativity has reached one of its highest peaks ever.

In the forthcoming pages we will do our best to present some bits of this theory, but we definitely invite
the reader to look to as many other sources as possible to collect the broadest possible amount of information and
develop their own point of view on this topic of paramount importance.

§.§ Population dynamics, chemotaxis and random walks

Another fascinating situation in which partial differential equations surface very often
occurs in the description of biological populations. This is a very captivating field of study,
which is truly crossdisciplinary in spirit and collects fundamental questions from (at least) mathematics, physics,
biology, ethology and social sciences.
For the purpose of these pages, we limit ourselves to a simple description of a
biological species exhibiting at a point $x\in\R^n$ and at a time $t\in\R$ a population density of the form $u(x,t)$.
We assume that a fraction $\rho\in[0,1]$ of
the population has the tendency of moving randomly, while the rest of the population
(corresponding to a fraction $\mu:=1-\rho$ of the totality of individuals) is driven
towards regions of higher concentration of a given attractant (e.g. a chemical signal)
distributed according to a function $w(x,t)$ (the
movement of an organism in response to a chemical stimulus and specifically
in dependence of the increasing or decreasing concentration of a particular substance
is called in jargon “chemotaxis”chemotaxis).
Additionally, we suppose that the population is subject to a driftdrift (e.g. due to the wind, or the stream, or the tide)
in a given direction $b(x,t)$.
The simpler cases in which there is no chemotactic effect (corresponding to $\mu:=0$),
or no random motion (corresponding to $\rho:=0$), or no drift (corresponding to $b:=0$)
are interesting special situations of the complex phenomena for which we now present a mathematical model
(the random component of this discussion will also be further generalized on page 0uojf29249-45kpkfdSmd11493839429efv
to more elaborated environments).

In outline, we consider a small spacial scale $h$ and a small time scale $\tau$.
We will choose in what follows a suitable relation between space and time scales to provide
coherent asymptotics in the formal limit.
In this setting,
at every unit of time the variation of the density of
the population is influenced by the previous factors in a rather explicit way. For the moment, let us focus on the simple case in which there is no chemotaxis and no drift. In this situation, $\rho:=1$, $\mu:=0$, $b$ vanishes identically
and the only drive for the population comes from a random walk.
To get to the bottom of this case,
we consider, at a given point $x$, at time $t+\tau$, the new random population $\rho u(x,t+\tau)$
(which actually coincides with the whole population $u(x,t+\tau)$, but let us keep the fraction $\rho$
for future use, even if for the moment $\rho=1$). Then, we have that $\rho u(x,t+\tau)$
is produced by the random population at the previous time $t$ which is located at some point $x+h\omega$, for a given direction $\omega\in\partial B_1$, times the probability
that these individuals move from $x+he$ to $x$ in the unit of time. If all directions are equally probable, this says that the new random population $\rho u(x,t+\tau)$ is produced by a term of the form $\rho\fint_{\partial B_1} u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega$, that is

\begin{equation}\label{JOSN-Oohpiklfoigheigohbo-gherpoguev}
\rho u(x,t+\tau)
=\rho\fint_{\partial B_1} u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega.
\end{equation}

Subtracting $\rho u(x,t)$ from both sides, dividing by $\tau$ and choosing $h:=\sqrt\tau$,

\begin{equation}\label{7tnfesewuxttau-x21s46}\begin{split}
\rho \frac{u(x,t+\tau)-u(x,t)}\tau\,&
=\,\frac1{\tau}\left[\rho\fint_{\partial B_1} u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega-\rho u(x,t)\right]\\&
=\,\frac{\rho}{h^2}\fint_{\partial B_1} \Big(u(x+h\omega,t)-u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega.\end{split}
\end{equation}

Now we aim at showing that, for small $h$,

\begin{equation}\label{KMS:RUIDKVODIG}
\fint_{\partial B_1} \Big(u(x+h\omega,t)-u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega=
\frac{c h^2}2 \Delta u(x,t) +o(h^2),
\end{equation}

for some $c>0$ depending on the dimension $n$ (see formula (<ref>) below).
For this, we
consider the formal Taylor expansion

$$ u(x+h\omega,t)=u(x,t)+h\nabla u(x,t)\cdot \omega+\frac{h^2}2\,D^2u(x,t) \omega\cdot \omega+o(h^2).$$

We observe, by symmetry, that

$$ \int_{\partial B_1} \nabla u(x,t)\cdot \omega\,d{\mathcal{H}}^{n-1}_\omega=0,$$

because[Spotting these types of cancellations will be often an essential ingredient of our arguments.
After all, one way to obtain significant information about a complex problem is to understand
which quantities do not play much of a role since they “average out”.
At this stage, we discuss these simplifications at an intuitive level. Later on, see e.g. (<ref>),
we will give more rigorous arguments to justify them, but for the moment the main goal is to try to visualize the
possibility of these cancellations and, especially, to appreciate their importance.]
any (positive or negative) contribution to such surface integral coming from a point $\omega\in\partial B_1$
is canceled precisely by the (negative or positive) contribution coming from $-\omega$.
From these observations we arrive at

\begin{eqnarray*} &&\fint_{\partial B_1} \Big(u(x+h\omega,t)-u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega
=\fint_{\partial B_1} \left(h\nabla u(x,t)\cdot \omega+\frac{h^2}2\,D^2u(x,t) \omega\cdot \omega\right)\,d{\mathcal{H}}^{n-1}_\omega+o(h^2)\\&&\qquad=
\frac{h^2}2\,\fint_{\partial B_1}D^2u(x,t) \omega\cdot \omega\,d{\mathcal{H}}^{n-1}_\omega+o(h^2)=\frac{h^2}2
\sum_{i,j=1}^n\fint_{\partial B_1}\partial_{ij} u(x,t) \omega_i\omega_j\,d{\mathcal{H}}^{n-1}_\omega+o(h^2)

We now identify an additional simplification by observing that if $i\ne j$ then

$$ \fint_{\partial B_1}\partial_{ij} u(x,t) \omega_i\omega_j\,d{\mathcal{H}}^{n-1}_\omega=0,$$

because whatever contribution comes from $\omega_i$ is canceled precisely by that coming from $-\omega_i$
(notice how helpful was for all these cancellations that we are integrating over such a symmetric domain as $\partial B_1$,
which remains invariant under all these reflections).

$$ \fint_{\partial B_1} \Big(u(x+h\omega,t)-u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega=\frac{h^2}2
\sum_{i=1}^n\fint_{\partial B_1}\partial_{ii} u(x,t) \omega_i^2\,d{\mathcal{H}}^{n-1}_\omega+o(h^2).

Now, once again, we understand the symmetries of the problem under consideration. For this, we notice that,
for every $i$, $j\in\{1,\dots,n\}$

$$ \fint_{\partial B_1} \omega_i^2\,d{\mathcal{H}}^{n-1}_\omega
=\fint_{\partial B_1} \omega_j^2\,d{\mathcal{H}}^{n-1}_\omega,$$

because[Once again, we aim here at developing arguments in the most intuitive way possible.
This idea with be more rigorously retaken later on, see (<ref>).]
the role played by the $i$th coordinate in the sphere is precisely the same as the one played by the $j$th coordinate.
Therefore, we can define

\begin{equation}\label{CDEGLDDEFC} c:=\fint_{\partial B_1} \omega_i^2\,d{\mathcal{H}}^{n-1}_\omega\end{equation}

and we stress that this quantity does not depend on $i$.

$$ \fint_{\partial B_1} \Big(u(x+h\omega,t)-u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega=\frac{c h^2}2
\sum_{i=1}^n \partial_{ii} u(x,t) +o(h^2),

from which we obtain (<ref>), as desired.

Hence we insert (<ref>) into (<ref>) to find that

\begin{equation}\label{y8oqiwhfesjd8923wfyug3iuweLuhonwfdpinrgNqikedOKS}\begin{split}&
\rho \frac{u(x,t+\tau)-u(x,t)}\tau
=\frac{c\rho}{2}\Delta u(x,t) +o(1).
\end{split}
\end{equation}

By sending $\tau\searrow0$ (and thus $h\searrow0$), we thereby find that

\begin{equation}\label{HJA:A89MAOL8789892}
\rho \partial_t u(x,t)
= \frac{c\rho}{2}\Delta u(x,t) ,
\end{equation}

which coincides, up to constant, with the heat equation[See
e.g. <cit.>
for more information about the strong connections between
the random walk and the heat equation.] presented in (<ref>)
(here, with no source term).

Hard-working scientists (Public Domain image from

This computation was not only instructive from the technical point of view, but it also unveiled one of the mainsprings
of the mathematical theory of diffusion, revealing the strong conceptual connection
of a substance (or a population) moving randomly according to a Brownian motionBrownian motion and the way in which
heat dissipates: after all, the spreading of temperature all over the region subject to diffusion is, in some sense, nothing else[The link between the random motion of molecules and the
macroscopic phenomenon of diffusion is indeed a deep feature. While the Brownian motion was discovered
in 1827 by botanist Robert Brown by looking at dust grains floating in water, Albert Einstein
in his “annus mirabilis” 1905
gave a quantitative model for the motion of floating particles as being moved by individual water molecules <cit.>. Not only this article
founded the statistical physics analysis of Brownian motion, but also it
provides a way to determine the mass and the dimensions of the atoms involved in the
process, thus changing atomic theory from a controversial set of conjectures
into an established fact of science. We will come back to the theory proposed by Einstein in Section <ref>.

See Figure <ref>
depicting Albert Einstein together with Niels Bohr.

By the way, Niels Bohr's brother was Harald Bohr, who was a mathematician, pioneering almost periodic functions, and a soccer player. The Denmark national soccer team in which he was playing took part in the  1908 Summer Olympics, where football was an official event for the first time. Harald Bohr scored two goals in the first match and played in the semifinal (Denmark 17 - France 1, which remains the Olympic record of the most goals scored by one team). After that, the Denmark team lost the final and won the silver medal. In all the subsequent editions of the Olympic Games, Denmark soccer team won two more silver medals (in 1912 and 1960) and a bronze (in 1948); hence Harald Bohr's team also holds the record (ex aequo) of best placement for a Denmark soccer team in the Olympics.

See Figure <ref>
for a photo of the Danish soccer team in the 1908 Olympic Games: Harald Bohr is the second player from the left in the top row.

Niels Bohr was also a passionate soccer player. He played goalkeeper in the Copenhagen-based team Akademisk Boldklub (at the time, one of the best clubs in Denmark; actually, the two brothers played several matches together in this team). One can wonder however why Niels never made it to the national team. Well, according to https://www.theguardian.com/football/2005/jul/27/theknowledge.panathinaikos
in a match against a German team, one of the midfielders of the opposite team launched a very long shot and Niels, leaning against the post, did not react at all, missing an easy save and letting the German team score. After the game, Niels admitted that on that occasion he had been distracted by a mathematical problem he was thinking about. Strangely enough, Niels did not play for Akademisk Boldklub after that season.]
but the tendency of random movements to distribute mass in average from regions of high density
towards regions of small density. In both cases, the process
has the strong tendency to “balance out” differences:
as time goes, hot spots lose their temperature in favor of cold spots which are heated up by
the temperature of their neighbors, as well as under a random walk
the regions with low density get occupied by the population
coming from the highly populated regions.

This “democratic” tendency of averaging out differences is typical of the Laplace operator and it
will be the underlying feature of all the “regularity theories” for elliptic equations that we will present in the forthcoming

Danish soccer team at the 1908 Olympic Games (Public Domain image from

Having well understood the case in which the biological species is subject solely on a random motion,
we now retake the interesting case in which a chemotactic agent
also comes into play, for a fraction $\mu$ of the population (for the moment, we still suppose that there is no additional drift, hence $b$ vanishes identically). In this situation,
the fraction $\rho$ of the population performing the random walk would be still described
by (<ref>) (that was the reason to include $\rho$ in the previous
computation, even if before $\rho$ was just equal to $1$) and we focus now on
the fraction $\mu$ of the population following the attractant $w$.
To appreciate the effect of the chemotactic factor, we suppose that such a fraction $\mu$ of the population with density $u$ does not move completely randomly, but, at each time step, picks a direction of motion with a probability $\lambda$ that is proportional to values of the attractant having density $w$.

More explicitly, given $\omega\in\partial B_1$,
we suppose for the chemotactic population that a jump from the point $x$ to the point $x+h\omega$
occurs with probability

\begin{equation}\label{RISUGIBNSUPJDMD} \lambda(x,\omega,t):=\frac1{{\mathcal{H}}^{n-1}(\partial B_1)}+w(x+h\omega,t)-w(x-h\omega,t).\end{equation}

In this setting, we are assuming that the oscillations of $w$ are sufficiently small[More precisely,
one could rewrite (<ref>) as

\[ \lambda(x,\omega,t):=\frac1{{\mathcal{H}}^{n-1}(\partial B_1)}+\e_0\Big(w(x+h\omega,t)-w(x-h\omega,t)\Big),\]

where $\e_0>0$ is a parameter that takes into account the “sensibility” of the chemotactic population
to the chemical attractant. Just to be consistent with the probability scenario, it is convenient
to think that $\lambda\ge0$ in view of the smallness of $\e_0$ with respect to the oscillations
of the attractor's density $w$.]
to make the above quantity positive.
In this sense, the probability $\lambda$ differs from that of the classical random walk (corresponding to $w$
being constant) since it increases in favor of the directions $\omega$ in which the concentration of the chemical
attractant is higher. We stress that the above setting provides a normalized probability,
since the total probability of
jumping to the sphere of radius $h$ around the given point $x$ is $1$, because

\begin{equation*}
\int_{\partial B_1}\lambda(x,\omega,t)\,d{\mathcal{H}}_\omega^{n-1}
=1+\int_{\partial B_1}\Big(w(x+h\omega,t)-w(x-h\omega,t)\Big)\,d{\mathcal{H}}_\omega^{n-1}=1,
\end{equation*}

due to odd symmetry cancellations.

We also point out that, in the above setting, the probability of a jump from a point $a$ to
a point $b$ with $|a-b|=h$ is given by

\begin{eqnarray*} \lambda\left(a,\frac{b-a}{h},t\right)&=&\frac1{{\mathcal{H}}^{n-1}(\partial B_1)}+w\left(a+h\frac{b-a}{h},t\right)-w\left(a-h\frac{b-a}{h},t\right)\\&=&
\frac1{{\mathcal{H}}^{n-1}(\partial B_1)}+w(b,t)-w(2a-b,t).\end{eqnarray*}

In consequence, given $\omega\in\partial B_1$,
taking $a:=x+h\omega$ and $b:=x$, we find that
the jump from $x+h\omega$ to $x$ occurs with probability

$$ \lambda (x+h\omega,-\omega,t)=
\frac1{{\mathcal{H}}^{n-1}(\partial B_1)}+w(x,t)-w(x+2h\omega,t).$$

This being so, we can detect the chemotactic counterpart of (<ref>).
Indeed, since the chemotactic population $\mu u(x,t+\tau)$ is produced by a term of the form $\mu\int_{\partial B_1} \lambda (x+h\omega,-\omega,t) u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega$, we have that

\begin{equation*} \begin{split}
\mu u(x,t+\tau)&=\mu\int_{\partial B_1} \lambda (x+h\omega,-\omega,t) u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega\\
&=\mu\int_{\partial B_1} \left(\frac1{{\mathcal{H}}^{n-1}(\partial B_1)}+w(x,t)-w(x+2h\omega,t)\right)u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega\\&=\mu\fint_{\partial B_1} u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega
+\mu\int_{\partial B_1} \Big(w(x,t)-w(x+2h\omega,t)\Big)u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega

Thus, we can subtract to both sides $\mu u(x,t)$, divide by $\tau=h^2$ and recall (<ref>)
and (<ref>) to find that

\begin{eqnarray*}
&&\mu \frac{u(x,t+\tau)-u(x,t)}\tau\\ &=&
\frac1{h^2}\left[\mu\fint_{\partial B_1}\Big( u(x+h\omega,t)-u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega+
\mu\int_{\partial B_1} \Big(w(x,t)-w(x+2h\omega,t)\Big)u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega
\right]\\&=& \frac{c\mu}{2}\Delta u(x,t) +
\frac{\mu}{h^2}\int_{\partial B_1} \Big(w(x,t)-w(x+2h\omega,t)\Big)u(x+h\omega,t)\,d{\mathcal{H}}^{n-1}_\omega
+o(1)\\&=& \frac{c\mu}{2}\Delta u(x,t)-
\frac{\mu}{h^2}\int_{\partial B_1} \Big(
2h\nabla w(x,t)\cdot\omega
\Big)\Big(u(x,t)+h\nabla u(x,t)\cdot\omega\Big)\,d{\mathcal{H}}^{n-1}_\omega\\&&\qquad\qquad\qquad
+o(1)\\&=& \frac{c\mu}{2}\Delta u(x,t) -
\frac{\mu}{h^2}\int_{\partial B_1} \Big(
2h^2\nabla w(x,t)\cdot\omega\,\nabla u(x,t)\cdot\omega
+2h^2D^2w(x,t)\omega\cdot\omega\, u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega
+o(1)\\&=& \frac{c\mu}{2}\Delta u(x,t) -
2\mu \int_{\partial B_1} \Big(
\nabla w(x,t)\cdot\omega\,\nabla u(x,t)\cdot\omega
+D^2w(x,t)\omega\cdot\omega\, u(x,t)\Big)\,d{\mathcal{H}}^{n-1}_\omega\\&&\qquad\qquad\qquad
+o(1)\\&=& \frac{c\mu}{2}\Delta u(x,t) -
2\mu \sum_{i=1}^n\int_{\partial B_1} \Big(\partial_i w(x,t)\partial_i u(x,t)\omega_i^2
+\partial_{ii}w(x,t) u(x,t)\omega_i^2\Big)\,d{\mathcal{H}}^{n-1}_\omega
\\&=& \frac{c\mu}{2}\Delta u(x,t) -
2\mu \sum_{i=1}^n\int_{\partial B_1} \partial_i\big(\partial_i w(x,t) u(x,t)\big)\omega_i^2
\,d{\mathcal{H}}^{n-1}_\omega
+o(1)\\&=& \frac{c\mu}{2}\Delta u(x,t) -
\widetilde{c}\mu \sum_{i=1}^n \partial_i\big(\partial_i w(x,t) u(x,t)\big)
+o(1)\\&=& \frac{c\mu}{2}\Delta u(x,t) -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)
\end{eqnarray*}

for a suitable $
\widetilde{c}>0$.

Notice that here above we have repeatedly taken advantage of odd symmetry cancellations.
Gathering this and (<ref>) we thereby find that

\begin{equation}\label{DRIFPAG2}\begin{split}&
\frac{u(x,t+\tau)-u(x,t)}\tau=
\rho \frac{u(x,t+\tau)-u(x,t)}\tau +\mu \frac{u(x,t+\tau)-u(x,t)}\tau\\&\qquad
=\frac{c\rho}{2}\Delta u(x,t) +\frac{c\mu}{2}\Delta u(x,t) -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)+o(1)\\&\qquad
=\frac{c}{2}\Delta u(x,t) -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)+o(1).\end{split}\end{equation}

Sending $\tau\searrow0$ (and thus $h\searrow0$) we obtain

\begin{equation}\label{CHEMLKSDEQ}
\partial_t u(x,t)=\frac{c}{2}\Delta u(x,t)  -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big).\end{equation}

For constant distributions $w$ of the chemical attractant, equation (<ref>)
can be seen as a special case (up to normalizing constants) of the heat equation in (<ref>).

Stationary solutions of (<ref>) solve

\begin{equation}\label{OJS-PJDN-0IHGDOIUGDBV02ujrf}
\frac{c}{2}\Delta u(x,t)  -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)=0.\end{equation}

This is the equation solved by the equilibria corresponding to a biological population in the presence of
a chemotactic factor.

Let us now introduce one last complication in the model by accounting for a possible drift
(as mentioned on page DRIFPAG). In this situation, we may suppose that, in the time unit,
the population is also moved by the velocity vector $b$ by a space displacement $b\tau$.
Regarding this, the population density $u(x+2\tau,t)$

has in fact moved to the location $x+b\tau$: correspondingly, to take into account this drift,
one can observe that $ u(x,t+2\tau)=u(x-b(x,t)\tau,t+\tau)$, and accordingly, recalling (<ref>)
and using that

$$u\big(x-b(x,t)\tau,t+\tau\big)=u(x,t )
-\tau b(x,t)\cdot\nabla u (x,t )+\tau\partial_t u(x,t)
we have that

\begin{eqnarray*}
&& \frac{u(x,t+2\tau)-u(x,t)}\tau= \frac{u(x,t+2\tau)-u(x,t+\tau)}\tau+ \frac{u(x,t+\tau)-u(x,t)}\tau
\\&&\qquad
\frac{c}{2}\Delta u(x,t) -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)+o(1)\\&&\qquad=
\frac{u(x,t )
-\tau b(x,t)\cdot\nabla u (x,t )+\tau\partial_t u(x,t)-
u(x,t)-\tau\partial_t u(x,t) +o(\tau)
\\&&\qquad\qquad+
\frac{c}{2}\Delta u(x,t) -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)+o(1)\\&&\qquad=
-b(x,t)\cdot\nabla u (x,t )+
\frac{c}{2}\Delta u(x,t) -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)+o(1)
\end{eqnarray*}

Thus, we can pass to the limit and obtain

\begin{equation}\label{OJS-PJDN-0IHGDOIUGDBV02ujrf220}
\partial_t u(x,t)=\frac{c}{2}\Delta u(x,t) -
\widetilde{c}\mu \div\big(u(x,t) \nabla w(x,t)\big)-b(x,t)\cdot\nabla u(x,t).\end{equation}

Stationary solutions of (<ref>) give rise to the equation

\begin{equation}\label{OJS-PJDN-0IHGDOIUGDBV02ujrf22}
\Delta u(x,t) -
\div\big(u(x,t) \nabla w(x,t)\big)-b(x,t)\cdot\nabla u(x,t)=0,\end{equation}

up to normalizing constants that we have omitted for the sake of simplicity.
In this spirit, (<ref>) models[The equation
in (<ref>), together with the diffusion of the attractant,
reveals the complexity
of the model under consideration which, depending on different parameter thresholds, exhibit a variety of
different stable geometrical patterns, colony formation, blowup mechanisms, flocculation and aggregation phenomena.

All these traits play a fundamental role for life as we know it, since these configurations can arise
from equations such as the ones presented here and possess quite a precise counterpart in nature, when
the corresponding behaviors occur, for example, in response to external stress,
or to lack of resources (sometimes the chemical attractant is produced by the biological species itself,
e.g. by amoebae in case of food scarceness) or in response to predation
(autoaggregation provides protection and defense against predators). These aggregation processes often produce slimy yet durable coatings, called biofilms. Even blowup phenomena
can be exploited by bacteria to enhance their possibility of surviving in hostile environments
(e.g., the individual on the top of the tower formed by the blowup have more chances to be picked up, say, by the wind
or by an external factor, and be deposited possibly in a less hostile environment), so all in all the variety
of patterns exhibited by the solutions of a relatively simple differential equation is in some sense the mathematical counterpart
of the variety of ways in which “life finds a way” (as Ian Malcolm utters in “Jurassic Park”).] the equilibria reached by a biological species possibly subject
to chemotaxis and to an external drift.

Notice that (<ref>) reduces to (<ref>)
in the absence of drift.

We have written (<ref>) in a form which emphasizes its “divergence structure”:
if instead one wants to highlight the presence of a second order
differential operator, it suffices to notice that $\div(u \nabla w)=u\Delta w+\nabla w\cdot\nabla u$ to recast (<ref>) into

\begin{equation}\label{OJS-PJDN-0IHGDOIUGDBV02ujrf22NDF}
\Delta u(x,t) -\big( b(x,t)+\nabla w(x,t)\big)\cdot\nabla u(x,t)- u(x,t)\Delta w(x,t)=0,\end{equation}

which is also a telling expression since it is revealing that the effect of the chemotaxis
is to alter the diffusion with an additional drift term $\nabla w$
which is favoring the movements towards regions with higher concentrations
of the the attractant.

We refer to <cit.> and the references therein for further information about chemotaxis and its biological implications.

§.§ Pattern formation, or how the leopard gets its spots

We know well that the natural world offers a variety of amazing
regular patterns, such as symmetries, tessellations, spirals, spots, etc.,
see e.g. Figure <ref>.
These patterns visible in nature
are determinant factors in the processes of natural and sexual selection,
providing organisms with structures recognizable by their conspecifics (thus favoring social
interactions and reproduction) and often realizing optimal configurations for mobility, hunt or camouflage (see e.g. <cit.> and the references therein).

Leopard mating dance (image from
Wikipedia, licensed under the Creative Commons Attribution 2.0 Generic license;
see Steve Jurvetson's website https://www.flickr.com/photos/jurvetson/5913330010/
for the complete series of shots).

The biological process leading to the development of the specific shape
of an organism is dubbed “morphogenesis” morphogenesis
and a scientific investigation of the formation of patterns in nature
is relatively recent.
One of the pioneer scientists interested in the mathematical analysis of the growth, form and evolution of plants and animals
was Sir D'Arcy Wentworth Thompson: his 1116 pages book <cit.>
(793 pages in the first edition of 1917) combined classical natural philosophy,
biology and mathematics to give insight on a number of biological shapes
and analyzed the differences in the forms in nature in the light of mathematical transformations.

A modern mathematical treatment of morphogenesis was initiated
by[Not only the father of mathematical morphogenesis,
computer science and artificial intelligence, Alan Turing was also an exceptional long-distance runner,
capable of world-class marathon standards. See https://kottke.org/18/04/alan-turing-was-an-excellent-runner
for a nice photo of Turing during one of his running performances.

During World War II, Turing worked as a code breaker for the British Intelligence,
devising an electromechanical machine, named “Bombe”,
to decipher the encrypted messages created by the German cipher device “Enigma”.
The headquarter of the British code breakers was located in a
mansion named Bletchley Park which hosts today a statue of Turing
made with slates (by Stephen Kettle). In this statue, Turing is depicted seated and looking at a German Enigma machine,
see Figure <ref>.

Turing was 39 years old when he was charged with “gross indecency” on the basis of a homosexual relationship
under Section 11 of the Criminal Law Amendment Act 1885.
The case, Regina v Turing and Murray, was brought to trial on 31 March 1952
and Turing pleaded guilty, insisting that he saw nothing wrong with his actions.

He was convicted and given a choice between imprisonment and chemical castration via synthetic estrogens.
As a result of choosing the second option, Turing suffered of severe mutations of his body,
was barred from his work for Government Communication Headquarters and
was denied entry into the United States.
On June 8 1954, Turing was found dead of cyanide poisoning with a half-eaten apple
next to his bed lay. The inquest into his death recorded a verdict of suicide, but
the apple was not tested for cyanide.

Turing's favorite fairy tale was Snow White and the Seven Dwarfs.]
Alan Mathison Turing in 1952, see <cit.>.
Turing's brilliant idea is that simple physical laws were sufficient to justify
the shaping of complex patterns, such as
animal markings and arrangement of leaves and florets in plants.
In a nutshell, Turing proposed that the process of morphogenesis
is regulated by the interaction of chemical substances, called morphogens, which
diffuse through a tissue: these substances could be
hormones, or genes, or any other essence which can act and react in presence of another one.
In practice, see e.g. <cit.>, one of the morphogens may act as an
“activator”, which is self-sustaining and introduces positive feedback, while the other
may play the role of an “inhibitor” which tends to suppress the self-amplification of the activator.
In this interplay, the pattern may be created by the different speeds of diffusion
of the two substances:
namely, the faster diffusion of the inhibitor
can catch up with the activator's  self-replication (that is, roughly speaking, on the one side the activator's capacity of self-replicating could be strong enough to produce
local patches, but the predominant speed of the inhibitor could avoid that these patches grow
incessantly). All in all, the whole process can be thereby considered as a “diffusion driven instability”. See <cit.> for a very clear explanation of the roles of activators and inhibitors, via an example with sweating grasshoppers.

Statue of Alan Turing (photo by Jon Callas; image from
Wikipedia, licensed under the Creative Commons Attribution 2.0 Generic license.).

The mathematical formulation of Turing's idea combines the notion of diffusion,
as modeled for instance by the heat equation (<ref>),
and that of reaction, taking into account that each morphogen can chemically react with the others
and the effect of this interaction can depend on the concentration of the diffusing substances in the tissue. The combination of reaction and diffusion in this type
of systems of partial differential equations justified the name
of “reaction-diffusion
equations”. reaction-diffusion equation

For instance, one can consider the case of two morphogens with density $u$ and $v$ respectively
and a system of reaction-diffusion equations of the form

\begin{equation}\label{SYS:READI}
\begin{dcases}
&\partial_t u=\mu\Delta u+f(u,v),\\
&\partial_t v=\nu\Delta v+g(u,v),
\end{dcases}
\end{equation}

where $\mu$ and $\nu$ are positive coefficients which model the speed of diffusion
of the chemical substances with density $u$ and $v$ respectively.

Now, the formation of patterns as an outcome of (<ref>)
is, in a sense, not completely intuitive: on the contrary, given the
“democratic” tendency of the Laplace operator (as discussed on page DEEFFNVST)
one may imagine that the solution $(u(t),v(t))$, as $t\to+\infty$, will evolve spontaneously towards
some constant values $(u_0,v_0)$, which are just the common zeros of $f$ and $g$
(that is, such that $f(u_0,v_0)=g(u_0,v_0)=0$). And this is indeed one of the possible
destination for the solutions of (<ref>). However, there is also a more intriguing possibility:
the constant state $(u_0,v_0)$ may well exist (and it may also be “stable” for small
perturbations when $\mu:=0$ and $\nu:=0$), but it may be triggered off by random disturbances.
In this situation, when $\mu>0$ and $\nu>0$, the solutions may end up drifting away
from the constant $(u_0,v_0)$. Thus (unless for some reasons the solution diverges)
several interesting patterns may arise, such as oscillations between equilibria,
stationary waves, moving fronts, etc.

Though a full understanding of Turing's theory of morphogenesis goes well beyond the scopes of this
set of notes, following <cit.> one can at least grasp some of the ideas involved. For instance,
one can consider a simple subcase of (<ref>) in which $f$ and $g$ are linear functions
(this simplification can also be inspiring to treat the general case, since, in the vicinity of
the constant equilibrium $(u_0,v_0)$,
one can try to “linearize the equation” to obtain information on its dynamics).
That is, let us consider the system of equations

\begin{equation}\label{52568-345256890hro3}
\begin{dcases}
&\partial_t u=\mu\Delta u+a(u-u_0)+b(v-v_0),\\
&\partial_t v=\nu\Delta v+c(u-u_0)+d(v-v_0),
\end{dcases}
\end{equation}

for some $a$, $b$, $c$, $d\in\R$. For simplicity, let us also suppose that the problem is set on a circle,
namely $x\in\R$ and $u$ and $v$ are periodic functions of period $2\pi$.
One can thus look for solutions of (<ref>) in Fourier Series of the form

\begin{equation}\label{52568-345256890hro4} u(x,t)=u_0+\sum_{j\in\Z} U_j(t)\,e^{ijx}\qquad{\mbox{and}}\qquad
v(x,t)=v_0+\sum_{j\in\Z} V_j(t)\,e^{ijx},\end{equation}

with $U_j$ and $V_j$ to be determined.

Substituting (<ref>) into (<ref>) we have that

\begin{eqnarray*}&&
\sum_{j\in\Z} \dot U_j(t)\,e^{ijx}=
\partial_t u=\mu\Delta u+a(u-u_0)+b(v-v_0)\\&&\qquad
= -\mu\sum_{j\in\Z} j^2U_j(t)\,e^{ijx}
+a\sum_{j\in\Z} U_j(t)\,e^{ijx}+b\sum_{j\in\Z} V_j(t)\,e^{ijx}
\end{eqnarray*}

and similarly

\[ \sum_{j\in\Z} \dot V_j(t)\,e^{ijx}=
-\nu\sum_{j\in\Z} j^2V_j(t)\,e^{ijx}
+c\sum_{j\in\Z} U_j(t)\,e^{ijx}+d\sum_{j\in\Z} V_j(t)\,e^{ijx}.\]

From these equations we arrive at

\begin{equation}\label{52568-345256890hro4-1} \left(\begin{matrix} \dot{U}_j(t)\\ \dot{V}_j(t)
\end{matrix}\right)=
\left( \begin{matrix}a-j^2\mu & b\\
c& d-j^2\nu\end{matrix}\right)
\left(\begin{matrix} {U}_j(t)\\ {V}_j(t)
\end{matrix}\right).
\end{equation}

This is a first order ordinary differential equation with constant coefficients.
Hence, we suppose for simplicity that the (possibly complex) eigenvalues of the matrix $\left( \begin{matrix}a-j^2\mu & b\\
c& d-j^2\nu\end{matrix}\right)$, which we denote by $\lambda_j$ and $\Lambda_j$, are distinct.
We also denote by $w_j$ and $W_j\in \C^2$ the corresponding eigenvectors.
With this notation,
we find (see e.g. <cit.>) that the solutions of (<ref>)
are of the form

\begin{equation}\label{52568-345256890hro4-2} \left(\begin{matrix} {U}_j(t)\\ {V}_j(t)\end{matrix}\right)
=\xi_j \,e^{\lambda_j t} w_j+\Xi_j\, e^{\Lambda_j t} W_j
\end{equation}

for some $\xi_j$, $\Xi_j\in\C$.

Strem plot of the system of ordinary differential equations in (<ref>).

It is convenient to use the vector notation

$$ \xi_j=\left(\begin{matrix} \xi_{j1}\\ \xi_{j2}\end{matrix}\right)\qquad{\mbox{and}}\qquad
\Xi_j=\left(\begin{matrix} \Xi_{j1}\\ \Xi_{j2}\end{matrix}\right)$$

and let

$$ A_j:=\xi_j w_{j1},\qquad B_j:=\Xi_j W_{j1},\qquad
C_j:=\xi_j w_{j2}\qquad{\mbox{and}}\qquad D_j:=\Xi_j W_{j2}.$$

In this way, (<ref>) yields that

\begin{equation}\label{UJ:lao-myND9ona9nyuton-0201}
U_j(t)=A_j e^{\lambda_j t}+B_j e^{\Lambda_j t}\qquad{\mbox{and}\qquad}
V_j(t)=C_j e^{\lambda_j t}+D_j e^{\Lambda_j t}.\end{equation}

\begin{equation*}\begin{split}&\left( \begin{matrix}a-j^2\mu & b\\c& d-j^2\nu\end{matrix}\right)
\left( \begin{matrix}A_j\\ C_j\end{matrix}\right)=
\left( \begin{matrix}a-j^2\mu & b\\c& d-j^2\nu\end{matrix}\right)
\left( \begin{matrix}\xi_jw_{j1} \\ \xi_jw_{j2}\end{matrix}\right)
\\&\qquad=\xi_j\left( \begin{matrix}a-j^2\mu & b\\c& d-j^2\nu\end{matrix}\right)w_j=
\xi_j\lambda_jw_j=\lambda_j\left( \begin{matrix}A_j\\C_j\end{matrix}\right),
\end{split}\end{equation*}

leading to

\begin{equation}\label{UJ:lao-myND9ona9nyuton-0202}
(a-j^2\mu-\lambda_j)A_j+b C_j=0.
\end{equation}

\begin{equation}\label{UJ:lao-myND9ona9nyuton-0203}
(a-j^2\nu-\Lambda_j)B_j+b D_j=0.\end{equation}

In jargon, (<ref>) and (<ref>)
are sometimes dubbed “dispersion relations”dispersion relation:
their interest lies in the fact that they relate the speed of oscillation in the time variable
(quantified in (<ref>)
by the eigenvalues $\lambda_j$ and $\Lambda_j$) with the spatial periodicity of the medium
(characterized by the eigenvalues $-j^2$ of the one-dimensional Laplacian and modulated
by the speeds of diffusion $\mu$ and $\nu$).

To recap briefly, from the system of reaction-diffusion equations in (<ref>)
one arrives at the solutions introduced in (<ref>),
with $U_j$ and $V_j$ as in (<ref>),

\begin{equation}\label{o2re2s3p4on21di21ng214ly}
{\mbox{$\lambda_j$ and~$\Lambda_j$ are (distinct)
complex eigenvalues of the matrix }}\left( \begin{matrix}a-j^2\mu & b\\
c& d-j^2\nu\end{matrix}\right),\end{equation}

and with the parameters satisfying (<ref>)
and (<ref>).

As detailed in <cit.>
and <cit.>, this explicit mathematical construction
has a number of important biological consequence and presents a sufficiently rich structure
to account for many patterns visible in nature. To see these features, one may focus
on the case in which one of the eigenvalues has the largest real part (roughly speaking,
one expects that the other modes are dominated by this one). Also,
it is convenient to distinguish between the case in which the dominant eigenvalue is real
from the one in which it is complex and with nonzero imaginary part: indeed, real eigenvalues will be related to stationary states
and complex eigenvalues to oscillatory cases.

More specifically, suppose that

\begin{equation}\label{CMOPSDEArajeke24w}
{\mbox{$\Lambda_{j_0}$ is the eigenvalue with largest real part.}}\end{equation}

We notice that also $\Lambda_{-j_0}=\Lambda_{j_0}$, since

\begin{equation}\label{INVAMAJJ}
{\mbox{the matrix }}\left( \begin{matrix}a-j^2\mu & b\\ c& d-j^2\nu\end{matrix}\right) {\mbox{ remains the same if we exchange~$j$ with~$-j$}}.\end{equation}

Hence, dropping the higher order terms, we can assume that the dynamics of the solutions in (<ref>) is governed by the following long-time asymptotics:

\begin{equation}\label{INVAMAJJ2}
\begin{split}&
u(x,t)\simeq u_0+ U_{j_0}(t)\,e^{i{j_0}x}+U_{-j_0}(t)\,e^{-i{j_0}x}
\simeq u_0+ e^{\Lambda_{j_0} t}\Big( B_{j_0} e^{i{j_0}x}+B_{-j_0} e^{-i{j_0}x}\Big)\\
v(x,t)\simeq v_0+ e^{\Lambda_{j_0} t}\Big( D_{j_0} e^{i{j_0}x}+D_{-j_0} e^{-i{j_0}x}\Big).\end{split}
\end{equation}

The invariance in (<ref>) also suggests that if $A_j$, $B_j$, $C_j$ and $D_j$ are solutions of (<ref>) and (<ref>), then so are $A_{-j}$, $B_{-j}$, $C_{-j}$ and $D_{-j}$: for this reason, we can suppose that $B_{-j_0}=B_{j_0}$
and $D_{-j_0}=D_{j_0}$
in (<ref>), obtaining that

\begin{equation}\label{INVAMAJJ3}
\begin{split}&
u(x,t)\simeq u_0+ B_{j_0} e^{\Lambda_{j_0} t}\Big( e^{i{j_0}x}+e^{-i{j_0}x}\Big)=
u_0+ 2B_{j_0} e^{\Lambda_{j_0} t}\cos( j_0x)\\
v(x,t)\simeq v_0+2D_{j_0} e^{\Lambda_{j_0} t}\cos( j_0x).\end{split}
\end{equation}

Without aiming at exhausting all the possible patterns included in (<ref>),
let us now show a concrete case of interest.

For instance, let $\vartheta\in\N$,
\begin{equation}\label{G1A2MbdRbnU2U8n4ngVA9}\gamma:=
\frac{33\,\vartheta^2}{13\sqrt{3}-9}
\end{equation}
and also

\begin{equation}\label{MSUOJS5UIHJ6OR6YUFH8YTFUYHGJFDTYG}
u_0:=0,\qquad v_0:=0, \qquad a:=\gamma, \qquad b:=-2\gamma,\qquad c:=2\gamma\qquad {\mbox{and}}\qquad d:=-2\gamma. \end{equation}

In this case,
the system in (<ref>) describes an activator with density $u$
which is an autocatalytic activator: that is, such a substance stimulates the production
of itself (since $a>0$) and also activates the production of the substance with density $v$ (since also $c>0$). Also, the substance with density $v$ corresponds to a self-degrading inhibitor:
indeed, higher concentrations of this reactant are noxious for itself (since $d<0$)
and for the activator with density $u$ (since $b<0$).

Interestingly, the origin, which corresponds to the equilibrium $(u_0,v_0)$, is a stable sink for the system of ordinary differential equations corresponding to (<ref>) when $\mu:=0$ and $\nu:=0$. See Figure <ref> for a sketch of the trajectories

\begin{equation}\label{STREAMPLO}
\begin{dcases}
&\partial_t u= \gamma u -2\gamma v,\\
&\partial_t v=2\gamma u-2\gamma v.
\end{dcases}
\end{equation}

Quite remarkably, as discovered by Turing, the stability of (<ref>) can be destroyed by random fluctuations arising from the diffusivity of the chemical reactant. To appreciate this, given $u_0$, $v_0$, $a$, $b$, $c$ and $d$ as in (<ref>),
we take $\mu:=1$ and $\nu:=12$ in (<ref>). Note that this corresponds to a situation in which the diffusion of the inhibitor is faster than the one of the activator.
This scenario gives that

$$ \left( \begin{matrix}a-j^2\mu & b\\ c& d-j^2\nu\end{matrix}\right)=
\left( \begin{matrix}\gamma-j^2 & -2\gamma\\ 2\gamma& -2\gamma- 12j^2\end{matrix}\right),$$

which possesses eigenvalues of the form

\begin{equation}\label{CMOPSDEArajeke24w2}-\frac12\left(\gamma+13 j^2  \pm \sqrt{(11j^2+7\gamma)(11j^2-\gamma)}\right).\end{equation}

By (<ref>), we aim at detecting the greatest possible real part in (<ref>). To this end, note that when $(11j^2+7\gamma)(11j^2-\gamma)\le0$
then the real part in (<ref>) is equal to $-\frac12(\gamma+13 j^2)<0$.
Instead, if $(11j^2+7\gamma)(11j^2-\gamma)>0$ the largest possible real part in (<ref>) is equal to

\begin{equation}\label{LESOCHEDIBMECOHowrj-0}\begin{split}&
\sup_{j\in\Z}\frac12\left(\sqrt{(11j^2+7\gamma)(11j^2-\gamma)}
-\gamma-13 j^2\right)\\=\,&\frac\gamma2\sup_{j\in\Z}\left(
\sqrt{\left(\frac{11j^2}\gamma+7\right)\left(\frac{11j^2}\gamma-1\right)
}-1-\frac{13 j^2}\gamma\right)\\=\,&
\frac\gamma2\sup_{j\in\Z}\Phi\left(\frac{j^2}\gamma\right),\end{split}\end{equation}

$$ \Phi(\tau):=\sqrt{(11\tau+7)(11\tau-1)}-1-13 \tau.$$

Using elementary calculus, one checks that

$$ \max_{\tau\ge0}\Phi(\tau)= \frac4{11} \left(7 -4\sqrt{3}\right)

thanks to (<ref>).

This observation and (<ref>) give that the eigenvalues with largest possible real part in (<ref>) correspond to the choice $j:=\vartheta$ and are of the form

$$ \frac{2\gamma}{11} \left(7 -4\sqrt{3}\right)
=\frac{66\left(7 -4\sqrt{3}\right)\,\vartheta^2}{11\left(13\sqrt{3}-9\right)}

Hence, up to constants, the corresponding setting in (<ref>) takes the form

\begin{equation}\label{LESOCHEDIBMECOHowrj}
\begin{split}&
u(x,t)\simeq  \exp\left(\frac{66\left(7 -4\sqrt{3}\right)\,\vartheta^2}{11\left(13\sqrt{3}-9\right)}\,t\right)\cos(\vartheta x)\\
v(x,t)\simeq \exp\left(\frac{66\left(7 -4\sqrt{3}\right)\,\vartheta^2}{11\left(13\sqrt{3}-9\right)}\,t\right)\cos (\vartheta x) .\end{split}
\end{equation}

See Figure <ref>
for an example with $\vartheta:=2$.

Plot of the function $(x,t)\mapsto\exp\left(\frac{264\left(7 -4\sqrt{3}\right)}{11\left(13\sqrt{3}-9\right)}\,t\right)\cos (2 x)$.

Of course, the asymptotics in (<ref>) are divergent as $t\to+\infty$,
which would correspond to the chemical substances to reach infinite density, which is certainly unfeasible in practice,
hence the meaning of (<ref>) has to be understood only in the vicinity of the equilibrium $(u_0,v_0)$,
which was set to be the origin for simplicity. Indeed, in practice the linear system in (<ref>)
must be considered as an efficient linearization only in the vicinity of the equilibrium $(u_0,v_0)$,
while the general situation is more accurately described by a nonlinear system as in (<ref>).
For practical purposes, the nonlinear sources $f$ and $g$ would force a bound on the densities $u$ and $g$
and possibly favor the convergence of the solution for large times to steady solutions $u=u(x)$ and $v=v(x)$ of

\begin{equation}\label{SYS:READI:STEAD}
\begin{dcases}
&\mu\Delta u(x)+f(u(x),v(x))=0,\\
&\nu\Delta v(x)+g(u(x),v(x))=0.
\end{dcases}
\end{equation}

That is, for small times, the linear mechanism identified in (<ref>) is helpful to detect how diffusion can drive
instability and place the kinetics of the system out of the “trivial” state $(u_0,v_0)$; then, at a longer time scale,
the nonlinear structure of (<ref>) becomes instrumental to confine the solution and lead it towards spatially
inhomogeneous patterns, as described by the steady solutions of (<ref>)
(as another option, the evolution of the equation may be stopped after a certain time in case
the release of the chemical substances stops; this could be the case in which
the pattern is formed at an embryonic stage for the animal due to chemical substances that are released only
during specific periods of the early stage development of an organism).

We also remark that, for the sake of simplicity, here we confined ourselves to the case in which the spatial domain
is a circle (i.e., the real line with periodic assumptions in $x$): in general, if one considers more complicated domains
(say, closer to biological situations of specific interest) then it is convenient to replace $e^{ijx}$ in (<ref>)
with the eigenfunctions
of the Laplacian in the domain of interest (with the corresponding boundary conditions).
In like manner, the dispersion relations in (<ref>) and (<ref>)
must take into account the corresponding eigenvalues in the place of $-j^2$.
These eigenvalues replace $-j^2$ in the matrix in (<ref>) too, and the diffusion eigenvalues $\lambda_j$ and $\Lambda_j$ must be modified accordingly.
The diffusion eigenvalues corresponding to diffusive instability still correspond to the ones with positive real part,
and one can focus for concreteness on the diffusive eigenvalue with largest real part in (<ref>).
The structural difference in the general case is however that the excited modes, that is the diffusive
eigenvalues with positive real part, depend on the domain. Since their region of positive and negative values
create a visible pattern in the dynamics of the solution (compare with the positive and negative values of the solution
depicted in Figure <ref>), it is conceivable that these regions have a connection with the visible
patterns in nature. For instance, two-dimensional
rectangular domains present eigenfunctions of the type $\cos\left(\frac{2\pi j_1 x_1}{\ell_1}\right)\cos\left(\frac{2\pi j_2 x_2}{\ell_1}\right)$, where $\ell_1$ and $\ell_2$ account for the lengths of the side of the rectangle
and $j_1$, $j_2\in\N$. In this scenario, excited states corresponding to $j_1=0$ give rise to horizontal stripes,
the ones corresponding to $j_2=0$ to vertical stripes, and the ones with $j_1\ne0$ and $j_2\ne0$ to maculate patterns, see
e.g. Figure <ref>.

$\quad$  $\quad$

Level sets of $(x_1,x_2)\mapsto2\cos x_2$, $(x_1,x_2)\mapsto2\cos x_1$ and $(x_1,x_2)\mapsto2\cos x_1\cos x_2$.

With respect to this, notice also that different excite modes correspond to different size of the pattern (such as width of the
or possible elongation of the spots). Of course, the analysis of simple rectangular regions is insufficient to capture the whole complexity of animal patterns: yet, it is suggestive to “approximate” an animal's coat with rectangular regions and
observe how for instance the orientation of stripes “locally” follow the proportion of the approximating rectangles,
see e.g. Figure <ref>
which clearly shows a change between vertical and horizontal stripes patterns in zebras
at the junctions between the body and the legs and between the body and the tail.

Mutually grooming zebras (photo by Duvignau Alain from
Wikipedia, available
under the Creative Commons CC0 1.0 Universal Public Domain Dedication).

Animal patterns may also bifurcate from stripes to spots in different regions of the organism, possibly in response of a variation
of the “local geometry”, compare e.g. with the tail markings of the cheetah in which the typical spots of the animal
become stripes at the tip of the tail, see Figure <ref>.

Caress of the Sphinx, by Fernand Khnopff (Public Domain image from

After <cit.>,
Turing's ideas about morphogenesis have become a cornerstone
in mathematical biology and led to a number of fantastic accomplishments, see e.g.
the mathematical reconstruction of very sophisticated animal patterns
presented in <cit.>.
Also, the name of Turing patternTuring pattern
is nowadays commonly used to denote biological structures
and specifically those arising by the imbalances between diffusion rates of different chemical agents which make a stable system sensitive to perturbations.

See <cit.>, <cit.>, <cit.>
and the references therein for a thorough analysis of diffusion driven pattern formation,
with several remarkable examples inspired by concrete natural phenomena.
See also <cit.> and the reference therein for more information about
reaction-diffusion equations.

§.§ Space invaders

A cane toad (bufo marinus) and the spread of cane toads in Australia from 1940 to 1980 in five-year intervals
(photo by Bill Waller and animated map by  Froggydarb; images from
Wikipedia, licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license).

A topical argument in mathematical biology consists in the study of invasive species and in their territorial colonization ability,
see e.g. Figure <ref>.
To introduce ourselves to this subject, it is first opportune to understand the so-called logistic equation
logistic equation. This model was introduced by Pierre François Verhulst,
see Figure <ref> for his portrait,
and describes the evolution in time of a biological population.
The number of individuals $N(t)$ is supposed to grow at an intrinsic growth rate, parameterized by a given $\rho>0$,
which somewhat accounts for the ideal birth rate of the population ($N$ individuals would give rise to $\rho N$ newborns in the unit of time). Additionally, the population undergoes an intraspecific competition which causes the death of some individuals
due to possible overcrowding (assuming, for instance
that the environmental resources only allow a maximum number of individuals $K$). The combination of these effects
lead to the logistic ordinary differential equation

\begin{equation}\label{HJA:A89MAOL8789892BIOEQLOG}
\frac{dN}{dt}(t)=\rho N(t)\left(1-\frac{N(t)}K\right).
\end{equation}

This equation considers all the population located basically at the same place, but, for many practical purposes,
it is also convenient to describe biological individuals in a spatial environment, in which case the function $N$
depends on the time variable $t$ and also on the space variable $x$. If the population performs some kind of random walk,
as discussed in Section <ref> (see in particular equation (<ref>)), combining (<ref>) with the random diffusive tendency of the population we arrive at the equation

\begin{equation}\label{HJA:A89MAOL8789892BIOEQLOG-2}
\partial_t N(x,t)=c\Delta N(x,t)+\rho N(x,t)\left(1-\frac{N(x,t)}K\right),
\end{equation}

for some diffusion coefficient $c>0$.

We can actually simplify (<ref>) via a simple rescaling: namely, defining

\begin{equation}\label{HJA:A89MAOL8789892BIOEQLOG-4}

equation (<ref>) reduces to

\begin{equation}\label{HJA:A89MAOL8789892BIOEQLOG-3}
\partial_t u=\Delta u+u(1-u),
\end{equation}

which is dubbed[Equation (<ref>)
is named after Sir Ronald Aylmer Fisher <cit.>,
Andrey Nikolaevich Kolmogorov, Ivan Georgievich Petrovsky and Nikolai Semenovich Piskunov <cit.>.

Though biological invasions and
the spread of a genetic trait were
the prime motivations for the study of (<ref>),
nowadays this equation is also intensively studied for its applications
in combustion and flame propagation: that's the unifying power of
the good mathematics!

Besides the introduction of equation (<ref>) and some questionable
studies in eugenics, Fisher is considered one of the founders of population genetics and of modern statistical science.
His name is also linked to the so-called
Fisherian runaway sexual selection mechanism, aiming at an explanation of
some exaggerated, costly and apparently maladaptive male ornamentation in nature (seemingly in conflict with natural selection)
as dictated by persistent female choice (possibly triggered initially by the
ornament signaling greater potential fitness, hence likelihood of leaving more descendants):
the classical example of Fisherian runaway is the elaborate peacock plumage, see Figure <ref>.

A classical story about Fisher's contribution towards the standardization of statistical experiments is the “lady tasting tea” <cit.>. In a nutshell, once Fisher offered a cup of tea to phycologist Blanche Muriel Bristol.
She politely declined it, saying that she preferred the flavor when the milk was poured into the cup before (not after) the tea. So we have an intriguing dilemma: can the order of pouring milk affect the flavor of the tea? and, specifically, can Bristol notice the difference?

To test this, Fisher provided Bristol with eight randomly ordered cups of tea: four of which were prepared by first pouring the tea and then adding milk, and four by first pouring the milk and then adding the tea. Bristol had to identify the different cups.
Fisher had to interpret the data to distinguish the case in which Bristol could actually tell the difference in the pouring order
and the one in which she could identify correctly some cups, but just by chance.

For this,
with a bit of combinatorics, Fisher computed that there are ${{8}\choose{4}}=\frac{8!}{4! (8-4)!}=70$ possible combinations of cups.
Since Bristol was aware that there were four cups of each type, her answer would have included four of each.

There is of course only one possibility of making precisely the right identification: hence achieving a complete identification
just by chance occurs with probability $\frac{1}{70}=0.01428571428...$;
there are instead 16 possibilities of making exactly one error (indeed, one error comes from a single swap of two cups of different type, say the $i$th cup of the first type with the $j$th cup of the second type, with $i$, $j\in\{1,2,3,4\}$).
Hence, there are $1+16=17$ possibilities of making at most one mistake
and therefore, if we adopt a conventional probability criterion of success to be below the $5\%$ threshold, the ability of Bristol to properly categorize the cups of tea would have been confirmed if, and only if, she was able to correctly identify all of them, without making any mistake.

This story is relevant since it became an example of a randomized experiment for checking a “null hypothesis” in statistics and it contributed to a scientific establishment of randomization analysis of experimental data.

For the record, Muriel Bristol allegedly succeeded in classifying all eight cups correctly.

See <cit.> for further reading on the lady tasting tea experiment.

See also Figure <ref>
for a picture of Kolmogorov and footnote <ref>
on page KOLMONOTE for more information about the prominent figure of Kolmogorov.]
in jargon the Fisher-Kolmogorov-Petrovsky-Piskunov equation
(or Fisher-KPP equation for short). Fisher-Kolmogorov-Petrovsky-Piskunov equation
Interestingly, this equation is a reaction-diffusion equation reaction-diffusion equation
in the setting introduced on page REDI6NFHARLROA7KOLM789GIJ7solFUMHDNOJHNFOJED231.

Pierre Verhulst (Public Domain image from

The question of determining the invasivity of biological species is very related to (<ref>).
In this setting, in light of (<ref>), the level $u=0$ corresponds to the absence of the population
and the level $u=1$ to the maximal density permitted by environmental resources.
Remarkably, for any given $v\in[2,+\infty)$

\begin{equation}\label{BEVEPAVKCiuneA3OMAkac}
\begin{split}&
{\mbox{equation~\eqref{HJA:A89MAOL8789892BIOEQLOG-3} admits a solution~$u$ describing the full environmental invasion}}\\ &
{\mbox{of the biological species with velocity~$v$.}}\end{split}\end{equation}

More precisely, given any $v\in[2,+\infty)$ and any $\omega\in \partial B_1$, there exists a smooth function $U:\R\to(0,1)$, with

\begin{equation}\label{BEVEPAVKCiuneA3OMAkac22} \lim_{\tau\to+\infty}U(\tau)=0\qquad{\mbox{and}}\qquad\lim_{\tau\to-\infty}U(\tau)=1,\end{equation}

and a solution $u$ of (<ref>) of the form

\begin{equation}\label{SOLIWA} u(x,t):=U(\omega\cdot x-vt).\end{equation}

The solutions[The concrete interest of the very special solutions in (<ref>) is given by the fact that
under suitable initial conditions, “reasonable” solutions of (<ref>) evolve into a solitary wave,
see <cit.>. Roughly speaking, wells in the initial conditions
are quickly reabsorbed by the evolving dynamics, since regions with $u$ close to $1$ tend to grow
and eat up adjacent regions with $u$ close to $0$, producing, for large times, a monotonic traveling wave with constant speed.
This phenomenon makes solitary waves an essential building block for understanding the dynamics of the
Fisher-Kolmogorov-Petrovsky-Piskunov equation.] of this form are often called “solitary waves”. solitary wave
See Figure <ref> for a representation of a solitary wave.

A blue peacock presenting its feathers (image by from
Wikipedia, available under the Creative Commons Attribution-Share Alike 3.0 Unported).

To establish (<ref>),
we consider[As remarked in <cit.>, when $v:=\frac{5}{\sqrt {6}}=2.04124145232...$,
one can find explicit solutions of (<ref>) in the form

$$ U(\tau)=\frac1{\displaystyle\left(1+C\exp\left(\frac\tau{\sqrt{6}}\right)\right)^{2}},$$

for all $C>0$.

Figure <ref>
was obtained with this function with $C:=1$ and the frames of the second image there are obtained
by taking times $t\in\{0,2,4,6,8\}$.

See also <cit.> for more information about this explicit solution.
We are not aware of other solutions of (<ref>) which
can be written explicitly in a nice and closed form.]
the ordinary differential equation

\begin{equation}\label{BEVEPAVKCiuneA3OMAkac2}
\frac{d^2U}{d\tau^2}(\tau)+v\frac{dU}{d\tau}(\tau)+U(\tau)\big(1-U(\tau)\big)=0.\end{equation}

By setting $V(\tau):=\frac{dU}{d\tau}(\tau)$, we obtain the system of equations

\begin{equation}\label{BEVEPAVKCiuneA3OMAkac3b}
\left(\begin{matrix}
\displaystyle\frac{dU}{d\tau} \\
\\
\displaystyle\frac{dV}{d\tau}
\end{matrix}\right)=\left(\begin{matrix}
V \\
\\
\end{matrix}\right).\end{equation}

The equilibria of (<ref>) are $E_0:=(0,0)$ and $E_1:=(1,0)$.
The linearized dynamics in the vicinity of $E_1$ is described by the matrix $M_1:=\left(
\begin{matrix}
0 & 1\\
1 & -v
\end{matrix}\right)$, which possesses eigenvalues $\frac12\left(-v \pm \sqrt{v^2 + 4}\right)$
with eigenvectors $\left(\frac12 \left(v\pm\sqrt{v^2+4}\right), 1\right)$.

In particular, $E_1$ is a hyperbolic saddle for all values of $v\in\R$. So we consider the unstable manifold
(see e.g. <cit.>) of $E_1$ in the direction $-\left(\frac12 \left(v+\sqrt{v^2+4}\right), 1\right)$.
This provides us with an orbit $(U(\tau),V(\tau))$ such that

\begin{equation}\label{JMNSRAPPORA1}
\lim_{\tau\to-\infty}(U(\tau),V(\tau))=E_1\end{equation}

\begin{equation}\label{JMNSRAPPORA}
\lim_{\tau\to-\infty}
\left(\frac{dU}{d\tau} (\tau),\frac{dV}{d\tau}(\tau)\right)= -
\left(\frac12 \left(v+\sqrt{v^2+4}\right), 1\right),\end{equation}

see Figure <ref>.

Given $v\in[2,+\infty)$, we now consider the planar region

$$ {\mathcal{R}}:=\left\{{\mbox{$(U,V)\in\R^2$ \;s.t.\; $U<1$, \;$V\in\left(-\displaystyle\frac1v,0\right)$\;
and \; $\displaystyle\frac{v+\sqrt{v^2-4}}{2}\,U+V>0$}}\right\}.$$

A traveling wave of a species invading the spacial domain from left to right.

Notice that, in light of (<ref>) and (<ref>),
there exists $\tau_\star\in\R$ such that $(U(\tau),V(\tau))\in{\mathcal{R}}$ for all $\tau\in(-\infty,\tau_\star)$.

We claim that

\begin{equation}\label{CONGFIuynU6mnemt}
{\mbox{$(U(\tau),V(\tau))\in{\mathcal{R}}$\; for all \;$\tau\in(-\infty,+\infty)$.}}\end{equation}

Indeed, suppose not. Then, there exists $\tau_0\ge \tau_\star$ such that $(U(\tau),V(\tau))\in{\mathcal{R}}$ for all $\tau\in(-\infty,\tau_0)$
but either $V(\tau_0)=-1/v$, or $V(\tau_0)=0$
or $U(\tau_0)=1$, or $\frac{v+\sqrt{v^2-4}}{2}\,U(\tau_0)+V(\tau_0)=0$. We exclude all these cases to prove (<ref>).

First, we show that $V(\tau_0)=-1/v$ cannot hold. This is because otherwise $
\frac{dV}{d\tau}(\tau_0)\le0$, but we know from (<ref>)

$$ \frac{dV}{d\tau}(\tau_0)=-vV(\tau_0)-U(\tau_0)\big(1-U(\tau_0)\big)\ge-vV(\tau_0)-\frac14=1-\frac14>0.$$

Similarly, $V(\tau_0)=0$ cannot happen, otherwise

$$ 0\le\frac{dV}{d\tau}(\tau_0)=-vV(\tau_0)-U(\tau_0)\big(1-U(\tau_0)\big)\le-vV(\tau_0)=0.$$

This entails that either $U(\tau_0)=0$ or $U(\tau_0)=1$, whence either $(U(\tau_0),V(\tau_0))=E_0$
or $(U(\tau_0),V(\tau_0))=E_1$, which is impossible, being both $E_0$ and $E_1$ equilibria for the system.

Also $U(\tau_0)=1$ cannot hold true, otherwise, by (<ref>),

$$ 0\le\frac{dU}{d\tau}(\tau_0)=V(\tau_0)\le0,$$

that gives that $(U(\tau_0),V(\tau_0))=E_1$, which is impossible.

Finally, the case $\frac{v+\sqrt{v^2-4}}{2}\,U(\tau_0)+V(\tau_0)=0$ cannot occur as well, otherwise we argue as follows.

$$ F(\tau):=\frac{v+\sqrt{v^2-4}}{2}\,U(\tau)+V(\tau).$$

We note that $F(\tau)>0$ for $\tau\in(-\infty,\tau_0)$ and $F(\tau_0)=0$. As a result, from (<ref>) we arrive at

\begin{eqnarray*} 0&\ge&\frac{dF}{d\tau}(\tau_0)\\&=&
\frac{v+\sqrt{v^2-4}}{2}\,\frac{dU}{d\tau}(\tau_0)+\frac{dV}{d\tau}(\tau_0)\\&=&
\frac{v+\sqrt{v^2-4}}{2}\,V(\tau_0)-vV(\tau_0)-U(\tau_0)+U^2(\tau_0)\\&=&
\frac{\sqrt{v^2-4}-v}{2}\,V(\tau_0)-U(\tau_0)+U^2(\tau_0)\\&=&-
\frac{\big(\sqrt{v^2-4}-v\big)\big(v+\sqrt{v^2-4}\big)}{4}\,U(\tau_0)-U(\tau_0)+U^2(\tau_0)\\&=&-
\frac{(v^2-4)-v^2}{4}\,U(\tau_0)-U(\tau_0)+U^2(\tau_0)\\&=&
\end{eqnarray*}

This gives that $U(\tau_0)=0$, whence $V(\tau_0)=0$, but this contradicts the fact that $E_0$ is an equilibrium for the system.

The proof of (<ref>) is thereby complete.

As a byproduct of (<ref>), we also have that, for all $t\in\R$,

\begin{equation}\label{CONGFIuynU6mnemtLI9}
\end{equation}

Now we claim that

\begin{equation}\label{CONGFIuynU6mnemtLI}
\lim_{\tau\to+\infty}(U(\tau),V(\tau))=E_0.
\end{equation}

Indeed, by (<ref>) we know that $V(\tau)<0$ for all $\tau\in\R$, therefore by (<ref>), $U(\tau)$ is decreasing.
Since also by (<ref>) we have that $U(\tau)\in(0,1)$ for all $\tau\in\R$, we infer that
there exists $U_\infty\in[0,1)$ such that $U(\tau)\to U_\infty$ as $\tau\to+\infty$.

Now, let $\delta>0$, to be taken as small as we wish here below. Using the “dots” to denote derivatives, recalling (<ref>) and (<ref>) we observe that

\begin{equation}\label{KAQUINLtyY93l2}\begin{split}
\le \frac1\delta\left|\int_\tau^{\tau+\delta}\big(V(\tau)-V(\theta)\big)\,d\theta\right|+
\frac1\delta\left|\int_\tau^{\tau+\delta} V(\theta)\,d\theta\right|\\&=
\frac1\delta\left|\int_\tau^{\tau+\delta}\left(\int_\theta^\tau\dot{V}(\zeta)\,d\zeta\right)\,d\theta\right|+
\frac1\delta\left|\int_\tau^{\tau+\delta} \dot{U}(\theta)\,d\theta\right|\\&=
\frac1\delta\left|\int_\tau^{\tau+\delta}\left(\int^\theta_\tau vV(\zeta)+U(\zeta)\big(1-U(\zeta)\big)\,d\zeta\right)\,d\theta\right|+
\frac1\delta\left|U(\tau+\delta)-U(\tau)\right|\\&\le
\frac2\delta\int_\tau^{\tau+\delta}  \left(\int^\theta_\tau \,d\zeta\right)\,d\theta+
\frac1\delta\left|U(\tau+\delta)-U(\tau)\right|\\&=
\delta+
\frac1\delta\left|U(\tau+\delta)-U(\tau)\right|.\end{split}
\end{equation}

From this we arrive at

$$ \limsup_{\tau+\infty}|V(\tau)|\le\limsup_{\tau+\infty}
\delta+
\frac1\delta\left|U(\tau+\delta)-U(\tau)\right|=\delta+
\frac1\delta\left|U_\infty-U_\infty\right|=\delta$$

and accordingly, by choosing $\delta$ arbitrarily small,

$$ \lim_{\tau+\infty} V(\tau)=0.$$

As a result, we have that $(U(\tau),V(\tau))\to(U_\infty,0)\in[0,1)\times\{0\}$.
Since the only equilibrium of (<ref>) on $[0,1)\times\{0\}$
is $E_0$, the proof of (<ref>) is complete.

Phase portrait and heteroclinc connection for the dynamical system in (<ref>) with $v:=5/2$.

Hence, in view of (<ref>), (<ref>), (<ref>) and (<ref>),
the trajectory $U(\tau)$ provides a solution of (<ref>) satisfying the asymptotics in (<ref>). Then, by (<ref>),

\begin{eqnarray*}&&
-\partial_t u(x,t)+\Delta u(x,t)+u(x,t)\big(1-u(x,t)\big)\\&&\qquad=
v \frac{dU}{d\tau}(\omega\cdot x-vt)
+\frac{d^2U}{d\tau^2}U(\omega\cdot x-vt)
+U(\omega\cdot x-vt)\big(1-U(\omega\cdot x-vt)\big)=0.
\end{eqnarray*}

This shows that $u$ is the desired solitary wave solution of the Fisher-Kolmogorov-Petrovsky-Piskunov equation in (<ref>).

For completeness, we show also that biologically interesting solitary waves of the Fisher-Kolmogorov-Petrovsky-Piskunov equation
do not exist when $v\in[0,2)$, since

\begin{equation}\label{1343g43678EPAVKCiuneA3OMAkI536}
\begin{split}&
{\mbox{when~$v\in[0,2)$ any solution~$U:\R\to(-\infty,1]$ of~\eqref{BEVEPAVKCiuneA3OMAkac22} and~\eqref{BEVEPAVKCiuneA3OMAkac2}}}\\&{\mbox{necessarily attains also negative values.}}\end{split}\end{equation}
We recall indeed that biologically interesting values of the solutions are the ones modeling a population density,
hence they are confined in the interval $[0,1]$, due to (<ref>).

To check the claim in (<ref>),
suppose by contradiction that a solution $U:\R\to(0,1]$ of (<ref>) and (<ref>)
exists for some $v\in[0,2)$ and let $V:=\frac{dU}{d\tau}$ to find the system of equations in (<ref>).

We claim that, for all $\tau\ge0$,

\begin{equation}\label{KAQUINLtyY93l}
|V(\tau)|< |V(0)|+\frac1v.
\end{equation}

Indeed, suppose not. Then, there exists $\underline\tau>0$ such that $|V(\tau)|<|V(0)|+\frac1v$
for all $\tau\in[0,\underline\tau)$ and $|V(\underline\tau)|=|V(0)|+\frac1v$. Two cases can occur,
either $V(\underline\tau)\ge0$ or $V(\underline\tau)<0$. In the first case, we have that

$$ 0\le\frac{dV}{d\tau}(\underline\tau)=-vV(\underline\tau)-U(\underline\tau)\big(1-U(\underline\tau)\big)
\le-1,$$

which is a contradiction.

If instead $V(\underline\tau)<0$, then

$$ 0\ge\frac{dV}{d\tau}(\underline\tau)=-vV(\underline\tau)-U(\underline\tau)\big(1-U(\underline\tau)\big)=

which is a contradiction too, hence (<ref>) is proved.

Using this and computing as in (<ref>), for all $\tau\ge0$ and $\delta>0$ we arrive at

\begin{eqnarray*}
\frac1\delta\left|\int_\tau^{\tau+\delta}\left(\int^\theta_\tau vV(\zeta)+U(\zeta)\big(1-U(\zeta)\big)\,d\zeta\right)\,d\theta\right|+
\frac1\delta\left|U(\tau+\delta)-U(\tau)\right|\\&\le&
C\delta +
\frac1\delta\left|U(\tau+\delta)-U(\tau)\right|,
\end{eqnarray*}

with $C:=v\left(|V(0)|+\frac1v\right)+1$.

|V(\tau)|\le C\delta +
\frac1\delta\left|0-0\right|=C\delta$$

and accordingly, taking $\delta$ as small as we wish, we find that $V(\tau)\to0$ as $\tau\to+\infty$.

In this way, we have established that $(U(\tau),V(\tau))\to(0,0)$ as $\tau\to+\infty$.
Hence, we pick $\eta>0$, to be chosen conveniently small in what follows and we find $\tau_\eta\in\R$ such that $|(U(\tau),V(\tau))|\le\eta$ for all $\tau\ge\tau_\eta$. Up to a time shift, i.e. up to replacing $(U(\tau),V(\tau))$ with $(U(\tau-\tau_\eta),V(\tau-\tau_\eta))$, we can assume from now on that $\tau_\eta=0$ and write that

\begin{equation}\label{L-X2343545y678we34eS}
\big|(U(\tau),V(\tau))\big|\le\eta\qquad{\mbox{for all }}\,\tau\ge0.
\end{equation}

Now we perform an argument related to the nonlinear stability of equilibria in planar dynamical systems.
Namely, we take $\varrho(\tau)>0$ and $\varphi(\tau)\in\R$ such that

$$ U(\tau)=\varrho(\tau)\,\cos(\varphi(\tau))\qquad {\mbox{and}}\qquad

Actually, since $U(\tau)>0$, we have that

\begin{equation}\label{922i3k2e944racrhfo}
{\mbox{we can pick~$\varphi(\tau)$ in the interval }}\,\left(-\frac\pi2,\frac\pi2\right).\end{equation}

Thus, denoting by the dot the derivative with respect to $\tau$,

$$ \dot U=\dot\varrho\cos\varphi- \varrho\dot\varphi\sin\varphi
\qquad {\mbox{and}}\qquad
\dot V=\dot\varrho\sin\varphi+\varrho\dot\varphi\cos\varphi.$$

From this we arrive at

\begin{eqnarray*}
\dot U\sin\varphi-\dot V\cos\varphi&=&
\big(\dot\varrho\cos\varphi- \varrho\dot\varphi\sin\varphi\big)\sin\varphi-
\big(\dot\varrho\sin\varphi+\varrho\dot\varphi\cos\varphi\big)\cos\varphi
\\&=&
\dot\varrho\sin\varphi\cos\varphi- \varrho\dot\varphi\sin^2\varphi-
\dot\varrho\sin\varphi\cos\varphi-\varrho\dot\varphi\cos^2\varphi
\\&=&-\varrho\dot\varphi.
\end{eqnarray*}

Since, by (<ref>),

\begin{eqnarray*}
\dot U\sin\varphi-\dot V\cos\varphi&=&
\varrho\cos^2\varphi-\varrho^2\cos^3\varphi\\&=&\varrho+\frac{v\varrho}2\sin(2\varphi)
\end{eqnarray*}

we conclude that

\[ -\varrho\dot\varphi=\varrho+\frac{v\varrho}2\sin(2\varphi)

and accordingly

\[ -\dot\varphi=1+\frac{v}2\sin(2\varphi)

In particular, using (<ref>) and the assumption that here $v\in[0,2)$,

\[ -\dot\varphi\geq 1-\frac{v}2-\varrho\ge 1-\frac{v}2-\eta
\ge\frac{1-\displaystyle\frac{v}2}{2}>0,\]

as long as $\eta$ is small enough.
This gives that there exists $\tau_\star>0$ such that $\varphi(\tau_\star)<-\frac{\pi}{2}$,
in contradiction with (<ref>).
This completes[Here is an alternative argument
to prove (<ref>), relying on the additional simplifying assumptions that

$$ \lim_{\tau\to+\infty}\dot{U}(\tau)=0$$

and the following limit exists

$$ \ell:=\lim_{\tau\to+\infty} \frac{\ddot{U}(\tau)}{\dot{U}(\tau)}.$$

From this, if $U(\tau)\in(0,1)$ for all $\tau\in\R$ and $U(\tau)\to0$ as $\tau\to+\infty$, we can apply L'Hôpital's Rule and deduce that

$$ \lim_{\tau\to+\infty} \frac{\dot{U}(\tau)}{U(\tau)}=
\lim_{\tau\to+\infty} \frac{\ddot{U}(\tau)}{\dot{U}(\tau)}=\ell.$$

As a result, dividing (<ref>) by $U(\tau)$ and taking the limit we find that

\begin{eqnarray*}
&=&\lim_{\tau\to+\infty}\frac{\ddot{U}(\tau)}{\dot U(\tau)}\;
\frac{\dot{U}(\tau)}{U(\tau)}+v\frac{\dot{U}(\tau)}{U(\tau)}+1\\
&=&\left( \ell+\frac{v}2\right)^2+1-\frac{v^2}4
\\&\ge&1-\frac{v^2}4,
\end{eqnarray*}

which gives that $v^2\ge4$, hence (in the convention that $v\ge0$) necessarily $v\ge2$ to allow invading fronts.

These computations are also interesting since for more general equations such as $\ddot{U}+v\dot{U}+f(U)=0$ with $f$ smooth, $f(r)>0$ for all $r\in(0,1)$, $f(0)=f(1)=0$ and $f'(1)<0$, a calculation as above leads to

\begin{eqnarray*}&& 0=\ell^2+v\ell+\lim_{\tau\to+\infty}\frac{f(U(\tau))}{U(\tau)}=
\ell^2+v\ell+f'(0)=\left( \ell+\frac{v}2\right)^2+f'(0)-\frac{v^2}4\ge f'(0)-\frac{v^2}4,
\end{eqnarray*}

leading to the necessary condition for invasion $v\ge 2\sqrt{f'(0)}$.

This is quite interesting, since it also points out that the speed of the invasive front is dictated by its low density fringes not by the high density regions (namely, the velocity of invasion is completely determined by the forcing term $f$
and specifically by the values of $f$ in the vicinity of $0$, not in the vicinity of $1$).

A function $f$ with the above structure is often called a bistable nonlinearity bistable nonlinearity.

For further details on bistable nonlinearities and plane wave solutions, see <cit.>.]
the proof of (<ref>).

For additional information on the Fisher-Kolmogorov-Petrovsky-Piskunov, see <cit.>,
<cit.>, <cit.>
and the references therein.

§.§ Equations from hydrodynamics

A classical arena for partial differential equations is provided by fluid mechanics and hydrodynamics.
Let us present some of the situations that naturally occur in this framework.
Let us consider a fluid with density $\rho(x,t)$, vectorial velocity $v(x,t)$ and pressure $p(x,t)$.
The position of a parcel of fluid at time $t$ is given by $x(t)$ and therefore such a parcel
travels according to the velocity prescription

\begin{equation}\label{0uojfSANMFOSJFEL}
\dot x(t)=v(x(t),t).
\end{equation}

The title page to Euler's original article “Principes généraux du mouvement des fluides”,
published in Mémoires de l'Académie des Sciences de Berlin in 1757 (Public Domain source from
Internet Archive).

Also, the total mass of a fluid element occupying a region $\Omega$ of the space at a given time $t_0$
is given by the quantity $\int_{\Omega} \rho(x,t_0)\,dx$. This amount of fluid will move around in a small time $\tau$
possibly occupying a different region, that we name $\Omega_\tau$ which collects all the evolution trajectories $x(t_0+\tau)$
of the fluid parcels such that $x(t_0)\in\Omega$. Since we are assuming that matter is neither created nor disappear in this process
we deduce a conservation law given by the equation

$$ \int_{\Omega} \rho(x,t_0)\,dx=\int_{\Omega_\tau} \rho(x,t_0+\tau)\,dx$$

or equivalently

\begin{equation}\label{0ojCOMAS} \frac{d}{d\tau}\int_{\Omega_\tau} \rho(x,t_0+\tau)\,dx=0.\end{equation}

Notice that we can change variable in the previous integral simply by flowing back the fluid.
Specifically, given a point $q\in\R^n$ we denote by $\Phi^t(q)$ the evolution according to the fluid vector field $v$
with initial position at time $t_0$ equal to $q$, or more formally we define $\Phi^t(q)$ as the solution of the following
Cauchy problem for ordinary differential equations:

$$ \begin{dcases}
\frac{d}{dt} \Phi^t(q)=v\big(\Phi^t(q),t),\\
\Phi^{t_0}(q)=q.
\end{dcases}$$

This is nothing else as a rewriting of (<ref>), but we are keeping track here explicitly of the initial position.
With this notation, we have that $\Omega_\tau=\{\Phi^\tau(y)$, $y\in\Omega\}$.
As a result, one can perform the change of variable

$$x=\Phi^\tau(y)=y+\tau v(y,t_0)+o(\tau),$$

which produces

$$ |\det \partial_y x|=1+\tau \div v(y,t_0)+o(\tau)$$

and therefore

\begin{eqnarray*}&&
\int_{\Omega_\tau} \rho(x,t_0+\tau)\,dx\\&=&
\int_{\Omega} \rho(y+\tau v(y,t_0),t_0+\tau)\,\big(1+\tau \div v(y,t_0)\big)\,dy+o(\tau)
\\&=&\int_{\Omega}
\Big(\rho(y,t_0)+\tau\nabla\rho(y,t_0)\cdot v(y,t_0)
+\tau\partial_t \rho(y,t_0)\Big)
\,\big(1+\tau \div v(y,t_0)\big)\,dy+o(\tau)\\&=&\int_{\Omega}
\Big(\rho(y,t_0)+\tau\nabla\rho(y,t_0)\cdot v(y,t_0)
+\tau\partial_t \rho(y,t_0)+
\tau\rho(y,t_0)\div v(y,t_0)\Big)\,dy+o(\tau)\\&=&\int_{\Omega}
\Big(\rho(y,t_0)
+\tau\partial_t \rho(y,t_0)+
\tau\div\big(\rho(y,t_0)\, v(y,t_0)
\big)\Big)\,dy+o(\tau).
\end{eqnarray*}

This and (<ref>) formally lead to

\begin{equation}\label{D23q4wy5OPACLPASIJMDNIAKS279uyt3M}
\begin{split}0=\left.
\frac{d}{d\tau}\int_{\Omega_\tau} \rho(x,t_0+\tau)\,dx\right|_{\tau=0}=
\int_{\Omega}
\Big(\partial_t \rho(y,t_0)+\div\big(\rho(y,t_0)\, v(y,t_0)\big)\Big)\,dy.\end{split}\end{equation}

Since this holds true for all domains $\Omega$ and all times $t_0$, we have obtained the equation

\begin{equation}\label{OJS-PJDN-0IHGDOIUGDBV02ujrfMTE}
\partial_t \rho+\div(\rho v)=0.
\end{equation}

This is called[It is interesting to observe that in the framework of equation (<ref>)
one can reinterpret the divergence term in (<ref>)
as a transport term with velocity $\nabla w$ (that is, the effect of chemotaxis
is to produce a transport of the biological population with a velocity proportional to the gradient of the chemical attractant).

Moreover, in chemistry one often models the spreading of some substrate with concentration $\rho$
in a fluid with velocity $v$: in this context equation (<ref>)
is often called “convection equation”convection equation.]
in jargon “mass transport equation”mass transport equation
(or simply “transport equation”) transport equation
or “continuity equation”continuity equation.

Of course, transport of mass, corresponding to the conservation of matter, is not the only ingredient to accurately
describe the dynamics of a fluid, hence one has to complement (<ref>) with additional information. In particular, one can take into account Newton's Law regarding momentum conservation.
In this setting, the change of fluid momentum must correspond to the forces acting on the fluid
(for simplicity, we assume here that the fluid is only subject to gravity and to its own pressure force
with no other external forces to be exerted).
Notice that the momentum corresponding to the mass of the fluid in a region of space $\Omega$
is given by the quantity $\int_\Omega \rho(x,t) v(x,t)\,dx$.
The corresponding gravity force is provided by $-ge_n\int_\Omega \rho(x,t)\,dx$, being $g$ the gravity acceleration constant
(and $e_n=(0,\dots,0,1)$, supposing the gravity acting downwards in the vertical direction).
As for the pressure force, it arises from $p$ in the normal direction
along the surface $\partial\Omega$, with a conventional
minus sign (so that a resistance to velocity increasing occurs when moving towards regions with high pressures, while
low pressure regions produce a sucking effect). These considerations translate
Newton's Law about momentum conservation into the formula

\begin{equation} \label{90uyoihfs83wEUSlNOE}
\frac{\partial}{\partial t}\int_\Omega \rho(x,t) v(x,t)\,dx=
-ge_n\int_\Omega \rho(x,t)\,dx-\int_{\partial\Omega} p(x,t) \nu(x)\,d{\mathcal{H}}^{n-1}_x,\end{equation}

being, as usual, $\nu(x)$ the unit external normal at $x\in\partial\Omega$.

Furthermore, in light of (<ref>)
(applied here to $\rho v_j$ instead of $\rho$) we have that, for all $j\in\{1,\dots,n\}$,

\begin{eqnarray*}
&&\frac{\partial}{\partial t}\int_\Omega \rho(x,t) v_j(x,t)\,dx\\
&=& \int_{\Omega}
\Big(\partial_t (\rho v_j)(x,t)+\div\big((\rho v_j)(x,t)\, v(x,t)\big)\Big)\,dx\\&=&
\int_{\Omega}
\Big(\rho (x,t)\partial_t v_j(x,t)+\rho(x,t) v(x,t)\cdot\nabla v_j(x,t)\Big)\,dx,
\end{eqnarray*}

where (<ref>) has been exploited in the latter step.

Thus, using the notation $v\cdot\nabla$ to denote the operator $\displaystyle\sum_{i=1} v_i\partial_i$,
possibly applied to a vector valued function componentwise, we can write that

\begin{eqnarray*}
&&\frac{\partial}{\partial t}\int_\Omega \rho(x,t) v(x,t)\,dx=
\int_{\Omega}
\Big(\rho (x,t)\partial_t v(x,t)+\rho(x,t) \big(v(x,t)\cdot\nabla\big) v(x,t)\Big)\,dx.
\end{eqnarray*}

From this and (<ref>)
we obtain that

\begin{equation} \label{90uyoihfs83wEUSlNOE-VXBNC}\begin{split}&
\int_{\Omega}
\Big(\rho (x,t)\partial_t v(x,t)+\rho(x,t) \big(v(x,t)\cdot\nabla\big) v(x,t)\Big)\,dx\\&\qquad\qquad=
-ge_n\int_\Omega \rho(x,t)\,dx-\int_{\partial\Omega} p(x,t) \nu(x)\,d{\mathcal{H}}^{n-1}_x.\end{split}\end{equation}

Also, from the Divergence Theorem, for all $j\in\{1,\dots,n\}$,

$$ \int_{ \Omega} \partial_j p(x,t)\,dx=
\int_{ \Omega} \div\big( p(x,t) e_j\big)\,dx
=\int_{\partial\Omega} p(x,t) e_j\cdot\nu(x)\,d{\mathcal{H}}^{n-1}_x$$

and in consequence

\[ \int_{ \Omega} \nabla p(x,t)\,dx=\int_{\partial\Omega} p(x,t) \nu(x)\,d{\mathcal{H}}^{n-1}_x.\]

From these considerations and (<ref>)
we arrive at

\Big(\rho (x,t)\partial_t v(x,t)+\rho(x,t) \big(v(x,t)\cdot\nabla\big) v(x,t)\Big)\,dx=
-ge_n\int_\Omega \rho(x,t)\,dx-\int_{ \Omega} \nabla p(x,t)\,dx$$

and accordingly

\begin{equation}\label{EUMSD-OS-D}
\rho \partial_t v+\rho (v\cdot\nabla) v=
-ge_n \rho-\nabla p.
\end{equation}

The system containing the mass conservation equation in (<ref>)
and the momentum balance equation in (<ref>) (plus possibly a constitutive relation
linking pressure and density) constitutes what in jargon[The fountainhead of these equations
came indeed from the work of Leonhard Euler,
see Figure <ref>.

See also Figure <ref>
for Euler's iconic portrait (by Jakob Emanuel Handmann).
Many of us wonder why Euler is depicted in a pijama with a towel on his head.
Arguably, it was not a pijama but a fashionable banyan, not a towel but a silk cloth
and possibly at that time it was believed that loose and informal
dresses (the banyan, no wig)
contribute to exercise of the faculties of the mind (and indeed we also find it much more comfortable
to do mathematics wearing very informal clothes).

Note that Euler is portrayed as facing left,
possibly to mitigate the impact on the spectator of a right eyelid ptosis and a divergent strabismus.
At age 31, Euler became almost blind in his right eye <cit.>
(and this might be a possible explanation for which Euler was nicknamed “the cyclops”
by Frederick II, King of Prussia).

This partial loss of vision did not discourage Euler from producing the finest possible
mathematics (actually, he allegedly stated “Now I will have fewer distractions”).
At age 59, a surgical restoration for a cataract in his left eye rendered Euler
almost totally blind. Yet, Euler's passion and talent for mathematics remained undefeated
and he kept producing a vast number of impressive works which changed the course of mathematics
till he died at age 76. In his eulogy, Marquis de Condorcet wrote
“il cessa de calculer et de vivre” (French for “he ceased to calculate and to live”).]
are called “Euler fluid flow equations”Euler equation.

Portrait of Leonhard Euler (Public Domain image from

In a sense, these equations are obtained for a rather ideal situation and several modifications
of the previous setting can be performed to account for more complex models. Among these
modifications, one of the most popular one consists in accounting for some friction between the fluid molecules
which creates a “viscosity”viscosity resisting to the change of fluid velocity. More specifically,
let us quantify[This idea of confronting the pointwise value of a function
with the average nearby will be extensively used in Chapter <ref> and it will be
one of the leitmotifs of the study of harmonic functions.]
the difference of velocity between a fluid parcel at $x$ and that of the parcels nearby
by the difference of $v(x,t)$ and its average in a small ball, say of radius $h$ centered at $x$, namely

$$ {\mathcal{D}}(x,t):= v(x,t)-\fint_{B_h(x)}v(y,t)\,dy.$$

By exploiting (<ref>) and polar coordinates, we find that

\begin{equation*}
\begin{split}
{\mathcal{D}}(x,t)&= \fint_{B_h(x)}\big( v(x,t)-v(y,t)\big)\,dy\\
&=-\frac{1}{|B_h|} \int_0^h \left( \int_{\partial B_1}
r^{n-1} \big(v(x+r\omega,t)-v(x,t)\big)\,d{\mathcal{H}}^{n-1}_\omega\right)\,dr\\
&=-\frac{|B_1|}{|B_h|} \int_0^h
r^{n-1} \left(
\frac{c h^2}2 \Delta v(x,t) +o(h^2)
\right)\,dr\\
&=-\mu h^2  \Delta v(x,t) +o(h^2),
\end{split}\end{equation*}

for some $\mu>0$.

Then, we suppose that the velocity of the fluid is reduced when ${\mathcal{D}}$ is positive (since the parcel at $x$
is faster than the ones in its vicinity, and therefore it is dragged back by them)
and is enhanced when ${\mathcal{D}}$ is negative (since the parcel at $x$
is slower than the ones in its vicinity, and therefore it is pulled forward by them).
For simplicity, we can therefore assume that such an additional acceleration term is proportional to ${\mathcal{D}}$
(say, to make it a finite quantity in the limit, proportional to ${\mathcal{D}}/h^2$). By taking $h\searrow0$,
we thus obtain an additional acceleration (or deceleration) given by a term of the form $\mu \Delta v$.
By incorporating this correction into (<ref>) we thus find the equation

\begin{equation}\label{EUMSD-OS-DNS}
\rho \partial_t v+\rho (v\cdot\nabla) v=\mu\Delta v
-ge_n \rho-\nabla p,
\end{equation}

which is called in jargon[Equation (<ref>) (possibly complemented with
other structural equations of the fluid)
come from Claude Louis Marie Henri Navier and Sir George Gabriel Stokes, 1st Baronet.
Our experience with the Laplacian as a “democratic” operator (recall the discussion on page DEEFFNVST)
may suggest that the Laplacian in (<ref>) is providing a smoothing effect.
This may be the case, but our present understanding of the solutions of (<ref>)
is unsatisfactory and highly incomplete. In particular, it is not yet known if smooth (say,
globally defined and meeting some natural conditions)
solutions always exist. If the reader finds an answer to this question,
they do not only earn immortal fame among the circle of PDEs enthusiasts, but also obtain a
substantial amount of money, since the problem is listed in the Millennium Prize Problems
and a $10^6$ US $ award is offered to the first person providing a solution
(the case of  incompressible fluids in $\R^3$ is hard enough). See e.g. <cit.> for further details about the most difficult way to become rich.]
the Navier-Stokes equationNavier-Stokes equation.

We stress that (<ref>) is a vectorial equation (or, equivalently, a system of equations in each component of $v$). There are also interesting situations in which (<ref>) reduces to a scalar equation,
such as the one in which the fluid is constrained within an infinite pipe with a given direction. For instance,
to make things as simple as possible,
let us suppose that the fluid is incompressibleincompressibility: hence, we have that the density of the fluid remains constant in the flow
(namely, a change in the density over time would imply that the fluid had either compressed or expanded). That is, for incompressible fluids the quantity $\rho(x(t),t)$ must be constant in time, leading to

\begin{equation*} 0=\frac{d}{dt}\rho(x(t),t)=\nabla\rho(x(t),t)\cdot \dot x(t)+\partial_t\rho(x(t),t)=
\nabla\rho(x(t),t)\cdot v(x(t),t)+\partial_t\rho(x(t),t).\end{equation*}

Combining this with (<ref>),

\begin{equation*} \rho\div v=\div(\rho v)-\nabla\rho\cdot v=\div(\rho v)+\partial_t\rho(x(t),t)=0,\end{equation*}

\begin{equation}\label{INCOMPRE}
\div v=0,\end{equation}

that is the velocity is a divergence-free vector field.

A horizontal pipe in direction $e_1$ with cross section $\Omega$.

For this reason, recalling (<ref>), stationary solutions of incompressible viscous fluids
are solutions of

\begin{equation}\label{EUMSD-OS-DNSSER}
\begin{dcases}&
\mu\Delta v-\rho (v\cdot\nabla)v-ge_n \rho-\nabla p=0,\\&
\div v=0.
\end{dcases}
\end{equation}

Suppose now that the above fluid moves in parallel streamlines
through a horizontal straight pipe of a given cross sectional form.
Up to rotations which leave invariant the vertical directions, we can assume that the pipe has the form $\R\times\Omega$
for some $\Omega\subset\R^{n-1}$ and $v=(u,0,\dots,0)$ for some scalar function $u$ (that is, the pipe directing the fluid is oriented towards the $e_1$-direction, see Figure <ref>). In this situation,
we have that $\Delta v=\Delta u \,e_1$,
$(v\cdot\nabla)v=u\partial_1u\,e_1$ and $\div v=\partial_1u$.
Consistently with these observations,
(<ref>) reduces to

\begin{equation}\label{EUMSD-OS-DNSSE234R}
\begin{dcases}&
\mu\Delta u\,e_1-\rho u\partial_1u\,e_1-ge_n \rho-\nabla p=0,\\&
\partial_1u=0.
\end{dcases}
\end{equation}

Taking the components $\{2,\dots,n-1\}$ of the first equation in (<ref>) we see that $\partial_i p=0$
for each $i\in\{2,\dots,n-1\}$, giving that the fluid pressure is constant in these directions.
Taking the last component instead, we find $\partial_n p=-g\rho$.

The most interesting information however comes from the first component of the first equation in (<ref>) since
this provides an elliptic equation for $u:\Omega\to\R$. Indeed, focusing on this aspect, we infer from (<ref>) that $u$ is a solution of

\begin{equation*}
\begin{dcases}&
\mu\Delta u-\rho u\partial_1u-\partial_1 p=0,\\&
\partial_1u=0.
\end{dcases}
\end{equation*}

Hence, plugging the second equation into the first one,

\begin{equation}\label{EUMSD-OS-32456i7DNSSE234R}
\mu\Delta u=\partial_1 p \qquad{\mbox{ in }}\,\Omega.
\end{equation}

A special case of interest, which constitutes one of the main motivations for the theory that will be presented
in Sections <ref> and <ref>, is when:

\begin{equation}\label{NK090SOJDKHJNZ}\begin{split}
&{\mbox{the pressure rate is constant,}}\\
&{\mbox{the velocity on the boundary of the pipe is constant}}\\
&{\mbox{and the tangential stress on the pipe is constant.}}
\end{split}\end{equation}

The problem discussed in detail in Sections <ref> and <ref>, and dating back
to James Serrin <cit.>, is precisely
to identify the possible shape of the pipe that allows such special configurations (that is, to determine all the cross sections $\Omega$ that give rise to solutions of (<ref>)
which are compatible with the prescriptions in (<ref>)):
this is certainly a problem of great practical importance, since ideally one would like to project pipes which
do not produce excessive tangential stress, to reduce the wearing of the material and the consequential formation of defects, holes
and licking (for this, maintaining a constant stress along the pipe could be a nice way to balance
the wearing somewhat uniformly along the whole surface of the pipe).

To efficiently address a question of this type, we first need to translate the physical
requests described in (<ref>)
into a mathematical framework. In order to achieve this,
we point out that the prescription that the pressure rate is constant simply says that $\partial_1 p=c_1$, for some $c_1\in\R$.

Also, the assumption that the velocity on the boundary of the pipe is constant means that $u=c_2$ along $\partial\Omega$,
for some $c_2\in\R$. As a matter of fact, since equation (<ref>) remains invariant
if we replace $u$ with $u-c_2$, we can simply assume[Equivalently, one can argue that the invariance of the form of the description of physical problems
among mutually translating reference frames allows us to choose an inertial frame of reference
moving at constant speed $c_2$. In this frame the previous prescription reduces to that
the velocity on the boundary of the pipe is equal to zero.]
that $c_2=0$, hence $u=0$ along $\partial\Omega$.

As for the stress on the pipe, we can suppose that this is due to the fact that fluid parcels near the pipe may move faster or slower
that the ones directly adhering to the pipe (that is, in the normalization above, the fluid parcel in the vicinity of the pipe
could have strictly positive or strictly negative velocity). We can additionally suppose that the wearing of the pipe is directly
influenced by this change of velocity (say, some particles of the pipe get removed by the fluid parcels moving either forwards or backwards). In this spirit, one can model the stress along the pipe as proportional to the normal derivative of the fluid velocity.
With this, the prescription of having a constant stress along the pipe is translated into $\partial_\nu u=c_3$ along $\partial\Omega$,
for some $c_3\in\R$.

In consideration of this, and normalizing constants for simplicity,
we write (<ref>) and (<ref>) as

\begin{equation}\label{EUMSD-OS-32456i7DNSSE234R-BNM}
\begin{dcases}
\Delta u=1&{\mbox{ in }}\Omega,\\
u=0&{\mbox{ on }}\partial\Omega,\\
\partial_\nu u=c&{\mbox{ on }}\partial\Omega.
\end{dcases}\end{equation}

This equation will be indeed[We observe that (<ref>)
does possess a solution when $\Omega$ is a ball. Specifically, given $x_0\in\R^n$ and $r>0$, the function $u(x):=\frac{|x-x_0|^2-r^2}{2n}$ is a solution of (<ref>)
in the ball $B_r(x_0)$
(and note that this solution corresponds to our physical intuition of the problem,
since it corresponds to a fluid reaching its maximal speed at the center of a circular pipe, with the velocity of the fluid being minimal at the boundary of the pipe, due to the viscous effects).

We will discuss in Sections <ref> and <ref>
whether this is the only type of possible solution
or whether there are domains different from the ball which allow different kinds of solutions.]
the starting point of the topics presented in
Sections <ref> and <ref>. A different model leading to the same mathematical
problem will be also presented in Section <ref>.

We end this section by recalling a simple, but helpful, use of the “principle of inertia” in the fluid dynamics setting.
Namely, a case of special interest is provided by a rigid body (say, with a given shape described by a set $\Omega$) moving in a fluid with constant speed. In this situation, the equations presented in this section hold outside the moving domain described, for instance, by $\Omega(t):=\{x\in\R^3$ s.t. $x+v_0 t\in\Omega\}$, for a given constant vector $v_0\in\R^3$.
Since fluid dynamics is already very difficult in a given domain, it is usually more agreeable to change the inertial frame to describe the motion by following the moving body (and we expect by Newton's First Law that this new choice of coordinates does not alter the physical description of the phenomena).

More explicitly, it comes in handy to define

\begin{equation}\label{ISUJONPRENsyhfngTNAbaujnfRGbhndfIKAR}
\widetilde v(x,t):=v(x-v_0 t,t)+v_0,\qquad
\widetilde p(x,t):=p(x-v_0 t,t)\qquad{\mbox{and}}\qquad\widetilde \rho(x,t):=\rho(x-v_0 t,t)\end{equation}
and observe that
equations (<ref>), (<ref>), (<ref>)
and (<ref>) are all preserved by this transformation, since $\div\widetilde v=\div v$,

\begin{eqnarray*}
\partial_t \widetilde\rho+\div(\widetilde\rho \widetilde v)&=&
\partial_t \rho-v_0\cdot\nabla\rho+\div\big(\rho( v+v_0)\big)\\&=&\partial_t \rho+\div(\rho v)
\end{eqnarray*}

\begin{eqnarray*}&&
\widetilde\rho \partial_t \widetilde v+\widetilde\rho (\widetilde v\cdot\nabla) \widetilde v-\mu\Delta \widetilde v+ge_n \widetilde\rho+\nabla \widetilde p\\&=&
\rho \big(\partial_t v-v_0\cdot\nabla v\big)
+\rho \big((v+v_0)\cdot\nabla) v-\mu\Delta v+ge_n \rho+\nabla p\\&=&
\rho \partial_t v
+\rho v\cdot\nabla v-\mu\Delta v+ge_n \rho+\nabla p.
\end{eqnarray*}

Also, the new domain of reference is now simply the complement of $\Omega$.

Notice that the transformation in (<ref>)
is simply replacing the situation of a rigid body in constant speed motion in a fluid at rest at infinity with the one of a still body in a fluid with constant speed at infinity. This is precisely the idea leading to the construction of wind tunnels: to replicate the aerodynamic interactions between a moving object and the surrounding air, it is common to construct large tubes with air blowing through them against a static model of the object.

In the forthcoming Sections <ref> and <ref>, where we will describe rigid bodies moving at a constant speed, we will tacitly assume to have performed the transformation in (<ref>) to reduce ourselves to the case of still objects and domains that do not vary with time.

§.§ Irrotational fluids

Among all possible fluids, a class deserving some special attention is provided by those which present “no vortex”.
In light of Stokes' Theorem, this notion is made precise by the mathematical concept of curl, that is
we say that a fluid is irrotational if the curl of its velocity field vanishes.

A natural question in this setting is whether “vortexes can be produced out of nothing”:
for instance if the initial conditions of a fluid present no vortexes, is it possible that vortexes arise at a later stage?

To increase our familiarity with multivariate calculus and partial differential equations,
we show that this is impossible, at least for inviscid, incompressible and barotropic flows
(in this setting, a flow is said to be barotropic if its density is a function of the pressure only, say $\rho=g(p)$
for some positive function $g$). More precisely, we consider a barotropic solution in the absence of gravity of
the momentum balance equation in (<ref>)
satisfying the incompressibility condition in (<ref>), that is

\begin{equation}\label{CURLDESE1}
\begin{dcases}
\rho \partial_t v+\rho (v\cdot\nabla) v=-\nabla p,\\
\div v=0,\\
\rho=g(p)
\end{dcases}\end{equation}

and we show that

\begin{equation}\label{CURLDESE}
{\mbox{if~$\curl v=0$ at time~$t=0$, then~$\curl v=0$ at every time~$t$.}}\end{equation}

To achieve this goal, we first recall the vector calculus identity, valid for all smooth vector fields $V=(V_1,V_2,V_3):\R^3\to\R^3$,

\begin{equation}\label{VCISTPOK}
(V\cdot\nabla )V=\frac{1}{2}\nabla |V|^2-
V\times (\curl V),\end{equation}

having denoted by $\times$ the vector product operation.
To prove (<ref>), up to exchanging the order of the coordinates, we can concentrate our attention on the first
coordinate: hence we aim at proving that

\begin{equation}\label{VCISTPOK-7}
\sum_{j=1}^3 V_j\partial_j V_1=\frac{1}{2}\partial_1 |V|^2-
\big(V\times (\curl V)\big)_1.\end{equation}

To check this, we calculate:

\begin{eqnarray*}&&
\big(V\times (\curl V)\big)_1-\frac{1}{2}\partial_1 |V|^2\\&=&V_2(\curl V)_3-V_3(\curl V)_2-
V\cdot\partial_1 V\\
&=& V_2\big( \partial_1 V_2-\partial_2 V_1\big) - V_3 \big( \partial_3 V_1-\partial_1 V_3\big)-V_1\partial_1 V_1
-V_2\partial_1 V_2-V_3\partial_1 V_3\\
&=& - V_2\partial_2 V_1 - V_3 \partial_3 V_1-V_1\partial_1 V_1,
\end{eqnarray*}

which proves (<ref>) and thus (<ref>).

We also recall another vector calculus identity, valid for all smooth vector fields $V=(V_1,V_2,V_3)$, $W=(W_1,W_2,W_3):\R^3\to\R^3$:

\begin{equation}\label{nhfhgtfredwM3S54E4y5C236terdhrCHE}
\curl(V\times W)
= \div W\;V- \div V \;W+(W\cdot\nabla)\, V- (V\cdot\nabla)\,W.\end{equation}

To prove this, we can write $V=V_1e_1+V_2e_2+V_3 e_3$ and notice that, since (<ref>) is linear in $V$, it suffices to prove it when $V$
reduces to each of the components $V_j e_j$.

Thus, up to reordering the coordinates, we can focus on the case in which $V$ is actually $V_1 e_1$.
Hence we calculate

\begin{eqnarray*}
&&\curl((V_1e_1)\times W)+ \div( V_1e_1) \;W-(W\cdot\nabla)\, V_1e_1+( (V_1e_1)\cdot\nabla)\,W
\\&=&
\curl(0,-V_1W_3,V_1W_2)
+ \partial_1 V_1\,W-\sum_{j=1}^3W_j\partial_j V_1e_1+V_1\partial_1 W\\
&=&\big( \partial_2(V_1W_2)+\partial_3(V_1W_3),\,-\partial_1(V_1W_2),\,-\partial_1(V_1W_3)\big)
+ \partial_1 V_1\,W-\sum_{j=1}^3W_j\partial_j V_1e_1+V_1\partial_1 W
\\&=&
\big(
\partial_2V_1\,W_2+V_1\partial_2 W_2+\partial_3 V_1\,W_3+V_1\partial_3 W_3+\partial_1 V_1\,W_1
0,\,0\big)-\sum_{j=1}^3W_j\partial_j V_1e_1\\&=&
\big(V_1\partial_2 W_2+V_1\partial_3 W_3
\\&=& \div W\;V_1e_1
\end{eqnarray*}

and this completes the proof of (<ref>).

The geometric argument used to prove (<ref>).

It is also useful to recall the vector calculus identity

\begin{equation}\label{YY08ojew9oik903hyto3298ythgfjb8ihyri-1urX} \div(\curl V)=0.\end{equation}

For this, we take any ball $B\subset\R^3$ and we consider the spherical surface $S:=\partial B$.
We split $S$ into two hemispheres $S^+$ and $S^-$
with a common equator $\eta$, which is the intersection between the boundaries of the surfaces $S^+$ and $S^-$.
Looking at $\eta$ as a curve, the natural direction of travel along $\eta$ when considered as the boundary of $S^+$
is opposite to the one obtained when considering it as the boundary of $S^-$, see Figure <ref>
(here, we are endowing $S^+$ and $S^-$ with the external unit vector field of $\partial B$).
As a result, the circulation of a vector field $V$ along $\eta$ considered as the boundary of $S^+$ (that we denote by ${\mathcal{C}}^+(V)$)
is opposite to the circulation of $V$ along $\eta$ considered as the boundary of $S^-$
(that we denote by ${\mathcal{C}}^-(V)$), that is

$$ {\mathcal{C}}^+(V)=-{\mathcal{C}}^-(V).$$

In addition, by Stokes' Theorem, we know that ${\mathcal{C}}^\pm(V)$ agrees with the flux
of $\curl V$ through the surface $S^+$.

Using these bits of information and the Divergence Theorem, we conclude that

\begin{eqnarray*}
\int_{S^+}\curl V\cdot\nu+\int_{S^-}\curl V\cdot\nu=\int_{\partial B}\curl V\cdot\nu=
\int_{B}\div(\curl V).
\end{eqnarray*}

From the arbitrariness of the ball $B$, we obtain (<ref>), as desired.

We also point out that, for any smooth function $\psi$,

\begin{equation}\label{08ojew9oik903hyto3298ythgfjb8ihyri-1urX} \curl\nabla\psi=0.\end{equation}

This can be checked by employing Stokes' Theorem. Indeed, if $\Sigma$ is any smooth oriented element of surface in $\R^3$
with boundary $\partial \Sigma $, the flux of $\curl\nabla\psi$ through the surface $\Sigma$
is equal to the circulation of $\nabla\psi$ over the loop $\partial\Sigma$. The latter object, if we describe $\partial\Sigma$
by a closed curve $\gamma:[0,L]\to\R^3$ with arc length parameterization, can be written as

$$ \int_0^L \nabla\psi(\gamma(\tau))\cdot\dot\gamma(\tau)\,d\tau=
\int_0^L \frac{d}{d\tau}\psi(\gamma(\tau))\,d\tau=\psi(\gamma(L))-\psi(\gamma(0))=0.$$

Since this is valid for an arbitrary surface $\Sigma$, the proof of (<ref>) is complete.

Now, to prove (<ref>), we let $\omega:=\curl v$, we use the momentum balance equation
and the barotropicity in (<ref>)
in combination with (<ref>) and we see that

\begin{equation}\label{dui39urteifjewgtlie4y9548796575903549873204327957608497685490549030594}
\begin{split}
0\,&=\curl\left( \partial_t v+ (v\cdot\nabla) v+\frac{\nabla p}{\rho}\right)\\&=\curl\left( \partial_t v+ \frac{1}{2}\nabla |v|^2-
v\times (\curl v)+\frac{\nabla p}{g(p)}\right)\\
&=\partial_t \omega+\curl\left(  \frac{1}{2}\nabla |v|^2-
v\times \omega+\frac{\nabla p}{g(p)}\right).
\end{split}\end{equation}

Now we observe that

\begin{eqnarray*} &&
\curl\left(\frac{\nabla p}{g(p)}\right)= \left(\partial_2 \left(\frac{\partial_3 p}{g(p)}\right)-
\partial_3 \left(\frac{\partial_2 p}{g(p)}\right),\partial_3 \left(\frac{\partial_1 p}{g(p)}\right)-
\partial_1 \left(\frac{\partial_3 p}{g(p)}\right),
\partial_1 \left(\frac{\partial_2 p}{g(p)}\right)-
\partial_2 \left(\frac{\partial_1 p}{g(p)}\right)
\right)\\&&\qquad=\left( \frac{\partial_{23} p}{g(p)}-\frac{\partial_3p\,g'(p)\,\partial_2p}{g^2(p)}
-\frac{\partial_{23} p}{g(p)} + \frac{\partial_2 p\,g'(p)\,\partial_3p}{g^2(p)},\right.\\ &&\qquad\qquad\qquad\left.
\frac{\partial_{13} p}{g(p)}-\frac{\partial_1 p\,g'(p)\,\partial_3p}{g^2(p)}
-\frac{\partial_{13} p}{g(p)}+\frac{\partial_3 p\,g'(p)\,\partial_1p}{g^2(p)},\right.\\ &&\qquad\qquad\qquad\left.
\frac{\partial_{12} p}{g(p)}-\frac{\partial_2 p\,g'(p)\,\partial_1p}{g^2(p)}
-\frac{\partial_{12} p}{g(p)}+\frac{\partial_1 p\,g'(p)\,\partial_2p}{g^2(p)}
\right)\\&&\qquad=(0,0,0).
\end{eqnarray*}

Hence, from this, (<ref>) and (<ref>),

\begin{equation*}
\begin{split}
0=\partial_t \omega-\curl\left( v\times \omega\right).
\end{split}\end{equation*}

This and (<ref>) lead to

\begin{equation*}
\partial_t \omega=
\curl(v\times \omega)
= \div \omega\;v- \div v \; \omega+( \omega\cdot\nabla)\, v- (v\cdot\nabla)\, \omega.\end{equation*}

Thus, combining the incompressibility condition in (<ref>)
and (<ref>),

\begin{equation*}
\partial_t \omega=( \omega\cdot\nabla)\, v- (v\cdot\nabla)\, \omega\end{equation*}

and therefore, following the fluid parcel $x(t)$,

\begin{equation*}\begin{split}&
\frac{d}{dt} \omega(x(t),t)=
\partial_t \omega(x(t),t)+(\dot x(t)\cdot\nabla)\omega(x(t),t)\\&\qquad
=\partial_t \omega(x(t),t)+\big(v( x(t),t)\cdot\nabla\big)\omega(x(t),t)
=\big( \omega(x(t),t)\cdot\nabla\big)\, v(x(t),t).
\end{split}\end{equation*}

As a result, if $Z(t)=(Z_1(t),Z_2(t),Z_3(t)):=\omega(x(t),t)$, we have that $Z$ is a solution of the Cauchy problem
for ordinary differential equations

$$ \begin{dcases}\displaystyle
\dot Z(t)=\sum_{k=1}^3 Z_k(t)\, \partial_k v(x(t),t),\\
\end{dcases}$$

By the uniqueness result for solutions of ordinary differential equations we thereby infer that $Z(t)=0$ for all $t\in\R$
and this completes[As a technical detail, we point out that we are freely assuming that if, for a given function $f$,
we know that $f(x(t),t)=0$ for all times $t$, then we can infer that $f(x,t)=0$
for all spatial positions $x$ and times $t$. This is because we are implicitly assuming that the flow
of the fluid exists for all times: hence, given any position $x$ we can “flow $x$ backwards for a time $t$”,
that is consider the fluid parcel evolution starting at some point $x_0$ at time $0$ arriving at $x$ at time $t$
(this gives us the possibility of choosing $x(t)=x$).]
the proof of (<ref>).

§.§ Propagation of sound waves

Now we retake the set of Euler fluid flow equations for the velocity $v$ of an inviscid fluid,
namely we look at the mass conservation equation in (<ref>), at
the momentum balance equation in (<ref>) and at a constitutive relation
linking the pressure $p$ and the density $\rho$, namely (neglecting the gravity effects)

\begin{equation}\label{EUMSD-OS-D-TUT}
\begin{dcases}&
\partial_t \rho+\div(\rho v)=0,\\&
\rho \partial_t v+\rho (v\cdot\nabla) v=-\nabla p,\\&
\end{dcases}
\end{equation}

for some function $f$. In practice, it is useful to take $f$ to be strictly increasing (the higher the density of the fluid, the higher the pressure produced; this would also give that the fluid is barotropic in the setting of Section <ref>).
We can think of (<ref>) as a very simple model for a gas
and we aim here at understanding the propagation of sound waves in such a medium.

For this objective, we first point out that a solution of (<ref>)
is provided by $(v,p,\rho)=(0,p_0,\rho_0)$, for every $\rho_0\in(0,+\infty)$ and $p_0:=f(\rho_0)$.
This configuration corresponds to the physical situation of a gas at rest, with constant density and pressure.

Suppose now that we perturb this configuration by creating a small variation of the gas pressure,
for instance by singing, or by playing[To play guitar, see the forthcoming Section <ref>.]
a guitar. The idea is thus to look for (at least approximate) solutions of the form

$$ v(x,t)=\e v_1(x,t),\qquad p(x,t)=p_0+\e p_1(x,t)\qquad{\mbox{and}}\qquad \rho(x,t)=\rho_0+\e\rho_1(x,t).$$