On the Linear Stability of the Lamb-Chaplygin Dipole

The Lamb-Chaplygin dipole (Lamb1895,Lamb1906,Chaplygin1903) is one of the few closed-form relative equilibrium solutions of the 2D Euler equation characterized by a continuous vorticity distribution. We consider the problem of its linear stability with respect to 2D circulation-preserving perturbations. It is demonstrated that this flow is linearly unstable, although the nature of this instability is subtle and cannot be fully understood without accounting for infinite-dimensional aspects of the problem. To elucidate this, we first derive a convenient form of the linearized Euler equation defined within the vortex core which accounts for the potential flow outside the core while making it possible to track deformations of the vortical region. The linear stability of the flow is then determined by the spectrum of the corresponding operator. Asymptotic analysis of the associated eigenvalue problem shows the existence of approximate eigenfunctions in the form of short-wavelength oscillations localized near the boundary of the vortex and these findings are confirmed by the numerical solution of the eigenvalue problem. However, the time-integration of the 2D Euler system reveals the existence of only one linearly unstable eigenmode and since the corresponding eigenvalue is embedded in the essential spectrum of the operator, this unstable eigenmode is also shown to be a distribution characterized by short-wavelength oscillations rather than a smooth function. These findings are consistent with the general results known about the stability of equilibria in 2D Euler flows and have been verified by performing computations with different numerical resolutions and arithmetic precisions.


Introduction
The Lamb-Chaplygin dipole is a relative equilibrium solution of the two-dimensional (2D) Euler equations in an unbounded domain R 2 that was independently obtained by Lamb (1895Lamb ( , 1906) ) and Chaplygin (1903); the history of this problem was surveyed by Meleshko & van Heijst (1994).The importance of the Lamb-Chaplygin dipole stems from the fact that this is a simple exact solution with a continuous vorticity distribution which represents a steadily translating vortex pair (Leweke et al., 2016).Such objects are commonly used as models in geophysical fluid dynamics where they are referred to as "modons" (Flierl, 1987).Interestingly, despite the popularity of this model, the stability properties of the Lamb-Chaplygin dipole are still not well understood and the goal of the present investigation is to shed some new light on this question.
We consider an unbounded flow domain Ω := R 2 (":=" means "equal to by definition").Flows of incompressible inviscid fluids are described by the 2D Euler equation which can be written in the vorticity form as where t ∈ (0, T ] is the time with T > 0 denoting the length of the interval considered, ω : (0, T ] × Ω → R is the vorticity component perpendicular to the plane of motion and u = [u 1 , u 2 ] T : (0, T ] × Ω → R 2 is a divergence-free velocity field (i.e., ∇ • u = 0).
The space coordinate will be denoted x = [x 1 , x 2 ] T .Introducing the streamfunction ψ : (0, T ] × Ω → R, the relation between the velocity and vorticity can be expressed as T and ∆ψ = −ω. (2) System (1)-(2) needs to be complemented with suitable initial and boundary conditions, and they will be specified below.
In the frame of reference translating with the velocity −U e 1 , where U > 0 and e i , i = 1, 2, is the unit vector associated with the ith axis of the Cartesian coordinate system, equilibrium solutions of system (1)-( 2) satisfy the boundary-value problem (Wu et al., 2006) where the "vorticity function" F : R → R need not be continuous.Clearly, the form of the equilibrium solution is determined by the properties of the function F (ψ). Assuming without loss of generality that it has unit radius (a = 1), the Lamb-Chaplygin dipole is obtained by taking where b ≈ 3.8317059702075123156 is the first root of the Bessel function of the first kind of order one, J 1 (b) = 0, and η ∈ (−∞, ∞) is a parameter characterizing the asymmetry of the dipole (in the symmetric case η = 0).The solution of (3)-( 4) then has the form of a circular vortex core of unit radius embedded in a potential flow.The vorticity and streamfunction are given by the following expressions stated in the cylindrical coordinate system (r, θ) (hereafter we will adopt the convention that the subscript "0" refers to an equilibrium solution) • inside the vortex core (0 < r ≤ 1, 0 < θ ≤ 2π): • outside the vortex core (r > 1, 0 < θ ≤ 2π): The vortical core region will be denoted A 0 := {x ∈ R 2 : ∥x∥ ≤ 1} and ∂A 0 will denote its boundary.The streamline pattern inside A 0 in the symmetric (η = 0) and asymmetric (η > 0) case is shown in figures 1a and 1b, respectively.Various properties of the Lamb-Chaplygin dipole are discussed by Meleshko & van Heijst (1994).In particular, it is shown that regardless of the value of η the total circulation of the dipole vanishes, i.e., Γ 0 := A 0 ω 0 dA = 0. We note that in the limit η → ±∞ the dipole approaches a state consisting of a monopolar vortex with a vortex sheet of opposite sign coinciding with the part of the boundary ∂A 0 above or below the flow centerline, respectively, for positive and negative η. Generalizations of the Lamb-Chaplygin dipole corresponding to differentiable vorticity functions F (ψ) were obtained numerically by Albrecht et al. (2011), whereas multipolar generalizations were considered by Viúdez (2019b,a).
Most investigations of the stability of the Lamb-Chaplygin dipole were carried out in the context of viscous flows governed by the Navier-Stokes system, beginning with the computations of dipole evolution performed by Nielsen & Rasmussen (1997); van Geffen & van Heijst (1998).While relations (5)-( 6) do not represent an exact steady-state solution of the Navier-Stokes system, this approximate approach was justified by the assumption that viscous effects occur on time scales much longer than the time scales characterizing the growth of perturbations.A first such study of the stability of the dipole was conducted by Billant et al. (1999) who considered perturbations with dependence on the axial wavenumber and found several unstable eigenmodes together with their growth rates by directly integrating the three-dimensional (3D) linearized Navier-Stokes equations in time.Additional unstable eigenmodes were found in the 2D limit corresponding to small axial wavenumbers by Brion et al. (2014).The transient growth due to the non-normality of the linearized Navier-Stokes operator was investigated in the related case of a vortex pair consisting of two Lamb-Oseen vortices by Donnadieu et al. (2009) and Jugier et al. (2020), whereas Sipp & Jacquin (2003) studied Widnalltype instabilities of such vortex pairs.The effect of stratification on the evolution of a perturbed Lamb-Chaplygin dipole in 3D was considered by Waite & Smolarkiewicz (2008); Bovard & Waite (2016).The history of the studies concerning the stability of vortices in ideal fluids was recently surveyed by Gallay (2019).
The only stability analysis of the Lamb-Chaplygin dipole in the inviscid setting we are aware of is due to Luzzatto-Fegiz & Williamson (2012); Luzzatto-Fegiz (2014) who employed methods based on imperfect velocity-impulse diagrams applied to an approximation of the dipole in terms of a piecewise-constant vorticity distribution and concluded that this configuration is stable.Finally, there is a recent mathematically rigorous result by Abe & Choi (2022) who established orbital stability of the Lamb-Chaplygin dipole (orbital stability implies that flows corresponding to "small" perturbations of the dipole remain "close" in a certain norm to the translating dipole; hence, this is a rather weak notion of stability).
As noted by several authors (Meleshko & van Heijst, 1994;Waite & Smolarkiewicz, 2008;Luzzatto-Fegiz & Williamson, 2012;Abe & Choi, 2022), the stability properties of the Lamb-Chaplygin dipole are still to be fully understood despite the fact that it was introduced more than a century ago.To the best of our knowledge, the present study is the first comprehensive investigation of the linear stability of the Lamb-Chaplygin dipole in the inviscid case, which is the only setting where it represents a true equilibrium solution of the governing equations.As a result, we find behavior that was not observed in any of the earlier studies.It is demonstrated that the Lamb-Chaplygin dipole is in fact linearly unstable, but the nature of this instability is quite subtle and cannot be understood without referring to the infinite-dimensional nature of the linearized governing equations.More specifically, both the asymptotic and numerical solution of an eigenvalue problem for the 2D linearized Euler operator suitably localized to the vortex core A 0 confirm the existence of an essential spectrum with the corresponding approximate eigenfunctions in the form of short-wavelength oscillations localized near the vortex boundary ∂A 0 .However, the time-integration of the 2D Euler system reveals the presence of a single exponentially growing eigenmode and since the corresponding eigenvalue is embedded in the essential spectrum of the operator, this unstable eigenmode is also found not to be a smooth function and exhibits short-wavelength oscillations.These findings are consistent with the general mathematical results known about the stability of equilibria in 2D Euler flows (Shvydkoy & Latushkin, 2003;Shvydkoy & Friedlander, 2005) and have been verified by performing computations with different numerical resolutions and, in the case of the eigenvalue problem, with different arithmetic precisions.
The structure of the paper is as follows: in the next section we review some basic facts about the spectra of the 2D linearized Euler equation and transform this system to a form in which its spectrum can be conveniently studied with an asymptotic method and numerically; a number of interesting properties of the resulting eigenvalue problem is also discussed, an approximate asymptotic solution of this eigenvalue problem is constructed in § 3, the numerical approaches used to solve the eigenvalue problem and the initialvalue problem (1)-(2) are introduced in § 4, whereas the obtained computational results are presented in § 5 and § 6, respectively; discussion and final conclusions are deferred to § 7; some more technical material is collected in three appendices.

2D Linearized Euler Equations
The Euler system (1)-( 2) formulated in the moving frame of reference and linearized around an equilibrium solution {ψ 0 , ω 0 } has the following form, where ψ ′ , ω ′ : (0, T ] × Ω → R are the perturbation variables (also defined in the moving frame of reference) in which ∆ −1 is the inverse Laplacian corresponding to the far-field boundary condition (7c) and w ′ is an appropriate initial condition assumed to have zero circulation, i.e., Ω w ′ dA = 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 2D Euler flows it was proved by Vishik & Friedlander (2003); Lin (2004) 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., 2006) predicts that equilibria satisfying system (3) are nonlinearly stable if F ′ (ψ) ≥ 0, which however is not the case for the Lamb-Chaplygin dipole, since using (4) we have F ′ (ψ 0 ) = −b 2 < 0 for ψ 0 ≥ η.Thus, Arnold's criterion is inapplicable in this case.

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 du/dt = Au on a Banach space X (in general, infinite-dimensional) with the state u(t) ∈ X and a linear operator A : X → X .Solution of this problem can be formally written as u(t) = e At u 0 , where u 0 ∈ X is the initial condition and e At the semigroup generated by A (Curtain & Zwart, 2013).While in finite dimensions linear operators can be represented as matrices which can only have point spectrum Π 0 (A), in infinite dimensions the situation is more nuanced since the spectrum Λ(A) of the linear operator A may in general consist of two parts, namely, the approximate point spectrum Π(A) (which is a set of numbers λ ∈ C such that (A − λ) is not bounded from below) and the compression spectrum Ξ(A) (which is a set of numbers λ ∈ C such that the closure of the range of (A − λ) does not coincide with X ).We thus have Λ(A) = Π(A) ∪ Ξ(A) and the two types of spectra may overlap, i.e., Π(A) ∩ Ξ(A) ̸ = ∅ (Halmos, 1982).A number λ ∈ C belongs to the approximate point spectrum Π(A) if and only if there exists a sequence of unit vectors {f n }, referred to as approximate eigenvectors, such that ∥(A − λ)f n ∥ X → 0 as n → ∞.If for some λ ∈ Π(A) there exists a unit element f such that Af = λf , then λ and f are an eigenvalue and an eigenvector of A. The set of all eigenvalues λ forms the point spectrum Π 0 (A) which is contained in the approximate point spectrum, Π 0 (A) ⊂ Π(A).If λ ∈ Π(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 λ ∈ C is referred to as the essential spectrum Π ess (A) := Π(A)\Π 0 (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 e At one is usually interested in understanding the relation between its growth abscissa γ(A) := lim t→∞ t −1 ln ∥e At ∥ X and the spectrum Λ(A) of A. While in finite dimensions γ(A) is determined by the eigenvalues of A with the largest real part, in infinite dimensions the situation is more nuanced since there are examples in which sup z∈Λ(A) ℜ(z) < γ(A), e.g., Zabczyk's problem (Zabczyk, 1975) also discussed by Trefethen (1997); some problems in hydrodynamic stability where such behavior was identified are analyzed by Renardy (1994).
In regard to the 2D linearized Euler operator L, cf.(7a), it was shown by Shvydkoy & Latushkin (2003) 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 λ max in the flow field and to the index m ∈ Z of the Sobolev space H m (Ω) in which the evolution problem is formulated (i.e., X = H m (Ω) above).The norm in the Sobolev space H m (Ω) is defined as ∥u∥ , where α 1 , α 2 ∈ Z with |α| := α 1 + α 2 (Adams & Fournier, 2005).More specifically, we have (Shvydkoy & Friedlander, 2005) In 2D flows Lyapunov exponents are determined by the properties of the velocity gradient ∇u(x) at hyperbolic stagnation points x 0 .More precisely, λ max is given by the largest eigenvalue of ∇u(x) computed over all stagnation points.As regards the Lamb-Chaplygin dipole, it is evident from figures 1a and 1b that in both the symmetric and asymmetric case it has two stagnation points x a and x b located at the fore and aft extremities of the vortex core.Inspection of the velocity field ∇ ⊥ ψ 0 defined in (5a) shows that the largest eigenvalues of ∇u(x) evaluated at these stagnation points, and hence the Lyapunov exponents, are λ max = 2 regardless of the value of η.
While characterization of the essential spectrum of the 2D linearized Euler operator 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 (Drazin & Reid, 1981;Chandrasekhar, 1961) or the cellular cat's eyes flows (Friedlander et al., 2000).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 (2004) 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.

Linearization Around the Lamb-Chaplygin Dipole
The linear system (7) is defined on the entire plane R 2 , however, in the Lamb-Chaplygin dipole the vorticity ω 0 is supported within the vortex core A 0 only, cf.(6b).This will allow us to simplify system (7) 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 (7d) is also supported in A 0 , then the initial-value problem (7) can be regarded as a free-boundary problem describing the evolution of the boundary ∂A(t) of the vortex core (we have A(0) = A 0 and ∂A(t) = ∂A 0 ).However, as explained below, the evolution of this boundary can be deduced from the evolution of the perturbation streamfunction ψ ′ (t, x), hence need not be tracked independently.Thus, the present problem is different from, e.g., the vortexpatch problem where the vorticity distribution is fixed (piecewise constant in space) and in the stability analysis the boundary is explicitly perturbed (Elcrat & Protas, 2013). Denoting the perturbation streamfunction in the vortex core and in its complement, system (7) can be recast as where n is the unit vector normal to the boundary ∂A 0 pointing outside and conditions (9d)-(9e) represent the continuity of the normal and tangential perturbation velocity components across the boundary ∂A 0 with f ′ : ∂A 0 → R denoting the unknown value of the perturbation streamfunction at that boundary.
The velocity normal to the vortex boundary ∂A(t) is u n := u•n = ∂ψ 1 /∂s = ∂ψ 2 /∂s, where s is the arc-length coordinate along ∂A(t), cf.(9d).While this quantity identically vanishes in the equilibrium state ( 5)-( 6), cf ( 17), in general it will be nonzero resulting in a deformation of the boundary ∂A(t).This deformation can be deduced from the solution of system (9) as follows.Given a point z ∈ ∂A(t), the deformation of the boundary is described by dz/dt = n u n | ∂A(t) .Integrating this expression with respect to time yields where z(0) ∈ ∂A 0 and 0 < τ ≪ 1 is the time over which the deformation is considered.Thus, the normal deformation of the boundary can be defined as ρ(τ We also note that at the leading order the area of the vortex core A(t) is preserved by the considered perturbations We notice that in the exterior domain R 2 \A 0 the problem is governed by Laplace's equation (9c) subject to boundary conditions (9d)-(9f).Therefore, this subproblem can be eliminated by introducing the corresponding Dirichlet-to-Neumann (D2N) map which is constructed in an explicit form in Appendix A. Thus, equation (9c) with boundary conditions (9d)-( 9f) can be replaced with a single relation 1 holding on ∂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 ∂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 ∂A 0 using relation (10).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 (RHS) in (9a) evaluated on the boundary ∂A 0 .
Noting that the base state satisfies the equation 3)-( 4), and using the identity • ∇ ⊥ ψ 0 , the vorticity equation ( 9a) can be transformed to the following simpler form where we also used (9b) to eliminate ω ′ in favor of ψ ′ 1 .Supposing the existence of an eigenvalue λ ∈ C and an eigenfunction ψ : A 0 → C, we make the following ansatz for the perturbation streamfunction ψ ′ 1 (t, x) = ψ(x) e λt which leads to the eigenvalue problem ∂∆ ψ ∂r = 0, at r = 0, where ∆ −1 M is the inverse Laplacian subject to the boundary condition ∂ ψ/∂n − M ψ = 0 imposed on ∂A 0 and the additional boundary condition (13b) ensures the perturbation vorticity is differentiable at the origin (such condition is necessary since the differential operator on the RHS in (13a) is of order three).Depending on whether or not the different differential operators appearing in it are inverted, eigenvalue problem (13) 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 ( 13) is preferred and we will focus on it hereafter.
We note that the proposed formulation ensures that the eigenfunctions ψ have zero circulation, as required where we used the divergence theorem, equations (9b)-(9c) and the boundary conditions (9e)-(9f).Since it will be needed for the numerical discretization described in § 5, we now rewrite the eigenvalue problem (13) explicitly in the polar coordinate system where ∂θ 2 and the velocity components obtained as u r 0 , u θ They have the following behavior on the boundary ∂A 0 , where "∼" means the norms on the left and on the right are equivalent (in the precise sense of norm equivalence), the essential spectrum (8) of the operator H will have m = −2, so that Π ess (H) is a vertical band in the complex plane with Operator H, cf.(15a) has a non-trivial null space Ker(H).To see this, we consider the "outer" subproblem Then, the eigenfunctions ψ C spanning the null space of operator H are obtained as solutions of the family of "inner" subproblems where C = 2, 3, . . . .Some of these eigenfunctions are shown in figures 2a-d, where distinct patters are evident for even and odd values of C. -1

Asymptotic Solution of Eigenvalue Problem (15)
A number of interesting insights about certain properties of solutions of eigenvalue problem ( 15) 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 (η = 0) and begin by introducing the ansatz where f m , g m : [0, 1] → C, m = 1, 2, . . ., are functions to be determined.Substituting this ansatz in (15a) with the Laplacian moved back to the left-hand side (LHS) and applying well-known trigonometric identities leads after some algebra to the following system of coupled third-order ordinary differential equations (ODEs) for the functions where the Bessel operator B m is defined via r f , whereas the coefficient functions have the form, cf. ( 16), The functions g m (r), m = 1, 2, . . ., satisfy a system identical to (21), which shows that the eigenfunctions ψ(r, θ) are either even or odd functions of θ (i.e., they are either symmetric or antisymmetric with respect to the flow centerline).Moreover, the fact that system (21) couples Fourier components corresponding to different m implies that the eigenvectors ψ(r, θ) are not separable as functions of r and θ.Motivated by our discussion in § 2.1 about the properties of approximate eigenfunctions of the 2D linearized Euler operator, we will construct approximate solutions of system (21) in the short-wavelength limit m → ∞.In this analysis we will assume that λ ∈ Π ess (L) is given and will focus on the asymptotic structure of the corresponding approximate eigenfunctions.We thus consider the asymptotic expansions where λ 0 , λ 1 ∈ C are treated as parameters and f 0 m , f 1 m : [0, 1] → C are unknown functions.Plugging these expansions into system (21) and collecting terms proportional to the highest powers of m we obtain It follows immediately from (24a) that f 0 m−1 = f 0 m+1 .Since this analysis does not distinguish between even and odd values of m, we also deduce that which is an inhomogeneous first-order equation defining the leading-order term f 0 m (r) in (23) in terms of f 1 m (r).Without loss of generality the boundary condition (21b) can be replaced with f 0 m (0) = 1.The solution of ( 25) is then a sum of two parts: the solution of the homogeneous equation obtained by setting the RHS to zero and a particular integral corresponding to the actual RHS.Since at this level the expression f 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 where The limiting (as r → 1) behavior of functions (27a)-( 27b) exhibits an interesting depen-dence on λ 0 , namely, In particular, the limiting value of I r (r) as r → 1 changes when ℜ(λ 0 ) = 4, which defines the right boundary of the essential spectrum in the present problem, cf. ( 8).Both I r (r) and I i (r) diverge as O(1/(1 − r)) when r → 1 which means that the integrals under the exponentials in ( 26), 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 0 m (r), the factor depending on I r (r) determines its growth as r → 1: we see that |f 0 m (r)| becomes unbounded in this limit when ℜ(λ 0 ) < 4 and approaches zero otherwise.The real and imaginary parts of f 0 m (r) obtained for different eigenvalues λ 0 are shown in figures 3a,b, where it is evident that both the unbounded growth and the oscillations of f 0 m (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 1 m (r) is rather difficult to compute and we do not attempt this here.If f 1 m−1 ̸ = f 1 m+1 , the solution of equation ( 25) will also include some extra terms in addition to (26)-( 27) which would correspond to another possible family of approximate eigenfunctions.However, as will be evident from the discussion below, the solutions given in ( 26)-( 27) capture the relevant behavior.Finally, in view of ansatz (20), the leading-order approximations to eigenfunctions are obtained multiplying the function f 0 m (r) by cos(mθ) or sin(mθ) with m → ∞ which introduces rapid oscillations in the azimuthal direction.
We thus conclude that when ℜ(λ 0 ) < 4, the solutions of eigenvalue problem (15) constructed in the form (20) include functions dominated by short-wavelength oscillations whose asymptotic, as m → ∞, structure involves oscillations in both the radial and azimuthal directions and are localized near the boundary ∂A 0 .Since as a result their pointwise values on ∂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 (21c)-(21d), which is consistent with the fact that they represent approximate eigenfunctions associated with the essential spectrum Π ess (H) of the 2D linearized Euler operator.In order to find solutions of eigenvalue problem (15) which do satisfy all the boundary conditions we have to solve this problem numerically which is done next.

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

Discretization of Eigenvalue Problem (15)-(16)
Eigenvalue problem (15)-( 16) is solved using the spectral collocation method proposed by Fornberg (1996), see also the discussion in Trefethen (2000), which is based on a tensor grid in (r, θ).The discretization in θ involves trigonometric (Fourier) interpolation, whereas that in r is based on Chebyshev interpolation where we take r ∈ [−1, 1] which allows us to avoid collocating (15a) at the origin when the number of grid points is even.Since then the mapping between (r, θ) and (x 1 , x 2 ) is 2-to-1, the solution must be constrained to satisfy the condition which is fairly straightforward to implement (Trefethen, 2000).
In contrast to (15a), the boundary condition (15b) 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 ∈ [−1, 1] and θ ∈ [0, 2π] are linked and both given by N which is an even integer.The resulting algebraic eigenvalue problem then has the form where ψ ∈ C N 2 is the vector of approximate nodal values of the eigenfunction and H ∈ R N 2 ×N 2 the matrix discretizing the operator H, cf.(15a), obtained as described above.Problem (30) is implemented in MATLAB and solved using the function eig.
The discretization of all operators in H, cf (15), 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 H and hence also the matrix H are nonnormal and singular, the numerical conditioning of problem (30) 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 H by performing projections on a certain number N C of eigenfunctions associated with the eigenvalue λ = 0 (they are obtained by solving problem (19) with different source terms ϕ C , C = 2, 3, . . ., N C + 1, cf. ( 47)).However, solutions of problem (30) obtained in this way were essentially unchanged as compared to the original version.Moreover, in addition to examining the behavior 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 (2017).Increasing the arithmetic precision up to O(10 2 ) significant digits was also not found to have a noticeable effect on the results obtained with small and medium resolutions N ≤ 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- § 2.2, we know the spectrum of the operator H includes essential spectrum in the form of a vertical band in the complex plane |ℜ(z)| ≤ 4, z ∈ 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 (30) is finitedimensional and therefore can only have 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 Π ess (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 eigenfunction, we also solve a regularized version of problem (15) in which ψ is replaced with ψ δ := R −1 δ ψ, where R δ := (Id −δ 2 ∆), δ > 0 is a regularization parameter and the inverse of R δ is defined with the homogeneous Neumann boundary conditions.The regularized version of the discrete problem (30) then takes the form where the subscript δ denotes regularized quantities and R δ is the discretization of the regularizing operator R δ .Since the operator R −1 δ can be interpreted as a low-pass filter with the cut-off length given by δ, the effect of this regularization is to smoothen the eigenvectors by filtering out components with wavelengths less than δ.Clearly, in the limit when δ → 0 the original problem (30) is recovered.An analogous strategy was successfully employed by Protas & Elcrat (2016) in their study of the stability of Hill's vortex where the eigenfunctions also turned out to be singular distributions.

Solution of the Time-Dependent Problem (1)-(2)
The 2D Euler system (1)-( 2) is transformed to the frame of reference moving with velocity −U e 1 and rewritten in terms of the "perturbation" vorticity ω(t, x) := ω(t, x) − ω 0 (x) and the corresponding perturbation streamfunction ψ(t, x), such that it takes the form, cf. ( 7), To facilitate solution of this system with a Fourier pseudospectral method (Canuto et al., 1988), we approximate the unbounded domain with a 2D periodic box Ω ≈ 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, 2001).We note that this is a standard approach and has been successfully used in earlier studies of related problems (Nielsen & Rasmussen, 1997;van Geffen & van Heijst, 1998;Billant et al., 1999;Donnadieu et al., 2009;Brion et al., 2014;Jugier et al., 2020).Since the instability has the form of shortwavelength 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 ω(t, x) , where M is the number of grid points in each direction in T 2 .Substitution of expansion (33) into (32a) yields a system of coupled ordinary differential equations describing the evolution of the expansion coefficients ω k (t), k ∈ V M , which is integrated in time using the RK4 method.Product terms in the discretized equations are evaluated in the physical space with the exponential filter proposed by Hou & Li (2007) used in lieu of dealiasing.We use a massively-parallel implementation based on MPI with Fourier transforms computed using the FFTW library (Frigo & Johnson, 2003).Convergence of the results with refinement of the resolution M and of the time step ∆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π.

Solution of the Eigenvalue Problem
In this section we describe solutions of the discrete eigenvalue problem (30) and its regularized version (31).We mainly focus on the symmetric dipole with η = 0, cf.figure 1a.In order to study dependence of the solutions on the numerical resolution, problems (30)-( 31) were solved with N ranging from 20 to 260, where the largest resolution was limited by the amount of RAM memory available on a single node of the computer cluster we had access to.The discrete spectra of problem (30) obtained with N = 40, 80, 160, 260 are shown in figures 4a-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 |λ| = O(N ) when the resolution is refined; they are not shown in figures 4a-d).We see that as N increases the cloud formed by the complex eigenvalues remains restricted to the band −2 ⪅ ℜ(λ) ⪅ 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 figures 5a-d we show the eigenvalue density defined as the number of eigenvalues in a small rectangular region of the complex plane, i.e., ) where ∆λ r , ∆λ i ∈ R are half-sizes of a cell used to count the eigenvalues with ∆λ i ≈ 500∆λ 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 µ(z) increases near the origin.
As discussed in § 2.1, a key question concerning the linear stability of 2D Euler flows is the existence of point spectrum Π 0 (L) of the linear operator L, cf. ( 7).However, the usual approach based on discretizing the continuous eigenvalue problem ( 15)-( 16), 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 (30) 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 Π ess (L) ever more densely in this limit.In this way we have found a single eigenvalue, denoted λ 0 , which together with its negative −λ 0 and complex conjugates ±λ * 0 , satisfy the above condition, see Table 1.As is evident from this table, the differences between the real parts of λ 0 computed with different resolutions N are very small and just over 1%, although the variation of the imaginary part is larger.Moreover, as will be discussed in § 6, λ 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, 2014).In particular, if the maximum period is bounded τ max < ∞, this implies the presence of a 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 Π ess (L) includes the entire imaginary axis iR.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.
Finally, we analyze eigenvectors of problem ( 30) and choose to present them in terms of vorticity, i.e., we show ω i = −∆ ψ i , where the subscript i = 0, 1, 2 enumerates the corresponding eigenvalues.First, in figures 7a,b,c we illustrate the convergence pattern of the eigenvector ω 0 corresponding to the eigenvalue λ 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 | ω 1 (r, θ)| of the eigenvector grows rapidly near the boundary, i.e., as r → 1, which is consistent with the behavior of the function f 0 m (r) describing the asymptotic solution, cf.expressions ( 26)-( 28) and figures 3(a,b).In addition, a rapid variation of | ω 1 (r, θ)| in the azimuthal coordinate θ is also evident for r ≲ 1.However, something that could not be discerned by the asymptotic analysis is that these oscillations are mostly concentrated near the azimuthal angles θ = ±π/4, ±3π/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 (15), 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 will be shown in § 6, even if they are not fully resolved, the eigenvectors computed here still contain useful information.
Next, in figures 8a,c,e we compare the real parts of the eigenvectors associated with different eigenvalues: the complex eigenvalue λ 0 corresponding to the exponentially growing mode (already shown in figure 7), a purely real eigenvalue λ 1 and a purely imaginary eigenvalue λ 2 .We see that while these eigenvectors are qualitatively similar and share the features described above, the eigenvector ω 0 is symmetric with respect to the flow centerline, whereas the eigenvectors ω 1 and ω 2 are antisymmetric.Another difference is that in the eigenvector ω 0 associated with the eigenvalue λ 0 the oscillations are mostly concentrated near the azimuthal angles θ = ±π/4, ±3π/4, cf.figure 8a; on the other hand, in the eigenvectors ω 1 and ω 2 the oscillations are mostly concentrated near the stagnation points x a and x b , cf. figure 8(c,e).
The numerical approximations of the eigenvectors are characterized by short-wavelength oscillations.Here, "short-wavelength" means that a significant variation of the magnitude | ω(r, θ)| of the eigenvector with respect to both r and θ 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 ω 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 ∂A 0 .Since due to their highly oscillatory nature the eigenvectors shown in figures 8a,c,e are not fully resolved, in figures 8b,d,f we show the corresponding eigenvectors of the regularized eigenvalue problem (31) where we set δ = 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 (31) are distributed following a similar pattern as revealed by the eigenvalues of the original problem (30), cf.figures 4b and 5b.In particular, for the main eigenvalue of interest λ 0 , cf.Table 1, the difference with respect to the corresponding eigenvalue of the regularized problem λ δ,0 is rather small and we have |ℜ(λ 0 − λ δ,0 )|/ℜ(λ 0 ) ≈ 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 Π ess (L).
Solution of the discrete eigenvalue problem (30) for asymmetric dipoles with η > 0 leads to eigenvalue spectra and eigenvectors qualitatively very similar to those shown in figures 4a-d and 8a,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 centerline.
Figure 7: Real parts of the eigenvector ω 0 corresponding to the eigenvalue λ 0 , cf.Table 1, and representing an exponentially growing mode obtained by solving the discrete eigenvalue problem (30) with different resolutions N .The grids covering the surface plots represent the discretizations of the domain A 0 used for different N .

Solution of the Evolution Problem
As in § 5, we focus on the symmetric case with η = 0.The 2D Euler system (1)-( 2) is solved numerically as described in § 4.2 with the initial condition for the perturbation vorticity ω(t, x) given in terms of the eigenvectors shown in figures 8a-f, i.e., Unless indicated otherwise, the numerical resolution is M = 512 grid points in each direction.By taking ε = 10 −4 we ensure that the evolution of the perturbation vorticity is effectively linear up to t ≲ 70 and to characterize its growth we define the perturbation enstrophy as The evolution of this quantity is shown in figure 9a 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( λt), where the growth rate (found via a least-squares fit) is λ ≈ 0.127 and is essentially equal to the real part of the eigenvalue λ 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 ω 0 and ω δ,0 are used as the initial conditions in (35) (in fact, in the latter case the transient is barely present).The reason for this behavior is that ω 0 is the sole true eigenvector of the operator L, whereas ω 1 and ω 2 are only approximate eigenvectors associated with the (approximate) eigenvalues λ 1 and λ 2 belonging to the essential spectrum Π ess (L) rather than to the point spectrum Π 0 (L).
As a result, ω 0 represents the only linearly growing mode, such that when ω 1 , ω 2 or any other approximate eigenvector is used as the initial condition in (35), a transient behavior ensues where the solution ω(t) of system (32) approaches the trajectory involving the growing mode ℜ e λ 0 t ω 0 .Hereafter we will focus on the flow obtained with the initial condition (35) given in terms of the eigenfunction ω 0 , cf. figure 8a.The effect of the numerical resolution N used in the discrete eigenvalue problem (30) is analyzed in figure 9b, where we show the perturbation enstrophy (36) in the flows with the eigenvector ω 0 used in the initial conditions (35) computed with different N .We see that refined resolution leads to a longer transient period while the rate of the exponential growth λ 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 (35) and of the perturbation vorticity ω(t, x) at different times t ∈ (0, 60] is shown in figure 10 as a function of the wavenumber k := |k|.It is defined as where σ is the azimuthal angle in the 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 k as a continuous rather than discrete variable).Since its enstrophy spectrum is essentially independent of the wavenumber k, the eigenvector ω 0 in the initial condition ( 35) turns out to be a distribution rather than a smooth function.
The enstrophy spectra of the perturbation vorticity ω(t, 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 ⪅ t ≤ 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 9a the relative change of the total energy Ω |u(t, x)| 2 dx and of the total enstrophy Ω ω(t, x) 2 dx, which are conserved quantities in the Euler system (1)-( 2), is at most of order 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 [ ω(t)] 0 = [ ω 0 (t)] 0 + [ ω i (t)] 0 , it remains zero by construction throughout the entire flow evolution.We now go on to discuss the time evolution of the perturbation vorticity in the physical space and in figures 11a and 11b we show ω(t, 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 ⪅ t ≤ 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 θ ≈ ±π/4, ±3π/4, i.e., in the regions of the vortex boundary where most of the short-wavelength oscillations evident in the eigenvector ω 0 are localized, cf.figure 8a.With thickness on 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 centerline for all times and since the vorticity ω 0 of the base flow is antisymmetric, the resulting total flow ω(t, x) does not possess any symmetries.The perturbation vorticity ω(t, x) realizing the exponential growth in the flows corresponding to the initial condition involving the eigenvectors ω 1 and ω 2 (and their regularized versions ω δ,1 and ω δ,2 ) is essentially identical to the perturbation vorticity shown in figure 11b, although its form during the transient regime can be quite different.In particular, the perturbation eventually becomes symmetric with respect to the flow centerline even if the initial condition ( 35) is antisymmetric.The same is true for flows obtained with initial condition corresponding to all approximate eigenvalues other than λ 1 and λ 2 (not shown here for brevity).We did not attempt to study the time evolution of asymmetric dipoles with η > 0 in (5a), since their vorticity distributions are discontinuous making computation of such flows using the pseudospectral method described in § 4.2 problematic.

Discussion and Final Conclusions
In this study we have considered an open problem concerning the linear stability of the Lamb-Chaplygin dipole which is a classical equilibrium solution of the 2D Euler equation in an unbounded domain.We have considered its stability with respect to 2D circulationpreserving perturbations and while our main focus was on the symmetric configuration with η = 0, cf.figure 1a, we also investigated some aspects of asymmetric configurations with η > 0. Since the stability of the problem posed on a unbounded domain is difficult to study both with asymptotic methods and numerically, we have introduced an equivalent formulation with all relations defined entirely within the compact vortex core A 0 , which was accomplished with the help of a suitable D2N map accounting for the potential flow outside the core, cf.Appendix A. The initial-value problem for the 2D Euler equation with a compactly supported initial condition is of a free-boundary type since the time evolution of the vortex boundary ∂A(t) is a priori unknown and must be determined as a part of the solution of the problem.This important aspect is accounted for in our formulation of the linearized problem, cf.relation (10).The operator representing the 2D Euler equation linearized around the Lamb-Chaplygin dipole has been shown to have an infinite-dimensional null space Ker(L) and the eigenfunctions ψ C , C = 2, 3, . . ., spanning this null space, cf.figures 2a-d, can potentially be used to search for nearby equilibrium solutions.
We have studied the linear stability of the Lamb-Chaplygin dipole using a combination of asymptotic analysis ( § 3) and numerical computation ( § 5) employed to construct approximate solutions of the eigenvalue problem (15) together with the numerical time-integration of the 2D Euler system (32) in § 6.These three approaches offer complementary insights reinforcing the main conclusion, namely, that the Lamb-Chaplygin dipole is linearly unstable with the instability realized by a single eigenmode ω 0 , cf. figure 8a, featuring high-frequency oscillations localized near the vortex boundary ∂A 0 and the corresponding eigenvalue λ 0 embedded in the essential spectrum Π ess (L) of the linearized operator L. In other words, there is no "smallest" length scale characterizing the unstable mode, which is why it cannot be accurately resolved using any finite numerical resolution.This is one of the reasons why this form of instability specific to the inviscid evolution is so fundamentally different from the mechanisms underlying the growth of perturbations during the viscous evolution of the dipole that were observed in all earlier studies (Nielsen & Rasmussen, 1997;van Geffen & van Heijst, 1998;Billant et al., 1999;Donnadieu et al., 2009;Brion et al., 2014;Jugier et al., 2020).
An approximate solution of eigenvalue problem (15) obtained in § 3 using an asymptotic technique reveals the existence of approximate eigenfunctions in the form of shortwavelength oscillations localized near the vortex boundary ∂A 0 .Remarkably, eigenfunctions with such properties exist when ℜ(λ 0 ) < 4, i.e., when λ 0 is in the essential spectrum Π ess (H) of the 2D linearized Euler operator and it is interesting that the asymptotic solution has been able to capture this value exactly.We remark that with exponential terms involving divergent expressions as arguments, cf. ( 26), this approach has the flavor of the WKB analysis.We note that while providing valuable insights about the structure of the approximate eigenvectors the asymptotic analysis developed in § 3 does not allow us to determine the eigenvalues of problem (15), i.e., λ 0 serves as a parameter in this analysis.Moreover, since the obtained approximate solution represents only the asymptotic (in the short-wavelength limit m → ∞) structure of the eigenfunctions, it does not satisfy the boundary conditions (21c)-(21d).To account for these limitations, complementary insights have been obtained by solving eigenvalue problem (15) numerically as described in § 4.1.
Our numerical solution of eigenvalue problem (15) obtained in § 5 using different resolutions N yields results consistent with the general mathematical facts known about the spectra of the 2D linearized Euler operator, cf.§ 2.1.In particular, these results feature eigenvalues of the discrete problem (30) filling ever more densely a region around the origin which is bounded in the horizontal (real) direction and expands in the vertical (imaginary) direction as the resolution N is increased, which is consistent with the existence of an essential spectrum Π ess (H) in the form of a vertical band with the width determined by the largest Lyapunov exponent of the flow, cf. ( 8).The corresponding eigenvectors are dominated by short-wavelength oscillations localized near the vortex boundary ∂A 0 , a feature that was predicted by the asymptotic solution constructed in § 3.However, solutions of the evolution problem for the perturbation vorticity with the initial condition (35) corresponding to different eigenvectors obtained from the discrete problems ( 30)-( 31) reveal that λ 0 (and its complex conjugate λ * 0 ) are the only eigenvalues associated with an exponentially growing mode with a growth rate equal to the real part of the eigenvalue, i.e., for which λ ≈ ℜ(λ 0 ).When eigenvectors associated with eigenvalues other than λ 0 or λ * 0 are used in the initial condition (35), the perturbation enstrophy (36) reveals transients of various duration followed by exponential growth with the growth rate again given by ℜ(λ 0 ).This demonstrates that ±λ 0 and ±λ * 0 are the only "true" eigenvalues and form the point spectrum Π 0 (H) of the operator associated with the 2D Euler equation linearized around the Lamb-Chaplygin dipole.On the other hand, all other eigenvalues of the discrete problems ( 30)-( 31) can be interpreted as numerical approximations to approximate eigenvalues belonging to the essential spectrum Π ess (H).More precisely, for each resolution N the eigenvalues of the discrete problems other than ±λ 0 and ±λ * 0 approximate a different subset of approximate eigenvalues in the essential spectrum Π ess (H) and the corresponding eigenvectors are approximations to the associated approximate eigenvectors.This interpretation is confirmed by the eigenvalue density plots shown in figures 5a-d and is consistent with what is known in general about the spectra of the 2D linearized Euler operator, cf.§ 2.1.
In figure 9a we noted that when the initial condition ( 35) is given in terms of the eigenvector ω 0 , the perturbation enstrophy E(t) also exhibits a short transient before attaining exponential growth with the rate λ ≈ ℜ(λ 0 ).The reason for this transient is that, being non-smooth, the eigenvector ω 0 is not fully resolved, which is borne out in figure 10 (in fact, due to the distributional nature of this and other eigenvectors, they cannot be accurately resolved with any finite resolution).Thus, this transient period is needed for some underresolved features of the perturbation vorticity to emerge, cf.figure 11a vs. figure 11b.However, we note that in the flow evolution originating from the eigenvector ω 0 the transient is actually much shorter than when other eigenvectors are used as the initial condition (35), and is nearly absent in the case of the regularized eigenvector ω δ,0 .We emphasize that non-smoothness of eigenvectors associated with eigenvalues embedded in the essential spectrum is consistent with the known mathematical results predicting this property (Lin, 2004).Interestingly, the eigenfunctions ψ C , C = 2, 3, . . ., associated with the zero eigenvalue λ = 0 are smooth, cf.figures 2a-d.
We also add that there are analogies between our findings and the results of the linear stability analysis of Hill's vortex with respect to axisymmetric perturbations where the presence of both the continuous and point spectrum was revealed, the latter also associated with non-smooth eigenvectors (Protas & Elcrat, 2016).
In the course of the linear evolution of the instability the vortex region A(t) changes shape as a result of the ejection of thin vorticity filaments from the vortex core A 0 , cf. figures 11a,b.However, both the area |A(t)| of the vortex and its total circulation Γ are conserved at the leading order, cf. ( 11) and ( 14).We reiterate that the perturbation vorticity fields shown in figures 11a,b were obtained with underresolved computations and increasing the resolution M would result in the appearance of even finer filaments such that in the continuous limit (M → ∞) some of the filaments would be infinitely thin.
In this study we have considered the linear stability of the Lamb-Chaplygin dipole with respect to 2D perturbations.It is an interesting open question how the picture presented here would be affected by inclusion of 3D effects.We are also exploring related questions in the context of the stability of other equilibria in 2D Euler flows, including various cellular flows.
where α k , β k ∈ R, k = 1, 2, . . ., are expansion coefficients to be determined and the constant term is omitted since we adopt the normalization ∂A 0 f ′ (s) ds = 0.The boundary value f ′ of the perturbation streamfunction on ∂A 0 serves as the argument of the D2N operator, cf.(9c).Expanding it in a Fourier series gives where f c k , f s k ∈ R, k = 1, 2, . . ., are known coefficients.Then, using the boundary condition ψ ′ 2 (1, θ) = f ′ (θ), θ ∈ [0, 2π], cf.(9c), the corresponding Neumann data can be computed as which expresses the action of the D2N operator M on f ′ .In order to make this expression explicitly dependent on f ′ , rather than on its Fourier coefficients as in (40) , we use the formulas for these coefficients together with their approximations based on the trapezoidal quadrature (which are spectrally accurate when applied to smooth periodic functions (Trefethen, 2000)) f ′ (θ l ) cos(kθ l ), (41a) where {θ l } N l=1 are grid points uniformly discretizing the interval [0, 2π].Using these relations, the D2N map (40) truncated at N/2 Fourier modes and evaluated at the grid point θ j can be written as where k [cos(kθ j ) cos(kθ l ) + sin(kθ j ) sin(kθ l )] (43) are entries of a symmetric matrix M ∈ R N ×N approximating the D2N operator.

Figure 1 :
Figure 1: Streamline pattern inside the vortex core A 0 of (a) a symmetric (η = 0) and (b) asymmetric (η = 1/4) Lamb-Chaplygin dipole.Outside the vortex core the flow is potential.The thick blue line represents the vortex boundary ∂A 0 whereas the red symbols mark the hyperbolic stagnation points x a and x b .

Figure 3 :
Figure 3: Radial dependence (a) of the eigenvectors f 0 m (r) associated with real eigenvalues λ 0 = 2 (red solid line) and λ 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 0 m (r) associated with complex eigenvalue λ 0 = 3 + 10i.Panel (b) shows the neighbourhood of the endpoint r = 1.

Figure 4 :
Figure 4: Eigenvalues obtained by solving the discrete eigenvalue problem (30) with different indicated resolutions N .The eigenvalues ±λ 0 and ±λ *0 which converge to well-defined limits as the resolution N is refined, cf.Table1, are marked in red.The eigenvalue λ 0 is associated with the only linearly unstable mode, cf.§ 6.

Figure 8 :
Figure 8: Real parts of the eigenvectors corresponding to the indicated eigenvalues obtained by solving (a,c,e) eigenvalue problem (30) and (b,d,f) the regularized problem (31) using the resolution N = 80.The grid shown on the surface represents the discretization of the domain A 0 used in the numerical solution of problems (30) and (31).

Figure 11 :
Figure 11: Perturbation vorticity ω(t, x) in the flow corresponding the initial condition (35) involving the eigenvector ω 0 during (a) the transient regime and (b) the period of exponential growth.

Table 1 ,
are marked in red.The eigenvalue λ 0 is associated with the only linearly unstable mode, cf.§ 6.

Table 1 :
Eigenvalue λ 0 associated with the linearly growing mode, cf.§ 6, obtained by solving the discrete eigenvalue problem (30) with different resolutions N .