## Abstract

The notion of a *finite dipole* is introduced as a pair of equal and opposite strength point vortices (i.e. a vortex dipole) separated by a finite distance. Equations of motion for *N* finite dipoles interacting in an unbounded inviscid fluid are derived from the modified interaction of 2*N* independent vortices subject to the constraint that the inter-vortex spacing of each constrained dipole, ℓ, remains constant. In the absence of all other dipoles and background flow, a single dipole moves in a straight line along the perpendicular bisector of the line segment joining the two point vortices comprising the dipole, with a self-induced velocity inversely proportional to ℓ. When more than one dipole is present, the velocity of the dipole centre is the sum of the self-induced velocity and the average of the induced velocities on each vortex comprising the pair due to all the other dipoles. Each dipole orients in the direction of shear gradient based on the difference in velocities on each of the two vortices in the pair. Several numerical experiments are shown to illustrate the interactions between two and three dipoles in abreast and tandem configurations. We also show that equilibria (multi-poles) can form as a result of the interactions, and we study the stability of polygonal equilibria, showing that the *N*=3 case is linearly stable, whereas the *N*>3 case is linearly unstable.

## 1. Introduction

This paper formulates the dynamics of multiple finite dipoles. The primary motivation is to develop models for fish schooling that account for the hydrodynamic interactions among school members. Fish schools are often examined using behaviour-based models that consider homogeneous particles (fish) interacting locally based on rules of repulsion, alignment and attraction to other fish (Couzin *et al.* 2002, 2005). While capable of exhibiting visual patterns similar to those seen in biological schools from disordered swarming to ordered parallel and circular motions (see Parrish *et al.* 2002), these behaviour-based models do not account for the fluid medium. Yet, it is widely believed that the hydrodynamic coupling plays an important role in fish schooling. In his seminal paper, Weihs (1973) uses von Kármán wake models to examine optimal positioning of an individual fish within the school, thus deducing that a diamond shape pattern was beneficial for schooling. These are *steady-state models*, where the von Kármán wakes, and thus the whole school, translate with a constant velocity. More recently, the effects of near-field vorticity, presumably produced by a neighbouring fish, on the dynamics of a single fish have been examined both experimentally (Liao *et al.* 2003, Muller 2003) and analytically/computationally (Kanso & Oskouei 2008; Eldredge & Pisani 2008; Alben 2010; Oskouei 2011). Other recent studies focused on the far-field effects of the hydrodynamic interactions, captured via the potential flow component, on the motion of two submerged bodies (Nair & Kanso 2007; Tchieu *et al.* 2010). These models are yet to be scaled to produce dynamical models for fish schooling.

In this paper, we propose the finite-dipole dynamical system as a model that captures the far-field fluid effects of multiply-interacting swimming bodies. To this end, consider the dynamics of a collection of several bodies propelling themselves in a two-dimensional inviscid fluid as shown, for example, in figure 1*b*. The smallest length scale, ℓ, is associated with the finite size of each body, which, to leading order, produces a dipolar velocity field as depicted in figure 1*a*. The largest length scale, *R*, is associated with the separation distance between bodies. For well-separated cases in which ℓ≪*R*, one might imagine that a dynamical system consisting of a collection of interacting finite-sized self-propelled dipoles would be a reasonable model governing the hydrodynamic interactions among the bodies. The proposed system is depicted in figure 1*c*, where the bodies have been effectively replaced by their finite, self-propelled dipole counterparts.

We derive a dynamical system based on a constrained set of 2*N* interacting point vortices, paired with equal and opposite strengths with each pair constrained to remain a distance ℓ apart (hence the term finite-sized), and we investigate some of the dynamical properties associated with the system. Note that interacting dipolar fields have been considered by other authors. Newton (2005) derives a *point dipole* dynamical system by considering the limit for each dipole (made up of equal and opposite strength point vortices), making the assumption that the self-induced velocity vector of each dipole aligns itself in the direction of the background velocity field at the dipole centre. The infinite self-induced velocity associated with the limit is dealt within two ways. First, it is subtracted out of the field, allowing each dipole to move according the field generated by all the others. Second, a regularized version is proposed and compared with the first approach. In the regularized version, a regularizing parameter is introduced into the denominator of each dipole, much like regularizations in point vortex systems such as Krasny (1986). Yanovsky *et al.* (2009) and Kulik *et al.* (2010) derive equations of motion by representing the vorticity field in terms of generalized functions and their derivatives. The most recent paper of Llewellyn Smith (2011) reviews previous methods, including those associated with point vortex systems and formulates a new method to derive equations of motion for general point singularities using generalized momentum arguments.

In addition to these studies on point dipole dynamics, earlier work such as that of Love (1893), Eckhardt & Aref (1988), Price (1993) and subsequently Tophoj & Aref (2008) analyse the complex and chaotic interaction of four point vortices paired up, initially, in two tight groups of equal and opposite pairs, hence two *unconstrained dipoles*. More specifically, these papers take a look at the scattering properties of the unconstrained interactions between the two pairs of vortices. The work of Newton & Shokraneh (2008) formulates the general problem of 2*N* point vortices grouped into equal and opposite pairs on a rotating sphere, hence an unconstrained *N*-dipole problem on the sphere, while that of Kimura (1999) formulates the problem of isolated dipolar motion on surfaces of constant curvature (e.g. spheres) and proves that they follow geodesic paths (e.g. great circles).

The organization of this paper is as follows. Section 2 addresses the formulation of the finite-dipole problem in the unbounded plane as a constrained 2*N* point vortex system. This is done by directly modifying the standard point vortex equations of motion to respect the constraint that each pair of point vortices of equal and opposite strength remains a fixed distance ℓ apart. Examples of two and three finite-dipole systems are addressed in §3. In §4, we investigate the formation of multi-pole equilibria and subsequently discuss the stability of such equilibria. We conclude with §5 and discuss the outlook of employing the finite-dipole system as a model of biological organisms swimming in potential flow as well as giving an example touching on aspects of formation motion.

## 2. Formulation

### (a) Equations of motion

We formulate the problem in an inviscid, unbounded, two-dimensional fluid. Consider *N* pairs of point vortices of equal and opposite strengths, placed a distance ℓ_{n} apart in the complex *z*-plane (where *z*=*x*+i*y* and ) for *n*=1,…,*N* as depicted in figure 2. A single pair is shown in the detailed schematic. The positive (+*Γ*_{n}) and negative (−*Γ*_{n}) vortices are located at
2.1a
respectively, where *z*_{n} denotes the midpoint lying between the pair, i.e. the vortex dipole centre
2.2
and *α*_{n} represents the orientation of the vortex dipole with respect to the *x*-axis. Corresponding to figure 2, we will refer to the positive vortex as the left vortex and the negative vortex as the right vortex as identified by the respective subscripts.

It is well known that the kinematic solution to the incompressibility condition in inviscid flow amounts to determining a complex potential, *F*(*z*)=*ϕ*+i*ψ*, that is analytical in the domain (except at singular points) and decaying at infinity. Here, *ϕ* and *ψ* are the so-called (real-valued) velocity potential and streamfunction, respectively. For the system provided in figure 2, *F*(*z*) is constructed from a linear superposition of point vortex potentials at their respective positions,
2.3
The complex velocity field, *w*=*u*−*iv*, is given by the derivative of (2.3) with respect to *z*,
2.4
It can also be confirmed that in the limit and , equation (2.3) reduces to a sum of point dipoles at *z*_{n} with scalar strength *μ*_{n} and orientation *α*_{n}, that is to say, one has
which is in agreement with Newton (2005, see eqn (6)).

We seek to determine equations of motion governing the interaction between all vortex dipoles as functions of time, while the vortex pair in each dipole is constrained so that . These constrained vortex pairs are aptly named *finite dipoles*, as they maintain a constant, finite distance between the vortex pair for all time. This amounts to deriving equations of motion for all vortex dipole centres, *z*_{n}, and the orientations, *α*_{n}. To this end, we consider the 2*N*-vortex problem (see Newton 2001) subject to an additional inter-dipole velocity to allow us to apply the constraint, . Each respective vortex of each individual finite dipole convects with the modified velocity
2.5a
and
2.5b
where and represent the complex conjugate and time derivative, respectively, and *λ*_{n} is real. The first term in equation (2.5) represents the self-induced velocity and is given by
2.6
that is, an isolated vortex dipole moves in a straight line in the direction defined by *α*_{n} at the speed *Γ*_{n}/(2*π*ℓ_{n}). The second term, *w*_{n,o}, represents the velocity induced by all other finite dipoles and has a similar form to (2.4), namely
2.7
The third term, ±i*λ*_{n}e^{−iαn}, is an additional, attractive (*λ*_{n}>0) or repulsive (*λ*_{n}<0) inter-dipole velocity that allows for a constraint on ℓ_{n} to be applied. This velocity is modelled so that *z*_{n,l} and *z*_{n,r} move towards or away from each other at the speed *λ*_{n} along the line joining these two vortices. Typically, in the unconstrained interaction of *N*-vortex dipoles in an inviscid fluid (or equivalently, 2*N* point vortices), the inter-dipole spacing ℓ_{n} is not constant, see, for example, Eckhardt & Aref (1988), and the equations of motion for each vortex are exactly given by equations (2.5) with *λ*_{n}=0. The introduction of non-zero *λ*_{n} introduces another degree of freedom to enforce the constraint that , much like applying Lagrange multipliers in constrained mechanics. This allows each finite dipole to retain its ‘particle-like’ identity throughout its time evolution.

We now rewrite (2.5) as a system of equations governing the motion of the centre *z*_{n} and the orientation *α*_{n} of each finite dipole. To this end, we take the conjugate and time derivative of equation (2.2) and substitute equations (2.5) to get the equation of motion for the dipole centre,
2.8
One can see that the velocity of the dipole centre is a sum of the self-induced velocity and the average of the velocity that each left and right vortex feels from all other dipole pairs. In instances ℓ_{n}/*R*≪1, where *R* is some character separation distance between dipoles, the motion of the dipole centre *z*_{n} can be approximated with the background velocity field evaluated at *z*_{n} (as in Newton 2005; Yanovsky *et al.* 2009; Llewellyn Smith 2011) but with an additional self-induced velocity term specified by *w*_{n,s}.

In similar fashion, taking the time derivative of the separation between the vortices gives
2.9
Imposing the constraint distinctly determines the newly added parameter
2.10
and enforcing (2.10) allows us to write (2.9) as
2.11
Equations (2.8) and (2.11) form a closed system of 3*N* real equations (*N* complex+*N* real) governing the motion of *N* finite dipoles interacting in an unbounded plane. We note, for this reason alone, the system cannot be put directly into canonical Hamiltonian form of 2*N* real equations. The equations are physically based on the interaction of point vortices and allow for a finite self-induced velocity of each dipole. It is also worthy to note that in this formulation, the dipole centre *z*_{n} is not singular. Classically, a singular point dipole induces infinite velocity through its centre. Without proper desingularization, neighbouring point dipoles have the tendency to advect through the singular points of other neighbouring dipoles, as can be simulated in the formulations given by Newton (2005), Yanovsky *et al.* (2009) and Llewellyn Smith (2011). Our formulation circumvents this possibility and the equations of motion remain well behaved in this circumstance.

### (b) Dynamic translation and alignment of the finite dipole

Equation (2.8) implies that the translational motion of the centre of a finite dipole depends on its self-induced velocity as well as the average of the velocities of the two vortices owing to other dipoles and background flows. That is to say, the centre of the dipole need not always translate in the direction of the self-induced velocity (i.e. along a straight line) because, in general, the background flow can induce a velocity on the centre in the direction and/or lateral to the self-propelled velocity.

Equation (2.11) implies that the instantaneous evolution of the orientation of a finite dipole depends on the difference in velocities between the two vortices that comprise the dipole. More specifically, one can see from equation (2.11) that the important quantity defining the instantaneous evolution of *α*_{n} is which is the projection of the difference in velocities onto the direction of self-induced motion. The difference between the velocities *w*_{n,o}(*z*_{n,r}) and *w*_{n,o}(*z*_{n,l}) defines the local shear of the fluid. We discuss these effects on the motion and alignment of the dipole through examples.

Consider the case of a single isolated dipole or a dipole moving in a locally uniform flow *U* such that *w*_{n,o}(*z*_{n,r})=*w*_{n,o}(*z*_{n,l})=*U* (where *U* is a constant, possibly zero), as shown in figure 3*a*. Since both vortices comprising the dipole feel the same background velocity, equation (2.11) implies that is precisely zero with no dependence on the direction of the background flow. Thus, even when the dipole is not aligned with the background flow direction, the dipole simply translates along a straight line combining its self-induced velocity and the velocity of the locally uniform flow. Because of the lateral components of the background velocity, the trajectory of the dipole centre will be altered to move in a direction different from its self-induced velocity, as shown, but the dipole does not alter its orientation. The model of Newton (2005), by contrast, forces the self-induced velocity vector to align itself with the uniform background flow.

Consider now the example of a single dipole in a shear flow as shown in figure 3*b*. Because one of the vortices of the dipole feels a larger component of background velocity than the other, the dipole aligns its orientation in the direction of the shear gradient. In this example, the top vortex feels a larger induced velocity than the bottom vortex and thus rotates in the clockwise direction while moving on the trajectory shown, owing to the lateral components of the background flow field. Note that one has only if the dipole moves perpendicular to the direction of local shear. It is clear through these two examples that dipoles turn in the direction of the shear gradient as they move on curved paths which take into account not only the self-induced velocity vector, but also the parallel and lateral projections of the background field on this vector.

It is worth emphasizing that this is the immediate behaviour and not the long-time behaviour of the coupled system. The dynamics of the coupled system governed by equations (2.8) and (2.11) is much more rich, complex and likely chaotic for *N*>1. Consider, for example, the dynamics of a finite dipole located at *z*_{1} with orientation *α*_{1} in the presence of the flow field generated by another finite dipole at *z*_{0}=0 with orientation *α*_{0}=*π*/2. We fix *z*_{0} and *α*_{0} so that the dipole at *z*_{0} exhibits no dynamics and acts as a background flow field on the *z*_{1}-dipole. Two sample trajectories are given in figure 4. Initially in figure 4*a*, and the dipole turns counterclockwise towards the velocity shear gradient because the left vortex (i.e. the vortex closer to the origin) feels a larger velocity in the negative *y*-direction than the right vortex. After the initial turn, the dipole makes a small adjustment before heading along a straight line at its self-induced velocity as the effect of the background flow gradually diminishes.

The second example given in figure 4*b* shows clearly how the translation and orientation equations couple to produce more complex dynamics. Initially, the dipole is oriented to travel in the direction of the gradient, thus . It appears that the dipole should always be oriented in the *α*_{1}=0 direction since it moves across streamlines, but the local velocity field moves *z*_{1} in the lateral direction and thus the dipole begins to move downward. Subsequently, the dipole begins to reorient itself since its equation of motion is coupled to the translational dynamics of the dipole. Eventually, the dipole moves far enough that the influence of the background flow becomes negligible and it continues to travel along a straight line with its self-induced velocity. It is in this sense that the long-time dynamics for a collection of finite dipoles is difficult to predict based solely on the instantaneous background flow.

### (c) Force acting on a finite dipole

A free-point vortex in an inviscid fluid advects force free so that momentum is conserved everywhere except at the singular point vortex location. This result is due to Kirchhoff (1876) and is discussed in Lamb (1945). In the constrained dipole model presented here, it is expected that each vortex of the finite dipole experiences an additional force to maintain its constrained motion. This is also reflected in the additional velocity component ±i*λ*_{n} e^{−iαn} in equations (2.5). In this section, we derive an expression for the constraint force and interpret its effect on the constrained dipole.

Consider the fluid momentum in the moving contour surrounding a single vortex of the dipole pair as pictured in figure 5. The force applied to a finite dipole satisfying the fixed-separation constraint is equivalent to the time rate of change of momentum for the fluid enclosed in contour . Following Michelin & Llewellyn Smith (2009), the complex force is given as
2.12
where *w*_{c} is the velocity of the contour. The *x*- and *y*-components of the force are given by and , respectively.

As shrinks infinitesimally small to surround only the vortex pair and the branch cut connecting the two vortices, the contribution from becomes identically zero and therefore . Now consider the integrals presented in (2.12) around . The complex potential and velocity can be decomposed as
2.13
and
2.14
where and are the desingularized potential and velocity field, respectively. That is to say, in the vicinity of *z*_{n,l}, the functions and remain analytical and single-valued. Noting that as shrinks to surround solely the point vortex, . Now the integrals in (2.12) can be exactly evaluated around ,
2.15a
2.15b
and
2.15c
Similarly, the integrals around can be evaluated and substituted, along with (2.5), into (2.12) to obtain
2.16
This expression can be simplified further by substituting (2.9) into (2.16) and noting that and . To this end, one gets
2.17
Each individual vortex contributes exactly to equation (2.17) thus there is no additional moment about *z*_{n}. As expected, when no constraint is applied (*λ*_{n}=0), the dipole pair convects force free. In the presence of other point vortices, a vortex belonging to the finite dipole, except under very special conditions, will attempt to separate from or be attracted to its partner. When *λ*_{n}>0, a force is applied in the same direction as the self-induced velocity. This is contrary to applying a force along the line joining the two, which may seem more intuitive since the constraint acts along this line. The reason that the force is applied in the direction of travel is that, if a vortex dipole begins to separate, the pair itself begins to lose its forward momentum. In this instance, to enforce the fixed-separation constraint, the dipoles are pushed and given extra momentum to maintain a constant speed and separation. Therefore, an additional force in the direction e^{iαn} is necessary whenever the constrained variable *λ*_{n}>0, and vice versa for *λ*_{n}<0.

## 3. Finite-dipole interactions

We present numerical examples of the dynamic interactions of finite dipoles in free-space with the emphasis placed on investigating the motion of dipoles initially travelling in the same direction. Each dipole is given a self-propelled speed |*w*_{s,n}|=1 by specifying *Γ*_{n}=1 and ℓ_{n}=ℓ=1/(2*π*). In what follows, the initial orientation is specified as *α*_{n}=*π*/2 (thus initially travelling upward). Keep in mind that these finite dipoles have no control or collision-avoidance mechanisms and act solely under the evolution of the equations given in equations (2.8) and (2.11). These equations are solved numerically in Matlab using the packaged solver ode45 (a variable-step, explicit time integrator) with relative and absolute tolerances of 10^{−6} and 10^{−8}, respectively.

We begin by considering the interaction of two dipoles. We initially place the dipoles a horizontal distance *b* and vertical distance *h* apart, as defined by figure 6, and allow the dipoles to interact freely. Several regimes of motion occur in this example and are depicted in figure 7 for a total integration time *t*=5.

When initializing the dipoles perfectly in tandem to one another (*b*=0), the dipoles remain in relative equilibrium with travelling speed
3.1
that is to say, the dipoles help each other in travelling forward at a rate faster than their self-induced velocity, a phenomena we call ‘dipole drafting.’ As , the dipoles approach a maximum speed of twice its self-propelled speed. A sample simulation where *h*=4ℓ is given in figure 7*a*. As we increase *b* while holding *h*=4ℓ, the dipoles initially assist each other in forward motion, but then diverge quickly (figure 7*b*). As *b* is further increased, the rate of divergence decreases (figure 7*c*). It is interesting to note that this diverging behaviour of the trailing dipole is also seen in numerical simulations of a drafting cylinder behind a forced cylinder (Tchieu *et al.* 2010). In contrast to Tchieu *et al.* (2010), here the two dipoles are dynamically coupled with one another.

We decrease *h* while holding *b*=4ℓ. For a critical value of *h*=2.25ℓ, the pair moves in a relative equilibrium but in an oblique direction to their self-induced velocity as depicted in figure 7*d*. There is no propulsive benefit for the dipoles because the travelling speed of this equilibrium is less than the self-propelled speed. Still, it is intriguing that two dipoles can influence each other to change their overall direction of motion without changing their orientations. Interestingly enough, as we decrease *h*<2.25, the two dipoles begin to oscillate in what seems close to a periodic trajectory (modulo translations). This is illustrated in figure 7*e*. The dipoles maintain this oscillation throughout the integration time. As we decrease *b*<0.83ℓ the dipoles enter a regime where they collide with one another. For example, with the dipoles initially starting abreast (*h*=0), the two dipoles eventually collide symmetrically and annihilate each other as shown in figure 7*f*. Indeed, the dipole dynamical system can be viewed as a low-order model for interacting rigid bodies in which case the collision of the two entities is not unexpected. Notice that although the dipoles collide in finite time, the evolution equations are well defined because the equations describe the evolution of the centre *z*_{n} which is a non-singular point. As the vortices of the respective pairs approach one another in the collision, they cancel the effect of the nearest neighbour and the dipoles are annihilated.

We now consider the case of three finite dipoles. The dipoles are initialized as schematically drawn in figure 6 with the horizontal and vertical spacings, *b* and *h*, respectively. In figure 8*a*, the dipoles are placed in abreast formation (*h*=0). The two outer dipoles induce a velocity in the negative *y*-direction on the centre dipole and travel faster than the centre dipole. The outer dipoles approach the centreline, turn inward owing to the influence of the centre dipole, and then system approaches a fixed equilibrium where the lines joining the vortex pairs of each respective dipole lie on the edges of an equilateral triangle. Although it may seem fortuitous, this equilibrium is attained for situations where *h*<0.36ℓ and seems to be an attracting state for this regime of initial conditions.

As we increase *h* to 0.36ℓ<*h*<1.5ℓ, the trailing pair of dipoles collide as seen in figure 8*b*. Interestingly enough, when 1.50ℓ<*h*<5.83ℓ, the dipoles enter a regime where the outer initially angle themselves towards the centre and then subsequently diverge (figure 8). When integrating the dynamics for longer time as done in figure 9, we see that the system continues to oscillate with growing amplitude for long times, which is characteristic for this regime. For *h*>5.83ℓ, owing to the separation between the lead and trailing vortices, the influence of the lead dipole is too small to counteract the tendency of the trailing pair of dipoles to collide as in figure 7*f*.

## 4. Multi-pole formation and stability

The occurrence of the three-dipole equilibrium is rather peculiar and suggests the question of whether other polygonal equilibria exist and whether or not they are stable. In similar fashion, more than a century ago, Lord Kelvin (Thomson 1878) discussed the relative equilibria of point vortices lying at the vertices of a polygon, and Thomson (1883) revisited the problem and investigated the stability of this configuration. One can consult Aref *et al.* (2003) for a review on point vortex equilibria. Here, we study the analogous case applied to the absolute equilibria, , of the finite-dipole dynamical system.

There is an important distinction between the *free* point vortex dynamical system and the finite dipole dynamical system discussed here. The free point vortex system, as a discretization of the Euler equations of ideal fluid flow, is Hamiltonian, as is well known. A characteristic property of Hamiltonian systems is that equilibrium configurations cannot be reached from generic initial conditions in finite time. In other words, vortices have to be initially placed in their equilibrium configuration and cannot dynamically attain an equilibrium. This is not true for the finite-dipole system since, as mentioned earlier, the force required to maintain the constraint breaks the Hamiltonian structure of the problem. This is seen very clearly in a Hamiltonian energy plot shown in figure 10*a* of the equivalent unconstrained point vortex system. In addition, the equilibrating constraint, which is directly proportional to the magnitude of the constraining force (2.17), is given in figure 10*b*. It is clear that the point vortex Hamiltonian is not conserved in general.

Indeed, in figure 8*a*, three dipoles starting abreast reached, in finite time, what appears to be an equilibrium configuration where all the vortices constituting the dipoles are placed equidistantly on a circle of radius *r*. We now examine the general case of fixed equilibria of *N* dipoles more closely.

### (a) Equilibria

We wish to find equilibrium configurations for *N* finite dipoles of strength *Γ*_{n}=*Γ* and characteristic separation ℓ_{n}=ℓ. We specifically look for equilibria where the dipoles are placed equidistantly on the circle of radius *r*, that is,
4.1
The orientation is initially prescribed to move towards the centre of the formation,
4.2
In this configuration, . Substituting equation (4.1) and (4.2) into (2.8) gives
4.3
where *A*_{N}(*r*;ℓ) and *B*_{N}(*r*;ℓ) are polynomial with the highest power of *r* being *r*^{N+1} for *N* odd and *r*^{N} for *N* even. The positive real roots of *A*_{N}(*r*;ℓ) are the radii that define the equilibria. For *N*=2, *A*_{2}(*r*;ℓ)=−2*r*^{2} and the only solution is the trivial case where *r*=0. For *N*=3, *A*_{3}(*r*;ℓ)=48*r*^{4}−24ℓ^{2}*r*^{2}−ℓ^{4} and the only real solution is . For *N*=4, *A*_{4}(*r*;ℓ)=−8*r*^{4}+10ℓ^{2}*r*^{2} bifurcates into two solutions; the trivial solution, *r*=0, and . For *N*≥5, owing to the increasing order of *A*_{N}(*r*;ℓ), equilibrium radii are found numerically. In figure 11, the non-trivial radii *r* are plotted versus *N*, the number of dipoles in the system. There are two non-trivial solutions for *r* for *N*=5, 6, three for *N*=7, 8, four for *N*=9, 10 (the fourth solution not shown in figure 11 for clarity) and so on.

Figure 12 depicts the configuration and associated streamlines of the first set of solutions for various *N*. These are the equilibria for which the radii are the largest (circle markers in figure 11). Characteristically, the lines joining the vortex pairs do not intersect in this configuration. In fact, each vortex is equidistant from its nearest neighbour. The growth of *r* for this set of equilibria is nearly linear with a slight departure at small *N*.

The additional equilibria that exist are defined with increasingly smaller radii (triangle and square markers in figure 11). These alternative equilibria can be viewed as ‘higher harmonics’ of the primary equilibria presented in figure 12. This is more clearly shown in figure 13, which depicts the hierarchy of the *N*=10 equilibria. The lines joining the dipoles in figure 13*a* now intersect with two other neighbouring dipoles in this configuration. Although the locations of the dipole centres are defined, the equilibrium solution requires that the vortices comprising the dipoles must also remain equidistant. In the subsequent equilibria presented in figure 13*c*, the line joining the dipoles intersects with four and six others, respectively, and all individual vortices remain equidistant.

We remark that in contrast to this finite-dipole problem, the free-point vortex dynamics of this 2*N* configuration of alternating strength vortices is not in equilibrium. When unconstrained, the vortices of adjoining finite dipoles in these formations eventually pair and travel towards (e.g. in the first set of equilibria solutions presented in figure 12, the right vortex of *n*=1 partners with the left vortex of *n*=2). The constraint modelled here allows for a finite force to be applied to each finite dipole to maintain its equilibrium in the system.

### (b) Stability

We now examine the linear stability of the equilibria presented in figure 12 numerically. The nonlinear system of equations (2.8) and (2.11) is rewritten as a system of real equations in standard state-space form as
4.4
where **x**=(Re[*z*]_{1},…,Re[*z*]_{N},Im[*z*]_{1},…,Im[*z*]_{N},*α*_{1},…,*α*_{N})^{T} and **F**(**x**) is a 3*N*-dimensional vector function defined by the right-hand sides of equations (2.8) and (2.11). Equations (4.4) are linearized about the equilibrium positions. The linearized equation of motion take the form:
4.5
where *δ***x** denotes an infinitesimal perturbation about the equilibrium state and **J**=**∇****F** is the 3*N*×3*N* Jacobian matrix evaluated at the equilibrium of interest. Owing to the number of states when *N*≥3, the **J** is computed numerically using finite differences. The maximum non-zero eigenvalues for several values of *N* are presented in table 1. It is observed for *N*>3 and odd, the maximum eigenvalues have double multiplicity.

We examine the eigenvalues and eigenvectors of **J** more closely. For *N*=3, the matrix **J** (9×9) is degenerated and has five zero eigenvalues, suggesting that there are five invariants of motion. The first three invariants are due to the invariance in *x*-translations, *y*-translations, and in-plane rotation of the initial configuration. The additional two invariants stem from the fact that reflections about the *x*- and *y*-axes do not change the system. The four other eigenvalues are negative and thus stable. The stable modes of the *N*=3 system are depicted in figure 14. There are four modes that are denoted as expansion/contraction, shear, opening/closing and tilting modes (figure 14*a*–*d*, respectively). The linear superposition of these four stable modes constitute any deformation for which the system is subsequently stable. The stability owing to small perturbations is numerically verified.

For *N*≥3, the polygonal configuration of dipoles also exhibits the same invariants. For *N*=4, **J** has one positive unstable eigenvalue that causes the system to collapse in a straining mode (which is defined by its associated eigenvector). The evolution of this unstable mode is given in figure 15. When perturbed in this fashion, the two dipoles separate, strain and eventually collide with each other. Numerically, it is confirmed that the other stable modes do not excite an immediate collapse, although the system eventually becomes unstable owing to numerical noise. For *N*≥4, it is verified that the equilibria are linearly unstable.

We note that this linear analysis does not guarantee stability owing to finite perturbations. A nonlinear stability analysis is necessary for this and beyond the scope of this paper. It is also interesting to note the classical result of Thomson (1883), which was later corrected by Havelock (1931) (also see, Dritschel 1985; Aref 1995, for modified forms) that the regular polygon configuration of identical point vortices is linearly stable when there are six or fewer vortices, neutrally stable when there are seven vortices and linearly unstable when there are eight or more vortices. In this analogous analysis for identical finite dipoles restricted to the polygon, the *N*=3 case, consisting of six vortices, is linear also stable, and the *N*≥4 cases, consisting of eight or more vortices, are also linearly unstable.

## 5. Conclusion and discussion

We derived equations of motion for the finite-dipole dynamical system by starting from the standard point vortex equations and imposing that the separation distance between the two vortices in each dipole remains constant. The enforcement of the constraint dictated that an additional force must be applied to each dipole, effectively increasing the momentum of the dipole when the unconstrained vortex pair begins to separate and decreasing the momentum of the dipole when the unconstrained pair begins to merge. The balance is maintained by the modified velocity term in equation (2.5) and results in a dynamical system that is not Hamiltonian in the canonical sense. The resulting equations of motion imply that the dipole centre convects with a self-induced velocity plus the average of the background velocity at the two vortices that comprise the finite dipole. The alignment of the dipole changes according to the projection of the difference in the velocities of the two vortices onto the direction of self-induced motion.

We examined the interaction of two and three finite dipoles through examples. We showed that two finite dipoles can exhibit steady forward and oblique motion, that is to say, move in relative equilibrium to each other. We also identified regimes where the two dipoles diverged or collided and others where they seemed to move in periodic or quasi-periodic motion. Interestingly, in a different context, the motion of a collection of vortex dipoles in the dilute-gas regime of a Bose–Einstein condensate was examined and quasi-periodic motion in the form of stable epicyclic orbits were observed (Middelkamp *et al.* 2011). It is, of course, well known that the equations governing quantized vortices in helium II and Bose–Einstein condensates are the same as that in ideal, incompressible fluids (see Donnelly 1991), so we certainly expect our model to be relevant in this context as well. Based on what we know from the previously mentioned works on the dynamics of unconstrained dipoles, we can anticipate that the long-time dynamics of a collection of finite dipoles will be chaotic, but that is outside the scope of the present paper. In the case of three finite dipoles, the emergence of dipole equilibria, or multi-poles, were observed and a class of stationary polygon equilibria were found. A linear stability analysis was then performed and it was found that the *N*=3 case (triangle) was linearly stable, and for *N*>3, the polygonal equilibria were linearly unstable to infinitesimal perturbations.

As mentioned in §1, our main motivation for considering this class of finite-dipole models stems from our interest in developing low-order models for fish schooling. To this end, we briefly examine the motion of dipoles initially placed in square and diamond formations, respectively. We consider these configurations because in the context of fish schooling, it has been argued that the diamond formation is favourable for schooling (see Weihs 1973). Weihs (1973) based his analysis on a stationary, infinite lattice and computed the locomotory benefits, a given fish in diamond formation gets from the vortical wakes (modelled as idealized vortex streets) of its neighbouring fish. Here, we view the dipole system as a model for a group of dynamically interacting, self-propelled swimmers in the absence of vortex shedding and only consider a finite number of dipoles (fish). In figure 16*a*, nine dipoles are placed initially in a square formation with spacing equal to . The dipoles are placed abreast and in tandem to one another. Figure 16*a* shows that the dipoles interact such that they break the square formation and form three of the equilibrium formations. In this case, the interactions between the dipoles is not beneficial for locomotion.

In figure 16*b*, the dipoles are placed in a staggered or a diamond formation and initially spaced with *b*=4ℓ and *h*=4ℓ (that is to say, same separation distance between dipoles as in the square formation). These finite dipoles have no control or collision-avoidance mechanisms and act solely under the evolution of the equations given in equations (2.8) and (2.11). Figure 16*b* shows that, for the same integration time as in 16*a*, though slightly diverging, the dipoles essentially retain their formation. The dipoles maintain much better spacing than in the case of the square formation, leading to more favourable conditions when travelling in a group than when travelling alone. Eventually, after some time that is much longer than the characteristic time it takes for the dipoles to reach a stationary equilibrium in the square formation, the dipoles either collide or diverge to the point where they act independently of one another.

These observations raise questions associated with the collective motion of self-propelled swimmers in the inviscid limit. Are there locomotory advantages to schooling? What is the role of the hydrodynamic coupling in motion coordination? How does the addition of control or decision-making features (e.g. attitude control on *α*_{n}) affect the collective behaviour of a group of dipoles? These questions will be addressed in a future study.

## Acknowledgements

P.K.N. was supported by The National Science Foundation under grant NSF-DMS-0804629. E.K. and A.T. are partially supported by the National Science Foundation under the CAREER award CMMI 06-44925 and the grants CCF 08-11480 and CMMI 07-5709.

- Received February 26, 2012.
- Accepted April 12, 2012.

- This journal is © 2012 The Royal Society