2. Two-dimensional linearized Euler equations

The Euler system (1.1)–(1.2) formulated in the moving frame of reference and linearized around an equilibrium solution $\{\psi _0,\omega _0\}$ has the following form, where $\psi ', \omega ' \, : \, (0,T] \times \varOmega \rightarrow {\mathbb {R}}$ are the perturbation variables (also defined in the moving frame of reference):

(2.1a) \begin{gather} \frac{\partial{\omega'}}{\partial{t}} ={-} (\boldsymbol{\nabla}^\perp\psi_0 - U\boldsymbol{e}_1 )\boldsymbol{\cdot} \boldsymbol{\nabla}\omega' - \boldsymbol{\nabla}^\perp\psi' \boldsymbol{\cdot} \boldsymbol{\nabla}\omega_0 \nonumber\\ ={-} (\boldsymbol{\nabla}^\perp\psi_0 - U\boldsymbol{e}_1 )\boldsymbol{\cdot} \boldsymbol{\nabla}\omega' + \boldsymbol{\nabla}\omega_0 \boldsymbol{\cdot} (\boldsymbol{\nabla}^\perp \varDelta^{{-}1})\omega' \nonumber\\ =: {\mathcal{L}} \omega', \quad \text{in} \ \varOmega, \end{gather}

(2.1b)$$\begin{gather} \Delta\psi' ={-} \omega', \quad \text{in} \ \varOmega, \end{gather}$$

(2.1c)$$\begin{gather}\psi' \rightarrow 0, \quad \text{for} \ |\boldsymbol{x}| \rightarrow \infty, \end{gather}$$

(2.1d)$$\begin{gather}\omega'(0) = w', \quad \text{in} \ \varOmega, \end{gather}$$

in which $\varDelta ^{-1}$ is the inverse Laplacian corresponding to the far-field boundary condition (2.1c) and $w'$ is an appropriate initial condition assumed to have zero circulation, i.e. $\int _{\varOmega } w'\, {\rm d}A = 0$. Unlike for problems in finite dimensions where, by virtue of the Hartman–Grobman theorem, instability of the linearized system implies the instability of the original nonlinear system, for infinite-dimensional problems this need not, in general, be the case. However, for 2-D Euler flows it was proved by Vishik & Friedlander (Reference Vishik and Friedlander2003) and Lin (Reference Lin2004) that the presence of an unstable eigenvalue in the spectrum of the linearized operator does indeed imply the instability of the original nonlinear problem.

Arnold's theory (Wu et al. Reference Wu, Ma and Zhou2006) predicts that equilibria satisfying system (1.3) are nonlinearly stable if $F'(\psi ) \ge 0$, which, however, is not the case for the Lamb–Chaplygin dipole, since using (1.4) we have $F'(\psi _0) = -b^2 < 0$ for $\psi _0 \ge \eta$. Thus, Arnold's criterion is inapplicable in this case.

2.1. Spectra of linear operators

When studying spectra of linear operators, there is fundamental difference between the finite- and infinite-dimensional cases. To elucidate this difference and its consequences, we briefly consider an abstract evolution problem ${\rm d}u/{\rm d}t = {\mathcal {A}} u$ on a Banach space ${\mathcal {X}}$ (in general, infinite-dimensional) with the state $u(t) \in {\mathcal {X}}$ and a linear operator ${\mathcal {A}} \, : \, {\mathcal {X}} \rightarrow {\mathcal {X}}$. Solution of this problem can be formally written as $u(t) = {\rm e}^{{\mathcal {A}} t} u_0$, where $u_0 \in {\mathcal {X}}$ is the initial condition and ${\rm e}^{{\mathcal {A}} t}$ the semigroup generated by ${\mathcal {A}}$ (Curtain & Zwart Reference Curtain and Zwart2013). While in finite dimensions linear operators can be represented as matrices which can only have point spectrum $\varPi _0({\mathcal {A}})$, in infinite dimensions the situation is more nuanced since the spectrum $\varLambda ({\mathcal {A}})$ of the linear operator ${\mathcal {A}}$ may in general consist of two parts, namely the approximate point spectrum $\varPi ({\mathcal {A}})$ (which is a set of numbers $\lambda \in {\mathbb {C}}$ such that $({\mathcal {A}} - \lambda )$ is not bounded from below) and the compression spectrum $\varXi ({\mathcal {A}})$ (which is a set of numbers $\lambda \in {\mathbb {C}}$ such that the closure of the range of $({\mathcal {A}} - \lambda )$ does not coincide with ${\mathcal {X}}$). We thus have $\varLambda ({\mathcal {A}}) = \varPi ({\mathcal {A}}) \cup \varXi ({\mathcal {A}})$ and the two types of spectra may overlap, i.e. $\varPi ({\mathcal {A}}) \cap \varXi ({\mathcal {A}}) \neq \emptyset$ (Halmos Reference Halmos1982). A number $\lambda \in {\mathbb {C}}$ belongs to the approximate point spectrum $\varPi ({\mathcal {A}})$ if and only if there exists a sequence of unit vectors $\{\,f_n\}$, referred to as approximate eigenvectors, such that $\|({\mathcal {A}} - \lambda ) f_n \|_{\mathcal {X}} \rightarrow 0$ as $n \rightarrow \infty$. If for some $\lambda \in \varPi ({\mathcal {A}})$ there exists a unit element $f$ such that ${\mathcal {A}} f = \lambda f$, then $\lambda$ and $f$ are an eigenvalue and an eigenvector of ${\mathcal {A}}$. The set of all eigenvalues $\lambda$ forms the point spectrum $\varPi _0({\mathcal {A}})$ which is contained in the approximate point spectrum, $\varPi _0({\mathcal {A}}) \subset \varPi ({\mathcal {A}})$. If $\lambda \in \varPi ({\mathcal {A}})$ does not belong to the point spectrum, then the sequence $\{\,f_n\}$ is weakly null convergent and consists of functions characterized by increasingly rapid oscillations as $n$ becomes large. The set of such numbers $\lambda \in {\mathbb {C}}$ is referred to as the essential spectrum $\varPi _{ess}({\mathcal {A}}) := \varPi ({\mathcal {A}}) \backslash \varPi _0({\mathcal {A}})$, a term reflecting the fact that this part of the spectrum is normally independent of boundary conditions in eigenvalue problems involving differential equations. It is, however, possible for ‘true’ eigenvalues to be embedded in the essential spectrum.

When studying the semigroup ${\rm e}^{{\mathcal {A}} t}$ one is usually interested in understanding the relation between its growth abscissa $\gamma ({\mathcal {A}}) := \lim _{t \rightarrow \infty } t^{-1} \ln \| {\rm e}^{{\mathcal {A}} t} \|_{\mathcal {X}}$ and the spectrum $\varLambda ({\mathcal {A}})$ of ${\mathcal {A}}$. While in finite dimensions $\gamma ({\mathcal {A}})$ is determined by the eigenvalues of ${\mathcal {A}}$ with the largest real part, in infinite dimensions the situation is more nuanced since there are examples in which $\sup _{z \in \varLambda ({\mathcal {A}})} {\rm Re}(z) < \gamma ({\mathcal {A}})$, e.g. Zabczyk's problem (Zabczyk Reference Zabczyk1975) also discussed by Trefethen (Reference Trefethen1997); some problems in hydrodynamic stability where such behaviour was identified are analysed by Renardy (Reference Renardy1994).

In regard to the 2-D linearized Euler operator ${\mathcal {L}}$, cf. (2.1a), it was shown by Shvydkoy & Latushkin (Reference Shvydkoy and Latushkin2003) that its essential spectrum is a vertical band in the complex plane symmetric with respect to the imaginary axis. Its width is proportional to the largest Lyapunov exponent $\lambda _{max}$ in the flow field and to the index $m \in {\mathbb {Z}}$ of the Sobolev space $H^m(\varOmega )$ in which the evolution problem is formulated (i.e. ${\mathcal {X}} = H^m(\varOmega )$ above). The norm in the Sobolev space $H^m(\varOmega )$ is defined as $\| u \|_{H^m} := [ \int _{\varOmega } \sum _{|\alpha | \le m} ({\partial ^{|\alpha |} u}/{\partial ^{\alpha _1} x_1 \partial ^{\alpha _2} x_2} )^2 \, {\rm d}A]^{1/2}$, where $\alpha _1,\alpha _2 \in {\mathbb {Z}}$ with $|\alpha | := \alpha _1+\alpha _2$ (Adams & Fournier Reference Adams and Fournier2005). More specifically, we have (Shvydkoy & Friedlander Reference Shvydkoy and Friedlander2005)

(2.2)\begin{equation} \varPi_{ess}({\mathcal{L}}) = \{ z \in {\mathbb{C}}, \ -|m| \lambda_{max} \le {\rm Re}(z) \le |m| \lambda_{max} \}. \end{equation}

In 2-D flows Lyapunov exponents are determined by the properties of the velocity gradient $\boldsymbol {\nabla }\boldsymbol {u}(\boldsymbol {x})$ at hyperbolic stagnation points $\boldsymbol {x}_0$. More precisely, $\lambda _{max}$ is given by the largest eigenvalue of $\boldsymbol {\nabla }\boldsymbol {u}(\boldsymbol {x})$ computed over all stagnation points. As regards the Lamb–Chaplygin dipole, it is evident from figure 1(a,b) that in both the symmetric and asymmetric cases it has two stagnation points $\boldsymbol {x}_a$ and $\boldsymbol {x}_b$ located at the fore and aft extremities of the vortex core. Inspection of the velocity field $\boldsymbol {\nabla }^\perp \psi _0$ defined in (1.5a) shows that the largest eigenvalues of $\boldsymbol {\nabla }\boldsymbol {u}(\boldsymbol {x})$ evaluated at these stagnation points, and hence the Lyapunov exponents, are $\lambda _{max} = 2$ regardless of the value of $\eta$.

While characterization of the essential spectrum of the 2-D linearized Euler operator ${\mathcal {L}}$ is rather complete, the existence of a point spectrum remains in general an open problem. Results concerning the point spectrum are available in a few cases only, usually for shear flows where the problem can be reduced to one dimension (Chandrasekhar Reference Chandrasekhar1961; Drazin & Reid Reference Drazin and Reid1981) or the cellular cat's eyes flows (Friedlander, Vishik & Yudovich Reference Friedlander, Vishik and Yudovich2000). In these examples unstable eigenvalues are outside the essential spectrum (if one exists) and the corresponding eigenfunctions are well behaved. On the other hand, it was shown by Lin (Reference Lin2004) that when an unstable eigenvalue is embedded in the essential spectrum, then the corresponding eigenfunctions need not be smooth. One of the goals of the present study is to consider this issue for the Lamb–Chaplygin dipole.

2.2. Linearization around the Lamb–Chaplygin dipole

The linear system (2.1) is defined on the entire plane ${\mathbb {R}}^2$; however, in the Lamb–Chaplygin dipole the vorticity $\omega _0$ is supported within the vortex core $A_0$ only, cf. (1.6b). This will allow us to simplify system (2.1) so that it will involve relations defined only within $A_0$, which will facilitate both the asymptotic analysis and numerical solution of the corresponding eigenvalue problem, cf. §§ 3 and 5. If the initial data $w'$ in (2.1d) are also supported in $A_0$, then the initial-value problem (2.1) can be regarded as a free-boundary problem describing the evolution of the boundary $\partial A(t)$ of the vortex core (we have $A(0) = A_0$ and $\partial A(t) = \partial A_0$). However, as explained below, the evolution of this boundary can be deduced from the evolution of the perturbation streamfunction $\psi '(t,\boldsymbol {x})$, hence need not be tracked independently. Thus, the present problem is different from, for example, the vortex-patch problem where the vorticity distribution is fixed (piecewise constant in space) and in the stability analysis the boundary is explicitly perturbed (Elcrat & Protas Reference Elcrat and Protas2013).

Denoting $\psi _1' \, : \, (0,T] \times A_0 \rightarrow {\mathbb {R}}$ and $\psi _2' \, : \, (0,T] \times {\mathbb {R}}^2\backslash \bar {A}_0 \rightarrow {\mathbb {R}}$ the perturbation streamfunction in the vortex core and in its complement, system (2.1) can be recast as

(2.3a)$$\begin{align} \frac{\partial{\omega'}}{\partial{t}} &={-} (\boldsymbol{\nabla}^\perp\psi_0 - U\boldsymbol{e}_1 )\boldsymbol{\cdot} \boldsymbol{\nabla}\omega' - \boldsymbol{\nabla}^\perp\psi_1' \boldsymbol{\cdot} \boldsymbol{\nabla}\omega_0, \quad \text{in} \ A_0, \end{align}$$

(2.3b)$$\begin{align}\Delta\psi_1' &={-} \omega', \quad \text{in} \ A_0, \end{align}$$

(2.3c)$$\begin{align}\Delta\psi_2' &= 0, \quad \text{in} \ {\mathbb{R}}^2\backslash \bar{A}_0, \end{align}$$

(2.3d)$$\begin{align}\psi_1' &= \psi_2' = f', \quad \text{on} \ \partial A_0, \end{align}$$

(2.3e)$$\begin{align}\frac{\partial{\psi_1'}}{\partial{n}} &= \frac{\partial{\psi_2'}}{\partial{n}} , \quad \text{on} \ \partial A_0, \end{align}$$

(2.3f)$$\begin{align}\psi_2' &\rightarrow 0 , \quad \text{for} \ |\boldsymbol{x}| \rightarrow \infty, \end{align}$$

(2.3g)$$\begin{align}\omega'(0) &= w', \quad \text{in} \ A_0, \end{align}$$

where $\boldsymbol {n}$ is the unit vector normal to the boundary $\partial A_0$ pointing outside and conditions (2.3d)–(2.3e) represent the continuity of the normal and tangential perturbation velocity components across the boundary $\partial A_0$ with $f' \, : \, \partial A_0 \rightarrow {\mathbb {R}}$ denoting the unknown value of the perturbation streamfunction at that boundary.

The velocity normal to the vortex boundary $\partial A(t)$ is $u_n := \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {n} = \partial \psi _1 / \partial s = \partial \psi _2 / \partial s$, where $s$ is the arc-length coordinate along $\partial A(t)$, cf. (2.3d). While this quantity identically vanishes in the equilibrium state (1.5)–(1.6), cf. (2.11), in general it will be non-zero resulting in a deformation of the boundary $\partial A(t)$. This deformation can be deduced from the solution of system (2.3) as follows. Given a point $\boldsymbol {z} \in \partial A(t)$, the deformation of the boundary is described by ${\rm d}\boldsymbol {z} / {\rm d}t = \boldsymbol {n} u_n|_{\partial A(t)}$. Integrating this expression with respect to time yields

(2.4)\begin{equation} \boldsymbol{z}(\tau) = \boldsymbol{z}(0) + \int_0^\tau \boldsymbol{n} u_n|_{\partial A_\tau} \, {\rm d}\tau' = \boldsymbol{z}(0) + \tau \boldsymbol{n} u_n|_{\partial A_0} + \mathcal{O}(\tau^2), \end{equation}

where $\boldsymbol {z}(0) \in \partial A_0$ and $0 < \tau \ll 1$ is the time over which the deformation is considered. Thus, the normal deformation of the boundary can be defined as $\rho (\tau ) := \boldsymbol {n} \boldsymbol {\cdot }[ \boldsymbol {z}(\tau ) - \boldsymbol {z}(0)] \approx u_n|_{\partial A_0} \tau$. We also note that at the leading order the area of the vortex core $A(t)$ is preserved by the considered perturbations:

(2.5)\begin{equation} \oint_{\partial A_0} \rho(\tau) \, {\rm d}s = \tau \oint_{\partial A_0} \frac{\partial{\psi}}{\partial{s}} \, {\rm d}s = \tau \oint_{\partial A_0} \, {\rm d}\psi = 0 \Longrightarrow\ |A(t)| \approx |A_0|. \end{equation}

We notice that in the exterior domain ${\mathbb {R}}^2\backslash \bar {A}_0$ the problem is governed by Laplace's equation (2.3c) subject to boundary conditions (2.3d)–(2.3f). Therefore, this subproblem can be eliminated by introducing the corresponding Dirichlet-to-Neumann (D2N) map $M \, : \, \psi _2'|_{\partial A_0} \rightarrow \partial \psi _2'/\partial n|_{\partial A_0}$ which is constructed in an explicit form in Appendix A. Thus, (2.3c) with boundary conditions (2.3d)–(2.3f) can be replaced with a single relation $\partial \psi _1'/\partial n = M\psi _1'$ holding on $\partial A_0$ such that the resulting system is defined in the vortex core $A_0$ and on its boundary only. It should be emphasized that this reduction is exact as the construction of the D2N map does not involve any approximations. We therefore conclude that while the vortex boundary $\partial A(t)$ may deform in the course of the linear evolution, this deformation can be described based solely on quantities defined within $A_0$ and on $\partial A_0$ using relation (2.4). In particular, the transport of vorticity out of the vortex core $A_0$ into the potential flow is described by the last term on the right-hand side in (2.3a) evaluated on the boundary $\partial A_0$.

Noting that the base state satisfies the equation $\Delta \psi _0 = - b^2 (\psi _0 - \eta )$ in $A_0$, cf. (1.3)–(1.4), and using the identity $(\boldsymbol {\nabla }^\perp \psi _1')\boldsymbol {\cdot } \boldsymbol {\nabla }\psi _0 = - (\boldsymbol {\nabla } \psi _1')\boldsymbol {\cdot } \boldsymbol {\nabla }^\perp \psi _0$, the vorticity equation (2.3a) can be transformed to the following simpler form:

(2.6)\begin{equation} \frac{\partial{\Delta \psi_1'}}{\partial{t}} ={-} (\boldsymbol{\nabla}^\perp \psi_0) \boldsymbol{\cdot} \boldsymbol{\nabla}(\Delta \psi_1' + b^2 \psi_1') \quad \text{in} \ A_0, \end{equation}

where we also used (2.3b) to eliminate $\omega '$ in favour of $\psi _1'$. Supposing the existence of an eigenvalue $\lambda \in {\mathbb {C}}$ and an eigenfunction ${\tilde {\psi }} \, : \, A_0 \rightarrow {\mathbb {C}}$, we make the following ansatz for the perturbation streamfunction $\psi _1'(t,\boldsymbol {x}) = {\tilde {\psi }}(\boldsymbol {x}) \, {\rm e}^{\lambda t}$ which leads to the eigenvalue problem:

(2.7a)$$\begin{align} \lambda {\tilde{\psi}} &= \varDelta_M^{{-}1} [ (\boldsymbol{\nabla}^\perp \psi_0) \boldsymbol{\cdot} \boldsymbol{\nabla}(\varDelta {\tilde{\psi}} + b^2 {\tilde{\psi}}) ] \quad \text{in} \ A_0, \end{align}$$

(2.7b)$$\begin{align}\frac{\partial{\Delta {\tilde{\psi}}}}{\partial{r}} &= 0, \quad \text{at} \ r=0, \end{align}$$

where $\varDelta _M^{-1}$ is the inverse Laplacian subject to the boundary condition $\partial {\tilde {\psi }} / \partial n - M{\tilde {\psi }} = 0$ imposed on $\partial A_0$ and the additional boundary condition (2.7b) ensures the perturbation vorticity is differentiable at the origin (such condition is necessary since the differential operator on the right-hand side in (2.7a) is of order three). Depending on whether or not the different differential operators appearing in it are inverted, eigenvalue problem (2.7) can be rewritten in a number of different, yet mathematically equivalent, forms. However, all these alternative formulations have the form of generalized eigenvalue problems and are therefore more difficult to handle in numerical computations. Thus, formulation (2.7) is preferred and we focus on it hereafter.

We note that the proposed formulation ensures that the eigenfunctions ${\tilde {\psi }}$ have zero circulation, as required:

(2.8)\begin{align} \varGamma' &:= \int_{A_0} \omega' \, {\rm d}A ={-} \int_{A_0} \Delta \psi_1' \, {\rm d}A ={-} \oint_{\partial A_0} \frac{\partial{\psi_1'}}{\partial{n}} \, {\rm d}s\nonumber\\ & ={-} \oint_{\partial A_0} \frac{\partial{\psi_2'}}{\partial{n}} \, {\rm d}s ={-} \int_{{\mathbb{R}}^2 \backslash \bar{A}_0} \Delta \psi_2' \, {\rm d}A = 0, \end{align}

where we used the divergence theorem, (2.3b)–(2.3c) and the boundary conditions (2.3e)–(2.3f).

Since it will be needed for the numerical discretization described in § 5, we now rewrite the eigenvalue problem (2.7) explicitly in the polar coordinate system:

(2.9a)$$\begin{align} \lambda {\tilde{\psi}} &= \varDelta_M^{{-}1}\left[ \left( u_0^r \frac{\partial{}}{\partial{r}} + \frac{u_0^\theta}{r} \frac{\partial{}}{\partial{\theta}} \right) (\varDelta + b^2 ) {\tilde{\psi}} \right] =: {\mathcal{H}} {\tilde{\psi}} \quad \text{for}\ 0 < r \le 1, \ 0 \le \theta \le 2{\rm \pi}, \end{align}$$

(2.9b)$$\begin{align}\frac{\partial{\Delta {\tilde{\psi}}}}{\partial{r}} &= 0, \quad \text{at} \ r=0, \end{align}$$

where $\varDelta = \partial ^2/\partial r^{2} + ({1}/{r})(\partial / \partial r) + ({1}/{r^2})(\partial ^2/\partial \theta ^2)$ and the velocity components obtained as $[ u_0^r, u_0^\theta ] := \boldsymbol {\nabla }^\perp \psi _0 = [ ({1}/{r})(\partial /\partial \theta ), -\partial /\partial r] \psi _0$ are

(2.10a)$$\begin{align} u_0^r &= \frac{2U J_1(br)\cos\theta}{b J_0(b) r}, \end{align}$$

(2.10b)$$\begin{align}u_0^\theta &={-} \frac{2U \left[ J_0(b r) - \dfrac{J_1(b r)}{b r} \right]\sin\theta + \eta b J_1(b r)}{J_0(b)}. \end{align}$$

They have the following behaviour on the boundary $\partial A_0$:

(2.11a,b)\begin{equation} u_0^r(1,\theta) = 0, \quad u_0^\theta(1,\theta) = 2U \sin\theta. \end{equation}

Since $\|\psi \|_{L^2} \sim \| \Delta \psi \|_{H^{-2}} = \| \omega \|_{H^{-2}}$, where ‘$\sim$’ means the norms on the left and on the right are equivalent (in the precise sense of norm equivalence), the essential spectrum (2.2) of the operator ${\mathcal {H}}$ will have $m = - 2$, so that $\varPi _{ess}({\mathcal {H}})$ is a vertical band in the complex plane with $|{\rm Re}(z)| \le 4$, $z \in {\mathbb {C}}$ (since $\lambda _{max} = 2$).

Operator ${\mathcal {H}}$ (cf. (2.9a)) has a non-trivial null space $\operatorname {Ker}({\mathcal {H}})$. To see this, we consider the ‘outer’ subproblem

(2.12a)$$\begin{align} {\mathcal{K}} \phi &:= \left( u_0^r \frac{\partial{}}{\partial{r}} + \frac{u_0^\theta}{r} \frac{\partial{}}{\partial{\theta}} \right) \phi = 0 \quad \text{for} \ 0 < r \le 1, \ 0 \le \theta \le 2{\rm \pi}, \end{align}$$

(2.12b)$$\begin{align}\frac{\partial{\phi}}{\partial{r}} &= 0, \quad \text{at} \ r=0, \end{align}$$

whose solutions are $\phi (r,\theta ) = \phi _C(r,\theta ) := B [J_1(br) \sin \theta ]^C$, $B \in {\mathbb {R}}$, $C = 2,3,\ldots$ (see Appendix B for derivation details). Then, the eigenfunctions ${\tilde {\psi }}_C$ spanning the null space of operator ${\mathcal {H}}$ are obtained as solutions of the family of ‘inner’ subproblems

(2.13a)$$\begin{gather} \left( \frac{\partial^{2}}{\partial{r}^{2}} + \frac{1}{r} \frac{\partial{}}{\partial{r}} + \frac{1}{r^2} \frac{\partial^{2}}{\partial{\theta}^{2}} + b^2 \right) {\tilde{\psi}}_C = \phi_C \quad \text{for} \ 0 < r \le 1, \ 0 \le \theta \le 2{\rm \pi}, \end{gather}$$

(2.13b)$$\begin{gather}\frac{\partial{{\tilde{\psi}}_C}}{\partial{r}} + M {\tilde{\psi}}_C = 0, \quad \text{at} \ r=0, \end{gather}$$

where $C = 2,3,\ldots$. Some of these eigenfunctions are shown in figure 2(a–d), where distinct patterns are evident for even and odd values of $C$.

Figure 2. Eigenfunctions ${\tilde {\psi }}_C$, $C=$ (a) 2, (b) 3, (c) 4, (d) 5, corresponding to the zero eigenvalue of problem (2.9).

3. Asymptotic solution of eigenvalue problem (2.9)

A number of interesting insights about certain properties of solutions of eigenvalue problem (2.9) can be deduced by performing a simple asymptotic analysis of this problem in the short-wavelength limit. We focus here on the case of the symmetric dipole ($\eta = 0$) and begin by introducing the ansatz

(3.1)\begin{equation} {\tilde{\psi}}(r,\theta) = \sum_{m=0}^{\infty} f_m(r) \cos(m\theta) + g_m(r)\sin(m\theta), \end{equation}

where $f_m, g_m \, : \, [0,1] \rightarrow {\mathbb {C}}$, $m=1,2,\ldots$, are functions to be determined. Substituting this ansatz in (2.9a) with the Laplacian moved back to the left-hand side and applying well-known trigonometric identities leads after some algebra to the following system of coupled third-order ordinary differential equations for the functions $f_m(r)$, $m=1,2,\ldots$:

(3.2a) \begin{gather} \lambda {\mathcal{B}}_m f_m = \frac{1}{2}P(r) \frac{{\rm d}}{{\rm d}r} ( {\mathcal{B}}_{m-1}f_{m-1} + b^2 f_{m-1} + {\mathcal{B}}_{m+1}f_{m+1} + b^2 f_{m+1} ) \nonumber\\ \times \frac{1}{2}Q(r) \frac{m}{r} ( {\mathcal{B}}_{m-1}f_{m-1} + b^2 f_{m-1} - {\mathcal{B}}_{m+1}f_{m+1} - b^2 f_{m+1} ), \quad r \in (0,1), \end{gather}

(3.2b)$$\begin{gather} f_m \quad \text{bounded}\quad \text{at} \ r = 0, \end{gather}$$

(3.2c)$$\begin{gather}\frac{{\rm d}}{{\rm d}r} f_m ={-} m f_m \quad \text{at} \ r = 1, \end{gather}$$

(3.2d)$$\begin{gather}\frac{{\rm d}}{{\rm d}r} {\mathcal{B}}_m f_m = 0 \quad \text{at} \ r = 0, \end{gather}$$

where the Bessel operator ${\mathcal {B}}_m$ is defined via ${\mathcal {B}}_m f := ({{\rm d}^2}/{{\rm d}r^2}) f + ({1}/{r}) ({{\rm d}}/{{\rm d}r}) f - ({m^2}/{r}) f$, whereas the coefficient functions have the following form, cf. (2.10):

(3.3a)$$\begin{gather} P(r) := \frac{2U J_1(br)}{b J_0(b) r}, \end{gather}$$

(3.3b)$$\begin{gather}Q(r) :={-} \frac{2U \left[ J_0(b r) - \dfrac{J_1(b r)}{b r} \right]}{J_0(b)}. \end{gather}$$

The functions $g_m(r)$, $m=1,2,\ldots$, satisfy a system identical to (3.2), which shows that the eigenfunctions ${\tilde {\psi }}(r,\theta )$ are either even or odd functions of $\theta$ (i.e. they are either symmetric or antisymmetric with respect to the flow centreline). Moreover, the fact that system (3.2) couples Fourier components corresponding to different $m$ implies that the eigenvectors ${\tilde {\psi }}(r,\theta )$ are not separable as functions of $r$ and $\theta$.

Motivated by our discussion in § 2.1 about the properties of approximate eigenfunctions of the 2-D linearized Euler operator, we construct approximate solutions of system (3.2) in the short-wavelength limit $m \rightarrow \infty$. In this analysis we assume that $\lambda \in \varPi _{ess}({\mathcal {L}})$ is given and focus on the asymptotic structure of the corresponding approximate eigenfunctions. We thus consider the asymptotic expansions

(3.4a,b)\begin{equation} \lambda = \lambda^0 + \frac{1}{m} \lambda^1 + {\mathcal{O}}\left(\frac{1}{m^2}\right), \quad f_m(r) = f_m^0(r) + \frac{1}{m} f_m^1(r) + {\mathcal{O}}\left(\frac{1}{m^2}\right), \end{equation}

where $\lambda ^0, \lambda ^1 \in {\mathbb {C}}$ are treated as parameters and $f_m^0, f_m^1 \, : \, [0,1] \rightarrow {\mathbb {C}}$ are unknown functions. Plugging these expansions into system (3.2) and collecting terms proportional to the highest powers of $m$, we obtain

(3.5a) \begin{align} &{\mathcal{O}}(m^3): \quad\quad\quad\quad f_{m-1}^0 - f_{m+1}^0 = 0, \end{align}

(3.5b) \begin{align} &{\mathcal{O}}(m^2): \quad \frac{1}{2} \frac{Q(r)}{r^3}(\,f_{m-1}^1 - f_{m+1}^1 ) =\frac{1}{2} P(r) \frac{{\rm d}}{{\rm d}r}\left[ \frac{1}{r^2} (\,f_{m-1}^0 + f_{m+1}^0) \right] \nonumber\\ & + \frac{Q(r)}{r^3}(\,f_{m-1}^0 + f_{m+1}^0 ) + \frac{\lambda^0}{r^2} f_m^0. \end{align}

It follows immediately from (3.5a) that $f_{m-1}^0 = f_{m+1}^0$. Since this analysis does not distinguish between even and odd values of $m$, we also deduce that $f_m^0 = f_{m-1}^0 = f_{m+1}^0$, such that relation (3.5b) takes the form

(3.6)\begin{equation} P(r) \frac{{\rm d}}{{\rm d}r}\left(\frac{1}{r^2} f_m^0 \right) - 2 \frac{Q(r)}{r^3} f_m^0 - \frac{\lambda^0}{r^2} f_m^0 = \frac{1}{2} \frac{Q(r)}{r^3}(\,f_{m-1}^1 - f_{m+1}^1 ), \quad r \in (0,1), \end{equation}

which is an inhomogeneous first-order equation defining the leading-order term $f_m^0(r)$ in (3.4) in terms of $f_m^1(r)$. Without loss of generality the boundary condition (3.2b) can be replaced with $f_m^0(0) = 1$. The solution of (3.6) is then a sum of two parts: the solution of the homogeneous equation obtained by setting the right-hand side to zero and a particular integral corresponding to the actual right-hand side. Since at this level the expression $f_{m-1}^1 - f_{m+1}^1$ is undefined, we cannot find the particular integral. On the other hand, the solution of the homogeneous equation can be found directly noting that this equation is separable and integrating which gives

(3.7)\begin{equation} f_m^0(r) = \exp\left[ {\rm i} \int_0^r I_i(r') \, {\rm d}r' \right] \exp\left[ \int_0^r I_r(r') \, {\rm d}r' \right],\quad r \in [0,1], \end{equation}

where

(3.8a)$$\begin{gather} I_r(r) := \frac{{\rm Re}(\lambda^0) b J_0(b) r^2 - 4 U b J_0(br) r + 8 U J_1 (br)}{2 U J_1(b r) r }, \end{gather}$$

(3.8b)$$\begin{gather}I_i(r) := \frac{{\rm Im}(\lambda^0) b J_0(b) r}{2 U J_1(b r)}. \end{gather}$$

The limiting (as $r \rightarrow 1$) behaviour of functions (3.8a)–(3.8b) exhibits an interesting dependence on $\lambda ^0$, namely:

(3.9a)$$\begin{gather} \lim_{r \rightarrow 1} I_r(r) =\begin{cases} + \infty, & {\rm Re}(\lambda^0) < 4 \\ \phantom+ 0, & {\rm Re}(\lambda^0) = 4 \\ - \infty, & {\rm Re}(\lambda^0) > 4 , \end{cases} \end{gather}$$

(3.9b)$$\begin{gather}\lim_{r \rightarrow 1} I_i(r) ={-} \operatorname{sign}[{\rm Im}(\lambda^0)] \infty. \end{gather}$$

In particular, the limiting value of $I_r(r)$ as $r \rightarrow 1$ changes when ${\rm Re}(\lambda ^0) = 4$, which defines the right boundary of the essential spectrum in the present problem, cf. (2.2). Both $I_r(r)$ and $I_i(r)$ diverge as ${\mathcal {O}}(1/(1 - r))$ when $r \rightarrow 1$ which means that the integrals under the exponentials in (3.7), and hence the entire formula, are not defined at $r = 1$. While the factor involving $I_i(r)$ is responsible for the oscillation of the function $f_m^0(r)$, the factor depending on $I_r(r)$ determines its growth as $r \rightarrow 1$: we see that $|f_m^0(r)|$ becomes unbounded in this limit when ${\rm Re}(\lambda ^0) < 4$ and approaches zero otherwise. The real and imaginary parts of $f_m^0(r)$ obtained for different eigenvalues $\lambda ^0$ are shown in figure 3, where it is evident that both the unbounded growth and the oscillations of $f_m^0(r)$ are localized in the neighbourhood of the endpoint $r = 1$. Given the singular nature of the solutions obtained at the leading order, the correction term $f_m^1(r)$ is rather difficult to compute and we do not attempt this here. If $f_{m-1}^1 \neq f_{m+1}^1$, the solution of (3.6) will also include some extra terms in addition to (3.7)–(3.8) which would correspond to another possible family of approximate eigenfunctions. However, as is evident from the discussion below, the solutions given in (3.7)–(3.8) capture the relevant behaviour. Finally, in view of ansatz (3.1), the leading-order approximations to eigenfunctions are obtained multiplying the function $f_m^0(r)$ by $\cos (m \theta )$ or $\sin (m \theta )$ with $m \rightarrow \infty$ which introduces rapid oscillations in the azimuthal direction.

Figure 3. Radial dependence (a) of the eigenvectors $f_m^0(r)$ associated with real eigenvalues $\lambda ^0 = 2$ (red solid line) and $\lambda ^0 = 6$ (blue dashed line) and (b) of the real part (red solid line) and the imaginary part (blue dashed line) of the eigenvector $f_m^0(r)$ associated with complex eigenvalue $\lambda ^0 = 3+10{\rm i}$. Panel (b) shows the neighbourhood of the endpoint $r=1$.

We thus conclude that when ${\rm Re}(\lambda ^0) < 4$, the solutions of eigenvalue problem (2.9) constructed in the form (3.1) include functions dominated by short-wavelength oscillations whose asymptotic, as $m \rightarrow \infty$, structure involves oscillations in both the radial and azimuthal directions and are localized near the boundary $\partial A_0$. Since as a result their pointwise values on $\partial A_0$ are not well defined, these solutions should be regarded as ‘distributions’. We remark that the asymptotic solutions constructed above do not satisfy the boundary conditions (3.2c)–(3.2d), which is consistent with the fact that they represent approximate eigenfunctions associated with the essential spectrum $\varPi _{ess}({\mathcal {H}})$ of the 2-D linearized Euler operator. In order to find solutions of eigenvalue problem (2.9) which do satisfy all the boundary conditions, we have to solve this problem numerically which is done next.

4. Numerical approaches

In this section we first describe the numerical approximation of eigenvalue problem (2.9)–(2.10) and then the time integration of the 2-D Euler system (1.1)–(1.2) with the initial condition in the form of the Lamb–Chaplygin dipole perturbed with some approximate eigenfunctions obtained by solving eigenvalue problem (2.9)–(2.10). These computations offer insights about the instability of the dipole complementary to the results of the asymptotic analysis presented in § 3.

4.1. Discretization of eigenvalue problem (2.9)–(2.10)

Eigenvalue problem (2.9)–(2.10) is solved using the spectral collocation method proposed by Fornberg (Reference Fornberg1996), see also the discussion in Trefethen (Reference Trefethen2000), which is based on a tensor grid in $(r,\theta )$. The discretization in $\theta$ involves trigonometric (Fourier) interpolation, whereas that in $r$ is based on Chebyshev interpolation where we take $r \in [-1,1]$ which allows us to avoid collocating (2.9a) at the origin when the number of grid points is even. Since then the mapping between $(r,\theta )$ and $(x_1,x_2)$ is 2-to-1, the solution must be constrained to satisfy the condition

(4.1)\begin{equation} {\tilde{\psi}}(r,\theta) = {\tilde{\psi}}({-}r,(\theta+{\rm \pi})(\text{mod} \ 2{\rm \pi})), \quad r \in [{-}1,1], \ \theta \in [0,2{\rm \pi}], \end{equation}

which is fairly straightforward to implement (Trefethen Reference Trefethen2000).

In contrast to (2.9a), the boundary condition (2.9b) does need to be evaluated at the origin which necessities modification of the differentiation matrix (since our Chebyshev grid does not include a grid point at the origin). The numbers of grid points discretizing the coordinates $r \in [-1,1]$ and $\theta \in [0,2{\rm \pi} ]$ are linked and both given by $N$ which is an even integer. The resulting algebraic eigenvalue problem then has the form

(4.2)\begin{equation} \lambda \boldsymbol{\psi} = \boldsymbol{H} \boldsymbol{\psi}, \end{equation}

where $\boldsymbol {\psi } \in {\mathbb {C}}^{N^2}$ is the vector of approximate nodal values of the eigenfunction and $\boldsymbol {H} \in {\mathbb {R}}^{N^2 \times N^2}$ the matrix discretizing the operator ${\mathcal {H}}$, cf. (2.9a), obtained as described above. Problem (4.2) is implemented in MATLAB and solved using the function eig. The discretization of all operators in ${\mathcal {H}}$, cf. (2.9), was carefully verified by applying them to analytic expressions and then comparing the results against exact expressions. Expected rates of convergence were observed as the resolution $N$ was increased.

Since the operator ${\mathcal {H}}$ and hence also the matrix $\boldsymbol {H}$ are non-normal and singular, the numerical conditioning of problem (4.2) may be poor, especially when the resolution $N$ is refined. In an attempt to mitigate this potential difficulty, we eliminated a part of the null space of $\boldsymbol {H}$ by performing projections on a certain number $N_C$ of eigenfunctions associated with the eigenvalue $\lambda = 0$ (they are obtained by solving problem (2.13) with different source terms $\phi _C$, $C = 2,3,\ldots,N_C+1$; cf. (B4)). However, solutions of problem (4.2) obtained in this way were essentially unchanged as compared with the original version. Moreover, in addition to examining the behaviour of the results when the grid is refined (by increasing the resolution $N$ as discussed in § 5), we have also checked the effect of arithmetic precision using the toolbox Advanpix (Reference Advanpix2017). Increasing the arithmetic precision up to ${\mathcal {O}}(10^2)$ significant digits was also not found to have a noticeable effect on the results obtained with small and medium resolutions $N \le 100$ (at higher resolutions the cost of such computations becomes prohibitive). These observations allow us to conclude that the results presented in § 5 are not affected by round-off errors.

In the light of the discussion in §§ 2.1 and 2.2, we know that the spectrum of the operator ${\mathcal {H}}$ includes essential spectrum in the form of a vertical band in the complex plane $|{\rm Re}(z)| \le 4$, $z \in {\mathbb {C}}$. Available literature on the topic of numerical approximation of infinite-dimensional non-self-adjoint eigenvalue problems, especially ones featuring essential spectrum, is very scarce. However, since the discretized problem (4.2) is finite-dimensional and therefore can only have a point spectrum, it is expected that at least some of the eigenvalues of the discrete problem will be approximations of the approximate eigenvalues in the essential spectrum $\varPi _{ess}({\mathcal {H}})$, whereas the corresponding eigenvectors will approximate the approximate eigenfunctions (we note that the term ‘approximate’ is used here with two distinct meanings: its first appearance refers to the numerical approximation and the second to the fact that these functions are defined as only ‘close’ to being true eigenfunctions, cf. § 2.1). As suggested by the asymptotic analysis presented in § 3, these approximate eigenfunctions are expected to be dominated by short-wavelength oscillations which cannot be properly resolved using any finite resolution $N$. Thus, since these eigenfunctions are not smooth, we do not expect our numerical approach to yield an exponential convergence of the approximation error. To better understand the properties of these eigenfunctions, we also solve a regularized version of problem (2.9) in which ${\tilde {\psi }}$ is replaced with ${\tilde {\psi }}_{\delta } := {\mathcal {R}}_{\delta }^{-1} {\tilde {\psi }}$, where ${\mathcal {R}}_{\delta } := (\operatorname {Id} - \delta ^2 \varDelta )$, $\delta > 0$ is a regularization parameter and the inverse of ${\mathcal {R}}_\delta$ is defined with the homogeneous Neumann boundary conditions. The regularized version of the discrete problem (4.2) then takes the form

(4.3)\begin{equation} \lambda_{\delta} \boldsymbol{\psi} = \boldsymbol{R}_{\delta} \boldsymbol{H} \boldsymbol{R}_{\delta}^{{-}1} \boldsymbol{\psi} =: \boldsymbol{H}_{\delta} \boldsymbol{\psi}, \end{equation}

where the subscript $\delta$ denotes regularized quantities and $\boldsymbol {R}_{\delta }$ is the discretization of the regularizing operator ${\mathcal {R}}_{\delta }$. Since the operator ${\mathcal {R}}_{\delta }^{-1}$ can be interpreted as a low-pass filter with the cut-off length given by $\delta$, the effect of this regularization is to smoothen the eigenvectors by filtering out components with wavelengths less than $\delta$. Clearly, in the limit when $\delta \rightarrow 0$ the original problem (4.2) is recovered. An analogous strategy was successfully employed by Protas & Elcrat (Reference Protas and Elcrat2016) in their study of the stability of Hill's vortex where the eigenfunctions also turned out to be singular distributions.

4.2. Solution of the time-dependent problem (1.1)–(1.2)

The 2-D Euler system (1.1)–(1.2) is transformed to the frame of reference moving with velocity $-U \boldsymbol {e}_1$ and rewritten in terms of the ‘perturbation’ vorticity ${\bar {\omega }}(t,\boldsymbol {x}) := \omega (t,\boldsymbol {x}) - \omega _0(\boldsymbol {x})$ and the corresponding perturbation streamfunction ${\bar {\psi }}(t,\boldsymbol {x})$, such that it takes the following form, cf. (2.1):

(4.4a)$$\begin{align} &\frac{\partial{{\bar{\omega}}}}{\partial{t}} + (\boldsymbol{\nabla}^\perp {\bar{\psi}} - U \boldsymbol{e}_1 ) \boldsymbol{\cdot} \boldsymbol{\nabla}\omega_0 + (\boldsymbol{\nabla}^\perp \psi_0 + \boldsymbol{\nabla}^\perp {\bar{\psi}} - U \boldsymbol{e}_1 ) \boldsymbol{\cdot} \boldsymbol{\nabla}{\bar{\omega}}= 0 \quad \text{in} \ \varOmega, \end{align}$$

(4.4b)$$\begin{align}- &\Delta {\bar{\psi}} = {\bar{\omega}} \quad \text{in} \ \varOmega, \end{align}$$

(4.4c)$$\begin{align}& {\bar{\psi}} \rightarrow 0 \quad \text{for} \ |\boldsymbol{x}| \rightarrow \infty. \end{align}$$

To facilitate solution of this system with a Fourier pseudospectral method (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988), we approximate the unbounded domain with a 2-D periodic box $\varOmega \approx {\mathbb {T}}^2 := [-L/2,L/2]^2$, where $L > 1$ is its size. While this is an approximation only, it is known to become more accurate as the size $L$ of the domain increases relative to the radius of the dipole which remains fixed at one (Boyd Reference Boyd2001). We note that this is a standard approach and has been successfully used in earlier studies of related problems (Nielsen & Rasmussen Reference Nielsen and Rasmussen1997; van Geffen & van Heijst Reference van Geffen and van Heijst1998; Billant et al. Reference Billant, Brancher and Chomaz1999; Donnadieu et al. Reference Donnadieu, Ortiz, Chomaz and Billant2009; Brion et al. Reference Brion, Sipp and Jacquin2014; Jugier et al. Reference Jugier, Fontane, Joly and Brancher2020). Since the instability has the form of short-wavelength oscillations localized on the dipole boundary $A_0$, interaction of the perturbed dipole with its periodic images does not have a significant effect.

The perturbation vorticity is then approximated in terms of a truncated Fourier series:

(4.5)\begin{equation} {\bar{\omega}}(t,\boldsymbol{x}) \approx \sum_{\boldsymbol{k} \in V_M} \hat{{\bar{\omega}}}_{\boldsymbol{k}}(t) \, {\rm e}^{{\rm i} \boldsymbol{k} \ \cdot \ \boldsymbol{x}}, \end{equation}

in which $\hat {{\bar {\omega }}}_{\boldsymbol {k}}(t) \in {\mathbb {C}}$ are the Fourier coefficients such that $\hat {{\bar {\omega }}}_{\boldsymbol {k}}(t) = \hat {{\bar {\omega }}}^*_{-\boldsymbol {k}}(t)$ and $V_M := \{ \boldsymbol {k} = [k_1,k_2] \in {\mathbb {Z}}^2 \, : \, -M/2 \le k_1, k_2 \le M/2 \}$, where $M$ is the number of grid points in each direction in ${\mathbb {T}}^2$. Substitution of expansion (4.5) into (4.4a) yields a system of coupled ordinary differential equations describing the evolution of the expansion coefficients $\hat {{\bar {\omega }}}_{\boldsymbol {k}}(t)$, $\boldsymbol {k} \in V_M$, which is integrated in time using the RK4 method. Product terms in the discretized equations are evaluated in physical space with the exponential filter proposed by Hou & Li (Reference Hou and Li2007) used in lieu of dealiasing. We use a massively parallel implementation based on MPI with Fourier transforms computed using the FFTW library (Frigo & Johnson Reference Frigo and Johnson2003). Convergence of the results with refinement of the resolution $M$ and of the time step $\Delta t$ as well as with the increase of the size $L$ of the computational domain was carefully checked. In the results reported in § 6 we use $L = 2{\rm \pi}$.

5. Solution of the eigenvalue problem

In this section we describe solutions of the discrete eigenvalue problem (4.2) and its regularized version (4.3). We mainly focus on the symmetric dipole with $\eta = 0$; cf. figure 1a. In order to study dependence of the solutions on the numerical resolution, problems (4.2)–(4.3) were solved with $N$ ranging from 20 to 260, where the largest resolution was limited by the amount of RAM available on a single node of the computer cluster to which we had access. The discrete spectra of problem (4.2) obtained with $N = 40$, 80, 160, 260 are shown in figure 4(a–d). We see that for all resolutions $N$ the spectrum consists of purely imaginary eigenvalues densely packed on the vertical axis and a ‘cloud’ of complex eigenvalues clustered around the origin (for each $N$ is there is also a pair of purely real spurious eigenvalues increasing as $|\lambda | = {\mathcal {O}}(N)$ when the resolution is refined; they are not shown in figure 4(a–d). We see that as $N$ increases the cloud formed by the complex eigenvalues remains restricted to the band $-2 \lessapprox {\rm Re}(\lambda ) \lessapprox 2$, but expands in the vertical (imaginary) direction. The spectrum is symmetric with respect to the imaginary axis as is expected for a Hamiltonian system. The eigenvalues fill the inner part of the band ever more densely as $N$ increases and in order to quantify this effect, in figure 5(a–d) we show the eigenvalue density defined as the number of eigenvalues in a small rectangular region of the complex plane, i.e.

(5.1)\begin{align} \mu(z) := \frac{\text{number of eigenvalues } \lambda \in \{ \zeta \in {\mathbb{C}} \, : \, |{\rm Re}(\zeta -z)| \le \Delta\lambda_r, \ |{\rm Im}(\zeta -z)| \le \Delta\lambda_i \}} {4 \Delta\lambda_r \Delta\lambda_i}, \end{align}

where $\Delta \lambda _r, \Delta \lambda _i \in {\mathbb {R}}$ are half-sizes of a cell used to count the eigenvalues with $\Delta \lambda _i \approx 500 \Delta \lambda _r$ reflecting the fact that the plots are stretched in the vertical direction. We see that as the resolution $N$ is refined the eigenvalue density $\mu (z)$ increases near the origin.

Figure 4. Eigenvalues obtained by solving the discrete eigenvalue problem (4.2) with different indicated resolutions $N$: (a) $N=40$, (b) $N = 80$, (c) $N = 160$, (d) $N= 260$. The eigenvalues $\pm \lambda _0$ and $\pm \lambda _0^*$ which converge to well-defined limits as the resolution $N$ is refined, cf. Table 1, are marked in red. The eigenvalue $\lambda _0$ is associated with the only linearly unstable mode, cf. § 6.

Figure 5. Eigenvalue densities (5.1) corresponding to the spectra shown in figure 4(a–d): (a) $N=40$, (b) $N=80$, (c) $N=160$, (d) $N=260$.

Table 1. Eigenvalue $\lambda _0$ associated with the linearly growing mode, cf. § 6, obtained by solving the discrete eigenvalue problem (4.2) with different resolutions $N$.

As discussed in § 2.1, a key question concerning the linear stability of 2-D Euler flows is the existence of point spectrum $\varPi _0({\mathcal {L}})$ of the linear operator ${\mathcal {L}}$, cf. (2.1). However, the usual approach based on discretizing the continuous eigenvalue problem (2.9)–(2.10), cf. § 4.1, is unable to directly distinguish numerical approximations of the true eigenvalues from those of the approximate eigenvalues. This can be done indirectly by solving the discrete problem (4.2) with different resolutions $N$ since approximations to true eigenvalues will then converge to well-defined limits as the resolution is refined; in contrast, approximations to approximate eigenvalues will simply fill up the essential spectrum $\varPi _{ess}({\mathcal {L}})$ ever more densely in this limit. In this way we have found a single eigenvalue, denoted $\lambda _0$, which together with its negative $-\lambda _0$ and complex conjugates $\pm \lambda _0^*$, satisfy the above condition (see table 1). As is evident from table 1, the differences between the real parts of $\lambda _0$ computed with different resolutions $N$ are very small and just over 1 %, although the variation of the imaginary part is larger. Moreover, as is discussed in § 6, $\lambda _0$ is in fact the only eigenvalue associated with a linearly growing mode.

We now take a closer look at the purely imaginary eigenvalues which are plotted for different resolutions $N$ in figure 6. It is known that these approximate eigenvalues are related to the periods of Lagrangian orbits associated with closed streamlines in the base flow (Cox Reference Cox2014). In particular, if the maximum period is bounded $\tau _{max} < \infty$, this implies the presence of horizontal gaps in the essential spectrum. However, as shown in Appendix C, the Lamb–Chaplygin dipole does involve Lagrangian orbits with arbitrarily long periods, such that the essential spectrum $\varPi _{ess}({\mathcal {L}})$ includes the entire imaginary axis $i{\mathbb {R}}$. The results shown in figure 6 are consistent with this property since the gap evident in the spectra shrinks, albeit very slowly, as the numerical resolution $N$ is refined. The reason why these gaps are present is that the orbits sampled with the discretization described in § 4.1 have only finite maximum periods which, however, become longer as the discretization is refined.

Figure 6. Purely imaginary eigenvalues obtained by solving the discrete eigenvalue problem (4.2) with different indicated resolutions $N$.

Finally, we analyse eigenvectors of problem (4.2) and choose to present them in terms of vorticity, i.e. we show ${\tilde {\omega }}_i = - \Delta {\tilde {\psi }}_i$, where the subscript $i=0,1,2$ enumerates the corresponding eigenvalues. First, in figure 7(a–c) we illustrate the convergence pattern of the eigenvector $\tilde {\omega }_0$ corresponding to the eigenvalue $\lambda _0$, cf. table 1, and representing an exponentially growing mode as the resolution is refined. We see that as $N$ increases the approximations of the eigenvector converge to a constant value within the domain $A_0$ and diverge near its boundary, in agreement with the distributional nature of these eigenvectors established by our asymptotic analysis in § 3. More specifically, we see that the magnitude $|{\tilde {\omega }}_1(r,\theta )|$ of the eigenvector grows rapidly near the boundary, i.e. as $r \rightarrow 1$, which is consistent with the behaviour of the function $f_m^0(r)$ describing the asymptotic solution, cf. expressions (3.7)–(3.9) and figure 3(a,b). In addition, a rapid variation of $|{\tilde {\omega }}_1(r,\theta )|$ in the azimuthal coordinate $\theta$ is also evident for $r \lesssim 1$. However, something that could not be discerned by the asymptotic analysis is that these oscillations are mostly concentrated near the azimuthal angles $\theta = \pm {\rm \pi}/ 4, \pm 3{\rm \pi} /4$. As expected, both the growth in the radial direction and the oscillations in the azimuthal direction become more rapid as the resolution $N$ is refined. Given the distributional nature of the solutions of the eigenvalue problem (2.9), the classical notion of ‘convergence’ of a numerical scheme is not entirely applicable here and instead one would need to refer to more refined concepts such as ‘weak convergence’, but since they are quite technical, we do not pursue this avenue here. However, as is shown in § 6, even if they are not fully resolved, the eigenvectors computed here still contain useful information.

Figure 7. Real parts of the eigenvector $\tilde {\omega }_0$ corresponding to the eigenvalue $\lambda _0$, cf. table 1, and representing an exponentially growing mode obtained by solving the discrete eigenvalue problem (4.2) with different resolutions $N$: (a) $N = 40$, (b) $N = 80$ and (c) $N = 120$. The grids covering the surface plots represent the discretizations of the domain $A_0$ used for different $N$.

Next, in figure 8(a,c,e) we compare the real parts of the eigenvectors associated with different eigenvalues: the complex eigenvalue $\lambda _0$ corresponding to the exponentially growing mode (already shown in figure 7), a purely real eigenvalue $\lambda _1$ and a purely imaginary eigenvalue $\lambda _2$. We see that while these eigenvectors are qualitatively similar and share the features described above, the eigenvector ${\tilde {\omega }}_0$ is symmetric with respect to the flow centreline, whereas the eigenvectors ${\tilde {\omega }}_1$ and ${\tilde {\omega }}_2$ are antisymmetric. Another difference is that in the eigenvector ${\tilde {\omega }}_0$ associated with the eigenvalue $\lambda _0$ the oscillations are mostly concentrated near the azimuthal angles $\theta = \pm {\rm \pi}/ 4, \pm 3{\rm \pi} /4$, cf. figure 8(a); on the other hand, in the eigenvectors ${\tilde {\omega }}_1$ and ${\tilde {\omega }}_2$ the oscillations are mostly concentrated near the stagnation points $\boldsymbol {x}_a$ and $\boldsymbol {x}_b$, cf. figure 8(c,e).

Figure 8. Real parts of the eigenvectors corresponding to the indicated eigenvalues obtained by solving (a,c,e) eigenvalue problem (4.2) and (b,d, f) the regularized problem (4.3) using the resolution $N=80$: (a) $\lambda _0=0.126 \pm i25.258$, (b) $\lambda _{\delta,0}=0.129 \pm i67.489$, $\delta =0.05$, (c) $\lambda _1=1.585$, (d) $\lambda _{\delta,1}=0.406$, $\delta =0.05$, (e) $\lambda _2=\pm i149.873$ and ( f) $\lambda _{\delta,2}=\pm i150.233$, $\delta =0.05$. The grid shown on the surface represents the discretization of the domain $A_0$ used in the numerical solution of problems (4.2) and (4.3).

The numerical approximations of the eigenvectors are characterized by short-wavelength oscillations. Here, ‘short wavelength’ means that a significant variation of the magnitude $|{\tilde {\omega }}(r,\theta )|$ of the eigenvector with respect to both $r$ and $\theta$ occurs on the length scale given by the grid size which shrinks as the resolution is refined. This feature is also borne out in figure 10 showing the enstrophy spectrum of the initial condition involving the eigenvector ${\tilde {\omega }}_0$. It is evident from this figure that significant contributions to the enstrophy come from a broad range of length scales, including the smallest length scales resolved on the numerical grid. The eigenvectors associated with all other eigenvalues (not shown here for brevity) are also dominated by short-wavelength oscillations localized near different parts of the boundary $\partial A_0$. Since due to their highly oscillatory nature the eigenvectors shown in figure 8(a,c,e) are not fully resolved, in figure 8(b,d, f) we show the corresponding eigenvectors of the regularized eigenvalue problem (4.3) where we set $\delta = 0.05$. We see that in the regularized eigenvectors oscillations are shifted to the interior of the domain $A_0$ and their typical wavelengths are much larger. The eigenvalues obtained by solving the regularized problem (4.3) are distributed following a similar pattern as revealed by the eigenvalues of the original problem (4.2), cf. figures 4(b) and 5(b). In particular, for the main eigenvalue of interest $\lambda _0$, cf. table 1, the difference with respect to the corresponding eigenvalue of the regularized problem $\lambda _{\delta,0}$ is rather small and we have $|{\rm Re}(\lambda _{0} - \lambda _{\delta,0})| / {\rm Re}(\lambda _{0}) \approx 0.024$ (both are computed here with the resolution $N = 80$). The remaining eigenvalues of the regularized problem also form a ‘cloud’ filling the essential spectrum $\varPi _{ess}({\mathcal {L}})$.

Solution of the discrete eigenvalue problem (4.2) for asymmetric dipoles with $\eta > 0$ leads to eigenvalue spectra and eigenvectors qualitatively very similar to those shown in figures 4(a–d) and 8(a,c,e), hence for brevity they are not shown here. The only noticeable difference is that the eigenvectors are no longer symmetric or antisymmetric with respect to the flow centreline.

6. Solution of the evolution problem

As in § 5, we focus on the symmetric case with $\eta = 0$. The 2-D Euler system (1.1)–(1.2) is solved numerically as described in § 4.2 with the initial condition for the perturbation vorticity ${\bar {\omega }}(t,\boldsymbol {x})$ given in terms of the eigenvectors shown in figure 8(a–f), i.e.

(6.1)\begin{equation} {\bar{\omega}}(0,\boldsymbol{x}) = \varepsilon \frac{\| {\tilde{\omega}}_i\|_{L^2(\varOmega)}} {\| \omega_0 \|_{L^2(\varOmega)}} {\tilde{\omega}}_i(\boldsymbol{x}) \quad \text{or} \quad {\bar{\omega}}(0,\boldsymbol{x}) = \varepsilon \frac{\| {\tilde{\omega}}_{\delta,i}\|_{L^2(\varOmega)}} {\| \omega_0 \|_{L^2(\varOmega)}} {\tilde{\omega}}_{\delta,i}(\boldsymbol{x}), \quad i=0,1,2. \end{equation}

Unless indicated otherwise, the numerical resolution is $M = 512$ grid points in each direction. By taking $\varepsilon = 10^{-4}$ we ensure that the evolution of the perturbation vorticity is effectively linear up to $t \lesssim 70$ and to characterize its growth we define the perturbation enstrophy as

(6.2)\begin{equation} {\mathcal{E}}(t) := \int_{\varOmega} {\bar{\omega}}(t,\boldsymbol{x})^2 \, {\rm d}\kern0.7pt \boldsymbol{x}. \end{equation}

The evolution of this quantity is shown in figure 9(a) for the six considered initial conditions and times before nonlinear effects become evident. In all cases we see that after a transient period the perturbation enstrophy starts to grow exponentially as $\exp ({\tilde {\lambda }} t )$, where the growth rate (found via a least-squares fit) is ${\tilde {\lambda }} \approx 0.127$ and is essentially equal to the real part of the eigenvalue $\lambda _0$, cf. table 1. The duration of the transient, which involves an initial decrease of the perturbation enstrophy, is different in different cases and is shortest when the eigenfunctions ${\tilde {\omega }}_0$ and ${\tilde {\omega }}_{\delta,0}$ are used as the initial conditions in (6.1) (in fact, in the latter case the transient is barely present). The reason for this behaviour is that ${\tilde {\omega }}_0$ is the sole true eigenvector of the operator ${\mathcal {L}}$, whereas ${\tilde {\omega }}_1$ and ${\tilde {\omega }}_2$ are only approximate eigenvectors associated with the (approximate) eigenvalues $\lambda _1$ and $\lambda _2$ belonging to the essential spectrum $\varPi _{ess}({\mathcal {L}})$ rather than to the point spectrum $\varPi _0({\mathcal {L}})$. As a result, ${\tilde {\omega }}_0$ represents the only linearly growing mode, such that when ${\tilde {\omega }}_1$, ${\tilde {\omega }}_2$ or any other approximate eigenvector is used as the initial condition in (6.1), a transient behaviour ensues where the solution ${\bar {\omega }}(t)$ of system (4.4) approaches the trajectory involving the growing mode ${\rm Re}({\rm e}^{\lambda _0 t} {\tilde {\omega }}_0 )$. Hereafter we focus on the flow obtained with the initial condition (6.1) given in terms of the eigenfunction ${\tilde {\omega }}_0$, cf. figure 8(a).

Figure 9. Time evolution of the normalized perturbation enstrophy ${\mathcal {E}}(t) / {\mathcal {E}}(0)$ in the flow with the initial condition (6.1) given in terms of (a) the different eigenvectors shown in figure 8(a–f) and obtained with a fixed resolution $N = 80$ and (b) the eigenvector ${\tilde {\omega }}_0$ computed with different resolutions $N$. In (a) the red solid lines correspond to the eigenvectors ${\tilde {\omega }}_0$ and ${\tilde {\omega }}_{\delta,0}$, black dotted lines to ${\tilde {\omega }}_1$ and ${\tilde {\omega }}_{\delta,1}$ and blue dashed lines to ${\tilde {\omega }}_2$ and ${\tilde {\omega }}_{\delta,2}$; thin and thick lines represent flows with initial conditions involving eigenvectors obtained as solutions of the discrete eigenvalue problem (4.2) and its regularized version (4.3), respectively. In (b) the blue dashed and red solid lines correspond to initial conditions involving eigenvectors ${\tilde {\omega }}_0$ obtained with the resolutions $N=40$ and $N=80$, respectively.

The effect of the numerical resolution $N$ used in the discrete eigenvalue problem (4.2) is analysed in figure 9(b), where we show the perturbation enstrophy (6.2) in the flows with the eigenvector ${\tilde {\omega }}_0$ used in the initial conditions (6.1) computed with different $N$. We see that refined resolution leads to a longer transient period while the rate of the exponential growth ${\tilde {\lambda }}$ is unchanged. This demonstrates that this growth rate is in fact a robust property unaffected by the underresolution of the unstable mode.

The enstrophy spectrum of the initial condition (6.1) and of the perturbation vorticity ${\bar {\omega }}(t,\boldsymbol {x})$ at different times $t \in (0,60]$ is shown in figure 10 as a function of the wavenumber $k := |\boldsymbol {k}|$. It is defined as

(6.3)\begin{equation} e(t,k) := \int_{S_{k}} | \hat{{\bar{\omega}}}_{\boldsymbol{k}}(t)|^2 \, {\rm d}\sigma, \end{equation}

where $\sigma$ is the azimuthal angle in wavenumber space and $S_{k}$ denotes the circle of radius $k$ in this space (with some abuse of notation justified by simplicity; here we have treated the wavevector $\boldsymbol {k}$ as a continuous rather than discrete variable). Since its enstrophy spectrum is essentially independent of the wavenumber $k$, the eigenvector ${\tilde {\omega }}_0$ in the initial condition (6.1) turns out to be a distribution rather than a smooth function. The enstrophy spectra of the perturbation vorticity ${\bar {\omega }}(t,\boldsymbol {x})$ during the flow evolution show a rapid decay at high wavenumbers which is the effect of the applied filter, cf. § 4.2. However, after the transient, i.e. for $20 \lessapprox t \leq 60$, the enstrophy spectra have very similar forms, except for a vertical shift which increases with time $t$. This confirms that the time evolution is dominated by linear effects as there is little energy transfer to higher (unresolved) modes. This is also attested to by the fact that for all the cases considered in figure 9(a) the relative change of the total energy $\int _{\varOmega } |\boldsymbol {u}(t,\boldsymbol {x})|^2\, {\rm d}\boldsymbol {x}$ and of the total enstrophy $\int _{\varOmega } \omega (t,\boldsymbol {x}) ^2\, {\rm d}\boldsymbol {x}$, which are conserved quantities in the Euler system (1.1)–(1.2), is at most of ${\mathcal {O}}(10^{-4})$ (this small variation of the conserved quantities is due to the action of the filter and the fact that the time-integration scheme is not strictly conservative, cf. § 4.2). Since in the numerical solution the total circulation is given by the Fourier coefficient $[\hat {\omega }(t)]_{\boldsymbol {0}} = [\hat {\omega }_0(t)]_{\boldsymbol {0}} + [\hat {\omega }_i(t)]_{\boldsymbol {0}}$, it remains zero by construction throughout the entire flow evolution.

Figure 10. Enstrophy spectra (6.3) of (blue squares) the initial condition (6.1) involving the eigenvector ${\tilde {\omega }}_0$ and (red circles) the corresponding perturbation vorticity ${\bar {\omega }}(t,\boldsymbol {x})$ at times $t = 10,20,\ldots,60$. The arrow indicates the trend with the increase of time $t$.

We now go on to discuss the time evolution of the perturbation vorticity in physical space and in figures 11(a) and 11(b) we show ${\bar {\omega }}(t,\boldsymbol {x})$ at the times $t = 4$ and $t = 21$, respectively, which correspond to the transient regime and to the subsequent period of an exponential growth. During that period, i.e. for $20 \lessapprox t \leq 60$, the structure of the perturbation vorticity field does not change much. We see that as the perturbation evolves a number of thin vorticity filaments is ejected from the vortex core $A_0$ into the potential flow with the principal ones emerging at the azimuthal angles $\theta \approx \pm {\rm \pi}/ 4, \pm 3{\rm \pi} / 4$, i.e. in the regions of the vortex boundary where most of the short-wavelength oscillations evident in the eigenvector ${\tilde {\omega }}_0$ are localized, cf. figure 8(a). With thickness of the order of a few grid points, these filaments are among the finest structures that can be resolved in computations with the resolution $M$ we use. The perturbation remains symmetric with respect to the flow centreline for all times and since the vorticity $\omega _0$ of the base flow is antisymmetric, the resulting total flow $\omega (t,\boldsymbol {x})$ does not possess any symmetries. The perturbation vorticity ${\bar {\omega }}(t,\boldsymbol {x})$ realizing the exponential growth in the flows corresponding to the initial condition involving the eigenvectors ${\tilde {\omega }}_1$ and ${\tilde {\omega }}_2$ (and their regularized versions ${\tilde {\omega }}_{\delta,1}$ and ${\tilde {\omega }}_{\delta,2}$) is essentially identical to the perturbation vorticity shown in figure 11(b), although its form during the transient regime can be quite different. In particular, the perturbation eventually becomes symmetric with respect to the flow centreline even if the initial condition (6.1) is antisymmetric. The same is true for flows obtained with initial condition corresponding to all approximate eigenvalues other than $\lambda _1$ and $\lambda _2$ (not shown here for brevity). We did not attempt to study the time evolution of asymmetric dipoles with $\eta > 0$ in (1.5a), since their vorticity distributions are discontinuous making computation of such flows using the pseudospectral method described in § 4.2 problematic.

Figure 11. Perturbation vorticity ${\bar {\omega }}(t,\boldsymbol {x})$ in the flow corresponding the initial condition (6.1) involving the eigenvector ${\tilde {\omega }}_0$ during (a) the transient regime ($t= 4$) and (b) the period of exponential growth ($t = 21$).