## Abstract

In this investigation, we revisit the question of linear stability analysis of two-dimensional steady Euler flows characterized by the presence of compact regions with constant vorticity embedded in a potential flow. We give a complete derivation of the linearized perturbation equation which, recognizing that the underlying equilibrium problem is of a free-boundary type, is carried out systematically using methods of shape-differential calculus. Particular attention is given to the proper linearization of contour integrals describing vortex induction. The thus obtained perturbation equation is validated by analytically deducing from it stability analyses of the circular vortex, originally due to Kelvin, and of the elliptic vortex, originally due to Love, as special cases. We also propose and validate a spectrally accurate numerical approach to the solution of the stability problem for vortices of general shape in which all singular integrals are evaluated analytically.

## 1. Introduction

In this work, we are interested in the linear stability of solutions of two-dimensional Euler equations that are steady in the appropriate frame of reference and feature compact regions with constant vorticity embedded in an otherwise potential flow. Such solutions arise as inviscid limits of actual viscous flows, and therefore find numerous applications in different areas of fluid mechanics. The first linear stability analysis of such a flow is attributed to Kelvin [1], see also Lamb [2], Batchelor [3] and Saffman [4], who identified the dispersion relation of neutrally stable travelling waves with *m*-fold symmetry moving along the perimeter of a circular vortex. This problem was revisited by Baker [5], who used a contour-dynamics approach involving integrals with singular kernels defined on the vortex boundary. The stability of a rotating ellipse, the so-called Kirchhoff vortex [6], was analysed in a seminal study by Love [7]. The effect of an ambient strain field was investigated by Moore & Saffman [8], whereas Mitchell & Rossi [9] explored the relation between linear instabilities and the long-term nonlinear evolution of elliptic vortices. In this context, we also mention the analytical study of Guo *et al.* [10] where an integro-differential perturbation equation similar to ours is proposed for the case of the elliptic vortex. The linear stability of more complex vortex configurations, such as polygonal arrays of corotating vortices and translating vortex pairs, was investigated by Dritschel [11–13], see also Dritschel & Legras [14]. A noteworthy feature of their approach is that they also used a *continuous* perturbation equation independent of a particular discretization employed to obtain the equilibrium solution. Problems related to the stability of polygonal vortex arrays were also studied by Dhanak [15]. The linear stability of the so-called V-states [16], corotating vortex patches and infinite periodic arrays of vortices was investigated in detail by Kamm [17], see also Saffman [4]. This approach was based on the representation of the solution in terms of the Schwarz function, and required discretization and numerical differentiation in order to obtain the perturbation equation. A discrete form of the perturbation equation was also used by Elcrat *et al.* [18] in their investigation of the linear stability of vortices in a symmetric equilibrium with a circular cylinder and a free stream at infinity. There have also been a number of analytical investigations concerning the stability, both in the linear and nonlinear setting, of flows obtained as perturbations of the Rankine vortex [19–24]. Another family of approaches is based on variational energy arguments going back to Kelvin. They were initially investigated by Saffman & Szeto [25], Dritschel [11,26], Cerrtelli & Williamson [27] and Fukumoto & Moffatt [28], and were more recently pursued by Luzzatto-Fegiz & Williamson [29–34]. They rely on global properties of the excess energy versus velocity impulse diagrams and provide partial information about the linear stability properties without the need to actually perform a full linear stability analysis.

To the best of our knowledge, there is no complete derivation and validation of the stability equations for two-dimensional vortex patches available in the literature. The main goal of this paper is thus to fill the gap by proposing a systematic approach to the study of the linear stability of vortices which is quite general in the sense that it is based on the continuous, rather than discrete, formulation of the governing equations and a rigorous treatment of boundary deformations. This aspect, namely, that linearization of differential and integral expressions with respect to perturbations of the shape of the domain requires special treatment with methods other than the traditional calculus of variations, does not seem to have received much attention in the vortex dynamics literature. Therefore, the proposed approach may serve as a template for studying more complicated problems such as the stability of vortex equilibria in three dimensions. Recognizing the Euler equation describing finite-area vortices as a *free-boundary* problem, we apply methods of shape-differential calculus to derive a perturbation equation in the form of an integro-differential equation defined on the vortex boundary. We emphasize that, in contrast to Kamm [17] and Elcrat *et al.* [18], numerical differentiation is not required at any stage of the derivation, but only in the evaluation of the coefficients of the perturbation equation. In this sense, the formulation we propose is a contour-dynamics analogue of the Orr–Sommerfeld equation used to study the stability of viscous parallel flows [35]. In this work, we also demonstrate using analytical calculations how one can reproduce the dispersion relations obtained by Kelvin in the stability analysis of the Rankine vortex [1] and by Love in the stability analysis of the elliptic vortex [7] from our general perturbation equation. Generalization of our approach to the stability analysis of more complicated vortex flows, including three-dimensional axisymmetric configurations, is deferred to future research.

We are interested in developing a general approach to analysing the linear stability of two-dimensional steady-state solutions of the Euler equation, which is usually written in the following form [36]:
1.1aand
1.1bwhere *ψ* is the streamfunction, is the flow domain, whereas *ψ*_{b} represents the boundary condition for the streamfunction. The function is *a priori* undetermined, and any choice of *F* will yield a steady solution of the two-dimensional Euler equation. One choice of the function *F*, motivated by the flow models arising from the Prandtl–Batchelor theory [37], is
1.2where *H*(⋅) is the Heaviside function, whereas are, respectively, the vorticity inside the vortex and the value of the streamfunction at its boundary. The profile given in (1.2) corresponds to finite-area vortices with constant vorticity *ω* embedded in an otherwise potential flow. This formulation can be made more general by including in (1.2) a term proportional to the Dirac delta function representing a vortex sheet on the vortex boundary [38], but we do not consider this in the present study. A salient feature of problem (1.1)–(1.2) is that the shape of the vortex region, which we will denote *A*, is *a priori* unknown and therefore must be determined as a part of the solution to the problem. System (1.1)–(1.2) is thus a *free-boundary* problem, and in order to emphasize this fact that will play a central role in our analysis, we rewrite it in the following form, which makes the free-boundary property more evident:
1.3a
1.3b
1.3c
1.3d
and
1.3ewhere and are the restrictions of the streamfunction *ψ* to, respectively, the rotational and irrotational parts of the flow domain and **n** is the unit vector normal to the vortex boundary (‘’ means ‘equal to by definition’). Solutions of system (1.3) depend on two parameters *ω* and *ψ*_{0} or, equivalently, the vortex area and its circulation . There is evidence that by continuing the solutions of (1.3) corresponding to different flow configurations, such that *γ*=const. and , one can obtain the corresponding point-vortex equilibria [39,41]. In agreement with earlier studies [13,17], we will assume in this investigation that both parameters |*A*| and *γ* are fixed, so that the perturbations will be allowed to modify the shape of the vortex region *A*, leaving its area and circulation unchanged. In other words, no vorticity is created by the perturbations. Derivation of the Jacobian, needed for the linear stability analysis of a given problem, requires differentiation of the governing equation with respect to the state variables. In the present problem, since governing system (1.3) is of the free-boundary type, this is a rather delicate issue because the shape of the vortex *A* can be regarded as a ‘state variable’. A proper way to deal with this problem is therefore to use *shape-differential calculus* [42,44], which is a suite of mathematical techniques for differentiation of various objects such as differential equations, integrals, etc., defined on variable domains. Application of the shape-differentiation methods makes it possible to derive a continuous perturbation equation for vortices with very general shapes in a rigorous and systematic manner, which is the main contribution of this work. A key element is a proper linearization of the contour-dynamics equations with respect to arbitrary boundary perturbations. For derivations of the different shape-differentiation identities we use in this study, the reader is referred to monographs [42–44] and a review paper [45] for a survey of applications of the shape calculus to different problems in computational fluid dynamics. The structure of the present paper is as follows: in the next section, we derive a general form of the perturbation equation; in §§3 and 4, we use this perturbation equation to obtain the dispersion relations for the circular and elliptic vortices; comments about an accurate numerical technique for the solution of the perturbation equation are presented in §5, together with some computational results; and conclusions are deferred to §6. Some technical results are given in Appendix A.

## 2. Derivation of the perturbation equation

In this section, we develop a systematic approach to derive the perturbation equation characterizing the stability of steady vortex configurations. For brevity of notation, in the following, we will use the complex number representation of vector quantities. This is indicated by using ordinary letters for vector quantities (e.g. for and similarly for other vector quantities). Juxtaposition (i.e. *z*_{1}*z*_{2} for some *z*_{1}, ) indicates complex multiplication of *z*_{1} and *z*_{2}. However, when we need equations written in the vector form, as we will in referencing results from shape calculus, we will denote vector quantities in bold face and inner products of vectors using the usual dot notation. Thus, for vectors and the corresponding complex numbers, we have the identity , where the overbar denotes complex conjugation. To fix attention, we will focus on the simplest, yet generic, case of a single vortex region *A* in equilibrium with an external velocity field *U*=*u*−*iv*, which may represent a free stream, the influence of other vortices, or may result from transformation of the problem to a moving (i.e. translating or rotating) coordinate system in which a relative equilibrium exists. We will consider a point *z*=*x*+*iy* on the vortex boundary ∂*A* (figure 1) that is assumed to have counterclockwise orientation. Using techniques of vortex dynamics [4], the velocity of a particle located at the point *z* can be expressed in the fixed frame of reference as^{1}
2.1where the integral on the right-hand side (r.h.s.) has a removable singularity (*w* is a complex integration variable). Starting from relation (2.1), one can pursue our objective with one of the two alternative approaches, namely (figure 1):

— The Lagrangian approach in which one follows the trajectory of a fluid particle on the vortex boundary; we note that displacement of a material particle in general also involves a component tangential to the vortex boundary, which does not lead to deformations of the contour [5]; and

— the geometric approach, in which one considers displacements of the points on the boundary in the direction normal to the boundary only [13].

Since the Lagrangian approach involves additional (tangential) degrees of freedom, the stability results obtained with the geometric approach represent a subset of the results obtained with the Lagrangian approach. While both approaches can be pursued using methods of shape differentiation, in this work, we will focus on the geometric approach as it leads to somewhat simpler calculations.

As a first step, we transform relation (2.1) to a suitable moving (i.e. rotating or translating) coordinate system, characterized by the rigid-body velocity *V* _{0}(*z*), in which a relative equilibrium is attained. Denoting *ζ*(*τ*) as the position of the point *z*(*τ*) in this moving coordinate system, the left-hand side (l.h.s.) of (2.1) becomes , so that relation (2.1) can be rewritten as
2.2

In equation (2.2) and hereafter, we use both *z*(*τ*) and *ζ*(*τ*), and note that, since *ζ*(0)=*z*(0), either coordinate can be used on the r.h.s. of (2.2) to derive the perturbation equation. Let *n*=*n*_{x}+i*n*_{y} be the unit vector normal to the boundary ∂*A* and pointing outside of *A*. The relative equilibrium is then characterized by the condition
2.3which, for given vorticity *ω* and vortex area |*A*|, can be interpreted as a condition on the shape of the vortex boundary ∂*A*. Next, we introduce a parametrization of the contour *z*=*z*(*τ*,*s*)∈∂*A*, where is the time, such that *s*∈[0,*L*], in which *L* denotes the length of the boundary ∂*A*. We will also adopt the convention that the superscript *ϵ*, where 0<*ϵ*≪1, will denote quantities corresponding to the perturbed boundary, so that *ζ*=*ζ*^{ϵ}|_{ϵ=0} and *n*=*n*^{ϵ}|_{ϵ=0} are the quantities corresponding to the (relative) equilibrium. Then, points on the perturbed vortex boundary can be represented as
2.4where *r*(*τ*,*s*) represents the ‘shape’ of the perturbation. We note that, without affecting the final result, *n*(*s*) in the last term in (2.4) could be replaced with its perturbed counterpart *n*^{ϵ}(*τ*,*s*). Using equation (2.1), we thus deduce
2.5from which equilibrium condition (2.3) is obtained by setting *ϵ*=0. In (2.5), the integral operator is defined in terms of a contour integral. In order to be able to use shape-calculus formulae, we need to write such integrals as integrals with respect to the arc length. If *F* is a complex-valued function and *Γ* a closed *C*^{1} curve, we have
2.6where *t* is the unit tangent vector in the counterclockwise direction.

The perturbation equation is obtained by linearizing equation (2.5) around the equilibrium configuration characterized by (2.3), which is equivalent to expanding (2.5) in powers of *ϵ* and retaining the first-order terms. Since equation (2.5) involves perturbed quantities defined on the perturbed vortex boundary ∂*A*^{ϵ}, the proper way to obtain this linearization is using methods of shape-differential calculus [43]. Differentiating the l.h.s. in (2.5) and setting *ϵ*=0, we thus obtain
2.7

It can be shown using shape calculus [43] that 2.8

Noting also that (d**n**/d*τ*)|_{τ=0}=−(**∇**_{Γ}**V**)^{T}**n**, where **∇**_{Γ} denotes the tangential gradient [43], has only a component tangential to the boundary ∂*A* and expanding d*r*/d*τ*=∂*r*/∂*τ*+*V* _{1t}(∂*r*/∂*s*), where and ∂/∂*s* is the derivative with respect to the arc length coordinate, we can rewrite (2.7) as
2.9

For the r.h.s. in (2.5), we note that the dependence on the perturbation *ϵrn* is more complicated here, as the perturbation also affects the *shape* of the contour ∂*A*^{ϵ} on which the integral is defined. To evaluate these expressions, we need the concept of the ‘shape derivative’. Given a function and a perturbation vector field , such that **Z**|_{∂Ω}=**0**, the shape derivative *f*′ of *f* with respect to perturbation **Z** is defined through the identity [44]
2.10

In our problem, 2.11and in order to apply relation (2.10), we need to translate the vector operations into complex notation. This can be carried out with the help of the following lemma.

### Lemma 2.1

*For* *a complex-valued function,***Z***=[Z*_{1}*,Z*_{2}*] and Z=Z*_{1}+i*Z*_{2}*, the directional derivative of f in the direction* **Z** *can be expressed as*
2.12*where*
*when w=x*+i*y*.

We can then write for the r.h.s. of (2.5)
2.13where we assume **Z**=*r***n** in (2.10).

As concerns the shape derivative *f*′ in (2.13), since the reference velocity of the moving coordinate system *V* _{0}(*z*) does not depend on the perturbation, we have [*V* _{0}(*z*)]′≡0. For the integral term, we need to use the following shape-calculus formula [43], p. 354, eqn (4.17):
2.14where is a differentiable function (∂*A*^{ϵ}⊂*Ω*), *g*′ is its shape derivative and *κ* is the curvature of the contour ∂*A*=*Γ*. While in the standard texts of Delfour & Zolésio [43] and Haslinger & Mäkinen [44], formula (2.14) is derived for real-valued functions *g*, an extension to complex-valued integrals is straightforward. For us, then, cf. (2.1) and (2.6)
2.15

Hereafter, we will adopt the convention that the subscripts *z* and *w* on different symbols, e.g. *n*, *t* or *r*, will indicate where the corresponding quantity is evaluated. For the shape derivative *g*′ in (2.14), we note that *g* depends on the perturbation **Z** through the independent variables only, so that by (2.10) we have and therefore *g*′≡0. Finally, assembling (2.9) and (2.13), using (2.14), and performing the differentiations required in (2.13) and (2.14), we obtain, after rearranging terms, the following perturbation equation:
2.16describing the linear evolution of the local amplitude *r*(*τ*,*s*) of the vortex boundary perturbation in the normal direction (cf. (2.4)). In the above, we have also used the identity *t*=i*n* that follows from the counterclockwise orientation of ∂*A*. In combination with an initial condition,
2.17system (2.16)–(2.17) should be interpreted as an initial-value problem in which the solution *r*(*τ*,*s*) is subject to the periodic boundary condition *r*(*τ*,0)=*r*(*τ*,*L*), ∀*τ*≥0. We also add that some of the integrals in (2.16) are to be understood in the principal-value sense. We remark that, although some terms are present in both equations, equation (2.16) appears different from the perturbation equation obtained by Dritschel [13]. The assumption that, to the leading order in *ϵ*, the perturbations do not change the vortex area, i.e. , leads to a restriction on the admissible initial perturbations *r*_{0}. Shape differentiating the expression for the vortex area, we obtain
2.18which means that, in order to satisfy (at the initial instant of time) the constant-area condition, we need to restrict attention to initial perturbations *r*_{0} with the vanishing mean with respect to the arc length *s*. We add that, assuming constant vorticity *ω*=const., condition (2.18) also ensures that the vortex circulation *γ* remains unchanged by perturbations. Construction of perturbations satisfying condition (2.18) is discussed in §5 (cf. equation (5.2)). In the absence of the constant-area condition (2.18), our preceding analysis would need to be slightly modified, as then the shape derivative *g*′ would not vanish.

In the following two sections, we show how the stability analyses of the circular and elliptic vortices can be derived as special cases from general perturbation equation (2.16).

## 3. Stability of the circular vortex

In this section, we demonstrate how the well-known results concerning the stability of the Rankine vortex originally obtained by Kelvin [1] and reviewed by Lamb [2], Batchelor [3] and Saffman [4], in particular, the dispersion relation for the instability waves travelling along the vortex boundary, can be analytically deduced directly from general perturbation equation (2.16). Even though these results are classical, we go through these calculations in some detail because they will also play an important role in the stability calculations for vortices with arbitrary shapes. More specifically, as will be shown in §5 of this work, numerical stability computations for vortices with general shapes can be reduced to a ‘perturbation’ of the corresponding analysis for the circular (Rankine) vortex. The main advantage of this approach is that it allows us to compute all singular (principal-value) integrals analytically. For the Rankine vortex with radius *a*, we thus set *z*(*p*)=*a* e^{ip}, *n*(*p*)=e^{ip}, and *κ*=1/*a*, where *p*=(2*π*/*L*)*s*∈[0,2*π*]. For the moving coordinate system, because of the rotational invariance of the problem, its velocity *V* _{0}(*z*) can be selected arbitrarily, and without loss of generality, we choose *V* _{0}(*z*)≡0 resulting in in (2.16). The perturbation equation then simplifies to
3.1where we have also written *w*=*a* e^{iq}. In view of (2.18), and the periodicity of *r*, the first integral term on the r.h.s. in (3.1) vanishes for area-preserving perturbations. The last integral in (3.1) corresponds to the second to last in (2.16), and the last integral in (2.16) vanishes (this follows from the explicit formula for principal-value integrals that we give below). We note that relation (3.1) is equivalent to the perturbation equation obtained with the Lagrangian approach by Baker in his analysis of the stability of the Rankine vortex [5]. As in earlier studies [1–5], we will adopt the ‘normal mode’ approach and suppose that
3.2where and . As will be demonstrated below, the normal mode approach is justified because equation (3.1) does not lead to mode coupling. Using ansatz (3.2), the second integral in (3.1) becomes
3.3Each of the integrals on the r.h.s. in (3.3) can be written in the form for some . This principal-value integral can be evaluated analytically using contour integration (see appendix A), and the result is
3.4

Inserting ansatz (3.2) into perturbation equation (3.1), and then using (3.4) to reduce (3.3), we obtain an equation in which all terms are proportional either to e^{iks} or to e^{−iks}. The fact that there are no other terms justifies *a posteriori* the normal mode approach (3.2). Finally, isolating the terms proportional to e^{iks} and e^{−iks}, we obtain
3.5and its complex conjugate copy, so that the solution is *ρ*(*τ*)=*ρ*_{0} e^{(iω/2)(1−k)τ}, where is the amplitude of the initial perturbation. From (3.2), we obtain the following expression for the evolution of the perturbation
3.6in which the phase velocity *v*_{p}=(*ω*/2)((*k*−1)/*k*) agrees with Kelvin's solution [1–4]. We add that the problem analysed by Lamb [2] is in fact formulated in a bit different way because he assumed the perturbation in the form of a potential modification of the velocity field, so that the shape of the region *A* where *ω*≠0 remains unchanged [2]. It is therefore interesting that the solutions obtained in these two problems are identical, despite subtle differences in the underlying assumptions.

## 4. Stability of the elliptic vortex

In this section, we demonstrate how the well-known results concerning the linear stability analysis of Kirchhoff's elliptic vortex [6] obtained initially by Love [7] can be analytically deduced from general perturbation equation (2.16). Let the vortex boundary be described as
4.1where *a* and *b* are, respectively, the major and minor axes of the ellipse, and *p* coincides with one of the coordinates in the elliptic coordinate system. The key result is that the flow becomes linearly unstable when the ellipse aspect ratio *c*=*a*/*b*>3. As the aspect ratio increases, consecutive eigenmodes are destabilized. They are given in the form [7,9,10]
4.2where
4.3

As can be verified, representation (4.2) satisfies the area-preservation condition (2.18). The associated eigenvalues are
4.4so that eigenmodes (4.2) grow proportionally to e^{−iλktτ}. We note that only modes with *k*≥3 can become unstable, and to every such mode, there corresponds a decaying mode with . These results are deduced from perturbation equation (2.16) by plugging *r*(*τ*,*p*)=e^{−iλkτ}*r*_{k}(*p*) into this equation, where *r*_{k}(*p*) is given in (4.2), and then performing the following steps.

1. Substitute representation (4.1) for

*z*=*z*(*p*) and*w*=*w*(*q*) in (2.16) and convert the contour integrals to definite ones with*q*∈[0,2*π*] as the integration variable.2. Since

*p*and*q*appear only as arguments of trigonometric functions, use the following identities (): 4.5and likewise for and using to transform the resulting expressions to contour integrals over the unit circle |*u*|=1 with d*u*=i*u*d*p*.3. Perform partial fraction decomposition of the obtained integrand expressions and evaluate the integrals using the residue theorem for the poles

*u*_{0}inside the unit circle (|*u*_{0}|<1) and formula (3.4) for the poles*u*_{0}on the unit circle (|*u*_{0}|=1).

For any given value of *k*, operations described in steps 1–3 above can be performed symbolically using a symbolic algebra package such as Maple (the corresponding code is available as electronic supplementary material). Combining all resulting terms and noting that for the elliptic vortex
cf. (2.1), and , cf. (2.2), we recover expression (4.4) for the eigenvalue. We add that in the special case when *a*=*b*, the results of §3 are obtained, except that in this limit, the velocity characterizing the rotation of the coordinate system is *V* _{0t}=*ωa*/4.

## 5. Accurate numerical solution of stability equations

In this section, we describe a spectrally accurate numerical technique that allows us to solve perturbation equation (2.16) numerically for vortices of arbitrary smooth shape. A key feature of this approach is that the singular parts of the integrals in (2.16) are in fact evaluated analytically. We remark that such an approach was also used in other stability analyses of two-dimensional vortex flows, e.g. in Dritschel [13]. Our method is validated by comparing the results obtained numerically using different resolutions for the elliptic vortex with the analytical results derived by Love [7], and discussed in §4. First, we assume that the vortex boundary ∂*A* admits a smooth periodic and otherwise arbitrary parametrization *z*=*z*(*p*), *p*∈[0,2*π*]. In order to simplify the notation, we rewrite perturbation equation (2.16) as
5.1where the dependence of the perturbation *r* on time *τ* was suppressed for brevity, and , whereas , and represent, respectively, the three integral operators on the r.h.s. of (2.16). Quantities appearing in equation (5.1) that are related to the contour geometry, i.e. |d*z*/d*p*|^{−1}, *t*, *n* and *κ*, can be evaluated using spectral differentiation. Equation (5.1) is discretized by using the following ansatz for the perturbation:
5.2depending on *N* real-valued coefficients *α*_{k} and *β*_{k} and where is the growth rate we seek to determine. The factor |d*z*/d*p*|^{−1} in (5.2) ensures that area-preservation condition (2.18) is satisfied by all terms in the series, except for the constant term , which must be introduced in order to ensure the well-posedness of the collocation problem (it will be, however, subtracted via projection from the resulting eigenvectors). After using ansatz (5.2), the independent variable *p* is discretized using *N* equispaced collocation points (*N* is an even number). Denoting , this discretization leads to the following eigenvalue problem:
5.3where
5.4in which
5.5is the (invertible) collocation matrix,
5.6corresponds to the first term on the r.h.s in (5.1) with **D** denoting the spectral differentiation matrix associated with ansatz (5.2), and
5.7in which the first term on the r.h.s. is obtained via straightforward evaluation, whereas the calculation of the integral expressions , and is described below.

All integrals with bounded integrands are evaluated using the trapezoidal rule, which on a periodic domain is a spectrally accurate quadrature. Since *p* and *q* are now the independent variables, in this section, these characters will be used as subscripts to indicate where a given variable is evaluated. The integrand expressions in *V* _{1t}, cf. (2.1), and in
have removable singularities and are therefore bounded [4]. For the integral operators and , which are defined in the principal-value sense, we use a standard approach [47], in which they are split into principal-value integrals corresponding to a circle, which can be evaluated analytically using formula (3.4), and regular ‘corrections’, which can be approximated numerically with the spectral accuracy using the trapezoidal rule. In order to be able to use (3.4), the basis functions in ansatz (5.2) are expressed in terms of the complex exponentials, so that
and
and analogously for . Then, we have
5.8where . Expanding this expression in a Fourier series with respect to *q*, i.e. , the second integral on the r.h.s. in (5.8) can be evaluated as follows with the help of formula (3.4):
5.9

For the remaining integral operator, we have
5.10where . Expanding this expression multiplied by |(d*z*/d*q*)(*q*)|^{−1} in a Fourier series with respect to *q*, i.e. , the second integral on the r.h.s. in (5.10) can be evaluated as follows with the help of formula (3.4):
5.11

The first integrals on the r.h.s. in (5.8) and (5.10) are regular, and can be evaluated numerically (they also vanish in the circular case). Expressions and are computed in an analogous manner. After collecting expressions (5.4)–(5.11), the algebraic eigenvalue problem (5.3) can be solved using standard tools.

We close this section with a numerical validation of the proposed method. As the test problem, we use the case of the elliptic vortex, also studied analytically in §4. We compute the eigenvalues by solving problem (5.3) using different numerical resolutions *N* and compare the results against the closed-form expression (4.4) obtained by Love [7], see also Mitchell & Rossi [9] and Guo *et al*. [10]. More specifically, we analyse the imaginary parts ℑ(*λ*) of the eigenvalues responsible for the instability. In figure 2, we show the relative errors between the eigenvalues ℑ(*λ*^{N}) computed numerically with resolution *N* and ℑ(*λ*_{k}) obtained analytically by Love (cf. (4.4)) for two different aspects ratios of the elliptic vortex (for an aspect ratio equal to eight, there are in fact four distinct unstable eigenvalues). In figure 2, we note that this error decays exponentially fast, dropping to the machine precision level already for modest resolutions *N*, thereby confirming the spectral accuracy of the proposed method. We add that, as expected, the resolutions required to achieve a given accuracy increase with the aspect ratio of the vortex.

## 6. Conclusions

While perturbation equations similar to (2.16) have already been used for linear stability analysis of two-dimensional vortex patches [11,13,14], to the best of our knowledge, the present study offers a first complete derivation of this approach. In particular, it relies on methods of shape calculus that are a general and mathematically consistent way of dealing with the free-boundary aspect of the problem. We add that, in the context of vortex dynamics, such techniques have already been used to study continuation of families of solutions [41] and vortex control problems [48]. The proposed numerical approach is demonstrated to be spectrally accurate, reducing all singular integrals to closed form (3.4). The two validation tests presented offer ‘the proof of the concept’ for this approach.

Generalization of our method to problems involving several vortices is straightforward, and requires that an equation of the form (2.16) be written for each individual vortex with the interaction between the vortices captured by the field *U*(*z*) (which vanishes identically in the single-vortex examples considered in §§3 and 4). This description will be simplified by symmetries of the vortex configuration. Another interesting and far less researched application area for our approach is the stability analysis of three-dimensional axisymmetric vortex flows, with and without swirl, with compact vorticity support [49]. All these problems are left to future research.

Regarding limitations of the proposed approach, we remark that shape calculus in its standard formulation [50] is applicable only to problems posed on smooth manifolds, hence, our method would need to be modified so that it can be applied to contours with singularities, such as, e.g. corners (flows with such features typically arise as terminal members of families of steady solutions [41]). Some relevant ideas were mentioned by Elcrat & Miller [51]. Likewise, analysis of stability with respect to perturbations leading to topology changes requires the use of different methods.

## Acknowledgements

The authors are grateful to Prof. L. Rossi for insightful comments about Love's stability analysis of the elliptic vortex [7]. B.P. was partially supported through an NSERC (Canada) Discovery Grant.

## Appendix A. Evaluation of the integral

The calculations presented in §§3–5 require the values of the integral
A1for integer *M*. If *M*>0, let us consider
in which *C* is now the contour connecting the points 0, 2*π*, 2*π*+*Y* i, *Y* i, and 0 with a small semi-circular indentation into the upper half-plane above *p*, where *p* ∈(0,2*π*) and *Y* >0. This contour integral vanishes by Cauchy's theorem. The integrals over the left and right lateral sides of *C* cancel by periodicity. For the top segment of contour *C*, writing *q*=*x*+i*Y* , we have
so for , the integral over that segment goes to zero. It remains to take the limit of the integral over the bottom segment of the contour *C* as the radius of the semicircular indentation vanishes. This is the original principal-value integral (A1) augmented by the contribution from the indentation, which can be evaluated using the Cauchy integral formula. We therefore write the integrand as
and find, using L'Hôpital's rule,

Thus, the limit of the integral over the indentation is this value multiplied by −*πi*, and, finally, we obtain
A2

When *M*<1, instead of contour *C*, we use its reflection into the lower half-plane. Then, since *Y* is negative, the integral over the bottom segment goes to zero as and, since the small semicircle is positively oriented,
A3

## Footnotes

- Received December 4, 2012.
- Accepted January 11, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.