  Hydrogen and helium demix when sufficiently cool, and this bears on the
evolution of all giant planets at large separations at or below roughly a
Jupiter mass. We model the thermal evolution of Jupiter, including its evolving
helium distribution following results of ab initio simulations for helium
immiscibility in metallic hydrogen. After 4 Gyr of homogeneous evolution,
differentiation establishes a thin helium gradient below 1 Mbar that
dynamically stabilizes the fluid to convection. The region undergoes overstable
double-diffusive convection (ODDC), whose weak heat transport maintains a
superadiabatic temperature gradient. With a generic parameterization for the
ODDC efficiency, the models can reconcile Jupiter's intrinsic flux, atmospheric
helium content, and radius at the age of the solar system if the Lorenzen et
al. H-He phase diagram is translated to lower temperatures. We cast the
evolutionary models in an MCMC framework to explore tens of thousands of
evolutionary sequences, retrieving probability distributions for the total
heavy element mass, the superadiabaticity of the temperature gradient due to
ODDC, and the phase diagram perturbation. The adopted SCvH-I equation of state
favors inefficient ODDC such that a thermal boundary layer is formed, allowing
the molecular envelope to cool rapidly while the deeper interior actually heats
up over time. If the overall cooling time is modulated with an additional free
parameter to imitate the effect of a colder or warmer EOS, the models favor
those that are colder than SCvH-I. In this case the superadiabaticity is modest
and a warming or cooling deep interior are equally likely.
