No evidence for radius inflation in hot Jupiters from vertical advection of heat

Elucidating the radiative-dynamical coupling between the upper photosphere and deeper atmosphere may be key to our understanding of the abnormally large radii of hot Jupiters. Very long integration times of 3D general circulation models (GCMs) with self-consistent radiative transfer are needed to obtain a more comprehensive picture of the feedback processes between dynamics and radiation. Here, we present the longest 3D nongray GCM study to date (86000d) of an ultra-hot Jupiter (WASP-76 b) that has reached a ﬁnal converged state. Furthermore, we present a method that can be used to accelerate the path toward temperature convergence in the deep atmospheric layers. We ﬁnd that the ﬁnal converged temperature proﬁle is cold in the deep atmospheric layers, lacking any sign of vertical transport of potential temperature by large-scale atmospheric motions. We therefore conclude that coupling between radiation and dynamics alone is not su ﬃ cient to explain the abnormally large radii of inﬂated hot gas giants.


Introduction
Characterizing and understanding the atmospheres of exoplanets is now more within reach than ever before thanks to the recent launch of the James Webb Space Telescope (JWST).Many of the observed hot gas giants exhibit very large radii and therefore low bulk densities.One of the long-standing questions in the field of hot and ultra-hot Jupiters (UHJs) pertains to the mechanism that inflates those planets.Interior models of the observed exoplanet population reveal that the state of inflation links the incident flux with the temperature in the deep atmosphere (Thorngren & Fortney 2018;Thorngren et al. 2019;Sarkis et al. 2021).These studies find that the intrinsic temperatures increase with instellation and peak at an instellation of ≈2 × 10 9 erg −1 s −1 cm −2 (Thorngren et al. 2019).Different explanations have been proposed to understand the mechanisms that lead to these hot intrinsic temperatures.There are two hypotheses that could explain inflation (Fortney et al. 2021).The first proposes that planets are heated by tidal interactions (e.g., Bodenheimer et al. 2001;Arras & Socrates 2010;Socrates 2013), whereas the second claims that incident stellar flux is deposited into the deep atmosphere.
A prominent theory that fits with the latter hypothesis is Ohmic dissipation (Perna et al. 2010;Batygin & Stevenson 2010;Menou 2012;Rauscher & Menou 2013;Helling et al. 2021;Knierim et al. 2022), which proposes that magnetic fields should couple to ionized flows in the upper atmosphere, slowing down zonal winds and therefore depositing heat by fric-tion.However, it is unclear whether or not this mechanism is sufficiently active at sufficiently deep layers to heat the deep atmosphere (Rauscher & Menou 2013).Using a 2D model with parameterized vertical transport, Tremblin et al. (2017) find that heat would be transported downwards by means of vertical transport of potential temperature without the need of magnetic fields.The authors claim that incident energy input from stellar instellation is self-consistently advected downwards.Idealized 3D general circulation models (GCM) with parametrized thermal forcing seem to reproduce this mechanism (Sainsbury-Martinez et al. 2019, 2021).
In recent work, we introduced expeRT/MITgcm, a 3D GCM with self-consistent radiative transfer for hot Jupiters without the limitations of the above-mentioned models (Schneider et al. 2022).In this latter Letter, we focussed on the effects of surface gravity and rotation rate on the temperature evolution.Here, we extend this work to the final state of the converged atmosphere, using the same 3D nongray GCM and integrating for long enough to converge the temperature in the deep atmosphere.
In order to find a planet that is suitable for studying the issue of temperature convergence in the deep layers of the atmosphere, we decided to go to the temperature extremes and simulate the ultra-hot Jupiter WASP-76 b.WASP-76 b has been observed with radial velocity (West et al. 2016) and transit spectroscopy (Seidel et al. 2019) and has both a Spitzer phase curve (May et al. 2021) and winds characterized by high-resolution spectroscopy (Ehrenreich et al. 2020;Kesseli & Snellen 2021;Kesseli et al. 2022).Additionally, the climate of WASP-76 b was recently modeled (Wardenier et al. 2021;Savel et al. 2022;Beltz et al. 2022).WASP-76 b has a very low bulk density and can therefore be considered to be inflated, while we also know that it receives a high incident stellar flux.We therefore expect there to be a mechanism at play that inflates the radius of this giant planet by transporting potential energy downwards.
We begin by describing the methods used in this Letter (Sect.2).We then explain how we reached temperature convergence in our models (Sect.3) before discussing the limitations of our methods (Sect.4).Finally, we present our conclusions in Sect. 5.

General circulation model
This work utilizes expeRT/MITgcm described in detail in Schneider et al. (2022) and Carone et al. (2020).In particular, expeRT/MITgcm builds on the dynamical core of the general circulation model (GCM) MITgcm (Adcroft et al. 2004) and solves the 3D hydrostatic primitive equations of fluid dynamics (e.g., Showman et al. 2009, Eqs.(1)-( 4)) on an Arakawa C-type cubed-sphere grid with a resolution of C32 assuming an ideal gas.The vertical grid contains 41 logarithmically spaced layers from 1 × 10 −5 bar to 100 bar, extended by six linearly spaced layers between 100 bar and 700 bar.The model includes Rayleigh friction at the bottom of the computational domain (below 490 bar), a sponge layer at the top of the atmosphere (above 1 × 10 −4 bar), and a fourth-order Shapiro filter.We caution that the choice of numerical methods, such as the choice of the dynamical core, often has a significant influence on the circulation and the temperature (Heng et al. 2011;Polichtchouk et al. 2014;Skinner & Cho 2021;Carone et al. 2021).For the scope of this work, we identify the use of deep drag as having possible repercussions for our findings, and therefore discuss the influence of the drag on the results of this work in Appendix A.
The temperature in the atmosphere is forced by radiative heating and cooling using a multi-wavelength radiative transfer scheme that operates during the runtime of the climate model.More specifically, the radiation field is updated with a radiative time step of ∆t rad = 100 s, which is four times the dynamical time step ∆t dyn = 25 s.The radiative transfer includes scattering and uses five correlated-k wavelength bins with 16 Gaussian quadrature points per wavelength bin.The implementation of radiative transfer in expeRT/MITgcm is based on petitRADTRANS (Mollière et al. 2019(Mollière et al. , 2020)).One advantage of using expeRT/MITgcm is its flexibility and performance, which enable long integration times while maintaining the accuracy of a multi-wavelength radiation scheme (Schneider et al. 2022).
As the goal of this Letter is to look at the conservation of energy and the interplay between dynamical energy transport and radiative energy transport, we updated expeRT/MITgcm to conserve energy lost by friction at the bottom boundary.Kinetic energy lost due to friction on the horizontal velocity u [m s −1 ] is converted to thermal energy using the following equation (Rauscher & Menou 2013): where dT dt [K s −1 ] is the change in the temperature T [K], τ deep is the friction timescale (as described in Schneider et al. 2022;Carone et al. 2020), and c p [J kg −1 K −1 ] is the heat capacity at constant pressure.Otherwise, we use the same parameters as in Schneider et al. (2022).

Setup of WASP-76 b simulations
We use expeRT/MITgcm to model WASP-76 b, an ultra-hot Jupiter orbiting WASP-76, where we assume T = 6250 K and R = 1.73 R as the stellar temperature and radius of WASP-76 (West et al. 2016).We model WASP-76 b on a tidally locked orbit of 0.033 AU with an orbital period of 1.8 d (West et al. 2016).The planetary radius and mass are taken as 1.83 R jup and 0.92 M jup respectively (West et al. 2016), resulting in a low gravity of g = 6.81 m s −2 .We use c p = 13784 J kg −1 K −1 and R = 3707 J kg −1 K −1 as the specific heat capacity at constant pressure and the specific gas constant, respectively.
We set up two models in this Letter (Fig. 1).The nominal model is a model of WASP-76 b, which has been initialized with a particularly hot initial state.We use the same method to initialize the temperature profile that was used in Schneider et al. (2022), incorporating an adiabatic temperature profile below 10 bar with Θ ad = 4000 K as the temperature that the adiabat would have at 1 bar.We chose such a hot initial state because planets are thought to form in hot conditions (e.g., Fortney et al. 2021).
The steroids model has been initialized from the nominal model at t = 40 000 d.The initialization of the steroids model is done by extrapolating the temperature evolution of the nominal model using log-linear regression.The extrapolation to find the steroids model resembles a hotter limit on the prediction from the temperature evolution of the nominal model up until t = 40 000 d, while we take over the dynamical state (e.g., the velocity field).We explain the details of the fit and initialization in Appendix B.

Results from nongray 3D GCM
Utilizing the performance of expeRT/MITgcm, we monitor the temperature convergence of WASP-76 b up to t = 86000 d.We show the final atmospheric state and observables in Appendix C. The horizontally averaged thermal profile for pressures above p = 0.1 bar is shown in Fig. 2. The nominal model cools down from an initially hot state to an almost stable solution within 30 000 d, where the temperature changes at the bottom boundary become smaller than 0.1 K d −1 , as shown in Fig. 3.However, the temperature of the model continues to decrease further and will reach a much colder state within reasonable timescales.Structure models such as those of Thorngren et al. (2019) and Sarkis et al. (2021) often give predictions about the location of the radiative convective boundary (RCB), defined as the first pressure layer from the bottom upwards, where the temperature becomes subadiabatic (Thorngren et al. 2019).For better comparison of our model to those predictions, we calculate the RCB in the same way, and show this as gray dots in Fig. 2. We fitted the temperature evolution around 40 000 day in order to estimate how much colder the final state would be.Assuming an exponential decay of the temperature changes with time, we estimated the final temperature at p = 650 bar to be T = 4567 K (Appendix B).We performed a second simulation (dubbed steroids) using this final state as an initial condition to monitor the stability of such a model.The resulting temperature evolution of this model is shown as dotted lines in Fig. 2. The new initial state seems to be a good initial guess in the deep layers, whereas some adjustment happens in intermediate layers between 0.1 bar and 100 bar.Looking at the temperature change rates (Fig. 3), we find that the temperature reaches a steady state, where the temperature stops changing, and the temperature changes start to swap signs in an oscillating manner with ≈5 × 10 −3 K d −1 .We therefore claim that temperature convergence has been reached in the steroids model, implying that the temperature of our WASP-76 b models self-consistently converges to a cold interior.
The planet gains energy from stellar irradiation and loses energy by thermal radiation.We can therefore formulate the energy balance as an integral over the net radiative flux F tot as where ϑ and ϕ are latitude and longitude respectively, T int [K] is the intrinsic temperature and R pla is the planetary radius.Here, F pla and F are the planetary and stellar fluxes, respectively.In radiative equilibrium, without an internal energy source, it should hold that and subsequently T int = 0 K from Eq. ( 2).We note that we do not impose a nonzero flux as the boundary condition of the radiative transfer in our GCM.Instead, the boundary condition for the thermal flux at the deepest layer of the atmosphere is zero, which corresponds to T int = 0 K.We display the effect of the intrinsic temperature on the temperature profile in Fig. 4, where we show temperature profiles with three different intrinsic temperatures (500 K, 600 K, and 700 K) along with the final time steps of the steroids and nominal model from Fig. 2. The temperature profiles for the different intrinsic temperatures have been calculated from the temperature evolution of the nominal model (Fig. 2) using Eq.(2) with bolometric fluxes, which was computed during the runtime of the climate simulations.The steroids model has a negative total net flux, rendering it impossible to calculate a physical meaningful intrinsic temperature.
Combining the predictions from Thorngren et al. (2019) and Sarkis et al. (2021), we find that the intrinsic temperature of WASP-76 b should be approximately 600 ± 100 K and the RCB location should be 1 ± 1 bar according to structure models.Thus, we can compare the intrinsic temperature to the predictions and find that the intrinsic temperature of the final nominal model is 434 K, which is significantly below the predicted value of 600 K from structure models.The right plot in Fig. 4 displays the RCB location for the same models.Comparing the RCB location of the final globally averaged nominal and steroids model to the predictions from the above-mentioned 1D structure models, we find that the locations of the RCB in our models are much deeper than 10 bar, and therefore not in line with predictions.Rauscher & Showman (2014) showed that the RCB location varies horizontally, because the equatorial regions have higher stellar fluxes than the polar regions.We therefore performed additional calculations of the RCB with temperature averages of the polar (|ϑ| > 45 • ) and equatorial (|ϑ| < 45 • ) regions.We can confirm the findings of Rauscher & Showman (2014), namely that the RCB location seems to vary horizontally, where we find higher and deeper levels for polar and equatorial regions, respectively.The different values of the RCB for different horizontal regions is caused by the fact that the temperature becomes horizontally homogeneous only at deeper levels in the atmosphere.We therefore conclude that the intrinsic model differences of 1D and 3D models make the comparison of RCB locations more difficult.
We show the evolution of the area-integrated radiative fluxes (F and F pla ) along with the total energy of the steroids and nominal in Fig. 5.We follow Polichtchouk et al. (2014) and calculate the total energy (TE) as a sum of the potential, kinetic, and internal energy by where Z is the geometric height.The plot shows that the high temperatures in the deep atmosphere do not lead to large flux differences between the stellar and planetary fluxes, which can be explained by the relation between flux and temperature, which goes with the fourth power of the temperature (Eq.( 2)).Intrinsic temperatures on the order of 1000 K therefore become insignificant against the temperature of the thermal emission, which is on the order of 2000 K.The nominal model loses energy during its evolution as the planet cools down, whereas the steroids model heats up again slightly.These trends manifest themselves in the right panel of Fig. 5, showing that the steroids model conserves the total energy, whereas the nominal model is still losing energy through radiation to space.Our model conserves the radiative energy to roughly 99.9%.Deitrick et al. (2022) recently reported that in practice their THOR GCM always tends toward a negative radiative flux balance, with discrepancies of as large as a few percent.The authors claim that the reason for that can be found in the numerical dissipation processes in their model such as sponge layers and hyper-diffusion2 .Similarly, Rauscher & Menou (2012) report that their total radiative heating rate is always positive, balancing a loss of kinetic energy due to numerics.The conservation of energy and angular momentum in MITgcm was benchmarked by Polichtchouk et al. (2014), revealing that MITgcm is particularly poor at conserving angular momentum, whereas it performs well at conserving energy, which can be confirmed for the steroids model in the right panel of Fig. 5.We therefore think that conservation of the radiative energy with an accuracy of 99.9% is already very good.It is therefore perfectly reasonable that the steroids model is indeed losing more energy than it is gaining, while still being in equilibrium.

Discussion
We conclude that we cannot find large-scale atmospheric motions able to inflate a planet significantly.Instead, we find that heating and cooling by radiation dominate the temperature evolution even in deep layers, forcing our nongray simulations toward a cool final state.We therefore propose that hydrodynamical models of hot gas giants need extra physics or parameterization (such as dissipative drag or diffusion) to recreate an inflated interior.Youdin & Mitchell (2010) and Tremblin et al. (2017) showed that vertical heat transport by turbulent mixing could possibly inflate the interior.However, their models do not model  vertical mixing self-consistently, instead relying on a parameter that sets the mixing strength.Using idealized GCM simulations, Sainsbury-Martinez et al. ( 2019) claimed that large-scale atmospheric motions could be a reason for a tendency toward a hot adiabatic interior in their simulations.There are two major differences between the model used in this work and the GCM used in the work of these latter authors that need to be considered when questioning the lack of a hot interior in our models in comparison to their work.One of the obvious differences is the treatment of radiation: whereas Sainsbury-Martinez et al. ( 2019) do not include radiative cooling and heating in the deep atmosphere, we apply a nongray radiative transfer scheme.
Perhaps equally as important, the second difference is the use of numerical dissipation, where Sainsbury-Martinez et al. ( 2019) use a hyper-diffusion scheme to stabilize their model, whereas we use a fourth-order Shapiro filter.Previous studies showed that in general, both dissipation schemes (filtering and diffusion) have a profound influence on the temperature and dynamics of a GCM (e.g., Polichtchouk et al. 2014;Koll & Komacek 2018;Hammond & Abbot 2022).Sainsbury-Martinez et al. (2019) showed that the strength of numerical diffusion in their model critically impacts the timescale of the tendency toward the hot adiabat, where a strong numerical diffusion relates to a faster evolution toward a hot interior.It could therefore be possible that their finding might have been caused by heat diffusion from their numerical dissipation scheme.
It is still unclear how the interior of hot and ultra-hot Jupiters interact with the upper atmosphere.There are currently few studies that link the convective interior with the photosphere dynamics (Showman et al. 2020).The first steps toward this goal were made by Lian et al. (2022), who model the effect of convective forcing on the circulation patterns of hot irradiated planets with a combination of Newtonian forcing in the upper atmosphere and perturbations in the deep atmosphere.A model able to combine 1D interior models with realistic 3D atmospheric models including radiative transfer would also be desirable and would reveal insights into the connection between the convection in the deep interior and the atmosphere.Similar approaches, connecting convective interiors with upper radiative atmospheres, exist in the stellar community (e.g., Skartlien et al. 2000;Freytag et al. 2010Freytag et al. , 2017Freytag et al. , 2019)).
The present work is the most complete GCM study to date, which is aimed at understanding the coupling of the upper and the deeper atmosphere in hot Jupiters.However, we note that we are so far omitting important effects for ultra-hot Jupiters: Ohmic dissipation, dissociation of H 2 and clouds.The upper atmospheres of ultra-hot Jupiters like WASP-76b are highly ionized due to their high temperatures (e.g., Helling et al. 2021).Magnetic fields, if present, would therefore couple strongly to the ionized winds.Rauscher & Menou (2013) and Beltz et al. (2022) showed that friction induced by magnetic fields acting on ionized winds may alter the dynamical state of the atmosphere and lead to additional heating in the upper atmosphere, where friction is converted to heat.It is still debated and unclear (Rauscher & Menou 2013;Showman et al. 2020) as to whether this extra heat could be sufficient to inflate a planet.
Recent simulations of ultra-hot Jupiters showed that H 2 dissociation can alter their thermal structure (Tan & Komacek 2019).The energy needed to dissociate molecular hydrogen on the day side cools the atmosphere, while recombination of atomic hydrogen heats the night side, leading to an overall reduced day-night temperature contrast.Similarly, the presence of clouds in a 3D GCM will alter the thermal and dynamical state of hot Jupiters (Roman et al. 2021;Deitrick et al. 2022).Realistic models of ultra-hot Jupiters should therefore take these effects into account, even though it is unlikely that these effects are important for the deeper parts of the atmosphere.

Summary and conclusions
In this work, we investigate one of the main hypotheses proposed by Guillot & Showman (2002) and Youdin & Mitchell (2010) to explain the abnormally large radii of inflated hot Jupiters: Energy transport from the irradiated atmosphere into deeper layers.We performed long-term (86 000 d) nongray 3D GCM simulations of WASP-76 b to search for signs of vertical advection of potential temperature from the upper irradiated atmosphere into the interior.Our simulations started from a hot initial state and cooled down toward a colder state.We then used an exponential fit on the temperature convergence to extrapolate the temperature toward a possible final converged state, which was reached soon after for the extrapolated simulation.The final converged L11, page 5 of 9 A&A 666, L11 (2022) simulation exhibits a cold temperature profile, lacking signs of vertical advection of heat from the upper atmosphere into the interior.We find instead that the atmosphere has cooled down to radiative equilibrium, conserving energy with an accuracy of 99.9%.
Our work strongly disfavors vertical downward advection of energy from the irradiated part of the atmosphere by large-scale atmospheric circulation as a possible explanation for inflation.We suggest instead that physical processes other than radiation and dynamics need to be taken into account in order to match the abnormal large radii of inflated hot Jupiters.Such processes could be the interaction between a convective interior and the atmosphere, small-scale turbulent transport, or magnetic field interactions.Future investigations are needed to quantify 3 their effects on the state of radius inflation 4 .Previous studies revealed the enormous impact of the chosen (often unjustified) dissipation schemes (Heng et al. 2011;Skinner & Cho 2021).We perform a test simulation with a 1000 times longer friction timescale (τ deep ) compared to the nominal and steroids simulations.This longer friction timescale weak-ens the strength of the Rayleigh drag applied to the deep atmosphere and makes sure that the prescription of the bottom drag does not affect the conclusions of this work.This additional simulation has been started from the nominal model at 40 000 d and has been carried out for 10000 d.
We compare the temperature evolution of the nominal simulation to the simulation with the low drag in Fig. A.1, where we find that the simulation with the lower drag cools down slower than the nominal simulation.The reason for the difference in the temperature evolution can be found in the Fig . A.2, where we show the temperature maps at 650 bar of the model with low drag and the nominal model.Both models start from an average temperature of 10 080 K at 40 000 d and cool down by a few hundred Kelvin.However, we find that the drag stops the equatorial jet, which otherwise descends toward the bottom boundary in the simulation with low drag.Simultaneously, stopping the equatorial jet leads to an efficient north-south redistribution of the temperature at the bottom boundary.Thus, in the model with much smaller drag, we find that the temperature starts diverging between the poles and the equator, where the equator cools down slower than the poles.
Considering this result, it is reasonable to conclude that the final converged solution of a model with lower drag would be hotter than the steroids model, possibly reducing the gap between outgoing and incoming radiation in the left panel of Fig. 5. Nevertheless, we also find that this simulation continues to cool down, still showing no sign of the vertical advection of heat that would, on the contrary, lead to a warming of the deep layers.
where λ is a norm that ensures that the value inside log stays positive and unitless.We set λ = j (dT/dt) j , (B.2) the subscripts i, j ∈ {1, 2, 3, 4} to denote the four support points of the fit.We can now use the least-squares method to solve where c = exp(c).The fit, as displayed in Fig. 3, matches the temperature change rates around t = 40 000 d.However, the real temperature change rates for later times are slightly below the fit, meaning that the temperature cools down even faster than predicted by the fit.
We then estimate the final temperature from integrating Eq. (B.4) over time and get We evaluate this equation at the age of the mother star, which is ≈ 5 Gyr (West et al. 2016), and find that the temperature at 650 bar should be T 650 bar = 4567 K.This temperature at the bottom of the simulation domain is then translated to an adiabatic temperature profile for the atmospheric layers below 10 bar and is used to kick off the steroids model.In order to force the simulation to continue with this new temperature profile for pressure levels below 10 bar, we developed a method that would smoothly force the temperature towards the predicted state during runtime.This method consists of a time-based smoothing, where the temperature is forced towards the final state by splitting the total change of temperature over a period of 10 d.We make sure that the temperature is changed by the expected amount, but still modulate the change rates by a sine function dT dt ∝ sin π t 10 d to smoothen the transition.A linear vertical smoothing between 1 bar and 10 bar ensures that the upper atmosphere above 1 bar stays untouched from the artificial forcing, which should only change the temperature in the deep layers.These measures are only used for the period of 10 day, and their only purpose is to force the temperature of the steroids model towards its new initial condition.The temperature forcing after 40010 d is then given by standard radiative heating and cooling as computed self-consistently within the GCM.nominal model, which can be explained by the strong winds in the steroids model that are more efficient in equilibrating the temperature in the steroids model.Similarly, we can see that the hot spot offset is greater in the steroids model, which is another outcome of the strong winds that transport energy further, before it can be reradiated.

Fig. 2 .Fig. 3 .
Fig. 2. Temperature evolution (color coded) of the horizontally averaged temperature in the deep layers of the atmosphere for the nominal model (solid lines) and the steroids model (dotted lines).A gray dot displays the position of the RCB at the final state.

Fig. 4 .
Fig.4.Temperature profiles and RCB locations.Left: temperature profiles for different intrinsic temperatures (gray) and final temperature profiles of the nominal (T int ≈ 434 K) and steroids simulation.Right: calculated RCB for the different temperature profiles of the left panel.Additional markers display the RCB calculated for the equatorial and the polar region.The gray profiles are calculated from the temperature evolution of the nominal model.We note that the steroids model loses energy due to dissipation and therefore has no physical intrinsic temperature.

Fig. 5 .
Fig. 5. Energy budget of the nominal and steroids simulations.Left: radiative energy budget FdA .Right: evolution of the total energy.Both simulations converge toward radiative equilibrium with approximately 0.1% accuracy.

Fig
Fig. A.1.Evolution of the temperature in three pressure layers (57 bar, 450 bar, 650 bar) in a simulation with low drag compared to the nominal simulation.The temperature in each layer is normed to the respective initial value at 40 000 d.

Fig. A. 2 .
Fig. A.2. Isobaric temperature slices at t = 50000 d at p = 650 bar for the nominal and the model with low drag.The temperature has been averaged over the last 100 d.The substellar point is located in the center of the plot.
3)in order to obtain the best-fit values for m and c, which is equiv-

Fig
Fig. C.2. Zonal mean eastward velocity for the nominal and the steroids model.
Fig. C.3.Phase curve and phase offset.Left: Phase curve at λ = 4.5 µm, computed using prt_phasecurve.The dashed gray line corresponds to the center of the dayside.Right: Difference between the dayside center and the maximum of the phase curve.
Overview of the simulations in this work.The steroids model starts at t = 40 000 d from the nominal model with the guess of the final state obtained via the method outlined in Appendix B.