2 The family of Norbury’s vortex rings

In this section we recall the family of inviscid vortex rings discovered by Norbury (Reference Norbury1973). Since these flows are axisymmetric, we will adopt the cylindrical–polar coordinate system $(r,\unicode[STIX]{x1D719},z)$ and define the operator ${\mathcal{L}}:=\unicode[STIX]{x1D735}\boldsymbol{\cdot }((1/r)\,\unicode[STIX]{x1D735})$ in which $\unicode[STIX]{x1D735}:=[\unicode[STIX]{x2202}/\unicode[STIX]{x2202}r,\unicode[STIX]{x2202}/\unicode[STIX]{x2202}z]^{\text{T}}$ (‘ $:=$ ’ means ‘equal to by definition’). Given the Stokes streamfunction $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}(r,z)$ , these flows satisfy the following system in the frame of reference moving in the direction of the $z$ axis with the translation velocity $W$ of the vortex

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{L}}\unicode[STIX]{x1D713}=-rf(\unicode[STIX]{x1D713}),\quad \text{in}~\unicode[STIX]{x1D6FA}, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}\rightarrow -{\textstyle \frac{1}{2}}Wr^{2},\quad \text{as}~|\boldsymbol{x}|:=\sqrt{r^{2}+z^{2}}\rightarrow \infty , & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FA}=\mathbb{R}^{2}$

$\boldsymbol{x}=[r,z]^{\text{T}}$

$f(\unicode[STIX]{x1D713})$

(2.2) $$\begin{eqnarray}f(\unicode[STIX]{x1D713})=\left\{\begin{array}{@{}ll@{}}{\mathcal{C}},\quad & \unicode[STIX]{x1D713}>\unicode[STIX]{x1D713}_{0}\\ 0,\quad & \unicode[STIX]{x1D713}\leqslant \unicode[STIX]{x1D713}_{0}\end{array}\right.,\end{eqnarray}$$

in which ${\mathcal{C}}\neq 0$ and $\unicode[STIX]{x1D713}_{0}\geqslant 0$ are constants. System (2.1)–(2.2) therefore describes a compact region ${\mathcal{D}}$ with azimuthal vorticity $\unicode[STIX]{x1D714}:=\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}$ varying proportionally to the distance $r$ from the flow axis embedded in a potential flow. The boundary of this region $\unicode[STIX]{x2202}{\mathcal{D}}:=\{(r,\unicode[STIX]{x1D719},z):\unicode[STIX]{x1D713}(r,z)=\unicode[STIX]{x1D713}_{0}\}$ is a priori unknown and must be found as a part of the problem solution. System (2.1)–(2.2) thus represents a free-boundary problem and, as will become evident below, this property makes the study of the stability of its solutions more complicated.

Figure 1. Schematic representation of the problem geometry in the meridional plane $(r,z)$ with a representative member of Norbury’s family of vortex rings (marked with a think blue line).

Solution of the stability problem requires a convenient representation of the boundary $\unicode[STIX]{x2202}{\mathcal{D}}$ of the vortex ring. The intersection of the vortex region ${\mathcal{D}}$ with the meridional plane $(r,z)$ , which we will denote ${\mathcal{A}}$ , is shown schematically in figure 1. We consider the following general parametrization of the vortex boundary $\boldsymbol{x}(\unicode[STIX]{x1D703})=[r(\unicode[STIX]{x1D703}),z(\unicode[STIX]{x1D703})]^{\text{T}}$ , where $\unicode[STIX]{x1D703}\in [0,2\unicode[STIX]{x03C0}]$ and $r:[0,2\unicode[STIX]{x03C0}]\rightarrow \mathbb{R}^{+}$ , $z:[0,2\unicode[STIX]{x03C0}]\rightarrow \mathbb{R}$ are smooth periodic functions, which can be viewed as a mapping from a unit circle to $\unicode[STIX]{x2202}{\mathcal{A}}$ . In his original work, Norbury (Reference Norbury1973) employed the following particular form of this parametrization adapted to solutions with the fore-and-aft symmetry and given in terms of some even periodic function $g:[0,2\unicode[STIX]{x03C0}]\rightarrow \mathbb{R}$

(2.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle r(\unicode[STIX]{x1D703})=g(\unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D703}+1, & \displaystyle\end{eqnarray}$$

(2.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle z(\unicode[STIX]{x1D703})=g(\unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$

$r=r(\unicode[STIX]{x1D703})$

$z=z(\unicode[STIX]{x1D703})$

$\unicode[STIX]{x2202}{\mathcal{A}}$

$P$

${\mathcal{A}}$

$|OP|=:L=(1/2)[r(0)+r(\unicode[STIX]{x03C0})]$

$L=1$

$|{\mathcal{A}}|=\unicode[STIX]{x03C0}L^{2}\unicode[STIX]{x1D6FC}^{2}=\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}^{2}$

$\unicode[STIX]{x1D6FC}$

$U$

${\mathcal{C}}=U/(L^{2}\unicode[STIX]{x1D6FC}^{2})$

${\mathcal{C}}=1/\unicode[STIX]{x1D6FC}^{2}$

$\unicode[STIX]{x1D6FC}$

$W$

$\unicode[STIX]{x1D713}_{0}$

$\unicode[STIX]{x1D6FC}\rightarrow 0$

$\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$

${\mathcal{A}}$

A consistent way to study the evolution and equilibria of inviscid flows with discontinuous vorticity distributions is based on the formalism of contour dynamics (Pullin Reference Pullin1992). Given a time-dependent region ${\mathcal{A}}(t)$ , where $t$ is time, its evolution can be studied by tracking the points $\boldsymbol{y}(t)$ on its boundary $\unicode[STIX]{x2202}{\mathcal{A}}(t)$ via the equation

(2.4) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{y}(t)}{\text{d}t}=\boldsymbol{v}(\boldsymbol{y}(t)):={\mathcal{C}}\,\int _{\unicode[STIX]{x2202}{\mathcal{A}}(t)}\boldsymbol{K}(\boldsymbol{y}(t),\boldsymbol{y}^{\prime })\,\text{d}s_{\boldsymbol{ y}^{\prime }},\quad \forall \boldsymbol{y}(t)\in \unicode[STIX]{x2202}{\mathcal{A}}(t),\end{eqnarray}$$

where $\boldsymbol{v}$ denotes the velocity, $\boldsymbol{K}(\boldsymbol{y}(t),\boldsymbol{y}^{\prime })$ is a suitable Biot–Savart kernel, $\boldsymbol{y}$ and $\boldsymbol{y}^{\prime }$ are defined in the absolute frame of reference and $\text{d}s_{\boldsymbol{y}^{\prime }}$ is an arclength element of the vortex boundary in the meridional plane. An equilibrium shape of the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ when the right-hand side in (2.1a ) has a step discontinuity, cf. (2.2), can be characterized by transforming the coordinates to the translating frame of reference $\boldsymbol{x}(t)=\boldsymbol{y}(t)-Wt\,\boldsymbol{e}_{z}$ and considering the normal component of (2.4)

(2.5) $$\begin{eqnarray}\boldsymbol{n}_{\boldsymbol{x}}\boldsymbol{\cdot }\frac{\text{d}\boldsymbol{x}(t)}{\text{d}t}=\boldsymbol{n}_{\boldsymbol{x}}\boldsymbol{\cdot }\left[{\mathcal{C}}\,\int _{\unicode[STIX]{x2202}{\mathcal{A}}}\boldsymbol{K}(\boldsymbol{x}(t),\boldsymbol{x}^{\prime })\,\text{d}s_{\boldsymbol{ x}^{\prime }}-W\,\boldsymbol{e}_{z}\right]=0,\end{eqnarray}$$

where $\boldsymbol{n}_{\boldsymbol{x}}:=[n_{r},n_{z}]^{\text{T}}$ denotes the unit vector normal to the contour $\unicode[STIX]{x2202}{\mathcal{A}}$ at the point $\boldsymbol{x}$ and directed outward. Equation (2.5) expresses the vanishing of the normal velocity component on the vortex boundary in relative equilibrium. Further below we will also use the unit tangent vector $\boldsymbol{t}_{\boldsymbol{x}}:=[t_{r},t_{z}]^{\text{T}}$ at the point $\boldsymbol{x}$ (we will adhere to the convention that a subscript on a geometric quantity, such as the tangent or normal vector and the curvature, will indicate where this quantity is evaluated). Since we have $\boldsymbol{x}(0)=\boldsymbol{y}(0)$ , the arguments of the kernel $\boldsymbol{K}$ could be changed from $\boldsymbol{y}$ and $\boldsymbol{y}^{\prime }$ in (2.4) to $\boldsymbol{x}$ and $\boldsymbol{x}^{\prime }$ in (2.5). As the equilibrium shape of the boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ is in general a priori unknown, relation (2.5) reveals the free-boundary aspect of the problem. The Biot–Savart kernel $\boldsymbol{K}$ appearing in (2.4) and (2.5) is obtained from the standard form of this kernel in three dimensions by integrating it with respect to the azimuthal angle $\unicode[STIX]{x1D719}$ , which is justified by the axisymmetry assumption. As a result, in relations (2.4)–(2.5) integration needs to be performed in the meridional plane only. Expressions for such an axisymmetric Biot–Savart kernel were derived by Wakelin & Riley (Reference Wakelin and Riley1996) with an alternative, but equivalent, formulation also given by Pozrikidis (Reference Pozrikidis1986). Here we will use this kernel in the form rederived and tested by Shariff, Leonard & Ferziger (Reference Shariff, Leonard and Ferziger2008)

(2.6) $$\begin{eqnarray}\displaystyle & & \displaystyle \boldsymbol{K}(r\,\boldsymbol{e}_{r}+z\,\boldsymbol{e}_{z},r^{\prime }\,\boldsymbol{e}_{r}+z^{\prime }\,\boldsymbol{e}_{z})=-r^{\prime }H(r,z,r^{\prime },z^{\prime })\,n_{z}(r^{\prime },z^{\prime })\,\boldsymbol{e}_{r}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,[(z^{\prime }-z)\,G(r,z,r^{\prime },z^{\prime })\,n_{z}(r^{\prime },z^{\prime })-r\,H(r,z,r^{\prime },z^{\prime })\,n_{r}(r^{\prime },z^{\prime })]\,\boldsymbol{e}_{z},\end{eqnarray}$$

where

(2.7) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle G(r,z,r^{\prime },z^{\prime }):=\frac{r^{\prime }}{\unicode[STIX]{x03C0}\sqrt{A+B}}K(\unicode[STIX]{x1D70E}),\\ \displaystyle H(r,z,r^{\prime },z^{\prime }):=\frac{1}{2\unicode[STIX]{x03C0}r}\left[\frac{A}{\sqrt{A+B}}K(\unicode[STIX]{x1D70E})-E(\unicode[STIX]{x1D70E})\,\sqrt{A+B}\right]\end{array}\right\} & & \displaystyle\end{eqnarray}$$

in which

(2.8a-c ) $$\begin{eqnarray}A:=r^{2}+{r^{\prime }}^{2}+(z-z^{\prime })^{2},\quad B:=2rr^{\prime },\quad \unicode[STIX]{x1D70E}:=\sqrt{\frac{2B}{A+B}},\end{eqnarray}$$

whereas $K(\unicode[STIX]{x1D70E})$ and $E(\unicode[STIX]{x1D70E})$ are the complete elliptic integrals of the first and second kind, respectively (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010). We note that $\unicode[STIX]{x1D70E}\rightarrow 1$ as $r^{\prime }\rightarrow r$ and $z^{\prime }\rightarrow z$ , and at $\unicode[STIX]{x1D70E}=1$ the function $K(\unicode[STIX]{x1D70E})$ has a logarithmic singularity (more details about the singularity structure of kernel (2.6) will be provided in § 3).

Combining the kinematic condition (2.5) with the normalization condition $L=1$ , we obtain the following system of equations defining solutions of problem (2.1)–(2.2), i.e. the vortex region ${\mathcal{A}}$ and the translation velocity $W$ , as a function of the parameter (the mean core radius) $\unicode[STIX]{x1D6FC}\in (0,\sqrt{2}]$

(2.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{\unicode[STIX]{x1D6FC}^{2}}\boldsymbol{n}_{\boldsymbol{x}}\boldsymbol{\cdot }\int _{\unicode[STIX]{x2202}{\mathcal{A}}}\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x}^{\prime })\,\text{d}s_{\boldsymbol{x}^{\prime }}-W\,\boldsymbol{n}_{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{e}_{z}=0,\quad \forall \boldsymbol{x}\in \unicode[STIX]{x2202}{\mathcal{A}} & \displaystyle\end{eqnarray}$$

(2.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle r(0)+r(\unicode[STIX]{x03C0})-2=0, & \displaystyle\end{eqnarray}$$

(2.9c ) $$\begin{eqnarray}\displaystyle & \displaystyle |{\mathcal{A}}|-\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}^{2}=0. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}_{0}$

$\unicode[STIX]{x2202}{\mathcal{A}}$

$\unicode[STIX]{x1D713}$

${\mathcal{L}}$

$g$

$g$

$\unicode[STIX]{x1D6FC}\in (0,\sqrt{2})$

$\unicode[STIX]{x1D6FC}$

$\unicode[STIX]{x1D6FC}$

$\unicode[STIX]{x2202}{\mathcal{A}}$

The idea of the approach is to find a deformation of some initial reference contour such that conditions (2.9a )–(2.9c ) are satisfied. System (2.9) is solved, in a suitable discretized form, using Newton’s method in which the required contour deformation (Newton’s step) is found approximately by solving a linear system (Jacobian) obtained by differentiating (2.9) with respect to the shape of the contour $\unicode[STIX]{x2202}{\mathcal{A}}$ . Taking $\boldsymbol{x}$ to represent the contour approximating the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ at a certain iteration, its perturbation $\boldsymbol{x}^{\unicode[STIX]{x1D716}}$ obtained by applying a small but otherwise arbitrary deformation can be described as

(2.10) $$\begin{eqnarray}\boldsymbol{x}^{\unicode[STIX]{x1D716}}(\unicode[STIX]{x1D703})=\boldsymbol{x}(\unicode[STIX]{x1D703})+\unicode[STIX]{x1D716}\,\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703})\,\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703})},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}>0$ and $\unicode[STIX]{x1D70C}:[0,2\unicode[STIX]{x03C0}]\rightarrow \mathbb{R}$ is a periodic function describing the shape of the perturbation. We will adopt the convention that the superscript $\unicode[STIX]{x1D716}$ will denote quantities corresponding to the perturbed boundary, so that $\boldsymbol{x}=\boldsymbol{x}^{\unicode[STIX]{x1D716}}|_{\unicode[STIX]{x1D716}=0}$ and $\boldsymbol{n}_{\boldsymbol{x}}=\boldsymbol{n}_{\boldsymbol{x}^{\unicode[STIX]{x1D716}}}|_{\unicode[STIX]{x1D716}=0}$ are the quantities corresponding to the reference shape. Since an arbitrary contour $\unicode[STIX]{x2202}{\mathcal{A}}$ does not represent an equilibrium, system (2.9) defining the equilibrium can be equivalently expressed as the following conditions on the deformation $\unicode[STIX]{x1D70C}$ and the translation velocity $W$ obtained by setting $\unicode[STIX]{x1D716}=1$ in ansatz (2.10)

(2.11) $$\begin{eqnarray}F\left(\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}\\ W\end{array}\right]\right):=\left[\begin{array}{@{}c@{}}\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703})+\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703})\,\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703})}}\boldsymbol{\cdot }\boldsymbol{v}(\boldsymbol{x}(\unicode[STIX]{x1D703})+\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703})\,\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703})})-W\,\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703})+\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703})\,\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703})}}\boldsymbol{\cdot }\boldsymbol{e}_{z}\\ \left[r(0)+\unicode[STIX]{x1D70C}(0)\,n_{r}(0)\right]+\left[r(\unicode[STIX]{x03C0})+\unicode[STIX]{x1D70C}(\unicode[STIX]{x03C0})\,n_{r}(\unicode[STIX]{x03C0})\right]-2\\ |{\mathcal{A}}(\unicode[STIX]{x1D70C})|-\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}^{2}\end{array}\right]=\left[\begin{array}{@{}c@{}}0\\ 0\\ 0\end{array}\right],\end{eqnarray}$$

where the first component of $F$ is to be interpreted as a relation holding for all $\unicode[STIX]{x1D703}\in [0,2\unicode[STIX]{x03C0}]$ and $|{\mathcal{A}}(\unicode[STIX]{x1D70C})|$ is the area inside the perturbed contour. Infinite-dimensional (continuous) system (2.11) can be converted to an algebraic form amenable to numerical solution with Newton’s method by discretizing the parameter $\unicode[STIX]{x1D703}$ uniformly as $\unicode[STIX]{x1D703}_{j}=(j-1)\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ , $j=1,\ldots ,M$ , where $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=2\unicode[STIX]{x03C0}/M$ and $M$ is an even integer. The perturbation $\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703})$ is approximated in terms of a truncated Fourier series

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703})\approx \unicode[STIX]{x1D70C}_{N}(\unicode[STIX]{x1D703}):=\mathop{\sum }_{k=0}^{N}a_{k}\cos (k\unicode[STIX]{x1D703})+\mathop{\sum }_{k=1}^{N}b_{k}\sin (k\unicode[STIX]{x1D703}),\end{eqnarray}$$

where $N:=M/2$ . The cosine- and sine-series coefficients of the shape perturbation $\unicode[STIX]{x1D70C}_{N}$ can then be combined into a vector

(2.13) $$\begin{eqnarray}\boldsymbol{z}=\left[a_{0},a_{1},\ldots ,a_{N},b_{1},\ldots ,b_{N}\right]^{\text{T}}\in \mathbb{R}^{M+1}.\end{eqnarray}$$

After collocating the first component of (2.11), cf. (2.9a ), at the discrete parameter values $\unicode[STIX]{x1D703}_{j}$ , $j=1,\ldots ,M$ , and using suitable discretizations (to be described below), the resulting algebraic system takes the form

(2.14) $$\begin{eqnarray}\boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right)=\left[\begin{array}{@{}c@{}}\mathbf{0}\\ 0\\ 0\end{array}\right],\end{eqnarray}$$

where $\boldsymbol{F}:\mathbb{R}^{M+2}\rightarrow \mathbb{R}^{M+2}$ and $\mathbf{0}\in \mathbb{R}^{M}$ is a column vector of zeros. Relation (2.14) represents a system of nonlinear algebraic equations which can be solved using Newton’s method. Denoting $\unicode[STIX]{x1D6AB}\boldsymbol{z}$ and $\unicode[STIX]{x0394}W$ the components of Newton’s step corresponding to the updates, respectively, of the series coefficients (2.13) of the shape perturbation and of the translation velocity $W$ , a single iteration of Newton’s method requires solution of the system

(2.15a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right)\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D6AB}\boldsymbol{z}\\ \unicode[STIX]{x0394}W\end{array}\right]=-\boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right), & & \displaystyle\end{eqnarray}$$

where

(2.15b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right):=\left[\begin{array}{@{}cc@{}}\left(\unicode[STIX]{x1D648}-W\operatorname{diag}\left(\boldsymbol{e}_{z}\boldsymbol{\cdot }\boldsymbol{t}_{\boldsymbol{x}(\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{M})}\right)\unicode[STIX]{x1D729}\right), & \left(\boldsymbol{e}_{z}\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{M})}\right)\\ -[0,2,\ldots ,0,2][0,\ldots ,0], & 0\\ \left[U_{0},U_{1},\ldots ,U_{N},V_{1},\ldots ,V_{N}\right], & 0\end{array}\right] & & \displaystyle\end{eqnarray}$$

is the Jacobian of the function $\boldsymbol{F}$ , cf. (2.14), in which the expressions in the first column represent the partial derivatives of the discretized kinematic condition (2.9a ), the normalization condition (2.9b ) and the fixed-area condition (2.9c ) with respect to the expansion coefficients (2.13) obtained using methods of the shape calculus (Delfour & Zolésio Reference Delfour and Zolésio2001), cf. § 3, whereas the expressions in the second column are the partial derivatives of these conditions with respect to the translation velocity $W$ . The vectors in the first row are to be understood as $\boldsymbol{e}_{z}\boldsymbol{\cdot }\boldsymbol{t}_{\boldsymbol{x}(\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{M})}:=[t_{z}(\unicode[STIX]{x1D703}_{1}),\ldots ,t_{z}(\unicode[STIX]{x1D703}_{M})]^{\text{T}}$ and $\boldsymbol{e}_{z}\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{x}(\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{M})}:=[n_{z}(\unicode[STIX]{x1D703}_{1}),\ldots ,n_{z}(\unicode[STIX]{x1D703}_{M})]^{\text{T}}$ . Since the matrix $\unicode[STIX]{x1D648}$ also arises in the discretization of the stability problem, evaluation of its entries is discussed in detail in § 4, cf. (4.25). The term $W\operatorname{diag}(\boldsymbol{e}_{z}\boldsymbol{\cdot }\boldsymbol{t}_{\boldsymbol{x}(\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{M})})\unicode[STIX]{x1D729}$ , where the matrix $\unicode[STIX]{x1D729}$ is defined in (4.7), results from shape differentiation of the second term on the left-hand side in condition (2.9a ). Noting that the shape differential of the vortex area $|{\mathcal{A}}(\unicode[STIX]{x1D70C})|$ is $|{\mathcal{A}}(\unicode[STIX]{x1D70C})|^{\prime }=(\int _{{\mathcal{A}}}\,\text{d}{\mathcal{A}})^{\prime }=\int _{\unicode[STIX]{x2202}{\mathcal{A}}}\unicode[STIX]{x1D70C}\,\text{d}s$ (Delfour & Zolésio Reference Delfour and Zolésio2001), the partial derivatives of the fixed-area condition (2.9c ) with respect to expansion coefficients (2.13) are obtained as $U_{k}:=\int _{0}^{2\unicode[STIX]{x03C0}}\cos (k\unicode[STIX]{x1D703}^{\prime })(\text{d}s/\text{d}\unicode[STIX]{x1D703})(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }$ , $V_{k}:=\int _{0}^{2\unicode[STIX]{x03C0}}\sin (k\unicode[STIX]{x1D703}^{\prime })(\text{d}s/\text{d}\unicode[STIX]{x1D703})(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }$ , $k=0,\ldots ,M$ , where $s(\unicode[STIX]{x1D703})$ is the arclength of the contour. Details concerning a spectrally accurate approach to evaluate such expressions are provided in § 4.

The Jacobian $\unicode[STIX]{x1D735}\boldsymbol{F}$ in (2.15b ) is singular with a one-dimensional kernel space corresponding to arbitrary shifts of the contour $\unicode[STIX]{x2202}{\mathcal{A}}$ in the $z$ direction. In the process of solving (2.15a ) to compute Newton’s step the singular-value decomposition (SVD) is used to ensure that the obtained Newton’s step is orthogonal to this subspace.

We emphasize that, as will be evident in § 3, the Jacobian $\unicode[STIX]{x1D735}\boldsymbol{F}$ also contains key information necessary to set up the generalized eigenvalue problem for the study of the stability of vortex rings with respect to axisymmetric perturbations. This fact constitutes the main practical advantage of the approach developed here with respect to the method originally devised by Norbury (Reference Norbury1973). Since the vortex rings are known to possess a fore-and-aft symmetry, in ansatz (2.12) only the cosine-series coefficients will be non-zero, i.e. $b_{1}=\cdots =b_{N}=0$ , leading, via relation (2.10), to functions $r(\unicode[STIX]{x1D703})$ and $z(\unicode[STIX]{x1D703})$ which are, respectively, even and odd. However, we retain the sine series in subsequent derivations, because their outcomes will be also used in the stability analysis in § 3 where general perturbations without the fore-and-aft symmetry will be considered.

Since for all $\unicode[STIX]{x1D6FC}\in (0,\sqrt{2})$ the contours $\unicode[STIX]{x2202}{\mathcal{A}}$ are smooth (Norbury Reference Norbury1973), we can expect spectral convergence of the expansion (2.12), i.e. $|a_{k}|\sim e^{-\unicode[STIX]{x1D6FD}k}$ for $k\rightarrow \infty$ and some $\unicode[STIX]{x1D6FD}>0$ (Trefethen Reference Trefethen2013). However, since the kinematic condition (2.9a ) involves a convolution integral which is discretized using a collocation approach, one may expect that aliasing errors will arise (Boyd Reference Boyd2001). They will manifest themselves through an increase of $|a_{k}|$ with $k$ when $k>k_{a}$ for a certain $k_{a}>0$ . This issue is efficiently remedied by using an implicit filtering approach in which Newton’s step (2.15a ) is redefined as

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right)\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D642}^{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6AB}\boldsymbol{z}\\ \unicode[STIX]{x0394}W\end{array}\right]=-\boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right),\end{eqnarray}$$

in which $\unicode[STIX]{x1D642}^{\unicode[STIX]{x1D6FF}}\in \mathbb{R}^{(M+1)\times (M+1)}$ is a smoothing operator given by a low-pass filter, here represented by the following diagonal matrix

(2.17) $$\begin{eqnarray}\unicode[STIX]{x1D642}^{\unicode[STIX]{x1D6FF}}:=\operatorname{diag}\left(0,\left[\frac{1}{1+(\unicode[STIX]{x1D6FF}k)^{2p}},~k=1,\ldots ,M\right],\left[\frac{1}{1+(\unicode[STIX]{x1D6FF}k)^{2p}},~k=1,\ldots ,M\right]\right).\end{eqnarray}$$

In computations we will take $p=4$ as the order of the filter, whereas the cutoff length scale $\unicode[STIX]{x1D6FF}$ will be refined with the resolution $M$ as $\unicode[STIX]{x1D6FF}=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FE}/M$ for some $\unicode[STIX]{x1D6FE}>1$ ensuring that filtering does not affect the spectral accuracy of the approach.

As the vortex rings become ‘fatter’ in the limit $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ , the shape of their boundary $\unicode[STIX]{x2202}A$ deviates from the circle and develops ‘corners’ characterized by an increasing curvature. In order to resolve these features, increased resolutions must be used varying from $M=32$ when $\unicode[STIX]{x1D6FC}=0.1$ to $M=4096$ when $\unicode[STIX]{x1D6FC}=1.4$ . We have carefully checked the convergence of the computed shapes of the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ and of the translation velocity $W$ with refinement of the resolution $M$ . The iterations of Newton’s method are initiated with Norbury’s original solutions documented in Norbury (Reference Norbury1973), or with neighbouring solutions for values of $\unicode[STIX]{x1D6FC}$ not considered by Norbury (Reference Norbury1973). The iterations are declared converged when the following stringent termination conditions are simultaneously satisfied

(2.18a,b ) $$\begin{eqnarray}\left\Vert \boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right)\right\Vert _{2}<\widetilde{\unicode[STIX]{x1D700}}\quad \text{and}\quad \left\Vert \boldsymbol{F}\left(\left[\begin{array}{@{}c@{}}\boldsymbol{z}\\ W\end{array}\right]\right)\right\Vert _{\infty }<\widetilde{\unicode[STIX]{x1D700}},\end{eqnarray}$$

where the tolerance was set as $\widetilde{\unicode[STIX]{x1D700}}=10^{-10}$ for $\unicode[STIX]{x1D6FC}\in (0,1.2]$ and was gradually relaxed to $\widetilde{\unicode[STIX]{x1D700}}=10^{-5}$ as $\unicode[STIX]{x1D6FC}$ increased to 1.4. The reason for relaxing the tolerance is that at the higher resolutions required to represent the shapes of the vortex boundary when $\unicode[STIX]{x1D6FC}>1.2$ , the accuracy of the solution is ultimately limited by the poor algebraic conditioning of the (desingularized) Jacobian matrix (2.15b ). More specifically, the value of $\widetilde{\unicode[STIX]{x1D700}}$ roughly reflects the accuracy with which Newton’s step $[\unicode[STIX]{x1D642}^{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6AB}\boldsymbol{z}^{T},\unicode[STIX]{x0394}W]^{\text{T}}$ can be evaluated by solving problem (2.16) at a given resolution. In all cases convergence was typically achieved in $O(10)$ iterations, although due to the singularity of the Jacobian (2.15b ) convergence was only linear. In the cases when Norbury’s original solutions were used as the initial guesses, the residuals (2.18) would typically decrease by 8–10 orders of magnitude.

Figure 2. (a) Vortex boundaries $\unicode[STIX]{x2202}{\mathcal{A}}$ obtained by solving system (2.14) for $\unicode[STIX]{x1D6FC}=0.01,0.2,0.4,0.6,0.8,1.0,1.2,1.4$ . (b) The amplitudes of the Fourier coefficients $\widehat{x}_{k}$ , $k=1,\ldots ,N$ , of the corresponding parametrizations $r(\unicode[STIX]{x1D703})$ obtained for the same values of $\unicode[STIX]{x1D6FC}$ (the trend with the increase of $\unicode[STIX]{x1D6FC}$ is indicated with an arrow). The resolutions used in the computations varied from $N=16$ for small $\unicode[STIX]{x1D6FC}$ to $N=2048$ for $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ .

The shapes of the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ obtained as solutions to (2.14) for representative values of $\unicode[STIX]{x1D6FC}$ are presented in figure 2(a) where a gradual evolution of the shape from a thin vortex filament to Hill’s ‘fat’ vortex is evident as $\unicode[STIX]{x1D6FC}$ increases from 0.01 to nearly $\sqrt{2}$ . To simplify notation, hereafter the equilibrium shapes of the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ are denoted $\boldsymbol{x}=[r,z]^{\text{T}}$ , with the coordinates given in terms of parametrizations $r=r(\unicode[STIX]{x1D703})$ and $z=z(\unicode[STIX]{x1D703})$ . The amplitudes of the corresponding Fourier coefficients $\widehat{z}_{k}$ , $k=1,\ldots ,N$ , are shown in figure 2(b) (the spectra of $r(\unicode[STIX]{x1D703})$ are essentially indistinguishable and are not shown). We see that for all values of $\unicode[STIX]{x1D6FC}$ the amplitudes $|\widehat{z}_{k}|$ of the Fourier coefficients decay exponentially with the wavenumber $k$ , although this property is gradually lost as $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ where the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ develops a sharp corner (such that $r(\unicode[STIX]{x1D703})$ and $z(\unicode[STIX]{x1D703})$ are no longer smooth functions of $\unicode[STIX]{x1D703}$ ). In figure 2(b) we also observe that, interestingly, the dependence of $|\widehat{z}_{k}|$ on the wavenumber $k$ is characterized by a certain recurrent pattern. The translation velocity $W$ is positive for all values of $\unicode[STIX]{x1D6FC}$ , meaning that the vortices shown in figure 2(a) move to the left. When $\unicode[STIX]{x1D6FC}\rightarrow 0$ , this translation velocity becomes unbounded as $(1/4)(\ln (8/\unicode[STIX]{x1D6FC})-1/4)$ , whereas when $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ , $W$ approaches $4/15$ , both of which are well-known results for the propagation velocity of an infinitely thin vortex ring and Hill’s vortex, respectively (Wu et al. Reference Wu, Ma and Zhou2006). As regards the integral diagnostics, the variation of the vortex area $|{\mathcal{A}}|$ , its circulation $\unicode[STIX]{x1D6E4}$ and its hydrodynamic impulse ${\mathcal{I}}$ with $\unicode[STIX]{x1D6FC}$ is consistent with the behaviour of these quantities discussed in Norbury (Reference Norbury1973), both as $\unicode[STIX]{x1D6FC}\rightarrow 0$ and $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ . The data $\{r_{j},z_{j}\}$ , where $r_{j}=r(\unicode[STIX]{x1D703}_{j})$ and $z_{j}=z(\unicode[STIX]{x1D703}_{j})$ , $j=1,\ldots ,M$ , describing the vortex boundaries $\unicode[STIX]{x2202}{\mathcal{A}}$ for a range of different $\unicode[STIX]{x1D6FC}\in (0,\sqrt{2})$ is provided in Supplementary Information (available at https://doi.org/10.1017/jfm.2019.473) together with the corresponding values of the propagation velocity $W$ . The linear stability of these equilibrium configurations is discussed in the next section.

3 Linear stability analysis

In this section we derive a singular integro-differential equation characterizing the linear response of the boundary $\unicode[STIX]{x2202}{\mathcal{A}}(t)$ of Norbury’s vortices, cf. figure 2(a), to axisymmetric perturbations. This will be done by applying the tools of the shape-differential calculus (Delfour & Zolésio Reference Delfour and Zolésio2001) to the contour-dynamics formulation (2.4) of the problem. Finally, we obtain a constrained eigenvalue problem from which the stability properties of Norbury’s vortices can be deduced.

In addition to the circulation, impulse and energy conserved by all classical solutions of Euler’s equations, axisymmetric inviscid flows also conserve Casimirs $\iint _{{\mathcal{A}}(t)}\unicode[STIX]{x1D6F9}(\unicode[STIX]{x1D714}/r)\,r\,\text{d}{\mathcal{A}}$ , where $\unicode[STIX]{x1D6F9}:\mathbb{R}\rightarrow \mathbb{R}$ is an arbitrary function with sufficient regularity (Mohseni Reference Mohseni2001). Since circulation is a particularly simple Casimir, defined by $\unicode[STIX]{x1D6F9}(\unicode[STIX]{x1D709})=\unicode[STIX]{x1D709}$ for some $\unicode[STIX]{x1D709}\in \mathbb{R}$ , we will focus on stability analysis with respect to perturbations preserving this quantity, which is the same approach as was also taken by Moffatt & Moore (Reference Moffatt and Moore1978), Protas & Elcrat (Reference Protas and Elcrat2016) in their studies of the stability of Hill’s vortex. By this we mean that the perturbed vortex ring will have the same circulation as the corresponding vortex ring in the equilibrium configuration (however, under the linearized Euler evolution, this circulation will in general change). The circulation $\unicode[STIX]{x1D6E4}$ of the flow in the meridional plane is defined as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}:=\iint _{{\mathcal{A}}(t)}\unicode[STIX]{x1D714}\,\text{d}{\mathcal{A}},\end{eqnarray}$$

and, since the flows considered here have the property $\unicode[STIX]{x1D714}={\mathcal{C}}r$ , cf. (2.1)–(2.2), conservation of circulation (3.1) implies the conservation of the volume of the vortical region. We add that restricting the vorticity of the perturbed vortex to the form $\unicode[STIX]{x1D714}={\mathcal{C}}r$ is justified by recognizing that such configurations can be obtained as rearrangements of the vorticity distribution in the reference configuration (Brooke Benjamin Reference Brooke1975).

We now describe how axisymmetric perturbations are applied to the equilibrium configurations discussed in § 2, cf. figure 2(a). Following the ideas laid out by Elcrat & Protas (Reference Elcrat and Protas2013), we now assume that the position of the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}(t)$ depends on time $t$ and, to simplify calculations, introduce another parametrization $\boldsymbol{x}=\boldsymbol{x}(t,s)\in \unicode[STIX]{x2202}{\mathcal{A}}(t)$ for it, now in terms of the arclength $s$ of the contour. Since the maps $r=r(\unicode[STIX]{x1D703})$ and $z=z(\unicode[STIX]{x1D703})$ are assumed invertible and $\text{d}s=\sqrt{(\text{d}r)^{2}+(\text{d}x)^{2}}$ , we can rewrite ansatz (2.10) as

(3.2) $$\begin{eqnarray}\boldsymbol{x}^{\unicode[STIX]{x1D716}}(t,s)=\boldsymbol{x}(s)+\unicode[STIX]{x1D716}\,\unicode[STIX]{x1D70C}(t,s)\,\boldsymbol{n}_{\boldsymbol{x}(s)},\end{eqnarray}$$

where $s=s(\unicode[STIX]{x1D703})$ and the point $\boldsymbol{x}$ belongs to the boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ of the equilibrium configuration. We remark that for smooth reference contours $\boldsymbol{x}$ ansatz (3.2) does not restrict the class of admissible contour deformations except for the assumption that the perturbed contour $\boldsymbol{x}^{\unicode[STIX]{x1D716}}$ be ‘close’ to the reference contour $\boldsymbol{x}$ (Delfour & Zolésio Reference Delfour and Zolésio2001), which is however consistent with how the linear stability problem is formulated. Using (2.4) we then deduce

(3.3) $$\begin{eqnarray}\boldsymbol{n}_{\boldsymbol{x}^{\unicode[STIX]{x1D716}}}\boldsymbol{\cdot }\frac{\text{d}\boldsymbol{x}^{\unicode[STIX]{x1D716}}(t)}{\text{d}t}=\boldsymbol{n}_{\boldsymbol{x}^{\unicode[STIX]{x1D716}}}\boldsymbol{\cdot }\boldsymbol{v}^{\unicode[STIX]{x1D716}}(\boldsymbol{x}^{\unicode[STIX]{x1D716}}(t))=\boldsymbol{n}_{\boldsymbol{x}^{\unicode[STIX]{x1D716}}}\boldsymbol{\cdot }\left[{\mathcal{C}}\,\int _{\unicode[STIX]{x2202}{\mathcal{A}}^{\unicode[STIX]{x1D716}}(t)}\boldsymbol{K}(\boldsymbol{x}^{\unicode[STIX]{x1D716}}(t),\boldsymbol{x}^{\prime })\,\text{d}s_{\boldsymbol{ x}^{\prime }}-W\,\boldsymbol{e}_{z}\right],\end{eqnarray}$$

from which the equilibrium condition (2.5) is obtained by setting $\unicode[STIX]{x1D716}=0$ . The equation describing the evolution of the perturbation $\unicode[STIX]{x1D70C}$ is obtained by linearizing relation (3.3) around the equilibrium configuration characterized by (2.5) which is equivalent to expanding (3.3) in powers of $\unicode[STIX]{x1D716}$ and retaining the first-order terms. Since (3.3) involves perturbed quantities defined on the perturbed vortex boundary $\unicode[STIX]{x2202}A^{\unicode[STIX]{x1D716}}(t)$ , the proper way to obtain this linearization is using methods of the shape-differential calculus (Delfour & Zolésio Reference Delfour and Zolésio2001). Below we state the main results only and refer the reader to our earlier study (Elcrat & Protas Reference Elcrat and Protas2013) for details of all intermediate transformations. Shape differentiating the left-hand side of relation (3.3) and setting $\unicode[STIX]{x1D716}=0$ we obtain

(3.4) $$\begin{eqnarray}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D716}}\left.\left[\boldsymbol{n}_{\boldsymbol{x}^{\unicode[STIX]{x1D716}}}\boldsymbol{\cdot }\frac{\text{d}\boldsymbol{x}^{\unicode[STIX]{x1D716}}(t)}{\text{d}t}\right]\right|_{\unicode[STIX]{x1D716}=0}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}.\end{eqnarray}$$

As regards the right-hand side in (3.3), we obtain

(3.5) $$\begin{eqnarray}\displaystyle & & \displaystyle \left.\frac{\text{d}}{\text{d}\unicode[STIX]{x1D716}}\left\{\boldsymbol{n}_{\boldsymbol{x}^{\unicode[STIX]{x1D716}}}\boldsymbol{\cdot }\left[{\mathcal{C}}\,\int _{\unicode[STIX]{x2202}{\mathcal{A}}^{\unicode[STIX]{x1D716}}(t)}\boldsymbol{K}(\boldsymbol{x}^{\unicode[STIX]{x1D716}}(t),\boldsymbol{x}^{\prime })\,\text{d}s_{\boldsymbol{x}^{\prime }}-W\,\boldsymbol{e}_{z}\right]\right\}\right|_{\unicode[STIX]{x1D716}=0}\nonumber\\ \displaystyle & & \displaystyle \qquad =-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}s}\,\boldsymbol{v}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{t}_{\boldsymbol{x}}+{\mathcal{C}}\,\unicode[STIX]{x1D70C}(s)\,\boldsymbol{n}_{\boldsymbol{x}}\boldsymbol{\cdot }\int _{\unicode[STIX]{x2202}{\mathcal{A}}}\frac{\unicode[STIX]{x2202}\boldsymbol{K}}{\unicode[STIX]{x2202}n_{\boldsymbol{x}}}\,\text{d}s^{\prime }\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,{\mathcal{C}}\,\boldsymbol{n}_{\boldsymbol{x}}\boldsymbol{\cdot }\int _{\unicode[STIX]{x2202}{\mathcal{A}}}\boldsymbol{K}^{\prime }+\left[\frac{\unicode[STIX]{x2202}\boldsymbol{K}}{\unicode[STIX]{x2202}n_{\boldsymbol{x}^{\prime }}}+\unicode[STIX]{x1D718}_{\boldsymbol{x}^{\prime }}\boldsymbol{K}\right]\unicode[STIX]{x1D70C}(s^{\prime })\,\text{d}s^{\prime },\end{eqnarray}$$

where $\unicode[STIX]{x1D718}_{\boldsymbol{x}}$ is the curvature of the contour $\unicode[STIX]{x2202}{\mathcal{A}}$ and $\boldsymbol{K}^{\prime }$ is the shape differential of the kernel $\boldsymbol{K}$ , cf. (2.6). The orientation of the unit vectors $\boldsymbol{t}_{\boldsymbol{x}}$ and $\boldsymbol{n}_{\boldsymbol{x}}$ , and the sign of the curvature $\unicode[STIX]{x1D718}_{\boldsymbol{x}}$ satisfy Frenet’s convention. The terms in the integrand expressions in (3.5) involving the derivatives of the kernel in the directions of the normal vectors are expressed as $\unicode[STIX]{x2202}\boldsymbol{K}/\unicode[STIX]{x2202}n_{\boldsymbol{x}}=\boldsymbol{n}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x}^{\prime })$ and $\unicode[STIX]{x2202}\boldsymbol{K}/\unicode[STIX]{x2202}n_{\boldsymbol{x}^{\prime }}=\boldsymbol{n}_{\boldsymbol{x}^{\prime }}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}^{\prime }}\boldsymbol{K}(\boldsymbol{x},\boldsymbol{x}^{\prime })$ , where the subscript on the operator $\unicode[STIX]{x1D735}$ indicates which variable differentiation is performed with respect to, whereas the shape differential of the kernel is evaluated as follows

(3.6) $$\begin{eqnarray}\displaystyle & & \displaystyle \boldsymbol{K}^{\prime }(r,z,r^{\prime },z^{\prime })=-r^{\prime }H(r,z,r^{\prime },z^{\prime })\left.\frac{\text{d}}{\text{d}\unicode[STIX]{x1D716}}n_{z}^{\unicode[STIX]{x1D716}}(r^{\prime },z^{\prime })\right|_{\unicode[STIX]{x1D716}=0}\,\boldsymbol{e}_{r}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad -\,\left[(z^{\prime }-z)\,G(r,z,r^{\prime },z^{\prime })\left.\frac{\text{d}}{\text{d}\unicode[STIX]{x1D716}}n_{z}^{\unicode[STIX]{x1D716}}(r^{\prime },z^{\prime })\right|_{\unicode[STIX]{x1D716}=0}-r\,H(r,z,r^{\prime },z^{\prime })\left.\frac{\text{d}}{\text{d}\unicode[STIX]{x1D716}}n_{r}^{\unicode[STIX]{x1D716}}(r^{\prime },z^{\prime })\right|_{\unicode[STIX]{x1D716}=0}\right]\,\boldsymbol{e}_{z}\nonumber\\ \displaystyle & & \displaystyle \qquad =r^{\prime }H(r,z,r^{\prime },z^{\prime })\,n_{r}(r^{\prime },z^{\prime })\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}s}\,\boldsymbol{e}_{r}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,[(z^{\prime }-z)\,G(r,z,r^{\prime },z^{\prime })\,n_{r}(r^{\prime },z^{\prime })+r\,H(r,z,r^{\prime },z^{\prime })\,n_{x}(r^{\prime },z^{\prime })]\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}s}\,\boldsymbol{e}_{z}\nonumber\\ \displaystyle & & \displaystyle \qquad =:\boldsymbol{D}\boldsymbol{K}^{\prime }(r,z,r^{\prime },z^{\prime })\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}s},\end{eqnarray}$$

where the second equality follows from the expression for the shape derivative of the unit normal vector, i.e. $\text{d}\boldsymbol{n}^{\unicode[STIX]{x1D716}}/\text{d}\unicode[STIX]{x1D716}\,|_{\unicode[STIX]{x1D716}=0}=-(\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}s)\,\boldsymbol{t}$ , and Frenet’s identities (Elcrat & Protas Reference Elcrat and Protas2013).

As explained by Elcrat & Protas (Reference Elcrat and Protas2013), the three terms on the right-hand side of (3.5) represent the shape sensitivity of the right-hand side of (3.3) to perturbations (3.2) applied separately to the normal vector $\boldsymbol{n}_{\boldsymbol{x}}$ , the evaluation point $\boldsymbol{x}$ and the contour $\unicode[STIX]{x2202}{\mathcal{A}}$ over which the integral is defined. Since the flow evolution is subject to constraint (3.1), this will restrict the admissible perturbations $\unicode[STIX]{x1D70C}$ . Indeed, noting that $\unicode[STIX]{x1D714}={\mathcal{C}}r$ , cf. (2.1)–(2.2), and shape differentiating relation (3.1) we obtain the following condition (Elcrat & Protas Reference Elcrat and Protas2013)

(3.7) $$\begin{eqnarray}\int _{\unicode[STIX]{x2202}{\mathcal{A}}}r\unicode[STIX]{x1D70C}(0,s^{\prime })\,\text{d}s^{\prime }=0,\end{eqnarray}$$

restricting the class of admissible perturbations to those which do not modify the circulation (3.1) (analogous conditions can be obtained when the perturbations are constructed to preserve other integral invariants mentioned above). We note that the constraint imposed in the formulation of the stability problem, cf. (3.1), is different from the constraints which were imposed during the computation of the equilibrium configurations, cf. (2.9b )–(2.9c ), which resulted from the assumed normalization.

It is now convenient to return to the original parametrization in terms of $\unicode[STIX]{x1D703}$ , i.e. assume again that $r=r(\unicode[STIX]{x1D703})$ , $z=z(\unicode[STIX]{x1D703})$ , $r^{\prime }=r^{\prime }(\unicode[STIX]{x1D703}^{\prime })$ and $z^{\prime }=z^{\prime }(\unicode[STIX]{x1D703}^{\prime })$ , such that combining relations (3.4)–(3.7) and replacing the line integrals with the corresponding definite ones we finally obtain the equation describing the deformation $\unicode[STIX]{x1D70C}(t,\unicode[STIX]{x1D703})$ of the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}(t)$ , cf. (3.2), in response to axisymmetric perturbations in the linear regime

(3.8a ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t} & = & \displaystyle -{\mathcal{C}}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\left(\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}\right)^{-1}\,\int _{0}^{2\unicode[STIX]{x03C0}}I_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }+{\mathcal{C}}\,\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703})\,\int _{0}^{2\unicode[STIX]{x03C0}}I_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,{\mathcal{C}}\,\int _{0}^{2\unicode[STIX]{x03C0}}I_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }+{\mathcal{C}}\,\int _{0}^{2\unicode[STIX]{x03C0}}I_{3}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime },\end{eqnarray}$$

(3.8b ) $$\begin{eqnarray}\displaystyle & =: & \displaystyle ({\mathcal{G}}\unicode[STIX]{x1D70C})(\unicode[STIX]{x1D703})\nonumber\\ \displaystyle & & \displaystyle \text{subject to:}\quad \int _{0}^{2\unicode[STIX]{x03C0}}r(\unicode[STIX]{x1D703}^{\prime })\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }=0,\end{eqnarray}$$

${\mathcal{G}}$

(3.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle I_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime }):=\boldsymbol{t}_{\unicode[STIX]{x1D703}}\boldsymbol{\cdot }\boldsymbol{K}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime }), & \displaystyle\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle I_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime }):=\boldsymbol{n}_{\unicode[STIX]{x1D703}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{K}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })}{\unicode[STIX]{x2202}n_{\unicode[STIX]{x1D703}}}, & \displaystyle\end{eqnarray}$$

(3.9c ) $$\begin{eqnarray}\displaystyle & \displaystyle I_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime }):=\boldsymbol{n}_{\unicode[STIX]{x1D703}}\boldsymbol{\cdot }\left[\frac{\unicode[STIX]{x2202}\boldsymbol{K}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })}{\unicode[STIX]{x2202}n_{\unicode[STIX]{x1D703}^{\prime }}}+\unicode[STIX]{x1D718}_{\unicode[STIX]{x1D703}^{\prime }}\boldsymbol{K}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\right], & \displaystyle\end{eqnarray}$$

(3.9d ) $$\begin{eqnarray}\displaystyle & \displaystyle I_{3}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime }):=\boldsymbol{n}_{\unicode[STIX]{x1D703}}\boldsymbol{\cdot }\boldsymbol{D}\boldsymbol{K}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime }). & \displaystyle\end{eqnarray}$$

$\boldsymbol{K}$

$r(\unicode[STIX]{x1D703})$

$z(\unicode[STIX]{x1D703})$

$r^{\prime }(\unicode[STIX]{x1D703}^{\prime })$

$z^{\prime }(\unicode[STIX]{x1D703}^{\prime })$

As regards the singularities of the kernels, one can verify by inspection that

(3.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle \forall ~\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime }\in [0,2\unicode[STIX]{x03C0}]\quad \displaystyle \lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\,\boldsymbol{n}_{\unicode[STIX]{x1D703}}\boldsymbol{\cdot }\boldsymbol{K}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })=0, & \displaystyle\end{eqnarray}$$

(3.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\,\frac{I_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })}{\ln \sin ^{2}\left({\displaystyle \frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}}\right)}=-\frac{1}{4\unicode[STIX]{x03C0}}r(\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(3.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\,\frac{I_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })}{\ln \sin ^{2}\left({\displaystyle \frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}}\right)}=-\frac{1}{4\unicode[STIX]{x03C0}}n_{z}(\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(3.10d ) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\,\frac{I_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })}{\ln \sin ^{2}\left({\displaystyle \frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}}\right)}=\frac{1}{4\unicode[STIX]{x03C0}}n_{z}(\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(3.10e ) $$\begin{eqnarray}\displaystyle & \displaystyle \lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\,\frac{I_{3}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })}{\ln \sin ^{2}\left({\displaystyle \frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}}\right)}=-\frac{1}{4\unicode[STIX]{x03C0}}r(\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}$

$2\unicode[STIX]{x03C0}$

$\unicode[STIX]{x1D703}$

$\unicode[STIX]{x1D703}^{\prime }$

$\boldsymbol{K}$

After introducing the ansatz

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(t,\unicode[STIX]{x1D703})=\text{e}^{\text{i}\unicode[STIX]{x1D706}t}\,u(\unicode[STIX]{x1D703})+\text{C.C.},\end{eqnarray}$$

where $i:=\sqrt{-1}$ and $\unicode[STIX]{x1D706}\in \mathbb{C}$ , system (3.8) takes the form of a constrained eigenvalue problem

(3.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}\,\unicode[STIX]{x1D706}\,u(\unicode[STIX]{x1D703})=\left({\mathcal{G}}u\right)(\unicode[STIX]{x1D703}),\quad \unicode[STIX]{x1D703}\in [0,2\unicode[STIX]{x03C0}], & \displaystyle\end{eqnarray}$$

(3.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{subject to:}\quad \text{periodic boundary conditions}, & \displaystyle\end{eqnarray}$$

(3.12c ) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}r(\unicode[STIX]{x1D703}^{\prime })\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }=0, & \displaystyle\end{eqnarray}$$

${\mathcal{G}}$

$\unicode[STIX]{x1D706}$

$u$

4 Numerical approach

In this section we describe our numerical approach focusing on the discretization of system (3.12) and the solution of the resulting algebraic eigenvalue problem. We will also provide some details about how these techniques have been validated. Our approach is an adaptation to the present problem of the methodology validated and thoroughly tested in earlier studies (Elcrat & Protas Reference Elcrat and Protas2013; Protas & Elcrat Reference Protas and Elcrat2016). It also shares a number of key ingredients with the methodology employed to approximate Newton’s step (2.15a ) in the solution of problem (2.11) when computing equilibrium configurations in § 2.

As a starting point, we recall that the equilibrium configurations of the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ are represented in terms of the parametrization $[r(\unicode[STIX]{x1D703}),z(\unicode[STIX]{x1D703})]$ , $\unicode[STIX]{x1D703}\in [0,2\unicode[STIX]{x03C0}]$ , discretized (with spectral accuracy) as $\{r_{j},z_{j}\}$ , $j=1,\ldots ,M$ . The eigenvectors $u$ are then approximated with a truncated series similar to the one already used in (2.12), i.e.

(4.1) $$\begin{eqnarray}u(\unicode[STIX]{x1D703})\approx u_{N}(\unicode[STIX]{x1D703}):=\mathop{\sum }_{k=0}^{N}a_{k}\cos (k\unicode[STIX]{x1D703})+\mathop{\sum }_{k=1}^{N-1}b_{k}\sin (k\unicode[STIX]{x1D703}),\end{eqnarray}$$

(without the risk of confusion, the expansion coefficients in (2.12) and (4.1) are denoted with the same symbols). Since the total number of coefficients parameterizing $u_{N}(\unicode[STIX]{x1D703})$ in (4.1) must be equal to $M=2N$ , we arbitrarily choose to truncate the sine series in (4.1) at $k=N-1$ . After substitution of ansatz (4.1), equation (3.12a ) is collocated on the uniform grid $\unicode[STIX]{x1D703}_{1},\ldots \unicode[STIX]{x1D703}_{M}$ , whereas constraint (3.12c ) takes the form

(4.2) $$\begin{eqnarray}\displaystyle & & \displaystyle a_{0}\,\int _{0}^{2\unicode[STIX]{x03C0}}r(\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }+\mathop{\sum }_{k=1}^{N}a_{k}\,\int _{0}^{2\unicode[STIX]{x03C0}}r(\unicode[STIX]{x1D703}^{\prime })\cos (k\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }\nonumber\\ \displaystyle & & \displaystyle \qquad +\mathop{\sum }_{k=1}^{N-1}b_{k}\,\int _{0}^{2\unicode[STIX]{x03C0}}r(\unicode[STIX]{x1D703}^{\prime })\sin (k\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }=0.\end{eqnarray}$$

The factor $\text{d}s/\text{d}\unicode[STIX]{x1D703}=\sqrt{(\text{d}r(\unicode[STIX]{x1D703})/\text{d}\unicode[STIX]{x1D703})^{2}+(\text{d}z(\unicode[STIX]{x1D703})/\text{d}\unicode[STIX]{x1D703})^{2}}$ can be approximated at the collocation points $\unicode[STIX]{x1D703}_{j}$ , $j=1,\ldots ,M$ , as $\widetilde{\text{d}s/\text{d}\unicode[STIX]{x1D703}}|_{j}$ with spectral accuracy applying a suitable differentiation matrix for periodic functions (Trefethen Reference Trefethen2000) to the discrete representations of the parametrizations $r(\unicode[STIX]{x1D703})$ and $z(\unicode[STIX]{x1D703})$ . Approximating the integrals in (4.2) using the trapezoidal quadrature, which for smooth integrand expressions defined on periodic domains achieves spectral accuracy, we obtain the constraint vector $\boldsymbol{c}=[c_{1},\ldots ,c_{M}]$ , where

(4.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \mathop{\left.c_{1}:=\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\mathop{\sum }_{j=1}^{M}r_{j}\widetilde{\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}\right|}\nolimits_{j},\quad \mathop{\left.c_{1+k}:=\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\mathop{\sum }_{j=1}^{M}r_{j}\cos (k\unicode[STIX]{x1D703}_{j}^{\prime })\widetilde{\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}\right|}\nolimits_{j},\quad k=1,\ldots ,N,\\ \displaystyle \mathop{\left.c_{N+1+k}:=\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\mathop{\sum }_{j=1}^{M}r_{j}\sin (k\unicode[STIX]{x1D703}_{j}^{\prime })\widetilde{\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}\right|}\nolimits_{j},\quad k=1,\ldots ,N-1.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

We also introduce the notation $\boldsymbol{w}=[w_{1},\ldots ,w_{M}]^{\text{T}}$ , where $w_{1}:=a_{0}$ , $w_{1+k}:=a_{k}$ , $k=1,\ldots ,N$ , and $w_{N+1+k}:=b_{k}$ , $k=1,\ldots ,N-1$ . As a result, we obtain the following discretization of the generalized eigenvalue problem (3.12)

(4.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}\unicode[STIX]{x1D706}\mathop{\sum }_{k=1}^{M}A_{jk}w_{k}=\mathop{\sum }_{k=1}^{M}(B_{jk}+C_{jk}+D_{jk}+E_{jk})w_{k},\quad j=1,\ldots ,M, & \displaystyle\end{eqnarray}$$

(4.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{k=1}^{M}c_{k}w_{k}=0, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}A_{jk}:=\left\{\begin{array}{@{}ll@{}}\cos [(k-1)\unicode[STIX]{x1D703}_{j}],\quad & 1\leqslant k\leqslant N+1\\ \sin [(k-N-1)\unicode[STIX]{x1D703}_{j}],\quad & N+1<k\leqslant M\end{array}\right.,\quad j=1,\ldots ,M,\end{eqnarray}$$

is an (invertible) collocation matrix related to ansatz (4.1) and the uniform discretization of $\unicode[STIX]{x1D703}$ , whereas the matrices $\unicode[STIX]{x1D63D}$ , $\unicode[STIX]{x1D63E}$ , $\unicode[STIX]{x1D63F}$ and $\unicode[STIX]{x1D640}$ correspond to the four terms in operator ${\mathcal{G}}$ , cf. (3.8a ). We note that, as stated above, problem (4.4) appears overdetermined because there are $M+1$ conditions and only $M$ degrees of freedom $w_{1},\ldots ,w_{M}$ . However, this apparent inconsistency will be remedied further below in this section.

The entries of matrix $\unicode[STIX]{x1D63D}$ , corresponding to the first term in operator ${\mathcal{G}}$ , are defined as follows (there is no summation when indices are repeated)

(4.6) $$\begin{eqnarray}B_{jk}:={\mathcal{C}}\,\left[\left.\widetilde{\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}\right|_{j}\right]^{-1}\,\left(\int _{0}^{2\unicode[STIX]{x03C0}}I_{0}(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }\right)\,\unicode[STIX]{x1D6EF}_{jk},\quad j,k=1,\ldots ,M,\end{eqnarray}$$

where

(4.7) $$\begin{eqnarray}\unicode[STIX]{x1D6EF}_{jk}:=\left\{\begin{array}{@{}ll@{}}-(k-1)\sin [(k-1)\unicode[STIX]{x1D703}_{j}],\quad & 1\leqslant k\leqslant N+1\\ (k-N-1)\cos [(k-N-1)\unicode[STIX]{x1D703}_{j}],\quad & N+1<k\leqslant M\end{array}\right.,\quad j=1,\ldots ,M.\end{eqnarray}$$

The improper integral in (4.6), representing the tangential velocity component, cf. (3.9a ), can be evaluated using property (3.10b ) to separate the singular part of the kernel as (Hackbusch Reference Hackbusch1995)

(4.8) $$\begin{eqnarray}\displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}I_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime } & =: & \displaystyle \underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\left[I_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })+\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{{\mathcal{T}}_{0}(\unicode[STIX]{x1D703})}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D703})}.\end{eqnarray}$$

Here ${\mathcal{T}}_{0}(\unicode[STIX]{x1D703})$ has a bounded and continuous integrand and can be therefore approximated with spectral accuracy using the trapezoidal quadrature noting also that for all $\unicode[STIX]{x1D703}\in [0,2\unicode[STIX]{x03C0}]$

(4.9) $$\begin{eqnarray}\lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\left[I_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })+\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]=-\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\left[\ln \left(\frac{{\displaystyle \frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}}{r}\right)^{2}+4(1-\ln 2)\right].\end{eqnarray}$$

This and analogous relations below, cf. (4.15), (4.18) and (4.23), are obtained symbolically using the software package Maple by considering generalized series expansions of the expression on the left-hand side.

In order to evaluate the factor $\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D703})$ in the singular part of (4.8), we expand $\text{d}s/\text{d}\unicode[STIX]{x1D703}$ in a truncated trigonometric series

(4.10) $$\begin{eqnarray}\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703})=\mathop{\sum }_{k=0}^{N}\unicode[STIX]{x1D702}_{k}\cos (k\unicode[STIX]{x1D703})+\unicode[STIX]{x1D701}_{k}\sin (k\unicode[STIX]{x1D703}),\end{eqnarray}$$

where the coefficients $\{\unicode[STIX]{x1D702}_{k},\unicode[STIX]{x1D701}_{k}\}_{k=0}^{N}$ can be evaluated based on the discrete values $\widetilde{\text{d}s/\text{d}\unicode[STIX]{x1D703}}|_{j}$ , $j=1,\ldots ,M$ , using the fast Fourier transform (FFT) (while $\unicode[STIX]{x1D701}_{0}=\unicode[STIX]{x1D701}_{N}\equiv 0$ , for compactness of notation, we will include these terms in different expressions below). Then

(4.11) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D703})=\int _{0}^{2\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }=\mathop{\sum }_{k=0}^{N}\unicode[STIX]{x1D702}_{k}C_{k}(\unicode[STIX]{x1D703})+\unicode[STIX]{x1D701}_{k}S_{k}(\unicode[STIX]{x1D703}),\end{eqnarray}$$

where

(4.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle C_{k}(\unicode[STIX]{x1D703}):=\int _{0}^{2\unicode[STIX]{x03C0}}\cos (k\unicode[STIX]{x1D703}^{\prime })\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\,\text{d}\unicode[STIX]{x1D703}^{\prime }=\left\{\begin{array}{@{}ll@{}}-4\unicode[STIX]{x03C0}\ln 2,\quad & k=0\\ -{\displaystyle \frac{2\unicode[STIX]{x03C0}}{k}}\cos (k\unicode[STIX]{x1D703}),\quad & k\geqslant 1\end{array}\right., & \displaystyle\end{eqnarray}$$

(4.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle S_{k}(\unicode[STIX]{x1D703}):=\int _{0}^{2\unicode[STIX]{x03C0}}\sin (k\unicode[STIX]{x1D703}^{\prime })\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\,\text{d}\unicode[STIX]{x1D703}^{\prime }=\left\{\begin{array}{@{}ll@{}}0,\quad & k=0\\ -{\displaystyle \frac{2\unicode[STIX]{x03C0}}{k}}\sin (k\unicode[STIX]{x1D703}),\quad & k\geqslant 1\end{array}\right.. & \displaystyle\end{eqnarray}$$

$C_{k}(\unicode[STIX]{x1D703})$

$S_{k}(\unicode[STIX]{x1D703})$

The entries of matrix $\unicode[STIX]{x1D63E}$ , corresponding to the second term in operator ${\mathcal{G}}$ , are defined as follows

(4.13) $$\begin{eqnarray}C_{jk}:=\left({\mathcal{C}}\,\int _{0}^{2\unicode[STIX]{x03C0}}I_{1}(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }\right)\,A_{jk},\quad j,k=1,\ldots ,M,\end{eqnarray}$$

where the coefficient is given by an improper integral evaluated as above, i.e. using property (3.10c ) to split the integral into a regular and a singular part

(4.14) $$\begin{eqnarray}\displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}I_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime } & =: & \displaystyle \underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\left[I_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })+\frac{n_{z}(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{{\mathcal{T}}_{1}(\unicode[STIX]{x1D703})}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{n_{z}(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D703})}.\end{eqnarray}$$

The regular part ${\mathcal{T}}_{1}(\unicode[STIX]{x1D703})$ can be approximated with the trapezoidal quadrature noting that

(4.15) $$\begin{eqnarray}\displaystyle & & \displaystyle \lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\left[I_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })+\frac{n_{z}(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \qquad =\frac{-1}{4\unicode[STIX]{x03C0}\left({\displaystyle \frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}\right)^{2}}\left\{\left[6n_{r}^{2}+\ln \left(\frac{{\displaystyle \frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}}{4r}\right)^{2}\right]\left(\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}\right)^{2}n_{z}\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,2\frac{\text{d}n_{z}}{\text{d}\unicode[STIX]{x1D703}}rn_{r}\left(\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D703}}n_{z}+\frac{\text{d}r}{\text{d}\unicode[STIX]{x1D703}}n_{r}\right)+2\frac{\text{d}n_{r}}{\text{d}\unicode[STIX]{x1D703}}rn_{r}\left(\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D703}}n_{r}-\frac{\text{d}r}{\text{d}\unicode[STIX]{x1D703}}n_{z}\right)-2\frac{\text{d}n_{r}}{\text{d}\unicode[STIX]{x1D703}}\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D703}}r\nonumber\\ \displaystyle & & \displaystyle \qquad \quad \left.-\,2\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D703}}\left[\frac{\text{d}z}{\text{d}\unicode[STIX]{x1D703}}n_{z}\left(2n_{r}^{2}-1\right)+2\frac{\text{d}r}{\text{d}\unicode[STIX]{x1D703}}n_{r}\left(n_{r}^{2}-1\right)\right]\right\}=:{\mathcal{Q}}(\unicode[STIX]{x1D703}),\end{eqnarray}$$

whereas approximation of the singular part is already given by relations (4.11)–(4.12).

The entries of matrix $\unicode[STIX]{x1D63F}$ , corresponding to the third term in operator ${\mathcal{G}}$ , are defined as follows

(4.16) $$\begin{eqnarray}D_{jk}:={\mathcal{C}}\,\left\{\begin{array}{@{}ll@{}}\displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}I_{2}(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D703}^{\prime })\cos [(k-1)\unicode[STIX]{x1D703}^{\prime }]\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime },\quad & 1\leqslant k\leqslant N+1\\ \displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}I_{2}(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D703}^{\prime })\sin [(k-N-1)\unicode[STIX]{x1D703}^{\prime }]\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime },\quad & N+1<k\leqslant M\end{array}\right.,\!\!\quad j=1,\ldots ,M,\end{eqnarray}$$

which represents the action of a weakly singular integral operator on trigonometric functions and can be evaluated as above using property (3.10d ) to split the integrals in (4.16) into regular and singular parts. For the integrals involving the cosine functions we thus obtain when $k\geqslant 0$

(4.17) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}I_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\cos (k\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }\nonumber\\ \displaystyle & & \displaystyle \qquad =:\underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\left[I_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })-\frac{n_{z}(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]\cos (k\unicode[STIX]{x1D703}^{\prime })\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{{\mathcal{T}}_{2,k}(\unicode[STIX]{x1D703})}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,\frac{n_{z}(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\cos (k\unicode[STIX]{x1D703}^{\prime })\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{\unicode[STIX]{x1D6F9}_{k}(\unicode[STIX]{x1D703})},\end{eqnarray}$$

where the regular part ${\mathcal{T}}_{2,k}(\unicode[STIX]{x1D703})$ can be approximated with the trapezoidal quadrature noting that

(4.18) $$\begin{eqnarray}\lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\left[I_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })-\frac{n_{z}(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]=-{\mathcal{Q}}(\unicode[STIX]{x1D703}),\end{eqnarray}$$

where ${\mathcal{Q}}(\unicode[STIX]{x1D703})$ was defined in (4.15).

The singular part $\unicode[STIX]{x1D6F9}_{k}(\unicode[STIX]{x1D703})$ can be evaluated using the trigonometric series expansion (4.10) for $\text{d}s/\text{d}\unicode[STIX]{x1D703}$ which leads to

(4.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}_{k}(\unicode[STIX]{x1D703}) & = & \displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}\cos (k\unicode[STIX]{x1D703}^{\prime })\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{l=0}^{N}\unicode[STIX]{x1D702}_{l}\int _{0}^{2\unicode[STIX]{x03C0}}\cos (k\unicode[STIX]{x1D703}^{\prime })\cos (l\unicode[STIX]{x1D703}^{\prime })\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\,\text{d}\unicode[STIX]{x1D703}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D701}_{l}\int _{0}^{2\unicode[STIX]{x03C0}}\cos (k\unicode[STIX]{x1D703}^{\prime })\sin (l\unicode[STIX]{x1D703}^{\prime })\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\,\text{d}\unicode[STIX]{x1D703}^{\prime }\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2}\mathop{\sum }_{l=0}^{N}\unicode[STIX]{x1D702}_{l}\left[C_{l+k}(\unicode[STIX]{x1D703})+C_{l-k}(\unicode[STIX]{x1D703})\right]+\unicode[STIX]{x1D701}_{l}\left[S_{l+k}(\unicode[STIX]{x1D703})+S_{l-k}(\unicode[STIX]{x1D703})\right],\end{eqnarray}$$

where the last equality follows from the trigonometric identities

(4.20) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \cos (k\unicode[STIX]{x1D703})\cos (l\unicode[STIX]{x1D703})={\textstyle \frac{1}{2}}[\cos [(l+k)\unicode[STIX]{x1D703}]+\cos [(l-k)\unicode[STIX]{x1D703}]],\\ \displaystyle \cos (k\unicode[STIX]{x1D703})\sin (l\unicode[STIX]{x1D703})={\textstyle \frac{1}{2}}[\sin [(l+k)\unicode[STIX]{x1D703}]+\sin [(l-k)\unicode[STIX]{x1D703}]].\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The integrals in (4.16) involving the sine functions are approximated analogously.

Finally, the entries of matrix $\unicode[STIX]{x1D640}$ corresponding to the last term in operator ${\mathcal{G}}$ , are defined as follows

(4.21) $$\begin{eqnarray}\displaystyle E_{jk} & := & \displaystyle {\mathcal{C}}\,\left\{\begin{array}{@{}ll@{}}\displaystyle -(k-1)\int _{0}^{2\unicode[STIX]{x03C0}}I_{3}(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D703}^{\prime })\sin [(k-1)\unicode[STIX]{x1D703}^{\prime }]\,\text{d}\unicode[STIX]{x1D703}^{\prime },\quad & 1\leqslant k\leqslant N+1\\ \displaystyle (k-N-1)\int _{0}^{2\unicode[STIX]{x03C0}}I_{3}(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D703}^{\prime })\cos [(k-N-1)\unicode[STIX]{x1D703}^{\prime }]\,\text{d}\unicode[STIX]{x1D703}^{\prime },\quad & N+1<k\leqslant M\end{array}\right.,\nonumber\\ \displaystyle & & \displaystyle \qquad j=1,\ldots ,M,\end{eqnarray}$$

whereas before the integrals can be evaluated using property (3.10e ) to split them into regular and singular parts. For the integrals involving the sine functions in (4.21), we thus obtain when $k\geqslant 1$

(4.22) $$\begin{eqnarray}\displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}I_{3}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })\sin (k\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime } & =: & \displaystyle \underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\left[I_{3}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })+\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]\sin (k\unicode[STIX]{x1D703}^{\prime })\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{{\mathcal{T}}_{3,k}(\unicode[STIX]{x1D703})}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\underbrace{\int _{0}^{2\unicode[STIX]{x03C0}}\sin (\unicode[STIX]{x1D703}^{\prime })\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\,\text{d}\unicode[STIX]{x1D703}^{\prime }}_{S_{k}(\unicode[STIX]{x1D703})},\end{eqnarray}$$

where the regular part ${\mathcal{T}}_{3,k}(\unicode[STIX]{x1D703})$ can be approximated using the trapezoidal quadrature noting that for $\unicode[STIX]{x1D703}\in [0,2\unicode[STIX]{x03C0}]$

(4.23) $$\begin{eqnarray}\lim _{\unicode[STIX]{x1D703}^{\prime }\rightarrow \unicode[STIX]{x1D703}}\left[I_{3}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D703}^{\prime })+\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\ln \sin ^{2}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime }}{2}\right)\right]=-\frac{r(\unicode[STIX]{x1D703})}{4\unicode[STIX]{x03C0}}\left[\ln \left(\frac{{\displaystyle \frac{\text{d}s}{\text{d}\unicode[STIX]{x1D703}}}}{4r}\right)^{2}+4\right],\end{eqnarray}$$

whereas the singular part is given in (4.12b ). The integrals in (4.21) involving the cosine functions are evaluated analogously.

We now offer some comments about the validation of the approach described above. The technique in which an improper integral is approximated by splitting it into the regular and singular part, cf. (4.8), (4.14), (4.17) and (4.22), is standard (Hackbusch Reference Hackbusch1995) and was thoroughly tested on similar problems by Elcrat & Protas (Reference Elcrat and Protas2013), Protas & Elcrat (Reference Protas and Elcrat2016). Here the consistency of this approach was additionally validated by comparing the results with the results obtained by approximating the integrals on the left-hand side in (4.8), (4.14), (4.17) and (4.22) with the trapezoidal quadrature based on a staggered grid, such that the kernels were not evaluated at their singular points, i.e. for $\unicode[STIX]{x1D703}^{\prime }=\unicode[STIX]{x1D703}$ . Finally, the ultimate validation of the described numerical approach is provided by the rapid convergence of Newton’s method used to compute the equilibrium vortex configurations, cf. § 2, where the Jacobian (2.15b ) is constructed with the help of the matrices $\unicode[STIX]{x1D63D}$ , $\unicode[STIX]{x1D63E}$ , $\unicode[STIX]{x1D63F}$ and $\unicode[STIX]{x1D640}$ , cf. (4.25), used here to set up the eigenvalue problem (4.4a ) for the stability analysis. With the high precision of the numerical quadratures thus established, the shape differentiation results in (3.5) were validated by comparing them against simple forward finite-difference approximations of the shape derivatives. For example, the consistency of the first term on the right-hand side in (3.5) was checked by comparing it (as a function of $\unicode[STIX]{x1D703}$ ) to

(4.24) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{-1}(\boldsymbol{n}_{\boldsymbol{ x}^{\unicode[STIX]{x1D716}}}-\boldsymbol{n}_{\boldsymbol{x}})\boldsymbol{\cdot }\left[{\mathcal{C}}\,\int _{\unicode[STIX]{x2202}{\mathcal{A}}}\boldsymbol{K}(\boldsymbol{x}(t),\boldsymbol{x}^{\prime })\,\text{d}s_{\boldsymbol{ x}^{\prime }}-W\,\boldsymbol{e}_{z}\right],\end{eqnarray}$$

in the limit of vanishing $\unicode[STIX]{x1D716}$ . In the same spirit, the consistency of the second and third terms on the right-hand side of (3.5) was verified by perturbing the evaluation point $\boldsymbol{x}$ and the contour $\unicode[STIX]{x2202}{\mathcal{A}}$ , respectively.

Finally, we discuss solution of the constrained eigenvalue problem (4.4). We note that the problem needs to be restricted to eigenvectors satisfying condition (4.4b ) and this is done with a projection approach (Golub Reference Golub1973) which will also resolve the issue with system (4.4) being overdetermined. Defining

(4.25) $$\begin{eqnarray}\unicode[STIX]{x1D648}:=\unicode[STIX]{x1D63D}+\unicode[STIX]{x1D63E}+\unicode[STIX]{x1D63F}+\unicode[STIX]{x1D640},\end{eqnarray}$$

system (4.4a ) can be written as $\unicode[STIX]{x1D706}\,\boldsymbol{w}=(-\text{i}\boldsymbol{A}^{-1}\,\unicode[STIX]{x1D648})\,\boldsymbol{w}$ . Condition (4.4b ) restricts eigenvectors $\boldsymbol{w}$ to a subspace given by the kernel space ${\mathcal{N}}(\boldsymbol{c}):=\{\boldsymbol{w}\in \mathbb{R}^{M},\boldsymbol{c}\boldsymbol{w}=0\}$ of the constraint operator $\boldsymbol{c}$ with the corresponding projection operator given by ${\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}:=\unicode[STIX]{x1D644}-\boldsymbol{c}^{\dagger }\boldsymbol{c}$ , where $\unicode[STIX]{x1D644}$ is the $M\times M$ identity matrix and $\boldsymbol{c}^{\dagger }:=\boldsymbol{c}^{T}(\boldsymbol{c}\boldsymbol{c}^{T})^{-1}=\boldsymbol{c}^{T}/||\boldsymbol{c}||_{2}$ is the Moore–Penrose pseudo-inverse of the operator $\boldsymbol{c}$ (Laub Reference Laub2005). Applying this projection operator to the eigenvalue problem and noting that for vectors satisfying constraint (4.4b ) we have $\boldsymbol{w}={\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}\boldsymbol{w}$ , we obtain $\unicode[STIX]{x1D706}{\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}\,\boldsymbol{w}=\unicode[STIX]{x1D706}\,\boldsymbol{w}={\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}(-\text{i}\boldsymbol{A}^{-1}\unicode[STIX]{x1D648})\boldsymbol{w}$ . Next, employing the identities $\operatorname{sp}({\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}\,(-\text{i}\boldsymbol{A}^{-1}\unicode[STIX]{x1D648}))=\operatorname{sp}({\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}^{2}\,(-\text{i}\boldsymbol{A}^{-1}\unicode[STIX]{x1D648}))=\operatorname{sp}({\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}\,(-\text{i}\boldsymbol{A}^{-1}\unicode[STIX]{x1D648})\,{\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})})$ , where $\operatorname{sp}(\cdot )$ denotes the spectrum of a matrix and which are a consequence of the properties of the projection operator, we finally obtain an $M\times M$ algebraic problem equivalent to (4.4) (Golub Reference Golub1973)

(4.26) $$\begin{eqnarray}\unicode[STIX]{x1D706}\,\boldsymbol{w}={\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}(-\text{i}\boldsymbol{A}^{-1}\unicode[STIX]{x1D648}){\mathcal{P}}_{{\mathcal{N}}(\boldsymbol{c})}\,\boldsymbol{w},\end{eqnarray}$$

which can be efficiently solved in MATLAB using the function eig. Our numerical tests confirm the convergence of $K$ leading eigenvalues obtained by solving problem (4.26), where $1\ll K<2N$ , to well-defined limiting values as the resolution $M$ is refined, with analogous behaviour observed for the corresponding eigenvectors $\boldsymbol{w}$ . For each resolution $M=2N$ , a certain number $(2N-K)$ of eigenvalues with the largest moduli are spurious, i.e. they are numerical artefacts related to discretization, which is a well-known effect (Boyd Reference Boyd2001). In addition, since for $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ Norbury’s vortices have to be determined with a relaxed accuracy (cf. § 2), spurious eigenvalues with smaller magnitudes also appear for $\unicode[STIX]{x1D6FC}\geqslant 1.35$ (such spurious eigenvalues are distinguished from the true ones by the fact that they do not exhibit the expected trends as the resolution $M$ is refined and/or when the parameter $\unicode[STIX]{x1D6FC}$ is varied (Boyd Reference Boyd2001)). All of those spurious eigenvalues are excluded from the analysis in § 5. Based on such tests we can conservatively estimate the accuracy of the resolved eigenvalues as 8 significant digits for $\unicode[STIX]{x1D6FC}\rightarrow 0$ and 3 significant digits when $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ .

5 Computational results

In this section we present and analyse the results obtained by numerically solving eigenvalue problem (4.26). Our objective is to understand how the stability properties of the family of Norbury’s vortices, cf. figure 2(a), depend on the parameter $\unicode[STIX]{x1D6FC}$ . In general, two distinct regimes can be observed: when $\unicode[STIX]{x1D6FC}\in (0,\unicode[STIX]{x1D6FC}_{0}]$ , where $\unicode[STIX]{x1D6FC}_{0}\approx 0.925$ is determined approximately with the accuracy of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}=0.0125$ , all eigenvalues are purely real, i.e. $\unicode[STIX]{x1D706}_{k}=\text{Re}(\unicode[STIX]{x1D706}_{k})$ , $k\in [-K/2,\ldots ,K/2]$ , indicating that the vortices are neutrally stable; on the other hand, when $\unicode[STIX]{x1D6FC}\in (\unicode[STIX]{x1D6FC}_{0},\sqrt{2})$ , the spectra also include conjugate pairs of eigenvalues with non-vanishing imaginary parts which correspond to linearly stable and unstable eigenmodes.

As regards the first regime with purely real eigenvalues, the spectra obtained for $\unicode[STIX]{x1D6FC}=0.2,0.4,0.8$ are presented in figure 3(a). As can be observed in this figure, the difference $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})=\unicode[STIX]{x1D706}_{k+1}-\unicode[STIX]{x1D706}_{k}$ between two consecutive eigenvalues $\unicode[STIX]{x1D706}_{k+1}$ and $\unicode[STIX]{x1D706}_{k}$ (when they are ordered such that $\text{Re}(\unicode[STIX]{x1D706}_{k+1})>\text{Re}(\unicode[STIX]{x1D706}_{k})$ ) is approximately constant with respect to the index $k$ , except for the eigenvalues closest to the zero eigenvalue $\unicode[STIX]{x1D706}_{0}=0$ . It is also evident from this figure that $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})$ increases as $\unicode[STIX]{x1D6FC}$ is reduced. This effect is quantified in figure 4(a), where we can see that $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})\rightarrow \infty$ in the thin-ring limit as $\unicode[STIX]{x1D6FC}\rightarrow 0$ . We now analyse the behaviour of the eigenvalues obtained in the thin-ring limit in relation to the stability properties of Rankine’s vortex, for which the eigenvalues are given by $\unicode[STIX]{x1D706}_{k}^{R}=(\unicode[STIX]{x1D714}/2)(k-1)$ , $k\in \mathbb{Z}$ , where $\unicode[STIX]{x1D714}$ is the vorticity within the vortical region (Lamb Reference Lamb1932). For Rankine’s vortex the difference between two nearby eigenvalues is $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{R}:=|\unicode[STIX]{x1D706}_{k}^{R}-\unicode[STIX]{x1D706}_{k-1}^{R}|=(\unicode[STIX]{x1D714}/2)$ , where the vorticity corresponding to the thin-vortex limit of the vortex ring can be expressed as $\unicode[STIX]{x1D714}={\mathcal{C}}\,r=1/\unicode[STIX]{x1D6FC}^{2}$ (because as $\unicode[STIX]{x1D6FC}\rightarrow 0$ , $r\rightarrow L=1$ ). The relative difference between the quantities $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})$ , cf. figure 4, and $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{R}$ is shown as a function of $\unicode[STIX]{x1D6FC}$ in figure 5 where we see that $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})$ approaches $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{R}$ faster than linearly in $\unicode[STIX]{x1D6FC}$ as $\unicode[STIX]{x1D6FC}$ vanishes (since the spectra of Rankine’s vortex and of the thin vortex rings consist of uncountably many purely real eigenvalues, the present argument is more conveniently formulated in terms of average differences between nearby eigenvalues as done in figure 5, rather than by analysing the dependence of individual eigenvalues on $\unicode[STIX]{x1D6FC}$ ). This demonstrates that the behaviour of the eigenvalues in the thin-ring limit is in fact consistent with the stability properties of Rankine’s vortex. In the opposite limit, when $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ , the difference $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})$ vanishes, albeit rather slowly, cf. figure 4(b), demonstrating that the purely real eigenvalues become more densely packed in this limit and thus approach a continuous spectrum, which is consistent with the stability properties of Hill’s vortex (Protas & Elcrat Reference Protas and Elcrat2016).

Figure 3. Spectra of eigenvalue problem (4.26) obtained for Norbury’s vortices with (a)  $\unicode[STIX]{x1D6FC}=0.2$ (blue diamonds), $\unicode[STIX]{x1D6FC}=0.4$ (red squares), $\unicode[STIX]{x1D6FC}=0.8$ (black circles) and (b)  $\unicode[STIX]{x1D6FC}=0.975$ (blue diamonds), $\unicode[STIX]{x1D6FC}=1.2$ (red squares), $\unicode[STIX]{x1D6FC}=1.39$ (black circles). In panel (a) all eigenvalues are purely real and, for clarity, the vertical axis is not scaled uniformly.

Figure 4. Dependence of the average difference $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})$ between two consecutive purely real eigenvalues on (a) the parameter $\unicode[STIX]{x1D6FC}$ with the linear scaling and (b) the quantity $(\sqrt{2}-\unicode[STIX]{x1D6FC})$ with the logarithmic scaling. The dotted vertical lines represent the values of $\unicode[STIX]{x1D6FC}$ for which the unstable and certain neutrally stable eigenvectors are analysed in figures 7, 8 and 9.

Concerning the second regime in which linearly stable and unstable eigenmodes appear, representative spectra are shown in figure 3(b). In all cases we note the presence of a conjugate pair of purely imaginary eigenvalues, which we will denote $\unicode[STIX]{x1D706}_{1}$ and $\overline{\unicode[STIX]{x1D706}}_{1}$ , and two conjugate pairs of eigenvalues with non-zero real parts of equal magnitude, which we will denote $\unicode[STIX]{x1D706}_{2}$ , $\overline{\unicode[STIX]{x1D706}}_{2}$ , $-\unicode[STIX]{x1D706}_{2}$ , $-\overline{\unicode[STIX]{x1D706}}_{2}$ , in addition to a large number of uniformly spaced purely real eigenvalues. As regards the dependence of the imaginary parts of the eigenvalues $\unicode[STIX]{x1D706}_{1}$ and $\unicode[STIX]{x1D706}_{2}$ on $\unicode[STIX]{x1D6FC}$ , in figure 6(a) we observe that for all $\unicode[STIX]{x1D6FC}\in (\unicode[STIX]{x1D6FC}_{0},\sqrt{2})$ we have $|\text{Im}(\unicode[STIX]{x1D706}_{1})|>|\text{Im}(\unicode[STIX]{x1D706}_{2})|$ , indicating that the purely imaginary eigenvalue $\unicode[STIX]{x1D706}_{1}$ is associated with the most unstable eigenmode. In figure 6(a) we also note that $\text{Im}(\unicode[STIX]{x1D706}_{1})$ grows monotonically as $\unicode[STIX]{x1D6FC}$ increases and $\lim _{\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}}\unicode[STIX]{x1D706}_{1}=i0.2$ , which is equal to the eigenvalue associated with the most unstable eigenmode of Hill’s vortex (Moffatt & Moore Reference Moffatt and Moore1978; Protas & Elcrat Reference Protas and Elcrat2016). Concerning the dependence of $|\text{Im}(\unicode[STIX]{x1D706}_{2})|$ on $\unicode[STIX]{x1D6FC}$ , we note that, interestingly, it is not monotonic and in fact the unstable eigenvalues $\unicode[STIX]{x1D706}_{2}$ are present only when $\unicode[STIX]{x1D6FC}\in [\unicode[STIX]{x1D6FC}_{0},\unicode[STIX]{x1D6FC}_{1}]\cup [\unicode[STIX]{x1D6FC}_{2},\unicode[STIX]{x1D6FC}_{3}]\cup [\unicode[STIX]{x1D6FC}_{4},\sqrt{2})$ , where $\unicode[STIX]{x1D6FC}_{1}\approx 1.0$ , $\unicode[STIX]{x1D6FC}_{2}\approx 1.125$ , $\unicode[STIX]{x1D6FC}_{3}\approx 1.2375$ and $\unicode[STIX]{x1D6FC}_{4}\approx 1.3375$ are determined approximately with the precision of $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}=0.0125$ . This demonstrates that the associated eigenmodes are unstable only for certain values of $\unicode[STIX]{x1D6FC}$ . In figure 6(b) we see that $\text{Re}(\unicode[STIX]{x1D706}_{2})$ decreases with $\unicode[STIX]{x1D6FC}$ whenever $\text{Im}(\unicode[STIX]{x1D706}_{2})\neq 0$ . In particular, we note that $\lim _{\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}}\unicode[STIX]{x1D706}_{2}=i0.1$ , which is the eigenvalue associated with the second unstable eigenmode of Hill’s vortex (Protas & Elcrat Reference Protas and Elcrat2016). We thus conclude that the eigenvalue spectra of Norbury’s vortex rings approach the spectra of Rankine’s vortex and Hill’s spherical vortex, respectively, in the limits $\unicode[STIX]{x1D6FC}\rightarrow 0$ and $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ corresponding to the thin and fat vortex rings.

Figure 5. Dependence of the relative difference between the quantity $\unicode[STIX]{x0394}\text{Re}(\unicode[STIX]{x1D706})$ , cf. figure 4, and the quantity $\unicode[STIX]{x0394}\unicode[STIX]{x1D706}^{R}$ characterizing the stability spectrum of Rankine’s vortex on the parameter $\unicode[STIX]{x1D6FC}$ . The dotted vertical lines represent the values of $\unicode[STIX]{x1D6FC}$ for which the unstable and certain neutrally stable eigenvectors are analysed in figures 7, 8 and 9.

Figure 6. Dependence of (a) the imaginary parts of the unstable eigenvalues $\unicode[STIX]{x1D706}_{1}$ (black circles) and $\unicode[STIX]{x1D706}_{2}$ (red squares) and (b) the real part of the unstable eigenvalue $\unicode[STIX]{x1D706}_{2}$ on the parameter $\unicode[STIX]{x1D6FC}$ . The blue diamonds in panel (a) represent the eigenvalues characterizing the stability of Hill’s vortex (Moffatt & Moore Reference Moffatt and Moore1978; Protas & Elcrat Reference Protas and Elcrat2016). Since due to the reasons discussed at the end of § 4 the eigenvalue $\unicode[STIX]{x1D706}_{2}$ could not be reliably determined when $\unicode[STIX]{x1D6FC}=1.4$ , it is not shown in these figures. The dotted vertical lines represent the values of $\unicode[STIX]{x1D6FC}$ for which the unstable eigenvectors are analysed in figures 7 and 8.

Figure 7. (a) Unstable purely real eigenvectors $u_{1}$ obtained for $\unicode[STIX]{x1D6FC}=0.975,1.2,1.39$ presented as functions of the normalized arclength coordinate $s/L$ along the vortex boundary $\unicode[STIX]{x2202}A$ (the arrow indicates the trend when $\unicode[STIX]{x1D6FC}$ increases). (b–d) Deformations of the vortex boundary induced by the eigenvectors $u_{1}$ for the indicated values of $\unicode[STIX]{x1D6FC}$ (the magnitude of the perturbation, given by $\unicode[STIX]{x1D716}$ in ansatz (3.2), is arbitrary). The thick lines represent the vortex boundaries $\unicode[STIX]{x2202}{\mathcal{A}}$ in their equilibrium configurations. In this and in subsequent figures the eigenvectors are scaled such that $\max _{s\in [0,L]}|u(s)|=1$ .

Figure 8. (a) Real parts (blue solid lines) and imaginary parts (red dashed lines) of the unstable eigenvectors $u_{2}$ obtained for $\unicode[STIX]{x1D6FC}=0.975,1.2,1.39$ presented as functions of the normalized arclength coordinate $s/L$ along the vortex boundary $\unicode[STIX]{x2202}A$ (the arrow indicates the trend when $\unicode[STIX]{x1D6FC}$ increases). (b–d) Deformations of the vortex boundary induced by the real and imaginary parts of the eigenvectors $u_{2}$ for the indicated values of $\unicode[STIX]{x1D6FC}$ (the magnitude of the perturbation, given by $\unicode[STIX]{x1D716}$ in ansatz (3.2), is arbitrary). The thick lines represent the vortex boundaries $\unicode[STIX]{x2202}{\mathcal{A}}$ in their equilibrium configurations.

Figure 9. Real parts (blue solid lines) and imaginary parts (red dashed lines) of representative neutrally stable eigenvectors $u_{k}$ obtained for (a)  $\unicode[STIX]{x1D6FC}=0.1$ and (b)  $\unicode[STIX]{x1D6FC}=0.975,1.2,1.39$ presented as functions of the normalized arclength coordinate $s/L$ along the vortex boundary $\unicode[STIX]{x2202}A$ (the arrows in panel (b) indicate the trends when $\unicode[STIX]{x1D6FC}$ increases). The eigenvectors are associated with purely real eigenvalues $\unicode[STIX]{x1D706}_{k}$ close to $\pm 0.05$ . (c–f) Deformations of the vortex boundary induced by the real and imaginary parts of these neutrally stable eigenvectors for the indicated values of $\unicode[STIX]{x1D6FC}$ (the magnitude of the perturbation, given by $\unicode[STIX]{x1D716}$ in ansatz (3.2), is arbitrary). The thick lines represent the vortex boundaries $\unicode[STIX]{x2202}{\mathcal{A}}$ in their equilibrium configurations.

We now move on to discuss the eigenvectors of problem (3.12) focusing on how their structure depends on the parameter $\unicode[STIX]{x1D6FC}$ . We begin with the unstable eigenvectors corresponding to the eigenvalues $\unicode[STIX]{x1D706}_{1}$ and $\unicode[STIX]{x1D706}_{2}$ , denoted $u_{1}$ and $u_{2}$ , which we will consider for representative values of the parameter $\unicode[STIX]{x1D6FC}$ chosen from each of the three subintervals of $\unicode[STIX]{x1D6FC}$ where $\text{Im}(\unicode[STIX]{x1D706}_{2})\neq 0$ , cf. figure 6(a). We remark that when ansatz (3.11) is used and the eigenvalue $\unicode[STIX]{x1D706}$ turns out to be purely imaginary, then the imaginary part of the corresponding eigenvector $u$ is irrelevant for the stability analysis (in fact, in the present computations it vanishes identically in such cases). On the other hand, if the eigenvalue $\unicode[STIX]{x1D706}$ is complex, then both the real and imaginary part of the eigenvector $u$ play a role in the stability analysis. We add that purely imaginary eigenvalues correspond to exponential in time growth of the normal perturbations $\unicode[STIX]{x1D70C}(t,s)$ to the vortex boundary $\unicode[STIX]{x2202}{\mathcal{A}}$ , cf. ansätze (3.2) and (3.11), whereas for complex eigenvalues this exponential growth is combined with harmonic oscillations.

The unstable eigenvectors $u_{1}$ obtained for $\unicode[STIX]{x1D6FC}=0.975,1.2,1.39$ are compared as functions of the normalized arclength coordinate $s/L$ measured counterclockwise along the vortex boundary $\unicode[STIX]{x2202}A$ in figure 7(a) (the arclength coordinate is measured starting from the outer intersection of the vortex boundary $\unicode[STIX]{x2202}A$ with the $r$ -axis, cf.  $\unicode[STIX]{x1D703}=0$ in figure 1, whereas $L:=\int _{\unicode[STIX]{x2202}A}\,\text{d}s$ is the length of the vortex boundary). Deformations of the vortex boundary $\unicode[STIX]{x2202}A$ induced by these eigenvectors are shown for the different values of $\unicode[STIX]{x1D6FC}$ in figure 7(b–d). We note that for all considered values of $\unicode[STIX]{x1D6FC}$ the deformation of the vortex boundary is concentrated near its rear (right) extremity and becomes more localized as $\unicode[STIX]{x1D6FC}$ increases which is accompanied by the formation of a high-curvature region in the vortex boundary. We also observe that away from the rear extremity of the vortex the eigenvectors $u_{1}$ deform the boundary $\unicode[STIX]{x2202}A$ in an opposite sense, which is a consequence of the constraint requiring the perturbations to preserve the circulation, cf. (3.7). The real and imaginary parts of the eigenvectors $u_{2}$ associated with the second eigenvalue $\unicode[STIX]{x1D706}_{2}$ obtained for different values of $\unicode[STIX]{x1D6FC}$ are shown in figure 8(a) as functions of the normalized arclength coordinate $s/L$ , with the corresponding deformations of the vortex boundary $\unicode[STIX]{x2202}A$ illustrated in figure 8(b–d). Here, for both the real and imaginary parts of the eigenvectors we observe properties qualitatively similar to what was already noted for the eigenvectors $u_{1}$ in figure 7 with the eigenvectors $u_{2}$ also becoming increasingly localized near the rear extremity of the vortex as $\unicode[STIX]{x1D6FC}$ increases. We emphasize that this behaviour is indeed consistent with the stability properties of Hill’s vortex, corresponding to $\unicode[STIX]{x1D6FC}=\sqrt{2}$ , where the unstable eigenvectors $u_{1}$ and $u_{2}$ were found to have the form of singular distributions localized at the rear stagnation point (Protas & Elcrat Reference Protas and Elcrat2016). The unstable eigenvectors $u_{2}$ also have the interesting property that while their real and imaginary parts deform the equilibrium contour in the same sense near the rear extremity of the vortex ring, they do so in an opposite sense near the front extremity. We add that since the eigenvalues with non-zero imaginary parts come in conjugate pairs, cf. figure 3(b), for each unstable mode there is a corresponding stable mode. These stable eigenvectors have the form $u_{1}(L-s)$ and $u_{2}(L-s)$ for $s\in [0,L]$ , i.e. they are obtained from the unstable modes by reflecting them with respect to the symmetry axis of the vortex region $A$ , such that they become localized near its front (left) extremity as $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ .

Finally, we move on to discuss the dependence of the neutrally stable eigenvectors, denoted $u_{k}$ , $k\neq 1,2$ , associated with the purely real eigenvalues on the parameter  $\unicode[STIX]{x1D6FC}$ . Extrapolating the data shown in figure 3(a), we note that for all values of $\unicode[STIX]{x1D6FC}\in (0,\sqrt{2})$ there is a countable infinity of neutrally stable eigenvectors, although only a finite subset of these eigenvectors may be resolved with computations performed with a finite resolution $M$ . Thus, to fix attention, we will focus here on a few representative eigenvectors associated with eigenvalues $\unicode[STIX]{x1D706}_{k}$ close to $\pm 0.05$ . The real and imaginary parts of these eigenvectors obtained for Norbury’s vortices with $\unicode[STIX]{x1D6FC}=0.1$ and $\unicode[STIX]{x1D6FC}=0.975,1.2,1.39$ are shown as functions of the normalized arclength coordinate $s/L$ along the vortex boundary $\unicode[STIX]{x2202}A$ in figures 9(a) and 9(b), respectively, whereas the deformations of the vortex boundaries induced by these eigenvectors are illustrated in figures 9(c–f). In figures 9(a) and 9(c) we see that for small values of $\unicode[STIX]{x1D6FC}$ , when the vortex region $A$ is nearly circular, cf. figure 2(a), the eigenvectors have the form of a nearly harmonic oscillation with a small modulation and with the real and imaginary parts differing only by a phase shift. In fact, in the thin-vortex limit $\unicode[STIX]{x1D6FC}\rightarrow 0$ these eigenvectors approach Kelvin waves (i.e. purely harmonic oscillations) known to characterize the neutral stability of Rankine’s vortex (Lamb Reference Lamb1932). On the other hand, for $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ the neutrally stable eigenvectors have the form of oscillations which are increasingly localized near the front (left) and rear (right) extremity of the vortex region, a behaviour which in the limit is also consistent with the stability properties of Hill’s vortex (Protas & Elcrat Reference Protas and Elcrat2016). In general, for a given $\unicode[STIX]{x1D6FC}$ the number of oscillations exhibited by the neutrally stable eigenvectors increases with $|\text{Re}(\unicode[STIX]{x1D706}_{k})|$ .

6 Summary and conclusions

In this final section we briefly summarize and discuss our results before pointing to some open questions. Recognizing that Norbury’s vortex rings are solutions to the free-boundary problem (2.1)–(2.2), their stability to axisymmetric perturbations was studied here using an approach based on the shape-differential calculus which allows one to rigorously account for the effect of perturbations to the shape of the vortex boundary in the linearization of the governing equations, cf. § 3. In addition to axisymmetry, the perturbations are assumed not the change the circulation $\unicode[STIX]{x1D6E4}$ of the vortex ring from its equilibrium configuration. Such perturbations may be therefore regarded as resulting from the application of conservative forces to the vortex ring. Information about the stability of the vortex rings to such axisymmetric and circulation-preserving perturbations is then encoded in the spectrum of a singular integro-differential operator (3.8a ) defined on the vortex boundary $\unicode[STIX]{x2202}A$ in the meridional plane. After discretizing it with a spectrally accurate numerical technique and imposing the constraint ensuring that perturbations preserve the circulation of the vortex ring, we obtain a constrained eigenvalue problem which is solved numerically. Developed and validated by Elcrat & Protas (Reference Elcrat and Protas2013), this approach was used earlier to solve an analogous stability problem for Hill’s vortex (Protas & Elcrat Reference Protas and Elcrat2016). The base states (Norbury’s vortex rings) were recomputed using Newton’s method based on shape differentiation to ensure that they are determined with a much higher accuracy and for a broader range of the parameter $\unicode[STIX]{x1D6FC}$ than in Norbury’s original work (Norbury Reference Norbury1973). The data for these base states are provided as Supplementary Information accompanying this paper.

The results presented in § 5 demonstrate that in terms of stability properties with respect to axisymmetric perturbations the continuous family of Norbury’s vortex rings represents a ‘bridge’ connecting Rankine’s columnar vortex with Hill’s spherical vortex. Indeed, for small $\unicode[STIX]{x1D6FC}\in (0,\unicode[STIX]{x1D6FC}_{0}]$ , Norbury’s vortex rings are neutrally stable with the eigenvectors and eigenvalues approaching Kelvin waves and their eigenvalues in the thin-vortex limit $\unicode[STIX]{x1D6FC}\rightarrow 0$ . The fact that the stability properties of Rankine’s vortex are recovered in this limit can be explained by the observation that in this case a thin vortex ring can be treated as an isolated columnar (Rankine) vortex and under the assumption of axisymmetry perturbations are allowed in the meridional plane only where they give rise to Kelvin waves. As $\unicode[STIX]{x1D6FC}$ increases and the vortex rings grow fatter, they become unstable at $\unicode[STIX]{x1D6FC}_{0}\approx 0.925$ which is marked by the emergence of two unstable eigenmodes. There seems to be no obvious change to the form of the vortex region ${\mathcal{A}}$ or to the corresponding flow pattern to which this change of stability could be attributed. However, as observed by O’Farrell & Dabiri (Reference O’Farrell and Dabiri2012), the change of stability can be explained by the fact that as $\unicode[STIX]{x1D6FC}$ increases parts of the vortex boundary approach the rear stagnation point (which is on the flow axis) such that the resulting shear becomes strong enough to amplify the perturbation. On the other hand, for thin vortex rings corresponding to small values of $\unicode[STIX]{x1D6FC}$ , the shear affecting the vortex boundaries is rather weak and perturbations are swept along the vortex perimeter without amplification. The most unstable mode $u_{1}$ is associated with a purely imaginary eigenvalue $\unicode[STIX]{x1D706}_{1}$ whose imaginary part increases in magnitude monotonically with $\unicode[STIX]{x1D6FC}$ , cf. figure 6(a). Interestingly, the second unstable eigenmode $u_{2}$ is present only for certain values of the parameter $\unicode[STIX]{x1D6FC}\in [\unicode[STIX]{x1D6FC}_{0},\sqrt{2})$ . It is associated with the eigenvalue $\unicode[STIX]{x1D706}_{2}$ with a non-vanishing real part whose magnitude deceases with $\unicode[STIX]{x1D6FC}$ , cf. figure 6(b). As a result, in the fat-vortex limit $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ the two unstable eigenvalues $\unicode[STIX]{x1D706}_{1}$ and $\unicode[STIX]{x1D706}_{2}$ converge to the eigenvalues associated with the two unstable eigenvectors of Hill’s vortex, while the corresponding eigenvectors $u_{1}$ and $u_{2}$ become increasingly localized near the rear extremity of the vortex. In the limit $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ this high-curvature region becomes the rear stagnation point of Hill’s vortex. The two unstable eigenmodes of Hill’s vortex have the form of singular distributions localized at that point (Protas & Elcrat Reference Protas and Elcrat2016) and can be therefore regarded as the limits the unstable eigenvectors $u_{1}$ and $u_{2}$ of Norbury’s vortices converge to as $\unicode[STIX]{x1D6FC}\rightarrow \sqrt{2}$ , cf. figures 7(a) and 8(a). We add that analogous results concerning the stability of Hill’s vortex to axisymmetric circulation-preserving perturbations had been earlier obtained using different techniques by Moffatt & Moore (Reference Moffatt and Moore1978). At the same time, the set of purely real eigenvalues converges to the continuous spectrum which was also reported for Hill’s vortex, although the convergence is rather slow, cf. figure 4(b). Finally, we add that the properties of the stability spectra obtained here for Norbury’s vortex rings with different values of $\unicode[STIX]{x1D6FC}$ , cf. figures 3(a) and 3(b), are consistent with the Hamiltonian structure of the governing Euler equations, in the sense that in all cases the sum of the eigenvalues is zero.

The stability results summarized above are consistent with the responses of Norbury’s vortex rings to axisymmetric perturbations computed in the nonlinear regime by O’Farrell & Dabiri (Reference O’Farrell and Dabiri2012), see also Ye & Chu (Reference Ye and Chu1995). More specifically, for thin vortex rings O’Farrell & Dabiri (Reference O’Farrell and Dabiri2012) found that small perturbations to the vortex boundary travel around the vortex perimeter, which is consistent with the neutral stability reported in this regime in the present study. On the other hand, for fat vortex rings they observed that prolate perturbations of the shape of the vortex boundary are amplified and evolve into protrusions near the rear extremity of the vortex ring which then move upstream to become thin filaments. This behaviour too is consistent with the results of our analysis which predicts that fat vortex rings are linearly unstable to axisymmetric perturbations and in fact the deformation of the vortex boundary near the rear extremity of the vortex ring characterizing the most unstable eigenmodes, cf. figure 7, can be viewed as a precursor of the filament observed at later stages of the evolution by O’Farrell & Dabiri (Reference O’Farrell and Dabiri2012). Moreover, the value of the parameter $\unicode[STIX]{x1D6FC}$ demarcating the two regimes reported by O’Farrell & Dabiri (Reference O’Farrell and Dabiri2012) was $\unicode[STIX]{x1D6FC}=0.7$ , not too far from $\unicode[STIX]{x1D6FC}_{0}=0.925$ found here.

We reiterate that the physical significance of the results reported here is clearly restricted by the assumption of axisymmetry introduced to simplify the problem, which limits the class of admissible perturbations. More specifically, while the present analysis indicates that thin vortex rings are neutrally stable, they are in fact known to be linearly unstable and develop a bending instability when subject to perturbations depending on the azimuthal coordinate (Fukumoto & Hattori Reference Fukumoto and Hattori2005, see Introduction for more discussion). In the same spirit, for Hill’s vortex the studies of Fukuyu et al. (Reference Fukuyu, Ruzi and Kanai1994), Rozi (Reference Rozi1999) demonstrated that perturbations with dependence on the azimuthal angle $\unicode[STIX]{x1D719}$ which are still localized near the rear stagnation point have in fact larger growth rates than the unstable modes with axial symmetry discussed by Moffatt & Moore (Reference Moffatt and Moore1978), Protas & Elcrat (Reference Protas and Elcrat2016). Thus, the ultimate goal of this research program is to address this limitation and provide a complete understanding of the stability of inviscid vortex rings with respect to general perturbations; the present results are a stepping stone in this direction. This more general stability problem can be studied using methods analogous to the approach developed in § 3, except that the perturbation $\unicode[STIX]{x1D70C}$ will now also depend on the azimuthal angle $\unicode[STIX]{x1D719}$ , such that the eigenvalue problem corresponding to (3.12) will have both $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D719}$ as independent variables. In addition, the contour-dynamics formulation (2.4) is not available in general 3-D flows and must therefore be replaced with an expression involving a volume integral of the Biot-Savart kernel. Solving this problem will be our next step and results will be reported in the near future.