A dynamo model of Jupiter’s magnetic ﬁeld

Jupiter’s dynamo is modelled using the anelastic convection-driven dynamo equations. The reference state model is taken from French et al. [2012]. Astrophys. J. Suppl. 202, 5, (11pp), which used density functional theory to compute the equation of state and the electrical conductivity in Jupiter’s interior. Jupiter’s magnetic ﬁeld is approximately dipolar, but self-consistent dipolar dynamo models are rather rare when the large variation in density and the effective internal heating are taken into account. Jupiter-like dipolar magnetic ﬁelds were found here at small Prandtl number, Pr ¼ 0 : 1. Strong differential rotation in the dynamo region tends to destroy a dominant dipolar component, but when the convection is sufﬁciently supercritical it generates a strong magnetic ﬁeld, and the differential rotation in the electrically conducting region is suppressed by the Lorentz force. This allows a magnetic ﬁeld to develop which is dominated by a steady dipolar component. This suggests that the strong zonal winds seen at Jupiter’s surface cannot penetrate signiﬁcantly into the dynamo region, which starts approximately 7000 km below the surface. 2014 Authors. Published Elsevier Inc. creativecommons.org/licenses/by/3.0/).


Introduction
Jupiter has the strongest magnetic field of any planet in the Solar System (Connerney, 1993). It is believed to be generated by convection-driven flows in the metallic hydrogen region of the planet (Parker, 1979;Stevenson, 1983Stevenson, , 2003. As the planet gradually cools, convection rather than radiation carries the heat flux out (Guillot et al., 2005), leading to an equilibrium reference state which is close to adiabatic. To model the dynamo, we use the self-consistent convection-driven dynamo equations in the anelastic approximation (Braginsky and Roberts, 1995;Lantz and Fan, 1999), which takes into account the large variation in density with depth. The model is based on an equilibrium reference state which uses an equation of state derived from density functional theory, and the electrical conductivity used here is also based on ab initio calculations (French et al., 2012). The model contains a small rocky core releasing less than 2% of Jupiter's intrinsic heat flux. As Jupiter cools, it releases an approximately uniform specific entropy everywhere outside the core, so the driving is different from the geodynamo, where the main buoyancy source is believed to be near the inner core boundary (basal heating).
In Boussinesq convection-driven dynamos, dipolar solutions occupy a large region of the numerically accessible parameter space (e.g. Olson et al., 1999;. Strong dipolar dominance is found for low E=Pm and moderate Rm, where E is the Ekman number, Pm the magnetic Prandtl number, and Rm the magnetic Reynolds number. Here dynamos are even more dipolar than the geomagnetic field . Jupiter's magnetic field is approximately dipolar, but strongly dipolar solutions for anelastic dynamos with large density ratios across the convecting shell are much harder to find (Gastine et al., 2012), a result confirmed here. Polytropic reference state models with uniform electrical conductivity only give dipolar solutions for density ratios less than 5 at Pr ¼ 1 (Gastine et al., 2012). At high density ratios the convective velocities are largest near the surface both in the linear (Glatzmaier and Gilman, 1981; and nonlinear ) regimes, and this enhances the helicity and zonal flow in the outer regions, leading to large small-scale fields there which dominate the dipolar component (Gastine et al., 2012). When the low electrical conductivity region in the non-metallic outer zone is taken into account, dipolar dynamos have been found in Boussinesq (Gómez-Pérez et al., 2010) and polytropic models , because the strong convection beyond the transition zone no longer generates disruptive small-scale fields. However, it was generally found that the transition zone between the electrically conducting region and the molecular insulating region must be in the range 0.7 to 0:8r jup  to get dipolar solutions. The recent ab initio calculations (French et al., 2012) suggest the transition zone is further out at $0:9r jup . These polytropic models were driven by basal heating; the more realistic uniform entropy source models compound the difficulty by enhancing the

Equations of the model
The Lantz-Braginsky-Roberts anelastic dynamo equations (Braginsky and Roberts, 1995;Lantz and Fan, 1999; were used in a spherical shell between the core radius r c and the cut-off radius r cut .
The usual form of the equation of motion in dimensional form is (e.g.  @u @t þ ðu Á rÞu ¼ 1 q j Â B À 2X Â u þ F m À rp 0 q À rU 0 À q 0 rU q ; ð2:1Þ and a perfect gas is often assumed, but the metallic hydrogen in the dynamo region of Jupiter means that the gas is far from perfect. Here u is the velocity in the rotating frame, j is the current density, B is the magnetic field, X is the rotational angular velocity, p is the gas pressure (including that coming from electron degeneracy in the metallic hydrogen region), q is the density and U the gravitational potential. The equilibrium state is assumed to be spherically symmetric for simplicity, but no symmetry is assumed for the disturbances (denoted by primes) produced by the convection. As usual in the anelastic approximation these disturbances are assumed not to alter the equilibrium density and pressure significantly (e.g. Lantz and Fan, 1999), and since the convective velocity in Jupiter is always much less than the sound speed this is a reasonable assumption. Hence the gravitational acceleration of the equilibrium state is g ¼ Àgr ¼ ÀrU and U is decomposed into U þ U 0 . We now put this equation into Lantz-Braginsky-Roberts form without assuming a perfect gas: see also (Ingersoll and Pollard, 1982;Braginsky and Roberts, 1995;Kaspi et al., 2009).
We can rewrite this in a more useful form using Maxwell's thermodynamic relations. The enthalpy H ¼ U þ p=q can be expressed in differential form as dH ¼ TdS þ dp=q:. So

@H @S
The equations are scaled using the units where r cut ¼ 6:7 Â 10 7 m is the cut-off radius, r c is the core radius, and the subscript m denotes values at the midpoint r ¼ r m ¼ ðr c þ r cut Þ=2. The small entropy drop across the layer is DS. The equation of motion (2.7) becomes, dropping the *, 1 Pm @u @t The entropy equation is ð2:14; 2:15Þ The dimensionless parameters are the Rayleigh number, the Prandtl number, the magnetic Prandtl number, the Ekman number and the radius ratio, which is the same in all runs described here.
In the entropy Eq. (2.11) S is the entropy per unit mass or specific entropy. We have assumed that turbulent entropy diffusion arising from small-scale motions dominates over radiative conductivity. Note that the viscous and ohmic dissipation terms appear in the anelastic heat equation. The source term PmH=Pr arises from the secular cooling of Jupiter which drives the convection. We assume that Jupiter evolves along a sequence of adiabats, so the specific entropy drops uniformly in space throughout the shell as time passes, i.e. H is independent of position. The value of H was chosen so that the heat flux coming out of the core is small compared to the internal heat flux coming out of the planet, as the mass of the core is so small it seems unlikely that it contributes significantly to the heat budget. The boundary conditions used were a stress-free outer boundary, a no-slip boundary at the core, constant entropy on both boundaries, and electrically insulating in the core and outside the cut-off radius. The small-scale convection that occurs in the low density region as the surface is approached (Glatzmaier and Gilman, 1981;, causes severe numerical difficulties. It is not practical to compute this small scale flow as well as integrating on the long time-scales needed for the dynamo. We therefore introduced a cut-off radius for the model which we took as r ¼ r cut ¼ 6:7 Â 10 7 m for the runs presented here, but different values were also tried. Some simulations were done using fixed flux rather than fixed entropy boundaries, as this can make a difference to Boussinesq dynamo models (Sakuraba and Roberts, 2009;Hori et al., 2010), but we found less difference in the anelastic cases.

The reference state
The reference state equilibrium model used is based on model J11-8a of French et al. (2012). Since the code we use is pseudospectral, the equilibrium quantities need to be smooth functions, so an interpolating analytic model was used for the reference state density q, temperature T and magnetic diffusivity g ¼ 1=l 0 r where l 0 is the permeability of free space and r is the electrical conductivity. The electrical conductivity above the cut-off level r ¼ r cut ¼ 6:7 Â 10 7 m is negligible, so ignoring this region will not hopefully influence the dynamo generation significantly. It is however possible that the zonal flow and the enhanced turbulence in the upper layers (Heimpel et al., 2005;Kuzanyan, 2009, Showman et al., 2011) could influence the dynamics below. The conductivity model used decays exponentially towards the surface. The electrical conductivity model (French et al., 2012) decays superexponentially, but the difference is only significant outside the cut-off region of our model, where the conductivity is negligible from a dynamo point of view. The reference state is assumed to be spherically symmetric. The quantities fed into the dynamo model are the temperature TðrÞ and temperature gradient T 0 ðrÞ, the density qðrÞ, the logarithmic density gradient nðrÞ ¼ d log q=dr; d n=dr, the magnetic diffusivity gðrÞ and its gradient d g=dr. The kinematic viscosity and m and entropy diffusivity j are taken to be constants.
The Leeds code uses dimensionless variables, and the unit of density is the value at the midpoint of the shell r m . Details of the reference state are given in Appendix A, where Table 2 gives the dimensional values of the core radius, midpoint radius, density, temperature and diffusivity units and length unit d. The density, temperature and diffusivity used are shown in Fig. 1a-c, together with the data points from the J11-8a model (French et al., 2012). The ratio of the density at the cut-off radius to that at the core boundary is 0.046. The discontinuity in density at r ¼ 4:4 Â 10 7 m in the J11-8a data arises from an assumed compositional change which has been smoothed out here.

Simulation results
The anelastic dynamo benchmark  showed that dipolar fields can be found for Pr ¼ 1, very large Pm and E ¼ 2 Â 10 À3 and Ra close to critical. However, very large Pm is not appropriate for giant planets, and we expect the convection to be strongly nonlinear. The most promising parameter regime found in our exploration of parameter space, where robust dipolar solutions occur, is at low Pr $ 0:1 and at Rayleigh numbers high enough to ensure a magnetic Reynolds number over 10 3 , and in consequence a strong magnetic field. Dynamos with long periods of dipolarity were found at Pr ¼ 0:15, but the dipole occasionally collapsed and regenerated with the opposite sign. At Pr ¼ 0:1 we found robust stable dipolar fields. The advantage of these nonreversing models is that they can be integrated to a statistically steady state with large but realistic computing resources.
The molecular thermal diffusivity and kinematic viscosity at r ¼ r m are given by French et al. (2012) as 1:7 Â 10 À5 m 2 s À1 and 3 Â 10 À7 m 2 s À1 respectively, so the Prandtl number there is 0.018. In the molecular H/He region it becomes order unity, but there radiative diffusion may be important keeping the effective Pr small. It is sometimes argued that turbulent diffusion by small-scale eddies will lead to a turbulent Pr of order unity, because the same small-scale eddies that transport momentum will also transport entropy. It is therefore likely that the effective Prandtl number will be greater 0.018, but how much greater is hard to estimate, so it seems sensible to explore Prandtl numbers in the range 0.018 to unity. The linear theory of rapidly rotating convection , Zhang, 1992 shows that at small Pr the pattern of convection at onset is somewhat more spread out, even at large density ratio, and also has a lower azimuthal wavenumber than at Pr ¼ 1 or larger. This helps generate dipole fields, because it spreads the dynamo action throughout the metallic hydrogen zone, rather than having it concentrated near the outer boundary.

Details of the runs
The dimensionless parameters used for the runs described here are given in Table 1. Runs A, B, C and D all have the same input parameters. Runs A, C and D all started from a strong dipole solution, and maintained the dipole throughout. The parameters Nr, Nth and Nphi refer to the number of grid points in the r; h and / directions respectively. The maximum spherical harmonic degree (h-direction) is 2Nth/3-1, and the maximum order (/-direction) is Nphi/3-1 because of de-aliasing. The h and / resolutions shown in Table 1 are defined as where E lm is the energy in the spherical harmonic of degree l and order m and E tot is the energy summed over all harmonics. h res is the energy in the last 5 spherical harmonics, i.e. those of degree L À 4 6 l 6 L, divided by the product of the total energy and the number of modes that contribute to the last 5 harmonics. The /-resolution is defined similarly but using the order of the spherical harmonics rather than the degree. Note that because m cannot exceed l, a much larger number of harmonics can contribute to the last five l harmonics than to the last five m harmonics. To compensate for this, in (4.1) we divide by the number of modes contributing as well as the total energy. This gives a convenient measure of how fast the higher harmonics are dropping off in the energy spectrum. Table 1 shows the results using the kinetic energy spectrum. In this region of the parameter space, the convergence of the magnetic energy spectrum and the entropy spectrum are similar in magnitude. Run C had the lowest resolution, and was run for over 2 diffusion times. The behaviour is very similar to that of run A. Run D has the highest resolution in h and was run for about 0.5 magnetic diffusion times. Note that the extra harmonics introduced in runs D and G by increasing Nth have dramatically improved the h-convergence, as we would expect, but have not reduced the /-convergence as much. The cumulative average of the dipole at t ¼ 1:2 was 0.40 for Run A and 0.41 for run D. This indicates satisfactory convergence in the spherical harmonic expansions. The differences between the Nr = 128 and Nr = 160 runs were minor. Checks were also made on the timestep controller, to ensure the timestep was small enough to make no major differences to long term averages. The timestep was normally around 10 À7 or slightly less in the presented runs. Fig. 2 shows results for the case Pr ¼ 0:1; Pm ¼ 3; E ¼ 2:5 Â 10 À5 ; Ra ¼ 1:1 Â 10 7 , details of the runs being given in Table 1. The magnetic energy, kinetic energy, and dimensionless heat flux, are defined as in Section 4 of . To show how key quantities are converging to a steady mean, cumulative averages are shown starting at t ¼ 0:2 magnetic diffusion times, to remove the effect of initial transients. It was noted in  that very long runs are needed to obtain very accurate values of the average energies, and the high resolution requirement and short timestep makes this impractical, but the cumulative average gives an idea of how the energies and dipole moment are approaching statistically steady values. The cumulative average over all available data for each run is given in Table 1, with a brief comment on the nature of each dynamo. Starting with a strong dipolar field, run A maintained the dipole for over 2 magnetic diffusion times, i.e. 6 viscous diffusion times and 60 thermal diffusion times. Run B has the same parameters as run A but the run started from a smaller dipolar field. There is bistability (Simitev and Busse, 2009;Gastine et al., 2012;Duarte et al., 2013) as the run B solution does not have a persistent dipolar field. However, at smaller E strongly dipolar solutions do emerge from a small seed field, for example in the case E ¼ 1:5 Â 10 À5 ; Ra ¼ 2 Â 10 7 as shown in Fig. 3a. A-dip is the Gauss coefficient g 10 at the outer boundary. A-quad is the g 20 Gauss coefficient, proportional to the quadrupole moment. It is always small and has zero long term average. Curve A-me shows the magnetic energy, and A-me its cumulative average, and similarly for B-me. The magnetic energy tracks the dipole coefficient, a result of the dipole dominance. Run B has a lower magnetic energy than run A. The cumulative averages indicate that the run is long enough for the time-averaged behaviour to emerge, though even longer runs would be needed to get very accurate values of these mean quantities (precise definitions given in Jones et al. (2011)). Fig. 2b shows the heat fluxes, kinetic energies and zonal flow . The internal heating coefficient has been chosen so the heat flux coming out of the core, A-coreflux, is very small as it will be in Jupiter. A-topflux, the heat coming out of the cut-off surface, is very similar in both runs, and is very stable (only run A shown). Run B has a much higher zonal flow and kinetic energy. This indicates that the magnetic field is controlling the zonal flow in run A. This control appears to be essential to get a dipolar dynamo for our model. In the dipolar regime the magnetic energy is about twice the kinetic energy but the scaling arguments below suggest a ratio of around 200 in Jupiter.

Time-dependence of the solutions
The results from runs E and F are shown in Fig. 3a and b. Run E has a lower Ekman number and higher Rayleigh number. The dipole was maintained for this run, though because the timestep is smaller it was only run for just over 0.6 diffusion times. Run F has the same parameter values but started from a small seed field. This run grew into a strong dipole, and ends up on the same solution as run E, so the bistability is lost by going to lower E. Run G was a higher resolution version of run E. Comparison in Table 1 shows convergence is still satisfactory, but not quite as good as in runs A-D, as might be expected from the smaller Ekman number.
It is possible that Jupiter's magnetic field reverses, like the geomagnetic field, but there is no observational evidence. A question of interest is whether there are Jupiter-like dynamo models which have reversals, like the geodynamo models (Glatzmaier and Roberts, 1995). The dipole coefficient time-series of Run B has reversals, but the surface radial magnetic field is never very Jupiter-like in run B. Run H, shown in Fig. 4, was the best reversing dynamo model found here. The Prandtl number is a little higher at 0.15, but the Ekman number is low enough to avoid bistability. This dynamo lies at the border of the steady dynamo window and the intervals between reversals are longer than in run B, and Table 1 Parameters for the runs. Ra; E; Pr and Pm are the Rayleigh, Ekman, Prandtl and magnetic Prandtl numbers respectively. H is the heat source, chosen in all runs so that the heat flux from the core is small. M.E. and K.E. are the magnetic energy and kinetic energy, Zon. Flow is the kinetic energy associated with the axisymmetric component of u/, and g 10 is the dipole coefficient of the magnetic field: these are all estimates of the time-averaged values in dimensionless units. The remaining quantities concerned with the numerical resolution are defined in Section 4.1. whenever the absolute value of the dipole moment is above around 0.15 the field is similar to that shown in Fig. 5 for run D. Since the magnetic diffusion time is several hundred thousand years (see Section 5 below), a Jupiter-like field could be maintained for many thousands of years according to this model before declining and reversing. However, the reversal behaviour in run H is somewhat different from that found in geodynamo reversals, where there are long periods of rather steady dipolarity with short intervals over which the dynamo reverses (e.g. Glatzmaier and Roberts, 1995;Kutzner and Christensen, 2004). Here the decline and build up are comparatively slow, and the dynamo has significant intervals of low dipolarity.

Run
4.3. The spatial structure of the solutions Fig. 5 shows snapshots taken from the highest resolution run D, with the same parameters as run A. Fig. 5a is the actual field of Jupiter (Connerney, 1993) on a spherical surface at the mean radius of Jupiter. Spherical harmonics of degree greater than 5 are not yet known reliably, and so are not included in Fig. 5a. Fig. 5c is a snapshot from the highest resolution run shown at full resolution, while Fig. 5b is the same solution but with all spherical harmonics of degree greater than 5 removed, to compare with Fig. 5a. The general similarity between Fig. 5a and b is notable. This is a snapshot, and there are times when the strong polar flux patches in the polar regions split in two, but they always recombine subsequently. Fig. 5d shows the axisymmetric part of the azimuthal field. Note that is essentially zero outside the metallic hydrogen region where the electrical conductivity is very small. If we take the ratio of the surface field strength to the internal axisymmetric azimuthal field strength as representative for Jupiter, then since the observed field strength at r cut is about 1.3 mT, the axisymmetric component of the azimuthal field would be about 6.5 mT. The highest strength field found anywhere in the model is about 12 times the surface value, about 16 mT. The antisymmetric nature of B / is also predicted for geodynamo models (Olson et al., 1999;. Fig. 6a shows the azimuthal flow on the outer boundary, from which we see that the axisymmetric part (the zonal flow) dominates the convective part at the equator, but not away from the equatorial belt. Fig. 6b shows the zonal flow (axisymmetric component of u / ) in a meridional plane. Strong zonal flow near the outer boundary at the equator is a robust feature of rapidly rotating convection (Zhang, 1992;Christensen, 2002;Heimpel et al., 2005; Jones and Kuzanyan, 2009) with a stress-free outer boundary. In our dipolar solution run D the zonal flow is mostly restricted to the low conductivity region, the flow in the metallic hydrogen region having little zonal flow, due to locking by the magnetic field.
There is a small transition region at the edge of the metallic hydrogen region where there is some zonal flow with non-negligible electrically conductivity, but this transition region seems to be too thin to affect the run D dynamo significantly. In Fig. 6c, the run B solution, the zonal flow is important in the metallic hydrogen region. Differential rotation then shears the convection columns, and a highly nonaxisymmetric nondipolar field pattern results. This separation of the zonal flow and toroidal field into distinct regions outside and inside the transition zone appears to be crucial to obtaining a dipolar Jupiter-like magnetic field. Fig. 6d shows the radial velocity in the equatorial plane. The convection columns outside the metallic hydrogen region are disconnected from those inside. Movies show that the small scale convection in the currentfree region is strongly advected by the zonal flow, whereas the convecting columns in the metallic hydrogen are not significantly sheared over their lifetime. At Pr ¼ 1 the convection outside the metallic hydrogen region is on a smaller scale than the convection in the magnetically influenced region (Gastine et al., 2012) but at Pr ¼ 0:1 the dominant azimuthal wavenumber is not so different. Gastine et al. (2012) argue that the vigorous small scale convection in the outer regions makes large scale dipolar dynamo action problematic, so the larger scales at Pr ¼ 0:1 may be connected with the existence of the dipolar window there. The meridional section Fig. 6e of u r shows that the flow is columnar, though the columns  do not reach right across the planet. Columnar convection appears to be essential for dynamo action dominated by the dipolar component (Olson et al., 1999;Sreenivasan and Jones, 2011). Fig. 7 shows a sequence of snapshots of the axisymmetric part of the azimuthal field at different times taken from run B, a non-dipolar run. The sequence shows evidence of a dynamo wave progressing from pole to equator, as happens in the solar dynamo. In Fig. 7a, the field is dominated by negative field in the northern hemisphere and positive field in the southern hemisphere, the antisymmetric form of the azimuthal field being consistent with a dipolar dynamo generated field. However, in Fig. 7b high-latitude reversed field starts to grow, and by the time of Fig. 7c the original azimuthal field is being squeezed by these reverse flux patches moving towards the equator. By the time of Fig. 7d, the original flux patches have gone, and the field is approximately reversed from Fig. 7a. The cycle then repeats, with faint patches of the original field parity now visible at high latitudes. The sequence shown in Fig. 7 has been chosen because the dynamo wave is quite clearcut, but in general the dynamo waves are rather erratic, as can be anticipated from the run B plots in Fig. 2a, superimposed on chaotic field fluctuations typical of high Rm numerical dynamos. However, the radial component of the field is consistent with a dynamo wave interpretation. Dynamo waves were seen by Duarte (2014) in a compressible Jupiter model dynamo at Pr ¼ 1, though these travelled from equator to pole.
The parameter space is large, and it has not yet been fully mapped. However, at Pr between 0.1 and unity with uniform heating, the run B behaviour was frequently found at relatively low Rm. At larger Rm, the flow is less columnar, and the dynamo is small scale, as found in Boussinesq models when the local Rossby number becomes too large (Sreenivasan and Jones, 2006;Olson and Christensen, 2006). As the Prandtl number is reduced towards 0.1, the weak high latitude reversed flux patches seen in Fig. 7 no longer grow, though there is a faint trace of them in Fig. 5d.
In Boussinesq dynamo models, changing from a fixed temperature outer boundary condition to a fixed flux outer boundary condition can make a significant difference to the form of the magnetic field and the convection (Sakuraba and Roberts, 2009;Hori et al., 2010). This seems to be less true in these anelastic dynamo models, but to explore this effect we present in Fig. 8 typical snapshots from runs E and I, which differ only in that in run I the flux is fixed at the outer boundary to the average value found in run E (see Fig. 3b). This means of course that the entropy is no longer fixed at zero there, and in run I the poles were slightly cooler and the equator slightly hotter. Since the convection is driven by very small entropy changes (Jupiter's interior is close to adiabatic) this corresponds to only a very small pole-equator temperature difference, well below any observational constraints. Fig. 8b, d and f correspond to the fixed flux case, Fig. 8a, c, and e to the fixed entropy case. There is no great difference between the two cases, but the dipole is slightly stronger in the fixed flux case, and this is the case for most times, though occasionally the dipole in run E will exceed that in run I. The zonal flow in run I is more confined to the equatorial region as we would expect from the stronger dipole leading to more efficient locking of the zonal flow. In Fig. 8e and f the axisymmetric radial magnetic fields are compared. There is more reversed flux near the equator in the fixed entropy case, and this leads to a belt of slightly weaker radial field near the equator in Fig. 8a compared to Fig. 8b; although this feature is not very striking, it is persistent.

Scaling of the dimensionless units
Estimates of the velocity of convection in the deep interior of Jupiter varying from 10 À3 m s À1 to 10 À2 m s À1 are derived from scaling laws (Christensen and Aubert, 2006;Showman et al., 2011) using the surface convective flux which is (Guillot et al., 2005) 5.4 W m À2 . The uncertainty arises because an extrapolation over many orders of magnitude must be made from the numerically accessible parameters to those obtaining in Jupiter and the exact power laws governing rapidly rotating convection remain uncertain. An alternative would be to use the secular variation of Jupiter's magnetic field, a successful technique for the Earth (Holme, 2007). There is still uncertainty over whether the secular variation has been unambiguously observed (Yu et al., 2010). Nevertheless, the magnetic field data acquired by NASA over the period 1973-2003 is better fitted by a model with a 0.012% per year increase in g 10 over this period (Ridley, 2010). The root mean square _ g 10 =g 10 for run A was 41 in our dimensionless units, and identifying this with 0:012% per year sets our dimensionless time unit s ¼ d 2 =g m at s ¼ 340; 000 years. Since d ¼ 6:055 Â 10 7 m, this means our effective turbulent magnetic diffusivity at the shell midpoint is d 2 =s ¼ 342 m 2 s À1 , about 700 times the laminar value. Our unit of velocity d=s is then 6 Â 10 À6 m s À1 . The kinetic energy for run A averages out at 1:5 Â 10 6 , suggesting a typical velocity of 1:7 Â 10 3 near the midpoint of the shell in our dimensionless units which corresponds to 10 À2 m 2 s À1 , at the top end of the range deduced from the scaling laws. Note that even if the convective velocity is at the bottom of the suggested range, 10 À3 m s À1 , the magnetic Reynolds number based on the laminar magnetic diffusivity is 10 5 , out of reach numerically. The zonal flow in run A rises to around 5000 in dimensionless units near the equator at our cutoff depth 3000 km below the surface, corresponding to only 3 Â 10 À2 m s À1 , much less than the surface wind speed at the equator. However, the reason that the zonal flow in the low conductivity region only rises to this comparatively small value is that the viscosity is much greater in our numerical model than in Jupiter. Lower values of E show the zonal flow increasing, but it is not yet numerically possible to get E very much lower than 10 À5 in a full dynamo simulation and this is not low enough to get zonal flows of the correct speeds, but simpler models (e.g. Jones et al., 2003) do show faster zonal flows. Our model strongly suggests that intense zonal flows cannot penetrate into the metallic hydrogen zone, and so the alternating jets poleward of latitudes ±26°must be ageostrophic. These ageostrophic jets may not penetrate deep into the interior (Vasavada and Showman, 2005), but the winds between ±26°(which are much stronger) are geostrophic, and so can be approximately determined by the surface flow. This predicts a gravitational signature due to Jupiter's internal rotation which may be measured by the Juno mission (Kong et al., 2013). Comparing Fig. 5a and b suggests that our dimensionless unit of magnetic field corresponds to 1.3 mT. The azimuthal field in the deep interior rises to about 5 times this value, 6.5 mT, not far from the 8 mT predicted by scaling law arguments (Christensen and Aubert, 2006;Yadav et al., 2013). If the turbulent value of g m ¼ 342 m 2 s À1 , is adopted then ðXlq m g m Þ 1=2 (the unit of nondimensionalisation) is about 13.9 mT, roughly ten times too high, but 1.3 mT lies between the values of ðXlq m g m Þ 1=2 (the unit of nondimensionalisation) given by the laminar and turbulent values of g, closer to the laminar value. The fact that the scaled field is too large with respect to the turbulent magnetic diffusivity is no great surprise, as scaling laws are needed to translate simulation estimates into planetary estimates.
The thermal diffusivity (more precisely the entropy diffusivity) and particularly the viscosity in the model are much greater than in Jupiter (French et al., 2012), so we have to interpret all our diffusive processes as being enhanced by sub-grid scale turbulence. This makes it difficult to compare the model heat flux with the actual heat flux of Jupiter, so the numerical models must be interpreted in the light of scaling laws extrapolating them to the realistic parameter regime. The maximum Rossby number in our calculation is about 0.01, which is actually about the right value for 100 m s À1 zonal flows, but this agreement is rather artificial; if we chose Xd as our unit of velocity the maximum zonal flow would indeed be 100 m s À1 , but the typical convective velocity in the model would be 20 m s À1 , much too high. It is not possible with current computer resources to capture realistic convective velocities and realistic surface zonal flow speeds in the same model Showman et al., 2011).

Conclusions
After many failed attempts, Jupiter-like magnetic fields have been found using self-consistent convection-driven dynamo models with a state-of-the-art reference model of Jupiter and a plausible internal heating model. If the Juno mission succeeds in gaining a more detailed representation of Jupiter's field, it will be interesting to compare the new data with the predictions of dynamo models such as this. Also, the data expected from the measurements of the gravity field will give information relevant to the dynamo process, because the models described here are incompatible with a differential rotation of the order of magnitude of the surface flow extending into the metallic hydrogen zone. Such strong differential rotation would break down the dipolar form of the field. On the basis of these models, we expect the geostrophic regime to be limited to the equatorial region, and only to extend to the region between latitudes ±26°at most. However, our model does suggest that the zonal flow starts to develop near the equator just below where the conductivity drops to low values, which is 7000 km below the surface, much deeper than in some shallow 'weather layer' models.
A clear message arising from the anelastic models, both those of Gastine et al. (2012) and ours, is that dipole-dominated solutions are not easy to find. It is possible that at very low Ekman numbers (less than 10 À5 ) dipolar solutions exist in a larger region of parameter space, but unfortunately exploring below 10 À5 is too expensive using the numerical techniques currently available. At present low Pr (typically around 0.1) and Pm > 1 seem the most promising values; even lower Pr $ 0:03 was used by Glatzmaier (Stanley and Glatzmaier, 2010). We also found that the Rayleigh number had to be large enough to give a magnetic Reynolds number of order 10 3 in order to generate fields strong enough to suppress the differential rotation in the dynamo region. However, we note that gas giant dynamos are as yet relatively unexplored compared with geodynamo models.
An important area for further investigation is the role of the cutoff radius. Although in the model the cut-off occurs in the semiconducting region above the transition zone, it is possible that if the depth of the cut-off was reduced below 3000 km the turbulence level would be enhanced, and this enhanced turbulence might find its way into the dynamo region and hence affect the dynamo process. It may prove to be important for there to be relatively little dynamical coupling between the molecular hydrogen region near the surface and the dynamo region below to get stable dipolar solutions. Low coupling would also help to explain how a strong geostrophic differential rotation in the equatorial region can build up independently of the low zonal flow in the dynamo region.