text:
\begin{document}

\title[Boundary effects on the magnetic Hamiltonian dynamics]{Boundary
  effects on the magnetic Hamiltonian dynamics in two dimensions}
\author[D. T. Nguyen, N. Raymond, S. V\~u Ng\d{o}c]{Th\d{o} Nguyen
  Duc, Nicolas Raymond, San V\~u Ng\d{o}c}
\keywords{Magnetic Hamiltonian, dynamics, confinement,  scattering, boundary}
\subjclass[2010]{70H05, 37N05}

\maketitle

\begin{abstract}
  We study the Hamiltonian dynamics of a charged particle submitted to
  a pure magnetic field in a two-dimensional domain. We provide
  conditions on the magnetic field in a neighbourhood of the boundary
  to ensure the confinement of the particle. We also prove a formula
  for the scattering angle in the case of radial magnetic fields.
\end{abstract}

\section{Introduction}
\subsection{Magnetic Hamiltonian dynamics}
This article is concerned with the dynamics of a charged particle in a
smooth bounded domain $\Omega\subset\RR^2$ in the presence of a non
homogeneous magnetic field $\textbf{B}$. The motion of a particle of
charge $e$ and mass $m$ under the action of the Lorentz force can be
expressed by Newton's equation
\begin{equation}\label{eq.Newton}
m\ddot{q}= e \dot{q}\times \textbf{B}\,,
\end{equation}
where $q =(q_1, q_2,q_3)^T \in \mathbb{R}^3$. To simplify our
discussion, we assume that $e=1$ and $m=1$. The vector field
$\textbf{B}$, defined on $\Omega$, is assumed to be smooth and to
satisfy the Maxwell equation $\nabla\cdot\textbf{B}=0$. For our target
problem in two dimensions, we suppose that $\textbf{B}$ is
perpendicular to the plane $\RR^2$, \textit{i.e.},
$\textbf{B}(q)=(0,0,b(q))$. This assumption forces particles lying in
the $\RR^2$ plane and whose initial velocities are in the plane to
stay in this same plane for all time. Since a vector field in
$\mathbb{R}^3$ can be identified with a 2-form, we write the magnetic
field as $\textbf{B} = b(q) \dd q_1\wedge \dd q_2$. Then, if there is
a 1-form $\textbf{A} =A_1 \dd q_1 + A_2 \dd q_2$ such that
$\dd \textbf{A} = \textbf{B}$, we can write~\eqref{eq.Newton} in
Hamiltonian form.  Consider, for all
$(q,p)\in \mathbb{R}^2 \times \mathbb{R}^2$,
\begin{equation}\label{Ham in standard}
\mathcal{H}(q,p) = \frac{\|p-\textbf{A}(q)\|^2}{2} \h\,,
\end{equation}
where $\|.\|$ denotes the Euclidean norm on $\mathbb{R}^2$.

The matrix representing the right cross product with $\mathbf{B}$ in
the canonical basis is
\[M_{\mathbf{B}}= J^T_{\textbf{A}}-J_{\textbf{A}}\,,\]
where $J_{\mathbf{A}}$ is the Jacobian matrix of $\mathbf{A}$. Hence
Newton's equation~\eqref{eq.Newton} becomes
\[\ddot{q}=M_{\mathbf{B}}\dot{q}\,,\]
so that
\begin{equation*}
\frac{d}{dt}\left( \dot{q}+\textbf{A}(q) \right) = J_{\textbf{A}}^T \dot{q}\,.
\end{equation*}
By introducing the momentum variable $p = \dot{q}+\textbf{A}(q) $, we
see that $\mathcal{H}(q,p)=\frac{1}{2}\norm{\dot q}^2$ is the kinetic
energy of the system, and $(q,p)$ evolves according to the Hamiltonian
flow associated with $\mathcal{H}$:
\begin{equation}\label{eq.Hamilonian}
\left\{\begin{split}
\dot{q}&=\partial_{p}\mathcal{H}(q,p)\\
\dot{p}&=-\partial_{q}\mathcal{H}(q,p)
\end{split}\right..
\end{equation}
We shall always assume that $q\mapsto b(q)$ is locally
Lipschitz-continuous, ensuring that the system~\eqref{eq.Hamilonian}
has a unique local maximal solution, thanks to the Cauchy-Lipschitz
theorem. Then, the vector potential $\mathbf{A}$ will always be chosen
to be $C^1$-smooth.

\subsection{Two questions}

From now on, we call $b$ the magnetic field and it is identified with the $2$-form
\[b(q_1,q_2) \dd q_1 \wedge \dd q_2=\dd \left(A_{1}\dd q_{1}+A_{2}\dd q_{2}\right)\,.\]
This article addresses two classical dynamical problems: confinement and scattering.
\begin{enumerate}[-]
\item (Confinement) Consider a charged particle in the magnetized region $\Omega$. A natural question is the following:

\centerline{\enquote{Will the particle reach the boundary in finite time?}}
\noindent We will provide a precise answer to this question, depending
on the behaviour of the magnetic field at the boundary and on the
initial conditions. Our results will improve recent results by Martins
in \cite{GM17}. In particular, we will see that, even if the magnetic
field is infinite at the boundary, some trajectories can escape from
$\Omega$. This kind of (open) problems is mentioned in \cite[Section
1.4]{VT11}.
\item (Scattering) Consider a charged particle outside the magnetized
  region $\Omega$. Before it reaches the region $\Omega$, the
  trajectory is a straight line. If it enters the region $\Omega$,
  does the particle escape from it in finite time? And, if it does so,
  what is the deviation angle between the ingoing and outgoing
  directions? We will explicitly answer these questions in the case of
  radial magnetic fields and when $\Omega$ is a disc. In this case,
  the angular momentum commutes with the Hamiltonian and allows a
  reduction to a one degree of freedom system.
\end{enumerate}
For both problems, we provide numerical illustrations of our results.

These questions have intrinsic physical motivations. Their answers
allow a better understanding of the classical dynamics of charged
particles in magnetic fields.  The description of the classical
trajectories has also many applications, for instance, at the quantum
level. The quantum aspect of the trapped trajectories can be related
to the essentially self-adjoint character of the magnetic Laplacian
(see \cite{VT11, GI09, GI11, RS75}). It is also a key point to
describe the spectrum/resonances of magnetic Laplacians. As far as the
authors know, whereas the description of the magnetic dynamics has
allowed to estimate the spectrum of magnetic Laplacians (see
\cite{RVN15, HKRVN16}), no result seems to exist to estimate their
resonances near the real axis. Investigating the trapped trajectories
is a necessary step in this direction.

In the regime of large magnetic field and small energy, a special
treatment of the confinement problem can be done and takes advantage
of the near-integrable structure of the Hamiltonian dynamics, either
via Birkhoff normal form~\cite{RVN15}, or KAM
theorems~\cite{castilho01}. On the contrary, our results here will
give more explicit initial conditions and allow regimes where the
guiding center motion is not necessarily meaningful.

\subsection{Organization of the article}
The article is organized as follows. In Section \ref{sec.main-results}, we state our main results about confinement and scattering. Section \ref{sec.proofs} is devoted to the proofs.

\section{Statements}\label{sec.main-results}

\subsection{Confinement problem}

\subsubsection{Tubular coordinates}
In order to state our results, it is convenient to introduce tubular coordinates near the boundary of $\Omega$, following the analysis of \cite{GM17}.

We assume that the connected components of $\partial\Omega$ are
$C^2$-smooth closed curves without self-intersections. Let
$\mathcal{C}$ be a connected component of $\partial \Omega$. It can be
parametrized by its arc length
$\gamma : \mathbb{R}/L\mathbb{Z} \to \mathcal{C}$ where $L$ is the
length of $\mathcal{C}$.

There exists $\delta>0$ such that
\begin{equation}\label{eq.psi}
\psi :\left\{
\begin{split}
  (0,\delta)\times \mathbb{R}/L\mathbb{Z} &\to\Omega_{\mathcal{C}}(\delta)\\
  (n,s)\qquad &\mapsto\gamma(s)+nN(s)=q
\end{split}\right.
\end{equation}
is a smooth diffeomorphism. $N(s)$ denotes the inward pointing normal at $\gamma(s)$ and
\[\Omega_{\mathcal{C}}(\delta)=\{q\in\Omega : \dd \left(x,\mathcal{C}\right)<\delta\}\,.\]
Note that
\begin{equation}\label{eq.Bns}
\textbf{B} = b(q) \dd q_1\wedge \dd q_2 = b(\psi(n,s))(1-n\kappa(s)) \dd s\wedge \dd n\,,
\end{equation}
where $\kappa(s)$ is the signed curvature of $\mathcal{C}$ at $\gamma(s)$.
In these coordinates, we can write
\begin{equation*}
\textbf{A} = A_n(n,s) \dd n + A_s(n,s) \dd s
\end{equation*}
with $A_n , A_s$ defined on $(0,\delta)\times \mathbb{R}/L\mathbb{Z}$ such that
\begin{equation}
  \frac{\partial A_s}{\partial n} - \frac{\partial A_n}{\partial s} =:
  B(n,s)=- b(\psi(n,s))(1-n\kappa(s))\,.
\end{equation}
Via the tubular coordinates, we can define the symplectic change of coordinates
\begin{equation}\label{eq.Psi}
\Psi :\left\{
\begin{split}
(0,\delta)\times\mathbb{R}/L\mathbb{Z}\times\mathbb{R}^2&\to\Omega_{\mathcal{C}}(\delta)\times\mathbb{R}^2\\
(n,s, p_{n}, p_{s})&\mapsto (\psi(n,s),((\dd  \psi)_{(n,s)}^{-1})^\mathrm{T}(p_{n}, p_{s}))=(q,p)
\end{split}\right.,
\end{equation}
where we have explicitly $p=(1-n\kappa(s))^{-1}p_{s}\gamma'(s)+p_{n}N(s)$.

The Hamiltonian takes the form (see Lemma \ref{lem.cartesian.to.normal}):
\begin{equation}\label{eq.Htubular}
H(n,s,p_n,p_s) = \frac{1}{2}(p_n - A_n(n,s) )^2 + \frac{(p_s-A_s(n,s))^2}{2(1-\kappa(s)n)^2}\,.
\end{equation}

\subsubsection{General confinement theorems}
We can now state our confinement results.

Our first theorem provides a sufficient condition on $\mathbf{B}$ so that no trajectory can escape from $\Omega$.
\begin{theorem}\label{theo.confinement0}
For every connected component $\mathcal{C}$ of $\partial \Omega$, we assume that
\begin{equation} \label{eq.lim.of.integral}
\lim_{n \to 0} \left\lvert \int_{n}^{\delta_{\mathcal{C}}}\int_{0}^{L_{\mathcal{C}}} B(\eta,\xi)\dd \xi \dd \eta\right\rvert  =+ \infty\,,
\end{equation}
and that there exists $M_{\mathcal{C}}\geq 0$ such that, for all $(n,s) \in (0,\delta_{\mathcal{C}})\times \mathbb{R}/L_{\mathcal{C}}\mathbb{Z}$,
\begin{equation}\label{boundedness}
\left| B(n,s)-\frac{1}{L_\mathcal{C}}\int_{0}^{L_\mathcal{C}} B(n,\xi) \dd \xi \right| \leq M_{\mathcal{C}}\,.
\end{equation}
Then the magnetic Hamiltonian dynamics is complete (i.e. no solution
of \eqref{eq.Hamilonian}, starting in $\Omega$, reaches
$\partial\Omega$ in finite time).
\end{theorem}

Of course, given a starting point $q\in\Omega$, only the components
$\mathcal{C}$ that bound the connected component of $q$ in $\Omega$
need to be taken into account. Actually, there is a more quantitative
version of the previous theorem.
\begin{theorem}\label{theo.confinement1}
Consider a connected component $\mathcal{C}$ of $\partial\Omega$. Let
\[K = \sup_{s \in \mathbb{R}/L\mathbb{Z}} \vert \kappa(s)\vert\,,\qquad K' = \sup_{s \in\mathbb{R}/L\mathbb{Z}} \vert\kappa'(s) \vert\,.\]
We assume that, for some $\epsilon\in(0,1)$, $\delta$ satisfies $0<\delta \leq \epsilon/K$. We assume that there exists $M\geq 0$ such that, for all $(n,s) \in (0,\delta)\times\mathbb{R}/L\mathbb{Z}$,
\begin{equation}\label{boundness}
\left| B(n,s)-\frac{1}{L}\int_{0}^{L} B(n,\xi) \dd \xi \right| \leq M\,.
\end{equation}
Consider $T>0$ and $q(t) = \psi(n(t),s(t)) $ a trajectory contained in
$\Omega_{\mathcal{C}}(\delta)$ for $t\in [0,T]$ with energy $H_0$. Let
\begin{equation}\label{eq.fn}
f(n) = -\frac{1}{L}\int_{n}^{\delta}\int_{0}^{L} B(\eta,\xi)\dd \xi \dd \eta
\end{equation}
and assume that
\begin{equation} \label{eq:f very large}
\liminf_{n \to 0}  |f(n)|  > C(T)\,,
\end{equation}
where
\begin{multline*}
C(T) = \left| \dot{s}(0)[1-\kappa(s(0))n(0)] + \int_{n(0)}^{\delta}\int_{0}^{L} B(\eta,\xi)\dd \xi \dd \eta \right|\\
+\sqrt{2H_0}(1+\epsilon) +\left( M\sqrt{2H_0} + \frac{2H_0 K' N}{1-\epsilon} \right)  T\,.
\end{multline*}
Let $g$\footnote{such a function g always exists.} be a continuous and strictly decreasing function  such that
\[\lim_{n\to0} g(n)=\liminf_{n\to0}|f(n)|\,,\qquad g\leq |f|\,\quad \mathrm{ on }\,\, [0,\delta]\,.\]
Then, $g$ takes the value $C(T)$ and, for all $t\in[0,T)$,
\begin{equation}
n(t)> g^{-1}(C(T))\,.
\end{equation}
\end{theorem}

\begin{remark}
  Theorems \ref{theo.confinement0} and \ref{theo.confinement1} are
  improvements of \cite[Theorems 1\&2]{GM17}. They tell us that a
  particle in $\Omega$ never reaches the boundary of $\Omega$. In
  \cite{GM17}, it is assumed that $\partial_{s}B$ is integrable:
\begin{equation}\label{eq.intMartins}
\sup_{s\in\mathcal{C}}\int_{0}^N |\partial_{s}B(m,s)|\dd m<+\infty\,,
\end{equation}
and the question of removing this
assumption was explicitly mentioned as important (\emph{op. cit.,
  section 3}). Our theorems give a partially positive answer to this question,
thus allowing for magnetic fields having wilder tangential behaviors.
\begin{enumerate}[-]
\item Theorem \ref{theo.confinement0} generalizes \cite[Theorem
  1]{GM17} by replacing the integrability assumption by
  \eqref{boundedness}. This allows in particular to consider a
  magnetic field (on the unit disc) of the form
  \[B(n,s)=\frac{1}{n}+\sin\left(\frac{\chi(s)}{n}\right)\,,\]
  where $\chi$ is a smooth function supported in $(-\pi,\pi)$ such
  that $\chi'(0)\neq 0$ and $\chi(0)=0$. For this magnetic field, it
  is easy to check that \eqref{eq.intMartins} is not satisfied.  In
  fact, the $C^\infty$ smoothness is actually not required; in order to
  draw Figure \ref{fig:Pic1}, we took, for simplicity, a small
  perturbation of $\chi(s)=\arcsin(\sin(s))$.

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth ,height = 6 cm]{P09}
\caption{A trajectory obtained with a magnetic field on the unit disc
  that is strong near the boundary with a non-integrable tangential derivative:\\
  $\displaystyle B(q) =\frac{1}{1-\sqrt{q_{1}^2+q_{2}^2}} + \sin\left(\frac{\arcsin(q_2)}{1-\sqrt{q_{1}^2+q_{2}^2}} \right) + 5q_1^3 -7q_2 \,.  $}
   \label{fig:Pic1}
\end{figure}

\item An explicit lower bound for the escaping time of a magnetized region is given in \cite[Theorem 2]{GM17} in the case when
\begin{equation}\label{eq.Mh}
B(n,s) = \frac{M}{n^\alpha} + h(n,s)\,, \qquad \alpha \geq 1\,.
\end{equation}
where $M \neq 0$ and $h$ is bounded and smooth in
$\Omega_{\mathcal{C}}(\delta)$, and so that \eqref{eq.intMartins}
holds. Theorem \ref{theo.confinement1} implies \cite[Theorem 2]{GM17},
and also provides an explicit lower bound for magnetic fields that are
not in the form \eqref{eq.Mh}, see Figure \ref{fig:Con7} where the
magnetic field changes sign infinitely many times.
\end{enumerate}

\end{remark}
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth, height = 6 cm]{Con7.eps}
\caption{A trajectory obtained with a magnetic field on the unit disc
  that strongly oscillates near the boundary:\\
  $\displaystyle B(q) = \frac{\frac{1}{2}-\sin\left(
      \frac{1}{1-\sqrt{q_1^2+q_2^2}} \right)
  }{(1-\sqrt{q_1^2+q_2^2})^2} +10q_1 -2q_1^2-10q_2^2\,.  $
  } \label{fig:Con7}
\end{figure}

\subsubsection{Confinement results in the radial case}
When $\Omega=D(0,1)$ and when $B$ is radial, the dynamics is
completely integrable, and hence can be entirely described by a one
degree of freedom Hamilonian; concerning the confinement problem, this
of course leads to stronger results.

\begin{proposition}\label{prop.of.confinement2}
  Let $q(t)=(q_1(t),q_2(t))$ be a solution to \eqref{eq.Hamilonian}
  starting at $t=0$ from inside the unit disc. If the initial data
  $(q(0), \dot{q}(0))$ satisfies either~\textbf{H1} or \textbf{H2}
  below:
\begin{enumerate}
\item[\textbf{H1}:]
\begin{equation}\label{eq.criterion.1}
\liminf_{r\to 1^-}\left|\frac{1}{2\pi}\int_{\|q(0)\|\leq \|q\|\leq r} B(q)\dd q-\det(q(0),\dot{q}(0))\right|>\|\dot{q}(0)\|\,,
\end{equation}
\item[\textbf{H2}:]
\begin{equation}\label{criterion2}
\liminf_{r\to 1^-}\left|\frac{1}{2\pi}\int_{\|q(0)\|\leq \|q\|\leq r} B(q)\dd q-\det(q(0),\dot{q}(0))\right|=\|\dot{q}(0)\|\,,
\end{equation}
and
\begin{equation}\label{criterion2'}
\limsup_{r\to 1^-}\frac{\left|\frac{1}{2\pi}\int_{\|q(0)\|\leq \|q\|\leq r} B(q)\dd q-\det(q(0),\dot{q}(0)) \right|-\Vert \dot{q}(0) \Vert}{r-1} <0\,,
\end{equation}
\end{enumerate}
then the solution exists for all $t\geq 0$, and there exists $\eta \in[0,1)$ such that
\begin{equation}
\forall t\geq 0\,,\quad\Vert q(t) \Vert < \eta\,.
\end{equation}
\end{proposition}
One can find situations where none of the hypothesis of Proposition
\ref{prop.of.confinement2} hold and the trajectory can be arbitrarily
close to the boundary.
(see Figure \ref{fig.escape-turning}: this unusual behavior can be
explained by a critical point of the radial Hamiltonian at $r=1$, see
\eqref{eq.red.H}).
\begin{figure}[h]
\centering
\includegraphics[width=0.4\textwidth]{Con4}
\caption{$B(r)= e^{-r}-\frac{2}{r} $ .} \label{fig.escape-turning}
\end{figure}

If the magnetic field is $L^1$-integrable near the boundary of
$\Omega$, we can prove that there exist trajectories escaping from
$\Omega$ in finite time.  In particular, even if the magnetic field is
infinite at the boundary, the confinement is not ensured.
\begin{proposition}\label{prop.confinement.3}
When
\begin{equation}
\label{equ:L1}
\limsup_{r\to 1^-} \left\lvert \int_{D(0,r)} B(q)\dd q \right\rvert < +\infty\,,
\end{equation}
there exists a trajectory starting in $\Omega$ and reaching the boundary in finite time.
\end{proposition}

Of course, even under assumption~\eqref{equ:L1}, some trajectory may
be confined, depending on initial conditions (see Figure \ref{fig:confine} where the simulations
are performed with $B(r)=\ln^2(1-r)$).

\begin{figure}[]
\centering
\begin{subfigure}[t]{0.5\textwidth}
\centering\includegraphics[width=0.7\textwidth]{Con32.eps}
\end{subfigure}
~
\begin{subfigure}[t]{0.5\textwidth}
\centering\includegraphics[width=0.7\textwidth]{Con31.eps}
\end{subfigure}
\caption{$B(r)=\ln^2(1-r) $: the particle is confined or not.}\label{fig:confine}
\end{figure}

\subsection{Scattering in the radial case}
Let us now describe our scattering result in the radial case. We
assume that $\mathbf{B}_{|\Omega}$ admits a locally Lipschitz
extension in a neighbourhood of $\Omega$.

In polar coordinates, we have
\[
\textbf{B} = B(r)r \dd r \wedge \dd \theta=\dd \left(G(r) \dd
  \theta\right)\,,
\]
where
\begin{equation*}
G(r) = \int_0^r \tau B(\tau)  \dd \tau\,.
\end{equation*}
Via the symplectic change of coordinates
\[
\begin{split}
  \mathbb{R}_{+}^*\times \RR/2\pi\ZZ\times\mathbb{R}^2&\to (D\setminus\{0\})\times\mathbb{R}^2\\
  (r,\theta,p_{r},p_{\theta})&\mapsto \left(r\cos\theta,r\sin\theta,
    \cos \theta p_r - \frac{\sin \theta}{r}p_\theta,\sin \theta p_r
    +\frac{\cos \theta}{r}p_\theta\right)=(q,p)
\end{split}\quad\,,
\]
the Hamiltonian becomes
\begin{equation}\label{eq.Ham.in.polar}
  \tilde H(r,\theta,p_r,p_\theta) = \frac{p_r^2}{2}+ \frac{(p_\theta - G(r))^2}{2r^2}\,,
\end{equation}
In particular, the angular momentum $p_{\theta}$ is constant along the
flow and we consider the reduced one dimensional Hamiltonian on
$T^*\RR^*_+$
\begin{equation}\label{eq.red.H}
  H(r,p_r) := \frac{p_r^2}{2}+V(r)\,,\qquad V(r) :=  \frac{(p_\theta - G(r))^2}{2r^2}\,,
\end{equation}
where $V\in C^1(\RR^*_+)$.  We notice that (see, for example, Lemma
\ref{lem.cartesian.to.normal})
\[v_{r}=p_{r}\,,\qquad v_{\theta}=r^{-1}(p_{\theta}-G(r))\,,\]
where $v_{r}$ and $v_{\theta}$ are the classical radial and tangential components of the velocity $v$.

We consider a charged particle with energy $H_0$ arriving into the
disk with velocity $v_{1}$. In particular,
$H_{0}=\frac{1}{2}\|v_{1}\|^2$. If the particle escapes from the disc
with velocity $v_{2}$ (see Figure \ref{Scattering}), we have
$\norm{v_2}=\norm{v_1}$, and a natural question is to compute the
(scattering) angle between these two vectors. Let
$\omega \in (-\pi,\pi]$ be the oriented angle between $v_1$ and $v_2$.

\begin{theorem}\label{theo.scattering}
  Consider a trajectory starting on $\partial\Omega$, with velocity
  $v_{1}\neq 0$ and entering $\Omega$. This means that either
  $v_{r}<0$, or $v_{r}=0$ and $\frac{B(1)}{v_{\theta}}<-1$. We define
  $\gamma$ the angle between the outward pointing normal and $v_{1}$.

We also assume
\begin{enumerate}[\rm i.]
\item either that the equation $V(r)=H_{0}$ has a solution for $r\in(0,1)$ and that the closest solution to $1$, denoted by $r^*$, satisfies $V'(r^*)<0$.
\item or, only when $p_{\theta}=0$, that the equation $V(r)=H_{0}$ has no solution.
\end{enumerate}

Then the trajectory escapes from $\Omega$ in finite time with velocity $v_{2}$, and we can compute the scattering angle $\omega \mod 2\pi$:
\begin{enumerate}[\rm i.]
\item either the trajectory does not pass through the origin and
\[\omega=\alpha+\pi-2\gamma\,,\]
where
\begin{equation}\label{alpha angle}
\alpha = 2 \int_{r^*}^{1}  \frac{p_\theta-G(r)}{r\sqrt{2H_0 r^2- (p_\theta-G(r))^2} } \dd r\,,
\end{equation}
\item or the trajectory passes through the origin (in this case $p_{\theta}=0$) and
\[\omega=\alpha-2\gamma\,,\]
where
\begin{equation}\label{alpha angle}
\alpha = 2 \int_{0}^{1}  \frac{-G(r)}{r\sqrt{2H_0 r^2- G(r)^2} } \dd r\,.
\end{equation}
\end{enumerate}

\end{theorem}

\begin{figure}[!h]
 \centering
\begin{tikzpicture}
\coordinate (O) at (0,0);
\coordinate (A) at (1,0);
\coordinate (B) at (-4,3);
\coordinate (C) at (-4,-3);
\coordinate (D) at (-5,3.75);
\coordinate (E) at (-3.2,2);
\coordinate (F) at (-4,-3);

 \draw[red,thick] (0,0) circle (5);

 \draw[red,thick,dashed] (0,0) circle (2);

 \draw[color=black,very thick, ->] (-4,3) -- (-3.2,2);

 \draw[color=black, very thick, ->] (-4,-3) -- (-4.8,-4);

 \draw [thick,dashed] plot[smooth, tension = 0.4] coordinates{ (-4,3)  (-2,0) (-4,-3) };

 \draw[color=gray, thick, -] (-5,3.75)--(0,0);
 \draw[color=gray, thick, -] (-5,-3.75)--(0,0);

 \draw[color=gray, thick, -] (0,0)--(5,0);

 \node[scale=1] at (-3.7,2) {$\vec{v}_1$};
 \node[scale=1] at (-4,-4) {$\vec{v}_2$};
 \node[scale=1] at (1,-0.5) {$r^*$};

\pic [draw, ->,angle radius = 0.3 cm, "$\theta_1$", angle eccentricity=2] {angle = A--O--B};
\pic [draw, ->,angle radius = 1 cm, "$\theta_2$", angle eccentricity=1.4] {angle = A--O--C};
\pic [draw, -,thick, angle radius = 0.3 cm, "$\gamma$", angle eccentricity=2] {angle = D--B--E};
\pic [draw, ->,thick, angle radius = 0.5 cm, "$\alpha$", angle eccentricity=1.4] {angle = B--O--F};

\draw [decorate,decoration={brace,amplitude=10pt},xshift=0pt,yshift=0pt] (2,0) -- (0,0) node [black,midway,xshift=0.5 cm] {};
\end{tikzpicture}
\caption{ The scattering arrows. }
\label{Scattering}
\end{figure}

\section{Proofs}\label{sec.proofs}
\subsection{Proof of Theorems \ref{theo.confinement0} and \ref{theo.confinement1}}

To reach the boundary, the particle has to be close to a connected
component $\mathcal{C}$ of $\partial \Omega$. Thus, we can assume
that, for all $t\in[0,T)$,
\begin{equation*}
q(t)\in\Omega_{\mathcal{C}}(\delta)\,.
\end{equation*}
Modifying the vector potential corresponds to a symplectic
transformation of the form $(q,p)\mapsto (q, p+\dd S(q))$, for some
smooth function $S$, and hence does not modify the trajectory of the
particle. Thus, we consider the function
\begin{equation*}
  \alpha(n,s) =   \frac{s}{L} \int_{0}^{L}  B(n,\xi) \dd \xi
  -\int_0^s B(n,\xi) \dd \xi\,.
\end{equation*}
Notice that $\alpha(n,\cdot)$ is $L$-periodic. Recalling \eqref{eq.fn}
and letting $\textbf{A} = \alpha(n,s) \dd n + f(n) \dd s$, we have
$ \textbf{B}=\dd\textbf{A}$.

By \eqref{eq.Htubular}, the corresponding Hamiltonian is
\begin{equation*}
  H(n,s,p_n,p_s) = \frac{(p_n - \alpha(n,s) )^2}{2} + \frac{(p_s-f(n))^2}{2(1-\kappa(s)n)^2}\,.
\end{equation*}
Concerning Hamilton's equations, we have in particular
\begin{equation*}
\dot{n} = p_n- \alpha(n,s)\,,\qquad \dot{p}_s =\tilde{B}(n,s)\dot{n} - \frac{(p_s-f(n))^2}{(1-\kappa(s)n)^3} \kappa'(s)n\,,
\end{equation*}
where
\begin{equation*}
  \tilde{B}(n,s) = \frac{1}{L}\int_{0}^{L}B(n,\xi) \dd \xi-B(n,s)\,.
\end{equation*}
We recall that, for all $t\in[0,T)$, $H(n(t),s(t),p_n(t),p_s(t)) = H_0$. We get
\begin{equation}\label{estimates 1}
\begin{split}
\vert \dot{n} \vert &\leq \sqrt{2H_0}\\
\abs{p_s-f(n)} &\leq \sqrt{2 H_0}(1+\epsilon)\\
\left\lvert \frac{(p_s-f(n))^2}{(1-\kappa(s)n)^3} \kappa'(s)n \right\rvert &\leq \frac{2H_0 K' \delta}{1-\epsilon}  \,,
\end{split}
\end{equation}
where in the last estimates we have used the notation of
Theorem~\ref{theo.confinement1} and in particular
$\abs{\kappa} n \leq K \delta \leq \epsilon$.
With our assumption~\eqref{boundness} on $\tilde{B}(n,s)$, we find,
for all $t\in[0,T)$,
\begin{eqnarray*}
\left\lvert p_s(t)\right\rvert  \leq   \vert p_s(0)\vert + \left( M\sqrt{2H_0} + \frac{2H_0 K' \delta}{1-\epsilon} \right) T\,,
\end{eqnarray*}
and thus
\begin{eqnarray}\label{eq.f <C}
\left\lvert  f(n(t)) \right\rvert &\leq& \left\lvert p_s(t)  \right\rvert + \left\lvert p_s(t) - f(n(t)) \right\rvert \leq C(T)\,,
\end{eqnarray}
with
\begin{equation*}
  C(T) = \vert p_s(0)\vert +\sqrt{2H_0}(1+\epsilon) +\left( M\sqrt{2H_0} + \frac{2H_0 K' \delta}{1-\epsilon} \right)  T\,.
\end{equation*}
If the trajectory reaches the boundary at $t=T$ , then
\begin{equation*}
\lim_{t\to T} n(t) = 0\,.
\end{equation*}
This, with \eqref{eq.f <C} and \eqref{eq.lim.of.integral}, gives a contradiction. This proves Theorem \ref{theo.confinement0}.

Now, consider a function $g$ as in Theorem \ref{theo.confinement1}. We have, for all $t\in[0,T)$,
\begin{equation*}
g(n(t)) \leq |f(n(t))| \leq  C(T)\,.
\end{equation*}
From~\eqref{eq:f very large}, we have $\lim_{n\to 0} g(n)>C(T)$; hence
$g$ must take the value $C(T)$ and the conclusion follows.

\subsection{Proof of Proposition \ref{prop.of.confinement2}}
Let us recall \eqref{eq.red.H}. The assumptions of Proposition \ref{prop.of.confinement2} can be written in terms of $V$.
\begin{enumerate}[\rm (H1)]
\item If
\begin{equation}\label{eq.liminfV1>H}
\liminf_{r\to 1^-} V(r) > H_0\,,
\end{equation}
we consider $\eta=\sup\{x\in(0,1) : V(x)=H_{0}\}\in(0,1)$. Consider a
trajectory $(q(t), p(t))$ with $q(0)\in D(0,1)$. We can assume that
$q(0)\neq 0$.  Let $T$ be the maximal time of existence in
$D(0,1)$. By energy conservation, we have, for all $t\in[0,T)$,
\[V(r(t))\leq H_{0}\,,\]
so that $r(t)\leq\eta$.

Note that \eqref{eq.liminfV1>H} means
\[\liminf_{r\to 1^-} \left\lvert G(r) - p_\theta \right\rvert>   \sqrt{2 H_0}\,.\]
Using the usual complex coordinate in the plane $\RR^2$, we can write
$\dot{q}=\left(\dot{r}+i\dot{\theta}r\right)e^{i\theta}$ and thus
\begin{equation*}
  \det (q(t), \dot{q}(t)) = r^2(t) \dot{\theta}(t) = p_\theta - G(r(t))\,.
\end{equation*}
Finally, we notice that $\|\dot{q}(0)\|=\sqrt{2H_{0}}$ and write
\[G(r) - p_\theta = G(r)-G(r(0)) - [p_\theta - G(r(0))]\,,\]
which gives~\eqref{eq.criterion.1}.

\item If
\begin{equation}\label{eq.liminfV1=H}
\liminf_{r\to 1^-} V(r) = H_0\,,
\end{equation}
and
\begin{eqnarray*}
\limsup_{r\to 1^-}\frac{V(r)- H_0}{r-1} <0\,,
\end{eqnarray*}
then we must again have
\[
\sup\{x\in(0,1) : V(x)=H_{0}\} < 1\,,
\]
and we can proceed as above.
\end{enumerate}

\subsection{Proof of Proposition \ref{prop.confinement.3}}
Consider $p_{\theta}=0$. Let
$\abs{V}_\infty := \sup_{r\in(0,1)}\abs{V(r)}$. By assumption,
$\abs{V}_\infty< +\infty$.

Let $r(0) \in (0,1)$ and choose $p_r(0)>0$ such that $ p_r^2(0) = 2\left(|V|_{\infty}-V(r(0))\right)+v^2$, with $v>0$. Since, for all $t\in[0,T)$,
\begin{equation*}
\frac{p_r^2(t)}{2}+V(r(t))=\frac{p_r^2(0)}{2}+V(r(0))\,,
\end{equation*}
we get $\dot{r}(t)= p_r(t) > v$ so that
\begin{equation*}
r(t) > vt+r(0)\,.
\end{equation*}
The particle escapes at $t=\frac{1-r(0)}{v}$.

\subsection{Proof of Theorem \ref{theo.scattering}}

We distinguish between the cases $p_{\theta}=0$ and $p_{\theta}\neq 0$.

\subsubsection{Case when $p_{\theta}\neq 0$}
In this case, $\lim_{r\to 0} V(r)=+\infty$; hence, due to energy
conservation, the trajectory does not approach the origin.
\begin{enumerate}[\rm i.]
\item Assume that $p_{r}(0)<0$. We have $V(1)<H_{0}$ and we can
  consider the right most turning point $r^*\in(0,1)$. By definition
  $V(r^*)=H_{0}$, and necessarily $V'(r^*)\leq 0$.

If $V'(r^*)<0$, it is easy to check that $r$ reaches $r^*$ in finite time, say $t=t^*$. This time is given by
\[t^*=\int_{r^*}^{1}\frac{\dd r}{\sqrt{2(H_{0}-V(r))}}\,.\]
By symmetry, the escape time is $2t^*$. Since $\dot{\theta}=\frac{p_{\theta}-G(r)}{r^2}$, we have
\[\theta(t^*)-\theta(0)=\int_{0}^{t^*} \frac{p_\theta-G(r)}{r^2} \dd t=\int_{0}^{t^*} \frac{(p_\theta-G(r)) \dot{r}}{r^2 p_r} \dd t=\int_{0}^{t^*}  - \frac{(p_\theta-G(r)) \dot{r}}{r^2\sqrt{2(H_0-V(r))}  }\dd t\,,\]
so that
\[\theta(t^*)-\theta(0)= \int_{r^*}^{1}  \frac{p_\theta-G(r)}{r^2\sqrt{2(H_0-V(r))} } \dd r\,.\]
By symmetry, we have
\[\theta(2t^*)-\theta(0)=2\int_{r^*}^{1}  \frac{p_\theta-G(r)}{r^2\sqrt{2(H_0-V(r))} } \dd r\,.\]

If $V'(r^*)=0$, $(r^*,0)$ is a critical point of the Hamiltonian and we get that $r$ reaches $r^*$ in infinite time (see Figure \ref{confine4}).

\item Assume that $p_{r}(0)=0$. Then $V(1)=H_{0}$. By assumption (the
  trajectory enters $D(0,1)$), we have $V'(1)\geq 0$, \textit{i.e.},
  $(p_{\theta}-G(1))B(1)+(p_{\theta}-G(1))^2\leq 0$ . If $V'(1)=0$,
  the particle sits at a fixed point of the Hamiltonian system, and
  hence $r(t)\equiv 1$ is constant. If $V'(1)>0$, it enters $D(0,1)$
  and the discussion is the same as previously.
\end{enumerate}
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{limitcircle.eps}
\caption{$B(r)= e^{-r}-\frac{2}{r}$ .} \label{confine4}
\end{figure}

\subsubsection{Case when $p_{\theta}= 0$}
In this case, since $G(0)=0$, $V(r)=\frac{1}{2r^2}G(r)^2$ admits a
continuous extension at $r=0$.

\begin{enumerate}[\rm i.]
\item Assume that $p_{r}(0)<0$. We have $V(1)<H_{0}$. The existence of
  $r^*$ such that $V(r^*)=H_{0}$ is not ensured. If $V(r)<H_{0}$ on
  $[0,1]$, the particle reaches $r=0$ in finite time $t=t^*$:
\[t^*=\int_{0}^{1}\frac{\dd r}{\sqrt{2(H_{0}-V(r))}}\,.\]
We get, by symmetry,
\[\theta(2t^*)-\theta(0)=2\int_{0}^{1}  \frac{-G(r)}{r^2\sqrt{2(H_0-V(r))} } \dd r+\pi\,.\]
If there exists $r^*\in(0,1)$ such that $V(r^*)=H_{0}$, the trajectory does not reach the origin and the discussion is the same as in the case $p_{\theta}\neq 0$.

\item Assume that $p_{r}(0)=0$. The discussion is the same as when $p_{\theta}\neq 0$.

\end{enumerate}

\subsubsection{Scattering angle}
We can now end the proof of Theorem \ref{theo.scattering}. In terms of complex numbers, we can write
\[v_{1}=(v_{r}(0)+iv_{\theta}(0))e^{i\theta_{1}}\,,\quad v_{2}=(-v_{r}(0)+iv_{\theta}(0))e^{i\theta_{2}}\,.\]
The scattering angle is
\[\theta_{2}-\theta_{1}+\mathrm{Arg}\left(\frac{-v_{r}(0)+iv_{\theta}(0)}{v_{r}(0)+iv_{\theta}(0)}\right)\,.\]
If $\gamma$ denotes the argument of $v_{r}(0)+iv_{\theta}(0)$, the scattering angle is thus
\[\theta_{2}-\theta_{1}+\pi-2\gamma\,.\]

\appendix
\section{Tubular coordinates}
\begin{lemma} \label{lem.cartesian.to.normal}
We write $\mathbf{A}=A_{1}\dd q_{1}+ A_{2}\dd q_{2}$. With \eqref{eq.psi}, we have
\[\mathbf{A}=A_{n}\dd n+A_{s}\dd s\,,\qquad \tilde A=(A_{n}, A_{s})^\mathrm{T}=(d\psi)^\mathrm{T}(A_{1}, A_{2})^\mathrm{T}\,.\]
We have
\begin{equation} \label{Ham in normal coor}
H(n,s,p_n,p_s) =\mathcal{H}\circ \Psi(n,s,p_{n}, p_{s})= \frac{(p_n - A_n(n,s))^2}{2} + \frac{(p_s- A_s(n,s))^2}{2(1-\kappa (s) n)^2}\,.
\end{equation}
Moreover, $v_{n}=p_n - A_n(n,s)$ and $v_{s}=(1-n\kappa(s))^{-1}(p_{s}-A_{s})$ are the normal and tangential component of $v$.
\end{lemma}
\begin{proof}
We write
\[2H(q,p)=\|p-A\|^2=\|(\dd\psi^{-1})^\mathrm{T}(\tilde p-\tilde A)\|^2=\langle (\dd\psi^{-1})(\dd\psi^{-1})^\mathrm{T}(\tilde p-\tilde A), \tilde p-\tilde A\rangle\,,\]
with $\tilde p=(p_{n}, p_{s})^{\mathrm{T}}$. Note that
\begin{equation}\label{eq.dpsi-1T}
(\dd\psi^{-1})^\mathrm{T}=[N(s)\,\,, (1-n\kappa(s))\gamma'(s)]\,.
\end{equation}
We get
\[(\dd\psi^{-1})(\dd\psi^{-1})^\mathrm{T}=\begin{pmatrix}1&0\\
0&(1-n\kappa(s))^{-2}\end{pmatrix}\,.\]
Concerning the velocity $v$, we write
\[v=p-A=(\dd\psi^{-1})^\mathrm{T}(\tilde p-\tilde A)\,,\]
and we use \eqref{eq.dpsi-1T}.
\end{proof}


TITLE: Is it correct? Prove that any fraction can be reduced
QUESTION [2 upvotes]: I want to know if my prove is correct.
My goal is proving:

Hypothesis: $a,b \in \mathbb Z$ and $a,b \notin \{-1, 0, 1\}$.
Thesis: for all $ a, b$, exist $a',b'\in \mathbb Z$ that verify $$\frac ab = \frac{a'}{b'}$$ and $a'/b'$ is irreductible, that is, $\gcd (a',b') = 1$.

I do it by reductio ad absurdum. I suppose that $\gcd(a', b') \neq 1$. That is, they've got some common factor.
Then I factorize $a$ and $b$ as a primer number product. So:
$$a = \varepsilon \cdot p_1^{a_1} p_2^{a_2} \cdots p_r^{a_r}$$
$$b = \gamma \cdot p_1^{b_1} p_2^{b_2} \cdots p_r^{b_r}$$
where $\varepsilon, \gamma = \pm 1$, $r \geq 1$, $a_i \geq 0$, $b_i \geq 0$ and every primer number $p_i$ is different.
So, I define $d$ as $\gcd(a,b)$, so, by definition:
$$d = p_1^{\min(a_1, b_1)} \cdots p_r^{\min(a_r, b_r)}$$
So I deduce:
$$\frac ab = \dfrac{\dfrac ad}{\dfrac bd} = \dfrac{\varepsilon \cdot\dfrac{ p_1^{a_1} p_2^{a_2} \cdots p_r^{a_r}}{p_1^{\min(a_1, b_1)} \cdots p_r^{\min(a_r, b_r)}}}{\gamma \cdot\dfrac{ p_1^{b_1} p_2^{b_2} \cdots p_r^{b_r}}{p_1^{\min(a_1, b_1)} \cdots p_r^{\min(a_r, b_r)}}} =  \frac\varepsilon\gamma \cdot
\dfrac{p_1^{a_1 - \min(a_1, b_1)} \cdots p_r^{a_r - \min(a_r, b_r)}}{p_1^{b_1 - \min(a_1, b_1)} \cdots p_r^{b_r - \min(a_r, b_r)}} $$
Then I set $$a' = \frac\varepsilon\gamma \cdot p_1^{a_1 - \min(a_1, b_1)} \cdots p_r^{a_r - \min(a_r, b_r)}$$
 $$b' = p_1^{b_1 - \min(a_1, b_1)} \cdots p_r^{b_r - \min(a_r, b_r)}$$
so $a/b = a'/b'$.
Now, as I have supposed that $a'$ and $b'$ have a factor in common, then $\exists i, a_i = b_i$, $a_i \neq 0$ and $b_i \neq 0$.
So, $\forall i$, we see the behaviour of $a_i$ and $b_i$ by cases.

If $a_i = b_i = 0$, then they haven't got factor $p_i$ in common.
If $a_i > b_i$, then $\min(a_i,b_i) = b_i$.
$a_i - \min(a_i,b_i) = a_i - b_i \neq 0$
$b_i - \min(a_i, b_i) = b_i - b_i = 0$
So $a_i \neq b_i$.
If $b_i > a_i$, by symmetry, also $a_i \neq b_i$.

We haven't found any $i$ that verify $\exists i, a_i = b_i$, $a_i \neq 0$ and $b_i \neq 0$, but we have supposed it previously! Contradiction.
So thesis is true.

REPLY [2 votes]: Simpler: Let $c=\gcd(|a|,|b|)$. Then $a/b=(a/c)/(b/c)$ and $\gcd(|a/c|,|b/c|)=1$.


TITLE: Contravariant components of metric tensor
QUESTION [0 upvotes]: I have a question regarding the contravariant components of the metric tensor in spherical coordinates.
I have calculated the covariant components as $g_{rr}=1$, $g_{\theta\theta}=r^2\sin\phi$, and $g_{\phi\phi}=r^2$ (all others are 0).
My question is how would I get the contravariant components from here?
I am told that $g^{ij}$ is such that $g^{ij}g_{jk}=\delta^i_k$ but have no clue how to proceed.
Thanks, in advance.
Eduardo

REPLY [2 votes]: First, I think you should have $g_{\theta\theta}=r^2\sin^2\phi$.
All this tensor notation is telling you to do is to find the inverse of the $3\times 3$ matrix $\big[g_{ij}\big]$. Since the matrix is diagonal, it's easily done!


TITLE: Transformation of random variable involving projectile distance
QUESTION [0 upvotes]: Problem:
Suppose a projectile is ﬁred at an angle $θ$ from the horizontal with a velocity $v$. The horizontal distance that the projectile travels, $R$, is given by
$$R = \frac{v^2}{g} \sin 2 \theta$$
where $g$ is the acceleration due to gravity ($9.8$ m/s2)
If $ \theta$ is uniformly distributed over the range $(0,π/4)$, ﬁnd the probability density function of $R$.
. . . . . . . . . . . . . . .
I know the answer involves the change of variable method. I get as far as solving the Jacobian derived to be $\dfrac{2v^2\cos2\theta}{g}$. At this point I am unsure how to proceed further to achieve the answer  which is given as $\dfrac{2g}{\pi\sqrt{v4-g^2r^2}}$ for $0 < r <v^2/g$
Any help would be appreciated.

REPLY [1 votes]: Hint:
Go the other way, i.e. write $\theta$ as a function of $r$: $\theta = v(r)$ the change of variable formula gives $f_R(r) = f_\theta (v(r)) |v'(r)|$ where $f_R$ is the distribution function for $R$ and $f_\theta = 4 / \pi$ and is uniform.


TITLE: guided proof that there are infinitely many primes on the arithmetic progression $4n + 3$
QUESTION [3 upvotes]: The exercise is to prove this proposition using the following lemmas:
1) Prove Lemma A: every prime ≠ 2 is either of the form $4n + 1$ or $4n + 3$
2) Prove Lemma B: the product of numbers of the form $4n + 1$ is of the form $4n + 1$.
I'll omit the proofs of the lemmas, since I'm not doubtful about their demonstrations. Again I ask if there is any flaw on the following proof.
Proof:
Suppose, by contradiction, that there are finitely many primes of the form $4n + 3$, say, $p_1, p_2, ..., p_k.$ Define $N = 4p_1...p_k - 1 = 4(p_1...p_k - 1) + 3.$ $N$ > $p_k$ and is of the form $4n + 3$, therefore is not prime. So $N$ must be divisible by some prime. Since $N$ is not divisible by any $p_i$ or 2, this prime divisor must be > $p_k$ and consequently, due to lemma A, of the form $4n + 1$. Now, the fundamental theorem of arithmetic states that any non-prime number must be factorable in prime numbers. Since $N$ is not divisible by any $p_i$, there must be primes $f_1, ..., f_l$ on which $N$ is factorable. From the hypothesis and from lemma A follows that these primes must be of the form $4n + 1$ and, by lemma B, so must be $N$, a contradiction.

REPLY [3 votes]: Your proof is good, however consider the following(slightly shorter) version:
Suppose there are finitely many primes of the form $4n + 3$, call them $p_1, p_2, \dots p_k$.  Then consider $$N = 4(p_1p_2\dots p_k) - 1$$ This is still of the form $4n + 3$ since it is congruent to $-1 \equiv 3(\mod 4)$.  Note also that none of the $p_i$ divide $N$. Thus, $N$ is either prime or composite.  If $N$ is prime, we reach a contradiction since it is a new prime of the desired form.  Suppose $N$ is composite, then there must be a prime of the form $4n+3$ dividing it by your lemma, also a contradiction since none of the $p_i$ divide $N$.
This proof is nearly identical to yours, however the thing to note is that you can just leave it in the form $4k - 1$ since this is equivalent to $4(k-1) + 3$.


TITLE: which is the cardinal of primes?
QUESTION [0 upvotes]: I know that $|\mathbb{N}|=|\mathbb{Z}|=|\mathbb{Q}|=\aleph_0$ but which is the cardinal of $|\mathbb{P}|$? since there is no formula for primes, then I can conclude that there is not bijection between $\mathbb{N}$ and $\mathbb{P}$ and so $|\mathbb{P}|\ne\aleph_0$?

REPLY [4 votes]: In short $|\mathbb{P}| = \aleph_0$. There are infinitely many primes, and every prime is an integer. So the primes are a countably infinite set, just like $\mathbb{N}$, $\mathbb{Z}$ and $\mathbb{Q}$. To put them in their place, consider the sum:
$$\sum_{k=1}^{\infty} \frac{1}{k^2} = \frac{\pi^2}{6} \, . $$
However, the sum of reciprocals $1/p$ where $p$ is prime, i.e. $1/2 + 1/3 + 1/5 + 1/7 + 1/11 + \cdots $ actually diverges! So, in some sense, there are "more" primes than there are square numbers. In short: there are "less" primes than there are natural numbers, but "more" primes than there are square numbers.


TITLE: How do you derive the eigenvalue and eigenvector of a linear combination of outer products?
QUESTION [2 upvotes]: I'm reading a paper that claims that it is easy to verify that the eigenvector-eigenvalue pair of
\begin{align*}
p_1\beta_1\beta_1^\top+p_2\beta_2\beta_2^\top,
\end{align*}
where $p_1$ and $p_2$ are just two non-zero and non-equal constants that sum to one, is
\begin{align*}
(v_1,\lambda_1)&=\left(\frac{\sqrt{\frac{1-\Delta_2}{1+\Delta_2}}\beta_1-\beta_2}{\sqrt{\frac{1+\Delta_2}{2}}+\sqrt{\frac{(1-\Delta_2)(1-\Delta_1)}{2(1+\Delta_2)}}},\frac{1+\sqrt{1-4(1-\langle\beta_1,\beta_2\rangle^2)p_1p_2}}{2}\right)\\
(v_2,\lambda_2)&=\left(\frac{\beta_2+\sqrt{\frac{1+\Delta_2}{1-\Delta_1}}\beta_1}{\sqrt{\frac{1-\Delta_2}{2}}+\frac{1+\Delta_2}{2}\sqrt{\frac{2}{1-\Delta_1}}},\frac{1-\sqrt{1-4(1-\langle\beta_1,\beta_2\rangle^2)p_1p_2}}{2}\right),
\end{align*}
where
\begin{align*}
\Delta_1&=\frac{(\lambda_1-\lambda_2)^2+p_1^2-p_2^2}{2(\lambda_2-\lambda_1)p_1} \\
\Delta_2&=\frac{(\lambda_2-\lambda_1)^2+p_2^2-p_1^2}{2(\lambda_1-\lambda_2)p_2}.
\end{align*}
I've been trying to verify the above statement but am unable to do so. Can anyone guide me as to how I should go about doing this? Thanks.
My attempt: I tried to show that
\begin{align*}
\lambda_1v_1&=(p_1\beta_1\beta_1^\top+p_2\beta_2\beta_2^\top)v_1 \\
\lambda_2v_2&=(p_1\beta_1\beta_1^\top+p_2\beta_2\beta_2^\top)v_2,
\end{align*}
by substituting the equations on the RHS and hoped to simplify to the LHS. However, the equations just got even more complex and I see no way of simplifying them.

REPLY [3 votes]: You can assume without loss of generality that $\beta_1$ and $\beta_2$ are unit norm (by integrating it into coefficients of $p_1$ and $p_2$).
Looking for a first eigenvalue-eigenvector relationship of this form:
$$(p_1\beta_1\beta_1^\top+p_2\beta_2\beta_2^\top)(a_1\beta_1+a_2\beta_2)=\lambda_1(a_1\beta_1+a_2\beta_2)$$
and denoting by $s$ the dot product of $\beta_1$ and $\beta_2$, we get, by expanding and decomposing onto basis $\beta_1,\beta_2$:
$$\begin{cases}p_1a_1+p_1sa_2&=&\lambda_1a_1\\
p_2a_1s+p_2a_2&=&\lambda_1a_2\end{cases}$$
which means that $\lambda_1$ can be found as (one of the two) eigenvalue(s) of the following matrix:
$$\begin{pmatrix}p_1&p_1s\\
p_2s&p_2\end{pmatrix}$$
and that $(a_1;a_2)$ will be an associated eigenvector to this eigenvalue.
(the second eigenvalue and eigenvector is therefore obtained in this way).
It remains do the last calculations whose resultsare very simple:
$$\begin{cases}\lambda_{1}&=&\dfrac12(p_1+p_2 - \sqrt{\Delta})\\
\lambda_{2}&=&\dfrac12(p_1+p_2 + \sqrt{\Delta})\end{cases} \ \text{with} \ \Delta:=(p_1-p_2)^2 + 4p_1p_2s^2$$
With resp. eigenvectors:
$$V_1=\binom{p_1 - p_2 -\sqrt{\Delta}}{2sp_2}$$
$$V_2=\binom{p_1 - p_2 +\sqrt{\Delta}}{2sp_2}$$
Why all these complicated formulas in the article ? A plausible reason is that they have taken raw results issued from a Computer Algebra System.


text:
\section{Conclusions}\label{sec:conclusions}

In this paper, we presented the design of sparse samplers for inverse problems with tensors. We have seen that by using samplers with a Kronecker structure we can overcome the curse of dimensionality, and design efficient subsampling schemes that guarantee a good performance for the reconstruction of multidomain tensor signals. We presented sparse sampling design methods for cases in which the multidomain signals can be decomposed using a multilinear model with a dense core or a diagonal core. For both cases, we have provided a near-optimal greedy algorithm based on submodular optimization methods to compute the sampling sets.

config: arxiv


TITLE: Torque on accelerating tire
QUESTION [0 upvotes]: So I have a problem understanding the torques on a tire as it accelerates from rest. Suppose I have a tire with mass $m = 10kg$, radius $r = 0.5m$ and I apply torque with magnitude $\tau = 30Nm$ to its axis. The force with which the tire pushes on the ground should be equal to $\tau / r$, so 60N. The maximum friction force $F_{f(max)} = \mu  m  g = 78.4N$, where $\mu = 0.8$ and $g = 9.8m/s^2$
so the force from the tire doesn't exceed the maximum friction force. And that means it won't slip. According to Newton's third law, the force produced by the tire is counteracted by an equal in magnitude and opposite in direction friction force $F_f = -F = 60N$. That friction creates torque about the axis of rotation which is in opposite direction to the torque I started with and has magnitude $\tau_f = F_f  r = 30Nm$ so they should completely cancel out, leaving 0 net torque on the tire and that means it won't rotate. I must be wrong somewhere because the tire will obviously accelerate when a torque is applied to it. I've read somewhere that some of the torque goes for the angular acceleration of the tire and not all of it produces the force at the bottom of the tire. Is that in any way correct?
Any help would be appreciated!

REPLY [0 votes]: I apply torque with magnitude =30Nm to its axis.

Sure.

The force with which the tire pushes on the ground should be equal to /, so 60N.

If the wheel were massless, that would be true.  100% of the torque from the engine would be "delivered" to the ground.  But the wheel is not massless.  Some of the torque is instead going into accelerating the wheel.
Using the sum of torques formula:
$$\tau_{net} = I\alpha$$
If the wheels rotation is accelerating, then there must be a non-zero net torque.  Therefore the torque from the axle and the torque from the ground cannot be equal.
There is a force pair where the force from the tire on the ground is exactly equal to the force from the ground on the tire.  But when examined as a torque, these will be less than the torque from the engine.
One way to think about this is the linear example of pushing two blocks forward.  You push on A, and A pushes on B.  If A is massless, B is pushed with an equal magnitude force.
But the more mass A has, the smaller the force A pushes on B.  In your car scenario, the larger the moment of inertia of the wheel, the smaller the torque that the ground creates.

And that would mean that if the wheel were massless, then the torque would be zero and it would not roll, right?

No.  A massless wheel can roll or accelerate.  It just cannot have any net torque applied.  Any torque placed on it must be countered by an equal torque elsewhere.  It's an ideal limit.

But is there a way to precisely calculate the force at the bottom from engine torque apart from the formula for net torque?

You have to solve simultaneous equations for the system.  You know the moment of inertia for the wheel, and you know the mass of the car.
The torque by the ground on the wheel will be an amount between zero and the engine torque, and that allows the wheel acceleration to match the vehicle acceleration.
If the torque is too high, the vehicle moves faster than the wheel turns.  If the torque is too low, the wheel spins faster than the vehicle moves.
$$\alpha = \frac{a}{t}$$
$$F_{road} = ma$$
$$\tau_{engine} - (F_{road} \times r) = I\alpha$$


TITLE: Change of variable. Integral vs continuous expectation
QUESTION [0 upvotes]: Suppose, the goal is to perform a change of variable in computing an expectation over continuous variables.
Now, I am wondering if the rules that can be applied are only the rules for change of variables in any integral or there are a set of rules also specific to expectations that one can take advantage of here?

REPLY [1 votes]: Evaluating an expectation by integration is ... an integration.   You may use the usual change of variables rules if you need to.


TITLE: Find function for multiplication operator such that there are at least 2 eigenvalues
QUESTION [1 upvotes]: Consider $V=C([a,b])$ and let $A:V\to V: f(t)\mapsto u(t)f(t)$ be the multiplication operator with fixed $u(t)\in V$.

Let $E$ denote the set of eigenvalues of operator $A$.

Let $S$ denote the spectrum of operator $A$. The spectrum is defined as the set of non-regular points; i.e. the points $\mu$ for which the resolvent $(A-\mu\mathbf{1})^{-1}$ is not well-defined or bounded on $V$.

It is clear that $S= \operatorname{im}(u)$, the image of $u$. By choosing $u(t)$ wisely, we can create 3 different situations: 1) $E=\emptyset$, 2) $E = S$ and 3) $\emptyset\ne E \subset S$. My question is regarding this third situation.
A possible example (let's work on $[0,1]$), would be to consider $u(t) = 0$ outside of $[1/3,2/3]$ and have it form a little baseless triangle within $[1/3,2/3]$ (to ensure $u$ is continuous); and then we can choose $f$ to be $0$ inside $[1/3,2/3]$, but non-zero outside of it (but still continuous). Then $u(t)f(t)=0 = 0f(t)$, so this way we have found $E=\{0\}\subset S=[0,1]$ (albeit without explicitly writing $f$).
Now, I want to find an example (just by drawing, not necessarily expressions) of $u(t)$ such that $|E|>1$ but $E\subset S$ strictly. I'm having trouble finding continuous functions for this to work. Any ideas?
Thanks.

REPLY [1 votes]: You can show that, in general,  $E$ will be the set of $x\in\mathbb{R}$ such that $\operatorname{int}(u^{-1}(x))\neq\varnothing$ (use an argument similar to the one you did with the little triangle). So just take a function whose graph looks like a stair to get more eigenvalues.
It will also follow that $E$ is always countable.


TITLE: Semilinear Parabolic PDE
QUESTION [5 upvotes]: I'm considering the following type of PDE:
$u_{t}=u_{xx}+u_{x}+u_{x}^2+u_{x}^3+\frac{u_{x}}{x(1-x)}+\left(\frac{u_{x}}{x(1-x)}\right)^3$
with periodic boundary conditions $u_{x}(0)=u_{x}(1)=0$.
Does anyone know of literature/articles I can look at to possibly show local existence(or blow up?). I've spent quite a bit of time looking through literature in particular Henry's Geometric theory of Semi-Linear Parabolic equations and papers by H.Amman on Quasi linear Parabolic equations.  Also due to the singularity I considered a weighted sobolev space approach but had no success.  I apologize if I'm too vague.

REPLY [2 votes]: I would try to find the Lie Symmetry groups for this equation first. A good book on this is Bluman, Cheviakov, Anco. The easiest to find them is to use some computer algebra packages such as Maple and e.g. the related GeM package from Cheviakov.


TITLE: Feynman's Integration technique, parameter finding
QUESTION [2 upvotes]: Recently I was studing the Feynman integration technique (differentiation under the integral sign), but I allways get stuck on the general function definition, and I was wondering if there is some kind of formula, trick, table, or something else that can help me evaluate such a funtion.
Thankfully liuzp.
Update:
For functions like:
$$\int_a^b \frac{f(x)}{g(x)} \hspace{1cm} f(x,\alpha )=e^{\alpha\ g(x) + c(x)} + d(x)$$
where $c(x)$ and $d(x)$ are special functions that make $f(x, \alpha ) = f(x)$ for some special value of $\alpha$.
But I've also founf that this won't work for all cases, and for cases where you haven't a denominator.
This also have a limitation because $g(x) \in \{\mathbb{P}^2,\ln[h(x)]\} $

REPLY [2 votes]: Here's a typical example.
$$ J = \int_0^\infty \frac{\arctan(x)}{x(1+x^2)}$$
We'd like to get rid of that nasty arctan.  We know the derivative of $\arctan$ is nice.
So we generalize to
$$ J(t) = \int_0^\infty \frac{\arctan(tx)}{x(1+x^2)}$$
take the derivative with respect to $t$, and the rest is easy...


text:
\begin{document}

\maketitle

\begin{abstract}
We derive both the finite and infinite population spatial replicator dynamics as the fluid limit of a stochastic cellular automaton. The infinite population spatial replicator is identical to the model used by Vickers and our derivation justifies the addition of a diffusion to the replicator. The finite population form generalizes the results by Durett and Levin on finite spatial replicator games. We study the differences in the two equations as they pertain to a one-dimensional rock-paper-scissors game.
\end{abstract}

\section{Introduction}
Evolutionary games using the replicator dynamic have been studied extensively and are now well documented \cite{ref:Fr91,ref:Wei95,ref:HS98,ref:HS03,T15,T19}. Variations on the classical replicator dynamic include discrete time dynamics \cite{ref:VRS11} and mutations \cite{ref:BTG12,ref:TS15}. Additional evolutionary dynamics, such as imitation  \cite{ref:Fr91,ref:Sc97,ref:FL98,ref:HS03} and exchange models \cite{ref:EGB16} have been studied. Alternatively evolutionary games have been extended to include spatial models by a number of authors \cite{ref:andrew,ref:kerr,ref:nowak,ref:roca,ref:RPSnetworks,ref:vickers1989,ref:vickers1991,ref:cressmanVickers}. Most of these papers append a spatial component to the classical replicator dynamics (see e.g., \cite{ref:vickers1989}) or discuss finite population replicator dynamics in which total population counts are used (see e.g., \cite{ref:andrew}). In the latter case, a spatial term is again appended to the classical replicator structure.

In \cite{DL94}, Durrett and Levin study discrete and spatial evolutionary game models and compare them to their continuous, aspatial analogs. For their study the authors focus on a specific class of two-player two-strategy games using a hawk-dove payoff matrix. Because their payoff matrix is $2 \times 2$, a single (spatial) variable $p(x,t)$ can be used to denote the proportion of the population playing hawk (Strategy 1) while a second variable $s(x,t)$ denotes the total population. Using the payoff matrix:
\begin{displaymath}
  \mathbf{A} = \begin{bmatrix}a & b\\c&d\end{bmatrix},
\end{displaymath}
the remarkable reaction-diffusion equation is analyzed:

\begin{align}
  \frac{\partial p}{\partial t} &= \Delta p + \frac{2}{s}\nabla p\cdot \nabla s + pq((a-c)p + (b-d)(1-p))
  \label{eqn:Durret1}
  \\
  \frac{\partial s}{\partial t} &= \Delta s + s\left(ap^2 + (b+c)p(1-p)+(1-p)^2d\right) - \kappa s^2.
  \label{eqn:Durret2}
\end{align}
Here $\kappa$ is a death rate due to overcrowding and $q = (1-p)$. Let $\mathbf{e}_i \in \mathbb{R}^n$ be the $i^\text{th}$ unit vector. Let $\mathbf{u} = \langle{p,q}\rangle$. When $\kappa = 0$, we can rewrite these equation as:
\begin{align}
  \frac{\partial p}{\partial t} &= \Delta p + \frac{2}{s}\nabla p\cdot \nabla s + p\left(\mathbf{e}_1^T-\mathbf{u}^T\right)\mathbf{A}\mathbf{u}\label{eqn:DL1}\\
  \frac{\partial s}{\partial t} &= \Delta s + s\cdot\mathbf{u}^T\mathbf{A}\mathbf{u}.
\end{align}
That is, Durrett and Levin have encoded the replicator dynamic into a finite population spatial partial differential equation, which differs from the one used by Vickers \cite{ref:vickers1989} because it assumes a finite population in its derivation. The second term in \cref{eqn:DL1} is an advection of species $p$ along a velocity vector given by the gradient of the whole population. However, this advection speed is inversely proportional to the population size. Thus small populations can have a dramatic effect on this term. As we show in \cref{sec:FluidLim}, this behavior  carries through for general games. To illustrate the advection clearly, we can define the potential function:
\begin{displaymath}
\phi(x) = 2\log\left[s(x,t)\right].
\end{displaymath}
We then have:
\begin{displaymath}
\mathbf{v} = \nabla\phi = \frac{2}{s}\nabla s.
\end{displaymath}
\cref{eqn:DL1} can then be written as:
\begin{displaymath}
\frac{\partial p}{\partial t} = \Delta p + \mathbf{v}\cdot \nabla p + p\left(\mathbf{e}_1^T-\mathbf{u}^T\right)\mathbf{A}\mathbf{u}\label{eqn:DL1a},
\end{displaymath}
giving a standard advection-diffusion equation. The logarithm in the potential function conveniently yields the per-capita advection rate.

In this paper we show that Durrett and Levin's finite population
spatial replicator is generalizable to an arbitrary payoff matrix. We
then focus our attention on the one-dimensional rock-paper-scissors
game under the replicator dynamics in both finite and infinite populations. This system has interesting properties in both the finite and infinite
population cases. In particular, we show: (i) the  model used by
Vicker's \cite{ref:vickers1989} arises naturally as the infinite
population limit of the generalization of Durrett and Levin's model,
which in turn can be derived from a stochastic cellular automaton
(particle) model as a fluid limit, and (ii) the one dimensional infinite population spatial rock-paper-scissors dynamic has a constant amplitude traveling wave solutions in biased RPS. However, the finite population version does not exhibit such solutions, but does seem to exhibit attracting stationary solutions. We illustrate the latter result numerically.

\section{Model}\label{sec:Model}
Let $\mathbf{A} \in \mathbb{R}^{n \times n}$ be a payoff matrix for a symmetric game \cite{ref:HS03}. All vectors are column vectors unless otherwise noted. Below we construct a stochastic cellular automaton and show that the fluid limit of this system yields a generalization of Durrett and Levin's specific finite population spatial replicator.

The state of cell $i$ at time index $k$ of the cellular automaton is a tuple $\langle{U_1(i,k),\dots,U_n(i,k)}\rangle$ where $U_j(i,k)$ provides the size of the population of species $j$ at position $i$ at time $k$. For simplicity, we assume that species interaction may only happen between cells and not within cells; i.e., the $U_1(i,k)$ members of species $1$ will not play  against the $U_n(i,k)$ members of species $n$. This assumption will become irrelevant in the limit.

During state update, an agent $A$ at cell $i$ chooses a random direction (cell $i'$ ) and a random member of the population (agent $A'$) within that cell. Assume Agent $A$ uses strategy $r$ while Agent $A'$ uses strategy $s$. After play, there are $\alpha\cdot\mathbf{A}_{rs}$ additional agents at cell $i$ using strategy $r$ and $\beta\cdot\mathbf{A}_{rs}$ additional agents playing strategy $r$ at cell $i'$, where $\alpha + \beta = 1$ is the probability of motion from cell $i$ to cell $i'$.
If $\mathbf{A}_{rs} < 0$, then agents are removed from their respective cells. To avoid computational issues when more members of a species die than are present, $\mathbf{A}$ can be modified so that $\mathbf{A}_{rs} \geq 0$ for all $r,s \in \{1,\dots,n\}$, without altering the evolutionary dynamics in proportion \cite{ref:Wei95}. The update rule for a single agent is illustrated in \cref{fig:Update}.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.4\columnwidth]{Figures/Figure.png}
\caption{A single interaction between two agent in a spatial game on
  a one-dimensional lattice is illustrated. The direction of interaction is chosen at random. Motion is random and governed by $\alpha$ and $\beta$ where $\alpha+\beta=1$.}
\label{fig:Update}
\end{figure}
The process described above and illustrated in \cref{fig:Update} is assumed to be happening simultaneously for each agent and we assume that the replication/death as a result of game play along with the migration are happening on (roughly) the same time scale.

Using these assumptions, we can construct mean-field equations for the population $U_r$ at position $i$ and time $k+1$ in a 1-D cellular automaton, assuming an equal likelihood that agents diffuse left or right. The mean number of agents playing strategy $r$ present at position $i$ at time $k+1$ is determined by:
\begin{enumerate}
\item The expected number of agents who remain at cell $i$: $\alpha U_r(i,k)$
\item The expected number of new agents created at cell $i$ who remain at cell $i$:
\begin{displaymath}
\frac{\alpha}{2}U_r(i,k)\left(\sum_{s}\mathbf{A}_{rs}\left(u_s(i+1,k) + u_s(i-1,k)\right)\right).
\end{displaymath}
\item The expected number of agents who migrate to position $i$ from neighboring cells:
\begin{displaymath}
\frac{\beta}{2}\left(U_r(i+1,k)+U_r(i-1,k)\right).
\end{displaymath}
\item The expected number of agents created in a neighboring cell who migrate to cell $i$:
\begin{displaymath}
\frac{\beta}{2}\left(U_r(i+1,k)+U_r(i-1,k)\right)\sum_{s}\mathbf{A}_{rs}u_s(i,k).
\end{displaymath}
\end{enumerate}
Here $u_s(i,k)$ is the \textit{proportion of the population} at cell $i$ playing strategy $s$ at time step $k$.

Let $\mathbf{u}(i,k) = \langle{u_1(i,k),\dots,u_n(i,k)}\rangle$. Re-writing sums as matrix products, the expected number of agents playing strategy $r$ at cell $i$ at time $k+1$ is:
\begin{multline}
U_r(i,k+1) =
\alpha U_r(i,k) + \frac{\beta}{2}\left(U_r(i+1,k)+U_r(i-1,k)\right) +\\
\frac{\alpha}{2}U_r(i,k)\left(\mathbf{e}_r^T\mathbf{A}\left(\mathbf{u}(i+1,k) + \mathbf{u}(i-1,k)\right)\right) + \\
\frac{\beta}{2}\left(U_r(i+1,k)+U_r(i-1,k)\right)\mathbf{e}_r^T\mathbf{A}\mathbf{u}(i,k).
\label{eqn:MatrEqn}
\end{multline}

Assume the cellular grid has lattice spacing $\Delta x$. Following \cite{ref:DGB+14} and using a Taylor approximation, we can write:
\begin{align*}
    U_r(x,t+\Delta t) &\approx U_r(x,t) +
        \Delta t \frac{\partial U_r(x,t)}{\partial t} + O(\Delta t^2)\\
    U_r(x+\Delta x,t) &\approx \sum_{j=0}^{2}\frac{\Delta x^j}{j!}
        \frac{\partial^j U_r(x,t)}{\partial x^j} + O(\Delta x^3).
\end{align*}

\section{Derivation of Fluid Limits}\label{sec:FluidLim}
We proceed to derive the mean-field approximation. Passing to the continuous case and assuming that interaction rates decrease linearly with $\Delta t$, we can write a second order approximation of \cref{eqn:MatrEqn} as:
\begin{multline}
\Delta t\frac{\partial U_r(x,t)}{\partial t} = \alpha U_r(x,t) +
\beta\left(U_r(x,t) + \frac{1}{2}\Delta x^2\frac{\partial^2 U_r}{\partial x^2}\right) + \\
\alpha\Delta t U_r(x,t)\cdot \mathbf{e}_r^T\mathbf{A}\left(\mathbf{u}(x,t) + \frac{1}{2}\Delta x^2\frac{\partial^2 \mathbf{u}}{\partial x^2} \right) + \\
\beta\Delta t\left(U_r(x,t) + \frac{1}{2}\Delta x^2\frac{\partial^2 U_r}{\partial x^2}\right)\mathbf{e}_r^T\mathbf{A}\mathbf{u}(x,t) - U_r(x,t).
\end{multline}
Where the $-U_r(x,t)$ on the right-hand-side arises from the formation of the Newton quotient on the left-hand-side.

Expanding and simplifying  yields:
\begin{multline}
\Delta t\frac{\partial U_r(x,t)}{\partial t} =
\frac{\beta}{2}\Delta x^2\frac{\partial^2 U_r}{\partial x^2} +
\Delta t U_r(x,t)\mathbf{e}_r^T\mathbf{A}\mathbf{u}(x,t) + \\
\frac{\alpha\Delta t}{2} \Delta x^2 U_r(x,t)\mathbf{e}_r^T\mathbf{A}\frac{\partial^2 \mathbf{u}}{\partial x^2} + \\
\frac{\beta\Delta t}{2}\Delta x^2\frac{\partial^2 U_r}{\partial x^2}\mathbf{e}_r^T\mathbf{A}\mathbf{u}(x,t).
\end{multline}
Assume $\beta \in (0,1]$. Dividing through by $\Delta t$ and assuming that $\lim_{\Delta t \rightarrow 0}\Delta x^2/\Delta t = 2D/\beta$ yields:
\begin{equation}
\frac{\partial U_r(x,t)}{\partial t} =  U_r(x,t)\mathbf{e}_r^T\mathbf{A}\mathbf{u}(x,t) +
D\frac{\partial^2 U_r}{\partial x^2}.
\label{eqn:righteqn}
\end{equation}
The constant $D$ is the diffusion constant and the assumption that
$$\lim_{\Delta t \rightarrow 0}\Delta x^2/\Delta t = 2D/\beta$$
is a variant of the assumption used to derive Fick's Law \cite{O09} and identical when $\beta = 1$.

These are the spatial dynamics used by Durrett and Levin, (in the first part of their paper), but are derived only by adding a diffusion term to the standard finite population growth equations. In \cite{DL94}, Durrett and Levin note that they derive a set of equations they feel are more appropriate for modeling finite spatial systems. Their derivation at the end of \cite{DL94} (for a specific hawk-dove system) rests on the assumption that migration happens ``on a much faster timescale'' than game interactions. Our model assumes that migration and game interactions occur on approximately the same time scale. Under this assumption, \cref{eqn:righteqn} is the correct spatial adaptation for finite populations; i.e., one simply adds a diffusion term. In the case where migration happens more quickly, then the derivation in \cite{DL94} should be used instead.

\section{Spatial Replicator with Finite Population}
The derivation of \cref{eqn:Durret1,eqn:Durret2} are not given in \cite{DL94}. They can be generalized for an arbitrary evolutionary game using \cref{eqn:righteqn} as the starting point. Let:
\begin{displaymath}
M(x,t) = \sum_{s}U_s(x,t).
\end{displaymath}
Differentiating we have:
\begin{displaymath}
\frac{\partial}{\partial t}\frac{U_r(x,t)}{M(x,t)} = \frac{1}{M(x,t)}\frac{\partial U_r(x,t)}{\partial t} - \\
u_r(x,t)\sum_{s}\frac{1}{M(x,t)}\frac{\partial U_s(x,t)}{\partial t}.
\end{displaymath}
Substituting from \cref{eqn:righteqn} we obtain:
\begin{equation}
\frac{\partial u_r}{\partial t} = u_r\cdot\left(\mathbf{e}_r^T\mathbf{A}\mathbf{u} -\mathbf{u}^T\mathbf{A}\mathbf{u}\right)+\\
\frac{D}{M}\left(\frac{\partial^2 U_r}{\partial x^2} - u_r\cdot\sum_{s}\frac{\partial^2 U_s}{\partial x^2} \right).
\label{eqn:PartialResult}
\end{equation}
Unlike in the derivation of the standard replicator dynamic, the rate of change of the population proportion is not solely a function of the proportions themselves.

We can remove dependence on the individual populations to derive an independent (coupled) system of differential equations that includes only the total population.
For arbitrary strategy $r$, we can apply the quotient rule to obtain:
\begin{displaymath}
M\frac{\partial u_r}{\partial x} = \frac{\partial U_r}{\partial x} - u_r\frac{\partial M}{\partial x}
\end{displaymath}
Differentiating again, multiplying by $1/M$ and re-arranging yields the expression:

\begin{displaymath}
\frac{1}{M}\frac{\partial^2 U_r}{\partial x^2} = \frac{u_r}{M}\frac{\partial^2 M}{\partial x^2} + \frac{2}{M}\frac{\partial u_r}{\partial x}\frac{\partial M}{\partial x} + \frac{\partial^2 u_r}{\partial x^2}.
\end{displaymath}
Using this we can write

\begin{displaymath}
\frac{D}{M}\left(\frac{\partial^2 U_r}{\partial x^2} - u_r\sum_{s}\frac{\partial^2 U_s}{\partial x^2} \right) =
D\left(\frac{2}{M}\frac{\partial M}{\partial x}\frac{\partial u_r}{\partial x}  + \frac{\partial^2 u_r}{\partial x^2} \right)
\end{displaymath}
Using this we can re-write \cref{eqn:PartialResult} as:
\begin{equation}
\frac{\partial u_r}{\partial t} = u_r\cdot\left(\mathbf{e}_r^T\mathbf{A}\mathbf{u} -\mathbf{u}^T\mathbf{A}\mathbf{u}\right)+\\
 D\left(\frac{2}{M}\frac{\partial M}{\partial x}\frac{\partial u_r}{\partial x}  + \frac{\partial^2 u_r}{\partial x^2} \right).
\end{equation}
The dynamics of $M$ can be derived (by addition) from \cref{eqn:righteqn}:
\begin{displaymath}
\frac{\partial M}{\partial t} = M\mathbf{u}^T\mathbf{A}\mathbf{u} + D\frac{\partial^2 M}{\partial x^2}.
\end{displaymath}
Thus, we have a coupled set of differential equations written entirely in terms of $\mathbf{u}$ and $M$, rather than $U_r$, $M$ and $\mathbf{u}$:
\begin{equation}
\left\{
\begin{aligned}
&\forall r \left\{\begin{aligned}
\frac{\partial u_r}{\partial t} = & u_r\cdot\left(\mathbf{e}_r^T\mathbf{A}\mathbf{u} -\mathbf{u}^T\mathbf{A}\mathbf{u}\right)+\\
& \hspace*{2em}D\left(\frac{2}{M}\frac{\partial M}{\partial x}\frac{\partial u_r}{\partial x}+ \frac{\partial^2 u_r}{\partial x^2} \right)
\end{aligned}\right.\\
&\hspace*{2em}\frac{\partial M}{\partial t} = M\mathbf{u}^T\mathbf{A}\mathbf{u} + D\frac{\partial^2 M}{\partial x^2}
\end{aligned}
\right..
\label{eqn:finitespatialrep}
\end{equation}
This is the spatial replicator equation for \textit{finite}
populations. Letting $M = s$ and $D = 1$, we recover the dynamics of
Durrett and Levin. In contrast to the aspatial replicator, the
inclusion of dynamics for $M$ yields a linearly independent system of
differential equations.

Allowing $M$ to approach infinity uniformly in $x$, we arrive at the
fluid limit in terms of $\mathbf{u}$ alone; this is the 1D nonlinear reaction-diffusion equation used by Vicker's \cite{ref:vickers1989,ref:vickers1991}:

\begin{equation}
\forall r \left\{
\begin{aligned}
\frac{\partial u_r}{\partial t} &= u_r\cdot\left(\mathbf{e}_r^T\mathbf{A}\mathbf{u} -\mathbf{u}^T\mathbf{A}\mathbf{u}\right) + D\frac{\partial^2 u_r}{\partial x^2}
\end{aligned}
\right..
\label{eqn:spatialrep}
\end{equation}
Generalization to $N$-dimensions is straightforward by replacing $\partial^2_x$ with the Laplacian $\Delta$. The $N$-dimensional spatial replicator with finite population is given by:
\begin{align*}
&\forall r \left\{\begin{aligned}
\frac{\partial u_r}{\partial t} = & u_r\cdot\left(\mathbf{e}_r^T\mathbf{A}\mathbf{u} -\mathbf{u}^T\mathbf{A}\mathbf{u}\right)+
\left(\frac{2}{M}\nabla M \cdot \nabla u_r +  \Delta u_r\right)
\end{aligned}\right.\\
&\hspace*{2em}\frac{\partial M}{\partial t} = M\mathbf{u}^T\mathbf{A}\mathbf{u} + D\Delta M.
\end{align*}

Thus, the finite population case adds a nonlinear advection term that forces $u_r$ to follow the population gradient. A similar system is studied by deForest and Belmonte in \cite{ref:andrew}, where the payoff gradient is followed instead of the population gradient.

In both \cref{eqn:finitespatialrep,eqn:spatialrep}, we see that the
aspatial replicator dynamics appear on the right hand side perturbed
by a spatial term. It is well known that the dynamics of the aspatial
replicator are confined to the $n$-dimensional simplex $\Delta_n$.
This remains true for the spatial replicator dynamics with finite populations.

Moreover, the solution $u_r = 1$ (i.e., there is only one population)
is a fixed point for the spatial replicator dynamic since the spatial
derivative of the probability distribution of the population
proportions is zero and the time derivative is identically zero as
expected. Thus, pure populations are constant stationary solutions for
these dynamics. Lastly, if $\tilde{\mathbf{u}}(x,t)$ is a constant
solution at a Nash equilibrium for the game defined by $\mathbf{A}$,
then the right-hand-side is again identically zero by the Folk Theorem
\cite{ref:HS03} of evolutionary game theory together with the fact
that there is no spatial variation. Thus every Nash equilibrium of the
matrix game corresponds to a spatially constant stationary solution of the spatial replicator dynamic in both the finite and infinite population cases.

\section{Example: One Dimensional Rock-Paper-Scissors}\label{sec:RPS}
Durrett and Levin's analysis of Hawk-Dove was aided by the fact that
one strategy can be eliminated, leaving a coupled system of two
partial differential equations. In the remainder of this paper, we analyze variations of rock-paper-scissors (RPS), which yield more interesting results because of its cyclic three-strategy nature and because it can be easily parameterized as discussed in \cite{ref:Wei95}. Substantial work has been done on (spatial) RPS without assuming a replicator dynamic or assuming a general replicator dynamic \cite{ML75,M10,PR19,SMR14,SMJS14,RMF07,RMF08,HMT10,PR17,SMR13}. In an example of a very recent generalization,  Kabir \& Tanimoto study pairwise evolution in RPS with noise \cite{KT21}.

A major focus of work in non-replicator spatial RPS has been on identifying spatial structures (e.g., spirals) that can emerge from these dynamics. For completeness, we quantify the relationship between the work in \cite{PR19,SMR14,SMJS14,RMF07,RMF08,HMT10,PR17,SMR13} and the replicator dynamic considered in this work in Appendix A. There, we show that the equation system used in these references does not subsume results on the replicator dynamic in general and hence one cannot extend these results automatically to this case.

The objective in the remainder of this work is to compare the finite and infinite population replicator for a more interesting three strategy game. For simplicity (and in contrast to much of the previous work in spatial RPS), we focus our attention on the one-dimensional PDE.

We consider a generalized RPS payoff matrix as given by:
\begin{displaymath}
\mathbf{A} = \begin{bmatrix}
0 & -1 & 1+a\\
1 +a & 0 & -1\\
-1 & 1+a & 0
\end{bmatrix}.
\end{displaymath}
When $a = 0$, this is the standard RPS game which has Nash equilibrium
$\langle{\tfrac{1}{3},\tfrac{1}{3},\tfrac{1}{3}}\rangle$. This is the
unique interior fixed point and the aspatial replicator exhibits an elliptic fixed point at this Nash equilibrium. This Nash equilibrium is preserved for $a \neq 0$;  and
corresponds to an asymptotically stable interior fixed point when $a
> 0$ and unstable fixed point when $a < 0$ \cite{Ze80} as a result of a degenerate Hopf bifurcation. This degenerate Hopf bifurcation is well understood and leads to traveling wave solutions, which we discuss in \cref{sec:TravelingWave}. We note $\mathbf{A}$ is not as fully general as the RPS matrix given in \cite{ref:HS03} or \cite{Ze80}, but it is substantially easier to work with.

Let $\mathbf{u} = \langle{u_r,u_p,u_s}\rangle$ and note:
\begin{displaymath}
\zeta(u_r,u_p,u_s) \overset{\Delta}{=} \mathbf{u}^T\mathbf{A}\mathbf{u} = a\left(u_ru_s + u_ru_p + u_s u_p\right).
\end{displaymath}

There are at least two classes of global solutions to \cref{eqn:finitespatialrep,eqn:spatialrep}:
\begin{description}
\item[Equilibrium Solution] Here, $u_r = u_p = u_s = \tfrac{1}{3}$ and $M$ solves:
\begin{displaymath}
0 = \frac{a}{3}M + M''
\end{displaymath}
or there is a single population (e.g., $u_r = 1$) and $M$ satisfies $M''=0$.

\item[Oscillating Solution] Here, $u_*(x,t) \equiv \upsilon_*(t)$ with $* \in \{r,p,s\}$ where $\upsilon_r$, $\upsilon_p$, $\upsilon_s$ are solutions to the standard RPS replicator and $M(x,t) = \mu(t)$ satisfies linear equation:
\begin{displaymath}
\dot{\mu} = \zeta(\upsilon_r,\upsilon_p,\upsilon_s)\mu.
\end{displaymath}
\end{description}
In either case, we may impose appropriate boundary conditions.

\subsection{Traveling Wave Solutions}\label{sec:TravelingWave}
Consider the infinite population model. Letting $z = x - ct$, we can re-write \cref{eqn:spatialrep} in compact form as:
\begin{equation}
\forall i\left\{
\begin{aligned}
Dv_i' &= -u_i\left(\mathbf{e}_i^T\mathbf{A}\mathbf{u} - \mathbf{u}^T\mathbf{A}\mathbf{u}\right) -c v_i\\
u_i' &= v_i
\end{aligned}\right..
\label{eqn:TravelingWave}
\end{equation}
For RPS we have the six dimensional non-linear system:
\begin{equation}
\left\{
\begin{aligned}
Dv_r' &= a u_r\zeta(u_r,u_s,u_p) - u_r(u_s - u_p) - a u_r u_s - c v_r\\
u_r' &= v_r\\
Dv_p' &= a u_p\zeta(u_r,u_s,u_p) - u_p(u_r - u_s) - a u_p u_r - c v_p\\
u_p' &= v_p\\
Dv_s' &= a u_s\zeta(u_r,u_s,u_p) - u_s(u_p - u_r) - a u_s u_p - c v_s\\
u_s' &= v_s.
\end{aligned}
\right.
\label{eqn:SysTravelWave}
\end{equation}
If $\mathbf{u}^*$ is a Nash equilibrium of $\mathbf{A}$, then the pair $\mathbf{u} = \mathbf{u}^*$ and $\mathbf{v} = \mathbf{0}$ is a fixed point of \cref{eqn:TravelingWave}. Linearizing about $u_r = u_p = u_s = \tfrac{1}{3}$, $v_r = v_p = v_s = 0$, we obtain the eigenvalues:
\begin{align*}
\lambda_{1,2} &= \frac{-3c \pm\sqrt{9c^2+12aD}}{6D}\\
\lambda_{3,4} &= \frac{-3c \pm \sqrt{9c^2 + 6aD + 6D\sqrt{3}\sqrt{-(a+2)^2}}}{6D}\\
\lambda_{5,6} &= \frac{-3c \pm \sqrt{9c^2 + 6aD - 6D\sqrt{3}\sqrt{-(a+2)^2}}}{6D}.\end{align*}
We can simultaneously show that for appropriate choice of wave speed, a Hopf bifurcation and hence a two dimensional center manifold exists and therefore a non-decaying traveling wave solution exists for the PDE. As a by-product, we compute the wave speed for a non-decaying traveling wave in terms of $a$ and $D$. For some constant $b$ (to be determined), let:
\begin{multline*}
(3c \pm bi)^2 = 9c^2 + 6aD \pm 6D\sqrt{3}\sqrt{-(a+2)^2} = \\
9c^2 + i6aD \pm 6D(a+2)\sqrt{3}.
\end{multline*}
Expanding the left hand side and relating real and imaginary parts we have:
\begin{align*}
&9c^2 - b^2 = 9c^2 + 6aD\\
&6bc  = 6D(a+2)\sqrt{3}.
\end{align*}
Solving for $b$ and $c$ yields:
\begin{align*}
b & = \sqrt{-6aD}\\
\pm\tilde{c} &= \pm\frac{(a+2)\sqrt{2k}}{2\sqrt{-a}}
\end{align*}
Without loss of generality, assume positive $\tilde{c}$. Using this information, we can rewrite the eigenvalues as:
\begin{align*}
\tilde{\lambda}_{1,2} &= \frac{-3\tilde{c} \pm\sqrt{9\tilde{c}^2+12aD}}{6D}\\
\tilde{\lambda}_{3}&= \frac{-6\tilde{c} - bi}{6D}\\
\tilde{\lambda}_{4}&= \frac{bi}{6D}\\
\tilde{\lambda}_{5}&= \frac{-6\tilde{c}+bi}{6D}\\
\tilde{\lambda}_{6}&= \frac{-bi}{6D}.
\end{align*}
For this solution to be physically realized, we must have $a < 0$. We also assume $a > -2$ or the dynamics changes (i.e., winning becomes losing). These requirements and our assumption that $\tilde{c} > 0$ implies that $\mathrm{Re}(\lambda_{1,2}) < 0$ for all choices of $D > 0$ and $a \in (-2,0)$. Therefore, this system has a four dimensional stable manifold and two pure imaginary eigenvalues, which satisfies the first requirement of Hopf's theorem (see \cite{GH13}, Page 152).

Consider $\lambda_{4,6}$ as functions of $c$ with:
\begin{displaymath}
\lambda_{4,6}(c) = \frac{-3c + \sqrt{9c^2 + 6aD \pm 6D\sqrt{3}\sqrt{-(a+2)^2}}}{6D}.
\end{displaymath}
For $c = \tilde{c} > 0$, we know that $\mathrm{Re}[\lambda_{4,6}(\tilde{c})] = 0$. Differentiating we have:
\begin{displaymath}
\lambda_{4,6}'(c) = - \frac{1}{2D} + \frac{3 c}{2D\sqrt{9c^2 + 6aD \pm 6D\sqrt{3}\sqrt{-(a+2)^2}}}.
\end{displaymath}
Evaluating at $c = \tilde{c}$ and simpilfying yields:
\begin{displaymath}
\lambda_{4,6}'(\tilde{c}) = - \frac{1}{2D} + \frac{3 \tilde{c}\left(3\tilde{c}\mp bi\right)}{2D\left(9\tilde{c}^2 - b^2\right)},
\end{displaymath}
by choice of $\tilde{c}$. We conclude:
\begin{displaymath}
\mathrm{Re}\left[\lambda_{4,6}'(\tilde{c})\right] = - \frac{1}{2D} + \frac{9\tilde{c}^2}{2D\left(9\tilde{c}^2 - b^2\right)} = \frac{b^2}{2D\left(9\tilde{c}^2 - b^2\right)}.
\end{displaymath}
This is non-zero since we assume $a < 0$ and $D > 0$. Thus the real parts of eigenvalues $\lambda_{4,6}$ cross the imaginary axis with non-zero speed, satisfying the second criterion of Hopf's theorem. Thus we conclude that the six dimensional traveling wave ODE has a solution and moreover exhibits a Hopf bifurcation, implying the existence of traveling wave solutions. In our numeric simulations, we show that fine tuning the parameters leads to a numerically stable traveling wave over the region of integration when $a = -\tfrac{4}{5}$ and $D = \tfrac{1}{12}$ (see \cref{fig:NegBiasPlots,fig:NegativeBiasTernary}).

\subsection{Illustrative Behavior in the Unbiased Games}
Consider \cref{eqn:finitespatialrep,eqn:spatialrep} and assume the periodic boundary conditions $u_{*}(-\pi,t) = u_{*}(\pi,t)$ for all time. When $a = 0$, then $\zeta(u_r,u_p,u_s) = 0$, so the population $M(x,t)$ satisfies the heat equation. For simplicity, we choose a solution to the heat equation that models the diffusion of a population outward:
\begin{displaymath}
M(x,t) = \frac{4000 e^{-\frac{x^2}{4 \beta  (t+10)}}}{\sqrt{\pi } \sqrt{\beta  (t+10)}}.
\end{displaymath}
The initial conditions:
\begin{align*}
u_r(x,0) &= \frac{1}{3}\left[1+\sin\left(x - \tfrac{4\pi}{3}\right)\right]\\
u_p(x,0) &= \frac{1}{3}\left[1+\sin\left(x - \tfrac{2\pi}{3}\right)\right]\\
u_s(x,0) &= \frac{1}{3}\left[1+\sin\left(x\right)\right]
\end{align*}
model three populations that are proportionally spread around a circle. The behavior of the three population proportions in both finite and infinite populations are shown in \cref{fig:NoBiasPlots}.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\columnwidth]{Figures/FiniteInfiniteRowFigNoBias.png}
\caption{Illustration of finite and infinite dynamics at various points in time on the circle for RPS with zero bias. The diffusing population causes finite population to converge to a stationary oscillating solution, while the infinite population converges to a stationary equilibrium solution.}
\label{fig:NoBiasPlots}
\end{figure}
In the infinite population case we observe that the three population proportions converge toward the equilibrium stationary solution. This is further illustrated in \cref{fig:NoBiasTernary} (left) where we show a ternary plot of the three strategies at $x = 0$. By contrast the finite population solution converges toward the oscillating stationary solution, as illustrated in \cref{fig:NoBiasTernary} (right). We show the corresponding solution to the ordinary RPS replicator dynamic to which the system converges at all points in space. The convergence of the infinite population system to the equilibrium at all points in space is illustrated in \cref{fig:NoBiasTernary} (left).
\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\columnwidth]{Figures/InfinitePopulationAsymptoticBehaviorNoBias.png} \quad
\includegraphics[width=0.45\columnwidth]{Figures/FinitePopulationAsymptoticBehaviorNoBias.png}
\caption{(left) The infinite population spatial replicator converges to an equilibrium point as illustrated from multiple starting points on the circle. (right) The finite population spatial replicator converges to an oscillating stationary solution as illustrated from multiple starting points on the circle. The oscillating solution is a solution to the replicator dynamic with unbiased RPS.}
\label{fig:NoBiasTernary}
\end{figure}
Before converging to an oscillating solution, we can see the one-dimensional proportions proportions are affected by the diffusion of the population at large. In \cref{fig:NoBiasPlots}, we note that the finite population plots are stretched with respect to their infinite population counterparts. This is particularly noticeable at time $t = 5$ and $t = 10$. This stretching, caused by the advection of the total population, leads to the difference in steady state solutions for the same initial conditions.

We contrast this behavior with the case when $a < 0$. Here, the population will collapse since $\zeta(u_r,u_p,u_s) \leq 0$ at all times. As noted, we set $a = -\tfrac{4}{5}$ and $D = \tfrac{1}{12}$. From this we compute a constant amplitude traveling wave speed of:
\begin{displaymath}
\tilde{c} = \frac{1}{2}\sqrt{\frac{3}{10}}.
\end{displaymath}
This traveling wave can be seen in \cref{fig:NegBiasPlots} in the infinite population plots.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\columnwidth]{Figures/FiniteInfiniteRowFigNegativeBias.png}
\caption{Illustration of finite and infinite dynamics at various points in time on the circle for RPS with negative bias. In the infinite population case, a stable traveling wave emerges. In the finite population case, population collapse causes the population proportions to swing with ever increasing amplitude.}
\label{fig:NegBiasPlots}
\end{figure}
In contrast the population collapse (shown in \cref{fig:PopulationDecay}) causes the finite population proportions to converge to an oscillating stationary strategy with greater and greater amplitude. This is precisely the behavior one expects to see from RPS under the replicator dynamics when $a < 0$. This is also shown \cref{fig:NegativeBiasTernary}(right), in which we show an example RPS trajectory with $a = -\tfrac{4}{5}$ and example trajectories at various points in space.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\columnwidth]{Figures/PopulationDecay.png}
\caption{The population collapses exponentially until $t \approx 30$. Additionally interactions with the individual species cause the population to become asymmetric as illustrated by the trajectories of $M(x^*,t)$ at $x^* = \pm \tfrac{\pi}{2}$.}
\label{fig:PopulationDecay}
\end{figure}
Interestingly, the population collapse slows dramatically as $t$ increases. This is caused by the fact that for (nearly) pure strategies, $\zeta(u_r,u_s,u_t) \approx 0$. This is illustrated in \cref{fig:PopulationDecay}, where the nearly exponential decay flattens  after $t \approx 30$. Additionally, \cref{fig:PopulationDecay} shows that the population becomes asymmetric as the system evolves, as a result of the various species dynamics.

Convergence to the stable traveling wave solution is illustrated in \cref{fig:NegativeBiasTernary} (left). Using the computed $\tilde{c}$, numerical analysis provided initial conditions which can be used to find a numerical solution to  \cref{eqn:SysTravelWave} that produces the closed curve (on the center manifold) that is guaranteed to exist by the eigenvalue analysis performed in the previous section. This is shown in \cref{fig:NegativeBiasTernary} (left).
\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\columnwidth]{Figures/InfinitePopulationAsymptoticBehaviorNegativeBias.png}\quad
\includegraphics[width=0.45\columnwidth]{Figures/FinitePopulationAsymptoticBehaviorNegativeBias.png}
\caption{(left) The infinite population spatial replicator converges to a limit cycle of the 6 dimensional traveling wave ODE as illustrated from multiple starting points on the circle. (right) The finite population spatial replicator converges to an oscillating stationary solution as illustrated from multiple starting points on the circle. The oscillating solution is a solution to the replicator dynamic with negative biased RPS and thus is converging to the boundary of the simplex.}
\label{fig:NegativeBiasTernary}
\end{figure}

\section{Conclusion}
In this paper we studied a finite and infinite population spatial
replicator. We showed how the finite population spatial replicator can
be derived from first principles from a stochastic cellular automaton
model and from there how the infinite population replicator used by
Vickers \cite{ref:vickers1989,ref:vickers1991} follows from this. This result generalizes the work of Durrett and Levin \cite{DL94} who first derived and studied the finite population spatial replicator for a specific game. We then compared the finite and infinite population spatial replicator for rock-paper-scissors on a circle ($S^1$). Our results are consistent with the work in \cite{RMF07,RMF08,HMT10,SMR13,SMR14,SMJS14,PR17,PR19}, which studies various characteristics of RPS, not necessarily in connection with spatial replicator. Consistent with the work in \cite{PR17,PR19} we show that for a certain rock-paper-scissors variant stable amplitude traveling waves can emerge as solutions to the infinite population spatial replicator by proving the existence of a Hopf bifurcation in the traveling wave ODE. These traveling waves are destroyed by population collapse in the finite population spatial replicator.

The finite population spatial replicator is intriguing because it is a highly non-linear reaction-advection-diffusion equation where advection is governed by the per capita bulk population motion. Identifying cases where complex behaviors emerge in the finite population case is a clear future direction. Additionally, studying stationary solutions may yield insights. For example, in our unbiased RPS, solutions to the stationary population equation are just the harmonic functions. Using this simplification may help identify interesting properties of the population proportion equations.

\section*{Acknowledgement}
Portions of CG's work were supported by the National Science Foundation Grant CMMI-1932991.

\appendix
\section{Relation to Other Rock Paper Scissors Games}\label{sec:A}
A substantial amount of work has been done on rock-paper-scissors outside of the context of the replicator dynamic. This work does not map conveniently to results on the replicator dynamic or its spatial variants
\cite{PR19,M10,SMR14,SMJS14,RMF07,ML75,RMF08,HMT10,PR17,SMR13,KT21}. The
earliest work to study competition among three species with cyclic
dominance may be \cite{ML75}, which is contemporary but does predate
the earliest work in evolutionary game theory (see e.g.,
\cite{TJ78,Ze80,Hof81,B83}). In the past twenty years, there has been
substantial work on the spatial dynamics of RPS that is independent of
the spatial replicator dynamic
\cite{ref:andrew,ref:kerr,ref:nowak,ref:roca,ref:RPSnetworks,ref:vickers1989,ref:vickers1991,ref:cressmanVickers}.
Work till 2014 is reviewed in \cite{SMJS14}. Mobility in cyclic
competition (RPS) is studied in \cite{RMF07,RMF08}. The impact of
reaction rates on spatial RPS is considered in \cite{HMT10}, while the
emergence of spiraling waves is studied in \cite{SMR13,SMR14}. More
recently, traveling waves, spirals and heteroclininc bifurcations have
been studied in \cite{PR17,PR19}. A discrete time model displaying
chaos is considered in \cite{SAF02}.

The results most closely related to those in this paper (specifically \cref{sec:TravelingWave}) can be found in the pair of papers by Postlethwaite and Rucklidge \cite{PR17,PR19}. Both these papers study a traveling wave solution for a specific spatial model of the rock-paper-scissors, with a more formal treatment given in \cite{PR19}. The motivating aspatial model is given by the system of equations:
\begin{align}
\dot{a} &= a (1-(a+b+c) -(\sigma + \zeta )b + \zeta c) \label{eqn:Leeds1}\\
\dot{b} &= b (1-(a+b+c) -(\sigma + \zeta )c + \zeta a) \label{eqn:Leeds2}\\
\dot{c} &= c (1-(a+b+c) -(\sigma + \zeta )a + \zeta b) \label{eqn:Leeds3}
\end{align}
The observation is made that when $\sigma = 0$, this system exhibits the conserved quantities $a+b+c = 1$ and $abc = K$ for some constant $K$. Letting $u_1 = a$, $u_2 = b$ and $u_3 = c$, then \crefrange{eqn:Leeds1}{eqn:Leeds3} can be written as the standard replicator:
\begin{equation}
\dot{u}_i = u_i\left(\left(\mathbf{e}_i - \mathbf{u}\right)^T\mathbf{A}\mathbf{u} \right),
\label{eqn:LeedsReplicator}
\end{equation}
where:
\begin{equation}
\mathbf{A} = \begin{bmatrix}
-1 & -\zeta -1 & \zeta -1 \\
 \zeta -1 & -1 & -\zeta  -1 \\
 -\zeta -1 & \zeta -1 & -1
\end{bmatrix}
\end{equation}
This is trivially diffeomorphic to the replicator dynamic with payoff matrix:
\begin{equation}
\tilde{\mathbf{A}} =
\zeta\begin{bmatrix}
0 & -1 & 1 \\
 1 & 0 & -1  \\
 -1 & 1 & 0
\end{bmatrix}.
\end{equation}
This is a scaled version of the standard (unbiased) RPS payoff matrix. However, when $\sigma \neq 0$, then \cref{eqn:Leeds1,eqn:Leeds2,eqn:Leeds3} can be written as:
\begin{equation}
\dot{u}_i = u_i\left(1 + \mathbf{e}_i^T\mathbf{B}\mathbf{u}\right),
\end{equation}
with:
\begin{equation}
\mathbf{B} = \begin{bmatrix}
-1 & -\sigma -\zeta -1 & \zeta -1 \\
 \zeta -1 & -1 & -\sigma -\zeta  -1 \\
 -\sigma -\zeta -1 & \zeta -1 & -1
\end{bmatrix}.
\end{equation}
Population evolution equations of this form are considered in \cite{PG16}.
When $\sigma = 0$, then $\mathbf{u}^T\mathbf{B}\mathbf{u} = \mathbf{u}^T\mathbf{A}\mathbf{u} = -1$, from which we derive \cref{eqn:LeedsReplicator}. For $\sigma \neq 0$, this is not the case and consequently the dynamics of \cref{eqn:Leeds1,eqn:Leeds2,eqn:Leeds3} are not trapped on the unit 2-simplex as will be the case when we study (un)biased spatial replicator dynamics in finite and infinite populations. To be clear, the authors of \cite{PR17,PR19} make no claim to this effect, however it does create an important distinction between the work in these papers and the work in this paper.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.45\columnwidth]{Figures/LeedsFigure.png}
\caption{We illustrate that the dynamics studied in \cite{PR17,PR19} are not trapped on the unit simplex for $\sigma \neq 0$, while the same dynamics are illustrative of an unbiased RPS game when $\sigma = 0$.}
\label{fig:Leeds}
\end{figure}
Additionally, as illustrated in \cref{fig:Leeds}, the ODE dynamics that give rise to the PDE are not trapped on the unit simplex (as they are in the replicator) thus suggesting distinct behaviors may be observed.

\bibliographystyle{unsrt}
\bibliography{mrabbrev,spatialreplicator,RevisionBib}

\end{document}


TITLE: Best known bounds for a product over primes in an interval
QUESTION [1 upvotes]: As the title ask, what is the best lower and upper bounds for the product below :
$$ \prod \limits_{x < p \leq y} \frac{p+1}{p}$$
such that $p$ denote the prime numbers in which fullfil the conditions under the product ?

REPLY [0 votes]: Though Dusart bounds are good ones for $\,\prod\limits_{p \le x}{\left(1 -1/p\right)}$, I think a better idea to the prime product asked by Ahmad is to explore the Dedekind psi function, used in the theory of modular functions, defined as $\,\psi(n)/n = \prod\limits_{p|n}\left( 1 +\frac{1}{p} \right)$.