Vorticity and divergence at scales down to 200 km within and around the polar cyclones of Jupiter

Since 2017 the Juno spacecraft has observed a cyclone at the north pole of Jupiter surrounded by eight smaller cyclones arranged in a polygonal pattern. It is not clear why this configuration is so stable or how it is maintained. Here we use a time series of images obtained by the JIRAM mapping spectrometer on Juno to track the winds and measure the vorticity and horizontal divergence within and around the polar cyclone and two of the circumpolar ones. We find an anticyclonic ring between the polar cyclone and the surrounding cyclones, supporting the theory that such shielding is needed for the stability of the polygonal pattern. However, even at the smallest spatial scale (180 km) we do not find the expected signature of convection—a spatial correlation between divergence and anticyclonic vorticity—in contrast with a previous study using additional assumptions about the dynamics, which shows the correlation at scales from 20 to 200 km. We suggest that a smaller size, relative to atmospheric thickness, of Jupiter’s convective storms compared with Earth’s, can reconcile the two studies. The extraction of the wind pattern, vorticity and divergence down to a scale of 200 km from the cluster of cyclones at Jupiter’s north pole shows evidence of an anticyclonic ring, which is needed to keep the system stable, around the central cyclone. No signatures of convection are observed at 200 km scales.

A t Jupiter's north pole there are eight cyclones that form an octagon, with one cyclone at each vertex and one additional cyclone in the centre 1,2 . The centres of the cyclones are at latitudes of 83 ± 10°, which is about 8,700 km from the pole. Jupiter's south pole is the same, except there are only five cyclones, which form a pentagon with one at the centre. The vertices are at latitudes of −83 ± 1°. The polygons and the individual vortices that comprise them have been stable for the 4 years since Juno discovered them [3][4][5] . The polygonal patterns rotate slowly, or not at all. The peak azimuthal wind speeds around each vortex range from 70 to 100 m s −1 , and the radial distance r from the peak to the vortex centre is about 1,000 km (ref. 6 ). In contrast, Saturn has only one vortex, a cyclone, at each pole 7 . The peak winds are 150 m s −1 , and the radius at the peak is 1,500 km (refs. 8,9 ). Saturn has a six-lobed meandering jet at 75°, but it has no cyclones associated with it. Both laboratory and theoretical models treat the hexagon as a stable wave-like pattern [10][11][12][13][14] .
There have been a handful of theoretical studies that specifically address the origin of polar cyclones on Jupiter and Saturn [15][16][17][18] . They comprise one-and two-layer models that introduce small-scale motions either as an initial condition or as continuous forcing balanced by dissipation. The small-scale vortices merge and become the large-scale vortices. The cyclones drift polewards, and the anticyclones drift equatorwards, as they do on Earth. In some cases the cyclones merge into one big cyclone at the pole. In other cases, with different parameter settings, the cyclones wander about without forming polygons. Only one theoretical study obtained stable polygons from random initial conditions, and only when the wavelengths of the initial random disturbances are less than 300 km (ref. 19 ). A Fourier analysis of Juno data reveals that flows with wavelengths larger than 215 km are gaining energy from smaller-scale flows-an example of an upscale energy transfer 20 . Therefore, one goal of this Article is to measure vorticity and divergence at scales much smaller than the main cyclones and determine how the upscale energy transfer takes place.
Another theoretical study 21 , which used shallow water equations, introduced cyclones that have the observed gross propertiesmaximum velocity and radius-and arranged them into different polygonal patterns around the pole to see which ones are stable. The stable ones have shielding (a ring of anticyclonic vorticity surrounding each of the cyclones) and the unstable ones do not. Some models with small-scale forcing develop shielding, but they do not organize into polygons [15][16][17][18]22 . So another goal of this Article is to measure the vorticity inside and outside the large cyclones and see whether they are shielded.
The small-scale forcing in the one-and two-layer models is a crude representation of convection. There are 3D models that treat convection explicitly, in some cases with the Boussinesq (quasi-incompressible) approximation [23][24][25] and in other cases with density varying vertically by up to five scale heights [25][26][27][28] . Some treat fluid in a box with periodic boundary conditions, and others use full spherical geometry. All the 3D models have small-scale convective plumes. The convective plumes produce large-scale vortices by mergers, an upscale transfer of kinetic energy, and some of the vortices arrange themselves into polygonal patterns 23,28 . Although a relation between divergence and vorticity is not discussed in any of these models, a negative correlation is expected for convection on a rotating planet. Therefore, a third goal of this Article is to measure divergence and vorticity at scales down to 180 km and search for this signature of convection.
Since 2017 the Juno spacecraft has observed a cyclone at the north pole of Jupiter surrounded by eight smaller cyclones arranged in a polygonal pattern. It is not clear why this configuration is so stable or how it is maintained. Here we use a time series of images obtained by the JIRAM mapping spectrometer on Juno to track the winds and measure the vorticity and horizontal divergence within and around the polar cyclone and two of the circumpolar ones. We find an anticyclonic ring between the polar cyclone and the surrounding cyclones, supporting the theory that such shielding is needed for the stability of the polygonal pattern. However, even at the smallest spatial scale (180 km) we do not find the expected signature of convection-a spatial correlation between divergence and anticyclonic vorticity-in contrast with a previous study using additional assumptions about the dynamics, which shows the correlation at scales from 20 to 200 km. We suggest that a smaller size, relative to atmospheric thickness, of Jupiter's convective storms compared with Earth's, can reconcile the two studies.
velocity and β = df/dy = 2Ωsinθ/a is the latitudinal gradient of the planetary vorticity f = 2Ωcosθ. Here θ is colatitude, y is the northward coordinate, and a and Ω are the planet's radius and angular velocity, respectively. L β plays a role in the stability of the zonal jets. On both Jupiter and Saturn, 2π/L β is approximately equal to the wavenumber of the zonal jet profile with respect to latitude when U is the root mean squared speed [29][30][31][32][33] . However, L β is infinite at the poles as β linearly approaches zero there. We therefore introduce a different scaling 19 , one based on the inverse gradient of β at the pole, γ = −dβ/dy = 2Ω/a 2 . The associated length scale is L γ = (U/γ) 1/3 , and for U = 80 m s −1 it is about 10,500 km (we ignore the oblateness and use Jupiter's equatorial radius throughout this paper). L γ represents the radius of the circle around the pole inside which the effect of the vortices-large-scale turbulence-is greater than the effect of β and the zonal jets. Note that L γ is the distance from a specific point (the pole) and L β is not. The value of L γ is close to the 8,700 km size of the polygons on Jupiter.
The radius of deformation L d is c/f, where c is the gravity wave speed of the gravest vertical mode-the one spanning Jupiter's weather layer, which extends from the base of the stratosphere down to the base of the water cloud. The value of c depends on the degree of stratification of the weather layer 34 , and is assumed to be independent of latitude. Different assumptions about the vertical stratification put the average L d at the poles in the range 350-1,300 km (refs. [34][35][36] ). This brackets the 1,000 km radius of the cyclones.
Originally L β was defined as the length scale where the flow transitions from turbulence to zonal jets as the scale of the flow increases 37 . Observations of Jupiter suggest that a similar transition occurs as the latitude of the flow decreases. The critical latitude, below which zonal jets dominate, was shown to be a decreasing function of L d (refs. 38,39 ). Here we argue that a critical latitude, π/2 − L γ /a, exists even for arbitrarily small values of L d .
We discuss the observations using parameters of the shallow water equations, where a single layer of fluid of thickness h floats hydrostatically on a much thicker fluid, which we assume is at rest 40,41 . The two dependent variables are the horizontal velocity v and the gravitational potential ϕ = g r h, where g r is the reduced gravity (the gravitational acceleration times the fractional density difference Δρ/ρ between the two layers 41 ). Here ϕ is the column density in the 2D continuity equation and √ϕ is the gravity wave speed c; ϕ controls vortex stretching and enters in the expression for potential vorticity (PV), which is a dynamical scalar that is conserved in fluid elements. For the shallow water equations PV is (ζ + f)/ϕ, where ζ = (∇ × v) · k is the relative vorticity (the curl of the horizontal velocity). Three-dimensional effects are not completely ignored: they enter through L d , which is proportional to the square root of h and the fractional density difference Δρ/ρ. Even these quantities are uncertain, so given the paucity of information about vertical structure, it is best to discuss our observations with the shallow water equations.
Vorticity and divergence. Figure 1 shows the octagon of cyclones surrounding the north pole 2 . Features in the clouds are visible at scales down to ~100 km, which is much smaller than the 1,000 km radius where the azimuthal velocity is greatest.  Figure 2 shows vorticity and divergence maps for two independent determinations of the wind field. The measurement required clouds in a pair of JIRAM images separated in time to be tracked to obtain the velocity, and then closed line integrals were taken to get vorticity and divergence. The magnitude of the vorticity is larger than that of the divergence. The persistence and movement of vorticity features, even those ~180 km in size, shows that the small-scale features are not measurement noise. The motion is visible when one toggles between the left and right vorticity maps in Fig Planetary signal and measurement noise. In Fig. 3, the top two panels show covariances between n0103 and n0204, which are the two independent determinations of the wind fields in Fig. 2. Vorticity is top left, and divergence is top right. If the two measurements gave exactly the same values at each point, the arrays of points would lie on a straight line running from the lower left corner to the upper right corner. Deviations from this line, that is, a correlation coefficient less than 1.0, come partly from measurement noise and partly from cloud motions on the planet. The vorticity measurements have a correlation coefficient η of 0.729. As the measurement noise is uncorrelated with cloud motions on the planet, the noise and planetary variances add up (equation (1)). And given that the noise from n0103 is uncorrelated with that from n0204, the covariations do not contain the variance of measurement noise (equation (2)):

Fig. 1 | Infrared image of the northern hemisphere as seen by JIRAM.
The circle at 80° latitude is about 12,000 km from the pole. The lines of constant longitude are 15° apart. The radiances have been corrected for nadir viewing, with bright yellow signifying greater radiance and dark red signifying lesser radiance. The average brightness temperature is in the range 215-220 K. Fig. 2 covers the central cyclone and the two cyclones at 135° and 315° east longitude, respectively. The two dark features at 120-150° E, 86° N whose filaments spiral towards their centres in a clockwise direction are anticyclones. Figure reproduced with permission from ref. 2 , Springer Nature Limited.
The quantities involving x and y are measurements and are known, so the measured correlation coefficient η allows one to separately determine the variance σ 2 p of vorticity on the planet and the variance σ 2 n of measurement noise. The same reasoning applies to the divergence (top right panel in Fig. 3). Table 1 shows the results for different values of the dimensions of the box used to measure vorticity and divergence. For the upper right panel of Fig. 3 and the 180 km × 180 km box, η = 0.299018, meaning that there is some divergence on the planet, but its variance is less than the measurement noise. Note that the noise values for vorticity and divergence are about the same for the same box size. The difference is that divergence decreases by a factor of 4.34, and vorticity decreases only by a factor of 1.61 from the 90 km × 90 km box to the 360 km × 360 km box. This difference is an indication that the divergence is a small-scale phenomenon that averages out for the larger box sizes. For the lower two panels of Fig. 3, divergence is plotted on the y axis with vorticity on the x axis, and η is essentially zero. Supplementary Table 2 shows that the noise estimates in Table  1 are best fitted by assuming that the measurement uncertainty for each component of velocity is 7.8 m s −1 .
Potential vorticity and shielding. Figure 4 shows the azimuthal mean v of the azimuthal velocity around the central cyclone as a function of radius r out to 6,000 km. Also shown are the mean relative vorticity ζ , the mean gravitational potential φ and the mean potential vorticity PV. In each panel there are three smooth curves. The middle one, coloured orange, was derived by a linear least squares fit to the velocity data. The basis functions are given in the Methods. The profile of v(r) agrees with earlier estimates 6 , including the profile at r > 2,000 km falling off faster than 1/r, implying negative vorticity in that region. Note that the fitted curve fits the data even where the velocity becomes negative (clockwise) at 4,000-6,000 km radial distance. The tabulated data are in Supplementary Data 3.
Bottom left: the ϕ values in Fig. 4 were computed from an integral and are therefore uncertain by an additive constant. However, ϕ = g r h is proportional to the thickness, and the thickness cannot be negative. The local maximum of ϕ is at r = 4,075 km, and Fig. 4 was computed with ϕ = 124 × 10 3 m 2 s 2 there. Having ϕ > 0 at the origin requires ϕ > 69 × 10 3 m 2 s 2 at r = 4,075 km. This gave L d > √φ /f = 749 km at r = 4,075 km, which is in the middle of the estimates obtained from lower latitudes when the variation in f with latitude was taken into account [34][35][36] .

Discussion
In a shallow water model that starts with cyclones of the observed size and velocity arrayed in polygonal patterns around the Jovian pole, stability requires an anticyclonic ring (shielding) around each cyclone 21 . With a peak azimuthal velocity of 80 m s −1 and radius at the peak of 1,000 km, a single parameter (b) controls the shape of the velocity profile and the depth of the shielding. The other free parameter is L d , but it has only a small effect on the results. The polygons are mainly stable in the range 1 < b < 2. Below this range, the shielding is too weak and the vortices merge. Above this range, the negative vorticity is too strong and the anticyclonic rings become two satellites orbiting around the cyclone 180° apart. At b > 3, these tripoles are unstable and the polygons fly apart chaotically. The blue curve in Fig. 4 has b = 1.35, which is safely in the stable zone The 200 km scale of vorticity and divergence is at least consistent with convection, Severe thunderstorms on Earth have diameters of 30-40 km, which is about five times the pressure scale height 43,44 . Granules, which are the convective elements in the solar photosphere, are about 1,000 km in diameter 45 , which is also about five times the scale height of the partially ionized hydrogen gas. The scale height on Jupiter is about 40 km at the water cloud base, so if the ratio of horizontal diameter to scale height were 5, as it is on Earth and the Sun, then convection elements on Jupiter would have diameters of 200 km. This is about the smallest scale we can measure.
The bottom row of Fig. 3 shows no correlation between divergence and vorticity, although both positive and negative values are At 80 m s −1 , which is about the maximum speed of the clouds, a feature moves 38 km. This is considerably less than the 180 km box size used to measure vorticity and divergence. Thus, to a good approximation, the two determinations show the same cloud features at the same time. The small motion is still visible, however, when it takes place on a large scale and includes many small-scale, high-contrast features, as with the slight rotation of the features visible in Fig. 2 and Supplementary Figs 1-4. bottom: points are plotted with divergence on the vertical axis and vorticity on the horizontal axis for n0103 (left) and n0204 (right). present. The implications for convection are uncertain. On the one hand, if a parcel conserves PV around a cycle of updrafts and downdrafts, the material derivative of ζ + f is proportional to the material derivative of ϕ throughout the cycle. But as ϕ is proportional to h, the material derivative of ϕ is also proportional to minus the divergence ∇ · v. Equivalently, the material derivative of negative (that is, anticyclonic) ζ is proportional to the divergence. As a result, negative vorticity lags the divergence by a quarter cycle and there is no measurable correlation. On the other hand, if a parcel has its vorticity reset to zero at the start of each updraft, then negative vorticity develops on rising trajectories, because they diverge at the top. In this case there is a correlation between divergence and negative vorticity, and that would be a sign of convection. Our feature-tracking approach yields vorticity and divergence at spatial scales of 200 km and larger. An entirely different approach 20 is to use infrared brightness itself as a dynamical variable, which extends the spatial scale down to wavelengths of ~15 km. The study 20 assumes that negative infrared brightness anomalies, which are related to cloud height, are upward displacements of pressure surfaces and therefore a measure of anticyclonic vorticity. This assumption is verified at scales from 250 to 1,600 km, which are smaller than the large cyclones but large enough that feature tracking is possible. It further assumes that the flow is quasigeostrophic, so the divergence is given by −1/f times the material rate of change of vorticity. This is the surface quasigeostrophic model 46,47 , which is used in meteorology and oceanography 48,49 . Applied to Jupiter 20 , one observes the signature of convection-a negative correlation between divergence and vorticity at 100 km scales and an upscale energy transfer from scales less than ~200 km to scales greater than ~200 km. These scales are just below the reach of our feature-tracking method. However, the quasigeostrophic approximation is based on Ro = ζ/f ≪ 1, where Ro is the Rossby number. The inequality is not strictly valid for Jupiter's polar cyclones (Fig. 4), and it is even less valid at smaller scales as Ro is predicted to increase with horizontal wavenumber 46 .
At mid-latitudes a downscale energy transfer is observed 50 , from wavelength scales of 2,000 km to the shortest scale measured, which is about 500 km. That study 50 used Cassini visible light imaging; ours used Juno infrared imaging, and the surface quasigeostrophic model uses infrared brightness. More work is needed to reconcile these three datasets. A parallel study 51 of Jupiter's south polar vortices, focusing on vorticity and stability, represents a step in the right direction.

Methods
A series of 12 images was started every 8 min to cover the same region at the north pole of Jupiter. Ideally the images in a series would fit together like tiles in a mosaic with no overlap and no spaces in between. The spacecraft was approaching Jupiter, and the image resolution changed from 22 km per pixel in the middle of the first series to 14 km per pixel in the middle of the fourth series, 24 min later. The two maps on the left of Fig. 2 were made by measuring cloud displacements between the first and third series, which are 16 min apart, and the two maps on the right were made from the second and fourth series, which are also 16 min apart. Therefore, the left and right maps are separated in time by only 8 min, but use entirely different images.
Supplementary Table 1 contains the archival filenames and our working names for the 48 images that were used in the analysis. The four series are named n01 to n04, each of which recorded roughly the same place on the planet 8 min after the one before. The first step in the processing was to determine the precise location on the planet of each resolution element in each image. This was done with NAIF/ SPICE data from the spacecraft and the precise geometric calibration of the JIRAM instrument. The second step was to map the brightness patterns onto a gridded reference plane tangent to the planet at the pole. We used 15 km per pixel for this mapping. These data are provided in Supplementary Data 1.   Under this assumption, the scaled PV at the pole is 4.18. The variations in PV are mostly due to variations in ζ and ϕ and less due to variations in f, which is the red line sloping gently down to the right. The other two curves (blue and green) were chosen to bracket the data and were used as initial conditions in a model study 21 that is described in the Discussion. The data points are given in Supplementary Data 3.
The third step was to measure cloud displacements in the reference plane using the Tracker3 software from JPL. The software automatically searches for the best correlation of brightness patterns between two images. This was done between images in series n01 and n03 and between images in series n02 and n04. Velocity was the displacement in kilometres divided by the time interval, which was always close to 16 min but depended on which image in each series was used. Correlation was done within a square box in the reference plane. After experimenting, we settled on a 15 × 15 pixel correlation box for the Tracker3 software. Thus, with 15 km per pixel in the reference plane, we were using squares 225 km on a side to define a feature. The resolution of the wind measurement was therefore ±112.5 km. We oversampled it by a factor of 2.5 to obtain wind vectors on a 45 km × 45 km grid. That dataset is Supplementary Data 2. We determined vorticity and divergence at every grid point by integrating around boxes of various sizes using Stokes's theorem and Gauss's law, respectively. Table 1 gives results for boxes 2, 4, 6 and 8 pixels on a side, corresponding to 90, 180, 270 and 360 km on a side, respectively.
The error in the velocity estimate σ v depended crucially on the granularity of the scene at the scale of the resolution element δ, which on average was about 18 km. Except for areas where there were no features at all, for which there were no estimates, the worst case was a single cloud feature of size ≤δ, for which the variance σ 2 v = 2δ 2 /Δt 2 , where Δt is the 16 min time step and the factor of 2 arises because we are subtracting position in two images. Then σ v = 2 1/2 δ/Δt, about 26.5 m s −1 . However, if the velocity measurement is the average of N statistically independent estimates of velocity, the variance is 2δ 2 /Δt 2 /N. The best-case scenario is when N is the number of resolution elements in the correlation box of L on a side such that N = (L/δ) 2 . Then σ v = 2 1/2 δ/Δt/N 1/2 = 2 1/2 δ 2 /L/Δt, which is 2.1 m s −1 for L = 225 km. Thus σ v is highly uncertain, but in Supplementary Table 2 we show that σ v = 7.81 m s −1 gives a good fit to the noise column in Table 1.
A quantitative measure of granularity is provided by image entropy H (ref. 52 ). We define it for each 15 × 15 correlation box from the histogram of brightness values in the box: The input data were 32 bit numbers, but we only had 225 pixels. We divided the range from the brightest to the darkest pixel into 256 grey levels, and we counted the number of times that each grey level appeared in the image. That number divided by 225 is p k , the frequency of occurrence of grey level k normalized so that ∑ p k = 1. The sum is over the 256 grey levels. If the brightness corresponding to a particular grey level k1 did not occur in the image, then p k1 = 0. At least 31 of the p k values must be zero. If all 225 pixels have brightnesses corresponding to grey level k2, then p k2 = 1 and all the other p k = 0, resulting in H = 0. If the brightness levels of all the 225 pixels are different, then H = log 2 (225) = 7.81. This is the maximum entropy for this problem. Low entropy is bad for feature tracking, and we experimented to find a value that eliminated the most suspicious data, such as the large pixel-to-pixel variations in the upper left and lower right corners of the divergence maps. We manually verified that the feature-tracking software was failing in those regions. Supplementary Fig. 5 is a histogram of entropy values, and Supplementary Figs 6 and 7 compare the vorticity and divergence maps with the low-entropy data present and with them masked out.
The data in Fig. 4 consisted of ~26,000 measured velocity vectors on the 45 km × 45 km grid of r < 6,010 km. Taking the azimuthal component v(r) of each vector, and knowing its r, we performed two separate least squares fits, one for a n and the other for b n in equation (2) to get analytic expressions for v(r) and ∂φ(r)/∂r, respectively. . (4) This choice of functions had no physical significance. The functions were chosen simply to fit the data and provide analytic expressions for integration and differentiation. For a good fit, the parameter r 0 must be close to the radius of the velocity maximum. From visual inspections, it was chosen to be 1,060 km for ∂φ/∂r and 1,200 km for v. We analytically integrated the expression for cyclostrophic balance in equation (2) to get φ (r) in Fig. 4, and we analytically differentiated the expression (1/r)∂(rv)/∂r =ζ to obtain vorticity. The measured azimuthal velocities are available in Supplementary Data 3.

Data availability
JIRAM data are available online at the Planetary Data System (PDS) at https:// pds-atmospheres.nmsu.edu/data_and_services/atmospheres_data/JUNO/jiram. The filenames of the images are listed in Supplementary Table 1. Calibrated, geometrically controlled radiance data mapped onto an orthographic projection centred on the north pole and velocity vectors derived from the radiance data are available in Supplementary Data 1-2.