In this setting, $\e$ is a small parameter and our goal is to determine the functions $(v_1,p_1,\rho_1)$ in order to satisfy the set of equations in (<ref>), formally up to negligible errors in $\e$. To this end, we observe that

$$ p_0+\e p_1=p=f(\rho)=f(\rho_0+\e \rho_1)=f(\rho_0)+\e f'(\rho_0) \rho_1+o(\e)=p_0+\e f'(\rho_0) \rho_1+o(\e),$$

leading to the choice

\begin{equation*}
p_1=f'(\rho_0) \rho_1,
\end{equation*}

up to higher orders in $\e$ that we here sloppily disregard.

A plane wave.

\begin{equation*}
0=\partial_t \rho+\div(\rho v)=\e\partial_t \rho_1+\e\rho_0\div v_1+o(\e)
\end{equation*}

\begin{equation*}
0= \rho \partial_t v+\rho (v\cdot\nabla) v+\nabla p=\e\rho_0 \partial_t v_1+\e\nabla p_1+o(\e)=
\e\rho_0 \partial_t v_1+\e f'(\rho_0)\nabla \rho_1+o(\e).
\end{equation*}

These observations give

\begin{equation*}
\begin{dcases}
&\partial_t \rho_1+ \rho_0\div v_1=0,\\
&\rho_0 \partial_t v_1+f'(\rho_0)\nabla \rho_1=0.
\end{dcases}\end{equation*}

As a consequence,

\begin{equation*}
\partial_{tt} \rho_1=-\partial_t\big(\rho_0\div v_1\big)
=-\div \big(\rho_0\partial_t v_1\big)=
\div\big(
f'(\rho_0)\nabla \rho_1
\big)
=f'(\rho_0)\Delta \rho_1

That is, the perturbed density $\rho_1$
satisfies the equation

\begin{equation}\label{WAYEBJJD121t416JH098327uyrhdhc832nbcM}
\partial_{tt}u=c^2 \Delta u
\end{equation}

with $c:=\sqrt{f'(\rho_0)}>0$.
Equation (<ref>) is called[In terms of the classification
mentioned in footnote <ref>
on page CLASSIFICATIONFOOTN, equation (<ref>) is hyperbolic.
Indeed, we can take here $N=n+1$, $X=(x,t)$ and

$$ a_{ij}=\begin{dcases}
c^2 & {\mbox{ if }}i=j\in\{1,\dots,n\},\\
-1 & {\mbox{ if }}i=j=n+1,\\
0 & {\mbox{ otherwise,}}
\end{dcases}$$

thus producing $n$ strictly positive, and one strictly negative, eigenvalues.]
in jargon the “wave equation”wave equation.

To have a feeling of the “propagation of waves” encoded in equation (<ref>),
we can consider a smooth function $u_\star:\R\to\R$ and a direction $\omega\in\partial B_1$,
and thus define $u_\omega(x,t):=u_\star(\omega\cdot x-ct)$. Notice that $u_\omega$ is indeed a solution of (<ref>),
physically corresponding to a traveling plane waveplane wave.
Indeed, its evolution in time corresponds to a translation of $u_\star$ along the direction $\omega$
with speed $c$. Also, for $\kappa\in\R$, the parallel hyperplanes $\{\omega\cdot x-ct=\kappa\}$
(again, traveling with constant speed $c$ in direction $\omega$)
correspond to the level sets $\{ u_\omega=\kappa_\star\}$, where $\kappa_\star:=u_\star(\kappa)$
(the level sets of a wave are often called “wavefronts”
and they have special importance since they correspond to the surfaces on which,
at a given moment of time, all particles of the medium undergo the same motion).

See Figure <ref> for a series of frames of a traveling plane wave
(the picture is obtained by choosing $u_\star(r):=e^{-r^2}$, $\omega:=\left(\frac{\sqrt2}{2},\frac{\sqrt2}{2}\right)$, $c=1$
and instants of time $t\in\{0,1,2,3\}$).
See also Figure <ref>
for the pictures of the corresponding level sets of the wave
(i.e., of the corresponding wavefronts).

Wavefronts of the plane wave in Figure <ref>.

§.§ Hydrodynamics doesn't always work right

Now that we built some confidence in the equations of hydrodynamics,
it is time to question and challenge our own knowledge, to appreciate its power as well as its limitation.
After all, it's not so important to know, but rather to know of not knowing.

In particular, the great Euler (“Master of us all” according to Laplace) allowing,
we will recall here a rather striking, and perhaps surprising, flaw of the hydrodynamic equations presented in Section <ref> in the description of aerodynamic forces.

We will focus here on the notion of “drag”drag, which is the
force acting in the opposite direction
to an object moving in a fluid (say, the air resistance, for example).
This is not just a mathematical curiosity: for instance, according to Wikipedia
“induced drag tends to be the most important component for airplanes during take-off or landing flight”
and, more importantly,
“in the physics of sports, the drag force is necessary to explain the motion of balls, javelins, arrows and frisbees and the performance of runners and swimmers”.

Portrait of Jean Le Rond d'Alembert (Public Domain image from

The issue about the drag that we discuss here was discovered by Jean le Rond d'Alembert in 1752.
Often, this remarkable discovery is called with the name of
“d'Alembert's paradox”d'Alembert's paradox;
see Figure <ref>
for a portrait of d'Alembert's intense glance (by Maurice Quentin de La Tour). In a nutshell, we will prove that

\begin{equation}\label{PARADDA}
\begin{split}
&{\mbox{if a body is moving with constant velocity}}\\&{\mbox{in an incompressible, inviscid, irrotational and steady flow,}}\\
&{\mbox{the corresponding drag force is zero.}}
\end{split}
\end{equation}

Of course, we may wonder how such a statement can be coherent with our own
everyday experience, in which the effect of air resistance is usually apparent and quite decisive
(at least when we catch a flight or we kick a ball). Though the full understanding
of hydrodynamics is likely way above the present possibilities of science,
one common explanation of the discrepancy between the claim in (<ref>) and real life situations
relies on remarking that the fluid's description in (<ref>) is too idealistic.
In particular, the occurrence of the paradox is usually attributed
to the neglected effects of viscosity[The importance of viscosity actually reflects the importance of the Laplace operator in the model: compare equations (<ref>) and (<ref>).

Though we do not explicitly use this here, for completeness we point out that
an incompressible and  irrotational fluid is automatically inviscid
(hence, the lack of viscosity is also a byproduct of two assumptions, namely the incompressibility and
the lack of vortexes, which, together, entail very restrictive byproducts): indeed, using the vector calculus

$$ \nabla (\div V)-\curl(\curl V)=\Delta V,$$

(see e.g. (<ref>) for a proof) we have that for incompressible and irrotational fluid it holds that

$$ \Delta v=\nabla (\div v)-\curl(\curl v)=0-0=0,$$

hence the Navier-Stokes equation in (<ref>)
boils down to the Euler equation in (<ref>).]
(roughly speaking, one of the effects of viscosity <cit.> is to produce
thin boundary layersboundary layer near the surface of the moving object, which may entail
friction, flow separation and a low-pressure wake, leading to pressure drag).
However, the official resolution of the question is possibly debatable in its full generality:
for instance, in relation to the paradoxical features of (<ref>),
Garrett Birkhoff <cit.> states[Garrett Birkhoff's statements
were in turn criticized by James J. Stoker <cit.>
(who was concerned about the possibility that readers were
“very likely to get wrong ideas about some of the important
and useful achievements in hydrodynamics”
being misled by the “negative aspects of the theory”). Some of the original statements in <cit.>
were revised in the second edition <cit.>.

All in all, much more knowledge has to be built before we reach a complete understanding of
fluid mechanics and hydrodynamics.]
that: “I think that to attribute them all to the neglect of viscosity is an unwarranted oversimplification. The root lies deeper, in lack of precisely that deductive rigor whose importance is so commonly minimized by physicists and engineers”.

Now we work out the mathematical details needed to prove the claim in (<ref>).
To achieve this goal, we recall
the Euler equation in (<ref>)
and we complement it with the mass transport equation in (<ref>),
the incompressibility condition in (<ref>)
and an irrotational condition. We also assume that the flow is “steady”:
this word usually refers to a situation in
which the fluid properties at a point in the system do not change over time
(they only depend on the point); concretely, to avoid ambiguities, we will suppose here that the density $\rho$
is only a function of $x$ (no dependence on the time $t$).

In this setting, the velocity $v$
of an incompressible fluid, subject to neither viscosity nor gravity effects, with steady density $\rho$ and pressure $p$
and presenting no vortexes can be described by the set of equations (valid outside the body $\Omega$)

\begin{equation}\label{VCISTPOK-2}
\begin{dcases}
\rho \partial_t v+\rho (v\cdot\nabla) v=-\nabla p,\\
\div(\rho v)=0,\\
\div v=0,\\
\curl v=0.
\end{dcases}
\end{equation}

Now, combining (<ref>) with the irrotationality condition in (<ref>) we deduce that

\begin{equation}\label{VCISTPOK-3}
(v\cdot\nabla )v=\frac{1}{2}\nabla |v|^2.\end{equation}

We also consider the differential form $\omega:=v_1\,dx_1+v_2\,dx_2+v_3\,dx_3$ and we observe that

\begin{eqnarray*} d\omega&=&-\partial_2 v_1\,dx_1\wedge dx_2-
\partial_3 v_1\,dx_1\wedge dx_3+\partial_1 v_2\,dx_1\wedge dx_2\\&&\qquad
-\partial_3 v_2\,dx_2\wedge dx_3+\partial_1v_3\,dx_1\wedge dx_3
+\partial_2v_3\,dx_2\wedge dx_3\\&=&(\curl v)_3\,dx_1\wedge dx_2-(\curl v)_2\,dx_1\wedge dx_3+
(\curl v)_1\,dx_1\wedge dx_2\\&=&0,
\end{eqnarray*}

thanks again to the irrotationality condition in (<ref>).

Accordingly, the differential form $\omega$ is exact, and therefore closed, thanks to the
Poincaré Lemma (see e.g. <cit.>). This gives that there exists a velocity potential function $\varphi$
satisfying $\omega=d\varphi$, that is

\begin{equation}\label{jiaknfujsuE0ujkjHAHLoefvujs98ou38jhPD98jhdjhfX8}
v=\nabla \varphi. \end{equation}

We stress that $\varphi$ is harmonic outside $\Omega$, since $0=\div v=\div\nabla\varphi=\Delta\varphi$.
Moreover, the first equation in (<ref>),
and (<ref>) entail that

\begin{equation}\label{kSUDSpREdfskljhREHGkozWI23f3aziw}\begin{split}&
0= \partial_t v+(v\cdot\nabla) v+\frac{\nabla p}\rho=
\partial_t v+\frac{1}{2}\nabla |v|^2+\frac{\nabla p}\rho\\&\qquad\quad=
\partial_t \nabla \varphi+\frac{1}{2}\nabla |\nabla\varphi|^2+\frac{\nabla p}\rho=
\nabla\left( \partial_t \varphi+\frac{1}{2}|\nabla\varphi|^2\right)+\frac{\nabla p}{\rho}.\end{split}
\end{equation}

Now we recall the assumption in (<ref>) that the object is moving with constant speed, say $v_0\in\R^3\setminus\{0\}$.
That is, if we denote by $x(t)$ the position at a given time $t$ of a parcel of the fluid, we have that $\dot x(t)=v_0$
and accordingly $v(x(t),t)=v_0$ for every time $t$.

From this observation we infer that

$$ 0=\frac{d}{dt}v(x(t),t)=\partial_t v(x(t),t)+(v_0\cdot\nabla) v(x(t),t).$$

and (<ref>) (recall also
footnote <ref> on page LAFOPRI)
give that

$$0= \partial_t \nabla\varphi+(v_0\cdot\nabla) (\nabla\varphi)=
\nabla\Big(\partial_t \varphi+(v_0\cdot\nabla)\varphi\Big)

Consequently, the function $\partial_t \varphi+(v_0\cdot\nabla)\varphi$ is constant in space
and therefore there exists a function of time, say $R=R(t)$, such that

\begin{equation*} \partial_t \varphi+(v_0\cdot\nabla)\varphi=R.\end{equation*}

From this and (<ref>) we arrive at

\begin{equation}\label{kSUDSpREdfskljhREHGkozWI23f3aziw2}\frac{\nabla p}{\rho}=
-\nabla\left( R-(v_0\cdot\nabla)\varphi+\frac{1}{2}|\nabla\varphi|^2\right)=
\nabla\left( (v_0\cdot\nabla)\varphi-\frac{1}{2}|\nabla\varphi|^2\right).
\end{equation}

Furthermore, using the conservation of mass and the incompressibility condition in (<ref>),

\begin{equation}\label{kSUDSpREdfskljhREHGkozWI23f3aziw7}
0=\div(\rho v)=\nabla\rho\cdot v+ \rho\div v=
\nabla\rho\cdot v=\nabla\rho\cdot v_0.
\end{equation}

With this, we can now compute the drag that the air exerts on the moving body.
For this, we recall that the force $F$ of the air is produced by the pressure, acting normally on the surface
of the body (see (<ref>)). Hence, calling the moving object $\Omega$ and $\nu$ its outer normal, we have that

\begin{equation}\label{FEUYT90hfbvjjxJ} F=-\int_{\partial\Omega} p\nu.\end{equation}

The minus sign here above is just recalling that positive pressures are assumed to go in the opposite direction
with respect to the exterior normal.
Since the drag $D$ is the force acting in the direction of motion of the body, using the Divergence Theorem we thus conclude that

\begin{equation}\label{EXDIVTHMQU} D=F\cdot\frac{v_0}{|v_0|}
=-\int_{\partial\Omega} p\frac{v_0\cdot\nu}{|v_0|}=\int_{\R^3\setminus\Omega}
\div\left(p\frac{v_0}{|v_0|}\right)=\int_{\R^3\setminus\Omega} \frac{\nabla p\cdot v_0}{|v_0|}

Notice that to apply the Divergence Theorem here to the infinite region $\R^3\setminus\Omega$
(that corresponds to the region occupied by the fluid) we are implicitly supposing that the pressure decays sufficiently fast at infinity
(that is, the pressure disturbance is essentially localized in the vicinity of the moving object, which is a reasonable[Once again, in this chapter, if an argument sounds convincing we buy it! We will be much more skeptical about heuristic reasonings from next chapter on. For instance, when we need to integrate on exterior domains in Section <ref>, the reader will appreciate how careful we will be in checking explicitly the appropriate decay properties of the functions involved (which is perhaps a delicate and annoying, albeit absolutely necessary, detail to be taken care of to avoid nonsensical computations).

In any case, let us mention that the possibility of utilizing the Divergence Theorem in an exterior domain
in (<ref>) is quite a delicate matter (see e.g. the thoughtful argument in <cit.>)
and it is a byproduct of the decay at infinity of the fluid velocity and pressure which
is modeled on that of harmonic functions.

Arguably, in the setting of these notes, alternative arguments to the ones in <cit.>
could be provided by using
the Green's Representation Formula in (<ref>)
(say, in the domain $B_R\setminus\overline{\Omega}$,
sending $R\to+\infty$): for this one would need to assume that the velocity potential
has a limit at infinity (thus, from Cauchy's Estimates
in Theorem <ref>, one can also bound the derivative of the potential along $\partial B_R$).

In any case, the decay estimates of the velocity fields rely
on the fact that the fluid cannot penetrate inside the object (as formalized in (<ref>)),
since this information allows one to get rid of one order of magnitude in the corresponding estimates. This is an interesting
physical feature since it reveals the “global” character of fluid dynamics, in which some features
in the proximity of the moving body have a significant influence on the velocity field at infinity, and vice versa.

In general, the decay analysis for solutions of fluid mechanics equations
is a delicate matter, see e.g. <cit.> and the references therein.]
physical assumption, though quantifying it precisely and checking it rigorously would require some technical skills).

Hence, recalling (<ref>),

\begin{equation*}
\begin{split} D\,&= \int_{\R^3\setminus\Omega} \frac{\rho v_0}{|v_0|}\cdot\nabla\left( (v_0\cdot\nabla)\varphi-\frac{1}{2}|\nabla\varphi|^2\right)\\
&=\int_{\R^3\setminus\Omega} \left\{\div\left[ \left( (v_0\cdot\nabla)\varphi-\frac{1}{2}|\nabla\varphi|^2\right)\frac{\rho v_0}{|v_0|}\right]-
\frac{\nabla\rho \cdot v_0}{|v_0|}\left( (v_0\cdot\nabla)\varphi-\frac{1}{2}|\nabla\varphi|^2\right)\right\}.
\end{split}\end{equation*}

Combining this and (<ref>),

\begin{equation}\label{inHTKMSue21567536uhtrf}
D=\int_{\R^3\setminus\Omega}\div\left[ \left( (v_0\cdot\nabla)\varphi-\frac{1}{2}|\nabla\varphi|^2\right)\frac{\rho v_0}{|v_0|}\right].
\end{equation}

We now stress that since the air cannot penetrate inside the object, necessarily

\begin{equation}\label{PSLmoLMSxHNS-AR}
{\mbox{the normal component
of the fluid velocity vanish along~$\partial\Omega$,}}\end{equation}
namely $v_0\cdot\nu=0$
on $\partial\Omega$.

As a result, using (<ref>) and the Divergence Theorem again,

\begin{eqnarray*}&&D=
\int_{\R^3\setminus\Omega} \div\left[ \left( (v_0\cdot\nabla)\varphi-\frac{1}{2}|\nabla\varphi|^2\right)\frac{\rho v_0}{|v_0|}\right]
=-\int_{\partial\Omega} \left( (v_0\cdot\nabla)\varphi-\frac{1}{2}|\nabla\varphi|^2\right)\frac{\rho v_0}{|v_0|}\cdot\nu=0.
\end{eqnarray*}

This completes the proof of d'Alembert's paradox in (<ref>).

See e.g. <cit.> and the references therein for additional information
on d'Alembert's paradox.

Portrait of Nikolai Yegorovich Joukowski,
Museum of Moscow Aviation Institute (image by Just from
Wikipedia, licensed under the Creative Commons Attribution-Share Alike 4.0 International license).

§.§ Lift of an airfoil

To safely recover from the shock received in Section <ref>
by d'Alembert's paradox, we give now a positive result
in terms of the aerodynamic lift of an airfoilairfoil. This result was established independently
by Martin Kutta and Nikolai Yegorovich Joukowski in the early twentieth century and, in its simplest formulation,
can be stated[Though the theory that
we present in these pages is way too crude to
comprise aerodynamics in its full complexity, the importance of a statement like (<ref>) in the theory of
flight appears to be paramount, as confirmed by the solid recognition bestowed to the scientists involved. For instance,
the Russian Air Force Academy and one of the airports in Moscow
were named after Joukowski. See also Figure <ref>.

In terms of real world applications, a common belief is indeed that the setting in (<ref>)
is too idealized to detect the intricate patterns produced by real fluids in the vicinity of a traveling object: however,
it is also believed that the circulation detected in (<ref>) does reflect significant physical
information since in many concrete situations
the flow around a thin airfoil is composed of a narrow viscous region near the body
(a sort of “boundary layer”) outside which the idealized
and inviscid description of the flow in (<ref>) turns out to be sufficiently realistic.

The gist is then to apply the setting in (<ref>) not quite to the traveling body but to the aggregate of the body and its own boundary layer:
in particular, the loop to compute the circulation in (<ref>) must be chosen outside this boundary layer.

With this, the boundary layer gives any traveling object an “effective shape” that may be different from its physical shape by accounting for the region in which the velocity changes from zero at the surface to the stream value away from the surface: the effectiveness of (<ref>) is thus to incorporate, as much as possible, the effects of turbulence and skin friction into the effective shape described by this boundary layer, remaining with a more mathematically treatable description away from it.

In this sense, the success of
the Kutta-Joukowski theory Kutta-Joukowski theory
consists in providing a simple, but not trivial, approach to aerodynamics which incorporates some aspects of viscous effects, while neglecting others.]
as follows:

\begin{equation}\label{JTTA}
\begin{split}&
{\mbox{an airfoil in relative motion with constant velocity~$-v_0$}}\\&
{\mbox{to an ambient inviscid homogeneous, irrotational fluid}}\\&{\mbox{has a lift force (that is the component of the force perpendicular to~$v_0$)}}\\&{\mbox{of magnitude }}
\rho \,|v_0| \Gamma, \\&{\mbox{where~$\Gamma$ is a circulation of~$v_0$ along the cross section of the airfoil.}}\end{split}\end{equation}

It is interesting to observe that (<ref>) explains for instance the generation of a lift on a wing as a result
of the contribution of the circulation $\Gamma$ of the velocity field around the wing.
Remarkably, the fact that (<ref>) takes into account the component of the force perpendicular to $v_0$
clearly states the importance of such a result for the theory of flight (that is, this lift force is what
allows, in principle
at least, an airplane to take off, provided they manage to create
a sufficiently large circulation $\Gamma$).

To model an airfoil, that is a “long wing” (actually, infinitely long for simplicity),
we consider a (nice, contractible) planar domain $\Omega\subset\R^2$
and the wing given by $\Omega_\star:=\Omega\times\R$, see Figure <ref>.
We reconsider the Euler equation in (<ref>) (in the absence of gravity) and the incompressibility condition in (<ref>) outside $\Omega_\star$
assuming that, by symmetry, the fluid parameters $v$ and $p$
are actually independent of $x_3$ and also on time (leading to a completely
steady solution of the problem).
In this way, we may consider
the following set of equations for $v=v(x,y)$ and $p=p(x,y)$, with $X=(x,y)\in\R^2\setminus\Omega$,
for a given constant $\rho\in(0,+\infty)$
(the constancy of the density being the mathematical translation
that the fluid is homogeneous):

\begin{equation}\label{DFRE}
\begin{dcases}&\rho (v\cdot\nabla) v=-\nabla p,\\
&\div v=0,\\&\curl v=0.\end{dcases}\end{equation}

Two-dimensional model of an airfoil, as inspired by a sailplane.

Notice that in this setting (thinking the vector field $(v_1(x,y),v_2(x,y),0)$ to be three-dimensional
but with a trivial last entry to compute the curl), we have that

\begin{equation}\label{CUR25t34yuD2}
0=\curl v=\big(\partial_xv_2-\partial_y v_1\big) e_3.
\end{equation}

Moreover, since, according to the principle of inertia
discussed on page PRINERY,
we are taking an inertial frame in which the body appears to be still,

\begin{equation}\label{ALLIVO0}
\lim_{|(x,y)|\to+\infty}v(x,y)=v_0.\end{equation}

Also, the Euler equation in (<ref>) and the assumption that the density is constant give that

\begin{eqnarray*}&& v\cdot\nabla\left(\frac{|v|^2}2+\frac{p}\rho\right)=
v\cdot\left(\sum_{k=1}^2 v_k\nabla v_k+\frac{\nabla p}\rho\right)=
v\cdot\left(\sum_{k=1}^2 v_k\nabla v_k-(v\cdot\nabla) v\right)
\\&&\qquad
=\sum_{j,k=1}^2 v_j v_k\partial_j v_k-
v\cdot\left(\sum_{j=1}^2 v_j\partial_j v\right)
=\sum_{j,k=1}^2 v_j v_k\partial_j v_k-
\sum_{j,k=1}^2 v_k v_j\partial_j v_k=0.
\end{eqnarray*}

As a consequence, if $\beta(X):= \frac{|v(X)|^2}2+\frac{p(X)}\rho$,

$$ v\cdot \nabla\beta=0,$$

\begin{equation*}\frac{d}{dt}\beta(X(t))=\nabla\beta(X(t))\cdot \dot X(t)=
\nabla\beta(X(t))\cdot v( X(t))=0.
\end{equation*}

As a result, $\beta(X(t))=\beta(X(0))$ for every time $t$, that is

\begin{equation}\label{BERNEQl}
\frac{|v(X(t))|^2}{2}+\frac{p(X(t))}{\rho}=
\frac{|v(X(0))|^2}{2}+\frac{p(X(0))}{\rho}.
\end{equation}

Equations such as (<ref>) are often referred to with the name of Bernoulli's Principle.

A topological argument used to take care of the fact that $\R^2\setminus\Omega$
has a hole.

Now we rely on the two-dimensional structure of the problem combined with the divergence free condition
in (<ref>) to see that the differential form $\zeta:= v_2\,dx-v_1\,dy$ is closed,
$$ d\zeta=-(\partial_y v_2+\partial_x v_1)\,dx\wedge dy=-\div v\,dx\wedge dy=0.$$

Now, the exterior of $\Omega$ is not simply connected in $\R^2$, however we can consider two
regions ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$ which are simply connected and their union is the exterior of $\Omega$, see Figure <ref>.
In this setting, we can apply the Poincaré Lemma (see e.g. <cit.>)
to each of the regions ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$, finding two potentials $\psi_j:{\mathcal{R}}_j\to\R$ such that $d\psi_j=\zeta$ in ${\mathcal{R}}_j$, with $j\in\{1,2\}$.
In particular,

\begin{equation}\label{56STRE78} \partial_x\psi_j=v_2\qquad{\mbox{and}}\qquad\partial_y\psi_j=-v_1.\end{equation}

We observe that ${\mathcal{R}}_1\cap{\mathcal{R}}_2$ consists of two connected components which reach infinity,
say ${\mathcal{S}}_1$ and ${\mathcal{S}}_2$. We observe that

\[ \nabla(\psi_1-\psi_2)=
\nabla\psi_1-\nabla\psi_2=(v_2,-v_1)-(v_2,-v_1)=0,\]

in ${\mathcal{S}}_1$, hence there exists $c_1\in\R$ such that $\psi_1-\psi_2=c_1$ in ${\mathcal{S}}_1$.

Similarly, there exists $c_2\in\R$ such that $\psi_1-\psi_2=c_2$ in ${\mathcal{S}}_2$.
Up to incorporating the constant into $\psi_2$ (which leaves (<ref>) invariant), we can suppose that $c_2=0$.

We claim that

\begin{equation}\label{Cik231ma3223y0-2tkhojnhoiuwhg8732gtpiuywfng2uahHBSne}
\end{equation}

For this, let $A\in(\partial\Omega)\cap\overline{{\mathcal{S}}_1}$ and $B\in(\partial\Omega)\cap\overline{{\mathcal{S}}_2}$
and let $\vartheta:[0,\ell]\to(\partial\Omega)\cap\overline{{\mathcal{R}}_1}$ be a curve, parameterized by its arc length, lying in $(\partial\Omega)\cap\overline{{\mathcal{R}}_1}$ and
joining $A$ to $B$, say with $\vartheta(0)=A$ and $\vartheta(\ell)=B$.
It follows from the impenetrability condition in (<ref>)
that the velocity field is tangential to $\partial\Omega$ and consequently

\begin{equation}\label{DEFSTTHiow9805e04938} v(\vartheta(\tau))
=\alpha(\tau)\,\dot\vartheta(\tau)\qquad{\mbox{ for all }}\tau\in[0,\ell],\end{equation}

for some scalar function $\alpha(\tau)$.

Thus, since $\vartheta$ lies in $\overline{{\mathcal{R}}_1}$, we can utilize (<ref>) with $j:=1$
and write that

$$ \Big( -\partial_y\psi_1\big(\vartheta(\tau)\big),\,\partial_x\psi_1\big(\vartheta(\tau)\big)\Big)=
\alpha(\tau)\,\dot\vartheta(\tau)$$

and therefore

$$ \nabla\psi_1\big(\vartheta(\tau)\big)=
\alpha(\tau)\,\big(\dot\vartheta_2(\tau),\,-\dot\vartheta_1(\tau)\big).$$

From this, we obtain that

$$ \frac{d}{d\tau}\psi_1\big(\vartheta(\tau)\big)=
\nabla\psi_1\big(\vartheta(\tau)\big)\cdot\dot\vartheta(\tau)
\alpha(\tau)\,\big(\dot\vartheta_2(\tau),\,-\dot\vartheta_1(\tau)\big)\cdot
\big(\dot\vartheta_1(\tau),\,\dot\vartheta_2(\tau)\big)=0

and consequently

$$ \psi_1(B)=\psi_1\big(\vartheta(\ell)\big)=\psi_1\big(\vartheta(0)\big)=\psi_1(A).$$

Similarly, by considering a curve lying in $(\partial\Omega)\cap\overline{{\mathcal{R}}_2}$, one obtains that

$$ \psi_2(B)=\psi_2(A).$$

As a result,

\begin{eqnarray*}
\end{eqnarray*}

and the proof of (<ref>) is complete.

Hence, we can now define, in the exterior of $\Omega$,

$$ \psi:=
\begin{dcases}
\psi_1 & {\mbox{ in }}{\mathcal{R}}_1,\\
\psi_2 & {\mbox{ in }}{\mathcal{R}}_2\\
\end{dcases}$$

and we stress that this is a fair definition, since $\psi_1=\psi_2$ in ${\mathcal{R}}_1\cap{\mathcal{R}}_2$,
owing to (<ref>).

Furthermore, by (<ref>), in the exterior of $\Omega$,

\begin{equation}\label{9uojhexterimega}
\partial_x\psi =v_2\qquad{\mbox{and}}\qquad\partial_y\psi =-v_1.\end{equation}

Note in particular that $|\nabla\psi|=v$. Moreover, $\psi$ is harmonic, thanks to the irrotationality
situation pointed out in (<ref>).

The function $\psi$ is often referred to with the name of
“stream function”.
This name comes from a physical intuition. Indeed, the streamlines, i.e. the trajectories of fluid particles, are described by the level sets of the stream function, because

\begin{equation}\label{STERTEAMSDN} \frac{d}{dt}\psi(X(t))=\nabla\psi(X(t))\cdot \dot X(t)=\nabla\psi(X(t))\cdot v(X(t))=0,\end{equation}

thanks to (<ref>).

Let us now reassess Bernoulli's Principle in (<ref>).
In our setting we need a refinement of (<ref>) along the boundary of our moving object, namely we claim that, for every $x\in\partial\Omega$,

\begin{equation}\label{BERNEQ}
\end{equation}

for some constant $p_0\in\R$.

To prove this, we take the simplifying assumption that the velocity field possesses at most finitely many zeros along $\partial\Omega$. We denote by $Z_1,\dots,Z_N$ these zeros (if any).
We pick $X_0\in(\partial\Omega)\setminus\{Z_1,\dots,Z_N\}$. Up to a rigid motion, we assume that $X_0=0$
and that the tangent vector of $\partial\Omega$ at $0$ is horizontal.
Say, we describe $\Omega$ near $0$ as the subgraph of a smooth function $f:\R\to\R$ with $f(0)=0$ and $f'(0)=0$.
We can also replace $\psi$ by $\psi-\psi(0)$ and thus assume additionally that $\psi(0)=0$.
In this setting, from the impenetrability condition in (<ref>), we have that the fluid vector field at the origin
is horizontal and thus, by (<ref>), we have that $\nabla\psi(0)$ is vertical and different from zero.
Let us suppose that

\begin{equation}\label{deltastreamline}
\partial_y\psi(0)>0,\end{equation}

the case $\partial_y\psi(0)<0$ being similar, just considering level sets
of $\psi$ corresponding to negative, instead of positive, values.

In this framework, near $0$, the level sets of $\psi$ can be written as graph of smooth functions: more precisely,
near the origin, for small $\delta>0$, the level set $\{\psi=\delta\}$ can be identified with the graph $\{y=f_\delta(x)\}$
for a suitable $f_\delta$. It follows from (<ref>) that $f_{\delta'}\ge f_\delta$ if $\delta'\ge\delta$.
Also, since, by (<ref>), $\{\psi=\delta\}$ describes the streamlines of the fluid,
we have that $\{\psi=\delta\}$ is contained in the complement of $\Omega$ for all $\delta>0$ and consequently,
near the origin, $\partial\Omega$ lies below the graph of $f_\delta$ for all $\delta>0$.

Thus, given the monotonicity of $f_\delta$, near the origin we can define

$$ f_0(x):=\lim_{\delta\searrow0}f_\delta(x)$$

and we have that $f\le f_0$. Moreover, since

$$ \psi(x,f_0(x))=\lim_{\delta\searrow0}\psi(x,f_\delta(x))=\lim_{\delta\searrow0}\delta=0=\psi(x,f(x)),$$

we deduce from (<ref>) that $f_0=f$, see Figure <ref>.

Foliation of streamlines in the vicinity of an airfoil.

This shows that the streamline emanating from the origin (that is, from every point of $\partial\Omega$ with nonzero velocity field)
remains on $\partial\Omega$.

Thus, if the velocity field does not vanish on $\partial\Omega$, then $\partial\Omega$ consists of a streamline
and then (<ref>) follows from (<ref>).

Partitioning $\partial\Omega$ as in (<ref>)
(not necessarily a realistic picture).

If instead
the velocity field vanishes at $Z_1,\dots,Z_N\in\partial\Omega$, we can decompose $\partial\Omega$
into open arcs $\Lambda_1,\dots,\Lambda_N$ such that

\begin{equation}\label{PTOP12DAILATENSIWASTbni9buj9unooREDItangeFIlAHOLEWFINASDJGIOSE2}\partial\Omega=\Lambda_1\cup\dots\cup\Lambda_N\cup\{Z_1,\dots,Z_N\},\end{equation}

in such a way that the velocity field does not vanish on each $\Lambda_j$ and $\Lambda_j$
is a streamline connecting $Z_j$ to $Z_{j+1}$ in infinite time (with $Z_{N+1}:=Z_1$), see Figure <ref>.
In this setting, we can consider a trajectory $X(t)$ starting at a given point of $\Lambda_j$ and deduce from (<ref>) that,
for every $x\in\Lambda_j$,

$$ |v(x)|^2+\frac{2p(x)}{\rho}=p_j,$$

for some $p_j\in\R$, and

$$ p_j=\lim_{t\to\pm\infty}|v(X(t))|^2+\frac{2p(X(t))}{\rho}$$

which leads to

$$ p_j=|v(Z_j)|^2+\frac{2p(Z_j)}{\rho}=|v(Z_{j+1})|^2+\frac{2p(Z_{j+1})}{\rho}.$$

This actually says that $p_1=\dots=p_N$ and thus (<ref>) plainly follows, as desired.

Now we describe $\partial\Omega$ as a smooth and closed curve, traveled counterclockwise and parameterized by its arc length, say $\gamma:[0,L]\to\C$, for some $L>0$ (in particular, notice that $\gamma(0)=\gamma(L)$).
We  identify points of $\R^2$ with complex numbers.
In this way, given $\tau\in[0,L]$, a unit tangent vector to $\partial\Omega$ at $\gamma(\tau)$
is given by $\dot\gamma(\tau)$ and the unit normal vector $\nu(\gamma(\tau))$ pointing outwards is given by a clockwise rotation
of $\dot\gamma(\tau)$ by an angle $\frac\pi2$, corresponding to the complex multiplication by $e^{-i\pi/2}=-i$
(that is, the outward unit normal can be written as $-i\dot\gamma(\tau)$).

For this reason, in light of (<ref>), we have that
the force acting on the airfoil (after identifying $\R^2$ with $\C$) takes the form

\begin{equation*} F=-\int_{\partial\Omega} p\nu
=-\int_0^L p(\gamma(\tau))\,\nu(\gamma(\tau))\,d\tau=i\int_0^L p(\gamma(\tau))\,\dot\gamma(\tau)\,d\tau.
\end{equation*}

From this and (<ref>) we arrive at

\begin{equation}\label{BERNEQ2} \frac{2F}\rho
=i\int_0^L \frac{2p(\gamma(\tau))}\rho\,\dot\gamma(\tau)\,d\tau
\left(p_0-|\nabla\psi(\gamma(\tau))|^2\right)
\,\dot\gamma(\tau)\,d\tau.
\end{equation}

$$ \int_0^L \dot\gamma(\tau)\,d\tau=\gamma(L)-\gamma(0)=0,$$

we deduce from (<ref>) that

\begin{equation}\label{BERNEQ3} \frac{2F}\rho
=-i\int_0^L |\nabla\psi(\gamma(\tau))|^2\,\dot\gamma(\tau)\,d\tau.
\end{equation}

We also recall that, in light of the impenetrability condition in (<ref>)
combined with (<ref>),
we have that $\nabla\psi$ is orthogonal to
the unit tangent vector $\dot\gamma$.
Thus, we write that

\begin{equation}\label{BERNEQ4} \nabla\psi(\gamma(\tau))=\pm i

Also, in the complex variable setting, that is using the notation $z=x+iy$ and identifying $x+iy$ with $(x,y)$, we set

\begin{equation}\label{ENGBSC2O232M4P3AND5D} w(z)=w(x+iy):=\partial_x\psi(x,y)-i\partial_y\psi(x,y),\end{equation}

and it follows that $w$ is the complex conjugated of $\nabla\psi$ and $|w|=|\nabla\psi|$.

Hence, taking the complex conjugation in (<ref>),

\begin{equation*} w(\gamma(\tau))=
\pm \overline{i|w(\gamma(\tau))|\,\dot\gamma(\tau)}=
\mp i|w(\gamma(\tau))|\,\overline{\dot\gamma(\tau)}.\end{equation*}

Taking the square of this identity

\begin{equation*} w^2(\gamma(\tau))=-|w(\gamma(\tau))|^2\,(\overline{\dot\gamma(\tau)})^2
\end{equation*}

and as a result

\begin{equation*} \begin{split}&w^2(\gamma(\tau)) \dot\gamma(\tau)=-|w(\gamma(\tau))|^2\,(\overline{\dot\gamma(\tau)})^2{\dot\gamma(\tau)}
\end{equation*}

We thus combine this information with the complex conjugation of (<ref>) and we conclude that

\begin{equation}\label{wdotgammatau}
\begin{split}&\overline{\frac{2i F}\rho}
=\overline{\int_0^L |\nabla\psi(\gamma(\tau))|^2\,\dot\gamma(\tau)\,d\tau}=
\int_0^L |\nabla\psi(\gamma(\tau))|^2\,\overline{\dot\gamma(\tau)}\,d\tau\\&\quad\qquad\qquad\qquad=-\int_0^L w^2(\gamma(\tau)) \dot\gamma(\tau)\,d\tau=-
\oint_\gamma w^2(z)\,dz.\end{split}
\end{equation}

Also, using the harmonicity of $\psi$, we observe that

$$ \partial_x(\Re w)=\partial_{xx}\psi=-\partial_{yy}\psi=\partial_y(\Im w)$$

and moreover

$$ \partial_y(\Re w)=\partial_{xy}\psi(x,y)=-\partial_x(\Im w).$$

These observations give that $w$ satisfy the Cauchy-Riemann equations
and therefore $w$ is holomorphic outside $\Omega$.

We can therefore take $R>0$ large enough such that $\Omega\Subset B_R$
and employ Laurent's Theorem (see e.g. <cit.>, applied here with $R_1:=R$, $R_2:=+\infty$
and $z_0:=0$): in this way we
obtain that, for all $z$ in the exterior of $B_R$,
the following series representation holds
true for suitable $a_k$, $b_k\in \C$:

$$ w(z)=w_0(z)+w_\star(z)\qquad{\mbox{with}}\qquad
w_0(z):=\sum_{k=0}^{+\infty} a_k z^k\qquad{\mbox{and}}\qquad
w_\star(z):=\sum_{k=1}^{+\infty}\frac{b_k}{ z^k},$$

where the power series representing $w_0$ converges in $\C$ and the one representing $w_\star$
converges in the exterior of $B_R$.

In particular, $w_0$ is holomorphic in the whole of $\C$
and the convergence of the power series defining $w_\star$ entails (see e.g. <cit.>) that

$$ \limsup_{k\to+\infty} |b_k|^{1/k}\le{2R}.$$

As a result, we can take $k_0\in\N\cap[2,+\infty)$ sufficiently large such that $|b_k|\le(3R)^k$ for every $k> k_0$
and therefore

\begin{eqnarray*}&& \limsup_{|z|\to+\infty}|w_\star(z)|\le
\limsup_{|z|\to+\infty}\sum_{k=1}^{k_0}\frac{|b_k|}{| z|^k}+\limsup_{|z|\to+\infty}\sum_{k=k_0+1}^{+\infty}\frac{|b_k|}{| z|^k}
\le0+\limsup_{|z|\to+\infty}\sum_{k=k_0+1}^{+\infty}\left(\frac{3R}{| z|}\right)^k\\&&\qquad\qquad\qquad\qquad
\left(\frac{3R}{| z|}\right)^{k_0+1}\frac{1}{1-\displaystyle\frac{3R}{| z|}}=0.
\end{eqnarray*}

For this reason, recalling (<ref>),

$$ |v_0|=\limsup_{|(x,y)|\to+\infty}|\nabla\psi	(x,y)|=\limsup_{|z|\to+\infty}|w(z)|=
\limsup_{|z|\to+\infty}|w_0(z)+w_\star(z)|=\limsup_{|z|\to+\infty}|w_0(z)|.$$

As a consequence, $w_0$ is bounded in the whole of $\C$ and therefore, by Liouville's Theorem
(see e.g. <cit.>) we have that $w_0$ is constant.

That is, by (<ref>), $w_0(z)=i \overline{v_0}$ for all $z\in\C$, whence, for all $z$ in the exterior of $B_R$,

\begin{equation}\label{EXTE6B-pjrmfR} w(z)=i \overline{v_0}+w_\star(z)=
i \overline{v_0}+\sum_{k=1}^{+\infty}\frac{b_k}{ z^k}.\end{equation}

\begin{eqnarray*}&&\limsup_{M\to+\infty}\left|\oint_{\partial B_M}\sum_{k=k_0+1}^{+\infty}\frac{b_k}{ z^k}\,dz\right|\le
\limsup_{M\to+\infty}\sum_{k=k_0+1}^{+\infty}\frac{2\pi|b_k|}{M^{k-2}}\\&&\quad
\le\limsup_{M\to+\infty}\sum_{k=k_0+1}^{+\infty}\frac{2\pi M^2(3R)^k}{M^{k}}
2\pi M^2\left(\frac{3R}{M}\right)^{k_0+1}\frac{1}{1-\displaystyle\frac{3R}{M}}=0.\end{eqnarray*}

Thus, by (<ref>) and Cauchy's Theorem (see e.g. <cit.>),

\begin{equation}\label{SunmgeDediamJofeMSy0r987wqugf7wegbDEDS}\begin{split}&
\oint_\gamma w(z)\,dz=
\lim_{M\to+\infty}\oint_{\partial B_M} w(z)\,dz\\&\quad=
\lim_{M\to+\infty}\oint_{\partial B_M} \left(i \overline{v_0}+\sum_{k=1}^{+\infty}\frac{b_k}{ z^k}\right)\,dz=
\lim_{M\to+\infty}\oint_{\partial B_M} \left(i \overline{v_0}+\sum_{k=1}^{k_0}\frac{b_k}{ z^k}\right)\,dz.\end{split}
\end{equation}

We stress that the function in the latter integrand is holomorphic in $B_M$ and we can thereby employ
Cauchy's Residue Theorem (see e.g. <cit.>), thus  deducing from (<ref>) that

\begin{equation}\label{SunmgeDediamJofeMSy0r987wqugf7wegbDEDS-2}
\oint_\gamma w(z)\,dz=
\lim_{M\to+\infty}\oint_{\partial B_M} w(z)\,dz=2\pi i b_1.
\end{equation}

In addition,

\begin{equation*}
\begin{split}&
\oint_\gamma w(z)\,dz=\int_0^L
\Big(\partial_x\psi(\gamma(\tau))-i\partial_y\psi(\gamma(\tau))\Big)
\Big( \dot\gamma_1(\tau)+i\dot\gamma_2(\tau)\Big)\,d\tau\\&\quad=
\int_0^L
\Big(\partial_x\psi(\gamma(\tau)) \dot\gamma_1(\tau) +\partial_y\psi(\gamma(\tau))\dot\gamma_2(\tau)-i\partial_y\psi(\gamma(\tau))\dot\gamma_1(\tau)
\Big)\,d\tau\\&\quad=
\int_0^L
\left(  \frac{d}{d\tau}\psi(\gamma(\tau)) -i\partial_y\psi(\gamma(\tau))\dot\gamma_1(\tau)
\right)\,d\tau\\&\quad=\psi(\gamma(L))-\psi(\gamma(0))-i
\int_0^L
\left(  \partial_y\psi(\gamma(\tau))\dot\gamma_1(\tau)-\partial_x\psi(\gamma(\tau))\dot\gamma_2(\tau)
\right)\,d\tau\\&\quad=
\int_0^L
\left(  \partial_y\psi(\gamma(\tau))\dot\gamma_1(\tau)-\partial_x\psi(\gamma(\tau))\dot\gamma_2(\tau)
\right)\,d\tau.
\end{split}\end{equation*}

Therefore, recalling the definition of $\Gamma$ in (<ref>) and the relations in (<ref>),

\begin{equation*}
\begin{split}&
\oint_\gamma w(z)\,dz=i
\int_0^L
\left( v_1(\gamma(\tau))\dot\gamma_1(\tau)+v_2(\gamma(\tau))\dot\gamma_2(\tau)
\right)\,d\tau=i\Gamma.
\end{split}\end{equation*}

As a result of this and (<ref>) we deduce that

$$ b_1=\frac{\Gamma}{2\pi }.$$

This and (<ref>) lead to

$$ w^2(z)=\left(
i \overline{v_0}+\frac{\Gamma}{2\pi z}+\sum_{k=2}^{+\infty}\frac{b_k}{ z^k}
\right)^2=-\overline{v_0}^2+\frac{i \overline{v_0} \Gamma}{\pi z}+\sum_{k=2}^{+\infty}\frac{c_k}{ z^k},$$

for suitable $c_k\in\C$.

This, Cauchy's Residue Theorem (see e.g. <cit.>) and (<ref>) give that

\begin{equation*}
\overline{\frac{2i F}\rho}
=2  \overline{v_0} \Gamma
\end{equation*}

and therefore

\begin{equation*}
{\frac{2i F}\rho}
=\overline{2\overline{v_0} \Gamma}=2v_0 \Gamma.
\end{equation*}

$$ F=-i\rho v_0 \Gamma$$

and this establishes the claim in (<ref>).

The level sets of the potential function in (<ref>).

A natural question now is however if it is possible to produce a nonzero circulation in (<ref>).
This is in general quite an intriguing problem: here we just provide a simple and explicit[Though
we do not really exploit this fact here, the example is inspired by a rotating cylinder producing
a Magnus effect, that is what soccer players use to make the ball curve during flight. See https://www.youtube.com/watch?v=XdL7EDKr_rk for a famous application of
the Magnus effect in soccer.] example.
For this, we let
\begin{equation}\label{HARGIRA4AXELHARLROA7789GIJ7soloDItangeFI2A}\Omega:=B_1\end{equation}

\begin{equation}\label{HARGIRA4AXELHARLROA7789GIJ7soloDItangeFI2} \psi(x,y):= y\left(1-\frac1{x^2+y^2}\right)+\ln(x^2+y^2).\end{equation}

See Figure <ref> for the level sets of $\psi$.
These level sets will correspond to streamlines, since we will consider the vector field $v=(v_1,v_2)$ with

\begin{equation}\label{HARGIRA4AXELHARLROAGIJ7soloDItangeFI2}
\begin{split}
\frac{2 x^2}{(x^2 + y^2)^2 }-\frac{1 + 2 y}{x^2 + y^2}-1
\\ {\mbox{and}}\qquad&
v_2(x,y):=\partial_x\psi(x,y)=\frac{2 x}{x^2 + y^2} +\frac{2 x y}{(x^2 + y^2)^2}.\end{split}
\end{equation}

See Figure <ref> for a sketch
of the vector field $v$.

The vector field in (<ref>).

The complex
analysis enthusiasts will also enjoy the fact that the setting in (<ref>)
is also available by posing

\begin{equation}\label{COMPLEXPOTE}w:=\frac2z + i \left(\frac1{z^2} - 1\right).\end{equation}

Notice that $w$ is holomorphic away of the origin. And the knowledge of this complex
map is pretty much all that is needed to produce the fluid velocity field in (<ref>).
The action of the above map $w$ is sketched in Figure <ref>.

We observe that $v$ is tangential to $\partial\Omega$ since

$$ v(\cos\theta,\sin\theta)\cdot(\cos\theta,\sin\theta)=0$$

and this is consistent with the impenetrability condition in (<ref>).
Additionally, we have that $\Delta\psi=0$ outside $\Omega$, which leads that $\div=0$
outside $\Omega$ as well.

Furthermore, recalling (<ref>),

\begin{equation}\label{8ihTGBAlFer3ONmdfZ345tyIjdmvREBSLv02}
\curl v=\big(\partial_xv_2-\partial_y v_1\big) e_3=\big(\partial_x(\partial_x\psi)+\partial_y \partial_y\psi\big) e_3=0
\end{equation}

\begin{eqnarray*}
\lim_{|(x,y)|\to+\infty} v(x,y)=
\lim_{|(x,y)|\to+\infty}\left(
\frac{2 x^2}{(x^2 + y^2)^2 }-\frac{1 + 2 y}{x^2 + y^2}-1,
\;\frac{2 x}{x^2 + y^2} +\frac{2 x y}{(x^2 + y^2)^2}\right)

The complex map $w$ in footnote <ref>.

Thus, to check that this example is indeed a solution of our fluid dynamics problem, it remains
to show that the first equation in (<ref>) holds true.
For this, we choose $\rho:=1$ and, inspired by the Bernoulli's Principle
in (<ref>),

$$ p(x,y):=\frac{|v_0|^2-|v(x,y)|^2}{2}=\frac{1-|v(x,y)|^2}{2}$$

and we combine the vectorial identity in (<ref>) with (<ref>)
to calculate that
\begin{eqnarray*}&&\rho (v\cdot\nabla) v+\nabla p=
(v\cdot\nabla) v-\nabla\frac{|v|^2}{2}=-
v\times (\curl v)
\end{eqnarray*}

This shows that the vector field in (<ref>)
is consistent with the fluid dynamics setting in (<ref>). To show the interest of (<ref>)
it thus remains to check that this vector field provides a nontrivial circulation $\Gamma$: to this end, we compute that

\begin{eqnarray*}
\Gamma&=&\int_0^{2\pi} v(\cos\theta,\sin\theta)\cdot(-\sin\theta,\cos\theta)\,d\theta\\
&=&2\int_0^{2\pi} (\cos^2\theta-\sin\theta-1,\;\cos\theta +\sin\theta\cos\theta)\cdot(-\sin\theta,\cos\theta)\,d\theta\\&=&2\int_0^{2\pi} (1+\sin\theta)\,d\theta\\&=&4\pi
\\&\ne&0,
\end{eqnarray*}

as desired.

But hold on a sec, the skeptical reader (who is always very welcome) will complain we've been cheating on them:
we have promised this section was devoted to airfoils, but then we tried to sell in (<ref>) the disk as an example of airfoils.
Come on, airfoils are objects as the ones in Figure <ref>, nothing in Figures <ref>, <ref> or <ref>
looks like an airfoil, just because the disk is not an airfoil at all!

Joukowski airfoil
(image by Krishnavedala from
Wikipedia, licensed under the Creative Commons Attribution-Share Alike 4.0 International license).

Well, in fact, to a certain degree, it is: this is the beauty of mathematics, and of complex analysis in particular.
Indeed, one of the brilliant ideas of the pioneers of aerodynamics is that the potential in (<ref>)
describes all the fluid dynamics outside the disk of the example that we have discussed explicitly in detail,
therefore other examples can be constructed via conformal transformations.
In particular, using dilations, translations and maps of the form $z+\frac1z$, one can transform a disk
into an airfoil shaped as in Figure <ref>.

The set of conformal transformations linking circles and airfoils can be visualized for instance via the Wolfram Demonstrations Projects

or via the GeoGebra application

or likely via a number of IT resources.

Interestingly, on
one can also plot the resulting flow around the airfoil: the result that we have obtained
using this application is reported in Figure <ref>,
which nicely shows how a disk is mapped into an airfoil.

From a disk to an airfoil, and the corresponding fluid flow (images produced
by the online application in

For further readings about the Kutta-Joukowski theory and the analysis of airfoils see e.g. <cit.>
and the references therein.
See also <cit.> for an extensive treatment of fluid dynamics from a complex analysis perspective.

We cannot avoid mentioning that one of the most spectacular applications of fluid dynamics is probably showcased
by the motion[The origin of the name boomerang is a bit uncertain: some references link it to an extinct Aboriginal language of New South Wales, others to the language of the Turuwal people (a sub-group of the Darug) of the Georges River (Tucoerah River).
Besides the well-known traditional employment by some Aboriginal Australian peoples for hunting, see Figure <ref>, it seems that boomerangs have been used also in ancient Europe, Egypt, and North America (ancient Egyptian boomerangs have been tested and seemed to work well as returning boomerangs and a boomerang discovered in the Carpathian Mountains in Poland dated back to about $30\times10^3$ years ago.

Interestingly, a boomerang
was used to set the Guinness World Record
for the longest throw of any object by a human: namely, in 2005, at Murarrie Recreation Ground, in Queensland,
David Schummy
performed a throw of 427.2 metres, see https://www.youtube.com/watch?v=ly3nCEbcQig

Currently, long distance boomerangs are mostly
shaped as a question mark and often have
a beveled edge, to facilitate the pitch
and lower the drift, since the boomerang in this case
is usually thrown almost horizontally (indeed, in the above mentioned
record, the objective
was not to make the boomerang come back to the throw location;
actually the boomerang ended up on a tree).] of a boomerang.

Otto Jungarryi Sims with a boomerang sitting inside a cave in the Northern Territory of Australia
(photo by Ed Gold; image from Wikipedia for free use under ticket #2020101210010454).

Roughly speaking, each wing of a boomerang is shaped as an airfoil section,
allowing the airflow over the wings to create a significant lift.
If the boomerang is thrown nearly upright, the rotating blades generate more lift at the top than the bottom,
because at the top the speed of the rotation adds up to the forward speed, while at the bottom
the speed of the rotation subtracts from the forward speed, see Figure <ref>
in which the forward speed is represented by the yellow arrow.
This additional lift from the top produces a torque (represented by the blue arrow in Figure <ref>)
whose effect is to make the rotation plane turn around: we point out that this torque is typically not sufficient
to tilt the boomerang around its axis of travel, given its high
angular momentum (the spinning of the boomerang being represented by the red arrow in Figure <ref>):
hence, the stability of the rotating plane ensured by the gyroscopic precession combines with the aerodynamic torque
and leads to the curved trajectory sketched by the green arrow in Figure <ref>.

Why do boomerangs return?

§.§ Surfing the waves

A topical problem in the dynamics of fluids consists in the description of waves in shallow water,
since this is a typical case arising in the proximity of land. Though a large number of different models are available for this goal,
in this pages we recall the classical equation

Students picnic on the Union Canal in 1922
(copyright The University of Edinburgh, available for public use; image from
licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license).

\begin{equation}\label{KDV-or}
a\partial _{t}\phi +b\partial _{x}^{3}\phi +c\,\phi \,\partial _{x}\phi =0
\end{equation}

for $a$, $b$, $c\in\R\setminus\{0\}$, with $\phi=\phi(x,t)$, $x\in\R$ and $t\in(0,+\infty)$.
Equation (<ref>) is called[Equation (<ref>)
is named after Diederik Johannes Korteweg and Gustav de Vries <cit.>, though it was
introduced by Joseph Valentin Boussinesq <cit.>.
De Vries completed his PhD under Korteweg's supervision and then worked all his life
as a high school teacher in Haarlem in the Netherlands.

The strong interest in the formation and propagation of waves in canals was possibly the outcome of a direct observation by
John Scott Russell <cit.> which took place in the Union Canal
(a canal in Scotland, running from Falkirk to Edinburgh).
Russell's own words have been repeated in virtually all papers and books which even remotely discuss water wave problems and we have no intention of breaking this consolidated tradition, hence here is his report on the astounding experience of meeting a traveling wave for the first time: “I was observing the motion of a boat which was rapidly drawn along a narrow channel by a pair of horses, when the boat suddenly stopped – not so the mass of water in the channel which it had put in motion; it accumulated round the prow of the vessel in a state of violent agitation, then suddenly leaving it behind, rolled forward with great velocity, assuming the form of a large solitary elevation, a rounded, smooth and well-defined heap of water, which continued its course along the channel apparently without change of form or diminution of speed. [...] Such, in the month of August 1834, was my first chance interview with that singular and beautiful phenomenon which I have called the Wave of Translation”.

See Figure <ref> for a historical picture
of the Union Canal.] the
Korteweg-de Vries equation (or KdV equation for short). Korteweg-de Vries equation

The physical content of equation (<ref>) is that $\phi$ models the shape of a wave in a canal.
The position on the canal is given by $x\in\R$ and $t$ stands for the time variable.
The amplitude of the wave and the depth of the channel are supposed to be small
(as will be discussed below in the approximations leading to the derivation of the equation),
hence the validity of (<ref>) has to be limited to small amplitude waves in shallow waters.

Though equation (<ref>) should be considered as a very simplified model, not able to capture the complexity
of oscillatory phenomena and waves in the real world, it is interesting to notice that the Korteweg-de Vries equation
can describe some interesting features such as traveling waves
and provide interesting information on the shape of the waves, as we now discuss.

The traveling wave solution for the Korteweg-de Vries equation equation, as found in
in (<ref>).

First of all (recalling the presentation on page 01ojehniTNasIplamnewa) one can seek for
traveling wave solutions of (<ref>) in the form

\begin{equation}\label{BSNEGIAaon}
\phi(x,t)=\phi_0(x-vt)\end{equation}

for some velocity $v$.
In this setting, equation (<ref>) reduces to

\begin{equation}\label{KAlyhchtvdhr4htFAn}
0=-av\phi_0' +b\phi'''_0 +c\,\phi_0 \,\phi_0' =-av\phi_0' +b\phi'''_0 +\frac{c}2\,(\phi_0^2)'=\left(-av\phi_0 +b\phi''_0 +\frac{c}2\,\phi_0^2\right)'.
\end{equation}

We integrate this equation by assuming that at infinity $\phi_0$
and its derivatives converge to zero, finding that

\begin{equation*}0=
-av\phi_0 +b\phi''_0 +\frac{c}2\,\phi_0^2.
\end{equation*}

As a result,

\begin{equation*}
0=\left( -av\phi_0 +b\phi''_0 +\frac{c}2\,\phi_0^2\right)\phi_0'=
\left(-\frac{av}2\phi_0^2 +\frac{b}2 (\phi'_0)^2 +\frac{c}{6}\,\phi_0^3 \right)'.
\end{equation*}

Integrating this as above, we deduce that

\begin{equation*} 0
=-\frac{av}2\phi_0^2 +\frac{b}2 (\phi'_0)^2+\frac{c}{6}\,\phi_0^3.\end{equation*}

This ordinary differential equation
can be solved by separation of variables (assuming $\phi_0\ge0$, $c>0$ and $abv>0$) leading,
up to a translation, to

\begin{equation}\label{DIKDVSOLIFI-EQ} \phi_0(r)=
\frac{6 a v}{c \left( \cosh\left(r \sqrt{\frac{a v}b} \right)+ 1\right)}.\end{equation}

See Figure <ref> for a sketch[See also http://lie.math.brocku.ca/$\sim$sanco/solitons/kdv_solitons.php
for several animations of (possibly interacting) traveling waves
of the Korteweg-de Vries equation.] of this traveling wave (recall (<ref>))
when $a:=1$, $b:=1$, $c:=1$ and $v:=1$.

It is also interesting to observe that the highest crest of the wave in (<ref>) is attained for $r:=0$
and it is equal to

\begin{equation}\label{d3ae9ikjm3f82KSMpf334}
\frac{3 a v}{c}.\end{equation}
This suggests that the fastest waves (i.e., waves with larger velocity $v$) correspond to the
highest ones.

Wave trains crossing in front of Île de Ré, in the Atlantic Ocean
(photo by Michel Griffon; image from
licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license).

Another nice feature captured by  the Korteweg-de Vries equation in (<ref>)
is a sufficiently accurate description of the shape of a set of waves: it is indeed a common experience that
“real” waves exhibit sharper crests and flatter troughs than those of sine and cosine functions,
see Figures <ref> and <ref>.
To understand the shape of these waves, it is convenient to define, given $\varphi\in\R$ and $m\in[0,1)$,

\begin{equation}\label{NokmsLPSo8s0} U(\varphi,m):=\int _{0}^{\varphi }{\frac{d\theta }{\sqrt {1-m\sin ^{2}\theta }}}.\end{equation}

We observe that $\frac{\partial U}{\partial \varphi}>0$, hence, given $m\in(0,1)$, the function $\varphi\mapsto U(\varphi,m)$ is invertible
and, with a slight abuse of notation, we call $\varphi(U,m)$ the inverse function.
In this setting, one defines the elliptic cosine (in Latin, cosinus amplitudinis) by elliptic cosine

\begin{equation}\label{NokmsLPSo8s0BIS} \operatorname{cn} (U,m):=\cos (\varphi(U,m)).\end{equation}

Now, to find periodic solutions of the Korteweg-de Vries equation in (<ref>)
and to compare their shapes with Figure <ref> we go back to (<ref>)
and we integrate it, now not assuming that $\phi_0$ decays at infinity: in this way, we have that

$$ -av\phi_0 +b\phi''_0 +\frac{c}2\,\phi_0^2=\kappa_1,$$

for some constant $\kappa_1\in\R$.

As a consequence,

$$0=\Big(-av\phi_0 +b\phi''_0 +\frac{c}2\,\phi_0^2-\kappa_1\Big)\phi_0'=
\left(-\frac{av}2\phi_0^2 +\frac{b}2 (\phi'_0)^2 +\frac{c}6\,\phi_0^3-\kappa_1\phi_0\right)'

from which a further integration produces

\begin{equation}\label{MS-pqewkrfciaub} -\frac{av}2\phi_0^2 +\frac{b}2 (\phi'_0)^2 +\frac{c}6\,\phi_0^3-\kappa_1\phi_0=\kappa_2,\end{equation}

for some constant $\kappa_2\in\R$.

We observe that we can write (<ref>) as

\begin{equation}\label{MS-pqewkrfciaub2} (\phi'_0)^2={\mathcal{P}}(\phi_0),\end{equation}

for a suitable polynomial ${\mathcal{P}}$ of degree $3$. The setting that we consider now
is that in which ${\mathcal{P}}$ possesses three distinct real roots $\beta_1>\beta_2>\beta_3$.
With this, we rewrite (<ref>) as

\begin{equation}\label{MS-pqewkrfciaub-x}

where $\kappa:=c/(3b)$, that we assume to be positive.

Now it turns out to be useful to seek solutions in the form

\begin{equation}\label{MS-pqewkrfciaub-x98ytr} \phi_0(r)=\beta_1\cos^2 (\Phi(r))+\beta_2\sin^2 (\Phi(r)), \end{equation}

where the function $\Phi$ is supposed to be increasing and has to be determined.

We observe that

$$ (\phi_0')^2=\Big( -2\beta_1\cos\Phi\,\sin\Phi\, \Phi'+2\beta_2\sin\Phi \,\cos\Phi\,\Phi'\Big)^2=
4(\beta_1-\beta_2)^2 \sin^2\Phi\,\cos^2\Phi\,(\Phi')^2.$$

\begin{eqnarray*}
\phi_0-\beta_1=\beta_1\big( \cos^2 \Phi -1\big)+\beta_2\sin^2 \Phi=-(\beta_1-\beta_2)\sin^2 \Phi
\end{eqnarray*}

\begin{eqnarray*}
\phi_0-\beta_2=\beta_1 \cos^2 \Phi +\beta_2\big(\sin^2 \Phi-1)=(\beta_1-\beta_2)\cos^2 \Phi.
\end{eqnarray*}

These observations and (<ref>) yield that

\begin{equation*}
\end{equation*}

Thus, since

$$ \phi_0-\beta_3=\beta_1\Big(1-\sin^2 \Phi\Big)+\beta_2\sin^2 \Phi-\beta_3=(\beta_1-\beta_3)-

we conclude that

\begin{equation}\label{MS-pqewkrfciaub-x89} \Phi'=\frac{\sqrt\kappa}2\,\sqrt{ (\beta_1-\beta_3)-(\beta_1-\beta_2)\sin^2\Phi}=
\kappa_0\,\sqrt{ 1-m\sin^2\Phi}

$$ \kappa_0:=\frac{\sqrt{\kappa\,(\beta_1-\beta_3)}}2\qquad{\mbox{and}}\qquad

We can thus utilize the separation of variable method in (<ref>)
and, in view of (<ref>), obtain that

$$ \kappa_0 r=\int_{0}^{\Phi(r)}\frac{d\vartheta}{\sqrt{ 1-m\sin^2\vartheta}}=U\big(\Phi(r),m\big),

where we normalized the picture so that $\Phi(0)=0$.

From this and (<ref>) it follows that

$$ \operatorname{cn}(\kappa_0 r)=\operatorname{cn}\big( U\big(\Phi(r),m\big)\big)=\cos(\Phi(r)).

Going back to (<ref>) we thereby conclude that

\begin{equation}\label{EQcnoidal-waves} \begin{split}&\phi_0(r)=\beta_1\cos^2 (\Phi(r))+\beta_2\big(1-\cos^2 (\Phi(r))\big)
=(\beta_1-\beta_2)\cos^2 (\Phi(r))+\beta_2\\&\qquad\quad\qquad=(\beta_1-\beta_2)\operatorname{cn}^2(\kappa_0 r )+\beta_2=
A\operatorname{cn}^2(\kappa_0 r)+B,\end{split}
\end{equation}

where $A:=\beta_1-\beta_2$ and $B:=\beta_2$.
The profile of these functions, for suitable choices of parameters, resembles the shape of waves
in Figures <ref> and <ref> (which are often
called “cnoidal waves” cnoidal wave due to the presence of the elliptic cosine $\operatorname{cn}$
in their expression). To have a feeling on how the different parameters
change the shape of a cnoidal waves see Figure <ref>.

Cnoidal wave solution to the Korteweg-de Vries equation (image by Kraaiennest from
licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license).

It is also interesting to take note of the fact that, in spite of its (relative) simplicity, the Korteweg-de Vries equation
in (<ref>) is capable of capturing complex phenomena such as the formation of oscillatory waves, with crests
traveling at different speeds according to their heights: these complicated situations may even arise from very simple
initial data, see Figure <ref>.

Furthermore, we point out that the Korteweg-de Vries equation often turns in handy
in the description of natural phenomena that seem to be well beyond its range of applicability.
For example (see e.g. <cit.> and the references therein), the Korteweg-de Vries equation
has been exploited in the description[Although it is debatable whether the artists
intended their works to be interpreted in that way, it is customary to
refer to the artworks in Figure <ref>
in connection with tsunamis.]
of tsunamis tsunami
even though the extreme depth of oceans may appear to be incompatible with an equation which was originally
designed for shallow waters. The reason for the success of the Korteweg-de Vries equation even in this setting
is, at least, twofold. First of all, though oceans are very deep, tsunami waves can reach a spacial
extent of even larger size: therefore, when the depth of the ocean is sufficiently small with respect to the breadth of the wave
the shallow water theory may still be applied, since what counts is the relative (and not absolute) sizes of the length parameters.
Secondly, while tsunamis in the deep ocean consists typically of very long waves with quite small amplitudes, as they approach shallow waters the nonlinear effects become predominant and they change significantly
the shape and velocity of the leading waves, for which the Korteweg-de Vries equation can also provide an approximate, but often
effective, model.

Now we give a brief motivation for the Korteweg-de Vries equation
in (<ref>). While its expression may appear rather mysterious
at first sight, especially in view of the third derivative appearing as a leading term, we will
see now that equation (<ref>) arises naturally from the equations of fluid dynamics
presented in Section <ref>.

More specifically, we consider the Euler fluid flow equation in (<ref>).
Here we take coordinates $(x,z)\in\R^2$ where $x\in\R$ represents the direction of a canal
and $z\in\R$ is the vertical displacement. The fluid velocity $v$ will be written in components as $v=(\overline u,\overline w)$,
with $\overline u$ and $\overline w$ scalar functions. Hence, being $\rho$ the density of the fluid (assumed to be constant)
and $p$ the pressure, equation (<ref>) reads

\begin{equation}\label{Cm-PAPMam}
\rho \partial_t (\overline u,\overline w)+\rho ((\overline u,\overline w)\cdot(\partial_x,\partial_z)) (\overline u,\overline w)=
-\rho g(0,1) -(\partial_x,\partial_z) p.
\end{equation}

Hence, setting

\begin{equation}\label{Cm-PAPMam20} \overline P:=\frac{p}\rho+gz,
\end{equation}

we rewrite (<ref>) as the following system of two equations:

\begin{equation}\label{Cm-PAPMam2}
\begin{dcases}
\partial_t \overline u+
\overline u\partial_x\overline u+\overline w \partial_z \overline u=- \partial_x\overline P,
\\
\partial_t \overline w+
\overline u \partial_x\overline w
+\overline w \partial_z \overline w=
-\partial_z\overline P.
\end{dcases}
\end{equation}

We suppose that the bottom of the canal is $z=0$ and the fluid at rest has height $h_0\in(0,+\infty)$.
The wave is thus described as a perturbation of the rest fluid surface $z=h_0$
described by the graph of a function $\overline\eta=\overline\eta(x,t)$, modulated by a small amplitude $\alpha\in(0,+\infty)$:
namely, at every instant of time $t$ the fluid in the channel corresponds to the region of the points $(x,z)\in\R^2$
with $0<z<h_0+\alpha\overline\eta(x,t)$. That is, the equations in (<ref>) are set in the domain

\begin{equation}\label{Cm-PAPMam30} \Omega^{(t)}:=\Big\{ (x,z)\in\R^2{\mbox{ s.t. }}x\in\R{\mbox{ and }}z\in\big(0,h_0+\alpha\overline\eta(x,t)\big)\Big\}.\end{equation}

We also assume that the fluid is incompressible, hence, by (<ref>),

\begin{equation}\label{Cm-PAPMam3}
0=\div (\overline u,\overline w)=\partial_x\overline u+\partial_z\overline w.\end{equation}

The equations in (<ref>) and (<ref>) are complemented by boundary conditions
corresponding to the bottom of the canal and the free surface of the fluid.
Indeed, at the bottom $z=0$ we suppose that the floor of the canal is impenetrable, hence the vertical component
of the fluid velocity vanishes, which reads

\begin{equation}\label{Cm-PAPMam4}
\overline w(x,0,t)=0.\end{equation}

As for the free surface, we assume that fluid parcels staying at the top remain at the top:
hence if $(x(t),z(t))$ represents the trajectory of a fluid particle on the upper surface of the fluid at time $t$ we have that

$$ z(t)=h_0+\alpha\overline\eta(x(t),t).$$

Thus, taking the time derivative of this expression and recalling that the fluid parcel velocity $(\dot x(t),\dot z(t))$ agrees with $(\overline u,\overline w)$,

\begin{equation}\label{Cm-PAPMam44}
\overline w\big(x(t),h_0+\alpha\overline\eta(x(t),t),t\big)=\dot z(t)=\alpha\partial_x\overline\eta(x(t),t)\dot x(t)+\alpha\partial_t\overline\eta(x(t),t).\end{equation}

The cnoidal waves in (<ref>)
with $A:=1$, $B:=3$ and $(\kappa_0,m)\in\big\{ (1/7,0.9),\,(1/6,0.99),\,(1/5,0.999),\,(1/3,0.999999)
\big\}$.

The free surface $z=h_0+\alpha\overline\eta$ also presents an additional boundary condition related to pressure.
Indeed, the pressure on the top surface of the fluid is balanced by the atmospheric pressure, that we denote by $p_0$
and assume to be constant. Hence, by (<ref>),

\begin{equation}\label{Cm-PAPMam45} \overline P\big(x,h_0+\alpha\overline\eta(x,t),t\big)=\frac{p_0}\rho+g\,\big(h_0+\alpha\overline\eta(x,t)\big).\end{equation}

We now consider a regime of shallow waters and small amplitudes.
That is,
we consider the quantities

\begin{equation}\label{EPSIDELKDV}
\e:=\frac{\alpha}{h_0}\qquad{\mbox{and}}\qquad\delta:=\frac{h_0}{g}.\end{equation}

We suppose that $\e$ and $\delta$ are small and also, up to normalizing constants, that

\begin{equation}\label{sawqewet54b5494pppppplllllkkkk}
\delta=\e.\end{equation}

We define

\begin{equation}\label{RTArasmsImPNSd-2}\begin{split}& u^\star(x,z,t):=\frac{\overline u\big(\sqrt{gh_0} x,h_0z,t\big)}{\sqrt{gh_0}},\qquad
w^\star(x,z,t):=\frac{\overline w\big(\sqrt{gh_0} x,h_0z,t\big)}{h_0},\\&\qquad
\eta^\star(x,t):= {\overline\eta(\sqrt{gh_0} x,t)}
\qquad{\mbox{and}}\qquad
P^\star(x,z,t):=\frac{\overline P\big(\sqrt{gh_0} x,h_0z,t\big)}{{gh_0}}.
\end{split}\end{equation}

In view of (<ref>), the above functions are defined in the region

\begin{equation}\label{NEWPD-do1}
\begin{split}
\Omega_{ t }^\star&:=
\Big\{ (x,z)\in\R^2{\mbox{ s.t. }}(\sqrt{gh_0} x,h_0z)\in\Omega^{(t)}\Big\}\\&
=\Big\{ (x,z)\in\R^2{\mbox{ s.t. }}x\in\R{\mbox{ and }} z\in\big(0,1+\e\eta^\star( x,t)\big)\Big\}.\end{split}
\end{equation}

We thus collects the equations in (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>), exploiting also (<ref>) and (<ref>), by writing

\begin{equation}\label{w78uhdbc-384irton2rve7g}
\begin{dcases}
\partial_t u^\star+
u^\star\partial_xu^\star+w^\star \partial_z u^\star=- \partial_x P^\star
& {\mbox{ in }}\Omega_{ t }^\star,
\\
\e\big( \partial_t w^\star+
u^\star \partial_x w^\star
+ w^\star \partial_z w^\star\big)=
-\partial_z P^\star & {\mbox{ in }}\Omega_{ t }^\star,\\ \partial_xu^\star+\partial_zw^\star=0& {\mbox{ in }}\Omega_{ t }^\star,\\
P^\star=\displaystyle\frac{p_0}{\rho gh_0}+{1+\e\eta^\star }
& {\mbox{ on }} z=1+\e\eta^\star ,\\
& {\mbox{ on }} z=1+\e\eta^\star ,\\
w^\star=0& {\mbox{ on }} z=0.
\end{dcases}
\end{equation}

We let $c_0:=\frac{p_0}{\rho gh_0}+1$.
It is now convenient to “surf the wave” and look at the new variable $\xi:=x-t$
(roughly speaking, using this variable one is moving with a wave that travels at unit[Interestingly,
moving at unit speed for $\eta^\star$ corresponds to moving at speed $\sqrt{gh_0}$
for the original profile $\overline\eta$, thanks to the change of spatial variables in (<ref>).
The quantity $\sqrt{gh_0}$ happens indeed to be one of the characteristic velocities of small waves
in shallow waters (relating the velocity to $\sqrt{gh_0}$
also tells us that the speed of these waves decreases when the height of the canal
is smaller and, conversely, their speed increases when moving from very shallow to slightly deeper water).] speed).
It is also appropriate to look at large times by setting $\tau:=\e t$. Thus, we rephrase our physical quantities with
respect to these new variables by defining

\begin{equation}\label{RTArasmsImPNSd-1}
\begin{split}&
u(\xi,z,\tau):=u^\star\left(\xi+\frac\tau\e,z,\frac\tau\e\right),\qquad w(\xi,z,\tau)=
\eta^\star\left(\xi+\frac\tau\e,\frac\tau\e\right)
\qquad{\mbox{and}}\qquad P(\xi,z,\tau):=
\end{split}\end{equation}

By (<ref>), these functions are defined in the region

\begin{eqnarray*} \Omega_\tau&:=&\Big\{ (\xi,z)\in\R^2{\mbox{ s.t. $\xi=x-t$, $\tau=\e t$ and~$(x,z)\in\Omega_{ t }^\star$}}\Big\}\\&
=&\Big\{ (\xi,z)\in\R^2{\mbox{ s.t. }}\xi\in\R{\mbox{ and }} z\in\big(0,1+\e\eta( \xi,\tau)\big)\Big\}

In this setting, one can rewrite (<ref>) in the form

\begin{equation}\label{KDV-0m-001}
\begin{dcases}
-\partial_\xi u+\e\partial_\tau u+u\partial_\xi u+w\partial_z u=-\partial_\xi P&{\mbox{ in }}\Omega_\tau,\\
\e\big( -\partial_\xi w+\e\partial_\tau w+u\partial_\xi w+w\partial_z w\big)=-\partial_z P&{\mbox{ in }}\Omega_\tau,\\
\partial_\xi u+\partial_z w=0&{\mbox{ in }}\Omega_\tau,\\
P=c_0+\e\eta & {\mbox{ on }}z=1+\e\eta,\\
w=\e(\e\partial_\tau\eta+(u-1)\partial_\xi\eta)& {\mbox{ on }}z=1+\e\eta,\\
w=0 & {\mbox{ on }}z=0.
\end{dcases}
\end{equation}

Numerical solution of the KdV equation $\partial_t u + u\partial_x u + (0.022)^2\partial_x^3u = 0$
with initial condition $ u(x, 0) = \cos(\pi x)$.
The initial cosine wave evolves into a train of solitary-type waves <cit.>
(by Ta2o; image from
licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license).

It is now convenient to consider formal Taylor
expansions of the physical quantities involved in powers of $\e$ by writing

\begin{eqnarray*}&& u(\xi,z,\tau)=\sum_{k=0}^{+\infty}\e^k u_k(\xi,z,\tau),\qquad
w(\xi,z,\tau)=\sum_{k=0}^{+\infty}\e^k w_k(\xi,z,\tau)\\&&
P(\xi,z,\tau)=\sum_{k=0}^{+\infty}\e^k P_k(\xi,z,\tau)
\qquad{\mbox{and}}\qquad
\eta(\xi,\tau)=\sum_{k=0}^{+\infty}\e^k \eta_k(\xi,\tau).\end{eqnarray*}

Substituting these expressions into (<ref>) we obtain that

\begin{equation}\label{KDV-0m-002}
\begin{dcases}
-\displaystyle\sum_{k=0}^{+\infty}\e^k\partial_\xi  u_k+\displaystyle\sum_{k=0}^{+\infty}\e^{k+1}\partial_\tau u_k+
\displaystyle\sum_{k,m=0}^{+\infty}\e^{k+m}
u_k\partial_\xi u_m+\displaystyle\sum_{k,m=0}^{+\infty}\e^{k+m}
w_k\partial_z u_m=-\displaystyle\sum_{k=0}^{+\infty}\e^k\partial_\xi P_k&{\mbox{ in }}\Omega_\tau,\\
-\displaystyle\sum_{k=0}^{+\infty}\e^{k+1}\partial_\xi w_k+
\displaystyle\sum_{k=0}^{+\infty}\e^{k+2}\partial_\tau w_k+
\displaystyle\sum_{k,m=0}^{+\infty}\e^{k+m+1}u_k\partial_\xi w_m+
\displaystyle\sum_{k,m=0}^{+\infty}\e^{k+m+1}w_k\partial_z w_m=-\displaystyle\sum_{k=0}^{+\infty}\e^k\partial_z P_k&{\mbox{ in }}\Omega_\tau,\\
\displaystyle\sum_{k=0}^{+\infty}\e^k(\partial_\xi u_k+\partial_z w_k)=0&{\mbox{ in }}\Omega_\tau,\\
P\left( \xi, 1+\displaystyle\sum_{k=0}^{+\infty}\e^{k+1}\eta_k,\tau\right)=c_0+\displaystyle\sum_{k=0}^{+\infty}\e^{k+1}\eta_k, & \\
w\left( \xi, 1+\displaystyle\sum_{k=0}^{+\infty}\e^{k+1}\eta_k,\tau\right)=
\displaystyle\sum_{k=0}^{+\infty}\e^{k+2}\partial_\tau\eta_k+
\displaystyle\sum_{k,m=0}^{+\infty}\e^{k+m+1}u_k\partial_\xi\eta_m-
\displaystyle\sum_{k=0}^{+\infty}\e^{k+1}\partial_\xi\eta_k,& \\
w(\xi,0,\tau)=0. &
\end{dcases}
\end{equation}

The strategy is now to approximate (<ref>) by neglecting the terms of size $\e^k$ with $k\ge3$.
For this, we consider the orders in $\e^k$ in (<ref>) with $k\in\{0,1,2\}$ and we formally equate the corresponding
coefficients. Namely, by formally considering the coefficients related to $\e^0$ in (<ref>) we see that

\begin{equation*}
\begin{dcases}
-\partial_\xi  u_0+u_0\partial_\xi u_0+w_0\partial_z u_0=-\partial_\xi P_0&{\mbox{ in }}\Omega_\tau,\\
0=\partial_z P_0&{\mbox{ in }}\Omega_\tau,\\
\partial_\xi u_0+\partial_z w_0=0&{\mbox{ in }}\Omega_\tau,\\
P_0( \xi, 1,\tau)=c_0, & \\
w_0( \xi, 1,\tau)=0,& \\
w_0(\xi,0,\tau)=0. &
\end{dcases}
\end{equation*}

Since we are interested in small perturbations of steady states, we also take $u_0:=0$ and $w_0:=0$, therefore

\begin{equation}\label{KDV-0m-003}
\begin{dcases}
0=-\partial_\xi P_0&{\mbox{ in }}\Omega_\tau,\\
0=\partial_z P_0&{\mbox{ in }}\Omega_\tau,\\
P_0( \xi, 1,\tau)=c_0. &
\end{dcases}
\end{equation}

From the second and third lines in (<ref>) it follows that

\begin{equation}\label{KDV-0m-005}P_0(\xi,z,\tau)=c_0.\end{equation}

The Great Wave off Kanagawa, famous woodblock print by the ukiyo-e artist Katsushika Hokusai
(H. O. Havemeyer Collection, Bequest of Mrs. H. O. Havemeyer, 1929, Metropolitan Museum of Art; Public Domain
image from Wikipedia). Right: Die Woge, sculpture by Tobias Stengel in Dresden, Germany (photo by Christoph Münch,
Press Section of the city of Dresden's official homepage; image from
licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license).

Now we consider the coefficients related to $\e^1$ in (<ref>) and we obtain that

\begin{equation}\label{KDV-0m-004}
\begin{dcases}
-\partial_\xi  u_1=-\partial_\xi P_1&{\mbox{ in }}\Omega_\tau,\\
0=-\partial_z P_1&{\mbox{ in }}\Omega_\tau,\\
\partial_\xi u_1+\partial_z w_1=0&{\mbox{ in }}\Omega_\tau,\\
P_1(\xi,1,\tau)+\partial_z P_0(\xi,1,\tau)\eta_0=\eta_0, & \\
w_1(\xi,1,\tau)=-\partial_\xi\eta_0,& \\
w_1(\xi,0,\tau)=0. &
\end{dcases}
\end{equation}
Owing to (<ref>),
the fourth line in (<ref>) can be written as $P_1(\xi,1,\tau)=\eta_0(\xi,\tau)$, which
combined with the second line gives that

\begin{equation}\label{KDV-0m-009}
Hence, from the first line,

\begin{equation}\label{TH:EVE}
\partial_\xi u_1=\partial_\xi\eta_0.\end{equation}
This and the third line give that

$$ \partial_z (w_1+z\partial_\xi\eta_0)=\partial_z w_1+\partial_\xi\eta_0=-\partial_\xi u_1+\partial_\xi\eta_0=0

and therefore, by the last line in (<ref>), we infer that

\begin{equation}\label{KDV-0m-006}w_1=-z\partial_\xi\eta_0,\end{equation}

which is also in agreement with the second last line in (<ref>).

Additionally, from (<ref>), we have that $u_1(\xi,\tau,z)=\eta_0(\xi,\tau)+\upsilon(\tau,z)$,
for an auxiliary function $\upsilon$.
As a matter of fact, if we have a solution of (<ref>)
and we add to $u_1$ a function of $(\tau,z)$ we obtain another solution. Therefore,
we choose the “simplest possible” solution by picking $\upsilon$ identically equal to zero, thus obtaining

\begin{equation}\label{TH:EVE:2}
\end{equation}

Now we consider the order related to $\e^2$ in the expansion in (<ref>)
(the higher orders being formally neglected). In this way, we find that

\begin{equation}\label{KDV-0m-008}
\begin{dcases}
-\partial_\xi  u_2+\partial_\tau u_1+
u_1\partial_\xi u_1+w_1\partial_z u_1=-\partial_\xi P_2&{\mbox{ in }}\Omega_\tau,\\
-\partial_\xi w_1=-\partial_z P_2&{\mbox{ in }}\Omega_\tau,\\
\partial_\xi u_2+\partial_z w_2=0&{\mbox{ in }}\Omega_\tau,\\
\partial_zP_0(\xi,1,\tau)\eta_1+
\displaystyle\frac12\partial_{zz}P_0(\xi,1,\tau)\eta_0^2+\partial_zP_1(\xi,1,\tau)\eta_0+P_2(\xi,1,\tau)=
\eta_1, & \\
\partial_zw_1(\xi,1,\tau)\eta_0+w_2(\xi,1,\tau)=
\partial_\tau\eta_0+u_1\partial_\xi\eta_0-\partial_\xi\eta_1,& \\
w_2(\xi,0,\tau)=0. &
\end{dcases}
\end{equation}

Recalling (<ref>) and (<ref>), the fourth line in (<ref>) becomes

$$ P_2(\xi,1,\tau)=
\eta_1(\xi,\tau).

Also, from the second line in (<ref>) and (<ref>),

$$ \partial_zP_2=-\partial_\xi(z\partial_\xi\eta_0)=-z\partial_{\xi \xi}\eta_0.

These observations entail that, for all $z\in(0,1)$,

\begin{eqnarray*}&&\eta_1(\xi,\tau)-P_2(\xi,z,\tau)=
\int_z^1\partial_zP_2(\xi,z',\tau)\,dz'\\&&\qquad=-\int_z^1 z'\partial_{\xi \xi}\eta_0(\xi,\tau)\,dz'=-\frac{1-z^2}2\partial_{\xi \xi}\eta_0(\xi,\tau)
\end{eqnarray*}

and therefore

$$ P_2(\xi,z,\tau)=\eta_1(\xi,\tau)+\frac{1-z^2}2\partial_{\xi \xi}\eta_0(\xi,\tau).$$

From this and the first equation in (<ref>) we infer that

\begin{eqnarray*}
\partial_\xi  u_2&=&
\partial_\tau u_1+u_1\partial_\xi u_1+w_1\partial_z u_1+\partial_\xi P_2\\
\partial_\tau u_1+u_1\partial_\xi u_1+w_1\partial_z u_1+\partial_\xi \eta_1+\frac{1-z^2}2\partial_{\xi\xi \xi}\eta_0.
\end{eqnarray*}

Notice the appearance of a third derivative for $\eta_0$. Thus, recalling
the third equation in (<ref>),

\begin{eqnarray*}
-\partial_z w_2=
\partial_\xi u_2=
\partial_\tau u_1+u_1\partial_\xi u_1+w_1\partial_z u_1+\partial_\xi \eta_1+\frac{1-z^2}2\partial_{\xi\xi \xi}\eta_0.
\end{eqnarray*}

Combining this with the last two equations in (<ref>) we conclude that

\begin{eqnarray*}&&
\partial_\tau\eta_0(\xi,\tau)+u_1(\xi,1,\tau)\partial_\xi\eta_0(\xi,\tau)-\partial_\xi\eta_1(\xi,\tau)-\partial_zw_1(\xi,1,\tau)\eta_0(\xi,\tau)\\&=&w_2(\xi,1,\tau)\\
&=& w_2(\xi,1,\tau)-w_2(\xi,0,\tau)
\\&=&\int_0^1 \partial_zw_2(\xi,z,\tau)\,dz\\&=&
-\int_0^1 \Bigg(
\partial_\tau u_1(\xi,z,\tau)+u_1(\xi,z,\tau)\partial_\xi u_1(\xi,z,\tau)+w_1(\xi,z,\tau)\partial_z u_1(\xi,z,\tau)\\&&
\qquad\qquad+\partial_\xi \eta_1(\xi,\tau)+\frac{1-z^2}2\partial_{\xi\xi \xi}\eta_0(\xi,\tau)
\Bigg)\,dz
\end{eqnarray*}

This and (<ref>) yield that

\begin{eqnarray*}&&
\partial_\tau\eta_0(\xi,\tau)+u_1(\xi,1,\tau)\partial_\xi\eta_0(\xi,\tau)-\partial_\xi\eta_1(\xi,\tau)+\partial_\xi\eta_0(\xi,\tau)\eta_0(\xi,\tau)\\&=&
-\int_0^1 \left(
\partial_\tau u_1(\xi,z,\tau)+u_1(\xi,z,\tau)\partial_\xi u_1(\xi,z,\tau)-z\partial_\xi\eta_0(\xi,\tau)\partial_z u_1(\xi,z,\tau)
\right)\,dz\\&&\qquad\qquad-\partial_\xi \eta_1(\xi,\tau)-\frac13\partial_{\xi\xi \xi}\eta_0(\xi,\tau).
\end{eqnarray*}

Then, we can simplify the term $-\partial_\xi\eta_1(\xi,\tau)$ and exploit (<ref>), finding that

\begin{eqnarray*}&&
\partial_\tau\eta_0(\xi,\tau)+\eta_0(\xi,\tau)\partial_\xi\eta_0(\xi,\tau)+\partial_\xi\eta_0(\xi,\tau)\eta_0(\xi,\tau)\\&=&
-\int_0^1 \left(
\partial_\tau \eta_0(\xi,\tau)+\eta_0(\xi,\tau)\partial_\xi \eta_0(\xi,\tau)
\right)\,dz-\frac13\partial_{\xi\xi \xi}\eta_0(\xi,\tau)\\&=&
-\partial_\tau \eta_0(\xi,\tau)-\eta_0(\xi,\tau)\partial_\xi \eta_0(\xi,\tau)
-\frac13\partial_{\xi\xi \xi}\eta_0(\xi,\tau),
\end{eqnarray*}

which is the Korteweg-de Vries equation

\begin{equation}\label{VJlaVDa8jr3} \partial_\tau\eta_0+ \frac32\partial_\xi\eta_0\eta_0+\frac16\partial_{\xi\xi \xi}\eta_0=0,\end{equation}

to be compared with (<ref>), as desired.

It is now instructive to recover equation (<ref>) in terms of the original variables. For this, owing to (<ref>) and (<ref>),

\begin{equation}\label{VJlaVDa8jr3dra6mo9}\eta(\xi,\tau)=
\eta^\star\left(\xi+\frac\tau\e,\frac\tau\e\right)=\overline\eta\left(\sqrt{gh_0} \left(\xi+\frac\tau\e\right),\frac\tau\e\right).\end{equation}

Moreover, since $\overline\eta(x,t)$ was the profile of the free surface of the fluid
in the canal and $\eta_0(\xi,\tau)$ the leading order of this profile in the new variables $(\xi,\tau)$
in a shallow water and low amplitude approximation, in light of (<ref>) we can consider the leading order of
the original profile $\overline\eta_0$ as defined implicitly by the relation

$$ \eta_0(\xi,\tau)=\overline\eta_0\left(\sqrt{gh_0} \left(\xi+\frac\tau\e\right),\frac\tau\e\right).$$

Therefore (<ref>) becomes

\begin{eqnarray*}\frac{
\sqrt{gh_0}}\e\partial_x\overline\eta_0+\frac1\e\partial_t\overline\eta_0+ \frac{3\sqrt{gh_0}}2\partial_x\overline\eta_0\overline\eta_0+\frac{ (gh_0)^{3/2}}6\partial_{x}^3\overline\eta_0=0
\end{eqnarray*}

and thus, recalling (<ref>),

\begin{eqnarray*}
\partial_t\overline\eta_0+\frac{ \sqrt{g} h_0^{5/2}}{6}\partial_{x}^3\overline\eta_0
+ \frac{3h_0^{3/2}}{2\sqrt{g}}\partial_x\overline\eta_0\overline\eta_0+\sqrt{gh_0} \partial_x\overline\eta_0=0.
\end{eqnarray*}

In this expression, one can get rid of the term $\partial_x\overline\eta_0$ by a vertical translation of $\overline\eta_0$,
that is, if

\begin{equation}\label{ZETATRAVEKDG32}\zeta_0(x,t):=\alpha\left(\overline\eta_0(x,t)+\frac{2g}{3h_0}\right)
=\frac{h_0^2}g \left(\overline\eta_0(x,t)+\frac{2g}{3h_0}\right),\end{equation}
leading to

\begin{eqnarray*}
\partial_t\zeta_0+\frac{ \sqrt{g} h_0^{5/2}}{6}\partial_{x}^3\zeta_0
+ \frac{3\sqrt{g}}{2\sqrt{h_0}}\partial_x\zeta_0\zeta_0=0,
\end{eqnarray*}

which corresponds[Retaking (<ref>), we know that
the velocity of the traveling wave $\phi_0$ is proportional to $\frac{c}a$ times the amplitude of the wave.
It is suggestive to compare this information with the numerology of the coefficients
and the relation with the physical shape of the wave. Namely, in this case we have that $\frac{c}a=c=\frac{3\sqrt{g}}{2\sqrt{h_0}}$
and, by (<ref>), the amplitude of the traveling wave would formally correspond to $
\frac{h_0^2}g \left(\max\overline\eta_0+\frac{2g}{3h_0}\right)=\frac{h_0^2}g \max\overline\eta_0+\frac{2h_0}{3}$.
The velocity of the traveling wave would thus be proportional to $\frac{3\sqrt{g}}{2\sqrt{h_0}}\left(
\frac{h_0^2}g \max\overline\eta_0+\frac{2h_0}{3}\right)=
\frac{3h_0^{3/2}}{2\sqrt{g}} \max\overline\eta_0+{\sqrt{gh_0}}
$. Once again, this suggests that the speed of small amplitude waves in a shallow canal is dictated by the quantity ${\sqrt{gh_0}}$,
plus an increasing term in the amplitude of the wave.

Of course, we do not aim here at exhausting the study of the
complicated topic of water waves' speed. See e.g. <cit.>
and the references therein for further details on this subject.]
to (<ref>) with $a:=1$, $b:=\frac{ \sqrt{g} h_0^{5/2}}{6}$ and $c:=\frac{3\sqrt{g}}{2\sqrt{h_0}}$.

For further information[For simplicity, we did not discuss the role of the surface tension
in the derivation of the Korteweg-de Vries equation, but we mention that this additional
term can be also easily included in the coefficient $b$ in (<ref>).] about the Korteweg-de Vries equation
and related topics see e.g. <cit.>
and the references therein.

§.§ Plasma physics

A state of matter is one of the distinct forms in which matter can exist in nature.
Classically, nature was classified
into three states of matter: solid, liquid, and gas.
Nowadays it is common[Also, a large number of
intermediate states are nowadays known to exist, such as
liquid crystals,
Bose-Einstein condensates,
Fermionic condensates,
neutron-degenerate matter, quark–gluon plasma, superfluids, superconductors, etc., but
we do not investigate these states here.

From the historical point of view,
plasma was first identified in laboratory by Sir William Crookes in 1879
(though more systematic studies of plasma only began in the twentieth century).

Besides discovering plasma, Crookes is also credited with discovering
thallium, inventing a radiometer and, quite conveniently, some $100\%$
ultraviolet blocking sunglass lens.

He was also interested in spiritualism and became president of the Society for Psychical Research,
whose purpose is to understand psychic and paranormal events.

See Figure <ref>
for a caricature of Sir William Crookes (by caricaturist Sir Leslie Ward).]
to consider also a fourth state
of matter, namely plasmaplasma.

Though our everyday experience with plasma appears to be rather limited
(we are experienced to plasma mostly only by
neon tubes, plasma lamps and lightnings),
plasma is allegedly the most abundant form of ordinary matter in the universe
(excluding dark matter and dark energy) and likely amounts to about $
99\%$ of the total mass of the visible universe.
It consists of a gas of ions and free electrons, usually produced by
very high temperatures ($10^4$ degrees or more)
which make electrons leave
their orbit around the corresponding nuclei. In this way,
if the initial status of the gas is overall in an electrical neutral state,
the plasma consists in a mixture of charged particles, namely the free electrons
and the corresponding ions.

Caricature of Sir William Crookes (Public Domain image from

The mathematical description of a plasma is a highly advanced subject but,
in a nutshell, there is a hierarchy of models accounting for the
the evolution of a plasma.

The most intuitive one would be to describe the plasma
by detecting positions and velocities of all its particles (ions and electrons) at a given time.

This approach is very precise but often impractical, given the high number of particles involved,
and therefore intermediate
models (often called “kinetic models”) have been introduced
to describe the plasma “in average” through the statistical analysis of the particle

Moreover, at a large scale, efficient models for a plasma leverage the knowledge of
fluid dynamics equations (somewhat close in the spirit to the ones that we
presented in Sections <ref> and <ref>):
namely, at a macroscopic level, close to thermodynamic equilibrium,
one can identify
each species of particles of a plasma (ions and electrons) with a fluid
and describe the corresponding density,
velocity and energy via a set of[To confirm the importance of partial differential
equations in our understanding of plasma, let us mention for instance that the
laboratory specialized in plasma and located in Toulouse, France,
is named LAboratoire PLAsma et Conversion d'Energie (LAPLACE).] partial differential equations.

To make this model concrete, we denote by an index $j\in\{1,2\}$
the species of particles of a given plasma (e.g., $j=1$ corresponding to ions and $j=2$ corresponding to
electrons). We suppose that the plasma, being neutral in average
but composed at a small scale by charged particles, is subject to its own magnetic and electric fields (denoted by $B$
and $E$, respectively).
We recall that charged particles are influenced by both the electric and the magnetic fields: more specifically,
the so-called
Lorentz force acting on a single charged particle (say, with charge $q_j$ and velocity $v_j$)
is of the form $q_j(E+v_j\times B)$. Therefore the total contribution of the Lorentz force acting on a portion of plasma with
particle density $\mu_j$ takes the form

\begin{equation}\label{TOTAFO-1}
\mu_jq_j(E+v_j\times B).\end{equation}

Additionally, the possible variations of density exert a force on the plasma similar to a pressure. Thus, we denote by
\begin{equation}\label{TOTAFO-2}
-\nabla p_j\end{equation}
this type of pressure (the minus sign
manifesting the fact that higher densities oppose the motion).

Differently from the case of classical fluids, in the description of a plasma we have also to account for
the effect of possible particle collisions. Since we are looking here at a macroscopic description
of the plasma, we suppose that the total momentum of the particles of species $j$
comes from the average over all particles of that species and therefore
like particle collisions do not change the total momentum (say, if two electrons collide,
one may slow down after collision and the other can be accelerated, but the total momentum
is preserved, thanks to Newton's Law; same if two ions collide).
We have therefore to focus on the collisions between unlike particles,
which allow momentum to be exchanged between the species:
in this situation, when an electron and a ion collide, the momentum loss of one particle entails
a corresponding momentum gain of the other particle:
that is if the first particle experiences a variation of momentum given by $P_1$,
the second particle experiences a variation of momentum given by $P_2=-P_1$.

Computing the momentum exchange of colliding particles can be in general a rather complicated task.
A simplifying assumption in this framework is to assume that the particles are reduced to a point.
An additional simplification arises if we are willing to suppose that the collisions are perfectly elastic (no dissipation
of energy due to particle collisions).
Furthermore, we can also reduce the problem to a simpler one if we are willing to use the fact that the mass
of the electrons is typically much smaller than the mass of the ions and take therefore
the mass of the ions as the leading term of a more complex approximation.
To make these ideas work, at least in this simplified setting, we recall that the elastic collision
of two point masses, say with masses $m_1$ and $m_2$,
initial (before collision) velocity $v_{1i}$ and $v_{2i}$
and final (after collision) velocity $v_{1f}$ and $v_{2f}$
is fully described by the equations

\begin{dcases}v_{1f}=\displaystyle{\frac {(m_{1}-m_{2})v_{1i}+2m_{2}v_{2i}}{m_{1}+m_{2}}},\\v_{2f}=\displaystyle{\frac {(m_{2}-m_{1})v_{2i}+2m_{1}v_{1i}}{m_{1}+m_{2}}},\end{dcases}

see e.g. <cit.>.

Hence, since we consider the index $j=1$ as corresponding to the electrons and we assume that $m_1\ll m_2$,
we have that the momentum exchange for the colliding electron is

\begin{eqnarray*}
m_1\left(\frac {(m_{1}-m_{2})v_{1i}+2m_{2}v_{2i}}{m_{1}+m_{2}}-v_{1i}\right)\\&=&
\frac {2m_1m_2(v_{2i}-v_{1i})}{m_{1}+m_{2}}\\
\end{eqnarray*}

In the framework of plasma collisions, this would give that $P_1=2\mu_1 m_1(v_2-v_1)$
and accordingly $P_2=2\mu_1 m_1(v_1-v_2)$.
Thus, if we denote by $\phi$ the collision frequency between the two species
(for simplicity, we are disregarding here the collisions between ions, or electrons,
with the remaining neutral gas),
the species corresponding to $j=1$ changes its momentum by a quantity $\phi\mu_1 m_1(v_2-v_1)$
and the species corresponding to $j=2$ changes its momentum by a quantity $\phi\mu_2 m_2(v_1-v_2)$
(where constants are omitted for the sake of simplicity).

The first durable colour photographic image (Public Domain image from Wikipedia).

From this, (<ref>) and (<ref>), by Newton's Law of momentum balance have that

\begin{eqnarray*} \mu_j m_j\partial_t v_j +\mu_j m_j(v_j\cdot\nabla)v_j&=& \mu_jm_j\frac{d}{dt} (v_j)\\&=&
\mu_jq_j(E+v_j\times B)-\nabla p_j+\phi\mu_j m_j(v_{k_j}-v_j),
\end{eqnarray*}

$$ k_j:=\begin{dcases}
1 & {\mbox{ if }}j=2,\\
2 & {\mbox{ if }}j=1.
\end{dcases}$$

This equation is usually complemented by a continuity equation (corresponding to the conservation
of the total number of particles, which can be seen as the counterpart here of the mass transport equation
in (<ref>)),
an equation of state (that is a constitutive law relating pressure and particle density,
which can be seen as the counterpart of the barotropic flow description given on page BAROFL,
usually taken of the form $p_j =\kappa \mu_j^\gamma$, for suitable positive constants $\kappa$ and $\gamma$)
and the classical Maxwell's equations[James Clerk Maxwell
founded the theory of electromagnetism, proving that electric and magnetic fields travel through space as waves moving at the speed of light and argued that light itself is an undulation causing electric and magnetic phenomena. This great conceptual unification of light and electrical phenomena made Maxwell one of the nineteenth century scientists having the greatest impact on the subsequent developments of relativity and quantum field theory.

Maxwell also worked on the kinetic theory of gases and pioneered the early days of control theory. Also, he established that the rings of Saturn were made of numerous small particles and he is credited for presenting the first durable colour photograph. For this, he had the idea of superimposing on a screen shots taken with red, green and blue filters: his result, depicting a tartan ribbon, is shown in Figure <ref>.
The photographic plates available at the time were almost insensitive to red and barely sensitive to green, hence the results in themselves were perhaps far from perfect, but Maxwell clearly knew this was after all just a rather minor detail in the Big Picture of Science
and he wrote
“by finding photographic materials more sensitive to the less refrangible rays, the representation of the colours of objects might be greatly improved”.] for the electric and magnetic fields:
omitting structural constants, this set of equations thus takes the form[Sometimes the setting in (<ref>)
is called in jargon “Euler-Maxwell equations”Euler-Maxwell equations, since it mixes the classical Euler's formulation
of fluid dynamics presented in Section <ref> (after suitable corrections
to deal with particle collision) with the Maxwell equationsMaxwell equations for electromagnetism.]

\begin{equation}\label{MAXWEUL}
\begin{dcases}
\mu_j m_j\partial_t v_j +\mu_j m_j(v_j\cdot\nabla)v_j=\mu_jq_j(E+v_j\times B)-\nabla p_j+\phi\mu_j m_j(v_k-v_j),\\
\partial_t\mu_j+\div(\mu_jv_j),\\
p_j =\kappa \mu_j^\gamma,\\
\div B=0,\\
\div E=\mu_1q_1+\mu_2q_2,\\
\curl E=-\partial_t E,\\
\curl B=\mu_1q_1v_1+\mu_2q_2v_2+\partial_t E,
\end{dcases}
\end{equation}

with $j\in\{1,2\}$.

For additional information about plasma physics, see e.g. <cit.> and the references therein.

§.§ Galaxy dynamics

Interestingly, a suitable modification of the model used to describe the motion of fluids
presented in Section <ref>
can represent the dynamics of galaxies and globular clusters
(see e.g. Figure <ref> for
a fascinating picture).
The ansatz of the model that we describe here is that the stars forming the galaxies interact only
by the gravitational field that they create collectively
and they can be described by parcels of a self-gravitating gas.
Neglecting collisional and relativistic effects, and disregarding the physical and chemical reactions
that continuously modify the internal structure of the stars, one
can provide a suitable set of equations describing the large scale galaxy motion.
These are named Jeans equationsJeans equations after Sir James Hopwood Jeans
(though were probably first derived by James Clerk Maxwell
and can be seen as a collisionless version of the Boltzmann equations
that describe the statistical behavior of a thermodynamic system not in equilibrium)
and are as follows:

\begin{equation}\label{STARS}
\begin{dcases}
\partial_t f+v\cdot\nabla_x f-\nabla u\cdot\nabla_v f=0,\\
\Delta u=\rho
\end{dcases}
\end{equation}

where $f=f(x,t,v)$ denotes the phase-space[By
phase-space here we mean the collection of position and velocity variables.]
density of stars, that is the density of stars corresponding to time $t$, position $x\in\R^3$ and velocity $v\in\R^3$,
$u=u(x,t)$ is the gravitational potential induced collectively by the stars,

$$ \rho=\rho(x,t):=\int_{\R^3}f(x,t,v)\,dv$$

is the star density corresponding to time $t$ and position $x$.

Hubble Space Telescope picture of the galaxies NGC 2207 (left) and IC 2163 (right) (Public Domain image from Wikipedia).

To understand the rationale of equation (<ref>) we argue as follows.
In analogy with the dynamics of fluid parcels described in (<ref>),
we assume that star dust trajectories follow the law $\dot x(t)=v(x(t),t)$
and that their acceleration follows the gravitational law $\ddot x(t)=-\nabla u(x(t),t)$.

One can therefore obtain an analogous of the continuity equation in (<ref>)
by assuming that, in absence of encounters, collisions
and collapses, the amount of stars is preserved by the flow.
Namely, the total quantity of stars occupying a region $Z$ of the phase-space $\R^3\times\R^3$ at a given time $t_0$
is given by the quantity $\int_{Z} f(x,t_0,v)\,dx\,dv$
and, in a small time $\tau$,
possibly occupying a different region, that we name $Z_\tau$ which collects all the evolution trajectories
in phase space given by

\begin{eqnarray*}
&& \Big(x(t_0+\tau),v\big(x(t_0+\tau),t_0+\tau\big)\Big)=
\big(x(t_0+\tau),\dot x(t_0+\tau)\big)=
(x(t_0),\dot x(t_0))+\tau (\dot x(t_0),\ddot x(t_0))+o(\tau)\\&&\qquad=
\Big(x(t_0),v\big(x(t_0),t_0\big)\Big)+\tau \Big(v\big( x(t_0),t_0\big),-\nabla u(x(t_0),t_0)\Big)+o(\tau).\end{eqnarray*}

Hence, the assumption that stars are conserved thus translates into

\begin{eqnarray*}
0&=&\left.\frac{d}{d\tau}\int_{Z_\tau} f(x,t_0+\tau,v)\,dx\,dv\right|_{\tau=0}\\
&=&\left.\frac{d}{d\tau}\left(\int_{Z} f\Big(x+\tau v(x,t_0), t_0+\tau,v-\nabla u(x,t_0)\Big)\,dx\,dv+o(\tau)\right)\right|_{\tau=0}\\
&=& \int_{Z}
\Big[ \nabla_x f(x,t_0,v)\cdot v(x,t_0)+\partial_t f(x,t_0,v)
-\nabla_v f(x,t,v)\cdot\nabla u(x,t_0)
\Big]\,dx\,dv.
\end{eqnarray*}

The arbitrariness of the phase-space domain $Z$ thus leads to

$$ \nabla_x f(x,t_0,v)\cdot v(x,t_0)+\partial_t f(x,t_0,v)
-\nabla_v f(x,t,v)\cdot\nabla u(x,t_0)=0,$$

which is the first equation in (<ref>).

Additionally, by Gauß's Flux Law for gravity,
we know that the flux of the gravitational field through a surfaces
balances the total mass enclosed within the surface: namely,
neglecting dimensional constants,
for every domain $\Omega\subset\R^3$,

$$ \int_{\partial\Omega} \nabla u(x)\cdot\nu(x)\,d{\mathcal{H}}^{n-1}_x=\int_\Omega \rho(x,t)\,dx.$$

From this and the Divergence Theorem we arrive at

$$ \int_{\Omega} \Delta u(x)\,dx=\int_\Omega \rho(x,t)\,dx,$$

and consequently, since $\Omega$ is arbitrary,
we find that $\Delta u=\rho$, which
is the second equation in (<ref>).

For additional information about galaxy dynamics see e.g. <cit.> and the references therein.

§.§ How to count what we cannot see

In 1812, Lorenzo Romano Amedeo Carlo Avogadro, Count of Quaregna and Cerreto,
hypothesized that the volume of a gas at a given pressure and temperature is proportional to the number of atoms or molecules.
That is, equal volumes of gases at the same temperature and pressure have the same number of molecules,
regardless of the nature of the gas.
In particular,
by Avogadro's Law, the number of molecules or atoms in a given volume of ideal gas is independent of their size.
While this prescription only holds true for ideal gases,
and real gases show instead small deviations from this ideal configuration,
Avogadro's Law is often a very useful approximation and it always provides a great conceptual tool,
since it detects a universal quantity that only depends on reasonable macroscopic parameters, such as
temperature and pressure, and is independent of the specific situation under consideration, namely the type of gas
that one is taking into account.
See Figure <ref> for a portrait of Avogadro (by C. Sentier).

Portrait of Amedeo Avogadro (Public Domain image from

It is perhaps worth stressing that Avogadro's Law is not so intuitive.
For instance, if one has three units of volumes of hydrogen
and one unit of volume of nitrogen, after they combine together and produce
ammonia, how many units of ammonia do we expect (assuming that
pressure and temperature are maintained constant)?
One unit? Three units? Four units?

The correct answer according to Avogadro's Law is two units of volume.
This is because three volumes of hydrogen contain $3k$ molecules of hydrogen $H_2$ (for some $k\in\N$)
and one volume of nitrogen contains $1k$ molecules of nitrogen $N_2$
(and the proportionality factor $k$ is the same for both hydrogen
and nitrogen, by Avogadro's Law).
Thus, from the reaction

\begin{equation}\label{REACCHI} 3H_2+1N_2=2NH_3,\end{equation}

the combination of $3k$ molecules of hydrogen
and $1k$ molecules of nitrogen produces $2k$ molecules
of ammonia $NH_3$. And again by Avogadro's Law
this corresponds to two volumes of ammonia. See Figure <ref>
for a sketch of this phenomenon.

Also, in chemistry, an obvious practical issue is given by the fact that molecules are small
while in the lab one has to deal with macroscopic quantities.
Thus, while it is desirable from the theoretical point of view to found the concept of
“amount of a given substance” on the number of “elementary entities”, i.e., atoms or molecules, present in the substance
(because these elementary entities are the ones which will play a role in the chemical reactions),
for practical purposes it is often convenient to relate the notion of “amount of a given substance”
to the ratio of measured macroscopic quantities (because measuring a huge number of microscopic entities
is typically unfeasible).

We remark that the equivalence between these two approaches is a notable consequence of
Avogadro's Law. Indeed, if, at a given pressure and temperature, we take, say, one gram of
atomic hydrogen $H$ (or, more simply, two grams of the common molecule of hydrogen $H_2$)
and we measure its volume, and then we consider the same volume of another gas, such as the carbon dioxide $CO_2$,
we know from Avogadro's Law that the latter contains the same number of molecules of $CO_2$
as the number of atoms of hydrogen $H$ (or the number of molecules of $H_2$) in the previous container.
If we measure the weight of the carbon dioxide $CO_2$ in the second container,
we find that it amounts to $44$ grams: we can accordingly conclude that
the mass of a molecule of carbon dioxide $CO_2$
is $44$ times bigger than the mass of an atom of hydrogen $H$
(or $22$ times bigger than the mass of a molecule of hydrogen $H_2$).
That is, we have measured, with a minimal effort, the microscopical weight of any molecule of any gases
in terms of a given “unit of measure”, such as the atomic weight of $H$ (or the molecular weight of $H_2$).

Hence, it has become a common practice to adopt
the notion of “mole” to denote the mass of a given substance
which contains the same number of molecules (or atoms in case of
pure atomic elements) as one gram of atomic hydrogen $H$.
In the previous example, since the mass of a molecule of carbon dioxide $CO_2$
happens to be $44$ times bigger than the mass of an atom of hydrogen $H$,
we infer that

\begin{equation}\label{LAMOLE}
{\mbox{a mole of carbon dioxide weights~$44$ grams.}}\end{equation}

The algebra in (<ref>)
with the corresponding gas volumes according to Avogadro's Law.

Moreover, by Avogadro's Law,
a mole of any gas at a given temperature and pressure occupies the same volume, making it a quite convenient
way to measure the (relative) amount of substance just by comparing volumes of gases.

The natural question is thus, given a certain chemical compound, how many molecules
of it correspond to a mole? How does this number depend on the specific chemical compound?
The answer is easy: by construction this number is the same for all substances and equals
the number of hydrogen atoms $H$ necessary to form one gram of hydrogen, and this
is called[We stress that
Avogadro's Number is not a “pure number” but a physical constant of dimension
“one over moles”. See also <cit.> for additional information on Avogadro's Number.]
Avogadro's Number ${\mathcal{N}}_{\!\mathcal{A}}$. Avogadro's Number

Notice that the knowledge of ${\mathcal{N}}_{\!\mathcal{A}}$ also permits to “pass from the ratios
to the real quantities”. For instance from (<ref>) we obtain that each molecule of
carbon dioxide weights $\frac{44}{ {\mathcal{N}}_{\!\mathcal{A}} }$ grams:
in general, the knowledge of Avogadro's Number entails the knowledge of atomic and molecular masses
as a byproduct of macroscopic measurements.

Therefore, it is highly desirable to have a precise[It turns
out (possibly with some mild approximation)
that ${\mathcal{N}}_{\!\mathcal{A}}=6.02214076\times 10^{23}$. Several people noticed
a similarity between this number and $2^{79}$ (which is half of the so-called yobibyte $2^{80}$):
indeed, it is a rather surprising coincidence that

$$ \frac{{\mathcal{N}}_{\!\mathcal{A}}-2^{79}}{{\mathcal{N}}_{\!\mathcal{A}}+2^{79}}$$

is almost zero.
There is nothing special in this, coincidences happen.
For instance, it is a mere, and quite remarkable, coincidence that the most important physical constant, namely the speed
of light, is almost equal to $300.000.000$ m/s.
It is also a coincidence that
the angular diameter of the Sun seen from Earth is quite close to that of the Moon
making it possible to see a full solar eclipse.
It is also a coincidence that

$$ \left(\frac{\pi^e-e^\pi}{\pi^e+e^\pi}\right)^{\pi+e}$$

turns out to be almost equal to zero.]
quantification of ${\mathcal{N}}_{\!\mathcal{A}}$. This is not a mathematical problem per se, since
to quantify ${\mathcal{N}}_{\!\mathcal{A}}$ one can only measure it through experiments.
However, one needs a brilliant piece of mathematics to contrive
a practical procedure to experimentally calculate ${\mathcal{N}}_{\!\mathcal{A}}$
and we recall here the method devised by Albert Einstein in <cit.>
(at that time, he was working as a patent clerk in Bern, Switzerland,
and the very notion of atoms and molecules was still a subject of controversy and intense scientific debate, since individual
atoms and molecules were, for the instruments available in 1905, simply
“too small” and “too fast”).

To measure Avogadro's Number, Einstein utilized a molecular theory of heat
more or less in the lines of the Brownian motion described here
in Section <ref> and leading to the heat equation in (<ref>).
For this, Einstein thought of some ideal spherical particles of radius $a$
subject to a Brownian motion in an ideal (for simplicity, one-dimensional and infinitely long) pipe
containing some fluid with viscosity coefficient equal to $\mu$.
In this setting, according to Einstein's calculation,

\begin{equation}\label{EINAVOGEAFG}
{\mathcal{N}}_{\!\mathcal{A}}=\frac{RTt}{3\pi \mu a\,X^2(t)},
\end{equation}

where $T$ denotes the temperature (or, more precisely, the thermodynamic temperature, that is the temperature measured in Kelvin), $t$ is any given time and $X^2(t)$ is the averaged square distance traveled by a Brownian particle in an interval
of time equal to $t$.

Equation (<ref>) is quite remarkable since it equates
Avogadro's Number with quantities which can be measured experimentally, thus leading
to a precise determination of ${\mathcal{N}}_{\!\mathcal{A}}$
(that's the power of mathematics!).

Jean Perrin (Public Domain image from

The experiment suggested by Einstein was indeed carried through by Jean Baptiste Perrin (see Figure <ref>)
and his team of research students
in 1909. The set up of the test consisted in a camera lucida with a microscope to
observe and record the motion of suspended gamboge[Gamboge
is a yellowish pigment.
The word gamboge comes from gambogium, the Latin word for the pigment, which in turns derives from Gambogia, the Latin word for Cambodia.

Besides its strong laxative properties, which did not play a role in Perrin's experiment,
gamboge also possessed the convenient feature of producing, with appropriate alcoholic additives and after a selective centrifuge,
almost perfect spherules with equal radius, thus providing the ideal suspension particles for the experiment
proposed by Einstein (this was not a cheap preparation,
about one kilo of gamboge and several months of work produced few decigrams of useful spherules).]
particles in
a liquid of a given viscosity and constant temperature.
Marking the particle's position on a piece of
graph paper at timed intervals, and repeating the test under different conditions
(such as different viscosity, temperature, radius of the gamboge grains, etc.), Perrin
obtained consistent results about the right hand side of (<ref>),
thus producing a rather accurate measurement of Avogadro's Number ${\mathcal{N}}_{\!\mathcal{A}}$
through a wealth of measurements that could not be contested.

This experiment also provided experimental confirmation of Einstein's equation (<ref>) and
raised atoms and molecules from the status of hypothetical objects to concrete and observable entities,
thus offering the ultimate confirmation to the atomic nature of matter
and concluding the struggle regarding the physical reality of molecules.
For these achievements, Perrin was awarded with the Nobel Prize for Physics in 1926.

It is now time to go back to Einstein's equation (<ref>)
to understand its roots in the theory of partial differential equations.
To this end, we recall that the first ingredient towards (<ref>) is Stokes' LawStokes' Law
about the viscosity force acting on a small sphere moving through a viscous fluid.
Namely, studying the incompressible steady flow of the Navier-Stokes equation (see Section <ref>), George Gabriel Stokes in 1851 had proposed the formula

\begin{equation}\label{SYO-EDCSJDKA235athcma}
F=6\pi \mu a v\end{equation}

to describe the frictional force $F$
acting on a spherical particle of radius $a$ with relative velocity $v$ in a fluid with viscous coefficient equal to $\mu$.

The second ingredient towards (<ref>) is
the van 't Hoff equationvan 't Hoff equation for[Equation (<ref>)
is named after Jacobus Henricus van 't Hoff Jr.,
first winner of the Nobel Prize in Chemistry, see Figure <ref>. Allegedly, van 't Hoff
chose to study chemistry
against the wishes of his father. Also, it seems he came up with equation (<ref>)
after a chance encounter and a quick chat with a botanist friend during a walk in a park in Amsterdam,
thinking about an analogy to the law for ideal gases.] the osmotic pressure $p$
of a solution containing $  n$  moles of solute particles in a solution of volume $  V $,
given by

\begin{equation}\label{BVNBSstVASNnsrNA}

in which $R$ is the universal gas constant.

Henry van 't Hoff (Public Domain image from

The third ingredient towards (<ref>) is the assumption that
the Brownian motion of the suspension particles can be effectively described by the heat equation,
as discussed in Section <ref>. Hence, we suppose that the density $u$ of the particles fulfills the equation

\begin{equation}\label{EINHEAT}
\partial_t u= c\partial_{xx} u,
\end{equation}

for some constant diffusion coefficient $c>0$ (here, we are assuming that the motion of the particles is confined
to a linear pipe, hence the variable $x$ is one-dimensional).

We can now retake van 't Hoff equation (<ref>)
by considering that the number of floating particles $N$ over the volume $V$
equals the density $u$ and, by definition, a mole contains ${\mathcal{N}}_{\!\mathcal{A}}$ particles.
These observations lead to

\begin{equation}\label{EDCSJDKA235S6IMfujNMimfr0}
p=\frac{NRT}{{\mathcal{N}}_{\!\mathcal{A}}\,V}=\frac{RT u}{{\mathcal{N}}_{\!\mathcal{A}}}.
\end{equation}

Also, the force $F$ acting on each floating particle arises from a pressure gradient, namely the gradient pressure force with respect to the unit of mass is $\frac{\partial_x p}u$.
This remark and (<ref>) give that

\begin{equation*}
F u=\partial_xp= \frac{RT\,\partial_x u }{{\mathcal{N}}_{\!\mathcal{A}}}.
\end{equation*}

At equilibrium, the force must correspond to the frictional force in (<ref>), therefore we obtain that

\begin{equation}\label{EDCSJDKA235S6IMfujNMimfr02}
6\pi \mu a v u=\frac{RT\,\partial_x u}{{\mathcal{N}}_{\!\mathcal{A}}}.\end{equation}

Now we compute the flux of floating particles through an ideal cross section of the pipe
by using two possible strategies. On the one hand, we can consider the model of particles traveling with velocity $v$,
hence we can quantify the flux by the product between $v$ and the particle density $u$.
On the other hand, the flux at some cross section $x$ can also be obtained by the variation in time
of the mass preceding $x$,
that is $\partial_t \int_{-\infty}^x u(\xi,t)\,d\xi$. These observations and (<ref>) give that

$$ v(x,t) u(x,t)=\partial_t \int_{-\infty}^x u(\xi,t)\,d\xi=
\int_{-\infty}^x \partial_t u(\xi,t)\,d\xi=c\int_{-\infty}^x \partial_{xx} u(\xi,t)\,d\xi=
c\partial_x u(x,t),$$

where, as usual, we have performed formal calculations and assumed some convenient decay of the solution at $-\infty$.

From this and (<ref>) we arrive at

\begin{equation}\label{EDCSJDKA235S6IMfujNMimfr03}
{\mathcal{N}}_{\!\mathcal{A}}=\frac{RT\,\partial_x u}{6\pi \mu a v u}=\frac{RT}{6\pi \mu a c}

Now we observe that

$$g (x,t):=\frac {1}{\sqrt {4\pi ct}}\exp \left(-{\frac {x^{2}}{4ct}}\right)$$

is a solution of (<ref>) concentrating at the origin when $t=0$, hence
we can suppose that the general solution $u$ of (<ref>) occurs as a superposition[The formalization
of this superposition method is the core of the notion of fundamental solution and, in the elliptic setting, it will
play a decisive role in Section <ref> and in particular in the forthcoming
Proposition <ref>.]
of traslations of $g$ of the form $g(x-x_0,t)$ (notice that each of these translations
would correspond to a particle located at $x_0$ when $t=0$).

The average squared displacement at time $t$ of the particle starting at the origin is therefore provided by the quantity

$$ \int_\R x^2 g(x,t)\,dx=\int_\R\frac {x^2}{\sqrt {4\pi ct}}\exp \left(-{\frac {x^{2}}{4ct}}\right)\,dx=2ct.

Since this displacement should not depend on the original position, we infer that it is the same for all
possibilities of initial location of the particles, namely

\begin{equation}\label{EINAVOGEAFG2} X^2(t)=2ct.\end{equation}

This is conceptually an important result since it suggests that velocities (namely ratios of $X(t)$ over $t$)
of Brownian motions behave sort of $\frac{\sqrt{2ct}}{t}=\frac{\sqrt{2c}}{\sqrt{t}}$,
hence quite irregularly[For a formalization of this heuristic and imprecise discussion see e.g. <cit.>.] for short times. That is,
measuring the velocity of the floating particles in an extremely short interval of time
just produces a result that approaches infinity (in a sense,
attempted experiments in this direction would simply end up measuring the wrong quantity).
One of the merit of Einstein's approach is thus to have bypassed this hindrance by detecting
the relevant physical quantities which were both consistent from the mathematical point of view and experimentally
measurable in laboratories. Notice indeed that by plugging (<ref>)
into (<ref>) we obtain (<ref>), as desired.

§.§ Vibrating strings

Another classical occurrence for partial differential equations arises in the description
of vibrating stringsvibrating string (e.g., the string of a guitar, see
Figure <ref>). Different models
are possible, taking into account different physical assumptions (see e.g. <cit.>),
but we focus here on the simplest possible (indeed linear) equation.
In our description, the string at rest is modeled by the segment $[0,L]\times\{0\}$
for some $L>0$ which represents the length of the string in the absence of further external forces.
One can assume that the string is constrained at the extrema
and can be deformed into the graph of a function $u$: namely, the position
of the vibrating string at each moment of time $t$ is described by the graph $(x,u(x,t))\in[0,L]\times\R$
and the constraints for the string correspond to the boundary prescription $u(0)=u(L)=0$.

Jimi Hendrix (photo by Steve Banks, image from
licensed under the Creative Commons Attribution-Share Alike 4.0 International license).

We suppose that the string is subject to gravity (acting downwards in the vertical direction)
and to its own “tension”tension (that is the internal force of the string acting between its elements
providing an elastic reaction to the external forces). The magnitude of this tension
at each point of the string $(x,u(x,t))$ will be denoted by $T(x,t)$.

We also consider the unit tangent vector to the string, as given by

$$ \tau(x,t):=\frac{\big(1,\partial_x u(x,t)\big)}{\sqrt{1+\big(\partial_x u(x,t)\big)^2}}.$$

The vibrating string model then consists in the assumption that, given a small $\delta>0$,
the tension forces acting at the string point $(x,u(x,t))$ are the byproduct of
a force $F_-$ acting on $(x-\delta,u(x-\delta,t))$ and a force $F_+$ acting on $(x+\delta,u(x+\delta,t))$
whose magnitude is proportional to the tension $T(x,t)$ and opposite tangential directions
(namely, the direction of $F_-$ is $-\tau(x-\delta,t)$ and
the direction of $F_+$ is $\tau(x+\delta,t)$), see Figure <ref>.

Therefore, the total tension force acting on the string at the point $(x,u(x,t))$ takes the form

\begin{equation}\label{TOTFOSTR}\begin{split}&F_+-F_-=\kappa T(x)\big(\tau(x+\delta,t)-\tau(x-\delta,t)\big)=
2\kappa \delta T(x)\partial_x\tau(x,t)+o(\delta)\\&\qquad\qquad\qquad=
2\kappa \delta T(x)\partial_x\left(
\frac{\big(1,\partial_x u(x,t)\big)}{\sqrt{1+\big(\partial_x u(x,t)\big)^2}}\right)

The constant $\kappa>0$ depends on the material of which the string is made
and accounts for the elastic properties of the string.

The total force acting on the string, taking into account gravity, is therefore

\begin{equation*} F=(F_1,F_2):=F_+-F_- -(0,m_\delta g)=
2\kappa\delta T(x)\partial_x\left(
\frac{\big(1,\partial_x u(x,t)\big)}{\sqrt{1+\big(\partial_x u(x,t)\big)^2}}\right)-(0,m_\delta g)

where $m_\delta$ is the mass of the string located between $x-\delta$ and $x+\delta$ and $g$ denotes
the gravity acceleration.

If we suppose that the density of the string is constant,  say equal to some $\rho>0$, the mass $m_\delta$
is the product of $\rho$ and the length of the string located between $x-\delta$ and $x+\delta$, namely

\begin{equation}\label{MEPSI} m_\delta= \rho \int_{x-\delta}^{x+\delta} \sqrt{1+\big(\partial_x u(\zeta,t)\big)^2}\,d\zeta
=2\rho\delta \sqrt{1+\big(\partial_x u(x,t)\big)^2}+o(\delta).\end{equation}

Furthermore, by Newton's Second Law, the vertical component of the force is equal to the mass times
the vertical acceleration of the string, that is

$$ F_2=m_\delta \partial_{tt} u(x,t).$$

In light of these considerations, we find that[From the geometric point
of view, it is interesting to observe that the term $\partial_x\left(
\frac{ \partial_x u}{\sqrt{1+\big(\partial_x u\big)^2}}\right)$ in (<ref>)
detects the curvature of the string
(remarks of this type will be formalized in a general setting
in Theorem <ref>).]
for small $\delta$

\begin{equation}\label{CUSTRJMSONGJMSD}
\begin{split}
\partial_{tt} u\,&= \frac{F_2}{m_\delta}
\\&=\frac1{m_\delta}\left[
2\kappa\delta T\partial_x\left(
\frac{ \partial_x u}{\sqrt{1+\big(\partial_x u\big)^2}}\right)- m_\delta g
\right]\\&=
\frac{2\kappa\delta T}{m_\delta}
\partial_x\left(
\frac{ \partial_x u}{\sqrt{1+\big(\partial_x u\big)^2}}\right)- g\\&=
\frac{\kappa T}{
\rho \sqrt{1+\big(\partial_x u(x,t)\big)^2}+o(1)}
\partial_x\left(
\frac{ \partial_x u}{\sqrt{1+\big(\partial_x u\big)^2}}\right)- g\\
\frac{\kappa T}{
\rho \sqrt{1+\big(\partial_x u(x,t)\big)^2}}
\partial_x\left(
\frac{ \partial_x u}{\sqrt{1+\big(\partial_x u\big)^2}}\right)- g+o(1).
\end{split}\end{equation}

By formally sending $\delta\searrow0$ we thus obtain the partial differential equation

\begin{equation}\label{CUSTRJMSONGJMSD-2}
\partial_{tt} u=
\frac{\kappa T}{
\rho \sqrt{1+\big(\partial_x u(x,t)\big)^2}}
\partial_x\left(
\frac{ \partial_x u}{\sqrt{1+\big(\partial_x u\big)^2}}\right)- g.
\end{equation}

Several different models can be consider to describe explicitly the tension $T$
and thus find an expression for (<ref>)
solely depending on the shape of the string, on its density and elastic constant
and on the gravity acceleration.
One commonly accepted model is to consider the curve as “inextensible”
and $T$ as a constant (the tension being distributed uniformly along the whole string).
Another possibility is to consider the string as a superposition of infinitesimal elastic springs and thus take $T$ as proportional to the string's elongation.
Both these models however share the common treat that for small elongations $T$
is constant. Therefore, in the small elongation approximation, we take $T:=1$
and we disregard quadratic terms in $u$ in (<ref>)
(the ansatz being that for small $u$ the linear terms
will prevail against the quadratic ones). With this simplification in mind,
one can reduce (<ref>) to

\begin{equation}\label{CUSTRJMSONGJMSD-3}
\partial_{tt} u=
\frac{\kappa }{
\rho }
\partial_{xx}u- g.
\end{equation}

In particular, in the absence of gravity, the linear approximation of the vibrating string can be seen
as a particular case of the wave equation in (<ref>).
Interestingly, the speed of propagation $c$ in (<ref>)
corresponds here to $\sqrt{\frac{\kappa }{
\rho }}$. Namely the higher the elastic parameter of the string, the higher the
speed of propagation; the higher the density (equivalently, the heavier the string),
the slower the speed of propagation.

A string subject to its own tension.

The dependence of the solution of (<ref>) upon the structural parameters
of the string is not a mere mathematical curiosity: on the contrary it has a deep impact on music,
since these parameters precisely determine the pitch of the sound that the string generates.
Let us see, for instance, how understanding partial differential equations may turn out
to be useful when tuning a guitar. For accomplishing this goal, we neglect the gravity effect and summarize (<ref>) and the boundary conditions of the string
into the following mathematical setting:

\begin{equation}\label{CwUSTRJMSONGJMSD-4}
\begin{dcases}
\partial_{tt} u(x,t)=\displaystyle{\frac{\kappa }{
\rho }}\partial_{xx}u(x,t)\qquad{\mbox{for all }}(x,t)\in(0,L)\times(0,+\infty),\\
u(x,0)=\eta\displaystyle\sin\frac{\pi x}{L},\\
\partial_t u(x,0)=0,\\
u(0,t)=u(L,t)=0\qquad{\mbox{for all }}t\in(0,+\infty).
\end{dcases}
\end{equation}

The prescriptions in (<ref>) model a guitar string
with rest length $L$ and null initial velocity which is initially displaced (say by
an expert fingerpicking) as a sinusoidal graph (for the sake of simplicity, though more complicated
initial situations could be taken into account). The small parameter $\eta>0$
has been introduced here above just to be consistent with the small elongation approximation.

It can be readily checked that the function

\begin{equation}\label{CwUSTRJMSONGJMSD-5} u(x,t)=\eta\cos\frac{\sqrt\kappa\,\pi t}{\sqrt\rho\,L}
\sin\frac{\pi x}{L}\end{equation}

solves (<ref>).
This solution is quite telling regarding the pitch played by the string.
Indeed, in light of (<ref>),
the string vibrates with a frequency

\begin{equation}\label{CwUSTRJMSONGJMSD-6}
\omega:=
\frac{\sqrt\kappa\,\pi}{\sqrt\rho\,L}.\end{equation}
In particular:

* Strings with a higher elastic coefficient produce a higher pitch (and we already know this by experience, because tightening the tuning pegs of a guitar enhances the string's tension and correspondingly raises the pitch),

* Longer strings produce lower pitches and shorter strings produce higher pitches (and we already know this by experience too, since, in playing, one changes the length of the vibrating string by holding it firmly against the fingerboard with a finger and shortening the string, that is stopping it on a higher fret, gives higher pitch),

* Strings with higher density produce lower pitches, that is more massive strings vibrate more slowly (and this is also in agreement with our experience since, for instance, on classical guitars, the low density nylon strings are used for the high pitches and the higher density wire-wound strings are employed for low pitches).
To appreciate even more the impact of mathematics on guitar playing,
we recall that an octave[The name “octave” comes from the Latin adjective “octava”, meaning “eighth”. This etymology however is possibly a bit confusing, since it only represents the interval between one musical pitch and another with double its frequency: hence, an octave is more related to the number $2$ than to the number $8$.
Also, in the common western musical composition, the octave is composed by $7$ notes, or $12$ semitones, hence the octave seems to relate more to the numbers $7$ and $12$ rather than $8$. The fact is that when talking about intervals of notes the tradition is to count from the first note to the final note included (which makes $8$ notes, somewhat justifying the name of octave).]octave
is the distance between one pitch and another with half or double its frequency:
thus, from (<ref>),
holding the string against the fingerboard halving its length has the effect of doubling the frequency of the root note (consistently with the fact that each fret on the guitar's fretboard corresponds[The choice of
dividing the octave in twelve semitonessemitone is also mathematically well grounded, since it corresponds
to frequency ratios of the form $2^{\frac{j}{12}}$ for $j\in\{0,\dots,12\}$.
On the one hand, these are irrational numbers (except for $j=0$ and $j=12$).
On the other, notes sound harmonious to our ear if the frequency of the notes is close to a simple interval:
for instance a frequency ratios such as $\frac32$ (the “perfect fifth”
in musical jargon), $\frac43$ (the “fourth”) and $\frac54$ (the “major third”).
These two observations seem contradictory, till we realize
how “close to rationals” the irrational numbers of the form $2^{\frac{j}{12}}$ are
(for $j\in\{1,\dots,11\}$ these numbers are “almost” $\frac{16}{15}$,
$\frac98$, $\frac65$, $\frac54$, $\frac43$, $\frac75$, $\frac32$, $\frac85$, $\frac53$, $\frac{16}9$
and $\frac{15}8$).
For example,

$$ \frac32-2^{\frac{7}{12}}=1.5-1.49830707688...=0.00169292312...$$

showing how close the seventh semitone is to the perfect fifth.

While other tuning scales are possible and are indeed employed in several contexts,
the equal temperament one (i.e., the one splitting the scale into equal intervals)
presents several advantages, such as the possibility of transposing a tune into other keys
thus allowing instruments to play in all keys
(or singers to sing a song in a key which is more congenial to their voice)
with minimal flaws in intonation.
The twelve-tone equal temperament based on
the division of the octave into twelve equally spaced parts
on a logarithmic scale appears the most widespread system in music today.
Likely, the common adoption of this specific temperament
also influenced the composition of music in order to accommodate the system and minimize dissonance
coming from irrational approximation of rational ratios.]
to one semitone,
the distance of twelve semitones is the octave, and the position of the twelfth fret
corresponds to a string half as long as the original one).

Gustav Kirchhoff
(Public Domain image from
Wikipedia, source Smithsonian Libraries).

For completeness, let us also mention a slightly different model for vibrating strings proposed[Besides the nonlocal equations for vibrating strings related to (<ref>),
Gustav Kirchhoff provided fundamental contributions
to the theory of electrical circuits, spectroscopy, and black-body radiation.
Which creates quite a confusion, because the names of
Kirchhoff's Law and Kirchhoff equation end up being used in all these topics, with different meanings.

By the way, Kirchhoff's home town was Königsberg, a
historic Prussian city that is now Kaliningrad, in Russia's Kaliningrad Oblast
(a small exclave of the vast Russian state
providing an interesting counterexample to the conjecture that political states are connected
regions, up to additional islands). The reason for which we mention Königsberg is because
one classical and very famous mathematical problem that was settled by Euler is related to the bridges of
Königsberg (see e.g. <cit.> to know more about the Königsberg bridges).] by Gustav Robert Kirchhoff in <cit.>, see Figure <ref>.
The gist of this model is that the elastic coefficient $\kappa$ in (<ref>)
is assumed to be a positive constant but, in reality, the elastic properties of the string may
depend on its elongation (according to the material, the string could either become “stiff”
or “slack off” for large
elongations). To account for such a possibility, Kirchhoff proposed that $\kappa$
instead depends on the length of the deformed string, namely on

\begin{equation}\label{leng2tht346y356ef7o8om-ely}
{\mbox{\footnotesize{\calligra{L}}}}\;\,(t):=\int_0^L \sqrt{1+\big(\partial_x u(x,t)\big)^2}\,dt.\end{equation}

In this model, $\kappa$ in (<ref>) must be interpreted as $\kappa\left(
{\mbox{\footnotesize{\calligra{L}}}}\;\,(t)\right)$, that is a function of the length of the string
(the corresponding equation is thus called Kirchhoff equation). Kirchhoff equation
Notice that not only in this situation $\kappa$ varies with time, but also
it depends on a “nonlocal” quantity. Indeed, the length of the deformed string in (<ref>) is a “global” object: to calculate it, it is not sufficient to
know the shape of the string at a given point, instead full information about its large scale geometry
is required. And, of course, mathematical problems requiring the knowledge of nonlocal
quantities become structurally harder (but often quite interesting, since they aim at
capturing the Big Picture).

§.§ Elastic membranes

The description of elastic membraneselastic membrane can be seen as a multi-dimensional version
of the one of vibrating strings presented in Section <ref>.
A simple model could be that of considering the membrane as constituted by
a web of infinitesimal strings in the coordinate directions.
In this way, one may think that the vertical force exerted
at a point $(x,u(x))\in\R^n\times\R$ of the membrane
is the sum of the forces produced by the tension of the infinitesimal strings
in each direction, plus gravity.
Recalling (<ref>), we may think that, for each $i\in\{1,\dots,n\}$,
the vertical force exerted by the tension of the infinitesimal string
located in direction $e_i$ (say, located from $x-\delta e_i$ and $x+\delta e_i$) is equal to

$$ 2\kappa \delta T(x)\partial_i\left(
\frac{\partial_i u(x,t)}{\sqrt{1+\big(\partial_i u(x,t)\big)^2}}\right),$$

up to higher order in $\delta$.

Hence the total vertical force acting on the membrane at $(x,u(x))\in\R^n\times\R$
is given by

$$ 2\kappa \delta T(x)\sum_{i=1}^n\partial_i\left(
\frac{\partial_i u(x,t)}{\sqrt{1+\big(\partial_i u(x,t)\big)^2}}\right)-m_\delta g,$$

where, in view of (<ref>),

=2\rho\delta \sum_{i=1}^n\sqrt{1+\big(\partial_i u(x,t)\big)^2},$$

up to higher orders in $\delta$, being $\rho$ the linear densities of the infinitesimal strings in each direction
(that we suppose to be constant and independent of the direction).
In this setting, Newton's Second Law gives that

\begin{equation}\label{CUSTRJMSONGJMSD-287729ou35t4jrg}
\begin{split}
\partial_{tt}u&=\frac{2\kappa \delta T}{m_\delta}\sum_{i=1}^n\partial_i\left(
\frac{\partial_i u}{\sqrt{1+(\partial_i u)^2}}\right)-g\\
&=\frac{\kappa T}{\displaystyle
\rho\sum_{i=1}^n\sqrt{1+(\partial_i u)^2}}
\sum_{i=1}^n\partial_i\left(
\frac{\partial_i u}{\sqrt{1+(\partial_i u)^2}}\right)-g,
\end{split}\end{equation}

which can be considered the higher dimensional[From the geometric point
of view, it is interesting to observe that the term $\displaystyle
\sum_{i=1}^n\partial_i\left(
\frac{\partial_i u}{\sqrt{1+(\partial_i u)^2}}\right)$ in (<ref>)
detects the mean curvature of the membrane
(remarks of this type will be formalized in a general setting
in Theorem <ref>).]
version of (<ref>).
In the small elongation approximation we thus consider $T:=1$ and disregard the terms that
are quadratic in the displacement $u$, thus reducing (<ref>)

\begin{equation}\label{CUSTRJMSONGJMSD-287729ou35t4jrg24qwerty} \partial_{tt}u=\frac{ \kappa }\rho\Delta u-g,\end{equation}

which is the higher dimensional version of (<ref>).

In particular, an elastic membrane in equilibrium is a stationary solution of (<ref>) and thus satisfies