Multipole Vortex Blobs (MVB): Symplectic Geometry and Dynamics

Vortex blob methods are typically characterized by a regularization length scale, below which the dynamics are trivial for isolated blobs. In this article, we observe that the dynamics need not be trivial if one is willing to consider distributional derivatives of Dirac delta functionals as valid vorticity distributions. More specifically, a new singular vortex theory is presented for regularized Euler fluid equations of ideal incompressible flow in the plane. We determine the conditions under which such regularized Euler fluid equations may admit vorticity singularities which are stronger than delta functions, e.g., derivatives of delta functions. We also describe the symplectic geometry associated with these augmented vortex structures, and we characterize the dynamics as Hamiltonian. Applications to the design of numerical methods similar to vortex blob methods are also discussed. Such findings illuminate the rich dynamics which occur below the regularization length scale and enlighten our perspective on the potential for regularized fluid models to capture multiscale phenomena.


Introduction
Vortices are important in hydrodynamics because they are the sources for the incompressible flow field. The vorticity distribution at any instant of time determines both the current state of the flow and its future evolution, for given boundary conditions. This property holds for any Hamiltonian system, and it can indeed be shown that the dynamics of vortices can be usefully expressed in Hamiltonian form. In the vorticity and stream function formulation of an ideal incompressible planar fluid, the evolution of the vorticity distribution ω(x, y, t) is given by where ω = − ψ is the vorticity, ψ is the stream function, and = ∂ x x + ∂ yy is the Laplace operator. The corresponding (x, y) components of the Eulerian velocity field are given by If one is willing to view the vorticity ω as a distribution, one can consider point vortex solutions. In particular, point vortices are obtained if one considers the vorticity solution ansatz where i (t) ∈ R, z = (x, y) ∈ R 2 and δ z i (t) is the Dirac delta distribution centered at the point z i (t) = (x i (t), y i (t)) ∈ R 2 at a given time t ∈ R. Substitution of this ansatz into (1) yields the following well-known finite dimensional system in the form of Hamilton's canonical equations, where G(z) = −(2π) −1 ln( z ) is the Green's function for the planar Laplacian. A point vortex approximation to a continuous distribution of vorticity for Euler's fluid equations is problematic, though a point vortex induces a flow velocity which becomes unbounded. However, when the point vortex is made smooth and bounded (regularized), the approximation becomes reasonable (Chorin 1973).
For example, one may consider the regularized form of the vorticity equation given by choosing a translationally and rotationally invariant smoothing kernel K δ of width δ > 0 and defining the regularized vorticity as K δ * ω = − ψ while continuing to use (1) to evolve ω in time. For example, K δ (z) = exp(− z 2 /δ 2 ) is considered in Beale and Majda (1985). In this case, the point vortex ansatz yields (2) again, except that the singular Green's function G is replaced by the smooth kernel where Ei(·) denotes the exponential integral function. The vorticity kernel G δ has no singularity at the origin for δ > 0 and is known as a vortex blob. This system is the starting point for the vortex blob method, introduced in Chorin (1973) (albeit with a different regularization). The economy of the vortex blob method derives from the property that Dirac delta distributions are hyper-local (i.e., parametrized by position), and the property that the vorticity equation (1) admits Dirac delta distributions as solutions. However, there are many distributions which are localized to a similar degree (e.g., derivatives of delta functions, ∂ x δ z i ).
In this paper, we study the more general vorticity solution ansatz, We find that this ansatz yields a closed finite dimensional system which generalizes vortex blobs. We call these new carriers of vorticity multipole vortex blobs or MVBs.

Main Contributions
(1) Section 2 briefly reviews the background for vortex methods in fluid modeling.
(2) Section 3 reviews the relationship between regularized fluids and vortex blob methods.
(3) Section 4 derives the equations of motion for point vortices and MVBs as exact solutions of a regularized vorticity equation. (4) Section 5 derives the conservation laws for these equations, such as energy, linear momentum, and angular momentum, and circulation. The derivation of these conserved quantities as symplectic momentum maps can be found in Appendix B. (5) Section 6 explains the relationship between the dynamical systems for MVBs and an implicitly defined closed dynamical system which governs the spatial moments of the vorticity distribution. (6) Section 7 discusses numerical aspects of using MVBs to model fluid dynamics, such as approximations of initial conditions (Sect. 7.2), and grouping of computational nodes (Sect. 7.1). (7) Section 8 presents the results of several numerical experiments involving small numbers of vortices, for N = 1, 2, and 3. (8) MVB dynamics are Hamiltonian. We present the symplectic and Hamiltonian structure of MVB dynamics in Sect. 9.

Background
Vortex methods for fluid modeling predate the computer age, and references to them can be found in the work of Helmholtz (Smith 2011, see the introductory section). For example, the use of point vortices as idealized solutions can already be found in a 1931 paper concerning a "line of discontinuity" in planar fluid flow (Rosenhead 1931). At the beginning of their development, the infinite velocities (and energies) associated with point vortices caused great difficulties, both numerically and theoretically. In fact, the point vortex approach did not produce a competitive numerical method until the 1970s, when the problems related to singularities were overcome by regularizing the singular vortex kernel to form a vortex blob. Stochastic perturbations were further included to model viscosity (Chorin 1973). These adjustments to the classical point vortex method yielded the vortex blob method, which quickly became of practical use for realistic fluid flow modeling. In particular, the regularized system proved more amenable to error analysis. It was shown that the solutions of the vortex blob method converge to solutions of the Navier-Stokes equations in Hald (1979). Later, stronger convergence rates were achieved by judicious choice of vortex kernels. By convolving the singular vortex kernel with sums of Gaussian smoothing kernels, a sequence of vortex blob kernels with faster convergence rates was found. Specifically, the convergence rate of the mth kernel was found to be of order h mq for any q ∈ (0, 1) where h = δ q is a grid-spacing parameter and δ > 0 is a length scale associated with the smoothing kernel Majda 1982, 1985). Simultaneously, the symplectic geometry of point vortices was clarified in Marsden and Weinstein (1983) by invoking Arnold's interpretation of ideal fluids (Arnold 1966). The findings of Marsden and Weinstein (1983) were developed further in Gay-Balmaz and Vizman (2012) to handle fluid flow on manifolds with nontrivial homology. While this theoretical development clarified the geometry of point vortices, vortex blobs were sometimes thought to be a numerical "trick" which violated the geometric interpretation. However, this thought was banished with the invention of the Euler-α model, a regularized model of ideal fluids with a parameter α representing the typical correlation length of fluctuations away from the mean of a Lagrangian fluid path (Foias et al. 2001). In particular, vortex blob solutions associated with a specific kernel serve as exact solutions to the Euler-α model (Oliver and Shkoller 2001). The Euler-α kernel is different from the kernels used in Chorin (1973) and Beale and Majda (1985). A comparison of the Euler-α kernel to the m = 1 kernel of Beale and Majda (1985) is given in Holm et al. (2006) for vortex filament and vortex sheet motion.
While vortex blobs performed well, they did not capture all of the qualitative richness observed in fluid vorticity dynamics. In particular, blobs of vorticity in real ideal fluids are known to change shape and deviate from initially circular distributions. A numerical method is proposed in Rossi (1997Rossi ( , 2005 to capture these shape dynamics by adding basis functions with nontrivial moments in the study of vortex merger (see, for example, Melander et al. 1988;Le Dizès and Verga 2002;Meunier et al. 2005). Another distinct model obtained by projection onto a Hermite basis is described in Nagem et al. (2009). This projection yielded a finite dimensional systems which modeled the (truncated) moments of the vorticity of an ideal incompressible fluid. The derivation of simplified combinatorial formulas invoked by the dynamics of this model was discovered in Uminsky et al. (2012), and these formulas have made the method numerically tractable for a large number of moments.
A dual approach to the moment-based methods of the previous paragraph (Rossi 1997(Rossi , 2005Nagem et al. 2009) is to consider multipole-based methods. This is the approach proposed in Nicolaides (1986), where an initial vortex ansatz consisting of sums of distributional derivatives of dirac delta distributions is considered. Such an idea has occurred intermittently in various forms in the literature, over many years.
For example, a regularized vortex blob model, in the spirit of Majda (1982, 1985), which considered vorticity distributions of the form ω Chiu and Nicolaides (1988). Here it was proven that this augmentation of the traditional vortex method will yield faster spectral convergence than that of a traditional vortex blob method. The current article considers higher order derivatives and can be seen as a natural follow-up to (Chiu and Nicolaides 1988). More recently, dynamics have been derived for interactions of pure vortices and pure dipoles. These come from vorticity distributions of the form ω i with the assumption that the locations of the dipoles, and the vortices never overlap and that their self-interaction terms may be neglected (Yanovsky et al. 2009;Tur et al. 2011). In a different approach, approximations of dipoles are created by holonomically constraining vortices of opposite strength to be a fixed distance from one another (Tchieu et al. 2012). The question remains, however, to what extent the dynamics of Tchieu et al. (2012) approximates those of Yanovsky et al. (2009) andTur et al. (2011) after self-interaction terms have been neglected. In summary, the removal of self-interaction terms is one of the primary obstacles to obtaining a multipole-based generalization of the point vortex method (Smith 2011). Moreover, the spectral error decay rates found in Hald (1979), Majda (1982, 1985) and Chiu and Nicolaides (1988) arise from the use of vortex blobs in place of (singular) point vortices. In this article, we will follow this regularization-based approach.

Vortex Blobs and Regularized Fluid Models
In this section, we review a class of regularized fluid models and their relationship with vortex blob methods (for a more detailed discussion see Holm et al. 2006). The sort of fluid models we consider take the form where L α is a SE(2) invariant linear psuedo-differential operator with a length scale parameter α > 0 such that lim α→0 Q op = 1. When L α is just the identity, the above "model" is Euler's equations of motion for an ideal fluid.
where is the Laplace operator, then we obtain the the Euler-α model, the solutions of which will converge to solutions of Euler's ideal fluid equations as α > 0 vanishes (Foias et al. 2001).
We may replace u with its stream function, ψ, in order to rewrite the above equations as This allows us to represent planar fluid dynamics in terms of scalar functions and distributions rather than vector fields.
The relationship between these regularized models and vortex blobs methods comes from first considering the point vortex ansatz If the operator, • L α , has a non-singular Green's function, G α , then substituting the ansatz into (5) implies that We should note that when L α is the identity (i.e., for an Euler fluid), then G α is singular, and an extra argument (perhaps a physical one) must be presented in order to allow the resulting singular velocity fields. In this paper, no such issues with singularity arise because we are modeling an Euler fluid with a regularized fluid where L α has a non-singular Green's function. Substitution of ψ into (4) then implies the following equations of motion for the vortex cores z i = (x i , y i ) and the strengths i (t): When α = 0 and L α = 1, this is nothing but the point vortex method. When α > 0, it is possible for ψ to be much more regular, and we obtain various vortex blob methods.
In particular, we obtain the smooth vortex blobs of Beale and Majda (1985). It is notable that (9) and (7) together form a finite dimensional ODE. The solutions of this ODE are exact solutions to the regularized fluid model. Again, this is valuable because the solutions of many regularized fluid models have been shown to converge to solutions of the ideal fluid equations as α vanishes. This paper seeks to generalize these point-like solutions to regularized fluid models to obtain a richer class of solutions with richer conservation properties.

Equations of Motion
In this section, we derive the equations of motion for the time-dependent parameters which specify multipole vortex blobs (MVBs). The zeroth-order MVBs are just standard vortex blobs, and the resulting equations of motion are those of the standard (non-stochastic) vortex blob algorithm (Chorin 1973). The first-order MVBs are regularized dipoles, and the equations of motion are those of Chiu and Nicolaides (1988). Here we will derive the equations of motion for N th-order MVBs following the approach of Chiu and Nicolaides (1988).
Consider the ansatz for the vorticity, . This form of the kernel produces one of the vortex blobs presented in Beale and Majda (1985), and the resulting numerical method is identical for spatially constant dynamical variables mn i (t) ∈ R for i ∈ S where S is some countable set. The stream function is The corresponding velocity field is given by Examples of the types of velocity fields produced are depicted in Figs. 1 through 3 on page 7. We seek equations of motion for the mn i (t)'s and z i (t)'s such that the velocity field (10) satisfies the vorticity equation (1). In the following calculations, we will not show the explicit time dependence of the dynamical variables.
We now find Fig. 2 The flow field around a first-order MVB with = 0, x = 1, y = 1 is that of a regularized dipole By invoking (24) of Appendix 1, we can rearrange the previous equation to obtain Similarly, we find Substitution of these expressions into (1) yields the vanishing of a linear combination of the distributions ∂ m x ∂ n y δ z i for m + n ≤ N + 1. Since each of these distributions is linearly independent of the others (assuming the z i 's are distinct), their individual coefficients must each vanish independently. If we focus on the terms of the sum where m + n = N , we obtain coefficients for ∂ x δ and ∂ y δ at the core locations. The vanishing of these coefficients yields the dynamics for MVB cores The vanishing of the coefficient of δ z i yields We observe that d k i /dt depends on ψ at the vortex cores z i , and the vortex core dynamics depend on ψ as well. Fortunately, we already found that ψ is purely a function of k i and z i , as stated in (9). Thus, (12) and (11) (with (9)) form a closed finite dimensional system. Most notably, by construction the vorticity equation (1) admits the N th-order MVB ansatz for the vorticity in (8) as a solution when the z i (t)'s and the i (t)'s satisfy the just derived finite dimensional system.

Remark 4.1
In the point vortex method (i.e., the unregularized case where L α = 1), this derivation of the dynamics requires an extra step. In particular, one must discard the self-interaction term, which we will describe here. For point vortices, ω = − ψ. Substituting the point vortex ansatz ω = i i δ(z − z i (t)) into the equations of motion (4) would then yield the nonsensical equatioṅ We say "non-sensical," because the right-hand side explodes when you evaluate the ith term in the sum, the self-interaction term. Historically, it is customary to discard this self-interaction term based on physical and symmetry principles (Marchioro and Pulvirenti 1994, Chapter 4). In contrast, for blob methods the logarithmic kernel is replaced with a differentiable kernel function, such as a Gaussian. This allows one to retain the self-interaction terms. In the case of standard vortex blobs (i.e., zeroth-order MVBs), this distinction makes no difference because the gradient of the kernel vanishes at the origin and the self-interaction term contributes nothing to the dynamics. However, the derivatives of the kernel of degree 2 and higher do not vanish at the origin. As a result, the self-interaction terms do contribute to the dynamics for MVBs of order 2 and higher. The choice to discard the self-interaction terms in Yanovsky et al. (2009), versus our choice to keep them explains one of the major discrepancies between our work and Yanovsky et al. (2009). In particular, Yanovsky et al. (2009) was concerned with generalizing the (unregularized) point vortex method in the same way that we have generalized the vortex blob method. Once the ansatz ω = mn i ∂ m x ∂ n y δ(z − z i ) was substituted into the equations of motion, they discarded the self-interaction terms in order to handle the singularities in the velocity field. They had no other choice. Except for the initial regularization step we took, this discarding of the self-interaction term is the primary place where the derivation of the equations of motion presented here diverges from the derivation in Yanovsky et al. (2009). Discarding the self-interaction term in Yanovsky et al. (2009) leads to contradictory compatibility equations for singularities of degree 2 and higher. This is one regime where the self-interaction terms have an impact on the dynamics in our regularized formulation. One of the major findings of Yanovsky et al. (2009) was that one could avoid these contradictory compatibility conditions by limiting one's self to combinations of point vortices and dipoles. Even in this limited scenario, our equations of motion do not match even in a regularized sense, as the vortex cores of the dipoles are not advected by the (singular) velocity field in Yanovsky et al. (2009). Additionally, as the regularization parameter goes to 0 in our framework, the velocity fields become singular, and the equations of motion for the 's will explode. So we cannot expect to observe any form of convergence to the finite valued ODEs of Yanovsky et al. (2009).

Conserved Quantities
In this section, we begin to touch upon some of the symplectic geometry of MVBs. To begin, let us consider a general vorticity distribution ω. The energy is defined as where ψ = G α * ω. The vorticity equation, (1), can be seen as an instance of Hamilton's equations on a Poisson manifold. In this case, the Poisson manifold is the space of vorticity distributions, and the Poisson bracket is the vorticity Poisson bracket derived in Marsden and Weinstein (1983). As the Hamiltonian is conserved by Hamilton's equations, we should expect H (ω) to be constant in time. Indeed, we find that if ω satisfies (1), then By integration by parts, we can remove the partial derivates from the ω's to find which vanished by the equivalence of mixed partials. As (1) is a Hamiltonian system, we can consider searching for symmetries to find other conserved quantities using Noether's theorem. We've relegated the discussion of the relevant symplectic structure to Appendix 2, where derivations and proofs of the following can be found. Here we can summarize the appendix.
It's simple to observe that the Hamiltonian H is translation invariant, and that H is rotationally invariant as long as the kernel G α has rotational symmetry. Thus, we should expect there to be conserved quantities tied to these symmetries. We find that the quantities are conserved. The relationship between these quantities and the symmetries of the system is explained in Appendix 2. Alternatively, one can observe the conservation of these quantities by direct calculation in the same way that conservation of energy was verified.
As the MVB ansatz is consistent with (1), we can substitute the MVB ansatz into the above conserved quantities, to obtain conserved quantities for the MVB evolution, (11) and (12). We obtain the following conserved quantities: Again, the first two quantities, J ang and J lin , are momenta derived from Noether's theorem for the rotational and translational symmetries of the fluid. The quantity H is the kinetic energy of the fluid. In Sect. 9, we will characterize the MVB dynamics as Hamiltonian systems, with Hamiltonian H . To each individual MVB, there are numerous conserved quantities which can be seen as a manifestation of the conservation of circulation. To show this, let u = (u, v) = (∂ y ψ, −∂ x ψ) be a time-dependent vector field which satisfies (1). The flow of u is the diffeomorphism, t : R 2 → R 2 , which sends particle labels at time 0 to their positions at time t. If ω t is the vorticity at time t, then ω t ( t (z)) = ω 0 is constant in time. This conservation law can be seen as a corollary of Kelvin's circulation theorem (Arnold and Khesin 1992). As a consequence, the quantity . By applying the change of variables formula and invoking the incompressibility condition, det(D ) = 1, we find This form of writing J (t) makes sense when ω t is a distribution. As a result, we find that for a vorticity of the form (8), the quantity is conserved for any f ∈ C ∞ 0 (R 2 ). While this conservation law holds for all functions with compact support, f , we do not obtain infinitely many conserved quantities when ω t satisfies the MVB ansatz and S is finite. This is because the expression on the righthand side only depends on the N th-order Taylor expansion of f at z i (0) ≡ −1 t (z i (t)), as is illustrated by the Faà di Bruno formula. We will not display the Faà di Bruno formula here because it requires nearly a page of notational definitions before to writing it down (Constantine and Savits 1996). Nonetheless, by computing the cardinality of jet spaces, one would obtain card(S) N (N +1)

Moments
In this section, we present how the moments of the vorticity distribution evolve in time. We will find that when the vorticity distribution is that of a MVB, and then the moments form a closed dynamical system at finite order.
The (a, b) th moment of the vorticity, ω, centered around the vortex position (x i , y i ) ∈ R 2 is given by We call the integer a +b the order of the moment. For a general vorticity, the evolution for the nth-order moments will depend on the (n + 1)th and higher order moments and so we cannot concoct a closed dynamical system for the moments of order n and less. However, this is not the case if ω satisfies the MVB ansatz, and the points (x i , y i ) are given by the locations of the jet vortices. If ω satisfies the MVB ansatz (8), then for a + b ≤ N and i ∈ S. Given the points z i ∈ S, one can write the moments in terms of the circulation strengths, the 's. For the moment μ ab i with a + b ≤ N with a, b ∈ N, we may invert this relationship to write mn i = mn i (μ), i.e., as a function of the moments. Invoking the motion equations for the 's and substituting, the relation between the 's and the μ's yields a closed dynamical system for the μ's.
Remark 6.1 This relation between the 's and the μ's may also be important in the context of plasma physics, especially when one recalls that (1) can be interpreted as a one-dimensional plasma model. Specifically, phase-space moments of the Vlasov probability distribution form an important dynamical link between Lagrangian-particle and Eulerian-continuum descriptions. The phase-space moments of the Vlasov probability distribution provide collective coordinates for the Hamiltonian dynamics of ensembles of particles. For more explanation of this property of Hamiltonian collectivization of the phase-space moments, see Guillemin and Sternberg (1990), Holm et al. (1990) and Gibbons et al. (2008a, b). In plasma dynamics, the phase-space moments arise from a Taylor expansion of the Vlasov particle distribution, taken around its centroid in phase-space. For planar incompressible flow of an ideal fluid, the phase-space comprises the (x, y) Lagrangian coordinates of a fluid particle, and the corresponding moments arise from Taylor expansions around the centroid of the (smooth) vorticity distribution. The duality between the resulting spatial moments of a smooth vorticity distribution and the MVBs corresponding to higher order singular vorticity distributions also obtained from a Taylor expansion raises the intriguing question of finding a relation between these two types of dynamical description. This question is particularly intriguing because the dynamics of moments beyond quadratic order in general does not close to form a finite dimensional Hamiltonian system, while the dynamics of MVBs closes at every order.
Remark 6.2 There exist other systems for approximating the dynamics of moments which differ from the one presented here. In particular, the equations of motion for the moments here form a closed system at order N , whereas other methods for deriving dynamical systems for moments (Uminsky et al. 2012;Nagem et al. 2009;Gibbons et al. 2008a, b) require truncations in order to form a closed system. For example, Uminsky et al. (2012) approximates the stream function as a sum of Hermite functions with evolving centroids and weights. In order to obtain the evolution for the weights and the centroids, they project the viscous vorticity equation onto this space via L 2 projection. The resulting formulas are explicit and efficient to compute, albeit more complex than the formulas found in this paper. The primary source of error for Uminsky et al. (2012) over long times is the discrepancy between the projected evolution equations and the true evolution equations. In contrast, we approximate an Euler fluid with a regularized fluid equation which we solve exactly. This is not to say that error is not accumulated in time. The primary source of error for our method over long times is the discrepancy between the regularized fluid equations and the true fluid equations.
Admittedly, the equations of motion for the moments in Uminsky et al. (2012) bear some resemblance to the equations of motion for the 's in our method. Both are quadratic in their respective variables, with coefficients involving combinatorial functions. A more precise relationship, if one exists, is difficult to discern. Philosophically, the methods share much in common. However, due to the fundamental approximation technique of projecting the equations of motion versus regularizing them, the methods are indeed distinct. This difference cascades throughout the study of both methods. For example, the convergence for Uminsky et al. (2012) is obtained via the convergence of spectral approximations, while the convergence of our method is a corollary of the convergence of a regularized fluid model (see Mumford and Michor 2013;Foias et al. 2001 for such convergence proofs).

Numerical Aspects
In this section, we discuss various numerical aspects of using MVBs to model fluid dynamics. We will observe how MVBs can be used to reduce the number of necessary pairwise computations without a drastic compromise in accuracy. We will also present an algorithm for constructing an initial condition of MVBs from a given stream function.
Remark 7.1 We refer to Chiu and Nicolaides (1988) for a convergence proof and error analysis of the first-order case. A convergence proof is beyond the scope of this article. Suffice it to say, such a proof would likely resemble Chiu and Nicolaides (1988).

Grouping and Reduction of Pairwise Computations
Let us consider the vorticity distribution ω = 1 δ z 1 + 2 δ z 2 .
If z 1 and z 2 are close, we can define the quantitiesz = (z 1 + z 2 )/2 and δz = z 1 − z 2 to obtain the approximation where h = δz . Therefore, the distributioñ serves as a o(h) approximation of ω in the sense of distributions. Moreover, the stream functionψ := G δ * ω is an o(h) approximation of ψ := G δ * ω in the traditional sense of analysis on functions.
We have just described the first case of grouping two N th-order MVBs concentrated at z 1 and z 2 into a single (N + 1)th-order MVB concentrated at the average position z. More generally, we can consider the ansatz The above computation implies that the quantitỹ serves as an o(h) approximation of ω. Of course, this again implies that the corresponding stream functions are approximated to order h as well. Note thatω is concentrated above a single point,z, while ω is concentrated above two points.
Remark 7.2 Such reductions are even more dramatic when considering higher order jets. In particular, 2 N zeroth-order MVBs can be approximated with a single N th-order MVB by applying the above approximations iteratively. The computation of pairwise interactions in the vortex method was once a major bottleneck in implementing the standard vortex method for real-world applications. It was not until the invention of the fast multipole method that it became tractable to compute millions of pairwise interactions by reducing the complexity from an O(n 2 ) calculation to an O(n log(n)) calculation, where n is the number of vortices (Greengard and Rokhlin 1987). However, in the case of viscous fluids with boundaries, vorticity is shed from the boundaries. As a result, the vortex blob method of Chorin (1973) created new vortices at the boundary by using the Kutta condition as a creation criteria. For these applications, n will grow in time without bound, and some means of discarding vortices must be invoked. It is here that the grouping of MVBs could be useful. If one merges two N th-order MVBs to obtained a (N + 1)th-order MVB, the amount of scalars and data typically increases. So one must still make a tough decision as to what data to discard (e.g., through some tolerance or by simply truncating at level M). Nonetheless, the analysis presented here could shed light on how best to implement this approach.

Remark 7.3
The merging of blobs of vorticity has been studied analytically (Melander et al. 1988) and numerically (Weiss and McWilliams 1993;Melander et al. 1988;Le Dizès and Verga 2002), as well as in the laboratory (Fine et al. 1991). All of this study has been in the slightly viscous (or nearly inviscid) regime. The grouping approach discussed here can be used to numerically resolve such collision events. In theory, there is no issue with collisions because we are considering regularized vortices where the induced velocity field from a single MVB is always finite. However, as δ becomes smaller, the velocity near the vortex core diverges. This should be of concern as the convergence analysis of the vortex method pre-supposes that δ 1. Typically such a near collision is handled by using a smaller time step (as the ODE is quite stiff). Grouping of MVBs suggests an alternative by avoiding this pairwise interaction altogether. Perhaps such an approach could be viewed as a variation of the punctuated dissipation events described in Weiss and McWilliams (1993) where an initial vorticity distribution is found to asymptotically approach a smoother axisymmetric vortex blob, and discrete vortex mergers are implemented to model this behavior.
Remark 7.4 There are qualitative questions which arise from mergers. For example, when two zeroth-order vortex blobs are near each other, they will typically scatter after some finite time. Merging these blobs into a single first-order blob will prohibit this scattering event from ever occurring. That both the zeroth-order MVB solution and the merged first-order MVB represent exact solutions of the fluid (after the merger event) are attributable to the long-term sensitivity to initial conditions near collision events. The scattering angle can be virtually anything since zeroth-order MVBs can waltz around each other many times before scattering. The amount of time two zerothorder MVBs can spend waltzing around each other, and perhaps the merged solution represent some sort of limiting solution. That is to say, the merged solutions can be interpreted as the "waltzing for eternity" solution.
The irreversibility of merging is disturbing when one takes it to its extreme, one massive high order MVB. In order to address this, a means of splitting high order MVBs into lower order ones should be considered. The primary difficulty here is in determining when to split. In the case of mergers, we can decide to merge MVBs when they are close. Such a criterion is not immediately apparent in the case of splitting MVBs.

A Numerical Experiment with Grouping
For illustrative purposes, we can numerically group four zeroth-order MVBs into two first-order MVBs, and then one second-order MVB. In particular, we can consider the initial condition The corresponding dynamics are depicted in the top row in Fig. 4. Next we group z 1 with z 0 and z 2 with z 3 in order to obtain two first-order MVBs with initial condition The corresponding dynamics are depicted in the middle row in Fig. 4. The dynamics appear qualitatively similar at the beginning of the evolution. Then the dynamics diverge around time t = 150 when the two first-order MVBs separate from one another, in contrast to the dynamics of the zeroth-order MVBs. Finally, we group the two first-order MVBs to obtain a single second-order MVB. Again, the dynamics appear qualitatively similar at the beginning of the evolution. Oddly, the dynamics of the second-order MVB appear qualitatively similar to the zeroth-order case even at t = 253. As there is only a single vortex, the separation of vortices mentioned in the first-order MVB experiment is not possible here. As a result, the dynamics of the original zeroth-order MVB dynamics appear to be approximately recovered.

Approximation of Initial Conditions
In this section, we will illustrate how initialize MVBs when given a stream function ψ at time 0. We can begin by defining an inner product on the space of distributions on R 2 , given by Consequently, the energy of the fluid is given by H (ω) = 1 2 ω 2 G δ = 1 2 ω, ω G δ . Let K be a compact set, and let 0 < h 1 be small so that we may define the finite grid h = {(ah, bh) ∈ K | (a, b) ∈ Z 2 }. 1 Given an ω ∈ D (R 2 ), we can attempt to approximate ω via Dirac deltas supported on hZ 2 . There is a natural way to do this with respect to the inner product ·, · G δ . We could define ω can be seen as a zeroth-order approximation to ψ = G δ * ω because ψ (0) h (z) = ψ(z) for all z ∈ h . Therefore, for smooth ω's, we obtain an error of order O( x) for a grid spacing of x using zeroth-order MVBs. The same reasoning applies if we consider ω for ψ = G δ * ω, z i ∈ h , and |β| ≤ k. Then ψ (k) h (z) = i,α (−1) m+n mn k ∂ m x ∂ n y G δ (z− z i ) serves as an order k approximation of ψ when ψ ∈ C k . In particular, for smooth ω's, we obtain an error of order O( x k+1 ) for a grid spacing of x using kth-order MVBs.
As an example, we numerically compute the corresponding approximations of the stream function The results are depicted in Fig. 5 where we observe sup-norm convergence on the interior of K . In particular, we measure the sup-norm error on the subregion (−3 < x, y < 3) with K = {(x, y) | −6 ≤ x, y ≤ 6}. We observe convergence using MVBs at orders zero, one, and two. In each case, a grid spacing is reached where the error plateaus (possibly due to machine precision). Nonetheless, higher order MVBs appear to out perform lower order ones for smaller grid spacings. In particular, we observe slopes in a log-log plot of magnitudes 1,2, and 3, suggesting that first-, second-, and third-order convergence rates for zeroth-, first-, and second-order MVBs, respectively. In terms of complexity, in order to achieve a desired error bound, e tol > 0, you would need to use a grid with O(e as k increases, one could object that a high order MVB is much more complex than a low order one. However, the number of degrees of freedom for a kth-order MVB is 2 + k j=0 (2 k /k!) which monotonically converges to a constant (roughly 9.39) as k → ∞. Therefore, the number of degrees of freedom is dominated by O(e −2/(k+1) tol ) as well. In other words, when ψ is highly differentiable we observe benefits in terms of complexity and storage to using a larger k regardless of whether one measures complexity by the number of parameters to keep track of, or the number of MVBs.

Numerical Experiments
In these section, we present the results of numerical experiments involving small numbers of vortices, for N = 1, 2, and 3.

Behavior of Isolated MVBs
Next, we will briefly explore the behavior of a single isolated kth-order MVB with mn = 0 with m + n < k for k = 0, 1, 2. This case allows us to investigate the dynamics induced by the higher order circulation variables in the absence of the lower order ones.

Order 0
The behavior of a single zeroth-order MVB is explicitly solvable because the dynamics are stationary.

Order 2
The behavior of a second-order vortex does not seem to be explicitly solvable. Here we consider initial conditions for which and all the other circulation variables are initially set to 0. The results are depicted in Fig. 7. We observe a structure which rigidly rotates counterclockwise.

A Scattering Experiment
Next we consider two MVBs. The first is a first-order MVB with an initial velocity pointed just slightly above origin. The second MVB is a standard zeroth-order vortex located at the origin. Specifically, we consider the initial conditions with mn i = 0 for m +n > 1 and i = 0, 1. The vortex at the origin appears to remain at the origin throughout the numerical run (t = 0 to t = 150). The first-order MVB starts by moving to the left in a straightline until it comes into proximity of the zeroth-order vortex. Then the first-order MVB swings around the zeroth-order vortex, traversing an angle of roughly 30 degrees before zooming off into the lower left quadrant of the plane in a straight line. These results are depicted in Fig. 8.  Fig. 8 A numerical run is shown with mirror-image initial conditions for two first-order MVBs, as given in (19). From left to right and top to bottom, these are snapshots at times t = 0, 25, 50, 75, 100, 125, respectively,

The Method of Images
Here we incorporate first-order MVBs into the method of images (Jackson 1975;Smith 2011). We consider the initial conditions consisting of two first-order MVBs which are mirror images of each other with respect to the x-axis. By symmetry, the resulting vector field should be tangential to the x-axis and provides a means of considering a boundary that satisfied the no-penetration condition. Specifically, we consider the initial condition: with mn i = 0 for m + n > 1 and i = 0, 1. The resulting dynamics depicted in Fig. 9 shows that as a first-order MVB approaches a boundary, it will turn its motion along the boundary and then move away so that its angle of reflection equals its angle of incidence.

Hamiltonians and Symplectic Structures
In modern Hamiltonian mechanics, as described in Abraham and Marsden (1978) and Arnold (2000), the Hamiltonian is a function on a symplectic manifold, which produces equations of motion. An important instance of a symplectic manifold is a Fig. 9 Numerical results are shown for two first-order MVBs with mirror-image initial conditions given by (20). From left to right and top to bottom, these are snapshots at times t = 0, 1.7, 3.4, 5, 6.7, 8.4, respectively. Apparently, a first-order MVB reflects elastically from a fixed boundary, so that its angle of reflection equals its angle of incidence coadjoint orbit (defined below). In this section, we compute the coadjoint orbit of a MVB as well as the associated symplectic structure. The coadjoint orbit of an initial vorticity distribution ω 0 comprises the set In fact, Orb(ω 0 ) inherits the structure of a smooth manifold, and a tangent vector on Orb(ω 0 ) at the pointω ∈ Orb(ω) is given by a distribution of the form £ u [ω] := u∂ xω + v∂ yω for some (nonunique) divergence-free vector field u = (u, v) ∈ X div (R 2 ). The symplectic structure is nothing more than a special case of the one derived via the Kirillov-Kostant-Souriau theorem (Abraham and Marsden 1978, see the boxed formula on p.303). In particular, the symplectic structure on Orb(ω) is given by Here we have used the change of variables formula and the fact that det(Dϕ) = 1. By the multivariate Faá di Bruno formula, the expression ∂ α | z=Z i ( f •ϕ)(z) is a sum of the partial derivatives of f at the points ϕ(Z i ) of order less than that of the multi-index α (Constantine and Savits 1996). Thus, ω 0 •ϕ −1 is contained in the finitely parametrized subset M (k) := { |α|≤k α i ∂ α δ Z i } for any ϕ ∈ SDiff(R 2 ). Therefore, Orb(ω 0 ) is a finite dimensional manifold when ω 0 satisfies the jet vortex ansatz.
Having identified a symplectic manifold, Orb(ω 0 ), we can then ask the question "are the dynamics Hamiltonian on Orb(ω 0 )?" Of course, the answer is "yes". This is the primary content of Marsden and Weinstein (1983). We provide our own explanation here for convenience.
For a general vorticity distribution ω, we may consider the kinetic energy Hamiltonian where ω may be of the form (8). In order to find Hamilton's equations on Orb(ω 0 ) choose some ω ∈ Orb(ω 0 ) and calculate the vector X H (ω) tangent to Orb(ω 0 ) given by Hamilton's equations. It must be the case that X H (ω) = £ u [ω] for some (nonunique) vector field u = (u, v) ∈ X div (R 2 ). Our goal is to solve for u. By the definition of the Hamiltonian vector field X H , we see that for any u = (u , v ) ∈ X div (R 2 ) If we let ψ := G δ * ω = G δ (· −z)ω(z)dz, then integration by parts implies We see that u = (−∂ y ψ, ∂ x ψ) is one possible solution. As is nondegenerate on the tangent spaces of Orb(ω), this is the unique solution. As a result, the evolution prescribed by X H is precisely (1). This proves that (1) can be seen as a Hamiltonian equation on Orb(ω 0 ) with respect to the symplectic structure (21) and the Hamiltonian (22).

The First-Order Case
Let us illustrate these Hamiltonian results for the case of the first-order MVB. Let z 1 , . . . , z n ∈ R 2 be distinct and define the initial vorticity distribution We desire to determine the coadjoint orbit, Orb(ω 0 ), and the symplectic structure. Indeed, we find that for any function f Collecting like terms we find By varying ϕ, we can obtain any collection of distinct points z 1 , . . . , z n ∈ R 2 and any collection of nonzero vectors 1 , . . . , n ∈ R 2 \{0}. This proves To derive the symplectic structure, recall the symplectic structure for a general vorticity (21). Now let ω = γ i δ z i + x i ∂ x δ z i + y i ∂ y δ z i . In this case, the left-hand side of (21) can be computed with respect to divergence-free vector field u = (u, v) and u = (u , v ) as Note that this is written entirely in terms of the first-order Taylor expansion of u and u evaluated at z i . Moreover, £ u [ω] = γ 0 u(z i )∂ x δ z i + . . . also has the property that it only depends on the first-order Taylor expansion of u and v at the points z 1 , . . . , z n . Therefore, both sides of (21) can be written as a function of the finite collection of numbers u(z i ), Du(z i ), v(z i ), Dv(z i ). The result then follows by identifying the scalars This proves that the symplectic structure on Orb(ω 0 ) is more concretely written as In essence, we have determined a finite dimensional Hamiltonian system whose solutions solve (1) when ψ is related to ω via an appropriate regularization.
Remark 9.1 The use of this symplectic structure shows that the map Remark 9.2 The corresponding Poisson bracket can be represented in tabular form by:

{·, ·}
x y x y The way to use this table is as follows. Let H = H (ξ ) be our Hamiltonian where ξ = (x 1 , . . . , x n , y 1 , . . . , y n , 1 , . . . , n , x 1 , . . . , x n , y 1 , . . . , y n ). Hamilton's equations are then given by for any function f , where B i j denotes the corresponding entry of the table. In particular, when f = ξ i , one recovers the equations of motion for the dynamics of the positions and strengths for a set of n first-order MVBs. Poisson geometers call B i j a Poisson tensor Abraham and Marsden (1978).

Conclusion
In this paper, we have considered a generalization of the standard vortex blob method, obtained by augmenting the vortices with higher order circulation variables and dubbing them multipole vortex blobs (MVBs). By viewing the vorticity equation as an advection equation, we have obtained equations of motion for these MVBs.
The extra degrees of freedom of MVBs resulted in richer dynamics near the vortex core. Moreover, these new vorticity carrying elements exhibited a variety of novel types of solution behavior. We also observed faster convergence rates in space using higher order MVBs. Moreover, we proposed a scheme to decrease the number of pairwise interactions, by grouping MVBs of lower order into a smaller number of MVBs of higher order. Lastly, the implications of Kelvin's circulation theorem were substantially richer in the case of MVBs than they were for the standard vortex blob method.
We have demonstrated the behavior of the MVBs with a sequence of simple numerical experiments consisting of small numbers of MVBs of various degrees. We found that first-order MVBs correspond to sums of vortex blobs and regularized dipoles which simply propagate themselves forward, while the second-order circulation variables activate richer (non-propagating) dynamics near the vortex core.
Finally, we derived the symplectic structure of MVBs using methods from Marsden and Weinstein (1983). The resulting structure turned out to be a direct generalization of the standard symplectic structure for vortex blobs.
The multiscale nature of ideal fluids is the principal obstacle to obtaining accurate models (Chorin 1994, Chapter 3). The use of MVBs augments the standard vortex blob method by allowing for singular vorticity distributions which model dynamics below the regularization length scale (i.e., at order δ k with δ 1 for a kth-order jet vortex). As the dynamics of MVBs are relatively easy to derive, and their analysis is tractable, we believe that MVBs will be of considerable value in understanding the place of regularized fluid models within the computational fluids community at large and they should provide renewed interest in the vortex blob method.
Future avenues of inquiry could include: • MVBs on manifolds, such as the sphere • The convergence properties of the MVB method • How does one choose the regularization length scale in relation to the grid resolution. This relationship is addressed quite well for zeroth-order MVBs in Beale and Majda (1982). It is not clear if higher order MVBs change those results. • An investigation of the kinetic theory of MVBs.
when ω is not smooth (e.g., a Dirac delta distribution), then one needs to invoke the mathematics of distributions as distinct from that of real-valued functions. Therefore, we have included this appendix to remind the reader of the basic theory of distributions. The main reference for this section is Hörmander (2003).
The space of distributions D (R 2 ) is the dual vector space to the space of smooth functions with compact supper C ∞ 0 (R 2 ). Therefore, a distribution is defined by how it maps functions to real numbers.
The distributional derivative of a given distribution ω ∈ D (R 2 ) in the ith coordinate direction may be defined as the distribution ∂ i ω obtained by For example, the Dirac delta distribution, δ 0 , is defined as the unique distribution such that The distributional derivative, ∂ i δ 0 , is given by Given a distribution ω ∈ D(R 2 ) and a function g ∈ C ∞ 0 (R 2 ), one can define the distribution g ω as (g ω)(z) f (z)dz = ω(z)g(z) f (z)dz, ∀ f ∈ C ∞ 0 (R 2 ).
For example, we find that g δ 0 = g(0)δ 0 . A slightly more involved, but standard, example is given by the computation of g ∂ i δ 0 . We find Therefore, On the left-hand side, note that g(0) and ∂ i g(0) are merely real numbers, which are multiplying the distributions ∂ i δ 0 and δ 0 . More generally, we find for a (perhaps nonunique) divergence-free vector field u = (u, v). Under this identification, the symplectic form at some ω ∈ S is given by an application of Kostant's formula Abraham and Marsden (1978). This is derived in Sect. 9 and found to be ω (£ u 1 [ω], £ u 2 [ω]) := ω(z)(u 1 (z) × u 2 (z))dz.
where u 1 and u 2 are divergence-free vector fields, and × denotes the planar crossproduct. Here we interpret the planar cross-product as taking values in the space of real numbers so that u 1 × u 2 is merely a smooth function. We now will translate formula (25) to this more specific scenario. Assume G acts upon R 2 , then G also acts upon distributions and upon S by symplectic group actions. In this context, a momentum map J associated with a G-action is defined by the equation where ξ · z denotes the action of ξ ∈ g on z ∈ R 2 , and "×" denotes the planar cross-product.

Translational Symmetry and J lin
The group R 2 acts upon R 2 by translation. That is to say, by sending anyz ∈ R 2 to z +z ∈ R 2 for anyz ∈ R 2 . This fact induces an action on smooth functions. In particular, there is a natural (right) action on C ∞ (R 2 ) sending the function φ(z) to the function (z) * φ(z) := φ(z −z). We denote the inverse operation by (z) * φ(z) := φ(z +z). This induces a (left) action on distributions which sends ω ∈ D (R 2 ) to the distribution (z) * ω(z) := ω(z +z). As a translation of R 2 byz is a volume-preserving diffeomorphism, we see that the coadjoint orbit S ⊂ D(R 2 ) is invariant under this action. Moreover, we observe the action, restricted to S, is symplectic because Given this symplectic action, we can seek a momentum map, J lin : S → (R 2 ) * . Consider an arbitrary element of the Lie-algebra δz = (δx, δỹ) ∈ R 2 and use (26)  Integrating by ω we find If ω satisfies the MVB ansatz The terms of the MVBs beyond the first order do not influence J lin .