Finite-core quasi-geostrophic circular vortex arrays with a central vortex

We numerically determine equilibrium states for three-dimensional quasi-geostrophic vortex arrays. The vortex arrays consist of five or eight equal volume and equal uniform potential vorticity peripheral eddies whose center lies equally spaced along a circle and of a similar vortex at the center of the array. We are specifically interested in the influence of the height-to-width aspect ratio of the vortices on the arrays’ properties. The linear stability of the vortex arrays is addressed and the vortex arrays are shown to be sensitive to instabilities when the vortices are close enough together. The onset of instability corresponds to a threshold for the distance between the peripheral vortices and the center of the array. Measured as a fraction of the mean vortex horizontal radius, the stability threshold increases as the height-to-width aspect ratio of the vortices increases. For a separation larger than the stability threshold between the vortices, the arrays are linearly stable, hence robust and long-lived in the nonlinear regime. We also show that prolate peripheral vortices exhibit a concave outer side when they are close to the center of the array, while it is convex otherwise. The nonlinear evolution of a selection of unstable vortex arrays is examined. In such cases, the vortices deform and some vortices merge, breaking the symmetry of the vortex array. The later evolution of the unstable vortex arrays can be convoluted.


I. INTRODUCTION
Vortices, or swirling masses of fluid, are ubiquitous features of fluid flows.Vortex arrays may exist in a state of equilibrium, where the vortices neither deform nor move relative to each other, leading to the notion of vortex statics introduced by Thomson in 1875. 1 For such vortex arrays, the flow induced by the vortices is steady in a frame steadily moving with the vortices.Vortex arrays with enough symmetries can often be in equilibrium.Inspired by the work by Mayer 2 on floating magnets, Thomson 3 first considered circular arrays of two-dimensional point vortices.The linear stability of such circular two-dimensional vortex arrays was first addressed by Thomson. 4 Thomson concluded that circular arrays of six or less identical vortices were linearly stable.A correction to Thomson's results later showed that a circular array of seven identical vortices was also stable. 5,6Point vortex arrays are only sensitive to displacement modes of instability, which results in the motion of the vortices away from their equilibrium location.The study of circular two-dimensional vortex arrays was later extended to finite-core vortices. 7,8Finite-core vortices can also be sensitive to deformation modes of instability, which primarily affects the shape of the vortices rather than their location.][11][12] Large-scale atmospheric and oceanic vortices are strongly influenced by the planetary rotation and stable background density stratification.When these two effects are dominant, the dynamics of the vortices is accurately captured by an asymptotic model, the quasi-geostrophic (QG) model.The study of circular vortex arrays was first extended to atmospheric and oceanic flows using the quasigeostrophic shallow water model [13][14][15][16] and then in a two-layer QG system. 16,17][20] In particular, Reinaud 21 investigated the stability of QG point vortex arrays and finite-core vortex arrays for vortices having a unit height-to-width ratio.It was shown that only arrays of less than sixpoint vortices are linearly stable in this regime.Arrays with more vortices can be stabilized if a like-signed vortex is placed in the center of the array.The present contribution is an extension to this study for arrays of vortices with various height-to-width aspect ratios.In a continuously stratified fluid, varying the vortex heightto-width aspect ratio is equivalent to varying the Burger number of the vortices.
3][24] An array of five vortices rotates around a central vortex (a "5 + 1"-vortex array) on Jupiter's south pole, and an array of eight vortices rotates around a central vortex (a "8 + 1"-vortex array) on its north pole.The exact nature and the vertical structure of the vortices over the poles of Jupiter are still debated, but elements of observations may suggest that the flow regime is similar to the Surface Quasi-Geostrophic (SQG) regime.Our purpose is not to explicitly model the Jovian polar vortex arrays but rather to explore the more general case of finite-core QG vortex arrays.Following Reinaud (2019), 21 we determine and characterize regular arrays of finite volume and uniform potential vorticity vortices in mutual equilibrium, and we address their linear stability.We focus on "5 + 1" and "8 + 1"-vortex arrays, as they provide generic examples of typical "m + 1"-vortex arrays.It should be noted that Reinaud (2019) only considered finite-core vortex arrays up to m = 7 but analyzed the stability of point vortex arrays up to m = 8.It should also be noted that within the point vortex approximation, the dynamics of co-planar QG vortices and SQG vortices are formally equivalent. 25Hence, the stability of SQG point vortex arrays can be directly recovered from the QG analysis. 21On the other hand, for finite-core vortices, an SQG vortex can be seen as the limit of QG finite-volume vortex whose height-to-width aspect ratio tends to zero. 26,27e show that the vortex arrays are sensitive to deformation modes when the vortices are close enough together.The impact of the instability is then investigated by numerically simulating the nonlinear evolution of a selection of unstable vortex arrays.Vortices first deform resulting in the partial merger of some of them.These mergers break the symmetry of the vortex arrays and eventually lead to a convoluted late evolution of the vortices.
This paper is organized as follows: The governing equations and the setup of the study are presented in Sec.II.Section III describes the equilibrium states for the vortex arrays and addresses their linear stability.The nonlinear evolution of a selection of unstable vortex arrays is presented in Sec.IV. Conclusions are drawn in Sec.V.

A. The quasi-geostrophic model
The quasi-geostrophic model used in this study derives from the three-dimensional Euler's equations for a rotating and stratified fluid under the Boussinesq approximation.Let U be a horizontal velocity scale, L a horizontal length scale, and H a vertical length scale.We denote f as the Coriolis frequency and N as the buoyancy (or Brunt-Väisälä) frequency, both assumed constant for simplicity.The flow regime in a rotating stratified fluid is characterized by two non-dimensional numbers, the Rossby number Ro = U/( f L) and the Froude number Fr = U/(NH).For Fr 2 ≪ Ro ≪ 1, an asymptotic expansion in Ro of the equations leads to the definition of the quasi-geostrophic potential anomaly (PV), q, as where φ = p/( f ρ 0 ) is the geostrophic stream function, p is the pressure, and ρ 0 is the mean reference density.The PV q is materially conserved through a layer-wise two-dimensional geostrophic velocity (u, v, 0), where Equations ( 1)-( 3) form a closed system.In Eq. ( 1), the physical vertical coordinate zp has been replaced by z = zpN/ f , where the constant ratio N/ f ≫ 1 in large parts of the atmosphere and oceans.It should be noted that in the QG regime, the vertical velocity w is not strictly speaking zero but too small to contribute to the advection of q.It, therefore, does not contribute to the dynamics but can be evaluated as a diagnostics.We will not discuss w in this paper.Equation ( 1) can be formally inverted using the appropriate Green's function, in an infinite fluid domain.Equation ( 4) can be formally differentiated to obtain the advection velocities u and v.The fluid domain occupied by the vortices in the vertical direction is discretized by n l horizontal layers of equal thickness, on which the Green's functions, for φ, u, and v, can be analytically integrated.For uniform PV vortices, the remaining horizontal integrals can be converted into contour integrals on the vortex boundary using Green's theorem. 28t this level of approximation, the flow is completely determined by the evolution of the materially conserved q.
The methods used to determine finite-core equilibrium vortex arrays, to address their linear stability and their nonlinear evolution, are all based on this contour representation of uniform PV vortices, and they all compute the necessary quantities by explicit contour integration, using local two-point Gaussian quadrature.All these approaches are Lagrangian in nature, and the fluid domain is explicitly unbounded.
The method used to determine the equilibrium states is described in Refs.21 and 29.It is based on making the vortexbounding contours converge to streamlines in the reference frame rotating with the equilibrium, i.e., a constant, along each horizontal vortex-bounding contour Ci, where φ is the stream function expressed in the frame rotating with the array.It uses an iterative method based on partial linearization of Eq. ( 5).The method was first used for two-dimensional vortices. 30he linear stability analysis is performed by analyzing deformation modes along the vortex-bounding contours, including a mode corresponding to the displacement of the full contour.
The nonlinear simulation is performed by advection explicitly the vortex-bounding contours using a fourth-order Runge-Kutta scheme.A contour surgery procedure 31 is used to control the accuracy and complexity of the contour description of the vortices.

B. General geometry
Each array consists of m = 5 and 8 peripheral vortices whose centroid lies on a circle and an additional vortex lying at the center of the array, which is set to the origin of the frame.Each vortex has the same uniform PV, q = 2π, the same volume V, and the same height-to-width aspect ratio = h/r h , where h is the half-height of the vortex and r h its mean horizontal radius.All vortices are centered on the plane z = 0.It should be noted that since the vertical direction has been rescaled by the constant ratio N/ f , if hp is the physical halfheight of the vortex (h = hpN/ f ), the vortex height-to-width aspect ratio is where Rov and Frv are the Rossby and Froude numbers based on the vortex dimensions and Bu is the vortex-based Burger number.
When determining the equilibrium states, symmetries are imposed to minimize the computational cost.First, we take advantage that the natural symmetry with respect to the plane z = 0, and corrections to the contours are only calculated on the upper half of the vortices.These corrections nevertheless depend on the influence of the full vortices explicitly.We also use the fact that the system has an m−fold azimuthal symmetry.Hence, corrections to the contours are calculated on 1/m th of the upper central vortex only and on only one of the m peripheral vortices.Finally, we also use the fact that the peripheral vortices are symmetric with respect to a vertical plane passing through the vortex centroid and the origin.This means that corrections to the contours are only determined on one-quarter of one of the peripheral vortices.On the other hand, no symmetry is explicitly imposed in the linear stability analysis or in the nonlinear simulations.
Vortices are mapped by n l = 83 horizontal layers.Contours bounding the vortices in each layer are discretized using np = 4m(n l /m) nodes, where (n l /m) is an integer division.This offers a high horizontal resolution of contours while ensuring the number of nodes is even and divisible by m.This is used in practice to impose the symmetries discussed above.We determine branches of equilibria.A branch of equilibria for a "5 + 1"-vortex array (respectively, "8 + 1"-vortex array) corresponds to a set value of .Equilibria within a given branch differ by the distance between the vortices.For both m = 5 and 8, we start each branch with a circular vortex array of typical radius ℓ = O (5 or 6 r h ).Then the iterative procedure is used to find the corresponding equilibrium.We use spheroidal vortices as a first guess for the shape of the vortices.This choice is justified because the vortices are well separated for ℓ ≃ 5, 6r h, and a spheroidal (axisymmetric) vortex standing alone is a steady solution of the QG equations.During the iterative procedure, the distance δi between the origin and the innermost edge of the peripheral vortices for m = 5 or the distance δo between the origin and the outermost edge of the peripheral vortices for m = 8 is fixed.These distances can be varied monotonously along each respective full branch of equilibria and thus can be used as a parameter to describe the full branch.The iterative procedure computes the corrections to the vortex-bounding contours as well as the correction to the angular velocity Ω at which the equilibrium steadily rotates.The iterative procedure is stopped when the correction to Ω becomes less than 10 −11 .Then the distance between the origin and the peripheral vortices is decreased by a small amount, δc, and the procedure is resumed for the next state along the branch.δc is not constant but is progressively decreased as the distance δi (respectively, δo) becomes small.In practice, we first set δc = O(0.1 r h ) and start to determine the branch of equilibria.When the method fails to converge, we restart the procedure from the last converged equilibrium while dividing δc by ten.This is necessary because the vortices at equilibrium become increasingly deformed as they are closer together, hence the corrections required between two equilibria increase as we decrease δi (respectively, δo).We stop the procedure when vortices are close enough to nearly touch or when the vortices form sharp corners, indicating the end of the branch of equilibria.For each steady state, we determine the location of the centroid of one of the peripheral vortices, which allow us to determine the radius ℓ of the circular vortex array.For both m = 5 and m = 8, we first consider = 0.1 and then from = 0.25 with increments of 0.25 in the aspect ratio.We were not able to determine numerically the full branches of equilibria for /r h > 2 for m = 5 and > 2.5 for m = 8.We label vortex 1 the peripheral vortex whose centroid lies on the semi-axis x > 0.

III. EQUILIBRIUM VORTEX ARRAYS
We first characterize the shape of the "5 + 1"-and "8 + 1"vortex arrays.We present examples of vortex arrays for two branches of equilibria in Fig. 1.The first branch of equilibria is for the "5 + 1"-vortex arrays with = 1.25, while the second branch corresponds to the "8 + 1"-vortex array with = 1.In both cases, as the radius of the array is decreased, the vortices are increasingly deformed.As already observed, vortices deform the most toward the nearest vortex. 18If m = 5 < 2π, the closest vortex to a peripheral vortex is the central vortex, and the innermost edge of the peripheral vortices points toward the center.This explains why δi decreases monotonously along the branch for the "5 + 1"-vortex arrays.On the other hand, for m = 8 > 2π, the closest vortices of given peripheral vortices are the two neighboring peripheral vortices in the circular array.The peripheral vortices develop a sharp edge in the azimuthal direction.A simple way to illustrate this deformation is to measure the radial extend Δv = δo − δi of the peripheral vortices, defined as the distance between the innermost and the outermost edges of the peripheral vortices.Results are presented in Fig. 2. For m = 5, we see the increase of Δv as ℓ decreases, indicating that the peripheral vortices elongate in the radial direction, pointing more and more toward the center of the array as the vortices get closer together.This trend accelerates near the end of the branch of equilibria, showing the increased vortex deformation as the strain induced by the vortices increases.The opposite is observed for m = 8, where the vortices flatten in the radial direction.Again, the deformation rapidly increases near the end of the branch of equilibria.This indicates that the vortices elongate in the azimuthal direction (the volume of the vortices remains the same along a branch of equilibria).Overall, we see that, for a given ℓ/r h , Δv/r h increases as increases for m = 5, while it decreases for m = 8.In both cases, Δv/r h increasingly differs from 2; the value for a circular horizontal vortex cross section as increases, indicating an increased deformation as increases for a given ℓ/r h .The angular velocity Ω at which the array steadily rotates is shown in Fig. 3.At leading order, Ω scales as Ω ∼ Q/ℓ 3 , where Q is the volume integrated PV of the vortices.Hence, Ω ∼ qr 2 h h/ℓ 3 = q (ℓ/r h ) −3 , explaining qualitatively the
trends observed: Ω increases as ℓ/r h decreases as the vortices interact more strongly, and Ω increases with for a given ℓ/r h .
pointing toward the neighboring peripheral vortices on either side, while their innermost side flattens.
For both m = 5 and m = 8, however, the most striking difference between the shape of the peripheral vortices as is varied in the vertical direction.Here, the PV distribution is a function of the vertical coordinate z; hence, the velocity field induced by the vortices also depends on z.This generates a vertical shear.We see that to remain in equilibrium and withstand the shear they create, the vortices bend.This is not unlike the deformation observed during the early nonlinear evolution of non-equilibrium pairs of interacting spheroidal vortices as they attempt to adapt to a nearby equilibrium. 32,33For h/r h ≤ 1, the outermost side of the peripheral vortices is convex while it becomes concave for h/r h > 1, near the end of the branches of equilibria.In the latter case, it should be nonetheless noted that the outermost side is still convex for large δi (respectively, δo).The vertical shape of the peripheral vortices shown in 6 where the vertical cross section of vortex 1 passing through its centroid is displayed in the − x * )/r h z/h)-plane.Here, x * is the midpoint between the innermost and outermost edges of the vortex.The figure also allows us to visualize the sharp innermost corner for m = 5 and the flattened innermost side for m = 8.
We next address the linear stability of the finite-core vortex equilibrium arrays.The linear stability analysis consists of the analysis of deformation modes of the vortex-bounding contours, taken proportional to exp σt, where σ = σr + iσi ∈ C. The real part σr of σ is therefore the mode growth rate while σi is its frequency.The detail of the approach is provided in Refs.21 and 29.Reinaud (2019) 21 has shown that both the arrays of "5 + 1" and "8 + 1" identical point vortices are linearly (neutrally) stable.This stability is independent of the distance between the vortices, and only displacement modes of instability can be studied within the point vortex approximation.It is therefore anticipated that the finite-core vortex arrays are linearly stable when the vortices are far apart, and in that case, well modeled by point vortices.The maximum growth rates σr are plotted vs ℓ/r h in Fig. 7.For m = 5, the equilibria are stable for large ℓ/r h , and a first mode of instability appears for a critical distance lc = (ℓ/r h )c.The distance lc increases with the vortex aspect ratio indicated that tall vortices are more sensitive to instabilities.This is related to the increased influence of the vertical shear as increases.The situation is similar for m = 8 and ≥ 2. For m = 8, smaller , and large ℓ, the numerical linear stability captures a small non-zero growth rate, slowly growing from large ℓ as ℓ decreases before the growth rate increases sharply.The small growth rates before the sharp increase are believed to be spurious.They are believed to result from numerical truncation errors and the fact that the equilibria are not exact but obtained numerically within a given accuracy.The latter is related to the finite spatial discretization of the vortex boundary and the tolerance on the correction to Ω used in the iterative method to obtain the equilibria.This may result in a residual unsteadiness in the "equilibria."Numerical nonlinear simulations, however, suggest that these equilibria (for ℓ large) are in fact stable.The sharp increase in σr, however, indicates the onset of a genuine instability.
The instability consists of deformation modes that primarily affect the regions of the vortices where sharp corners are forming: that is the innermost corner of the peripheral vortices for m = 5 and the side-corners for m = 8 as shown in Sec.IV, which addresses the nonlinear evolution of a selection of unstable equilibria.

IV. NONLINEAR EVOLUTION
We next perform nonlinear simulations of unstable equilibria using Contour Surgery. 28The initial conditions consist of the equilibria obtained in Sec.III.The modes of instability are not forced, but the perturbations triggering the instability arise from the numerical noise.The latter is mostly associated with redistributing the nodes on the vortex-bounding contours, decreasing their number while retaining an accurate description of the vortex boundary.In practice, the numerical method used to obtain the equilibria typically requires a high spatial resolution, beyond what is needed and practical in the nonlinear simulations.Reducing the initial number of nodes discretizing the vortex-bounding contours allows us to speed up the initial stages of the nonlinear simulation.In all cases, the time scale is implicitly defined by the vortex PV, q = 2π.For comparison, a single sphere of PV has a turnover period T = 6π/q.The length scale of the problem is set by the vortex half-height, h = 0.5.
We first consider the "5 + 1"-vortex array for = 0.25 and δi/r h = 2.583 corresponding to the last equilibrium of the branch.Recall that the last equilibrium of the branch is the most unstable one.The shape of the vortices is shown in Fig. 8 I.
in the evolution, and one of the peripheral vortices shares more material with the central vortex compared to the other peripheral vortices.This further increases rapidly the asymmetry of the flow.The later evolution is convoluted with more of the peripheral vortices merging.A large compound structure forms at the center, which exerts further complex shear and strain on itself and the surrounding vortices and filaments.The evolution is generic of the interaction for ≤ 1 and m < 7. 21 We next consider the last "5 + 1"-vortex array for = 2 and corresponding to δi/r h = 1.371.Results are shown in Fig. 9 at t = 12, 5, 15, and 25.The initial phase of the interaction is similar to the one in the previous case.The innermost edges of the peripheral deform, FIG. 5. Shape (vortex-bounding contours) of the equilibrium "8 + 1"-vortex arrays for the last state numerically obtained for indicated to the left of the images.Two views are presented: a top view (0 ○ from the vertical axis) and a view at 60 ○ from the vertical axis.The corresponding values for ℓ/r h are given in Table I   = 2 (solid red), = 2.25 (solid black), and = 2.5 (dashed yellow).The curves for = 0.1 and = 0.25 are almost undistinguishable.For each case, the corresponding value of ℓ/r h is given in Table I.

FIG. 7.
Maximum growth rates σr /q vs ℓ/r h for (a): "5 + 1"-vortex arrays and (b): "8 + 1"-vortex arrays and for = 0.1 (dashed green), = 0.25 (solid green), = 0.5 (dotted blue), = 0.75 (dashed blue), = 1 (solid blue), = 1.25 (dotted red), = 1.5 (dashed red), = 1.75 (solid red), = 2 (dotted black), = 2.25 (dashed black), and = 2.5 (solid black).shedding filaments of PV that wrap around and merge with the central vortex.In the QG regime, the angular impulse, is conserved.The migration of PV toward the center as the innermost edges of the peripheral vortices merge with the central vortex must be compensated by the migration of some PV away from the center.In the initial equilibrium vortex array, the top and bottom of the peripheral vortices are horizontally further away from the center, due to the vortex curvature.This is accentuated during the nonlinear evolution by the requirement of conserving J as the middle section of the peripheral vortices interacts strongly with the central vortex.
As the top and bottom of the peripheral vortices move from the center of the they are influenced by the PV accumulating near the center and rotate more slowly than the center.The differential rotation stretches the vortex top and bottom into long filaments, which spiral around the merged structure in the center as seen on the right panel of Fig. 9.We next turn our to selected "8 + 1"-vortex arrays.We again start by considering an array of very oblate vortices with = 0.25.We consider last state along the equilibrium branch corresponding to δo/r h = 4.408.Results are shown in Fig. 10.The instability results in the deformation of the peripheral vortices, in the region of the side corners where the peripheral vortices are close to their neighbors.Some of the vortices merge with their neighbors, forming structures.Similar results for = 0.5 (not shown) show that such structures can break and reform, resulting in a convoluted late evolution and the build-up of asymmetries.
This can also be seen in the next example for = 1.Again, we consider the last and most unstable equilibrium, corresponding to δo/r h = 4.528.Results are shown in Fig. 11.Mergers of some peripheral vortices create large compound structures.These large asymmetric structures surround and shear the central vortex, which elongates and is eventually absorbed by merger by one of the peripheral compound structures.The late evolution is convoluted and resembles a chaotic motion of a collection of vortices of various sizes and shapes, which occasionally strongly interact.
The final example presented is for the tallest "8 + 1"-vortex array with = 2.5.We consider again the most unstable state with δo/r h = 4.635.Results are presented in Fig. 12.As in the two other "8 + 1"-vortex arrays discussed above, the instability results in the deformation of the side corners of the peripheral vortices and the merger of some of them their neighbors.Like the previous example, one of the merged peripheral structures gets close to the central vortex and merges with it.As the flow evolves, filaments are shed away from the center to conserve angular impulse J, mostly the vortices top and bottom.

V. CONCLUSION
We have determined and studied branches of equilibria for three-dimensional quasi-geostrophic vortices "m + 1"-vortex arrays for m = 5, 8 and for various vortex aspect ratios .It should be noted that as decreases, the dynamics becomes similar to the dynamics of surface density anomalies of the surface quasi-geostrophic theory.The vortex aspect ratio is related to the vortex, hence the flow's, Burger number.Tall vortices tend to bend to reach equilibrium as increases to withstand the vertical shear they induce onto each other.When the vortices are far apart, the vortex arrays are linearly neutrally stable.When the distances between the vortices are less than a threshold, depending on the number m of peripheral vortices and the vortex height-to-width aspect ratio (or equivalently depending on the vortex Burger number), the vortex arrays are linearly unstable.The threshold is, however, typically close to the end of the branch of equilibria, meaning that, in practice, most vortex arrays are stable.This threshold, measured as a critical value of ℓ/r h for the onset of instability, increases as increases.The instability deforms the vortices in the regions where they are the closest to a neighboring vortex.This corresponds to the region of the vortices subjected to the strongest strain.This means that for m = 5, the instability primarily affects the innermost edges of the peripheral vortices and precipitates their merger with the central vortex.On the other hand, for m = 8, the instability triggers the merger of some peripheral vortices with their neighbors along the circular array.This re-organization of the peripheral vortices breaks the symmetry of the configuration.Eventually, the peripheral vortex structures can merge with the central vortex.For both m = 5 and m = 8, tall vortices bend, with their middle section closer to the center of the array than their top and bottom sections.During the nonlinear evolution, spiraling filaments are shed from the top and bottom of the peripheral vortices.
A natural extension to this work is to investigate threedimensional finite-core staggered vortex arrays.A staggering array can be seen as a combination of two concentric circular vortex arrays.For example, the vortices of the "8 + 1"-vortex array observed on the north pole of Jupiter are in fact slightly staggered.The staggering of the peripheral vortices is a natural mode of instability of circular vortex arrays with an even number of peripheral vortices. 21,34taggered point vortex arrays have already been studied, 20 but it would be interesting to determine finite-volume staggered vortex array and to investigate whether an unstable circular vortex array can re-organize itself into a stable staggered array.

FIG. 4 .
FIG. 4.Shape (vortex-bounding contours) of the equilibrium "5 + 1"-vortex arrays for the last state numerically obtained for indicated to the left of the images.Two views are presented: a top view (0 ○ from the vertical axis) and a view at 60 ○ from the vertical axis.The corresponding values for ℓ/r h are given in TableI.

TABLE I .
Value of ℓ/r h for the last equilibrium states found along the branch and shown in Figs.4 and 5.