Synthetic Evolution Tracks of Giant Planets

Giant planet evolution models play a crucial role in interpreting observations and constraining formation pathways. However, the simulations can be slow or prohibitively difficult. To address this issue, we calculate a large suite of giant planet evolution models using a state-of-the-art planetary evolution code. Using these data, we create the python program planetsynth that generates synthetic cooling tracks by interpolation. Given the planetary mass, bulk&atmospheric metallicity, and incident stellar irradiation, the program calculates how the planetary radius, luminosity, effective temperature, and surface gravity evolve with time. We demonstrate the capabilities of our models by inferring time-dependent mass-radius diagrams, estimating the metallicities from mass-radius measurements, and by showing how atmospheric measurements can further constrain the planetary bulk composition. We also estimate the mass and metallicity of the young giant planet 51 Eri b from its observed luminosity. Synthetic evolution tracks have many applications, and we suggest that they are valuable for both theoretical and observational investigations into the nature of giant planets.


INTRODUCTION
Studying the internal structure and bulk compositions of gas giants is crucial for our understanding of giant planet formation and evolution. However, because giant planets are mostly composed of compressible hydrogen and helium, they contract and their observable properties and interiors change as the planets evolve (Hubbard 1977). Accounting for the time-dependent nature properly requires comprehensive thermal evolution models. In order to better characterize giant planets, their the basic physical parameters such as their mass, bulk & atmospheric composition, and stellar irradiation must be known (Burrows et al. 2007; Thorngren et al. 2016). While observations provide valuable information -such as occurrence rates, orbital distances, masses and radii -combining the measurements with theory is critical for the data interpretation and for taking full advantage of the measurements.
One key objective of giant planet science is to determine the planetary bulk composition from mass-radius measurements Miller & Fortney 2011;Thorngren et al. 2016;Teske et al. 2019). This is important because the predicted composition depends on how planets form (Mousis et al. 2009;Johansen & Lambrechts 2017;Hasegawa et al. 2018;Ginzburg & Chiang 2020;Shibata et al. 2020) and can be used to narrow the range of possible formation pathways. An additional constraint is the planetary atmospheric composition that will be measured with high precision with the upcoming James Webb Space Telescope (JWST; Gardner et al. (2006)) and ESA's Ariel mission (Tinetti et al. 2018) , which can be used to further constrain the planetary interior and to provide ★ E-mail: simon.mueller7@uzh.ch information on the link between bulk and atmospheric composition ). More precise radius and mass measurements, as well as an accurate determination of stellar ages are expected from the ESA Plato mission (Rauer et al. 2014). Improved measurements are also expected from ground-based facilities such as HARPS (Mayor et al. 2003), NIRPS (Bouchy et al. 2017) and ESPRESSO (Pepe et al. 2021).
The upcoming data are expected to significantly improve our understanding of gaseous planets. However, the interpretation of the observations crucially depends on predictions from evolution models, which come with their own uncertainties: Previous studies have shown that model parameters such as equations of state and atmospheric models (Burrows et al. 2007; Baraffe et al. 2008;Vazan et al. 2013; Thorngren et al. 2016;Poser et al. 2019) influence the model predictions to a degree that is important as measurements become more accurate. Therefore, unyielding the full potential of these new observations requires careful modeling of the evolution of giant planets (Müller et al. 2020b).
Another increasingly important application for evolution models is the mass determination of directly imaged exoplanets by, for example, the Hubble Space Telescope or SPHERE (Beuzit et al. 2019). In the near future, JWST is expected to provide many more detections. Direct imaging is unique in the sense that it measures the planetary luminosity, which then requires evolution models to infer the planet's mass (e.g., Baraffe et al. 2003;Marley et al. 2007; Mordasini et al. 2017). For giant exoplanets, inferring the mass from luminosity measurements is commonly done by using pre-computed tables of planetary isochrones (Baraffe et al. 2003). The disadvantage of this method is that it over-simplifies the planetary evolution suggesting that it is mostly the planetary mass that controls the luminosity while in fact there are several other important parameters that affect the planetary evolution (Müller et al. 2020b).
Considering the wealth of new observations expected in the upcoming years, access to comprehensive evolution models is essential. With that in mind, we have calculated an extensive set of giant planet evolution models using a modified version of Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. (2011Paxton et al. ( , 2013Paxton et al. ( , 2015Paxton et al. ( , 2018Paxton et al. ( , 2019). These evolution models were then used to build planetsynth (https://github.com/ tiny-hippo/planetsynth), a publicly available python program that generates synthetic evolution tracks that account for the planetary mass, bulk & atmospheric metallicity and incident stellar irradiation. Given these parameters, planetsynth calculates how the planet's radius, luminosity, effective temperature and surface gravity evolve with time.
This paper is organised as follows: in §2, we present the planetary evolution calculations and describe the most important model assumptions and their influences. Then, we present the parameter space covered by our suite of evolution models in §2.1. We explain how the synthetic cooling tracks are generated in §2.2, and validate our approach in §3.1. We then propose key possible implications for the synthetic cooling tracks, and demonstrate their use on some real-world examples in §3.2. Key assumptions and limitations of our models as well as future work are discussed in §4. A summary of our work is given in §5.

THERMAL EVOLUTION OF GIANT PLANETS
Giant planets start their lives hot and inflated (Hubbard 1977;Cumming et al. 2018;Valletta & Helled 2020). Since their interiors are compressible, they contract as they cool down. To calculate the cooling of giant planets, we use a modified version of the stellar evolution code MESA. It solves the 1D hydrostatic equilibrium evolution equations (e.g., Kippenhahn et al. 2012) and calculates the evolution of the interior structure and of observable quantities such as the planetary radius, luminosity, and effective temperature.
The four parameters we use to create the planetary models are mass , bulk metallicity , envelope metallicity and incident stellar irradiation * (see §A1 for details on how the models were created with MESA). The bulk composition is defined by the hydrogen, helium, and heavy-element mass fractions. We assume a proto-solar hydrogen-helium ratio ( = 0.705, = 0.275), which defines as the free parameter that sets the planetary bulk composition.
For planets smaller than 5 and with ≠ , we further assume a core-envelope structure, i.e., the heavy-elements are mostly contained in a compact core in the deep interior. It should be noted, however, that formation models clearly show that the deep interiors of giant planets consist of composition gradients or more extended and dilute cores (Lozovsky et al. 2017;Helled & Stevenson 2017;Valletta & Helled 2020). Nevertheless, this region is typically relatively small. Furthermore, evolution models suggest that composition gradients are quickly destroyed by large-scale convection (Müller et al. 2020a) unless the internal temperatures are relatively low (Vazan et al. 2018).
For planets beyond 5 or with = a homogeneous composition is used. The heavy-element distribution directly affects the radius and luminosity of the planet (Baraffe et al. 2008;Vazan et al. 2013). The effect is, however, small compared to the choice of the equation of state (hereafter EoS) or the material that is used to represent the heavy elements (Müller et al. 2020b), in particular for massive giant planets. Table 1. Masses , bulk metallicities , atmospheric metallicities and incident stellar fluxes * covered by the suite of models. In total, the simulation grid is composed of about 45,000 planets.
[M ] * [erg/s/cm −2 ] 0.1 − 1 ≤ 0.80 ≤ 0.10 10 − 10 9 1 − 3 ≤ 0.50 ≤ 0.10 10 − 10 9 3 − 5 ≤ 0.20 ≤ 0.10 10 − 10 9 5 − 30 ≤ 0.04 ≤ 0.04 10 − 10 9 In order to self-consistently account for the heavy elements in the evolution calculations, we made substantial changes to the EoS of MESA as presented in Müller et al. (2020a). The modifications to the EoS make it suitable to model giant planets consisting of arbitrary mixtures of hydrogen, helium and a heavy element. For this study, we have added a 50-50 water-rock mixture, which we use in all our models to represent the heavy elements. More details on the EoS are given in §A2. Other than the EoS, both stellar irradiation and atmospheric metallicities both have a large influence on the evolution (e.g., Burrows et al. 2007;Müller et al. 2020b); we include both of these as model parameters (see §A3 for details).
Our calculations assume that giant planets form by core accretion (Mizuno 1980;Pollack et al. 1996) in the hot-start framework. In hotstart models the planets cool adiabatically from an initially large, but essentially arbitrary radius (e.g., Marley et al. 2007). As described above, we model the planets with either homogeneous compositions or core-envelope structures. In that case, it is reasonable to assume that the planets cool adiabatically for most of their lives, and different initial conditions converge quickly. It must be noted, however, that current studies of Jupiter and Saturn suggest that their interiors areat least in part -non-adiabatic (Debras & Chabrier 2019;Mankovich & Fuller 2021), and we hope to address this in future research. We discuss the implications of using adiabatic hot-start models in more detail in §4 and §A4.

Suite of Giant Planet Evolution Models
Giant exoplanets have large ranges of masses and compositions and are irradiated by their host stars with fluxes that can range many order of magnitudes (Thorngren et al. 2016;Teske et al. 2019), and the suite of planetary evolution models presented here needs to account for this diversity. The evolution parameters are sampled from a non-uniform four-dimensional grid spanned by the planetary mass , bulk heavyelement fraction , atmospheric metallicity and incident stellar flux * . The planetary masses we consider range from 0.1 − 30 . For planets more massive than about a Jupiter mass, the planetary radius is only changing slightly as a function of mass. Therefore, we sample the masses more coarsely beyond 1 . Note that the mass limit for the ignition of deuterium is generally believed to be at about 13 (e.g., Boss et al. 2007). However, since this would mostly affect very young and massive planets (Spiegel et al. 2011) for which uncertainties in the formation luminosity are high (e.g., Bowler 2016), we do not include deuterium burning in this work.
Depending on the mass, the grid covers different ranges of because i) the bulk metallicity tends to decrease with planetary mass and ii) there are limitations in the EoS. Nonetheless, this provides coverage for most of the inferred compositions of observed giant exoplanets population (Thorngren et al. 2016).
For the stellar flux, we consider the range from 10 − 10 9 erg s −1 cm −2 . This includes cool to warm Jupiters, but excludes extremely irradiated hot Jupiters, since their inflation mechanism remains un-known (e.g., Fortney & Nettelmann 2010;Weiss et al. 2013;Baraffe et al. 2014). The conventional cut-off irradiation for hot Jupiters is estimated to be 2 × 10 8 erg s −1 cm −2 (Miller & Fortney 2011). At the maximum flux that we consider, observations suggest that exoplanets can inflate to 1.5 (Thorngren & Fortney 2018). Therefore our evolution tracks with stellar fluxes beyond 2 × 10 8 erg s −1 cm −2 should be treated with caution. We do not include photo-evaporation, which is expected to become important for the lowest masses and highest stellar fluxes we consider (e.g., Jin et al. 2014). All the parameters covered by our suite of models are summarised in Table 1. In total, our simulation grid consists of about 45,000 planets.
The thermal evolution of the planets in the sample is calculated as described in §2 for 10 10 years. The relevant outputs are the planetary radius , luminosity , effective temperature and surface gravity as they evolve with time. The luminosity corresponds to the heat from the interior without any re-radiated stellar irradiation; the irradiation heats the outer layers of the planet and increases the luminosity at the surface (see §4 and §A3 for details). As discussed above, at early times the evolution could be significantly affected by the initial condition, and hence we exclude data for times below 10 7 years.
For a few parameters the evolution calculations crashed before reaching the end, because some temperatures and densities inside the planet are outside of the EoS boundaries. This occurs mostly for planets at the edges of our simulation grid towards the end of their cooling. However, after about 1 Gyr the planetary radius and luminosity change only very slowly with time. To be on the safe side, we discarded models that crashed before 2 Gyr. If the calculations fail between 2-10 Gyr, we extrapolate the subsequent luminosity and radius evolution with a function of the form ( ) = log( ) + , where and are parameters that are individually fitted to the planetary parameters using its evolution history. The effective temperature and surface gravity are then calculated from the extrapolated luminosity and radius. This approach leads to errors typically smaller than 0.1% in all output quantities, which is negligible compared to both observational and theoretical uncertainties.
To summarise: we simulate a suite of planetary models covering a large range of masses, bulk & envelope metallicities and incident stellar fluxes. The calculations determined the radius, luminosity, effective temperature and gravitational acceleration at the surface ( ( ), ( ), ( ), ( )) as a function of the planetary age .

Generating Synthetic Cooling Tracks
The suite of giant planet evolution models enables the generation of synthetic cooling tracks for a given set of planetary parameters (within the appropriate limits, see Table 1). We have thoroughly tested several approaches, such as directly interpolating on the irregular grid or using machine-learning regression with XGBoost (Chen & Guestrin 2016) and neural networks. The former is too slow when making a large number of predictions. The latter are fast and generally yield accurate outputs. However, ensuring that the predictions were physically consistent is difficult.
We have therefore settled on the following approach: (i) In order to generate data of equal length, a cubic spline interpolation is used to calculate the outputs (radius, luminosity and effective temperature) for each planet at discrete times between 10 7 − 10 10 years. The interpolated data from all planets is then combined into a single data set.
(ii) While it would have been possible to directly interpolate this data to generate new cooling tracks, it is quite slow for a four dimensional parameter space on an irregular grid. Therefore, we first map our data onto a four-dimensional regular, unevenly spaced grid using piece-wise linear interpolation. For the regions in parameter space where we do not have data (e.g., high-mass and metallicity planets), we use nearest neighbour interpolation. While this may be a decent estimate for parameters that are just outside the range covered by our suite of models, the errors will quickly grow large. Therefore, these data points should never be used to calculate predictions, but rather they were a convenient way to create a regular grid. The default behaviour when attempting to generate synthetic cooling tracks for parameters not covered by our model suite is to warn the user and return NaNs.
(iii) The synthetic cooling tracks are then generated by linearly interpolating on the regular grid. For each set of input parameters ( , , , log * ), the interpolation yields ( [ ], log [ ]) at discrete times between log = 7 − 10 [yrs]. In order to provide self-consistent results, we calculate the effective temperature from the generated radius and luminosity, assuming a black body: 4 = /4 2 where is the Stefan-Boltzmann constant. The gravitational acceleration is calculated from the generated radius as log = log G / 2 [cm/s 2 ].
(iv) If the prediction is required at specific time log [yr], a cubic spline interpolation is calculated from the synthetic cooling track.
Our approach is similar to that of the MESA Isochrones and Stellar Tracks (MIST) (Choi et al. 2016;Dotter 2016) project, which uses a grid of stellar evolutionary tracks and isochrones to generate cooling tracks for stars. The last two steps described above were implemented as the open-source python module called planetsynth, with the goal of it being both easy to use and update. Another advantage is that it is very fast: a million synthetic cooling tracks can be generated in a just a few seconds (for vectorized inputs).

planetsynth: A Fast and and Easy-to-Use Python Module for Generating Synthetic Cooling Tracks
The main class of the module is PlanetSynth, which contains the two most important methods synthesize and predict.
• synthesize: takes the planetary parameters ( , , , log * ) as arguments and calculates the cooling track at discrete times between 10 7 − 10 10 years.
• predict: uses single or several points in time to calculate ( , log , , log ) at those specific times for a given set of planetary parameters.
The module is available at https://github.com/tiny-hippo/ planetsynth and includes examples in some detail on how to use it.

RESULTS
In this section we first show that the synthetic cooling tracks yield accurate results when compared to the direct calculations with MESA. Then, we demonstrate the capabilities of our code by calculating time-dependent mass-radius diagrams, inferring metallicities from mass-radius measurements, and inferring the mass and metallicity of a directly imaged planet.

Validation
In order to validate our approach, we create a sample of 200 planets (see §B). The planetary parameters are generated with the Latin Hypercube sampling method (McKay et al. 1797), which is a statistical method commonly used in high-dimensional numerical experiments. Since the validation sample was not included when fitting the interpolants, it can be used as a non-biased sample to validate the performance of the synthetic models.
For each planet in the validation sample, we calculate the thermal evolution and then compare the output to the predictions from the synthetic cooling tracks. For a qualitative comparison, we plot the radius and luminosity evolution of 50 randomly chosen planets from the validation sample in the top panels of Fig. 1. It can be seen that the synthetic cooling tracks agree extremely well with those from the evolution calculations. In the bottom panels of Fig. 1, we plot the radii and luminosities from the evolution models against the synthesised predictions for the entire validation sample: The agreement is excellent, with only small deviations from the actual values. For a more detailed discussion of the model performance and the error distributions see §B.

Time-Dependent Mass-Radius Diagrams
Mass-Radius diagrams are often used to interpret the possible compositions of exoplanets, and to qualitatively compare exoplanets to our solar system planets. While the radii of planets made of refractory materials are not expected to change much with time, this is very different for gas-rich planets that contract significantly with time. The planetary age cannot be ignored when constructing mass-radius diagrams of giant planets since the radius depends on both mass and time.
Using our synthetic evolution tracks it is simple to calculate massradius diagrams that account for the time dependency of the radius. We show two examples in Fig. 2 with = 0.05 (left) and = 0.20 (right). In both cases, the incident stellar irradiation was * = 10 5 erg/s/cm 2 and the atmospheric metallicity was assumed to be solar. Plotting the figures side-by-side makes it clear that planets with large differences in composition can have the same radius at a different time. This will be particularly important when interpreting measurements from Plato, since it will provide accurate stellar ages with an uncertainty of ∼ 10% (Rauer et al. 2014).

Inferring the Bulk Metallicity based on Mass-Radius Measurements
Constraining the bulk compositions of giant exoplanets is critical for understanding giant planet formation models. While the bulk composition cannot be measured directly, it can be inferred from mass-radius measurements (from radial velocity and transit observations) with the help of evolution models (Guillot et al. 2006;Miller & Fortney 2011). Here, instead of using evolution models directly, we use the generated evolution tracks. We follow the approach described in Thorngren et al. (2016): for a given mass, age, stellar flux, and atmospheric metallicity we search for the bulk metallicities that fit the observed radius. We demonstrate this procedure on two giant exoplanets: i) Kepler-16 b with = 0.33 ± 0.02 , = 0.754 ± 0.003 , * = 4.19 × 10 5 erg/s/cm 2 , with an age between 0.5 − 10 Gyr (Doyle et al. 2011), and ii) HAT-P-54b with = 0.76 ± 0.03 and = 0.94 ± 0.03 , * = 1.04 × 10 8 erg/s/cm 2 and an age of 1.8 − 8.2 Gyr (Bakos et al. 2015).
In order to infer the bulk metallicity, samples are drawn from the probability distributions for the measured planetary mass, radius and age. We assume normal distributions for the mass and radius, and a uniform distribution for the age. For each draw, we use a root-finding algorithm and our synthetic cooling tracks to find the metallicity that reproduces the observed radius. This is done twice for each planet with different atmospheric metallicities: one or five time solar. We repeat the procedure 100,000 times, which yields the estimates and their associated uncertainties of the heavy-element contents. The Gaussian kernel density estimates (KDEs) based on these samples are shown in Fig. 3. Note that we limit the x-axis on both plots to exclude low-probability compositions for better visual clarity.
For Kepler-16b the mean inferred metallicity is = 0.34 (SD = 0.03) and = 0.37 (SD = 0.03) for the solar and 5x solar atmospheric compositions, respectively. For HAT-P-54b, the values are = 0.14 (SD = 0.08) and = 0.18 (SD = 0.05). The mean and standard deviations are calculated directly from the samples rather than from the KDEs. While the results for both atmospheric models are within their respective uncertainties, there is a clear trend towards higher heavy-element contents when the atmosphere is en-  riched in heavy elements. This is because the atmosphere becomes more opaque as heavy elements are added, and therefore the planet is hotter at a given age and can be more metal-rich (Burrows et al. 2007;Müller et al. 2020b). Note that our inferred metallicities for the solar atmosphere are significantly lower than previously reported estimates: = 0.39 +0.01 −0.02 for Kepler-16b and = 0.20 ± 0.04 for HAT-P-54b (Thorngren et al. 2016). This can mostly be attributed to the different hydrogen-helium EoSs (Müller et al. 2020b).

Further Constraining the Bulk Composition using Atmospheric Measurements
Both Kepler-16b and HAP-P-54b have uncertainties in mass and/or radius of up to a few percent, while the age is rather poorly constrained. In this case, measuring an atmospheric metallicity would not help to constrain the bulk composition. However, it is interesting to investigate whether this remains true if measurements were to improve. Therefore, we repeat the above calculations, but with mass and radius uncertainties of Δ = 0.005 and Δ = 0.005 , and ages constrained between 2 − 5 Gyr. The resulting distributions are shown in Fig. 4. Now, the inferred metallicities are = 0.34, 0.38 (SD = 0.02, 0.02) for Kepler-16b and = 0.16, 0.19 (SD = 0.02, 0.02) for HAT-P-54b, respectively.
For both planets the peaks in the probability density for the metallicity are now well-separated. This demonstrates that with improved observational uncertainties, measuring the atmospheric metallicity becomes a key ingredient in the interpretation of observations and can significantly improve the characterisation of giant exoplanets.

Inferring the Mass and Bulk Metallicity from Luminosity Measurements
In the examples so far, we have assumed that the planet's radius and mass are known. Alternatively, a giant planet can be observed by detecting its thermal emissions. In this case, the planet's luminosity (and atmospheric composition) are known, and the planetary mass can be inferred by comparing the observed luminosity with theoretical, mass-dependent evolutionary tracks (e.g., Baraffe et al. 2003).
Here, we use synthetic cooling tracks to infer the mass and metallicity of the young giant planet 51 Eri b. Its observed luminosity is log [ ] = −5.8 to −5.4, its distance from its star is ∼ 13 au and its mass was estimated to be ∼ 2 (assuming a hot start) (Macintosh et al. 2015). We follow a similar Monte Carlo approach as described in §3.2.4. We draw luminosities from a normal distribution with a mean value of log [ ] = −5.6 and a standard deviation of Δ log [ ] = 0.2. The age is inferred from a uniform distribution between 17 − 23 Myr. For the atmospheric metallicity, we use the measured host star [Fe/H], which is roughly solar.
For each sample, we use a multivariate optimization scheme to find the mass and bulk metallicity that closely matches the observed luminosity. We repeat this procedure 10,000 times and use a 2d Gaussian kernel density estimate to infer the probability distribution in the mass-metallicity plane. The resulting two-dimensional distribution is shown in Fig. 5. In Fig. 6, we show the 1d Gaussian kernel density estimate for mass of 51 Eri b.
The mass distribution is a slightly right-skewed Gaussian, with a mean value of = 2.3 (SD = 0.4 ) for the mass and = 0.11 (SD = 0.05) for the bulk metallicity. The mean mass estimate is higher than the 2 reported in Macintosh et al. (2015), which can be explained by different model assumptions, such as equations of state and atmospheric models.
It is interesting to note that the most probable mass 2.3 suggests a rather high bulk metallicity of 0.11, which corresponds to a heavy-element mass of roughly 70 ⊕ . Large heavy-element masses are difficult to explain with current planet formation models since cores of gas giants are currently thought to be limited in size (e.g., Bitsch et al. 2018). Additional super-stellar metallicities could be explained by the accretion of planetesimals (Mousis et al. 2009;Shibata et al. 2020) or pebbles (Johansen & Lambrechts 2017), or collisions between primordial cores (Ginzburg & Chiang 2020). However, it is still unclear whether these processes could explain highly-enriched gas giants. It is also important to restate that our models are based on a hot-start formation scenario. Assuming a cold-start, luminosities at a given mass could be orders of magnitude lower (see §A4), which would yield a mass estimate of 10 (Macintosh et al. 2015).

DISCUSSION
While this study represents an advancement in evolution calculations, several aspects in the models can clearly be improved. To start with, future work could provide a higher resolution for the planetary parameters, which would further reduce the interpolation errors. For the massive planets, our models neglect deuterium burning, which could affect the luminosity and radius of young and massive planets (Spiegel et al. 2011). We hope to include this in a future update of our evolution models. We also plan to include various opacity and EoS tables since they can influence the exoplanetary radius predictions (Guillot 1999). Another aspect is the atmospheric model: while the atmospheric opacity accounts for heavy-element enrichment and the effects of stellar irradiation are included, we do not consider the effect of clouds or grains. Atmospheres that are cloudy or enriched in grains are more opaque, and trap the primordial heat more efficiently. This leads to a delayed cooling, which can significantly affect the evolution of a planet (Vazan et al. 2013;Poser et al. 2019), in particular at young ages. Another aspect that can be improved is the treatment of stellar irradiation: while the method we used in this study is in agreement with more sophisticated treatments for planets that are not extremely irradiated (Valsecchi et al. 2015) it can clearly be improved. Future work should also include semi-grey atmospheric models (e.g., Guillot 2010;Parmentier & Guillot 2014) or full atmospheric models that include irradiation (e.g., Fortney et al. 2007) as well as different planetary albedos.
Due to limitations of the heavy-element EoS we did not model very massive metal-rich planets. While our suite of models covers most currently discovered planets (in terms of mass and metallicity), there a few extremely enriched planets (e.g., Barragán et al. 2018) that cannot currently be modelled with planetsynth. An additional limitation of the EoS is that the heavy-elements are represented by a 50-50 rock-ice mixture. However, giant exoplanets are diverse and the heavy elements can differ between planets. Giant planets form at various locations in the protoplanetary disk, which themselves can differ in elemental ratios, and the formation location determines the composition of the accreted solids due to different condensation temperatures. The exact heavy-element composition can affect the planetary radius, by a few percent (Baraffe et al. 2008).
The evolution calculations in this work correspond to the hot-start framework (Baraffe et al. 2003;Marley et al. 2007). If planets cool adiabatically, different initial conditions in hot-start models converge quickly . However, this is not true for non-adiabatic contraction, for example, due to the presence of large composition gradients (Vazan et al. 2013(Vazan et al. , 2015 or doublediffusive instabilities (Leconte & Chabrier 2012;Wood et al. 2013). While currently it seems that hot-starts are expected in the core accretion framework Cumming et al. 2018), future calculations could include cold Fortney et al. 2008) and warm-start initial configurations (Spiegel & Burrows 2012;Marleau & Cumming 2014).
Recent studies of the solar system gas giants, Jupiter and Saturn, suggest that their interiors are not fully convective/adiabatic (e.g., Debras & Chabrier 2019;Mankovich & Fuller 2021). Giant planets that are not fully convective cool down more slowly, because radiative/conductive or double-diffusive energy transport is less efficient. As a result, they would have lower luminosities early in their lifetimes because of the less efficient cooling. At some point, however, their luminosities would be higher, because their interiors stay hot for longer. Additional models should therefore ideally also consider non-adiabatic interiors.

CONCLUSIONS
In this work we calculate a comprehensive suite of giant planet evolution models with MESA for a large variety of planetary masses, bulk & atmospheric metallicities and incident stellar fluxes. We present the freely available python module planetsynth (https: //github.com/tiny-hippo/planetsynth) that interpolates on the model grid to generate synthetic evolution tracks. These cooling tracks converge well to those calculated by MESA, with accuracies generally well beyond observational and theoretical uncertainties.
These synthetic cooling tracks have many scientific applications, including: • Time-dependent mass-radius can easily be generated, which can aid the characterisation of giant planets.
• By using mass-radius measurements of exoplanets, cooling tracks can be used to infer a planet's bulk composition.
• We show that if observational uncertainties on mass, radius and age diminish further, measurements of the atmospheric composition can be used to constrain the bulk composition.
• We demonstrate how the the bulk metallicity and mass of the young directly imaged planet 51 Eri b can be inferred from synthetic cooling tracks.
Many accurate and new measurements of giant exoplanets are expected in the upcoming years. For example, JWST will detect young planets in the mid-infrared, which will allow it to image planets with smaller masses than currently possible. Another great advancement will come with the observation of many giant planet atmospheresin particular their compositions -with the Ariel mission. Combined with the improved determination of stellar ages from Plato and improved measurements from ground-based facilities such as HARPS, NIRPS and ESPRESSO, this will further constrain the compositions of giant exoplanets, opening up compelling avenues to improve our understanding of how giant planets form and evolve. To make of the most of these exciting opportunities, it is essential to consider giant planets as evolving objects. Giant planet evolution models will therefore play a crucial role, and it will be important that these models continue to be refined. The framework for creating synthetic evolution tracks presented in this work is flexible, and it will be possible to include a new generation of giant planet evolution models that account for non-adiabatic cooling histories, composition gradients and other physical mechanisms that can affect the thermal evolution of giant exoplanets.

A1 Creating Planetary Evolution Models with MESA
The evolution models are calculated in three steps: (i) First, we create initial models for the range of planetary masses that we consider using the create_initial_model namelist option. We use similar initial radii for all masses, except that we slightly scale them upwards with increasing mass. In this step, we use the same composition for all the planets, except when = 0. The models are then evolved briefly for a few thousands years for numerical stability.
(ii) As the next step, we use the relax_initial_composition namelist option to create the core-envelope profile for the assigned bulk and atmospheric metallicity. This sets the size of the core as well as the metallicity in the envelope and atmosphere. This step was skipped if = 0. Again, the models are evolved for 10 3 years.
(iii) Lastly, the relaxed models are evolved for 10 10 years with the stellar irradiation added via the column_depth_for_irradiation and irradiation_flux options. Since we consider stable coreenvelope or homogeneous structures, mixing of elements is disabled by setting mix_factor = 0, and convective regions are determined with the Schwarzschild criterion (use_Ledoux_criterion = .false.).

A2 Equation of State
Depending on the planetary mass, we use two different versions of MESA. For lower-mass planets ( ≤ 5 ), we use release 10108 with a modified EoS suitable for a wide range of planetary conditions and compositions (Müller et al. 2020a). Our EoS uses the recent Chabrier et al. (2019) hydrogen-helium EoS (hereafter CMS) and combines it with QEoS (More et al. 1988;Vazan et al. 2013) for rock (SiO 2 ) or water to create ideal mixture of hydrogen, helium and a heavy element. For this work, we implement an additional EoS that uses a 50-50 rock-water mixture, calculated as an ideal mixture. Our heavy-element EoS is valid for densities up to = 100 [g/cc], which is not enough for massive giant planets. Therefore, for highmass planets ( > 5 ) we use the MESA version 15140 with the built-in MESA EoS, which has the option to use the CMS EoS for H-He. For the conditions relevant in this work, the MESA EoS can be used for metallicities up to = 0.04, which sets the upper limit of metallicities that we consider for these massive planets.

A3 Atmosphere
The metallicity in the outer envelope can have a large effect on the evolution via the opacity. This effect is strongest for highly irradiated planets: higher atmospheric metallicities slow the cooling, and yield larger radii and luminosities at a given time (Burrows et al. 2007;Müller et al. 2020b). Therefore, the envelope metallicity is an important free parameter, which ranges between = 0 up to min(0.1, ) in this work. For the atmospheric model, we use the low-temperature opacities of Freedman et al. (2014) with the simple_photosphere MESA boundary condition. This defines the photosphere at an optical depth of = 2/3.
The irradiated surfaces of the planets are treated with the * − Σ * method: the stellar irradiation * is applied to a mass column Σ < Σ * , where Σ( ) = ∫ ( ) . The treatment of the stellar irradiation used in these models can be interpreted in the following way: Irradiation penetrates into the outer layers of the planet, which increases the transported luminosity at the surface. Excluding highly irradiated hot Jupiters, this method provides good agreement with semi-analytical models for irradiated planetary atmospheres (Guillot 2010;Guillot & Havel 2011). We use Σ * = 3×10 2 g cm −2 in order to match these models, however the planetary evolution is not sensitive to this choice (Valsecchi et al. 2015).

A4 Initial Conditions
Calculating the evolution of planets requires suitable initial conditions. We assume that the planets form by core accretion (Mizuno 1980;Pollack et al. 1996). In this framework a key assumption is the initial primordial thermal state, which defines the subsequent evolutionary path. Here, we consider so-called hot-start models (e.g., Baraffe et al. 2003;Marley et al. 2007). In a hot-start, the initially adiabatic planet is created with an arbitrarily large radius. Alternative core accretion formation models are cold- Fortney et al. 2008) or warm-start models (Spiegel & Burrows 2012;Marleau & Cumming 2014). For very young planets, predictions from the different models can vary greatly (up to several orders of magnitude in luminosity) (see e.g., Bowler 2016). Therefore, we focus on planets older than 10 7 years in this work.

APPENDIX B: VALIDATION SAMPLE & ERROR DISTRIBUTIONS
In §3.1 we demonstrate that the synthetic cooling tracks are a good approximation of the actual evolution calculations. In Fig. B1 we show the validation sample in the log − space; the size and colour of the scatter dots depicting their envelope metallicity and incident stellar irradiation, respectively.
We next present the distributions of mean absolute percentage errors (MAPE) for the radius and luminosity in Fig. B2. The distributions are approximately Gaussian, slightly skewed to the right. The cumulative distributions show that for almost all the validation samples the MAPE is less than 0.01%. As visible in Fig. 1, there are some rare cases where the errors are significant, since they start to become comparable to observational uncertainties. This may occur if the planetary parameters are at or near to the edge of the simulation grid (see Table 1). When interpreting the errors, it is important to keep in mind that there is always an inherent uncertainty in the evolution models due to the assumptions that have to be made. This paper has been typeset from a T E X/L A T E X file prepared by the author.