text:
\section{Introduction}
\label{sec:intro}

One approach to the pedigreed quest for determining the existence or nonexistence of finite-time singularity in the Euler and Navier-Stokes equations of fluid dynamics (see, e.g., \cite{doering} for an overview)  is to analyze them with very specialized  initial conditions.  Because  the vorticity at a singularity must diverge \cite{B_etal94},   effort has been spent on understanding the behavior of interacting localized de-singularized vortex tubes, leading to the study of vortex reconnection and its role in turbulence.    Various configurations have been proposed and investigated numerically for reconnection in classical turbulence (e.g., \cite{Kt94,K05,HD11}), quantum turbulence (e.g.,  \cite{B_etal08,Z_etal12,VPK17}),  and for the existence of singularity (e.g., \cite{SP85,PGGB,BMH16,KM18,K18}).

The present work investigates the reduced model of Moffatt and Kimura  \cite{MK19,MK19b} that describes the interaction of two circular vortex rings.  To the extent this model accurately represents  Euler's ideal fluid equations of motion, we have shown the existence of  finite-time singularity.  Explicitly, we have shown that within the proposed limits of applicability of this model there exist solutions that blow up in finite time.

Section \ref{sec:eom} describes the Moffatt and Kimura (MK) model.  This is followed by Sec.\ \ref{sec:HS} where  the Hamiltonian structure of the MK model is given, which is essential for our analysis.  Here we discover two constants of motion for the  MK model. One invariant  serves as  the  Hamiltonian for its noncanonical Hamiltonian formulation (flow on a Poisson manifold; see  \cite{pjm98}), while  the other turns out to be a Casimir invariant.  The Hamiltonian formulation allows us, in Sec.\ \ref{sec:quad},  to obtain geometrical intuition about the solution space by examining the intersection of the level sets of the two invariants, akin to the visualization afforded by the constancy of the energy and angular momentum magnitude for the Euler equations that describe the free rigid body.
Also in  this section we show how to reduce the MK system to quadrature and obtain for special initial conditions explicit solutions that have exact Leray scaling, which is representative of the finite-time singularity. Conditions for singularity within the range of applicability of the derivation given  in Refs.~\cite{MK19,MK19b} are presented.
In Sec.\ \ref{sec:conclu} we summarize  our results and mention some  future avenues.  Appendices are provided that exhibit additional features of our results.

\section{The Moffatt-Kimura system}
\label{sec:eom}

The MK system is a three-dimensional system of ordinary differential equations that describes the evolution of two initially circular vortices of radius $R$ and circulations $\pm\Gamma$, located symmetrically on planes $x=\pm z\tan \alpha$, with pitch angle $\alpha$,  in an $(x,y,z)$ Cartesian coordinate system.  The system is written in dimensionless form using the space scale $R$,  time scale  $R^2/\Gamma$, and effective Reynolds number $R_\Gamma= \Gamma/\nu=: 1/\varepsilon >> 1$.  Here $\nu$  is the usual
kinematic viscosity of the fluid. It is assumed that the vortices have Gaussian cores of radius measured by  $\delta$, a separation measured by $s$, and a curvature given by $\kappa$,  which evolve in terms of the dimensionless time according to
\bal
\dot{(\delta^2)}&=\varepsilon - c_2 \frac{\kappa \delta^2}{s}\,,
\label{dotde0}\\
\dot{s}&=-c_2\kappa \left[\ln\left(\frac{s}{\delta}\right) +\beta_1\right]
\label{dots0}\\
\dot{\kappa}&=c_1\frac{\kappa}{s^2}\,,
\label{dotka0}
\eal
subject to the  inequality constraints  for applicability
\begin{equation}
\delta< s< 1/{\kappa}\,,
\label{ineq}
\end{equation}
which are assumed in the derivation.  Roughly speaking, these inequalities assure that the vortex cores, with size measured by $\delta$,  are small enough so that cores do not overlap as the rings merge, i.e., as their separation given  by $s$ decreases, and that portions of the rings are sufficiently far away, as measured by $1/\kappa$,  to allow far field expansion.
The constants appearing in Eqs.~\eqref{dotde0}, \eqref{dots0}, and \eqref{dotka0}    are given by
\begin{equation}
c_1=  \frac{\cos\alpha \sin\alpha}{4\pi} \qquad\mathrm{and}\qquad
c_2=  \frac{\cos\alpha }{4\pi} \,,
\label{alpha}
\end{equation}
and the parameter $\beta_1$ depends on the vortex core,  with  the value  $\beta_1=0.4417$ for a gaussian core profile.

If the MK system is to be a  reduction that  inherits the noncanonical Hamiltonian structure of Euler's fluid equations (see e.g.\ \cite{pjm98}) then it will have a Hamiltonian form upon setting the  $\varepsilon=0$.  Thus we investigate
 \begin{eqnarray}
   \dot{\delta}&=&-\frac{c_2}{2}\,   \frac{\kappa \delta }{s}
  \label{dotde}\\
 \dot{s}&=&- c_2\,  \kappa \left[ \ln\left(\frac{s}{\delta}\right) +\beta_1\right]
   \label{dots}\\
 \dot{\kappa}&=&c_1\,  \frac{\kappa}{s^2}\,.
  \label{dotka}
\end{eqnarray}
For  later use we record here the equation of motion for $x:=s/\delta$,
\begin{equation}
\dot{x}= - \frac{c_2 \kappa}{\delta}(\Lambda -1/2)\,,
\label{dotx}
\end{equation}
which is easily verified.  Here, for convenience,  we have defined
 \begin{equation}
 \Lambda := \ln\left(\frac{s}{\delta}\right) +\beta_1 = \ln(x) +\beta_1 \,.
 \label{La}
 \end{equation}
  The smallest value of $\Lambda$ consistent with inequality \eqref{ineq} occurs when $\delta=s$; thus,
  \begin{equation}
 \Lambda    \geq  \beta_1 \qquad \mathrm{or}\qquad  \Lambda -1/2 \geq \beta_1-1/2 = -0.0583 \,,
 \label{critLa}
 \end{equation}
 where the gaussian value of $\beta_1=0.4417$ is used in the  second expression.  We note here that
 \begin{equation}
 \Lambda =1/2\qquad \mathrm{ at}\qquad  x^*= e^{1/2-\beta_1}\approx 1.0600\,.
 \end{equation}

\section{Hamiltonian Structure}
\label{sec:HS}

\subsection{Generalities}

Hamiltonian systems are usually written in terms of canonically conjugate sets of variables, a configuration space coordinate and its conjugate momentum.  The noncanonical Hamiltonian description is one where the form in terms of canonical variables is not necessary and replaced by algebraic properties of the Poisson bracket.   The terminology noncanonical Hamiltonian was introduced in the context of the ideal fluid and magnetohydrodynamics in  \cite{pjmG80}, but the ideas date back to the work of Sophus Lie.  (See, e.g.,  \cite{pjm98} for review.)

Given an $n$-dimensional  phase space with  coordinates $z=(z^1, z^2,\dots, z^n)$, a system of ordinary differential equations has noncanonical  Hamiltonian form  if there exists a Poisson bivector $J$,  an antisymmetric  second rank contravariant  tensor,
and a conserved phase space function $H(z)$ such that  the equations  can be written as follows:
\begin{equation}
\dot{z}^i = \{z^i,H\}= J^{ij} \frac{\partial H}{\partial z^j}  \qquad i,j=1,2,\dots n\,,
\label{genHamform}
\end{equation}
where repeated indices are summed and the Poisson bracket defined on functions of the coordinate $z$,
\begin{equation}
\{f,g\}= \frac{\partial f}{\partial z^i} J^{ij}  \frac{\partial g}{\partial z^i} \,
\end{equation}
is bilinear, antisymmetric, and most importantly satisfies the Jacobi identity,
\begin{equation}
\{f,\{g,h\}\} +\{f,\{g,h\}\} +\{f,\{g,h\}\} =0\,,
\label{jacobi0}
\end{equation}
for all functions $f,g,h$.  Unlike for the canonical description,  the tensor $J$,  the Poisson tensor, may depend on the coordinate $z$. It generates the Hamiltonian vector field as depicted on the righthand side of the second equality of \eqref{genHamform}. In  coordinates, \eqref{jacobi0} is equivalent to  the vanishing of  the following purely antisymmetric three-tensor:
\begin{equation}
S^{ijk}= J^{i\ell}\frac{\partial J^{jk}}{\partial z^\ell} + J^{j\ell}\frac{\partial J^{ki}}{\partial z^\ell} + J^{k\ell}\frac{\partial J^{ij}}{\partial z^\ell}  \equiv 0\,,
\label{jacobi0c}
\end{equation}
a quantity that is checked in practice.

When $\det J\neq 0$, an old theorem of Darboux, based on the algebraic properties of the Poisson bracket, implies that there is a coordinate change from the noncanonical coordinates $z$ to a set of canonically conjugate coordinates.  However, when $\det J=0$ this is not possible because of degeneracy, i.e., the existence of special functions $C$, called Casimir invariants,  that satisfy $\{C,f\}=0$ for all phase space functions $f$.   Thus, Casimir invariants are built-in to the phase space, for they will be conserved by a system generated by any Hamiltonian function. In the coordinates $z$, $C$ is a Casimir invariant if it satisfies
\begin{equation}
J^{ij}\frac{\partial C}{\partial z^j}=0\,.
\label{cascond}
\end{equation}

Noncanonical Hamiltonian systems possess the rich geometrical structure of so-called Poisson manifolds, where through every point of the phase space manifold is a conserved canonical Hamiltonian subspace, i.e., the manifold is foliated by symplectic leaves (see, e.g., the seminal reference \cite{weinstein} and the recent contribution \cite{pjmY20}).   We will see in practical terms how an interesting  Poisson  manifold emerges from the MK system.

For three-dimensional systems, like the MK system,  the Poisson tensor  $J$ has the form
\begin{equation}
J =
 \begin{bmatrix}
0& V_3& -V_2\\
-V_3& 0&  V_1\\
V_2& - V_1& 0\\
\end{bmatrix}
\end{equation}
for some vector $\mathbf{V}(z)=(V_1,V_2,V_3)$,  and it can be  shown easily that \eqref{jacobi0},  for the Jacobi identity,  is satisfied if  $\mathbf{V}(z)$ satisfies
\begin{equation}
\mathbf{V}\cdot \nabla \times \mathbf{V} =0\,.
\label{jacobi1}
\end{equation}
Thus, for three-dimensional systems there is a  convenient way to check the Jacobi identity.  Because antisymmetric matrices have even rank, $J$ must have rank 2 or 0.  The latter of course would generate trivial dynamics -- of interest is the case of rank 2 where there is a single Casimir invariant, and condition \eqref{cascond} can be written compactly as
\begin{equation}
\mathbf{V}\times \nabla C=0\,.
\label{3dcas}
\end{equation}
Relations \eqref{jacobi1}  and \eqref{3dcas} play central roles in our discovery of the Hamiltonian structure of the MK system.

\subsection{The Hamiltonian and Poisson bracket}
\label{ssec:mkham}

Given the equations of motion of a system, like \eqref{dotde}, \eqref{dots}, and \eqref{dotka}, and a constant of motion, one can seek a Poisson tensor by matching to the equations of motion while enforcing the Jacobi identity.   Thus we seek a suitable invariant, one that physically we expect to be an energy-like quantity.   Using the structure of the MK system  and some insight we find   \eqref{dotde}, \eqref{dots}, and \eqref{dotka} conserve the following:
 \begin{equation}
H=\frac1{\delta^2}\left[\ln\left(\frac{s}{\delta}\right) +\beta_1-\frac12\right]= \frac1{\delta^2}\includegraphics[height=5cm]{Hcloseup.pdf}[ \Lambda - 1/2\includegraphics[height=5cm]{Hcloseup.pdf}]\,,
 \label{H}
\end{equation}
 which can be shown directly. From \eqref{H} we obtain
 \begin{equation}
 s=\delta\, e^{\delta^2 {H} -\beta_1 +1/2}\,.
 \label{sde}
 \end{equation}
 Again, as with \eqref{critLa}, the threshold  for inequality \eqref{ineq} occurs when $\delta=s$, yielding
  \begin{equation}
 H > (\beta_1-1/2)/\delta^2=-0.0583/\delta^2\,,
 \label{critH}
 \end{equation}
 where again the gaussian value of $\beta_1$ is used in the  equality.

In terms of the  coordinates $z=(\delta, s, \kappa)$,  the analog of equation \eqref{genHamform} for the MK system takes the form
\begin{equation}
 \begin{bmatrix}
\dot{\delta}  \\
\dot{s} \\
\dot{\kappa}\\
\end{bmatrix}
=
 \begin{bmatrix}
0& V_3& -V_2\\
-V_3& 0&  V_1\\
V_2& - V_1& 0\\
\end{bmatrix}
 \begin{bmatrix}
\partial H/\partial \delta \\
\partial H/\partial s \\
\partial H/\partial \kappa \\
\end{bmatrix}\,,
\label{form}
\end{equation}
which upon making use of the MK equations of motion  \eqref{dotde}, \eqref{dots}, and \eqref{dotka}, and the candidate Hamiltonian \eqref{H}, we obtain the following equation for the Poisson tensor:
\begin{equation}
 \begin{bmatrix}
-{c_2}\, {\delta \kappa}/{(2s)} \\
- c_2\,  \kappa \, \Lambda  \\
c_1\,  {\kappa}/{s^2}\\
\end{bmatrix}
=
 \begin{bmatrix}
0& V_3& -V_2\\
-V_3& 0&  V_1\\
V_2& - V_1& 0\\
\end{bmatrix}
 \begin{bmatrix}
-2 \Lambda  /\delta^3\\
1/( s \delta^2)\\
0 \\
\end{bmatrix}
\label{form20}\,.
\end{equation}
Thus, the goal is to solve \eqref{form20} for a $\mathbf{V}=(V_1,V_2,V_3)$ that satisfies \eqref{jacobi1}.

It follows immediately from \eqref{form20} that
\begin{equation}
\dot{\delta}=-{c_2}\, {\delta \kappa}/{(2s)} =V_3/(s\delta^2)
\quadd\mathrm{and}\quadd
\dot{s}= - c_2\,  \kappa  \Lambda   = 2 V_3  \Lambda  /\delta^3\,;
\label{V3a}
\end{equation}
thus,
\begin{equation}
V_3= -c_2 \kappa \delta^3/2
\label{V3}
\end{equation}
works for both equations  of \eqref{V3a}.  The remaining equation yields the expression,
\begin{equation}
\dot{\kappa}= c_1\kappa/s^2= -2 V_2 \Lambda  /\delta^3  - V_1/(s\delta^2)\,.
\end{equation}
Thus, it appears there is  freedom in the choices of $V_1$ and $V_2$ to satisfy \eqref{jacobi1},
the Jacobi identity.  Upon setting
\begin{equation}
V_1 = -c_1\kappa\delta^2 /s -2 V_2 \Lambda  s/\delta \,.
\label{V1}
\end{equation}
we seek to find a $V_2$ that ensures   \eqref{jacobi1} is satisfied.  A direct calculation implies
 \bal
\mathbf{V}\cdot\nabla\times \mathbf{V}
&= - \frac{c_2}{2}\,  \kappa \delta^3 \, \frac{\partial  V_2 }{\partial \delta}
-  c_2\,  \kappa s\delta^2 \Lambda   \,  \frac{\partial V_2 }{\partial s}
+  c_1\,  \frac{\kappa\delta^2}{s}  \, \frac{\partial  V_2}{\partial \kappa}
\nonumber\\
& + V_2\left(-c_1\, \frac{\delta^2}{s} + \frac{c_2}{2}\, \kappa \delta^2
-c_2\,  \kappa \delta^2   \Lambda
\right)
+   \frac{c_1 c_2}{2}\, \frac{\delta^5 \kappa^2}{s^2}\,.
\label{jacobi2}
\eal
and the goal is to find a $V_2$ such that \eqref{jacobi2} vanishes.

Upon  inserting the following into \eqref{jacobi2},
\begin{equation}
V_2= \kappa\delta A(x)\,,
\label{V2}
\end{equation}
where recall $x:=s/\delta$ and  using
\begin{equation}
\partial_k V_2= \delta A(x)\,, \qquad \partial_s V_2 = \kappa A'(x)\,,
  \qquad \partial_\delta V_2= \kappa  A(x) - \kappa x A'(x)
\end{equation}
\eqref{jacobi2} becomes
\begin{equation}
\mathbf{V}\cdot\nabla\times\mathbf{V}=-c_2\kappa^2 \delta^3\left[x A' \left(\Lambda -\frac12\right) + \Lambda A -\frac{c_1}{2x^2}\right]\,.
\end{equation}
Thus,  the Jacobi identity is satisfied provided  we can find a solution $A(x)$ to
\begin{equation}
x A' \left(\Lambda -\frac12\right) + \Lambda A -\frac{c_1}{2x^2}=0\,,
\label{Aeq}
\end{equation}
where recall $\Lambda (x)$ is given by \eqref{La}.   If such a function $A$ is found, then the MK system is Hamiltonian with a bracket defined by
\begin{equation}
V_1= -c_1\frac{\kappa\delta^2}{s} -2 \kappa \Lambda  s \, A(x)  \,,\quad
V_2=\kappa\delta \,A(x)\,,
\quad\mathrm{and}\quad
V_3= -\frac{c_2}{2} \kappa \delta^3\,.
\label{Vs}
\end{equation}

Thus, we proceed to solve for $A$.   In the light of the inequality \eqref{ineq},   $\Lambda $ is seen to be a positive monotonic function with the  inverse
\begin{equation}
x=e^{\Lambda -\beta_1}\,.
\label{Lainv}
\end{equation}
Consequently, we can use   $\Lambda $ as the independent variable and  rewrite \eqref{Aeq} as
\begin{equation}
\frac{d}{d\Lambda } \left( e^\Lambda \sqrt{\Lambda -1/2\,}\, A\right) = \frac{c_1}{2x^2} \frac{e^\Lambda }{\sqrt{\Lambda -1/2\,}}
 =\frac{c_1e^{2\beta_1}}{2} \frac{e^{-\Lambda }}{\sqrt{\Lambda -1/2\,}}\,.
 \label{Lageq12}
 \end{equation}
 Here we have assumed $\Lambda -1/2\geq 0$.  In light of  \eqref{critLa} this may not be true; thus, we  will return and consider the case where $\Lambda -1/2\leq 0$.

 Antidifferentiating both sides gives
\begin{equation}
A=\frac{e^{-\Lambda }}{\sqrt{\Lambda -1/2\,}} \frac{c_1e^{2\beta_1-1/2}}{2 }  \int^{\Lambda -1/2}  \frac{e^{-\ell}}{\sqrt{\ell\,} }\, \,d\ell\,.
\end{equation}
Then, with the definition of  the incomplete gamma function
\begin{equation}
\gamma(\sigma,u)=\int_0^u  u'^{\sigma-1} e^{-u'} du'\,,
\label{ligf}
\end{equation}
and the identity
\begin{equation}
\gamma(1/2,u)=\sqrt{\pi} \mathrm{erf}(\sqrt{u})=2 \int_0^{\sqrt{u}} e^{-u'^2} d u'\,,
\label{gamaerf}
\end{equation}
with erf being the error function, we obtain

\begin{equation}
A=c_0\,  \frac{e^{-\Lambda }}{\sqrt{\Lambda -1/2\,}}
+  \frac{\sqrt{\pi} c_{1} \,  }{2}\,  e^{2\beta_1-1/2} \, \frac{e^{-\Lambda }}{\sqrt{\Lambda -1/2\,}}\,  \mathrm{erf}\left(\sqrt{\Lambda -1/2}\right)
\,,
 \label{A1}
\end{equation}
where $c_0\in\mathbb{R}$ is the integration constant.   Upon inserting \eqref{A1} into the expressions of  \eqref{Vs}, we have a  one-parameter family of Poisson brackets.
It is convenient to choose $c_0=0$,  giving
\begin{equation}
A=\frac{\sqrt{\pi} c_{1} \,  }{2\, x}\, \frac{e^{\beta_1-1/2} }{\sqrt{\Lambda -1/2\,}}\,  \mathrm{erf}\left(\sqrt{\Lambda -1/2}\right)
\,,
\label{Af}
\end{equation}
which is valid for $\Lambda -1/2\geq 0$, but we will see it is also valid for $\Lambda -1/2<0$. Observe, the choice $c_0=0$  gives us regularity at $\Lambda =1/2$.   Inserting \eqref{Af} into the equations of \eqref{Vs} defines  the  noncanonical Poisson bracket that we will use  for $\Lambda  \geq 1/2$.

Now consider the case where  $\Lambda  \leq 1/2$, which occurs for
\begin{equation}
0\leq 1/2- \Lambda \leq 1/2-\beta_1 - \ln x\leq 1/2-\beta_1 \,,
\label{Larange}
\end{equation}
where we assume $x\geq 1$ consistent again with \eqref{ineq}.  Instead of \eqref{Lageq12}, consider
\begin{equation}
\frac{d}{d\Lambda } \left( e^\Lambda \sqrt{1/2-\Lambda \,}\, A\right) = - \frac{c_1 e^{2\beta_1}}{2} \,\frac{e^{-\Lambda }}{\sqrt{1/2-\Lambda } }\,,
 \label{L leq12}
 \end{equation}
where \eqref{Aeq} has been used.   Integrating both sides   from $\Lambda =1/2$  leads to
  \begin{equation}
 A=   {c_1 \,e^{2\beta_1 -1/2}\, e^{-\Lambda }}  \, \,  M(1/2,3/2,1/2-\Lambda )\,, \qquad  \Lambda  \leq 1/2\,,
 \label{Af2}
 \end{equation}
 where
 \begin{equation}
 M(1/2,3/2,z)= \frac12\int _0^1 \frac{e^{zt}}{\sqrt{t}}\, dt
 = \frac{1}{2\sqrt{z}}\int _0^z \frac{e^{u}}{\sqrt{u}}\, du\,,
 \end{equation}
is the   Kummer function \cite{a&s},  which is  sometimes  called the confluent hypergeometric function of the first kind and denoted by ${}_1F_1(1/2;3/2;z)$.  Inserting \eqref{Af2} into the equations of \eqref{Vs} defines  the  noncanonical Poisson bracket defined  for $\beta_1\leq \Lambda  \leq 1/2$.

A comparison of \eqref{Af} and \eqref{Af2} follows from the identity
\begin{equation}
\mathrm{erf}(z)= \frac{2z}{\sqrt{\pi}}\,M(1/2,3/2,-z^2)\,,
\end{equation}
which implies
\begin{equation}
\frac{\sqrt{\pi}}{2}\frac{1}{\sqrt{\Lambda -1/2}}\, \mathrm{erf}\left(\sqrt{\Lambda -1/2}\right) = M(1/2,3/2,1/2-\Lambda ))\,.
\end{equation}
Because  $\mathrm{erf}(ix)$ for $x\in\mathbb{R}$  is pure imaginary,  like $\sqrt{\Lambda -1/2}$ for $\Lambda <1/2$, we can analytically continue the expression in terms of the error function through zero  to pure imaginary values yielding a real quantity.  Thus, expression \eqref{Af} can be used for both  $\Lambda \geq 0$ and $\Lambda <0$.

\subsection{The Casimir Invariant}
\label{ssec:cas}

Given the Hamiltonian structure of Sec.~\ref{ssec:mkham} we know immediately that the MK system  must possess  another constant of motion, the Casimir invariant.  A little thought reveals
the  Casimir must have the form
\begin{equation}
C= \frac{c_2}{2} \kappa + \frac{\lambda(x)}{\delta}\,.
\label{C}
\end{equation}
In the present context,  \eqref{3dcas} is equivalent to
\begin{equation}
 \begin{bmatrix}
0& V_3& -V_2\\
-V_3& 0&  V_1\\
V_2& - V_1& 0\\
\end{bmatrix}
 \begin{bmatrix}
\partial C/\partial \delta \\
\partial C/\partial s \\
\partial C/\partial \kappa \\
\end{bmatrix}
=
\begin{bmatrix}
0& V_3& -V_2\\
-V_3& 0&  V_1\\
V_2& - V_1& 0\\
\end{bmatrix}
 \begin{bmatrix}
-\lambda/\delta^2-s\lambda'/\delta^3 \\
\lambda'/\delta^2\\
c_2/2 \\
\end{bmatrix}
=0
\label{cas}
\end{equation}
The first equation  gives
\begin{equation}
0= -\frac{c_2}{2} \kappa  {\delta^3}\frac{\lambda'}{\delta^2} - \kappa\delta A \frac{c_2}{2}
\quad\Rightarrow\quad\lambda'=-A\,,
\end{equation}
while the second gives after some manipulation
\begin{equation}
0=-\frac{c_2\kappa\delta}{2}\left(\lambda - xA +\frac{c_1}{x}+2Ax\Lambda
\right)\,.
\end{equation}
Because of the Jacobi identity,  the third equation must automatically  be solved by the above, which can be  verified directly.
Thus we have
\bal
\lambda'(x)&=-A(x)
\label{laprime}
\\
\lambda(x)&= Ax(1-2\Lambda )-\frac{c_1}{x}=-2Ax(\Lambda -1/2) - \frac{c_1}{x} \,,
\label{lam}
\eal
where $A$ solves \eqref{Aeq}.   Upon differentiating \eqref{lam} we see that if $A$ satisfies \eqref{Aeq}, then \eqref{laprime} is automatic for a solution of   \eqref{lam}. In particular, using \eqref{Af} we obtain
\begin{equation}
 \lambda(x) = -  \sqrt{\pi} c_{1} \, e^{\beta_1-1/2} \, {\sqrt{\Lambda -1/2\,}}\,  \mathrm{erf}\left(\sqrt{\Lambda -1/2}\right)
 - \frac{c_1}{x} \,,
 \label{laf}
\end{equation}
where use has been made of  $xe^{-\Lambda }= e^{-\beta_1}$. Figure \ref{fig:Alambda} displays plots of both $A/c_1$ of \eqref{Af} (upper)
and $\lambda/c_1$ \eqref{laf} (lower) as  functions of $x$,  beginning with its critical value of $x=s/\delta=1$, revealing that both are  relatively simple monotonic functions consistent with inequality \eqref{ineq}.
\begin{figure}[htb]
\includegraphics[scale=.5]{ALambda.pdf}
\caption[sm]{Plots of $A/c_1$ (blue) and  $\lambda/c_1$ (orange) vs.\ $x=s/\delta$.  Recall $\lambda'=-A(x)$ [cf.\ Eq.~\eqref{laprime}].}
\label{fig:Alambda}
\end{figure}

Thus to summarize,  the Casimir, a second invariant, is given by
\bal
C&= \frac{c_2}{2} \kappa + \frac{\lambda(x)}{\delta}
\nonumber\\
&=  \frac{c_2}{2} \kappa   -  \sqrt{\pi} c_{1} \, e^{\beta_1-1/2} \, \frac{1}{\delta} {\sqrt{\Lambda -1/2\,}}\,  \mathrm{erf}\left(\sqrt{\Lambda -1/2}\right)
 - \frac{c_1}{s}\,,
\label{Cfinal}
\eal
where we have inserted  $\lambda$ from \eqref{laf} in \eqref{Cfinal}  and used $x=s/\delta$.  We remind the reader that \eqref{Cfinal} is valid for $\Lambda <1/2$ as well as for $\Lambda \geq1/2$.

Because of inequality \eqref{ineq}, not all values of \eqref{Cfinal} are permissible.  To understand the permissible range, we rewrite  \eqref{Cfinal}  as follows:
\begin{equation}
\frac{s}{c_1}C-\frac{c_2}{2c_1}\kappa s= F(\Lambda)\,,
\end{equation}
where
\begin{equation}
F(\Lambda )=  \sqrt{\pi}   \, e^{\Lambda -1/2} \,  {\sqrt{\Lambda -1/2\,}}\,  \mathrm{erf}\left(\sqrt{\Lambda -1/2}\right)   + 1 >0\,.
\label{FLa }
\end{equation}
As noted in \eqref{critLa},   the smallest allowable value of $\Lambda $ is $\beta_1$, and it is not hard to show that $F(\Lambda )$ obtains its minimum value at $\beta_1$, at  which it is positive.  Thus we see, the Casimir must satisfy
\begin{equation}
C<\frac{c_2}{2} \kappa
\label{Cineq}
\end{equation}
to be consistent with the inequality \eqref{ineq}.

\section{The nature of the solution, reduction to quadrature, and analysis}
\label{sec:quad}

\subsection{Geometrical solution}

Given that we have a three-dimensional system with two constants of motion, the solution space can be visualized by examining the  intersection of the level sets  of the  Hamiltonian $H$ of \eqref{H} with those of the Casimir $C$ of \eqref{Cfinal}.   This is a direct analog of how the stable and unstable trajectories of the  the free rigid body, as governed by Euler's equations,  are understood in terms of  the intersection of the angular momentum spheres with the energy ellipsoids.  The same situation occurs for other noncanonical Hamiltonian systems such as the Kida problem of fluid mechanics \cite{pjmMF97} and  the rattleback toy \cite{pjmYT17}, and indeed a large class of flows on Poisson manifolds \cite{pjmY20}.  Thus, by plotting level sets of  \eqref{H} and \eqref{Cfinal}  the nature of trajectories is revealed and, in addition, one can delineate the accessible phase space consistent with \eqref{ineq}.

 Figure~\ref{fig:Hcontours} displays contours of the Hamiltonian $H$ of \eqref{H}. Since $H$ is independent of the variable $\kappa$ it has  translational   symmetry along the $\kappa$ axis and, as can be seen in Fig.~\ref{fig:Hcontoursa}, has a sheet-like topology for positive  values of the variables.  Level sets corresponding to different  signs of $H$ have opposite curvature in the $s-\delta$ plane,  with there being a region of negative values of $H$ consistent with \eqref{critH} and \eqref{ineq}.   This `negaive energy' interval is shown in Fig.~\ref{fig:Hcontoursb} for $\delta<1$.  For all values of $H$,  the sheets become tangent to the $H=0$ plane with a slope given by $s= \delta e^{1/2-\beta_1}$.   Consequently, if $\delta$ and $s$ approach zero, they  do so in a clear and universal way independent of $\mathring{H}$,  the initial  value of the Hamiltonian constant of motion.

 \begin{figure}[htb]
\centering
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.43\textwidth]{Hcontours2.pdf}
\label{fig:Hcontoursa}
}
\hspace{.5cm}
\subfigure[{\footnotesize \  }]{
\def\includegraphics[height=5cm]{Hcloseup.pdf}{\includegraphics[height=5cm]{Hcloseup.pdf}}
\def\includegraphics[height=1.5cm]{des.pdf}{\includegraphics[height=1.5cm]{des.pdf}}
\defr{r}
\bottominset{\includegraphics[height=1.5cm]{des.pdf}}{\includegraphics[height=5cm]{Hcloseup.pdf}}{20pt}{5pt}
\label{fig:Hcontoursb}
}
\caption[sm]
{Contour plots of the Hamiltonian $H$ of \eqref{H}. (a) $\mathring{H}={-.50}$ (orange), $\mathring{H}={0}$ (blue), $\mathring{H}={1}$ (green), $\mathring{H}={10}$ (magenta), and $\mathring{H}={100}$ (lavender).  (b) Closeup plots of $s$ vs.\ $\delta$ at arbitrary $\kappa$.  Note the  region to the left of  $\delta=1$ for the case where $\mathring{H}=-0.0583$ (orange)  compared to the case $s=\delta$ (blue) with the inset showing the difference  (orange-blue), showing that  $s>\delta$ to the left.  The other curve  with   $\mathring{H}=0.005$ (green) is added for comparison.
}
\label{fig:Hcontours}
\end{figure}

In Figs.~\ref{fig:Ccontours4},
   \ref{fig:Ccontours9},   and \ref{fig:Ccontours90}
we show two views of contour plots of the Casimir of \eqref{Cfinal} for three values of the  vortex ring tilt  angle $\alpha\in\{\pi/4, \pi/6,\pi/90\}$.
Observe that these surfaces are again sheet like but with more interesting structure, no longer having the $\kappa$ independence of the $H$-surfaces.
\begin{figure}[htb]
\centering
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.45\textwidth]{CaPi-4a2.pdf}
\label{fig:Ccontours4a}
}
\hspace{.5cm}
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.45\textwidth]{CaPi-4b2.pdf}
\label{fig:Ccontours4b}
}
\caption[sm]
{(a) View of the level sets of the  Casimir $C$ for $\alpha=\pi/4$. (b)    Same as (a) but rotated.
}
\label{fig:Ccontours4}
\end{figure}

\begin{figure}[htb]
\centering
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.44\textwidth]{CaPi-9a2.pdf}
\label{fig:Ccontours9a}
}
\hspace{.5cm}
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.48\textwidth]{CaPi-9b2.pdf}
\label{fig:Ccontours9b}
}
\caption[sm]
{(a) View of the level sets of the  Casimir $C$ for $\alpha=\pi/9$. (b)    Same as (a) but rotated.
}
\label{fig:Ccontours9}
\end{figure}

\begin{figure}[htb]
\centering
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.45\textwidth]{CaPi-90a3.pdf}
\label{fig:Ccontours90a}
}
\hspace{.5cm}
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.46\textwidth]{CaPi-90b4.pdf}
\label{fig:Ccontours90b}
}
\caption[sm]
{(a) View of the level sets of the  Casimir $C$ for $\alpha=\pi/90$. (b)    Same as (a) but rotated.
}
\label{fig:Ccontours90}
\end{figure}

 The case where $\alpha=\pi/90$ corresponds to nearly parallel vortices.  For this case  $c_1\approx 5\times 10^{-5}<<1$ and \eqref{dotka0} implies $\kappa$ is nearly constant.  For this case the trajectory  lingers in the flat regions of Fig.~\ref{fig:Ccontours90} and the dynamics is approximately governed  by
\bal
   \dot{\delta}&=-\frac{c_2}{2}\,   \frac{\mathring{\kappa} \delta }{s}
  \label{dotde00}\\
 \dot{s}&=- c_2\, \mathring{\kappa}  \left[ \ln\left(\frac{s}{\delta}\right) +\beta_1\right]\,,
   \label{dots00}
\eal
with $c_2\approx 0.0796$.  Because $H$ is independent of $\kappa$, it is not a surprise that \eqref{H}  is still conserved by \eqref{dotde00}  and \eqref{dots00}.  This leads to the quadrature discussed in Appendix \ref{sec:kac}, where the `$\kappa$-clock' is proportional to ordinary time. As $s$ gets small,   $\kappa$  becomes activated,  demonstrating the importance of the  local induction velocity
for developing curvature.

Now consider the roles played by  the $H$ and $C$ surfaces. As a specific example, consider the case with the initial conditions
\begin{equation}
\mathring{\delta}=0.01 <  \mathring{s}=0.10 <  \mathring{\kappa}=1\,.
\label{SIC}
\end{equation}
With the values of \eqref{SIC} and a choice for the vortex ring tilt angle $\alpha$,    the solution lies on the intersection of the level sets with
\begin{equation}
\mathring{H}= 22442.9 \quadd\mathrm{and}\quadd \mathring{C}=-6.40628\,,\qquad \mathrm{for}\quad \alpha=\pi/9\,,
\end{equation}
where again we use $\beta_1=0.4417$.  Thus, these  initial conditions and corresponding initial values of
$\mathring{H}$ and $\mathring{C}$ are  consistent wtih  \eqref{ineq}.   This particular intersection is displayed in Fig.~\ref{fig:HC9_SIC}.
\begin{figure}[htb]
\centering
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.47\textwidth]{HC9-SIC-100-2.pdf}
\label{fig:SIC_100}
}
\hspace{.5cm}
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.47\textwidth]{HC9-SIC-1000-2.pdf}
\label{fig:SIC_1000}
}
\caption[sm]
{Plot of the intersection of the level set $\mathring{H}= 22442.9$ (orange) with that of  $\mathring{C}=-6.40628$ (blue) for $\alpha=\pi/9$.  These values correspond to the initial conditions of \eqref{SIC}.  In  (a) $\kappa$ ranges to 100, while in (b) we zoom out to $\kappa$ ranging to 1000  to give a more global perspective on the shape of the surfaces. }
\label{fig:HC9_SIC}
\end{figure}
Observe how the curve of intersection approaches increasingly  large values of $\kappa$ as both $\delta$ and $s$ approach zero.

In Fig.~\ref{fig:HCrange} we plot multiple intersections of the $H$ and $C$ contours.
\begin{figure}[htb]
\centering
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.47\textwidth]{C0H5,15,20,30,45k-2.pdf}
\label{fig:C0H}
}
\hspace{.5cm}
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.47\textwidth]{H0C-10--4-3.pdf}
\label{fig:H0C}
}
\caption[sm]
{Intersections of the level sets of $H$ and $C$  for the initial conditions of \ref{SIC} with $\alpha=\pi/9$.  In (a) we set
$\mathring{C}= -6.40628$ and show contours for $\mathring{H}\in\{5,000; 15,000; 20,100; 30,000;45,000\}$, while in (b) we set  $\mathring{H}=22442.9$ and show contours for $\mathring{C}\in\{-10,-9,-8,-7,-6,-4\}$.
}
\label{fig:HCrange}
\end{figure}
From this figure we see how $\kappa$ diverges for a variety of initial conditions,  behavior that is in fact generic.

Yet another picture of the singularity emerges if we eliminate $\delta$ between \eqref{H} and \eqref{Cfinal}, giving the expression
\bal
\frac{c_2\kappa}{2C}  &= 1- \frac{\sqrt{H}}{C} \frac{\lambda(x)}{\sqrt{\Lambda (x)-1/2}}
\nonumber\\
&=  1+  \frac{c_1\sqrt{H}/{C}} {\sqrt{\Lambda (x)-1/2}}\left(
 \sqrt{\pi} \, e^{\beta_1-1/2} \, {\sqrt{\Lambda -1/2\,}}\,  \mathrm{erf}\left(\sqrt{\Lambda -1/2}\right)
 + \frac{1}{x}
 \right)
 \nonumber\\
 &=:1+ \frac{c_1\sqrt{H}}{C}\, K[x]\,.
 \label{Keq}
\eal
This formula is well defined for $H<0$ because  $\delta= \sqrt{H/(\Lambda -1/2)}>0$.   Figure \ref{fig:K} shows the function $K$ with the singularity occurring at $x^* \approx 1.0600$.
\begin{figure}[htb]
\includegraphics[scale=.4]{K.pdf}
\caption[sm]{Plots of  the function $K$ of \eqref{Keq} vs.\ $x$, depicting the singularity that occurs for $\kappa$ at $x^* \approx 1.0600$.}
\label{fig:K}
\end{figure}

The reader may wonder what would happen to the Casimir  if we retained  the term of \eqref{A1} with  the integration constant $c_0$.  It turns out that this  merely adds a term proportional to $\sqrt{H}$ to the Casimir and thus just shifts the value of $C$,  and thus has no consequence other than changing the numerical value of $C$ for the same plots.

\subsection{Reduction to quadrature}

Let us now consider exact integration of the system.  Suppose $(\mathring{\delta},\mathring{s}, \mathring{\kappa})$ represent an arbitrary initial condition at time $\mathring{t}$ with corresponding initial values of our two invariants $\mathring{H} = H(\mathring{\delta},\mathring{s}, \mathring{\kappa})$ and $\mathring{C} = C(\mathring{\delta},\mathring{s}, \mathring{\kappa})$.  Because $H =\mathring{H}$ for all time,  \eqref{H} can be used in \eqref{Cfinal}  to obtain
\begin{equation}
\mathring{C}=\frac{c_2}{2}\kappa - c_1 e^{\beta_1-1/2} \sqrt{\pi\mathring{H}} \, \,
  \mathrm{erf}\left(\delta\sqrt{\mathring{H}}\right)    - \frac{c_1}{s} \,,
  \label{C1}
\end{equation}
where recall  this formula is well-defined when $\mathring{H}<0$ because it  analytically continues to the expression of \eqref{Af2}.   Next,  using \eqref{sde} to eliminate $s$  we  obtain
\begin{equation}
\mathring{C}=\frac{c_2}{2}\kappa - c_1 e^{\beta_1-1/2} \sqrt{\pi\mathring{H}} \, \,
  \mathrm{erf}\left(\delta\sqrt{\mathring{H}}\right)    - \frac{c_1}{\delta} \,  e^{-\delta^2\mathring{H} +\beta_1 - 1/2}\,.
  \label{intermed}
\end{equation}
Solving \eqref{intermed} for $\kappa$ is immediate
\begin{equation}
\frac{c_2}{2}\kappa= \mathring{C}  + c_1 e^{\beta_1-1/2} \sqrt{\pi\mathring{H}} \, \,
  \mathrm{erf}\left(\delta\sqrt{\mathring{H}}\right)   +  \frac{c_1}{\delta} \,  e^{-\delta^2\mathring{H} +\beta_1 - 1/2}\,.
  \label{kade}
\end{equation}
With  \eqref{sde} and \eqref{kade}, we obtain from   \eqref{dotde}, the following:
\bal
\dot{\delta}&=-\frac{c_2}{2}\frac{\kappa \delta }{s} =-
\frac{c_2}{2}\kappa \, e^{-\delta^2\mathring{H} +\beta_1 -1/2}
\nonumber\\
&=- \left( \mathring{C}  + c_1 e^{\beta_1-1/2} \sqrt{\pi\mathring{H}} \, \,
  \mathrm{erf}\left(\delta\sqrt{\mathring{H}}\right)   +  \frac{c_1}{\delta} \,  e^{-\delta^2\mathring{H} +\beta_1 - 1/2}
\right) \, e^{-\delta^2\mathring{H} +\beta_1 -1/2}\,,
\label{ddede}
\eal
which leads immediately to  the quadrature
\begin{equation}
\mathring{t} - t= e^{-2 \beta_1 +1}  \,\,\int_{\mathring{\delta}}^\delta
 \frac{ e^{\delta'^2\mathring{H}}\, \,\delta'\,   d \delta'}
{\delta'\,e^{-\beta_1 +1/2}\mathring{C}   + c_1  \sqrt{\pi\mathring{H}}\delta' \, \,
  \mathrm{erf}\left(\delta'\sqrt{\mathring{H}}\right)   +    c_1 e^{-\delta'^2\mathring{H}}
} \,.
\label{dequad}
\end{equation}
Although unwieldy, integration of \eqref{dequad} gives $\delta$ as a function of time, and  via  \eqref{sde} and \eqref{kade}  we obtain $s$ and $\kappa$ as functions of time, with the latter diverging.

 From \eqref{dequad} we see that if there is a finite-time singularity where $\delta\rightarrow 0$, then it occurs at a time $t_\infty$ in accordance with the following formula:
\begin{equation}
t_\infty=\mathring{t} +  e^{-2 \beta_1 +1}  \,\,\int^{\mathring{\delta}}_0
 \frac{ e^{\delta'^2\mathring{H}}\, \,\delta'\,   d \delta'}
{\delta'\,e^{-\beta_1 +1/2}\mathring{C}  +  c_1  \sqrt{\pi\mathring{H}}\delta' \, \,
  \mathrm{erf}\left(\delta'\sqrt{\mathring{H}}\right)   +   c_1  e^{-\delta'^2\mathring{H}}
} \,.
\label{tinfty}
\end{equation}
For $\mathring{H}\neq 0$,  we see from \eqref{H} that at  a singularity where $\delta\rightarrow 0$, we must have $\Lambda -1/2\rightarrow 0$.  As noted in \eqref{critH},
 the smallest  initial value of $H$ that satisfies   \eqref{ineq} is $\mathring{H}^*= (\beta_1-1/2)/\mathring{\delta}^2$, which is negative;
 if  $\mathring{H}> \mathring{H}^*$, then
\eqref{ineq} is satisfied initially.  If $\delta$ decreases  and it initially satisfies   $s>\delta$,  then it will satisfy it throughout its evolution  because  $\delta=s\, e^{-\delta^2\mathring{H} +\beta_1 -1/2}$.

\subsection{Exact  solutions with Leray scaling}
\label{ssec:exleray}

Evidently, the integral of \eqref{dequad} is dramatically simplified with the choice  $\mathring{H}=0$.
With this choice,  \eqref{sde} implies
\begin{equation}
s=\delta e^{-\beta_1 + 1/2}
\label{salde}
\end{equation}
and equations \eqref{dotde} and \eqref{dots} become identical upon setting
\[
\Lambda  =\ln(s/\delta) +\beta_1=1/2\,,
\]
a choice  consistent with  the inequality of \eqref{ineq}.
Given that for gaussian core profiles $\beta_1=0.4417$  (see  \cite{MK19}), we have
\begin{equation}
  e^{-\beta_1 + 1/2}  = 1.0600\quad\Rightarrow \quad s>\delta \,.
\end{equation}
Thus inequality \eqref{ineq} is satisfied, although only barely.   Proceeding, the  MK system reduces for $\mathring{H}=0$  to
 \begin{equation}
 \dot{s}=- \frac{c_2}{2}\,  \kappa
\quadd\mathrm{and}\quadd
 \dot{\kappa}=c_1\,  \frac{\kappa}{s^2}\,,
  \label{nudotka}
\end{equation}
with the Casimir becoming
\begin{equation}
C= \frac{c_2}{2} \kappa-\frac{c_1}{s}\,,
\label{CH0}
\end{equation}
which follows from \eqref{C1}.

Before considering arbitrary initial $\mathring{C}$, we consider the easily tractable case where both  $\mathring{H}=0$ and $\mathring{C}=0$.  This   implies
\begin{equation}
s\kappa= 2c_1/c_2 = 2\sin\alpha\,.
\label{ska}
\end{equation}
Therefore if we choose
\begin{equation}
2\sin\alpha <1\qquad \mathrm{or}\qquad \alpha<\frac{\pi}{6}\,,
\end{equation}
then the other part of the  inequality   \eqref{ineq} is satisfied, viz.\  $s<1/\kappa$.
 Note, the value of $\alpha=\pi/4$ chosen in \cite{MK19,MK19b} does not satisfy this inequality.
Using \eqref{ska} in the  $\dot{\kappa}$ equation of  \eqref{nudotka} gives
\begin{equation}
\dot{\kappa}=\frac{c_2^2}{4 c_1}  \kappa^3 \,,
\end{equation}
which is easily integrated to obtain
\begin{equation}
- \frac{1}{2\kappa^2} +  \frac{1}{2\mathring{\kappa}^2}  = \frac{c_2^2}{4 c_1}(t-\mathring{t})
\end{equation}
and  the exact solution
\begin{equation}
\kappa^{-1} = \sqrt{\frac{c_2^2}{2c_1}(t_\infty-t)}= \sqrt{\frac{\mathrm{cot}\alpha}{8\pi}(t_\infty-t)}\,,
\label{lerayE}
\end{equation}
where
\begin{equation}
t_{\infty}=
\mathring{t} + \frac{2c_1}{c_2^2\mathring{\kappa}^2}
=\mathring{t} + \frac{8\pi \tan\alpha }{\mathring{\kappa}^2}
\,.
\label{tinfytC0}
\end{equation}
As expected, smaller values of $\kappa$ take longer to diverge.

The rest of the solution is obtained from \eqref{ska} and \eqref{salde}, i.e.,
\bal
s&=2 \kappa^{-1} \, \sin\alpha= 2\sin\alpha \sqrt{\frac{\mathrm{cot}\alpha}{8\pi}(t_\infty-t)}
\nonumber\\
&= \sqrt{\frac{\sin(2\alpha)}{4\pi}(t_\infty-t)}=\delta e^{1/2-\beta_1} \,.
\label{sdeleray}
\eal
 It is a simple matter to insert the solutions of \eqref{lerayE} and \eqref{sdeleray} into \eqref{dotde}, \eqref{dots}, and  \eqref{dotka} to verify directly that they are indeed an exact solution, one that satisfies the inequalities of \eqref{ineq}.

Now consider the more general case $\mathring{C}\neq 0$.   Solving \eqref{CH0} for $s$ and inserting into \eqref{nudotka} gives
\begin{equation}
\dot{\kappa}=\frac{\kappa}{c_1}  \left(\frac{c_2}{2}\kappa - \mathring{C} \right)^2
= \frac{c_2^2}{4c_1}\kappa \left(\kappa - 2\mathring{C}/c_2 \right)^2\,,
\end{equation}
which is easily integrated to obtain the exact solution,
\begin{equation}
  \ln\Big(1-\frac{2\mathring{C}/c_2}{\kappa}\Big) +  \frac{2 \mathring{C}/c_2 }{\kappa\, (1 -   \frac{2 \mathring{C}/c_2}{\kappa} )}
  - \ln\Big(1-\frac{2\mathring{C}/c_2}{\mathring{\kappa}}\Big) - \frac{2 \mathring{C}/c_2 }{\mathring{\kappa}\, (1 -   \frac{2 \mathring{C}/c_2}{\mathring{\kappa}})}
= - \frac{\mathring{C}^2}{c_1} \,(t - \mathring{t})\,.
\label{kasol}
\end{equation}
 Because the model only makes sense if $\delta, s$, and $\kappa$ are all greater than or equal to zero, we must have $c_1/s=\kappa{c_2}/{2} -C\geq 0$ or $1-2 {C}/( {\kappa} c_2)\geq 0$, which  in fact according to \eqref{Cineq}  is true for all allowable values of $\mathring{H}$.

Assuming the physical initial conditions  satisfy $\mathring{\kappa}>  2   \mathring{C}/c_2$, we see the $\dot\kappa>0$, so $\kappa$ continues to grow,    with divergence occurring at the finite time
\begin{equation}
t_{\infty}=\mathring{t} + \frac{c_1}{\mathring{C}^2}\left(
   \ln\Big(1-\frac{2\mathring{C}/c_2}{\mathring{\kappa}}\Big) +  \frac{2 \mathring{C}/c_2 }{\mathring{\kappa}\, (1 -   \frac{2 \mathring{C}/c_2}{\mathring{\kappa}} )}
\right)\,,
\label{blowupH0}
\end{equation}
an approximation to \eqref{tinfty} but a generalization to \eqref{tinfytC0}.  This is depicted in Fig.~\ref{fig:blowupH0},

\begin{figure}[htb]
\centering
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.45\textwidth]{blowupH02.pdf}
\label{fig:blowupH0}
}
\hspace{.5cm}
\subfigure[{\footnotesize \  }]{\includegraphics[width=0.45\textwidth]{blowup-again.pdf}
\label{fig:blowup_again}
}
\caption[sm]
{Plots of the blowup times for $ \mathring{H}=0$, $ \mathring{C}\neq 0$, (a) according to Eq.~\eqref{blowupH0} where   $\tau= (t_\infty- \mathring{t})\mathring{C}^2/c_1$  vs.\  $2 \mathring{C}/(c_2\mathring{\kappa})$ and  (b)  according to Eq.~\eqref{blowupH0} where    $\tau= (t_\infty- \mathring{t})\mathring{C}^2/c_1$ vs.\  $X={2\sin\alpha}/({\mathring{\kappa}\mathring{s}})$.
}
\label{F:blowups}
\end{figure}

The quantity $t_{\infty}$ can also be written as a function of $\kappa s$, which makes it convenient for assessing initial conditions compatible with \eqref{ineq}.  Using  \eqref{CH0} we obtain
\begin{equation}
t_{\infty}=  \mathring{t} + \frac{c_1}{\mathring{C}^2} \left(\ln(X) +\frac{1}{X} - 1 \right)\,,
\end{equation}
where
\begin{equation}
X= \frac{2c_1}{c_2\mathring{\kappa}\mathring{s}} =\frac{2\sin\alpha}{\mathring{\kappa}\mathring{s}} \,.
\end{equation}
Observe from Fig.~\ref{fig:blowupH0}  there is blow up for all values of $\mathring{C}$, but only values consistent with
$\mathring{s}<1/\mathring{\kappa}$ are acceptable. Thus,
\begin{equation}
\quad \frac{1}{\mathring{s}\mathring{\kappa}} >1\quad\Rightarrow \quad
 \frac{2\sin\alpha}{\mathring{s}\mathring{\kappa}} >2\sin\alpha\,.
\end{equation}
Therefore after choosing  $2\sin\alpha<1$, only blowup times to the right of this value in Fig.~\ref{fig:blowup_again} are consistent with \eqref{ineq}.

As for the case  $\mathring{C}=0$, given the solution of \eqref{kasol} for $\kappa$,  we immediately obtain the solutions for $s$ and $\delta$
 from \eqref{CH0} and \eqref{salde}, respectively,
\begin{equation}
s= \frac{2c_1/c_2}{\kappa - 2\mathring{C}/c_2}= \delta\, e^{-\beta_1 + 1/2} \,.
\label{sdeexact}
\end{equation}

At late times, expansion of  \eqref{kasol}  again gives
\begin{equation}
\kappa^{-1} \sim \sqrt{\frac{c_2^2}{2c_1}(t_\infty-t)}\,,
\label{lerayE2}
\end{equation}
in agreement with the  Leray scaling  suggested in \cite{MK19} and proven above in \eqref{lerayE} and \eqref{sdeleray}.   Note, in \eqref{lerayE},   the constant $\mathring{C}$ only appears via  $t_\infty$.  The zero energy case is special in that it is tractable, but the Leray scaling is ubiquitous.    Inserting \eqref{lerayE} into \eqref{sdeexact} we obtain for late times
\begin{equation}
s =  \delta\, e^{-\beta_1 + 1/2}  \sim  \sqrt{2c_1(t_\infty -t)}
\label{sleray}
\end{equation}
and so  $\kappa s\sim 2c_2/c_1$ as for the case with $\mathring{C}=0$.

 \subsection{General analysis}
 \label{ssec:genanal}

Let us now return to  \eqref{dequad} and consider the general solution.  Because  $\mathring{C}$ is arbitrary, it follows immediately from \eqref{ddede}  that there is a family of equilibrium  points of $\delta$ when the righthand side vanishes.   However, by  \eqref{kade},   all of these correspond to $\kappa=0$ (see  Appendix \ref{sec:NHS} for an analysis).  Also note, for $\kappa>0$, $\dot\delta <0$ and so even when $\mathring{H}\neq 0$,  we  expect a solution where $\delta$ approaches  zero.  If indeed $\delta\rightarrow 0$, then \eqref{sde} implies
\begin{equation}
s\approx \delta e^{ -\beta_1 +1/2} +\mathcal{O}(\delta^3)
\end{equation}
for all $\mathring{H}$; i.e., as for $\mathring{H}=0$, $\delta$ must go to zero with $s$ along the line
\begin{equation}
s= \delta e^{ -\beta_1 +1/2} \,,
\end{equation}
while from \eqref{kade} $\kappa$ must diverge as
\begin{equation}
\frac{c_2}{2}\kappa  = \mathring{C}   +  \frac{c_1}{\delta} \,  e^{\beta_1 - 1/2}\,.
  \label{kade2}
\end{equation}
Note,   the sign of the righthand side of \eqref{ddede} is definite provided $\mathring{C}>0$ and $\delta>0$.

Let us explore further the behavior for small $\delta$.  To this end we use
\[
\mathrm{erf}(\sqrt{y})
= \frac{2}{\sqrt{\pi}} \sum_{n=0}^\infty \frac{(-1)^n \, y^{n+1/2}}{n!\,(2n+1)}\,,
\]
for small $y\in\mathbb{R}_+$,   $\sqrt{y}\, \mathrm{erf}(\sqrt{y})\in C^\omega$,  and the  identities are
\begin{equation}
\frac{d}{dx}  \mathrm{erf}(x)=  \frac{2}{\sqrt{\pi}}\,e^{-x^2}
\quad\mathrm{and}\quad
\mathrm{erf}(x)=\frac{d}{dx} \left(x\, \mathrm{erf}(x) + \frac{e^{-x^2}}{\sqrt{\pi}}\right)\,.
\label{useerfid}
\end{equation}
Upon multiplying \eqref{ddede} by $\delta$ and defining $u=\delta^2/2$, it  takes the form
\begin{equation}
\dot{u}= u^{\frac12}g(u) - f(u)
\end{equation}
where $f,g\in C^\omega(\mathbb{R})$  are defined by
\bal
f(u)&= e^{-2\mathring{H}u+2\beta_1 - 1}c_1\left( \sqrt{2\pi\mathring{H} u\,} \,
  \mathrm{erf}\left(\sqrt{2\mathring{H}u\,}\right)   +    e^{-2\mathring{H}u}
\right)
\label{fu}
\\
g(u)&=- e^{-2\mathring{H}u+\beta_1 - 1/2} \sqrt{2} \, \mathring{C}\,,
\label{gu}
\eal
and  satisfy
\begin{equation}
f(0)= c_1e^{2\beta_1-1}\qquad\mathrm{and}\qquad
g(0)=-\sqrt{2} \, \mathring{C}\, e^{\beta_1 -1/2}\,.
\end{equation}
Using \eqref{useerfid}  we see that the function  $f(u)$ is positive and monotonically decreasing, while the function $g(u)$ has  a sign determined by the sign of $\mathring{C}$,  and is at least  monotonic for small argument.

The functions $f$ and $g$ are depicted in Fig.~\ref{f&g}.
\begin{figure}[htb]
\includegraphics[scale=.6]{f+g.pdf}
\caption[sm]{Plots of  $\bar{f}=  e^{-x} (\sqrt{\pi x}\, \mathrm{erf}( \sqrt{x}) +e^{-x})$ (blue), $\bar{g}= 0.10 \times e^{-x}$ (orange),  and  $f/g=\sqrt{\pi x}\, \mathrm{erf}(\sqrt{x}) +e^{-x}$ (green) vs.\  $x=2\mathring{H}u$. Note,  $\sqrt{u_\infty} \sim f(u_\infty)/g(u_\infty)$ are the equilibria points corresponding to $\kappa=0$.}
\label{f&g}
\end{figure}

Thus, near $u=0$ the MK system behaves as
\begin{equation}
\dot{u}= u^{\frac12}g(0) - f(0)
\label{dotsqu}
\end{equation}
Because $u^{\frac12}g(0) - f(0)$ is not Lipschitz on $u\in[0,\epsilon]$ with $\epsilon>0$, we do not have the usual ODE existence theorem to rely on (see e.g.\  \cite{c&l}).  If $u(0)=\mathring{u}>0$ is small, then $\dot{u}(0)= \mathring{u}^{\frac12}g(0) - f(0)<0$, because $c_1>0$, and  $u$ should decrease.
However, the Lipschitz condition for uniqueness is  a sufficient but not a necessary condition.  Thus, further analysis is necessary, but indeed
\eqref{dotsqu} does have a unique solution. To see this let  $u_1$ and $u_2$ be  two solutions that coincide at some time.   A measure of their  difference
\begin{equation}
D = \left(\sqrt{u_1(t)} -\sqrt{u_2(t)}\right)^2
\end{equation}
satisfies $\dot D=0$ for all time,  so they must coincide.

\bigskip

From \eqref{dotsqu} we see there is a family of equilibrium points given by
\begin{equation}
 \sqrt{u_0} = \frac{f(0)}{g(0)}=  \frac{c_1e^{\beta_1-1/2}}{ \sqrt{2}\, \mathring{C} }
 \end{equation}
which only exists for $\mathring{C}>0$.  However, these are merely the $\kappa=0$ solutions mentioned above.
In Appendix \ref{sec:NHS} it is seen that linearization about any of these equilibrium solutions  yields a spectrum with   two zero eigenvalues and one unstable (positive) eigenvalue.   This implies $\kappa$ will grow so as to decrease the radius of curvature.

In any event,  it is a simple matter to integrate \eqref{dotsqu}, yielding
\begin{equation}
 \frac{2}{g(0)}\left[
 u^{\frac12} + \frac{f(0)}{g(0)} \ln\left(1-  {u^{\frac12} g(0)}/{f(0)}\right)
 \right]= t + \mathrm{const}\,,
\label{Cgeq0}
\end{equation}
 which, upon expanding,  gives for small $u$
\begin{equation}
\delta\sim   \sqrt{2f(0)({t_\infty-t)}}=e^{\beta_1-1/2}\sqrt{2 c_1(t_\infty-t)}\,.
\end{equation}
the  solution with Leray scaling consistent with \eqref{sleray} and \eqref{lerayE}.

\section{Conclusions}
\label{sec:conclu}

We have shown that the system put forth  by  Moffat and Kimura \cite{MK19,MK19b} for the interaction of two tilted vortex rings has solutions with  finite-time singularity.  This was  achieved by finding the noncanonical Hamiltonian structure of the equations, which naturally led to a geometrical depiction and explicit forms for the solutions by making use of the newly discovered Hamiltonian and Casimir invariants.   Exact Leray divergence was demonstrated, within the inequalities proposed in \cite{MK19,MK19b} for the  model.

Future publications will consider the effect of viscous dissipation, further physical interpretation of the results,  and various bounds and perturbation expansions.

\section*{Acknowledgment}
\noindent   PJM was supported by U.S. Dept.\ of Energy Contract \# DE-FG05-80ET-53088 and a Forschungspreis from the Alexander von Humboldt Foundation.  YK acknowledges support from JSPS KAKENHI grant \# JP19122056 and \# JP16747321.


text:
\section{Spectra and wave functions of Edge States}

In $K$ and $K'$ valleys of sTMD monolayers electrons and hole around band edges are described by the following Hamiltonian \cite{bib:Xiao}:
\small
\begin{equation}\label{SM_H_0}
H_{\tau} = v\left(\tau\sigma_xp_x + \sigma_yp_y\right) +\sigma_z\left[m - \tau \hat{s}_z(\Delta_v-\Delta_c)\right]+ \sigma_0 \tau \hat{s}_z(\Delta_v+\Delta_c),
\end{equation}
\normalsize
where $\tau=\pm 1$ is the valley index, $\sigma_{x,y,z}$ and $\sigma_0$ are Pauli and identity matrixes, respectively, acting in band subspace, $\hat{s}_z$ is $z$ component of the spin-1/2 operator, $v$ is matrix element of momentum between the band extremum wave functions, $2\Delta_{c,v}$ is atomic spin-orbit splitting in conduction and valence bands, respectively, $2m$ is band gap without the splitting.

Below in this section we calculate spectra and wave functions of edge states (ESs) for two different types of spin-orbit interaction (SOI) mechanisms.

\subsection{Bulk SOI}\label{section_BSOI}

First, we consider bulk mechanisms of SOI, such as Bychkov-Rashba SOI with substrate, SOI induced by perpendicular electric field, and SOI induced by in-plane magnetic field. In general, these interactions can be described by a perturbation Hamiltonian:
\begin{equation}\label{SM_H_1}
\delta H_{\tau} =
\left(
\begin{array}{cc}
0 & V_{\tau} \\
V^{+}_{\tau} & 0
\end{array}
\right).
\end{equation}
Here $V_{\tau}$ is square matrix of the second order, which in case of Bychkov-Rashba interface SOI or SOI caused by perpendicular electric field, has the form \cite{Roldan,Basko}:
\begin{equation}
V_{\tau} =
\left(
\begin{array}{cc}
0 & \lambda^*\, \delta_{\tau,-1} \\
-\lambda\, \delta_{\tau,1} & 0
\end{array}
\right),
\end{equation}
where $\lambda$ is parameter that describes value of the interactions ($\lambda\sim 1$ meV). In case of in-plane magnetic field $V_{\tau}=[g\mu_B/2](B_x-iB_y)\sigma_0$ (\mbox{$g\approx 2$} is effective $g$-factor and $\mu_B$ is the Bohr magneton). For bulk SOI we employ spin-decoupled boundary condition (BC) at the edge:
\begin{equation}\label{SM_BC_1} \left[\psi_{c,\tau}^{\uparrow (\downarrow)} + a_{1,2;\tau} \psi_{v,\tau}^{\uparrow (\downarrow)}\right]_{y=0}=0,
\end{equation}
where $a_{1,\tau}$ ($a_{2,\tau}$) is real phenomenological parameter that characterizes edge properties for spin up (spin down) states. Time reversal symmetry impose the following relation: $a_{1,\tau}=a_{2,-\tau}$ \cite{Enaldiev}. In zero approximation on $\delta H$ (\ref{SM_H_1}), Hamiltonian (\ref{SM_H_0}) with BC (\ref{SM_BC_1}) gives rise to two subbands of ESs in each valley \cite{Enaldiev}:
\begin{equation}\label{SM_ES_spectra}
\varepsilon_{n,p}^{(\tau)} = -\tau v_{n}^{(\tau)}p + \Delta_{n}^{(\tau)},
\end{equation}
here \mbox{$n=1,2$} is the ES subband index, \mbox{$v_{n}^{(\tau)} = 2a_{n,\tau}v/(1 + a^2_{n,\tau})$} is ES velocity and
\begin{equation}
\Delta_{n}^{(\tau)} = m\frac{a^2_{n,\tau}-1}{a^2_{n,\tau} + 1} + (-1)^{n+1}\tau\frac{\Delta_v + a^2_{n,\tau} \Delta_c}{1 + a^2_{n,\tau}}
\end{equation}
is value of e-h asymmetry. ES subbands possess linear dispersion and are splitted in energy as $\Delta^{(\tau)}_{1}\neq\Delta^{(\tau)}_{2}$. Zero approximation wave functions of the ESs have the following form:
\begin{eqnarray}\label{SM_ES_WF}
\psi^{(0)}_{1,p,\tau} = C_{1,p,\tau}
\left( 1, -1/a_{1,\tau}, 0, 0 \right)^{\rm T} e^{-\kappa_{1,p,\tau}y + ipx}, \\ \psi^{(0)}_{2,p,\tau} = C_{2,p,\tau}
\left( 0, 0, 1, -1/a_{2,\tau}\right)^{\rm T} e^{-\kappa_{2,p,\tau}y + ipx},
\end{eqnarray}
where $C_{n,p,\tau} = \left(2\kappa_{n,p,\tau}a^2_{n,\tau}/(1 + a^2_{n,\tau} )L_x\right)^{1/2}$ is normalization factor, and
\begin{equation}\label{SM_decay_length}
\hbar\kappa_{n,p,\tau} = -\tau p \dfrac{1 - a_{n,\tau}^2}{1 + a_{n,\tau}^2} +  \frac{ 2a_{n,\tau} }{1 + a^2_{n,\tau} } \left(\frac{m}{v} -(-1)^{n+1}\tau \frac{\Delta_v - \Delta_c}{2v} \right)
\end{equation}
is inverse decay length of ESs. ESs exist for those quasimomenta while their inverse decay length (\ref{SM_decay_length}) is positive. In Ref.\cite{Enaldiev} was shown that ESs (\ref{SM_ES_WF}),(\ref{SM_ES_spectra}) are in the gap at $a_{1,2;1}>0$ ($a_{1,2;-1}<0$). We also assume that $a_{\pm 1,\tau}\sim 1$ for ESs to be in the gap at small enough quasimomenta. Because of quite large band gap, the inverse decay length (\ref{SM_decay_length}) weakly depends on the quasimomentum of the ESs. Therefore, below we take its value at Fermi-level at calculation of plasmon spectra.

Perturbation (\ref{SM_H_1}) couples ESs belonging to different subbands. Since ES subbands are non degenerate, only their wave functions are changed in the first approximation on $\delta H_{\tau}$ (\ref{SM_H_1}):
\begin{equation}\label{SM_WF_bulk_perturb_1}
\psi_{1,p,\tau} = \psi^{(0)}_{1,p,\tau} + \frac{\langle\psi_{2,p,\tau}^{(0)}|\delta H_{\tau}|\psi_{1,p,\tau}^{(0)}\rangle}{\varepsilon_{1,p}^{(\tau)} - \varepsilon_{2,p}^{(\tau)}}\psi^{(0)}_{2,p,\tau},
\end{equation}

\begin{equation}\label{SM_WF_bulk_perturb_2}
\psi_{2,p,\tau} = \psi^{(0)}_{2,p,\tau} + \frac{\langle\psi_{1,p,\tau}^{(0)}|\delta H_{\tau}|\psi_{2,p,\tau}^{(0)}\rangle}{\varepsilon_{2,p}^{(\tau)} - \varepsilon_{1,p}^{(\tau)}}\psi^{(0)}_{1,p,\tau}.
\end{equation}

\subsection{Edge SOI}\label{section_ESOI}

In this section we consider the second type of SOI, that is caused by breaking of crystal potential in perpendicular to edge direction. Such SOI can be described by a spin-coupled BC \cite{Enaldiev}:
\begin{equation}\label{SM_BC_2}
\left[\Psi^{\uparrow}_{\tau} - M_{\tau}\Psi^{\downarrow}_{\tau} \right]_{y=0} = 0,
\end{equation}
where $\Psi^{\uparrow(\downarrow)}_{\tau} = (\psi_{c,\tau}^{\uparrow (\downarrow)},\psi_{v,\tau}^{\uparrow (\downarrow)})$ is two-component envelope functions, and $M$ is the second order matrix with elements:
\begin{align}
m_{11} = i(\sinh\xi + \tau\cosh\xi\cosh\eta\cos\nu)\\
m_{22} = i(\sinh\xi - \tau\cosh\xi\cosh\eta\cos\nu)\\
m_{12} = i\tau\cosh\xi(\sinh\eta - \cosh\eta\sin\nu) \\
m_{21} = -i\tau\cosh\xi(\sinh\eta + \cosh\eta\sin\nu),
\end{align}
where $\xi,\eta,\nu$ are phenomenological parameters that describe SOI at the edge. In case of the edge SOI, there also exists two subbands of ESs which are described by the following wave functions:
\begin{equation}\label{SM_WF_ESOI}
\psi_{n,p,\tau} = C_{n,p,\tau}
\left[
\begin{array}{c}
M_{\tau}
\left(
\begin{array}{c}
1 \\
G^{\downarrow}_{n,p,\tau}
\end{array}
\right)e^{\kappa^{\uparrow}_{n,p,\tau}y} \\
\left(
\begin{array}{c}
1 \\
G^{\downarrow}_{n,p,\tau}
\end{array}
\right)e^{\kappa^{\downarrow}_{n,p,\tau}y}
\end{array}
\right]e^{ipx},
\end{equation}
where
\begin{equation}
G^{\uparrow(\downarrow)}_{n,p,\tau} = \frac{v(\tau p - \hbar\kappa^{\uparrow(\downarrow)}_{n,p,\tau})}{\varepsilon_{n,p}^{(\tau)}+m\mp\tau\Delta_v}\nonumber
\end{equation}
and inverse decay lengthes are expressed as follows:
\begin{equation}\label{SM_decay_length_ESOI}
\hbar\kappa^{\uparrow(\downarrow)}_{n,p,\tau} = \sqrt{p^2-(\varepsilon_{n,p}^{(\tau)} - m \mp\tau\Delta_c)(\varepsilon_{n,p}^{(\tau)} + m \mp\tau\Delta_v)/v^2}.
\end{equation}
Energies of ES subbands are determined by a dispersion equation:
\begin{align}\label{SM_disp_eq_ESOI}
\left(1 + a_{1,\tau}G^{\uparrow}_{n,p,\tau}\right)\left(1 + a_{2,\tau}G^{\downarrow}_{n,p,\tau}\right) + \nonumber\\
\frac{\tau v(\tanh\xi - 1)}{\sinh\eta + \cosh\eta\sin\nu}\left(G^{\uparrow}_{n,p,\tau} - G^{\downarrow}_{n,p,\tau}\right)=0,
\end{align}
where
\begin{equation}\label{SM_a_12}
a_{1,2;\tau}=\pm\tau \frac{1 \pm\tau \cosh\eta\cos\nu}{\sinh\eta+\cosh\eta\sin\nu}
\end{equation}
are expressions for the introduced previously section boundary parameters via $\eta,\nu$. Coupling of the ES subbands by the edge SOI is provided by the last term in Eq.(\ref{SM_disp_eq_ESOI}). Taking this term as a perturbation, we obtain the same spectra of ESs as in case of the bulk SOI (\ref{SM_ES_spectra}), but with renormalized e-h asymmetry term:
\begin{align}
\widetilde{\Delta}_{n}^{(\tau)} = \Delta_{n}^{(\tau)} +
\frac{(-1)^{n+1}2\tau\Delta_v(\tanh\xi - 1)}{(1+a_{n,\tau}^2)\left(1 + \frac{\Delta_v (\tau(-1)^{n+1}+\cosh\eta\sin\nu)}{\varepsilon_{n,0}^{(\tau)} + m - (-1)^{n+1}\tau\Delta_v}\right)}
\end{align}
As it follows from Eq.(\ref{SM_a_12}) \mbox{$a_{1,1}>a_{2,1}$}. This means that ESs from the upper subband always has greater absolute value of velocity than from the lower one. Therefore, ES subbands do not intersect each other at any quasimomentum as without SOI as well taking the latter into account.

\section{Derivation of plasmon dispersion equation in self consistent Random phase approximation}

In this section we derive general dispersion equation for 1D plasmons in self-consistent random phase approximation (RPA), taking into account only transitions between ESs. First, we solve quantum kinetic equation for density matrix $\rho$ in the first order on electric potential:
\begin{equation}\label{SM_QKE}
i\hbar \frac{\partial \rho}{\partial t} = \left[H_0  - e\left(\varphi_q(y)e^{iqx-i(\omega+i0)t}+ c.c.\right), \rho\right]
\end{equation}
where $H_0=H_{\tau} + \delta H_{\tau}$ with BC (\ref{SM_BC_1}) for ES wave function set (\ref{SM_WF_bulk_perturb_1}), (\ref{SM_WF_bulk_perturb_2}), or $H_0=H_{\tau}$ with BC (\ref{SM_BC_2}) for ES wave function set (\ref{SM_WF_ESOI}). The first order amendment to density matrix reads as follows:
\begin{equation}\label{SM_DM_1}
\rho^{(1)}_{\lambda,\lambda'} = -e\langle\lambda|\varphi_q(y)e^{iqx}|\lambda'\rangle\delta_{\tau_{\lambda},\tau_{\lambda'}} \frac{f(\varepsilon_{\lambda'}) - f(\varepsilon_{\lambda})}{\varepsilon_{\lambda'} - \varepsilon_{\lambda} + \hbar\omega+ i0},
\end{equation}
where $e$ is absolute value of electron charge, $\lambda$ is set of ES quantum numbers, $f(\varepsilon)$ is Fermi-Dirac distribution function. Fourier harmonics of induced electron density is expressed in the following way:
\begin{equation}\label{SM_density}
n^{(1)}_{q}(y_e) = \frac{-e}{L_x}\sum_{\lambda,\lambda'}\rho^{(1)}_{\lambda,\lambda'}\langle \lambda'|\delta(y-y_e)e^{-iqx}|\lambda\rangle
\end{equation}

Inserting (\ref{SM_density}) in the Poisson equation we obtain a closed equation for the electric potential harmonics:
\begin{equation}\label{SM_Poisson}
\varphi_q(y)e^{iqx} = \frac{2}{\epsilon^*}\int_{0}^{+\infty}K_0\left(\left|q(y-y')\right|\right)n_q^{(1)}(y')e^{iqx}dy',
\end{equation}
where $K_0(x)$ is the Macdonald function describing Coulomb interaction in one dimension, $\epsilon^*$ is dielectric constant of environment. Substituting (\ref{SM_density}) in the previous equation and taking matrix element $\langle \mu'|\dots|\mu\rangle \delta_{\tau_{\mu}\tau_{\mu'}}$ from both side of the equation, we arrive to the final expression used in the main text:
\begin{eqnarray}\label{SM_RPA_gen}
\langle \mu'|\varphi_{q}(y)e^{iqx}|\mu\rangle\delta_{\tau_{\mu'}\tau_{\mu}} = \qquad\qquad\qquad\qquad\qquad\qquad\qquad\nonumber\\
 \frac{2e^2}{\epsilon^*L_x}\sum_{\lambda\lambda'}\frac{f(\varepsilon_{\lambda'}) - f(\varepsilon_{\lambda}) }{ \varepsilon_{\lambda'} - \varepsilon_{\lambda} + \hbar\omega +i0 }\langle \lambda|\varphi_{q}(y)e^{iqx}|\lambda'\rangle\delta_{\tau_{\lambda'}\tau_{\lambda}}\times \nonumber\\
\langle\lambda'|\langle\mu'|e^{iq(x-x')}K_0\left(|q(y-y')|\right)|\mu\rangle|\lambda\rangle \qquad\quad
\end{eqnarray}

\subsection{Dispersion equation for intrasubband plasmons}

Leaving only diagonal matrix elements on the subband number in Eq.(\ref{SM_RPA_gen}) (f.e. \mbox{$|\mu\rangle=|1,p',1\rangle $} and \mbox{$|\mu'\rangle =|1, p'+q,1\rangle$}) we arrive to the following equation:
\small
\begin{widetext}
\begin{equation}\label{SM_RPA_intra}
\begin{split}
\langle 1,p'+q,1\left|\varphi_{q}(y)e^{iqx}\right|1,p',1\rangle = \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
\\
\frac{2e^2}{\epsilon^*L_x}\sum_{p,\tau}
\frac{ f(\varepsilon_{1,p}^{(\tau)}) - f(\varepsilon_{1,p+q}^{(\tau)}) }{ \varepsilon_{1,p}^{(\tau)} - \varepsilon_{1,p+q}^{(\tau)} +\hbar\omega + i0} \langle 1,p+q,\tau\left|\varphi_{q}(y)e^{iqx}\right|1,p,\tau\rangle
\langle 1,p'+q,1|\langle 1,p,\tau\left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|1,p+q,\tau\rangle|1,p',1\rangle +  \\
\frac{2e^2}{\epsilon^*L_x}\sum_{p,\tau}
\frac{ f(\varepsilon_{2,p}^{(\tau)}) - f(\varepsilon_{2,p+q}^{(\tau)}) }{ \varepsilon_{2,p}^{(\tau)} - \varepsilon_{2,p+q}^{(\tau)} +\hbar\omega + i0 } \langle 2,p+q,\tau\left|\varphi_{q}(y)e^{iqx}\right|2,p,\tau\rangle
\langle 1,p'+q, 1|\langle 2,p,\tau\left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|2,p+q,\tau\rangle|1,p',1\rangle.
\end{split}
\end{equation}
\end{widetext}
\normalsize
As ESs are exponentially localized in vicinity of the edge, the diagonal matrix elements of electric potential in the long wavelength limit $q\kappa^{-1}\ll 1$ ($\kappa$ is inverse decay length of ESs at Fermi-level) can be estimated as:
\begin{eqnarray}\label{SM_potential_intra}
\langle 1,p+q,\tau|\varphi_{q}(y)e^{iqx}|1,p,\tau\rangle = \nonumber \qquad\qquad\qquad\\
\langle 2,p+q,\tau|\varphi_{q}(y)e^{iqx}|2,p,\tau\rangle \approx \varphi_{q}(0).
\end{eqnarray}
Taking asymptotic expression for 1D Coulomb potential $K_0(x)\approx \ln(2/x\zeta)$ in the long wavelength limit ($\zeta=e^{\gamma}$, $\gamma=0.577\dots$ is the Euler constant), we evaluate its matrix elements as follows:
\small
\begin{equation}\label{SM_Coulomb_intra_1}
\begin{split}
\langle 1,p'+q,1|\langle 1,p,\tau\left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|1,p+q,\tau\rangle|1,p',\tau\rangle \approx \qquad\\
 \ln\left[\frac{2\left(\kappa_{1,p,\tau} +\kappa_{1,p',1}\right) }{q\zeta}\right]\approx \ln\left[\frac{4\kappa_{1,p_{F_1},1}}{q\zeta}\right], \qquad\qquad\qquad\qquad\qquad
\end{split}
\end{equation}
\begin{equation}\label{SM_Coulomb_intra_2}
\begin{split}
\langle 1,p'+q,1|\langle 2,p,\tau\left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|2,{\small p+q},\tau\rangle|1,p',1\rangle \approx \qquad\qquad\qquad\\
\ln\left[\frac{2\left(\kappa_{2,p,\tau} +\kappa_{1,p', 1}\right) }{q\zeta}\right]\approx \ln\left[\frac{2\left( \kappa_{2,p_{F_2},1} + \kappa_{1,p_{F_1},1}\right)}{q\zeta}\right] \approx \qquad\qquad\qquad\qquad\quad \\ \ln\left[\frac{4\kappa_{2,p_{F_2},1}}{q\zeta}\right], \qquad\qquad\qquad\qquad\quad
\end{split}
\end{equation}
\normalsize
Since the inverse decay lengths (\ref{SM_decay_length}) or (\ref{SM_decay_length_ESOI}) is of weak dependence on quasimomentum, we take its value at corresponding Fermi-level. Further, due to smoothness of logarithm function, we will put $\ln(4\kappa_{1,p_{F_1},1}/q\zeta)\approx\ln(4\kappa_{2,p_{F_2},1}/q\zeta)$ in the matrix elements (\ref{SM_Coulomb_intra_1}),(\ref{SM_Coulomb_intra_2}).  Inserting Eqs.(\ref{SM_potential_intra})-(\ref{SM_Coulomb_intra_2}) in the system (\ref{SM_RPA_intra}) we finally obtain the dispersion equation for intrasubband plasmons that is used in the main text:
\begin{equation}
\begin{split}
1-\frac{2e^2}{\epsilon^*}\ln\left(\frac{4\kappa_{2,p_{F_2},1}}{|q|\zeta}\right)\sum_{\tau}\left[\Pi^{(\tau)}_{11}(\omega, q) + \Pi^{(\tau)}_{22}(\omega, q)\right]=0,\\
\Pi^{(\tau)}_{n'n}(\omega, q)=\frac{1}{L_x}\sum_{p}\frac{f\left(\varepsilon_{n',p}^{(\tau)}\right) - f\left(\varepsilon_{n,p+q}^{(\tau)}\right) }{ \varepsilon_{n',p}^{(\tau)} - \varepsilon_{n,p+q}^{(\tau)} + \hbar\omega +i0 }, \, n,n'=\{1,2\} \\
\\
\\
\\
\end{split}
\end{equation}
Therefore, in the main text $\kappa_{\rm intra}=\kappa_{2,p_{F_2},1}\approx \kappa_{1, p_{F_1},1}$.

\subsection{Dispersion equation for intersubband plasmons}

Taking \mbox{$|\mu'\rangle=|1,p'+q,1\rangle$}, \mbox{$|\mu\rangle=|2,p',1\rangle$} and keeping only intersubband matrix elements in Eq.(\ref{SM_RPA_gen}) we obtain the following equation:
\small
\begin{widetext}
\begin{equation}\label{SM_RPA_inter}
\begin{split}
\langle 1,p'+q, 1\left|\varphi_{q}(y)e^{iqx}\right|2, p', 1\rangle = \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad  \\
\frac{2e^2}{\epsilon^*L_x}\sum_{p,\tau}
\frac{ f(\varepsilon_{1,p}^{(\tau)}) - f(\varepsilon_{2,p+q}^{(\tau)}) }{ \varepsilon_{1,p}^{(\tau)} - \varepsilon_{2,p+q}^{(\tau)} +\hbar\omega + i0 } \langle 2,p+q,\tau\left|\varphi_{q}(y)e^{iqx}\right|1,p,\tau\rangle
\langle 1,p'+q, 1|\langle 1,p,\tau\left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|2,p+q,\tau\rangle|2,p', 1\rangle +  \\
\frac{2e^2}{\epsilon^*L_x}\sum_{p,\tau}
\frac{ f(\varepsilon_{2,p}^{(\tau)}) - f(\varepsilon_{1,p+q}^{(\tau)}) }{ \varepsilon_{2,p}^{(\tau)} - \varepsilon_{1,p+q}^{(\tau)} +\hbar\omega + i0 } \langle 1,p+q,\tau\left|\varphi_{q}(y)e^{iqx}\right|2,p,\tau\rangle
\langle 1,p'+q,1|\langle 2,p,\tau\left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|1,p+q,\tau\rangle|2,p',1\rangle.
\end{split}
\end{equation}
\end{widetext}
\normalsize
In long wavelength limit intersubband matrix elements of electric potential are evaluated as follows:
\small
\begin{equation}\label{SM_potential_inter_1}
\begin{split}
\langle 1,p+q,\tau\left|\varphi_{q}(y)e^{iqx}\right|2,p,\tau\rangle \approx \varphi_{q}(0)\langle 1,p+q,\tau\left|e^{iqx}\right|2,p,\tau\rangle \equiv \\
 \varphi_{q}(0)J_{12}^{(\tau)}(p,q),
\end{split}
\end{equation}
\begin{equation}\label{SM_potential_inter_2}
\begin{split}
\langle 2,p+q,\tau\left|\varphi_{q}(y)e^{iqx}\right|1,p,\tau\rangle\approx \varphi_{q}(0)\langle 2,p+q,\tau\left|e^{iqx}\right|1,p,\tau\rangle \equiv \\
 \varphi_{q}(0)J_{21}^{(\tau)}(p,q),
\end{split}
\end{equation}
\normalsize
where the last equality is in (\ref{SM_potential_inter_1}), (\ref{SM_potential_inter_2}) is a definition of $J_{12,21}^{(\tau)}(p,q)$. In the long wavelength limit intersubband matrix elements of 1D Coulomb potential read as:
\small
\begin{equation}\label{SM_Coulomb_inter}
\begin{split}
\langle 1,p'+q, 1|\langle 1,p,\tau\left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|2,p+q,\tau\rangle|2,p', 1\rangle \approx\\
\ln\left[\frac{2(\kappa_{1,p_{F_2},1}+\kappa_{2,p_{F_2},1})}{q\zeta}\right]J_{12}^{(1)}(p',q)J_{21}^{(\tau)*}(p,q) \qquad\qquad\qquad \\
\langle 1,p'+q, 1|\langle 2, p, \tau \left|K_0(\left|q(y-y')\right|)e^{iq(x-x')}\right|1,p+q,\tau\rangle|2, p', 1\rangle \approx\\
\ln\left[\frac{2(\kappa_{1,p_{F_2}, 1}+\kappa_{2,p_{F_2},1})}{q\zeta}\right]J_{12}^{(1)}(p',q)J_{12}^{(\tau)*}(p,q) \qquad\qquad\qquad \\
\end{split}
\end{equation}
\normalsize
where we took the inverse decay lengths at $p=p_{F_2}$. In the long wavelength limit matrix elements $J_{12,21}^{(\tau)}(p,q)$ are proportional to wave vector $q$, as at \mbox{$q=0$} ES wave functions belonging different subbands are orthogonal to each other. In the case of edge SOI the matrix elements satisfy symmetry condition: \mbox{$J_{12}^{(1)}(p,q,1)=J_{21}^{(-1)}(-p,-q,-1)$}. But, for the bulk SOI (\ref{SM_H_1}) the condition is following one: \mbox{$J_{12}^{(1)}(p,q,1)=-J_{21}^{(-1)}(-p,-q,-1)$}. Since at small wave vectors the absolute values of the matrix elements approximately equal to each other (i.e. $|J_{12}^{1}(p,q)\approx|J_{21}^{1}(p,q)|\approx|J_{12}^{-1}(-p,q)|\approx|J_{21}^{(-1)}(p,q)|$) and additionally possess weak dependence on quasimomentum (see Figs.\ref{SM_J_BSOI},\ref{SM_J_ESOI}) we treat them as constants taken at $p=p_{F_2}$ in the Eq.(\ref{SM_RPA_inter}). Finally , we obtain the following dispersion equation for intersubband plasmons:
\begin{equation}\label{SM_inter_disp}
\begin{split}
1-\frac{2e^2}{\epsilon^*L_x}\left|J_{12}(p_{F_2},q)\right|^2\ln\left ( \frac{2\left(\kappa_{1,p_{F_2}, 1}+\kappa_{2,p_{F_2},1}\right)}{|q|\zeta} \right )\times \\
\sum_{\tau}\left[\Pi^{(\tau)}_{12}(\omega,q) + \Pi^{(\tau)}_{21}(\omega,q)\right] =0,
\end{split}
\end{equation}
Therefore, in the main text $\kappa_{\rm inter}=(\kappa_{2,p_{F_2},1} + \kappa_{1,p_{F_2},1})/2$.

\begin{figure}
\includegraphics{SM_Fig2.eps}
\caption{Quasimomentum dependence of intersuband matrix element absolute values $|J_{12}^{(\tau)}(\tau p,q)|$ (blue), $|J_{21}^{(\tau)}(\tau p,q)|$ (red) in case of bulk SOI at three values of wave vectors $q=0.05p_{F_2}0.1p_{F_2},0.2p_{F_2}$. Solid (dashed) lines respond to the $K$ ($K'$) valley. Vertical line shows quasimomenta corresponding to Fermi-level in each edge subband (see, Fig.\ref{Fig1} of the main text). Calculation was carried out at: $2m=1.8$ eV, $v=2.5$ eV$\cdot$A, $2\Delta_c=3$ meV, $2\Delta_v=0.148$ eV, $a_{1,1}=a_{2,-1}=0.5$, $a_{2,1}=a_{1,-1}=0.35$, $|\lambda|=1$meV.
\label{SM_J_BSOI}}
\end{figure}

\begin{figure}
\includegraphics{SM_Fig1.eps}
\caption{Quasimomentum dependence of intersuband matrix element absolute values $|J_{12}^{(\tau)}(\tau p,q)|$ (blue), $|J_{21}^{(\tau)}(\tau p,q)|$ (red) in case of edge SOI at three values of wave vectors $q=0.05p_{F_2}0.1p_{F_2},0.2p_{F_2}$. Solid (dashed) lines respond to the $K$ ($K'$) valley. Vertical line shows quasimomenta corresponding to Fermi-level in each edge subband (see, Fig.\ref{Fig1} of the main text). Calculation was carried out at: $2m=1.8$ eV, $v=2.5$ eV$\cdot$A, $2\Delta_c=3$ meV, $2\Delta_v=0.148$ eV, $\xi=0.2$, $\eta=2.76$, $\nu=0.77$. \label{SM_J_ESOI}}
\end{figure}

\end{document}

redpajama_set_name: RedPajamaArXiv