Polar vortex crystals: Emergence and structure.

Vortex crystals are quasiregular arrays of like-signed vortices in solid-body rotation embedded within a uniform background of weaker vorticity. Vortex crystals are observed at the poles of Jupiter and in laboratory experiments with magnetized electron plasmas in axisymmetric geometries. We show that vortex crystals form from the free evolution of randomly excited two-dimensional turbulence on an idealized polar cap. Once formed, the crystals are long lived and survive until the end of the simulations (300 crystal-rotation periods). We identify a fundamental length scale, Lγ=(U/γ)1/3, characterizing the size of the crystal in terms of the mean-square velocity U of the fluid and the polar parameter γ=fp/a2p, with fp the Coriolis parameter at the pole and ap the polar radius of the planet.

The Juno spacecraft revealed a distinctive form of organized turbulence in the polar atmosphere of Jupiter (1). The North Pole exhibits a central polar cyclone surrounded by eight circumpolar cyclones; the South Pole is characterized by five cyclones surrounding a central polar cyclone. These cyclones are organized into a vortex crystal: a symmetric array of strong vortices in solid-body rotation within a background of weaker vorticity. The Jovian vortex crystals are located within 8 of the poles and have endured since their discovery in 2016 with very little change to their overall structure (2, 3). Jovian crystals are approximately stationary in the System III reference frame (4). The essential physical mechanisms responsible for the emergence and persistence of this peculiar polar dynamical regime are poorly understood.

Although observations of polar vortex crystals began with Juno, the most basic model of this phenomenon dates back to Kelvin’s 19th century theory of vortex atoms, i.e., the dynamics of equally spaced point vortices arranged in concentric rings (5–7). Laboratory experiments with magnetized pure-electron plasmas are another antecedent (8–11). This experimental system is isomorphic to a near-ideal two-dimensional (2D) fluid contained within a circular domain, with electron density equivalent to vorticity. Through vortex nucleation, merger, and finally self-organization, regular vortex arrays, in solid-body rotation (9), spontaneously crystallize from a filamented initial vorticity field. These electron-plasma experiments are simpler than polar planetary dynamics in several respects:

There is no analog of the variation of the Coriolis parameter with latitude (the β-effect);

The vorticity has strictly one sign (electrons are negatively charged);

The plasma is contained within a circular domain with a free-slip boundary;

Plasma dynamics are “barotropic” (the deformation length is infinite); and

There is effectively no dissipation and there is no energy source (small-scale convection) required to sustain the vortices.

Nonetheless many authors have remarked on the striking resemblance between electron-plasma vortex crystals and Jovian observations (1, 4, 12, 13). While these two systems have many features in common, it is crucial to relax simplifications 1 and 2 in any model of polar atmospheres. In particular, the small variation of the Coriolis parameter in the vicinity of the poles is key.

Recent observations show that moist convection drives an upscale energy transfer at Jovian high latitudes (14), consistent with the regime of rapidly rotating Rayleigh–Bénard convection (15). The barotropic component of the flow (i.e., the circumpolar and polar cyclones) exhibits kinetic energy spectra with a spectral slope, consistent with quasi-geostrophic (QG) dynamics.

These results motivate the hypothesis that barotropic QG dynamics, specialized to an idealized “polar-cap” geometry, might explain basic features of the polar planetary regime. Here we show that useful dynamical information, such as the radius of the polar vortex crystal, follows from dimensional analysis of the polar-cap QG dynamics, prepared with an initial condition having a length scale characteristic of convective cells (14).

A Barotropic QG Model of Polar-Cap Dynamics

An important characteristic of the QG polar-cap regime is that planetary parameters occur only in the combination , where f is the Coriolis parameter at the pole and a is the polar radius of the planet. The parameter γ is significant because close to the poles the Coriolis parameter is , where r is the distance from the pole. We discuss the North Pole so that f is positive and cyclones are positive vortices.

The barotropic polar-cap QG system is

where is the stream function, is the relative vorticity, and q is the potential vorticity (PV). We employ a Cartesian coordinate system centered on the pole so that . For the planetary PV η in [], it is instructive to compare two configurations. The first one is the polar-cap PV

where is the indicator function of the disk with radius ; i.e., and . Inside the specified radius is equal to the planetary PV, . The second configuration is the “flat-trap” PV

The flat-trap model , with no planetary PV gradient within the trap, is close to the electron-plasma experiments. Both and have a discontinuity of strength at the trap radius . This discontinuity results in an enclosure that largely confines the QG turbulence within the disk ().

The model is solved as an unforced initial value problem using a doubly periodic pseudospectral method provided by GeophysicalFlows.jl (16); the domain size is . We use a trap radius that is as large as possible, while still ensuring that there is minimal interaction between neighboring periodic domains. The vortex trap is an alternative to the dissipative sponge layer used in earlier studies of polar-cap dynamics (17–19). The trap approximately conserves energy throughout the 8-y duration of the simulation. For further numerical considerations, and discussion of the trap boundary, see .

The planetary radius a does not appear in the formulation above (nor in the initial condition below). Therefore, just as the planetary radius a is irrelevant for midlatitude β-plane dynamics, a is not a relevant length scale for polar-cap dynamics. The irrelevance of a can be appreciated intuitively because polar vortex crystals are confined in regions that are within of both poles. These caps occupy less than 1% of the total area of Jupiter. Within the QG framework, the Coriolis parameter f is also irrelevant: f and a are inextricably combined into the only relevant planetary parameter .

Initial Condition.

The initial relative vorticity is a random monoscale field, multiplied by a taper that confines ζ0 to the center of the trap: See Fig. 1 and for further details of the initial condition. The untapered initial relative vorticity is concentrated within a narrow annulus in wavenumber space; the central radius of the annulus is and the width of the annulus is much less than . In addition to , the initial condition provides a velocity scale U defined by

Fig. 1.

(A–I) Evolution of the relative vorticity ζ of the reference run with parameters in []. The central circle in each panel has radius km; in G–I the radius of the crystal is very close to . See Movies S1–S3 and showing the evolution of .

where is a domain average and ψ0 the initial stream function. The kinetic energy is approximately conserved throughout the evolution and thus U is the characteristic flow velocity. (More precisely, 5% of the initial kinetic energy is lost in the first 19 h of evolution and another 5% in the remaining 8 y; see and Fig. S1 for details.)

The Intrinsic Polar-Cap Length Scale.

From U and γ one can form the length.

is the polar analog of the Rhines length (20) on a midlatitude β-plane. To develop this analogy in more detail, note that the magnitude of the planetary PV gradient is and therefore the Rhines length at a distance r from the pole is

Well-known arguments identifying in midlatitudes all assume that the curvature of the planetary PV is negligible over the scale itself, i.e., that one can approximate the planetary PV by a linear function of latitude over a distance . Eq. shows that this assumption must fail as one approaches the poles: Curvature becomes important. We proceed heuristically to identify the new relevant length that replaces as : Equating r to in [] gives . Solving for then results in the length on the left of []. Thus is identified as an intrinsic polar length scale. We show below that is the radius of emergent QG polar vortex crystals.

Disclaimer.

Despite the adoption of Jupiter-like numbers in [] below, we do not claim quantitative Jovian accuracy for this barotropic QG model. For example, we use barotropic dynamics so that the deformation length is infinite (as it is in the plasma experiments). Observations summarized in ref. 12 indicate that the at the Jovian poles is somewhere in the range 350 to 1,300 km. This range encompasses various polar length scales, such as the vortex diameter defined as the radius of maximum azimuthal velocity ( to km). Thus baroclinic effects associated with finite are likely to be quantitatively important at the Jovian poles. While an unforced and weakly dissipative model cannot claim strict application to the continuously forced Jovian system, it does reveal the dominant evolutionary tendencies resulting from inertia and planetary PV gradients. In other words, the barotropic QG model, with the addition of γ, is sufficient to capture the key vortex-dynamical processes responsible for the genesis of polar vortex crystals.

Results: The Polar Cap

Fig. 1 shows a solution of [] and [] with in []. This “reference run” uses Jupiter-like parameters

corresponds roughly to the length scale of the shallow atmospheric convective cells (14). We used U = 80 m ⋅s−1 based on the latest observations of Jovian polar winds (14). We acknowledge that a lower value might have been closer to the mean-square velocity in Jovian high latitudes. However, a dimensional analysis in shows that all velocities in the polar-cap QG system are proportional to the external parameter U and thus it is easy to make adjustments.

We strive to achieve a large value of , so that the vortex crystal with radius (see below) is well separated from complications at the periphery of the trap. With available computational resources we use (see for a discussion of sensitivity to ). The PV discontinuity at is s−1. The domain length scale is . The initial condition in Fig. 1 shows random small-scale () vorticity concentrated in the center of the trap, with PV extrema s.

The development of a vortex crystal from random initial conditions is shown in detail in Movies S1–S3. The process is summarized by nine snapshots of the relative vorticity ζ in Fig. 1. There are several distinct evolutionary stages corresponding to the three rows of Fig. 1.

Vortex Nucleation.

The nucleation of small-scale vortices from the random initial condition in Fig. 1 is underway in Fig. 1 and, after 7 d of evolution, is almost complete in Fig. 1. Nucleation produces roughly equal numbers of cyclones and anticyclones (Fig. 1). Consistent with assumptions of vortex-gas scaling theory (21–23), the PV extrema in Fig. 1 form the cores of the nucleated vortices in Fig. 1 and C. Some of these vortices pair to form rapidly propagating dipoles, resulting in occasional collisions with the trap boundary at . Most freshly nucleated vortices, however, remain within the trap.

Vortex Segregation and Merger.

Subsequent evolution from the state in Fig. 1 involves merger (24–26) of like-signed vortices to form larger vortices and vortex segregation guided by γ. Vortex segregation refers to the systematic migration of cyclones toward the poles and anticyclones toward the equator (17, 18, 27–29). Segregation is noticeable in Fig. 1 and is largely complete after 75 d in Fig. 1 The pole is occupied by cyclones and anticyclones have accumulated up against the trap boundary at .

In view of PV mixing arguments, invoked below in connection with [], it is interesting that vortex segregation is “antidiffusive”: Cyclones () move poleward and reinforce the initial r = 0 maximum in the planetary PV . And the outward radial motion of anticyclones () reinforces the initial minimum in at .

Emergence of the Crystal.

After segregation is complete there are eight large cyclonic vortices moving chaotically within a cap of radius centered on the pole; the snapshot in Fig. 1 is representative of this state. The configuration in Fig. 1 is not a crystal; the vortex octet is not in approximate solid-body rotation. Moreover, the cyclones compete for occupancy of the central position—one cyclone will command the center until jostled out of position by a colleague (Movie S1). This jostling happens repeatedly until the configuration in Fig. 1 is finally achieved after about 1.5 y of evolution. The central cyclone in Fig. 1 remains in place until the end of the run in Fig. 1 at 8 y.

Once the permanent central cyclone is established, the octet crystallizes into solid-body rotation (Movie S3). The solid-body rotation is not perfect: Vortices have small epicyclic oscillations about their mean position, reminiscent of the oscillations of the cyclones about their mean position at the South Pole of Jupiter (3, 4). The oscillations may result from differences in strength between the eight cyclones; e.g., the smallest cyclones have the largest excursions.

Zonal-Mean Flow, Cooling, and PV Homogenization.

Although the polar cyclones in the lowest row of Fig. 1 are a prominent feature of the flow, one should not ignore the zonal-mean flow and the associated zonal mean relative vorticity . (See for a summary and discussion of zonal mean dynamics.) The eight cyclones interact with one another by mutual advection, but they also interact with the zonal-mean flow . For example, and in analogy with the electron-plasma experiments, we argue that the chaotic state in Fig. 1 is “cooled” as the irregular motion of the polar cyclones mixes the zonal-mean PV, . The vortex crystal is formed as the chaotic motion of the circumpolar cyclones is damped by transfer of angular momentum and energy to (10).

Following arguments along the lines of Taylor (30) and Bretherton (31), a main result of the zonal-mean equations is that the zonal-mean acceleration, , is related to radial vorticity transport, , by

In the central circle chaotic vortex motion directly results in and mean-flow acceleration. But there is also far-field and, through [], far-field zonal-mean acceleration. (The “far field” of the circumpolar cyclones is defined by .) Chaotic vortex motion within the circle produces a fluctuating far-field radial velocity, , corresponding to the unsteady dipole moment of the circumpolar cyclones. There is also an azimuthal component, , resulting from the constant-in-time total circulation. But only the weaker radial component, , results in radial transport and therefore far-field mean-flow acceleration via [].

Mixing is so strong that by the end of the reference run there is striking PV homogenization (32). Fig. 2 shows that the PV is mixed into discrete steps (33, 34). Thus, the vortex crystal coexists with an axisymmetric version of the β-plane PV staircase. Because of the axisymmetric geometry, this homogenized background PV is more accurately described as a “PV terrace.” Once the PV terrace is formed, γ-driven vortex segregation is no longer operative: The planetary γ-effect is annulled by the relative vorticity of the mean azimuthal flow . And there is no further mean-flow acceleration because ; i.e., PV fluctuations cannot be created by stirring homogeneous PV.

Fig. 2.

(A) PV as a function of along the diagonal cuts show in Fig. 4. (B) Azimuthally averaged azimuthal velocity, as a function of . The “numerical” curve is calculated from the final state (8 y) of the reference run; the “theoretical” curve is inferred from an analytic fit to the azimuthally averaged PV. The fit uses Gaussian polar cyclones embedded in a background of homogenized PV and does not attempt to model cyclonic vorticity outside the trap where (see for details).

Fig. 4.

The PV, , of six final states. (A) The reference solution shown in Fig. 1. Fig. 2 shows the PV along the diagonal sections in A. (B–F) The other five examples differ from the reference solution in A only by changing the seed in the random number generator used to construct the initial condition. We show the PV, , rather than the relative vorticity (e.g., ζ in Fig. 1) to make the PV steps visible.

Because of PV homogenization the eight cyclones in Fig. 1 sit on top of a flat PV terrace. Throughout the first PV terrace in Fig. 2, is less than the planetary vorticity, . Thus, the background relative vorticity, , is negative, i.e., opposite in sign to the circumpolar cyclones. A background of homogeneous weaker vorticity, of opposite sign to that of the polar cyclones, is an observed feature of the Jovian poles (14, 35).

Fig. 2 shows a comparison of the azimuthally averaged azimuthal velocity, , of the reference run with an analytic fit based on homogenized PV steps and Gaussian vorticity profiles for the polar cyclones (see for details of this PV terrace model). Although throughout the first terrace, the azimuthal velocity in Fig. 2 is dominated by the irrotational far-field velocity of the central cyclone and is therefore strongly positive. The secondary local maximum in at is the signature of the seven circumpolar cyclones.

PV terraces emerge in all our unforced simulations. But terrace formation in a continuously forced-dissipative system is less certain: We do not claim that the high-latitude Jovian atmosphere is PV terraced.

Anatomy of the Polar Cyclones.

Fig. 3 shows that the emergent polar cyclones are unshielded; i.e., in Fig. 3 the vorticity has one sign (positive) and each cyclone has a nonzero circulation. Therefore, the far-field azimuthal velocity around each polar cyclone is . (Here is the circulation of the cyclone; r is the distance from the vortex center; and “far-field” means that r is much greater than the vortex radius, defined as the radius of maximum velocity in Fig. 3.) Slow algebraic decay, , differs qualitatively from the exponential decay of u assumed in recent models: Those models employ completely shielded vortices with (12, 13).

Fig. 3.

(A) Zoom-in on the PV field of the reference-run polar vortex crystal. This crystal is composed of eight cyclones (numbered). The black circle has radius . (B) Diagonal cuts through the center of four of the cyclones (numbered) showing relative vorticity (r indicates the distance from the vortex center). The solid line is one diagonal cut and the dashed line is an orthogonal cut. The two diagonal cuts overlap almost perfectly, highlighting the axisymmetry of the cyclones. (C) Azimuthal velocity of the eight cyclones as a function of distance from their center, r, divided by . The thick curve is the average azimuthal velocity of the eight cyclones; the average profile has a maximum azimuthal velocity at radius or 500 km. In C the slow decrease of the far-field azimuthal velocity is evident.

In Fig. 3 the maximum cyclone velocity is between 5 and 10 × U, where U is the root-mean-square velocity defined in [] and set to 80 m ⋅s−1. The velocity 80 m ⋅s−1 is more typical of the maximum cyclone velocity of Jovian cyclones—see figure 6 of ref. 36. Therefore, the cyclones in Fig. 3 have peak velocities very much faster than those of Jovian polar cyclones. All velocities in the polar-cap system are, however, proportional to U (). Thus we can make an a posteriori adjustment so that the maximum cyclone velocity in the polar-cap model matches the Jovian observations: Reduce U by a factor of 8, which reduces by a factor of 2. With this reduction in U, the polar cyclones now have typical Jovian velocities, but the polar-cap crystal has a radius km, which is somewhat smaller than the 8,700-km radius of the Jovian crystals (36).

In addition to being too fast, the cyclones in Fig. 3 are also too small relative to those of Jupiter. The average cyclone radius indicated in Fig. 3 is about or, with U = 80 m ⋅s−1, a cyclone radius of 500 km. Jovian polar cyclones range in radius from 2,000 to 3,500 km (36).

In summary, if we use U = 80 m ⋅s−1, so that the crystal radius is km, then the reference-run cyclones are much faster, and much smaller, than those of Jupiter. We expect that these results depend on the use of a barotropic model with infinite deformation radius.

Sensitivity to Initial Conditions.

Although the emergence of as the radius of the crystal is robust and reproducible, the exact number of vortices in the crystal is an accident of initial conditions. For example, with the reference values in [] fixed, if we change the seed in the random number generator responsible for initializing ζ0, we produce a diverse set of final configurations; see Fig. 4 for six examples. In Fig. 4 there is a dipole consisting of a large polar cyclone and a smaller cyclone; this dipole orbits the pole. In Fig. 4, a small anticyclone is paired with the large polar cyclone, reminiscent of the Jovian observations. Small anticyclones are also seen in between the large cyclones at the North Pole of Jupiter; see figure 1 of Adriani et al. (1). In Fig. 4 there is a lone polar cyclone. Electron-plasma experiments exhibit an analogous sensitivity to initial conditions (8, 9).

The structure of the zonal-mean PV is less sensitive to initial conditions: The six runs in Fig. 4 all have two major PV terraces (see also the diagonal cuts in and S8). The reference run in Fig. 4 has two impressively flat PV steps in Fig. 2. The strong PV mixing results from the chaotic vortex motion in Fig. 1 . Other runs in Fig. 4 and have two major PV terraces and one or two additional minor terraces.

The PV terraces in Fig. 4, with a lone polar cyclone, are not very flat; i.e., in the PV is imperfectly mixed. Imperfect mixing in Fig. 4 may result from the early establishment of approximate axisymmetry, and therefore small , in this particular run. In the other examples in Fig. 4 there is protracted radial mixing by the chaotic nonaxisymmetric motion of the vortices as they crystallize.

Crystal Radius and .

In Fig. 1 the central circle has radius and this circle passes rather closely through the seven circumpolar cyclones in Fig. 1 and I. Dozens of solutions, with different initial conditions and variations in , U, γ, and away from those of the reference solution in [], confirm that is invariably the radius of the central cyclonic crystal (additional examples are in Fig. 4 and ). Perhaps one might have anticipated that the radius of the crystal is —but numerical solutions show that the dimensionless constant c is close to one. Identification of as the radius of the vortex crystal is the most important result of this work.

Vortex Radius, , and Stability of a Multivortex Crystal.

In Fig. 1 the radii of the eight polar cyclones are less than the radius of the crystal, , and also greater than the initial excitation wavelength . Our simulations indicate that formation of a multicyclone crystal requires (in the reference run ). The final cyclone radii fall in the large gap between and .

Varying , with all other parameters fixed at the reference values, changes the size and number of cyclones in the final state. Because of sensitivity to initial conditions, the following discussion of is qualitative. Increasing results in fewer, larger cyclones in the final state (). We find vortex crystals with a central polar cyclone surrounded by five or more circumpolar cyclones for . We show results only for km because the three values produce qualitatively similar crystals. If 400 km km, we find either a single central cyclone or crystals consisting of two to four vortices. Six runs with km all resulted in a single central cyclone.

We speculate that to form a stable multicyclone crystal the constituent vortices must be separated from each other with an intervortex spacing that is significantly greater than a vortex radius; e.g., a rough criterion for vortex merger is that two like-signed vortices come within a critical separation distance equal to 3.3 times the average of the two radii (24–26). If is too large, then after a few vortex mergers the resultant big vortices cannot avoid merger and also fit into the circle. And with rather large —for example km in the bottom row of —only one or two cyclones form during the nucleation phase.

Results: The Flat Trap

In Fig. 5 we use the same initial condition and parameters as those in the polar-cap reference run in Fig. 1. But instead of , we employ the flat-trap in []. This configuration resembles the electron-plasma experiments, except that the initial vorticity has both signs. The rapid initial vortex nucleation in Fig. 5 is similar to that of the reference run in Fig. 1 and produces equal numbers of cyclones and anticyclones. Note, however, that in Fig. 5 the emergent vortex ensemble consists mainly of anticyclones moving through a uniform background of weakly positive vorticity.

Fig. 5.

(A–I) This flat-trap run uses the same parameters as the reference run in Fig. 1. But with the flat-trap planetary PV in in [] there is no vortex segregation. In contrast to Fig. 1the final state in G–I consists mainly of anticyclones distributed throughout the trap.

Because there is no planetary PV gradient, vortex segregation does not occur in the flat-trap configuration. In Fig. 5 the only process discriminating between cyclones and anticyclones is the anticyclonic retention and cyclonic expulsion (): The final crystal is anticyclonic because most of the cyclones have escaped from the trap. In Fig. 5 there is an ensemble of anticyclones distributed much more uniformly throughout the trap than in its polar-cap cyclonic cousin (Fig. 1). There are also a few small stray cyclones protectively paired with much larger anticyclones. Without vortex segregation, many more dipoles strike the boundary and so at long time there are many more cyclonic escapees in Fig. 5 than in Fig. 1.

The large ensemble in Fig. 5 is not as beautifully regular and crystalline as some of the electron-plasma examples (8, 9). Moreover, in Fig. 5 only the larger vortices in the ensemble are in approximate solid-body rotation. We speculate that this difference between the flat trap and the electron-plasma experiments may result from different initial conditions. The plasma initial condition is a spiral vorticity filament with a well-defined width: See figure 2 of Fine et al. (8). In our numerical solutions increasing results in fewer and larger vortices, and a more regular crystal, persisting in time and in quasi–solid-body rotation ().

The flat-trap results in Fig. 5 and highlight the importance of γ in forming polar vortex crystals. Even though the vortex-core PV extrema are much greater than the trap PV discontinuity in both configurations, the relatively weak polar-cap planetary PV gradient, , is essential to vortex segregation and to the formation of a polar crystal with radius .

Discussion

The polar-cap dynamical regime differs qualitatively from 1) plain and simple 2D turbulence, which, in doubly periodic geometry, leads eventually to two opposite-signed final vortices (22), and 2) β-plane turbulence, which produces zonal jets characteristic of Jovian low- and midlatitudes. In this respect, polar-cap dynamics are a different regime of geostrophic turbulence, with polar vortex crystals being the analog of β-plane zonal jets.

The variation of the Coriolis parameter in the vicinity of the poles (the γ-effect) is crucial, as can be seen from comparison of Fig. 1 with the flat-trap solution in Fig. 5: Even though the total variation of planetary PV resulting from γ is much smaller than the relative vorticity of individual cyclones, γ is nonetheless essential for the formation of a crystal such as those observed at the poles of Jupiter. γ is essential to the definition of the intrinsic length scale , which emerges as the radius of the vortex crystals in the polar-cap configuration. Crystals form without γ, but these anticyclonic flat-trap crystals occupy the entire domain.

In Earth’s oceans a length scale analogous to may be important for the dynamics of long-lived anticyclones contained within roughly axisymmetric topographic bowls (37). The topographic PV gradient of the bowl greatly exceeds β. In these ocean situations, however, only solitary trapped anticyclones (the Loften eddy, the Mann eddy, etc.) have been observed. Another example is the suggestion of Bouchet and Sommeria (38) that the emergence of propagating oval-shaped midlatitude vortices, such as on Jupiter, requires compensation of the β-effect by vortex motion relative to the background flow. Because the β-effect is compensated, remaining quadratic terms in the environmental PV, analogous to γ, are important in determining the size of these midlatitude ovals.

Although γ is necessary for polar crystal formation, γ is not sufficient. Several authors have studied γ-driven polar vorticity accumulation and documented the formation of a lone polar cyclone (17–19, 29); these studies did not find multivortex crystals. We find that multivortex crystals most readily form with small values of . If is too large, then we also find only lone polar cyclones. We speculate that earlier studies may have involved cyclones that were too large relative to to form crystals.

For simplicity, and by analogy with the electron-plasma case, we have used a barotropic model. A finite deformation radius complicates our model because the conversion between kinetic and potential energies resulting from finite implies that the initial root-mean-square velocity, U defined in [], is not constant during the evolution. We take comfort in recent observations and numerical results (3, 14, 39, 40) indicating that polar cyclones might have a strong barotropic component.

Recent discussions (12, 13) of the stability and structure of polar cyclones model the crystal using shielded vortices; e.g., azimuthal velocity decreases exponentially with distance from the vortex center (much faster than ). These authors emphasize the role of a shielded central polar cyclone and the planetary PV, , in forming a circumpolar PV ridge, i.e., a local radial maximum in PV. Cyclones move up the PV gradient and thus this local maximum of q is a stable radial location at which circumpolar cyclones might sit (12, 13).

The scenario above does not apply to the numerical solutions in this paper. Within the polar cap, PV homogenization eliminates the poleward drift of cyclones. The vortex crystal, consisting of unshielded vortices, sits on top of a flat PV terrace.

Thus, although shielded vortices and the presence of a shielded central polar cyclone may help form a stable crystal (12, 13), it is clear from the minimal model presented here that these two ingredients are by no means necessary: The unshielded vortices in Fig. 3 and C indicate that shielding is not necessary to crystal formation and stability. No shielding is also consistent with Jovian observations (14) and with electron-plasma experiments. Moreover, in a few of our numerical solutions emergent vortex crystals are “hollow”; i.e., there is no central cyclone. There are five examples of hollow crystals in and three examples in figure 1 of Schecter et al. (9). These results are consistent with Kelvin’s point-vortex crystal, which consists of unshielded point vortices that form stable hollow crystals, with as many as seven vortices arranged in a ring (7).

Several differences between our crystals and Jovian crystals need to be mentioned. Jovian crystals have a preferred rotational frame of reference: They are quasi-stationary with respect to the System III reference frame (2), whereas the rotational frame of our crystals varies between realizations; i.e., the five multivortex crystals in Fig. 4 are in solid-body rotation with different rotation rates in the different realizations. Our simple model purposely omits processes that probably contribute to the stability and persistence of Jovian crystals and to the selection of a preferred rotational reference frame. Such processes may include the deep convective origin of the cyclones (39) and dissipation via deep magnetohydrodynamic coupling to the Jovian magnetic field (41).

Jovian vortices are more closely packed than those emerging in our simulations. This may reflect the absence of forcing, as suggested by the forced-dissipative nonhydrostatic polar-cap model of Cai et al. (40), which produces deep, closely packed and long-lived polar cyclones. A future avenue of research is examination of crystal formation and rotation and the vortex spacing in forced-dissipative configurations. It is hoped that sensitivity to initial conditions is less important in the forced-dissipative case.

Our integrations lasted only 300 crystal rotation periods with no indication of significant changes in the crystal on that time scale. Is this final state actually final, or are there spontaneous transitions to other states in the unexamined future? For example, one experimental electron-plasma crystal endured for about 104 rotation periods, after which one vortex disappeared and a new crystal reformed with one less vortex (11). Barotropic β-plane turbulence also presents examples of spontaneous transitions out of long-lived metastable states (42); e.g., on extremely long time scales, turbulent solutions abruptly switch to new configurations with a different number of zonal jets. It is possible that 300 crystal rotation periods are not long enough to reveal analogous rare transitions in the structure of a polar vortex crystal.