JUPITER MODELS WITH IMPROVED AB INITIO HYDROGEN EQUATION OF STATE (H-REOS.2)

Our new hydrogen EOS (H-REOS.2) is constructed in the same way as the former version H-REOS.1 (N08). The H EOS is assembled from different contributions at the low-density range (ρ < 0.1 g cm⁻³) and the high-density range (0.2 g cm⁻³ ⩽ρ ⩽ 9 g cm⁻³). At low densities, pressure dissociation and pressure ionization do not play an important role for the EOS. For these densities, models within the chemical picture agree well with experimental data. As for H-REOS.1, we here apply the Fluid Variational theory (FVT⁺) model (Holst et al. 2007) where ionization is neglected for temperatures below 3000 K.

In the high-density range, our EOS table is obtained from finite-temperature density-functional theory molecular dynamics (FT-DFT-MD) simulations using the code VASP (Kresse & Hafner 1993, 1994; Kresse & Furthmüller 1996). There are three improvements of our new H EOS. First, due to the enormous increase in computing power we were able to perform the simulations with 256 particles in a box compared to 64 particles (Holst et al. 2008). Second, we performed them up to 50,000 K, in order to cover the interior of a warm, young Jupiter, compared to 20,000 K (Holst et al. 2008; N08). In the simulations, a quantum mechanical treatment of the molecular vibrations of the ions is not included. Thus third, after completion of the simulations, we add the energy of a quantum mechanical oscillator, u_vib = k_BΘ_vib(0.5 + (exp(Θ_vib/T) − 1)⁻¹), to the internal energy of each molecule and subtract its classical energy of vibration, k_BT. The degree of dissociation is estimated via the coordination number (see Holst et al. 2008). This contribution is non-negligible if the temperature of the system is below the vibration temperature Θ_vib (French & Redmer 2009). For hydrogen, we use the experimental value Θ_vib = 6338.2 K.

For H-REOS.1, we obtained convergence of the EOS data within 3% and estimated the overall uncertainty due to both statistical and systematic errors to be below 5%. The new ab initio data in H-REOS.2 have an uncertainty of 1% for most of the density–temperature space. They are within the uncertainty of the former data for ρ ⩾ 0.3 g cm⁻³ but shifted to 2% higher pressures in a T–ρ region where pressure dissociation and ionization occurs. Figure 1 shows selected isotherms. With the onset of pressure dissociation at 0.2 g cm⁻³, the difference rises to 10%. This is a rather small density for FT-DFT-MD simulations that requires a large computational effort to reach converged results. The displayed region of 5000–10,000 K and 0.2–1 g cm⁻³ is relevant for the Jupiter adiabat. Here, the SCvH-i hydrogen EOS (Saumon et al. 1995) shows systematically higher pressures. At densities of 0.1 g cm⁻³ ⩽ρ ⩽ 0.2 g cm⁻³, the difference in pressure is a consequence of a revised interpolation between the low-density regime and the data points at 0.2 g cm⁻³. This difference is particularly pronounced for the internal energy (not shown). Along the Jupiter adiabat (Figure 2), the revised internal energy requires higher temperatures in order to keep the entropy constant at the level defined by Jupiter's outer boundary condition of 170 K at 1 bar (see Section  2.3), and also leads to a decrease in the density at low pressures of 0.1–10 GPa. However, due to little mass prevalent at those pressures in Jupiter, this region does not affect our Jupiter interior models. Instead, it is the 10–50 GPa pressure region where the lower densities of H-REOS.2 map onto the new Jupiter models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our EOS data tables provide the thermal EOS P(ρ, T) and the caloric EOS u(ρ, T), but not the specific entropy s(ρ, T). In Appendix A, we describe our method to derive an adiabat through a given point (P, T).

Careful inspection of our code used for the calculation of Jupiter's shape and rotation in N08 revealed that two coefficients (one of second order and the other one of third order) contained errors. Their effect on Jupiter models with H-REOS.1 is twofold: (1) the layer boundary shifts from ∼8 to ∼4 Mbar when the outer envelope is made to have solar metallicity and (2) the metallicity increases by ∼0.02 (in absolute value) in the outer envelope and ∼0.035 in the inner one so that the maximum total mass of metals increases from 32 to 41 M_⊕, where the latter value was presented in N08. All other results are almost unchanged. The two errors were found by a comparison of the resulting gravitational coefficients J_2n, n = 1–4, of linear and polytropic density Jupiter and Saturn models with analytic solutions given in Zharkov & Trubitsyn (1978), and with numerical solutions for these density distributions by T. Guillot (2008, personal communication). For details, see Nettelmann (2009). Another improvement of the code is the application of a Newton-Raphson scheme to find the central conditions that meet the envelope pressure and mass at the core–mantle boundary. For the current version of our code, we here introduce the label MOGROP-11 (MOdellierungsprogramm für GROße Planeten, 2011).

Apart from the code improvements mentioned above, our Jupiter interior and thermal evolution models are generated by essentially the same code as introduced in N08, Nettelmann (2011), and Fortney & Nettelmann (2010), where in the latter case the code was applied to the structure and evolution calculations of Uranus and Neptune. We here recall the properties of our three-layer structure model for giant planets, review the observational constraints for the new Jupiter models, and describe the method of structure and evolution modeling.

Three-layer structure model. A sufficient number of parameters to meet the observational constraints and boundary conditions are available if we assume a three-layer structure with two envelopes and a core. The envelopes are adiabatic, homogeneous, and consist of hydrogen, helium, and heavy elements, while the core is made of rocks. Heavy elements in the envelopes are represented by the water EOS H₂O-REOS (French et al. 2009; N08), for helium we use He-REOS (Kietzmann et al. 2007; N08), and for rocks the formula given in Hubbard & Marley (1989). The outer envelope has a lower He abundance (Y₁) than the inner envelope (Y₂). We allow for different heavy element mass fractions in the outer envelope (Z₁) and in the inner envelope (Z₂). The transition pressure between the envelopes P_{1 − 2} is a free parameter. The density and entropy are discontinuous at layer boundaries.

Observational constraints. The models are required to meet the observational constraints for the total mass M_J, equatorial radius at the 1 bar pressure level R_J, 1 bar temperature T₁, solid-body rotation rate ω, gravitational moments J₂, J₄, J₆, atmospheric helium abundance Y_atm, and have the same He mass fraction² Y as the protosolar cloud once had. As long as the oxygen abundance below Jupiter's water cloud deck is not measured, only the lower limit for the heavy element abundance Z_atm in Jupiter's atmosphere is known. This value is 1 × Z_☉ (see Section  4), where Z_☉ is the protosolar cloud's metallicity of 1.5% (Lodders 2003).

These observational constraints for the structure models can be divided into four groups depending on whether they enter the modeling procedure explicitly or need to be fitted by adjusting model parameters, and whether or not they are varied within the 1σ error bars. The first set (explicit, varied) is empty. The second set (explicit, not varied) consists of M_J = 317.8338 M_⊕, 2π/ω = 9^h55^m30^s (the period of the magnetic field), the upper limit T₁ = 170 K, and Y_atm = 0.238. While the small uncertainties in M_J and ω would not affect the resulting models, lowering T₁ by 5 K would reduce Z₁ by ∼0.3 Z_☉ (N08), in the opposite direction to what we are aiming for. Modifying Y_atm by 1σ (±0.007) changes Z₁ by about ∓0.3 Z_☉. Because Y_atm was measured by the Galileo entry probe at a depth where vertical convection ensures homogeneous distribution of those elements that do not form clouds (such as the noble gases) deeper inside, we set Y₁ = Y_atm. The parameters of the third group (to be fitted, not varied) are R_J = 7.1492 × 10⁷ m and Y = 0.275, where R_J is fitted by choosing the mean radius R_m of the equipotential surface that coincides with the 1 bar pressure level and Y by adjusting Y₂. We find R_m = 10.955 R_⊕ and Y₂ = 0.285–0.325. The gravitational moments J₂, J₄, and J₆ build the fourth group (to be fitted, varied). In particular, Z₂ is used to reproduce J₂/10⁻⁶ within the tight observational bounds 14,697(1), and Z₁ is used to adjust J₄/10⁻⁶ to either −587 (the observed central value) or −589 (the observed lower limit) (Jacobson 2003). When J₂ and J₄ are matched, J₆/10⁻⁶ is found to lie within 34.0 and 37.5, consistent with the observed limits (29.1–39.5; Jacobson 2003). Our computation of Jupiter's thermal evolution makes use of the observed effective temperature T_eff = 124.4 ± 0.3 K.

Modeling procedures. Internal profiles of m, P, T, and ρ which depend on the radial coordinate l are obtained by numerically integrating the structure equations

$$\begin{equation} \rho ^{-1}\:dP/dl=-d(V+Q)/dl \end{equation} \tag{ 1 }$$

and

$$\begin{equation} dm/dl = 4\pi \: l^2\:\rho . \end{equation} \tag{ 2 }$$

Equation (1) states that the force acting on a mass element with density ρ at radius l is balanced by a pressure gradient, where the force −d(V + Q)/dl arises from the gravity field V and the centrifugal force with potential Q due to rigid-body rotation. Equation (2) expresses the amount of mass dm included in a spherical shell of density ρ and thickness dl at radius l. These two ordinary, first-order differential equations are integrated inward starting at the surface l = R_m with the two boundary conditions P(R_m) = 1 bar and m(R_m) = M_J. In order to satisfy in addition the inner boundary condition m(0) = 0, a further parameter is needed. For this we choose the core mass M_c. Thus, for given values of the observational constraints and of the transition pressure P_1-2, the result is one single Jupiter model, obtained in terms of the internal profile and values for Z₁, Z₂, and M_c. To calculate the planetary shape and the gravitational moments, we apply the theory of figures (Zharkov & Trubitsyn 1978) up to third order as in N08 and also to fourth order. Such structures describe Jupiter at present time t₀ = 4.56 Gyr.

For modeling Jupiter's thermal evolution, we apply a procedure similar to that of Saumon et al. (1992). To describe Jupiter's interior at earlier times t  <  t₀, we select a few representative structure models and generate for each of them a sequence of ∼80 warmer interior profiles by increasing T₁ up to 800 K while m(P_1-2), M_c, and the element abundances are kept at constant values. Further models with intermediate T₁ values are then generated by interpolation. While for the structure models at the present time with known R_J(t₀) the core mass was chosen to ensure mass conservation, we here invert the problem to find the planet radius R_m(t) for the given core mass. Unlike Saumon et al. (1992), rotation is included in the approximation of spherical symmetry where the centrifugal force reduces to the zero-order term −dQ/dl = (2/3)ω²l. We calculate the cooling curves in the usual approximation of constant angular velocity (Saumon & Guillot 2004; Fortney et al. 2011), and then also by including angular momentum conservation (Hubbard 1970), which requires just one to three additional iterations for a given value of T₁, and the corresponding change in the energy of rotation, dE_rot. The time dt passed between profiles with different surface temperatures and hence different internal entropies is given by the equation of energy balance,

$$\begin{equation} L_{\rm eff}- L_{\rm eq}= -dE_{\rm int}/dt, \end{equation} \tag{ 3 }$$

where L_eff is the planet's observable luminosity in the infrared and L_eq is the luminosity the planet would have in equilibrium with the insolation. Using the Stefan–Boltzmann law for the energy emitted by a blackbody with uniform temperature, these luminosities can be used to define the commonly used temperatures T_eff, T_eq, and T_int. This way, T_eq(t₀) is derived from the L_eq value given in Guillot (2005). It is set constant over time or, alternatively, varied with time according to a linear increase of the irradiation flux F_eq := L_eq/4πR²_m, starting with F_eq(0) = 0.7F_eq(t₀) as predicted by theoretical standard models for the Sun (Bahcall & Pinsonneault 1995; Sackmann & Boothroyd 2003). Note that the solar standard evolution model predicts an increasing bolometric luminosity with time, whereas the short-wavelength (chromospheric and coronal) emissions of the young Sun may have been up to 1000 times stronger than those of the present Sun as indicated by the measured X-ray and ultraviolet (XUV) fluxes of young solar-like stars (Ribas et al. 2005). Absorbed XUV irradiation can lead to heating of and mass loss from the upper planetary atmosphere (Lammer et al. 2003) but absorption at high altitudes is not suspected to influence the energy balance of the interior (Guillot & Showman 2002). Therefore, we ignore the activity-related, time-dependent solar XUV flux in our evolution calculations for Jupiter.

The difference between the luminosities radiated away and absorbed equals the loss of the planet's intrinsic energy dE_int per time, see Equation (3). In the envelope, the energy lost per mass shell simply is the heat δq = T(m) ds(m). Guillot et al. (1995) offer a convenient closure relation by relating T_eff to T₁, namely, T₁ = K T^1.244_eff g^−0.167 according to the model atmosphere grid by Graboske et al. (1975), where the parameter K can be used to reproduce the observed T_eff(t₀) and g is the surface gravity. Also like Guillot et al. (1995), we take into account the core's contribution to the planet's intrinsic luminosity. The heat loss of the core is M_c c_v dT_core, where we use c_v = 1 J K⁻¹ g⁻¹ as in Guillot et al. (1995) and dT_core is the temperature difference of the isothermal cores of subsequent interior profiles with different T₁ values, as is ds(m) the entropy difference at mass shell m between those interior profiles. Although the abundances of radioactive elements decrease exponentially with time so that the energy production from radioactive decay in rock with meteoric composition may have been an order of magnitude larger in the past than it is today, we here use the present Earth's value L_radio  =  2 × 10¹³J s⁻¹ M⁻¹_⊕ to account for this contribution. This treatment seems justified for Jupiter because we find the prolongation of the cooling time due to the core's heat loss to be 0.01 Gyr only. Thus, the cooling equation reads

$$\begin{equation} dt = -\frac{\int _{M_{\rm c}}^{M_{\rm J}}dm\:T\,ds + M_{\rm c}c_v dT_{\rm core} + {d}E_{\rm rot}}{4\pi R^2\sigma \big(T_{\rm eff}^4-T_{\rm eq}^4\big) - L_{\rm radio}}. \end{equation} \tag{ 4 }$$

After integrating Equation (4) backward in time, we obtain the cooling time τ that the particular Jupiter model needs to cool down from an arbitrarily hot initial state to the present state.

Results for the core mass, the outer envelope metallicity, and the inner envelope metallicity as functions of the depth of the layer boundary between the envelopes are shown in Figure 3. If the layer boundary is located higher in the planet (lower transition pressures), Z₁ and Z₂ become smaller. The lighter envelopes then require a larger core mass to conserve the total mass. When Z₁ decreases down to zero, the core mass cannot grow any further. However, we do not consider models with Z₁  <  Z_☉ acceptable as this is inconsistent with Jupiter's observed atmospheric abundances. We find Z₂ ≫ Z₁ and M_c  <  10 M_⊕ for all our acceptable three-layer models. As the gravitational moments are most sensitive to the density in the outer part of the planet at pressures of a few Mbars (J₂) or even below 1 Mbar (J₄) and increase with the local density, Z₁ must decrease when the denser inner envelope extends farther out. Otherwise, |J₄| would become too large. On the other hand, as we go to higher pressures deeper inside the planet, the sensitivity of J₄, and later also of J₂ approaches zero. Therefore, Z₁ changes only weakly with P_1-2 for deep internal layer boundaries. However, since J₂ is adjusted by Z₂, Z₂ must rise strongly with P_1-2 in order to provide a sufficiently high mass density that is able to keep J₂ on the high level of the observed value. With rising metallicity the envelope becomes denser, so the core mass must decrease. When M_c = 0 is reached, the layer boundary is as deep as it can be, and Z₁ adopts a maximum value. For models with H-REOS.1, this maximum is 1.0 Z_☉ if J₄/10⁻⁶ is adjusted to −587, and 1.5 Z_☉ if J₄/10⁻⁶ is adjusted to −589. For models with our improved H EOS, this maximum is 2.5 Z_☉ (J₄/10⁻⁶ = −587) and 2.7 Z_☉ (J₄/10⁻⁶ = −589), respectively. This enhancement in Z₁ arises from the ∼2% higher pressures at given hydrogen mass density (Figure 1), implying ∼2% lower hydrogen mass density at given pressure. The lower partial density of hydrogen in the mixture requires it to be compensated for by a similar amount of metals.

Zoom In Zoom Out Reset image size

Including fourth-order coefficients to the calculations of Jupiter's shape and rotation systematically shifts the solutions to ∼0.8 Z_☉ lower outer envelope metallicities, inducing an additional uncertainty on Z₁ that is of the same size as the error bars of single observables such as T₁ and J₄; see below. For M_c and Z₂, the inclusion of higher-order terms has an even greater effect than the improvement in the H EOS. Obviously, a convergence check of the solution with increasing order of the theory is highly important but beyond the scope of this paper. The accuracy of our current computer code is insufficient to reliably study the effect of higher than fourth-order coefficients. From preliminary calculations, solutions with fifth-order coefficients are found to lie between those of third and fourth order.

Figure 3 also shows the observed element abundances in Jupiter's atmosphere. Assuming that O is as equally enriched as C, these abundances indicate an average enrichment of 2 (noble gases) to 4 (C, N, O). Taking two times solar as the lower limit for Z₁, models with P_1-2 < 4 Mbar drop out of the realm of acceptable models, reducing the upper limit for the core mass down to 8 M_⊕. The total mass of metals is M_Z = 29–32 M_⊕.

As representatives for our acceptable models we recommend models J11-4a and J11-8a, highlighted in Figure 3. Model J11-4a is the one with the outermost possible layer boundary, i.e., at 4 Mbar. Model J11-8a has P_1-2 = 8 Mbar for which Z₁ is ⩾2 Z_☉ for all of the above considered uncertainties. A machine-readable table for model J11-4a is provided as supplemental material; see Table 1. For these two models, we also vary Y, Y₁, T₁, and J₄ within their 1σ error bars. Modifying Y by ΔY = ±0.05 influences Y₂ in the same direction, and thus leads to a change in Z₂ of ∓8%. Modifying Y₁ by ΔY₁ = ±0.05 causes a response in Z₁ of ∓10%. Similarly, a 5 K colder 1 bar level leads to a 50% decrease in Z₁. In both cases, Z₂ consequently changes by ±3%–4% in the opposite direction to Z₁, and M_c changes by ∓2%–3%, again opposite to Z₂. The colder outer envelope adiabat also cools the interior, resulting in a core that is 100 K colder. The slightly denser H–He subsystem that results reduces the amount of metals in the envelope required to match J₂, but the net effect of ΔT₁ on Z₂ remains positive. Setting |J₄|/10⁻⁶ down to 585 lowers Z₁ by ∼15%. Clearly, the response of Z₂ (increase) and thus of M_c (decrease) in this case is stronger than in case of ΔT₁ and ΔY₁. Including these uncertainties, the total mass of metals changes but insignificantly to 28–32 M_⊕.

The striking difference between the Militzer et al. (2008, hereafter M08) predictions for the core mass of Jupiter (16 ± 2 M_⊕) and that of our group (0–8 M_⊕), where both groups applied similar ab initio simulations to generate the EOS data, has prompted speculations about causes for that difference (Militzer & Hubbard 2009). We here investigate how this difference relates to the modeling assumptions of two-layer (M08) or three-layer (N08; this work) structures. For that purpose, we have calculated two-layer models following M08. As in M08, the homogeneous envelope has Y₁ = 0.238, and the missing helium to obtain an average Y of 0.275 is not explicitly accounted for but may be part of the core mass. The core is either pure rocks or rock–ice, where we switch from ice to rock at 70 Mbar. As a result, only the innermost 4 M_⊕ are rocks. Because these models do not have the degree of freedom Z₂, we can only match either J₂ or J₄. As in M08 we choose to match J₂ by varying the envelope metallicity Z. We apply our usual numerical procedure (see Nettelmann 2011 for details) and vary Z, starting from 0.01, until J₂(Z) meets the observed value. The core mass must decrease with increasing envelope metallicity (denser envelopes) in order to ensure total mass conservation. If the number of layers assumed for the structure of the model were responsible for the different Jupiter models, we would expect to obtain the same two-layer model as in M08. This is not the case, as illustrated in Figure 4. While the M08 models require 4 ± 2 M_⊕ heavy elements in the envelope to meet J₂, corresponding to Z = 1.0 ± 0.5 Z_☉ for a 16 M_⊕ core, our two-layer models require Z = 4.3 Z_☉.

Zoom In Zoom Out Reset image size

This significantly larger envelope metallicity of our two-layer models enables us to make a transition to three-layer models where metals are shuffled from the outer part of the envelope to the inner part in order to get |J₄|/10⁻⁶ down to the observed value (from 614 to 587). The M08 models do not have this degree of freedom as the already low envelope metallicity would become zero long before J₄ is matched. Hence, the M08 models stay at J₄/10⁻⁶ = −614 and a large core, while our two-layer models can undergo a transition to three-layer models where Z₁ becomes smaller and Z₂ becomes larger than in the homogeneous envelope case. As a consequence of large Z₂ values, the core mass decreases further. In contrast, the core mass of the two-layer models is large because this mass also contains the missing mass of helium dM_He = (0.275–0.238)(M_J − M_c) ∼ 11 M_⊕ needed for an overall abundance Y = 0.275. Thus, the ice–rock core mass of two-layer models is about 5–10 M_⊕, just at the upper limit of our three-layer models. This estimate would also hold for the two-layer model of M08. However, we cannot reproduce their Jupiter model, the reason for which remains unexplained at the moment. We conclude that the difference between the M08 two-layer models and our three-layer models does not originate from assumptions about the number of layers, but instead is a symptom of already different models within the two-layer frame.

We have calculated the moment of inertia I as described in Appendix B. The non-dimensional form (λ) is obtained as λ = I M⁻¹_J R_m⁻². For all of our three-layer models, we find an almost invariant value λ = 0.27605 ± 0.03%. If we ignore Jupiter's shape deformation, the resulting λ decreases by 4% down to 0.26539 ± 0.03%. This is close to the predictions of 0.2629−0.2645 of Helled et al. (2011), who considered a broad set of interior models with core masses between 0 and 40 M_⊕ as allowed by the use of six-order polynomials to represent the pressure–density relation in Jupiter's envelope. From our calculations, we conclude that λ is not an appropriate parameter to constrain Jupiter's core mass further. In contrast, Helled et al. (2011) found a significant variation of 0.6% over the full range of their models and still a 0.2% variation when excluding models with M_c larger than 10 M_⊕. Juno's measurement of Jupiter's λ is expected to have an accuracy of 0.2% (Helled et al. 2011), hence sufficient to distinguish between the different predictions.

Jovian seismology has long been recognized as a unique opportunity to infer interior structure properties (Vorontsov et al. 1976). In particular, the predicted low-frequency acoustic free oscillations depend sensitively on the core mass and internal layer boundaries of theoretical Jupiter models (Gudkova & Zharkov 1999a). However, the expected small amplitudes (∼10 cm) and Jupiter's rapid rotation make it difficult to measure acoustic modes. Recently, Gaulme et al. (2011) used the SYMPA spectrometer to detect the lowest frequency modes (∼1000 μHz) of Jupiter from spatially resolved radial velocity measurements. In the approximation of low degrees l and overtones of high radial order n ≫ l, the acoustic modes ν_{n, l} are proportional to the characteristic frequency ν₀ (also called equidistance) and are equally spaced by ν₀ for a given l: ν_{n, l} ≃ (n + 1/2)ν₀. Gaulme et al. (2011) observed ν₀ = 155.3 ± 2.2 μHz. The inverse of ν₀ is twice the time a sound wave needs to travel from the center to the surface (Jupiter's troposphere),

$$\begin{equation} \nu _0=\left[2\int _0^{R_{\rm J}}dr\:c_s^{-1}\right]^{-1}, \end{equation} \tag{ 5 }$$

a useful parameter that interior models can easily be compared to. In Equation (5), c_s is the adiabatic sound velocity. Gudkova & Zharkov (1999a) calculate ν₀ = 152–155 μHz for SCvH-i EOS-based Jupiter models with five layers and core masses of 3–10 M_⊕. Gaulme et al. (2011) obtain the same range of values for ν₀ with three-layer models of Jupiter that have a small internal He discontinuity and core masses of 0–6 M_⊕. For our three-layer models with M_c = 0–8 M_⊕, we find ν₀ = 155.7–156.3 μHz. Since this range is rather narrow and slightly above that of former Jupiter models, a more accurate determination of the global mode spacing may help to rule out some Jupiter models but not constrain Jupiter's core mass further within the current uncertainty of ∼0–10 M_⊕.

For a better comparison with published calculations of the cooling of Jupiter, we first assume a constant irradiation flux and constant angular velocity over time. The latter assumption is justified by the rather small change of the rotational energy during evolution compared to Jupiter's intrinsic energy loss (Hubbard 1977) and the small change in planet radius over most of the contraction timescale (Guillot et al. 1995). Our H-REOS.2-based Jupiter model J11-4a has a calculated cooling time of 4.66 ± 0.04 Gyr, see Figure 5, and model J11-8a of 4.68 ± 0.04 Gyr, where the uncertainty mostly arises from the observational error bar of T_eff. This result can be considered in good agreement with the age of the solar system (τ_☉ = 4.56 Gyr) given that none of the alternative proposed Jupiter models that rely on free-energy models for the EOS (Saumon & Guillot 2004), such as LM-H4 (4.0 Gyr), SCvH-i (4.7 Gyr), and LM-SOCP (4.8 Gyr), gives a better agreement. On the other hand, Fortney et al. (2011) investigated the uncertainty arising from the use of the Graboske et al. (1975) model atmosphere grid for Jupiter's evolution. If instead a self-consistent model atmosphere grid was applied as is the common approach for exoplanets, Jupiter's cooling time increases by ∼0.5 Gyr. Sinking He droplets from H–He phase separation in Jupiter would also prolong the cooling time further. A process working in the opposite direction could be rising material from core erosion. During the first 500 Myr, Jupiter has shrunk from 1.4–1.5 R_J down to 1.1 R_J (Figure 5). Because the initial conditions for the long-term evolution of a 1 M_J planet are forgotten after 0.01 Gyr (Marley et al. 2007), we consider the displayed radius evolution at young ages realistic.

Zoom In Zoom Out Reset image size

In a second step, we keep the angular momentum conserved to Jupiter's current value (within a numerical accuracy of 0.4% which gives sufficiently smooth evolution curves) and include the subsequent change of rotational energy dE_rot(t), as also done in Hubbard (1970). Such second-order effects are important for estimates of the size of the correction factors that are necessary to bring the cooling time in agreement with τ_☉. It is clear that both effects (λ and dE_rot) will prolong the cooling time: (1) for the same total mass, a larger planet (the young Jupiter) rotates slower, implying a weaker centrifugal force that pushes matter outward, and the consequently smaller radius will allow less energy to be radiated away from the surface and (2) angular momentum conservation (L = Iω) then implies a lower energy of rotation E_rot = 1/2 Lω at young ages. The increase in E_rot with time must be compensated for by a reduced luminosity, implying again a longer cooling time. For model J11-4a, the first effect increases τ by 0.1 Gyr, and the second one by 0.2 Gyr, so that we find τ = 4.96 Gyr. In a third step, we consider a time-dependent solar irradiation by a linear approximation of the Sun's luminosity evolution, which is assumed to have started with 70% of the current value. We find the lower irradiation to speed up the cooling by ∼0.6 Gyr, so that we end up with τ = 4.41 ± 0.04 Gyr. Figure 6 shows the relations between the structure parameter T₁ that determines the planet radius, P = 2π/ω and λ. Their evolution with time can be read from the cooling curve T₁(t).

Zoom In Zoom Out Reset image size

The resulting cooling times of 4.4 Gyr (linearly increasing insolation) to 5.0 Gyr (constant insolation) suggest a reasonable understanding of Jupiter's interior and evolution, where remaining uncertainties seem to be attributable to uncertainties in the luminosity evolution of the Sun. This may be counterintuitive as the present Sun is, by means of accurate helioseismology, far better constrained than Jupiter. For instance, uncertainties in the sound velocity profile of the solar standard model are of the order of 0.1% only (Boothroyd & Sackmann 2003). On the other hand, non-standard solar evolution models that predict a bright and more massive young Sun (Sackmann & Boothroyd 2003) cannot completely be ruled out by observationally derived stellar mass-loss rates; see Güdel (2007) for a review.

However, homogeneous, adiabatic evolution models for Saturn and Uranus that are equivalent to those presented here for Jupiter keep failing to reproduce the observed luminosities (Fortney et al. 2011), indicating that the standard three-layer model assumption may be too simplistic for some giant planets. Therefore, Leconte & Chabrier (2012) investigated the possibility of layered convection in Jupiter, where heat is transported inefficiently by diffusion across thin stable layers, separated vertically by convective cells. They determined the maximum superadiabaticity that would give density distributions in agreement with the gravity field data. The deduced presence of roughly 10⁴ diffusive interfaces in Jupiter will qualitatively reduce the heat flux out of the deep interior and shorten the cooling time (Stevenson 1985). Quantitative estimates of (at least) this effect are necessary before one can conclude a reasonable understanding of Jupiter's evolution.