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 uid and the polar parameter γ = f p /a 2 p , with f p the Coriolis parameter at the pole and a p the polar radius of the planet.

The Juno spacecraft revealed a distinctive form of organized turbulence in the polar atmosphere of Jupiter (1). e North Pole exhibits a central polar cyclone surrounded by eight circumpolar cyclones; the South Pole is characterized by ve cyclones surrounding a central polar cyclone. ese cyclones are organized into a vortex crystal: a symmetric array of strong vortices in solid-body rotation within a background of weaker vorticity.

e 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). e 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)(6)(7). Laboratory experiments with magnetized pure-electron plasmas are another antecedent (8)(9)(10)(11). is experimental system is isomorphic to a near-ideal two-dimensional (2D) uid contained within a circular domain, with electron density equivalent to vorticity. rough vortex nucleation, merger, and nally self-organization, regular vortex arrays, in solid-body rotation (9), spontaneously crystallize from a lamented initial vorticity eld.

ese electron-plasma experiments are simpler than polar planetary dynamics in several respects: 1) ere is no analog of the variation of the Coriolis parameter with latitude (the β-e ect); 2) e vorticity has strictly one sign (electrons are negatively charged); 3) e plasma is contained within a circular domain with a free-slip boundary; 4) Plasma dynamics are "barotropic" (the deformation length is in nite); and 5) ere is e ectively 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 electronplasma vortex crystals and Jovian observations (1,4,12,13). While these two systems have many features in common, it is crucial to relax simpli cations 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). e barotropic component of the ow (i.e., the circumpolar and polar cyclones) exhibits kinetic energy spectra with a k -3 spectral slope, consistent with quasigeostrophic (QG) dynamics.

ese 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).

Vortex crystals, geometric arrays of like-signed vortices, are observed in natural systems with vastly different space and time scales: at the poles of Jupiter (∼10,000-km radius and lifetime of at least 5 y) and in laboratory experiments with pure-electron plasma (∼3.5-cm radius, lifetime of about 1.7 s). We follow the adage "less is more" and show that minimal physics is required for polar vortex crystals formation and persistence. Crystals, resembling those of Jupiter, form from the free evolution of an unstratified and rapidly rotating fluid in an axisymmetric geometry. An essential ingredient in this minimal model is the decrease of the vertical component of the Coriolis force with distance from the pole. Once formed, the crystal seems to survive indefinitely.

An important characteristic of the QG polar-cap regime is that planetary parameters occur only in the combination γ = f p /afoot_0 p , where f p is the Coriolis parameter at the pole and a p is the polar radius of the planet. e parameter γ is signi cant because close to the poles the Coriolis parameter is f (r ) ≈ f p -1 2 γr 2 , where r is the distance from the pole. We discuss the North Pole so that f p is positive and cyclones are positive vortices.

e barotropic polar-cap QG system is

where ψ(x , y, t) is the stream function, ζ = ψ xx + ψ yy is the relative vorticity, and q is the potential vorticity (PV). We employ a Cartesian coordinate system centered on the pole so that r = (x 2 + y 2 ) 1/2 . For the planetary PV η in [2], it is instructive to compare two con gurations. e rst one is the polar-cap PV

where 1 (r ) is the indicator function of the disk with radius r ; i.e., 1 (0 < r < r ) = 1 and 1 (r < r ) = 0. Inside the specied radius r , η pc (r ) is equal to the planetary PV, -1 2 γr 2 . e second con guration is the " at-trap" PV

e at-trap model η ft (r ), with no planetary PV gradient within the trap, is close to the electron-plasma experiments. Both η pc (r ) and η ft (r ) have a discontinuity of strength q = 1 2 γr 2 at the trap radius r .

is discontinuity results in an enclosure that largely con nes the QG turbulence within the disk 0 < r < r (SI Appendix, section 2). e model is solved as an unforced initial value problem using a doubly periodic pseudospectral method provided by Geophys-icalFlows.jl (16); the domain size is L dom × L dom . We use a trap radius r that is as large as possible, while still ensuring that there is minimal interaction between neighboring periodic domains.

e vortex trap is an alternative to the dissipative sponge layer used in earlier studies of polar-cap dynamics (17)(18)(19). e trap approximately conserves energy throughout the 8-y duration of the simulation. For further numerical considerations, and discussion of the r trap boundary, see SI Appendix, section 2.

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

e initial relative vorticity ζ 0 (x , y) is a random monoscale eld, multiplied by a taper that con nes ζ 0 to the center of the trap: See Fig. 1A and SI Appendix, section 1 for further details of the initial condition. e untapered initial relative vorticity is concentrated within a narrow annulus in wavenumber space; the central radius of the annulus is 2π/L init and the width of the annulus is much less than 2π/L init . In addition to L init , the initial condition provides a velocity scale U de ned by

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

The Intrinsic Polar-Cap Length Scale. From U and γ one can form the length

L γ 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 β(r ) = γr and therefore the Rhines length at a distance r from the pole is

Well-known arguments identifying L R in midlatitudes all assume that the curvature of the planetary PV is negligible over the scale L R itself, i.e., that one can approximate the planetary PV by a linear function of latitude over a distance L R . Eq. 7 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 L R as r → 0: Equating r to L R in [7] gives L R = (U /γL R ) 1/2 . Solving for L R then results in the length L γ on the left of [6]. us L γ is identi ed as an intrinsic polar length scale. We show below that L γ is the radius of emergent QG polar vortex crystals.

Disclaimer. Despite the adoption of Jupiter-like numbers in [8] below, we do not claim quantitative Jovian accuracy for this barotropic QG model. For example, we use barotropic dynamics so that the deformation length L def is in nite (as it is in the plasma experiments). Observations summarized in ref.

12 indicate that the L def at the Jovian poles is somewhere in the range 350 to 1,300 km. is range encompasses various polar length scales, such as the vortex diameter de ned as the radius of maximum azimuthal velocity (∼1,000 to 1,200 km). us baroclinic e ects associated with nite L def 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 su cient to capture the key vortex-dynamical processes responsible for the genesis of polar vortex crystals.

Fig. 1 shows a solution of [1] and [2] with η pc (r ) in [3]. is "reference run" uses Jupiter-like parameters

[8] L init 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 SI Appendix, section 4 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 r , so that the vortex crystal with radius L γ (see below) is well separated from complications at the periphery of the trap. With available computational resources we use r = 5L γ (see SI Appendix, section 7 for a discussion of sensitivity to L γ /r ). e PV discontinuity at r = 5L γ is q = 1 2 γr 2 = 9.9451 × 10 -5 s -1 . e domain length scale is L dom = 12L γ . e initial condition in Fig. 1A shows random small-scale (L init L γ ) vorticity concentrated in the center of the trap, with PV extrema q ext ∼ 6 × 10 -3 s -1 ∼ 60 q . e development of a vortex crystal from random initial conditions is shown in detail in Movies S1-S3. e process is summarized by nine snapshots of the relative vorticity ζ in Fig. 1. ere are several distinct evolutionary stages corresponding to the three rows of Fig. 1.

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

Vortex Segregation and Merger. Subsequent evolution from the state in Fig. 1C involves merger (24)(25)(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)(28)(29). Segregation is noticeable in Fig. 1E and is largely complete after 75 d in Fig. 1F : e pole is occupied by cyclones and anticyclones have accumulated up against the trap boundary at r .

In view of PV mixing arguments, invoked below in connection with [9], it is interesting that vortex segregation is "antidi usive": Cyclones (ζ > 0) move poleward and reinforce the initial r = 0 maximum in the planetary PV -1 2 γr 2 . And the outward radial motion of anticyclones (ζ < 0) reinforces the initial minimum in η pc (r ) at r .

Emergence of the Crystal. After segregation is complete there are eight large cyclonic vortices moving chaotically within a cap of radius L γ centered on the pole; the snapshot in Fig. 1G is representative of this state. e con guration in Fig. 1G 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). is jostling happens repeatedly until the con guration in Fig. 1H is nally achieved after about 1.5 y of evolution. e central cyclone in Fig. 1H remains in place until the end of the run in Fig. 1I at 8 y.

Once the permanent central cyclone is established, the octet crystallizes into solid-body rotation (Movie S3). e 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). e oscillations may result from di erences 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 ow, one should not ignore the zonal-mean ow ū(r , t) and the associated zonal mean relative vorticity ζ = (r ū) r /r . (See SI Appendix, section 5 for a summary and discussion of zonal mean dynamics.) e eight cyclones interact with one another by mutual advection, but they also interact with the zonal-mean ow ū. For example, and in analogy with the electron-plasma experiments, we argue that the chaotic state in Fig. 1G is "cooled" as the irregular motion of the polar cyclones mixes the zonal-mean PV, q = ζ -1 2 γr 2 . e 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, ūt , is related to radial vorticity transport, v ζ , by ūt + v ζ = 0.

[9]

In the central L γ circle chaotic vortex motion directly results in v ζ and mean-ow acceleration. But there is also far-eld v ζ and, through [9], far-eld zonal-mean acceleration. ( e "far eld" of the circumpolar cyclones is de ned by r L γ .) Chaotic vortex motion within the L γ circle produces a uctuating far-eld radial velocity, v ∼ r -2 , corresponding to the unsteady dipole moment of the circumpolar cyclones. ere is also an azimuthal component, u ∼ r -1 , resulting from the constant-in-time total circulation. But only the weaker radial component, v , results in radial transport v ζ and therefore far-eld mean-ow acceleration via [9].

Mixing is so strong that by the end of the reference run there is striking PV homogenization (32). Fig. 2A shows that the PV is mixed into discrete steps (33,34). us, 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: e planetary γ-e ect is annulled by the relative vorticity of the mean azimuthal ow ū. And there is no further mean-ow acceleration because ζ = 0; i.e., PV uctuations cannot be created by stirring homogeneous PV.

Because of PV homogenization the eight cyclones in Fig. 1 G-I sit on top of a at PV terrace. roughout the rst PV terrace in Fig. 2A, q ≈ -q /4 is less than the planetary vorticity, -1 2 γr 2 . us, 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. 2B shows a comparison of the azimuthally averaged azimuthal velocity, ū(r ), of the reference run with an analytic t based on homogenized PV steps and Gaussian vorticity pro les for the polar cyclones (see SI Appendix, section 6 for details of this PV terrace model). Although ζ(r ) < 0 throughout the rst terrace, the azimuthal velocity ū(r ) in Fig. 2B is dominated by the irrotational r -1 far-eld velocity of the central cyclone and is therefore strongly positive. e secondary local maximum in ū(r ) at r /r ≈ 0.25 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. 3B the vorticity has one sign (positive) and each cyclone has a nonzero circulation. erefore, the far-eld azimuthal velocity around each polar cyclone is u = κ/2πr v . (Here κ > 0 is the circulation of the cyclone; r v is the distance from the vortex center; and "far-eld" means that r v is much greater than the vortex radius, de ned as the radius of maximum velocity in Fig. 3C.) Slow algebraic decay, u ∝ r -1 v , di ers qualitatively from the exponential decay of u assumed in recent models: ose models employ completely shielded vortices with κ = 0 (12, 13).

In Fig. 3C the maximum cyclone velocity is between 5 and 10 × U, where U is the root-mean-square velocity de ned in [5] and set to 80 m • s -1 . e velocity 80 m • s -1 is more typical of the maximum cyclone velocity of Jovian cyclones-see gure 6 of ref. 36. erefore, the cyclones in Fig. 3C have peak velocities very much faster than those of Jovian polar cyclones. All velocities in the polar-cap system are, however, proportional to U (SI Appendix, section 4). us we can make an a posteriori adjustment so that the maximum cyclone velocity in the polarcap model matches the Jovian observations: Reduce U by a factor of 8, which reduces L γ 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 L γ = 5,000 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. e average cyclone radius indicated in Fig. 3C curve is ū(r) calculated from the final state (8 y) of the reference run; the "theoretical" curve is ū(r) 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 r > r (see SI Appendix, section 6 for details).

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 L γ = 10 4 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 in nite deformation radius. Sensitivity to Initial Conditions. Although the emergence of L γ 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 [8] xed, if we change the seed in the random number generator responsible for initializing ζ 0 , we produce a diverse set of nal con gurations; see Fig. 4 for six examples. In Fig. 4E there is a dipole consisting of a large polar cyclone and a smaller cyclone; this dipole orbits the pole. In Fig. 4C, 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 gure 1 of Adriani et al. (1). In Fig. 4F there is a lone polar cyclone. Electron-plasma experiments exhibit an analogous sensitivity to initial conditions (8,9). e structure of the zonal-mean PV is less sensitive to initial conditions: e six runs in Fig. 4 all have two major PV terraces (see also the diagonal cuts in SI Appendix, Figs. S5 andS8). e reference run in Fig. 4A has two impressively at PV steps in Fig. 2A. e strong PV mixing results from the chaotic vortex motion in Fig. 1 D-F. Other runs in Fig. 4 and SI Appendix, Fig. S3 have two major PV terraces and one or two additional minor terraces.

e PV terraces in Fig. 4F, with a lone polar cyclone, are not very at; i.e., in SI Appendix, Fig. S5 the PV is imperfectly mixed. Imperfect mixing in Fig. 4F may result from the early establishment of approximate axisymmetry, and therefore small v ζ , 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 Lγ . In Fig. 1 the central circle has radius L γ and this circle passes rather closely through the seven circumpolar cyclones in Fig. 1 H and I. Dozens of solutions, with di erent initial conditions and variations in r , U, γ, and L init away from those of the reference solution in [8], con rm that L γ is invariably the radius of the central cyclonic crystal (additional examples are in Fig. 4 and SI Appendix, Fig. S3). Perhaps one might have anticipated that the radius of the crystal is c L γ -but numerical solutions show that the dimensionless constant c is close to one. Identi cation of L γ as the radius of the vortex crystal is the most important result of this work.

Vortex Radius, L init , and Stability of a Multivortex Crystal. In Fig. 1 G-I the radii of the eight polar cyclones are less than the radius of the crystal, L γ , and also greater than the initial excitation wavelength L init . Our simulations indicate that formation of a multicyclone crystal requires L γ L init (in the reference run L γ /L init = 50). e nal cyclone radii fall in the large gap between L init and L γ .

Varying L init , with all other parameters xed at the reference values, changes the size and number of cyclones in the nal state. Because of sensitivity to initial conditions, the following discussion of L init is qualitative. Increasing L init results in fewer, larger cyclones in the nal state (SI Appendix, Fig. S3). We nd vortex crystals with a central polar cyclone surrounded by ve or more circumpolar cyclones for L init ∈ {100 km, 200 km, 300 km}. We show results only for L init = 200 km because the three values produce qualitatively similar crystals. If 400 km ≤ L init ≤ 1,200 km, we nd either a single central cyclone or crystals consisting of two to four vortices. Six runs with L init = 2,000 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 signi cantly greater than a vortex radius; e.g., a rough criterion for vortex merger is that two likesigned vortices come within a critical separation distance equal to 3.3 times the average of the two radii (24)(25)(26). If L init is too large, then after a few vortex mergers the resultant big vortices cannot avoid merger and also t into the L γ circle. And with rather large L init -for example L init = 2,000 km in the bottom row of SI Appendix, Fig. S3-only one or two cyclones form during the nucleation phase.

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 η pc (r ), we employ the at-trap η ft (r ) in [4]. is con guration resembles the electron-plasma experiments, except that the initial vorticity has both signs. e rapid initial vortex nucleation in Fig. 5 A-C 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 G-I the emergent vortex ensemble consists mainly of anticyclones moving through a uniform background of weakly positive vorticity.

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

e large ensemble in Fig. 5 G-I 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 di erence between the at trap and the electron-plasma experiments may result from di erent initial conditions. e plasma initial condition is a spiral vorticity lament with a well-de ned width: See gure 2 of Fine et al. (8). In our numerical solutions increasing L init results in fewer and larger vortices, and a more regular crystal, persisting in time and in quasi-solid-body rotation (SI Appendix, Fig. S4).

e at-trap results in Fig. 5 and SI Appendix, Fig. S4 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 con gurations, the relatively weak polar-cap planetary PV gradient, γ r , is essential to vortex segregation and to the formation of a polar crystal with radius L γ .

e polar-cap dynamical regime di ers qualitatively from 1) plain and simple 2D turbulence, which, in doubly periodic geometry, leads eventually to two opposite-signed nal vortices (22), and 2) β-plane turbulence, which produces zonal jets characteristic of Jovian low-and midlatitudes. In this respect, polar-cap dynamics are a di erent regime of geostrophic turbulence, with polar vortex crystals being the analog of β-plane zonal jets. e variation of the Coriolis parameter in the vicinity of the poles (the γ-e ect) is crucial, as can be seen from comparison of Fig. 1 with the at-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 poles of Jupiter. γ is essential to the de nition of the intrinsic length scale L γ , which emerges as the radius of the vortex crystals in the polar-cap con guration. Crystals form without γ, but these anticyclonic at-trap crystals occupy the entire domain.

In Earth's oceans a length scale analogous to L γ may be important for the dynamics of long-lived anticyclones contained within roughly axisymmetric topographic bowls (37).

e 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 β-e ect by vortex motion relative to the background ow. Because the β-e ect 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 su cient. Several authors have studied γ-driven polar vorticity accumulation and documented the formation of a lone polar cyclone (17)(18)(19)29); these studies did not nd multivortex crystals. We nd that multivortex crystals most readily form with small values of L init /L γ . If L init /L γ is too large, then we also nd only lone polar cyclones. We speculate that earlier studies may have involved cyclones that were too large relative to L γ to form crystals.

For simplicity, and by analogy with the electron-plasma case, we have used a barotropic model. A nite deformation radius L def complicates our model because the conversion between kinetic and potential energies resulting from nite L def implies that the initial root-mean-square velocity, U de ned in [5], 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 r -1 ). ese authors emphasize the role of a shielded central polar cyclone and the planetary PV, -1 2 γr 2 , 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). e scenario above does not apply to the numerical solutions in this paper. Within the polar cap, PV homogenization eliminates the poleward drift of cyclones. e vortex crystal, consisting of unshielded vortices, sits on top of a at PV terrace. us, 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: e unshielded vortices in Fig. 3 B 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. ere are ve examples of hollow crystals in SI Appendix, Fig. S3A and three examples in gure 1 of Schecter et al. (9). ese 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 di erences between our crystals and Jovian crystals need to be mentioned. Jovian crystals have a preferred rotational frame of reference: ey 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 ve multivortex crystals in Fig. 4 A-E are in solid-body rotation with di erent rotation rates in the di erent 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 eld (41). Jovian vortices are more closely packed than those emerging in our simulations. is may re ect 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 longlived polar cyclones. A future avenue of research is examination of crystal formation and rotation and the vortex spacing in forceddissipative con gurations. 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 signi cant changes in the crystal on that time scale. Is this nal state actually nal, or are there spontaneous transitions to other states in the unexamined future? For example, one experimental electron-plasma crystal endured for about 10 4 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 con gurations with a di erent 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. Data Availability. There are no data underlying this work. The code used to generate the runs have been deposited in GitHub, https://github.com/ FourierFlows/GeophysicalFlows.jl.