Hot Jupiter Diversity and the Onset of TiO/VO Revealed by a Large Grid of Non-Grey Global Circulation Models

The population of hot Jupiters is extremely diverse, with large variations in their irradiation, period, gravity and chemical composition. To understand the intrinsic planet diversity through the observed population level trends, we explore the a-priori scatter in the population created by the different responses of atmospheric circulation to planetary parameters. We use the SPARC/MITgcm 3D global circulation model to simulate 345 planets spanning a wide range of instellation, metallicity, gravity and rotation periods typical for hot Jupiters, while differentiating between models with and without TiO/VO in their atmosphere. We show that the combined effect of the planetary parameters leads to a large diversity in the ability of atmospheres to transport heat from day-side to night-side at a given equilibrium temperature. We further show that the hot-spot offset is a non-monotonic function of planetary rotation period and explain our findings by a competition between the rotational and divergent parts of the circulation. As a consequence, hot-spot offset and phase curve amplitude are not necessarily correlated. Finally, we compare the observables from our grid to the population of Spitzer and Hubble observations of hot Jupiters. We find that the sudden jump in brightness temperature observed in the Spitzer secondary eclipse measurements can be naturally explained by the cold-trapping of TiO/VO at approximately 1800K. The grid of modelled spectra, phase curves and thermal structures are made available to the community, together with a python code for visualization of the grid properties, at https://doi.org/10.5281/zenodo.10785321 and http://sim3d.oca.eu/.


INTRODUCTION
As the number of high-quality exoplanet observation increases, with the James Webb Space telescope (Greene et al. 2016) and upcoming Ariel mission (Tinetti et al. 2017) providing wide wavelength coverage and ground-based telescopes supplying high resolution data (Ehrenreich et al. 2020), hot Jupiters will continue to provide the best targets for atmospheric characterization.These planets are primarily classified based on their short-period orbits and high temperatures, meaning they offer optimal targets for transit observations.Many have been observed during secondary eclipse (Garhart et al. 2020;Baxter et al. 2020;Deming et al. 2023), with a smaller number viewed over full orbits to provide phase curves (Parmentier & Crossfield 2018).Strong gravitational interactions, due to the short-period orbits, are thought to cause tidal locking (Rasio et al. 1996;Lubow et al. 1997;Showman & Guillot 2002), leading to rotation periods synchronised with orbital periods.As all the incident stellar irradiation is deposited into the day-side hemisphere, steep day-night temperature contrasts form, making these planets intrinsically 3-dimensional in nature and driving strong winds (Showman & Guillot 2002).These winds facilitate energy transport from the hot day-side onto the cold night-side (Komacek & Showman 2016;Komacek & Tan 2018), in the form of ★ E-mail: Alexander.roth@physics.ox.ac.uk either a super-rotating equatorial jet which can reach speeds of up to several km/s or direct day-to-night divergent flow Hammond & Lewis (2021).The winds and day-to-night flow have been measured in these atmospheres (Snellen et al. 2010;Brogi et al. 2016).The jets, first theorized by modelling (Showman & Guillot 2002;Cooper & Showman 2005;Showman et al. 2008;Menou & Rauscher 2009;Heng et al. 2011;Perna et al. 2012;Rauscher & Menou 2012;Mayne et al. 2014;Showman et al. 2015;Dobbs-Dixon & Agol 2013;Tan & Komacek 2019), have also been inferred through observation (Knutson et al. 2007).However, these strong winds are often not strong enough to transport energy from day-to-night before it is radiated back to space, leading to large horizontal contrasts in temperature (Knutson et al. 2007;Stevenson et al. 2010), chemistry (Parmentier et al. 2018;Helling et al. 2019) and cloud coverage (Parmentier et al. 2021;Helling et al. 2019).
The formation pathways responsible for the observed exoplanet population is one of the largest remaining unanswered questions.Planets form from the local material in the stellar disk.The location of formation, or subsequent migration after formation has occurred, could affect the atmospheres of giant planets (Bitsch et al. 2022;Fortney et al. 2013).Consequently, measuring chemical abundances in planetary atmospheres can inform us on their formation and migration history (Cridland et al. 2019;Öberg et al. 2011;Madhusudhan et al. 2014).At the population level, however, it is not only the trend of chemical abundances with parameters, but also their scatter around the mean population value, that can inform us on the planetary formation mechanisms and their stochastic nature (Fortney et al. 2013).However, abundance measurements in exoplanet atmospheres have to rely on a thorough understanding of their 3D thermo-chemical structure to avoid bias in atmospheric retrievals (Arcangeli et al. 2018;Taylor et al. 2020).Particularly, at the population level, we wish to separate the part of the observed trends and scatter that is due to the diversity of the planetary parameters (e.g.temperature or rotation period) and the one that is due to intrinsic, formation related, parameters (such as metallicity).
The population of hot Jupiters is extremely diverse.As shown in Figure 1, their equilibrium temperatures vary from ≈1000 to 2400K, which corresponds to irradiation varying by more than an order of magnitude.Cooler planets are often called "warm Jupiters" whereas hotter ones are called "ultra-hot Jupiters".Furthermore, gravity and period also vary by an order of magnitude within the population.Composition-wise, observations of their spectra have shown that atmospheric metallicity can also have an order of magnitude variation within the population (Welbanks et al. 2019).There are multiple known correlations between these parameters at the population level which significantly boosts the spectral diversity we observe and introduce trends in their atmospheres.For example, because of tidal locking, equilibrium temperature and rotation are intrinsically linked.Additionally, the hottest hot Jupiters have an inflated radius (Fortney & Nettelmann 2010;Baraffe et al. 2014), leading to a correlation between equilibrium temperature and surface gravity.Finally, a correlation between atmospheric metallicity and planet mass is observed in the solar-system and could also be present in giant exoplanets (Welbanks et al. 2019).
The chemical composition of an atmosphere can strongly alter its thermal structure.In particular, the presence of TiO and VO, two strong optical absorbers that exist in M dwarfs, have been postulated to create a thermal inversion on the planet day-side (Hubeny et al. 2003;Fortney et al. 2008), meaning that the temperature would increase with decreasing pressure.Such thermal inversions have been observed by the presence of emission features in the thermal emission spectra of the hottest planets (Coulombe et al. 2023;Mikal-Evans et al. 2020;van Sluĳs et al. 2023;Cont et al. 2022).For the cooler ones, clear absorption features have shown the definitive lack of a thermal inversion (Line et al. 2016;Kreidberg et al. 2014), which is expected as TiO and VO are supposed to condense out of the atmosphere at low temperatures.However, the exact point at which the population transitions from TiO/VO rich atmospheres with thermal inversions towards TiO/VO poor atmospheres without inversions is still unclear.Direct detection of TiO and VO through observation is complex as its spatial distribution is shaped by multiple processes.Condensation of these chemicals in the deep layers (Spiegel et al. 2009) or on the night-side (Parmentier et al. 2013;Merritt et al. 2020) could deplete these chemicals from the planet's day-side atmosphere (Powell et al. 2018).Furthermore, atmospheres that are hot enough to suppress condensation can be too hot for TiO and VO to be stable against thermal dissociation (Parmentier et al. 2018).
Additionally, the 3D structure of hot Jupiters makes the measurement of chemical abundances from their spectra challenging.Indeed, observations, either in transit or eclipse, are usually sensitive to the spatially averaged spectra.Information from different parts of the atmosphere, which can have very different properties, are mixed together.1D retrieval models can have strong biases when applied to hot Jupiter data sets, for both emission (Feng et al. 2016;Blecic et al. 2017;Taylor et al. 2020) and transmission (Line & Parmentier 2016;MacDonald & Madhusudhan 2017;Caldas et al. 2019;Pluriel et al. 2020;MacDonald et al. 2020;Lacy & Burrows 2020).This is due to the spatially varying cloud coverage (Line & Parmentier 2016), thermal structure (Feng et al. 2016;Taylor et al. 2020;Caldas et al. 2019) or chemical abundance maps (MacDonald et al. 2020;Pluriel et al. 2020).This problem becomes even greater at higher spectral resolutions, as the spectral line shapes, depths and positions, which are heavily dependent on the chemical, thermal and dynamical structure in the atmosphere, become increasingly resolved (Showman et al. 2013;Miller-Ricci Kempton & Rauscher 2012;Kempton et al. 2014;Zhang et al. 2017;Beltz et al. 2021;Keles 2021;Wardenier et al. 2021).This has been shown to cause heavily asymmetric features in transmission (Ehrenreich et al. 2020;Kesseli & Snellen 2021;Wardenier et al. 2021;Kesseli et al. 2022;Savel et al. 2022).These effects can lead to biases when inferring abundances because of the non-linear nature of radiative transfer: the spatially averaged spectra is different from the spectra of the spatially averaged thermal structure of the planet.This is particularly true when spatial correlations between physical quantities are present in the planet (e.g.clouds forming in the cooler parts of the atmosphere, chemical abundances being determined by the local temperature etc...).Consequently, to obtain robust chemical abundances from observed spectra, we must understand how planetary parameters influence the physical processes in hot Jupiter atmospheres.Therefore, to draw conclusions on observational trends arising from planet formation processes, we must first understand expected trends due to planetary properties.With rapidly increasing numbers of high quality exoplanet observations, large scale population studies (e.g (Baxter et al. 2020;Garhart et al. 2020;May et al. 2022;Deming et al. 2023)) are now possible.The goal of this paper is to estimate how the 3D atmospheric dynamics can shape the observational properties of hot Jupiters across the observed population of planets.
Such work has been done in 1D and 2D by the use of radiative/convective model grids (Baeyens et al. 2021(Baeyens et al. , 2022;;Goyal et al. 2020Goyal et al. , 2021;;Mansfield et al. 2021;Molaverdikhani et al. 2019).Baeyens et al. (2021Baeyens et al. ( , 2022) ) explore how chemical mixing and photochemistry affect atmospheres in a pseudo-2D framework.Goyal et al. (2020Goyal et al. ( , 2021) ) generate a 1D-2D model grid with a range of re-circulation factors, metallicities, and C/O ratios and use this to explore population trends in the observed Spitzer thermal emission.In a similar way, Mansfield et al. (2021) use a 1D model grid to explore trends in the water feature strength of HST/WFC3 observations.
A number of 3D Global Circulation Model (GCM) grids have also been used to identify the importance of planetary parameters in shaping the atmospheric circulation in hot Jupiters.Day-to-night temperature contrast has been shown to depend heavily on irradiation received, rotation period and atmospheric composition (Showman et al. 2015;Kataria et al. 2016;Komacek & Showman 2016;Komacek et al. 2017;Perna et al. 2012;Tan et al. 2024).Super-solar metallicities are also found to increase the day-to-night temperature contrast at the photosphere and reduce eastward hot-spot offset (Showman et al. 2009;Kataria et al. 2015Kataria et al. , 2016;;Drummond et al. 2018a).Additionally, Kataria et al. (2016); Komacek et al. (2022) show how surface gravity plays an important role when determining the photospheric pressure in the atmosphere.However, until now no study has utilised non-grey 3D modelling to investigated the impact of systematically and jointly varying all of these intrinsic planetary properties on an atmosphere, which is the ultimate goal for this study.
We begin by describing the dynamical, radiative transfer and chemical modelling in Section 2. In Section 3, we then detail the setup of the model grid.This includes a brief discussion on each of the free parameters, why certain values have been selected and how we expect them to affect an atmosphere.Then, in Section 4 we show how the wide range of planetary free parameters chosen influence the atmosphere of hot Jupiters.Starting with an overview of the thermal structure, then a discussion on the wide variation in day-to-night heat redistribution within the grid.We then move on to an overview of the secondary eclipse spectra and phase curves, before discussing the phase curve offset in more detail.Comparisons with current planetary observations are then detailed in Section 5. We provide an explanation for the sharp jump in brightness temperature that is seen by Spitzer observations of the population of hot Jupiter atmospheres (Deming et al. 2023), inferences from investigating spectral indices for water in HST/WFC3 observations (Mansfield et al. 2021), and compare predictions from our grid to the Spitzer phase curve observations before summarizing our conclusions in Section 6.

Dynamics
We use the non-grey SPARC/MITgcm global circulation model (Showman et al. 2009), which has been widely applied to mod-elling the atmospheric dynamics for hot, Jupiter sized planets in tidally locked orbits (Showman et al. 2009(Showman et al. , 2015;;Kataria et al. 2015Kataria et al. , 2016;;Steinrueck et al. 2019;Lewis et al. 2017;Parmentier et al. 2013Parmentier et al. , 2016Parmentier et al. , 2018Parmentier et al. , 2021;;Tan et al. 2024).This model couples the dynamical code from the MITgcm (Adcroft et al. 2004) with the plane-parallel radiative transfer code of Marley & McKay (1999) and solves the primitive equations on a cube-sphere grid.
The setup of the GCM here is very similar to that described in Parmentier et al. (2018Parmentier et al. ( , 2021)).All models assume an internal heat flux corresponding to a temperature of 100K.Although it has been shown that internal heat flux probably increases with equilibrium temperature (Thorngren et al. 2019), the effect of the internal heat on the wind and thermal structure at the photosphere is not yet fully understood (Komacek et al. 2022).We therefore leave such study for future work and choose to keep the internal heat constant here.All of our models also use a fixed planetary radius of 1.3  , where   is the equatorial radius of Jupiter at the 1 bar pressure level (71, 492).Despite the radii of hot Jupiter's varying by around a factor of 2 across the population, this was done in order to minimize the number of parameter dimensions within this study.We assume a specific heat capacity of the atmosphere   = 1.3 × 10 4   1  −1 , with a specific heat capacity ratio of  = 1 + 2/7 and a mean molecular weight of  = 2.3which is valid for the H2-dominated atmospheres we are simulating.Other prescribed model parameters are changed on a model to model bases, including: equilibrium temperature, atmospheric metallicity, rotation period (changed with stellar parameters), surface gravity and the inclusion of certain molecular species, see Section 3 for details.Our models purposely do not add a Rayleigh drag parameterization as has been done in the past by e.g.Komacek & Showman (2016) and most dissipation is due to the 4th order Shapiro filter applied to the temperature and velocity fields, which dissipates kinetic energy by smoothing the horizontal gradients between grid points (Koll & Komacek 2018).We use a fixed dissipation timescale of  shap =50s for all models, meaning that the numerical dissipation is proportional to the kinetic energy of each model (Koll & Komacek 2018).This provides us with the least dissipative models that can be obtained at the spatial resolution we are using while being numerically stable.
The atmosphere was allowed to evolve for 150 days, 300 days for the 1000-1200K models, with all quantities being averaged over the final 50 days for the output.This is enough integration time for a pseudo-steady state to be reached at the photospheric pressure level (Showman et al. 2009), although not enough to reach a fully converged state which would take several thousand days (Mendonça et al. 2018a;Mendonça 2020;Wang & Wordsworth 2020;Mayne et al. 2017).As shown in Appendix E, running for much longer times does not significantly change the observable properties of the planets we are interested in.Given that the interactions between the deep interior and upper atmosphere are still not fully understood for hot Jupiter exoplanets, we are assuming that the deep layers are quiet as they don't have time to develop winds, giving us the same solution as a photosphere equilibriated with a quiet deep interior and meaning that our choice of deep boundary condition has no effect on the final result.The models were run at a horizontal resolution of C32, which corresponds to an approximate resolution of 128 longitudinal and 64 latitudinal cells.The atmosphere is split into 53 pressure levels, between 2×10 −6 -200 bars.For low metallicity models the time-step was set to 25 s, which was decreased to 12.5 s at higher metallicities in order to stabilize the simulations.

Radiative Transfer
In both the 3D GCM simulations and subsequent spectral calculations, radiative transfer is handled via the plane-parallel radiative transfer code of Marley & McKay (1999).Although this code was first developed to study Titan's atmosphere (McKay et al. 1989), it has previously been used to study both hot Jupiters (Fortney et al. 2005(Fortney et al. , 2008;;Parmentier et al. 2021) and ultra-hot Jupiters (Parmentier et al. 2018).
For the opacities, we use those detailed in Freedman et al. (2008), including updates from Freedman et al. (2014).When coupled to the GCM a lower spectral resolution including only 11 frequency bins is used in order to speed up calculations (Kataria et al. 2013).Within each of these bins, thousands of spectral lines are modelled through the correlated-k distribution method (Goody & Yung 1989;Lacis & Oinas 1991), allowing us to compress the information down to just 8 k-coefficients.
We then use the 3D thermal structure from the GCM and post process using a much higher spectral resolution of 196 wavelength bins, ranging from 0.2679 -227.5313.This is done using the twostream radiative transfer equations along the line of sight for each atmospheric column in the cube-sphere grid at each planetary phase, taking into account emission, absorption and scattering.By using this method we naturally take into account geometrical effects such as limb darkening.For more details see Parmentier et al. (2021Parmentier et al. ( , 2016)); Fortney et al. (2006).For the stellar spectra we use the NextGen stellar models (Hauschildt et al. 1999).Spectra from the grid produced using this method can be seen in Section 4.1.
To additionally speed up post-processing of the grid, the P-T structure from the models was recalculated with a lower spatial resolution of C16, equivalent to 64 longitude and 32 latitude cells.This is done as an insignificant portion of the spectral information is lost when decreasing the spatial resolution in this manner, whilst a significant amount of time is saved in the calculations (Robbins-Blanch et al. 2022).
The opacity affects the radiative timescale, and therefore temperature structure, of an atmosphere.Opacity itself is highly dependent on the temperature structure due to pressure and temperature broadening of spectral lines and changing chemistry (Parmentier et al. 2018;Arcangeli et al. 2018).Using non-grey opacity is therefore important for realistic modelling to ensures that energy is conserved when we calculate the spectrum of the planet from the modelled temperature field.

Chemistry
In our models we assume local chemical equilibrium with local rainout of condensate materials (Visscher et al. 2006(Visscher et al. , 2010) ) and a solar C/O ratio of 0.46 (Lodders et al. 2009).Both molecular and atomic chemical abundances are calculated using a version of the NASA CEA Gibbs minimisation scheme (Gordon & McBride 1994), assuming solar elemental abundances and local chemical equilibrium, developed to study gas and condensate equilibrium chemistry under different atmospheric conditions (Moses et al. 2013;Skemer et al. 2016;Kataria et al. 2015;Wakeford et al. 2017;Marley et al. 2017;Parmentier et al. 2018).Chemical equilibrium is far from a safe assumption in hot Jupiter atmospheres, where the extreme conditions cause effects such as photochemistry and chemical quenching which smooths out the very large day to night chemical gradients (Cooper & Showman 2006;Drummond et al. 2018b;Mendonça et al. 2018b;Steinrueck et al. 2019).However, although non-equilibrium abundances (through quenching) can change the overall heat transport and the observables in specific band-passes (i.e temperature changes of ∼ 100-200K can lead to a ∼ 20% change in the spectra (Steinrueck et al. 2019)), the effects of disequilibrium chemistry on the overall atmospheric energy balance is small compared to the population trend (Steinrueck et al. 2019;Drummond et al. 2016).
For the models with increased metallicity, we multiply the elemental abundances of all elements apart from hydrogen and helium by the metallicity factor and adjust the H/He abundance to ensure that the sum of all mixing ratios is equal to 1.To the first order, atmospheric metallicity increases the abundance of all molecules present in the atmosphere, overall increasing the opacity of the atmosphere.At second order, atmospheric metallicity can change the balance between molecules in chemical equilibrium and, for example, favour CO instead of CH4 in cool planets (JWST Transiting Exoplanet Community Early Release Science Team et al. 2023).Additionally, we include models both with and without TiO and VO in the grid, more details on these species can be found in Section 3.5.

GRID PARAMETERS
In the grid, we have five free planetary parameters which are varied: equilibrium temperature ( eq ), atmospheric metallicity (Log(M/H)), rotation period (determined via the stellar parameters), surface gravity (Log(g)) and the presence of strongly absorbing chemicals TiO and VO.The full range of values used can be seen in Table 1.Any other necessary model parameters (i.e C/O ratio, planet radius, internal temperature etc..) remain consistent across the grid.The values used for these fixed parameters can be seen in Section 2. Combinations of these free parameters were used to run a total of 345 simulations.The extent of our model grid when compared to the observed exoplanet population can be seen in Figure 1.From this, we note that the population scatter is well covered within the grid, with the vast majority of planets within our equilibrium temperature range having a closely analogous model.

Equilibrium Temperature
The equilibrium temperature,  eq , is a measure of the stellar irradiation received by an atmosphere.It is defined as the brightness temperature that a homogeneous planet with zero albedo would have.This temperature is incredibly important for the atmospheric structure, as it affects the dynamical, radiative and chemical timescales significantly.In our grid we have selected a wide range of equilibrium temperatures, from 1000K (which is on the lower end for hot Jupiter exoplanets) to 2400K.We sample every 200K giving us eight values.We don't use temperatures higher than 2400K as above this molecular dissociation (Parmentier et al. 2018), transport of heat through H2 dissociation (Bell & Cowan 2018;Roth et al. 2021;Komacek & Tan 2018) and atomic iron opacities (Lothringer et al. 2018) start to cause significant complications to the planetary energy balance.Magnetic drag may also start to disrupt the equatorial jet (Beltz et al. 2022;Rogers 2017;Rogers & Komacek 2014).On the other end, we don't use temperatures lower than 1000K as for these planets the assumption of tidally locked orbits starts to break down (Showman et al. 2014;Parmentier et al. 2015a).
The distribution of observed equilibrium temperatures is displayed in Figure 1.From this we can see that the majority of hot Jupiters have equilibrium temperatures between ∼ 1000 − 1800K, with a lower number at higher temperatures, and that our grid range covers most of the observed population.

Orbital Period
In order to stay close to the observed population, we decided to run models for planets orbiting three different stellar types.A K2 dwarf, similar to HD189733, a G0V dwarf, similar to HD209458 and a F6 dwarf similar to HAT-P-7.The influence of different stellar spectra on planet atmospheres, and particularly the sensitivity of temperature inversions to the spectral type, have been previously studied using 1D models for both hot (Mollière et al. 2015) and ultra-hot (Lothringer & Barman 2019) Jupiters.Associated stellar spectra come from the NextGen stellar model archive (Hauschildt et al. 1999).As such,  star ,  star and  star all vary between the different model tracks.Even though the exact stellar spectrum surely has an effect on the radiative energy balance of the atmosphere, the main effect of changing stellar parameters that we are interested in is the related change of orbital period.Because we assume that our planets are tidally locked, their rotation period is directly determined by the equilibrium temperature and stellar parameters: eq 1 3 days. (1) As shown in Table 2, our three stars correspond to three model tracks with periods shifted by 0.37 and 2.23 compared to our nominal case.Thus, at a given equilibrium temperature our grid has rotation period varying by a factor of 6.The goal of varying the orbital period in this way is to compare how the atmosphere varies when looking at different planets around the same stellar type.The distribution of observed orbital periods for the hot Jupiter population is displayed in Figure 1.We note that for high temperature models (>2000K), we do model some orbital periods that are lower than observed in the population, which was done to maintain the completeness of the grid.

Metallicity
We define the planetary metallicity here as the atomic molar mixing ratio of all other chemical elements to Hydrogen and Helium in the atmosphere.Elemental abundances are scaled from solar composition (Lodders et al. 2009), with all atoms other than H and He multiplied by the metallicity value.For the grid we have selected 3 values to sample in the metallicity: Log(M/H)=0.0, 0.7 and 1.5.These range from solar to roughly thirty-times solar metallicity.
We do not explore very high values of metallicity, ensuring that our approximation of a constant mean molecular weight and heat capacity throughout the population stay correct (see Zhang & Showman (2017)).The main effect of metallicity is therefore to increase the molecular opacities of the atmosphere, which leads to lower photospheric pressures.

Surface Gravity
The surface gravity of a planet is dependent on the mass and radius, and displays order of magnitude variations in the hot Jupiter mass regime.For example WASP-18b (Hellier et al. 2009) is a very high density planet with a surface gravity of log(g)=∼2.2 −2 , whilst there are also very low density planets, such as WASP-39b (Maciejewski et al. 2016) which has a surface gravity of log(g)=∼0.63 −2 .We therefore sample three values for the surface gravity, representing close to the extremes in the hot Jupiter population: log(g[SI]) =0.8, 1.3, 1.8 (g=6.31, 19.95, 63.10  −2 ).As shown in Zhang & Showman (2017), gravity naturally cancels out in the dynamic equations.However, the role of gravity is important in setting the pressure of the infrared photosphere and thus the opacities at the photosphere.

TiO & VO
If strong UV/optical absorbers are present in the upper atmosphere, a portion of the incident stellar irradiation will be intercepted creating an area of local heating which causes a temperature inversion.In hot Jupiter atmospheres, the primary chemicals thought to be responsible for these temperature inversion are gaseous titanium oxide (TiO) and vanadium oxide (VO) (Hubeny et al. 2003;Fortney et al. 2008).These species are predicted by chemical equilibrium to be present in sufficient quantities to produce a strong thermal inversion above temperatures of ∼ 1400, creating two distinct regimes in which hot Jupiter atmospheres can be classified (Fortney et al. 2008;Parmentier et al. 2015b;Piette et al. 2020) TiO and thermal inversions have indeed been detected in the hottest hot Jupiters, such as WASP-33b ( eq = 2781K) (Cont et al. 2022), WASP-121b ( eq = 2450K) (Evans et al. 2016), WASP-189b ( eq = 2636K) (Prinoth et al. 2022) and WASP-18b (  = 2505K)  ( eq = 2182K) (Pelletier et al. 2023).However, TiO, VO nor thermal inversions have been detected in cooler planets, such as HD209458b ( eq = 1477K) (Line et al. 2016) or WASP-43b ( eq = 1379K) (Kreidberg et al. 2014), even though these planets are hot enough to have gaseous TiO in their atmospheres.It is therefore not yet clear at which equilibrium temperature the population of planets transitions from TiO/VO dominated atmospheres to atmospheres devoid of TiO/VO (Mansfield et al. 2021).The main mechanisms through which to deplete TiO/VO from the day-side atmosphere are condensation processes that can either happen in the deep layers of the atmosphere (Spiegel et al. 2009;Parmentier et al. 2016) or in the planet's nightside (Parmentier et al. 2016;Beatty et al. 2017;Pelletier et al. 2023).
The actual strength of these cold-trap processes depends on poorly constrained microphysical parameters (such as the condensates size) and the strength of the atmospheric mixing (Powell et al. 2019).
Here, we decide to be agnostic about the transition from an atmosphere containing TiO/VO to one without.We therefore run both a set of models at chemical equilibrium (and thus containing TiO/VO) between 1000-2400K, and a set of models where the elemental abundances of Ti and V have been artificially set to zero, mimicking the disappearance of TiO/VO due to non-local condensation processes, between 1400-2400K.

Atmospheric Structure
For much of the analysis in this investigation we will be boiling the atmosphere down into a series of characteristic observable features.However, we show in Figure 2 the wide range of diversity displayed by our models.This shows the temperature structure, normalised by the equilibrium temperature, at the pressure level closest to the 1.4 micron photospheric pressure for all the models within our grid.The 1.4 micron photospheric pressure can be calculated as (Parmentier et al. 2018), where  H 2 O,1.4m is the absorption cross section of water at 1.4, which we assume to be constant,  H 2 O is the volume mixing ratio of water,  is the gravity and  is the mean molecular weight.this can be scaled by our model parameters as, (3) This method of calculating photospheric pressure works well for lower temperature planets, however it does become less accurate for higher temperatures as it does not include any treatment for water dissociation.Each model is sorted into columns and rows based on their grid parameters.Gaps in the grid shown in Figure 2 are due to models that are unable to converge.This usually happens at high temperature, high metallicity and fast rotation rates due to numerical instabilities.
The temperature structure in our models varies drastically based on the prescribed model parameters.Surface gravity and metallicity alter the location of the photospheric pressure level, orbital period changes the width and strength of equatorial jet circulation and TiO/VO create temperature inversions.From these temperature maps we can already start to see some trends among the models.The equatorial jet structure is clearly visible in many.This jet becomes more narrow in latitude for models with higher temperatures, which are rotating much faster than their cooler counterparts in the grid, and for models with lower metallicity.High latitude Rossby waves are also visible for many of the models.The day-to-night temperature contrast can clearly be seen to increase with the equilibrium temperature.
Figure 3 shows the zonal-mean zonal wind speed with pressure for the grid arranged in the same format.Again, there are some clear trends that emerge with our model parameters.The strength of the super-rotating jet is clearly increased with the equilibrium temperature.The jet also extends to deeper pressures in the atmosphere as the surface gravity is increased.As the normalised rotation period is Rotation Period P Norm = 0.37 P Norm = 2.23 increased, the jet also clearly becomes wider, and the flow structure transitions from multiple jets to a single one for higher temperature models.
To show how the secondary eclipse spectra vary with our model parameters, we include Figure 4.Here we can see that increasing metallicity and surface gravity act to increase the emission flux, with increasing normalised period having the opposite effect.This is caused by the effect of parameters on the heat redistribution, which is discussed in detail throughout Section 4.2.It is interesting to note here that the metallicity appears to have a larger affect on the flux within the emission bands, surface gravity has a larger affect outside of the bands and the normalised period has a fairly constant affect across the spectra.The bottom right panel of Figure 4 shows how the presence of TiO and VO affect the emission spectra.

Heat Redistribution
The incident stellar flux on a planet's atmosphere has, essentially, two main pathways.It is either reflected back into space, the portion of which can be determined by the bond albedo, or it is deposited in the day-side atmosphere where it can either be re-radiated or partially advected onto the night-side before being re-radiated (Cowan & Agol 2011).This can be described as the heat redistribution efficiency.Ideally, measurements of the emitted light from a planet at all phases is necessary to estimate the heat redistribution efficiency.However, for planets with a small albedo, this efficiency can be determined by comparing the incoming radiation ( 4 eq ) to the outgoing longwave radiation from the planetary day-side ( 4 day ).Thus we define the heat redistribution efficiency, or  -factor, as (Arcangeli et al. 2018), day is the temperature of the equivalent blackbody emitting as much energy in our direction as the planetary day-side.Following this definition, the redistribution factor (  -factor) has a minimum value of  = 1, corresponding to a case of maximum planet wide heat redistribution where the atmospheric temperature is fully homogenized between the day and night-sides, and a maximum value of  = 2.66 where there is zero heat redistribution between the two hemispheres.These limits only apply when considering the bolometric flux.A value of  = 2 corresponds to a theoretical case where there is dayside only heat redistribution.These  -factor values are related to the  values defined in (Cowan & Agol 2011) by Equation 4 is used to calculate the heat transport in our cloudless models, using the same method as described in Parmentier et al. (2021), since the albedos are always small (see Appendix D for our model albedos).This will further enable a direct comparison with secondary eclipse data.

Expected scaling for the heat redistribution factor
Heat is transported in hot Jupiter atmospheres via two mechanisms: direct advection of heat by the winds and wave adjustment mechanisms.As the atmosphere tries to move heat from day to night-side, energy is lost by the gas through radiation.Hammond & Lewis (2021) show that these two regimes can both contribute to the heat circulation efficiency equally for certain planets.As the jet circulation is highly dependent on a host of external factors such as the orbital period of the planet, these components will also have variable relative magnitudes.This means that their are essentially three ways to transport heat: wave adjustment mechanisms, rotational circulation and divergent circulation.Because of this, heat transport varies widely over the observed hot Jupiter population as it is highly dependent on atmospheric thermal and chemical structure.Fortney et al. (2008), for example, suggested that the presence of stratospheric temperature inversions could explain some variations in observed redistribution efficiency, which has driven much of the investigation into low pressure optical absorbers in hot Jupiter atmospheres.Additionally, heat redistribution is thought to be negatively effected with increasing temperature (Komacek & Showman 2016; Perez-Becker & Showman 2013) and metallicity (Charnay et al. 2015).
The wave timescale can be defined as as the timescale for gravity waves to travel across a hemisphere of the planet in the isothermal case (Perez-Becker & Showman 2013;Parmentier et al. 2021): where   is the planetary radius,  the mean molecular weight,   the Boltzmann constant and  photo the photospheric temperature.
The advective timescale can be defined as the timescale it takes for  a parcel of gas to be advected across one hemisphere by the equatorial zonal jet (Parmentier et al. 2021): where  jet is the mean wind speed at the equator and  eq the equilibrium temperature.Parmentier et al. (2021) find  jet to increase linearly with  eq in their cloudless GCM models, which we also find for our models (see Appendix A).  jet also increases linearly with increasing surface gravity.Increasing atmospheric metallicity slightly increases  jet at lower  eq , but reduces it at higher  eq .Additionally, increasing rotation period decreases  jet , although this relationship is non-linear.However, order of magnitude variations in these parameters change  jet by at most a factor of 2. The radiative timescale, the time it takes for a parcel of gas to lose its energy to space, at the photosphere of an atmosphere can be approximated as (Showman & Guillot 2002): where  is the Stefan-Boltzmann constant,  the surface gravity,  photo the photospheric pressure and   the specific heat capacity at constant pressure.2002).Given the radiative timescale has a much stronger dependency on the temperature and gravity parameters than the advective and wave timescales, we expect adjustments to the radiative timescale to dominate changes in the heat redistribution.The opacity will increase with atmospheric metallicity, as a higher ratio of metals to hydrogen will increase the atomic weighting of the atmosphere.From Equation 7, the radiative timescale at the photosphere will therefore reduce with increasing metallicity (Kataria et al. 2015;Drummond et al. 2018a), leading to less efficient heat redistribution.However, as  phot also depends on gravity, the scaling with  is not clear until we rewrite Equation 7. The radiative timescale can equally be expressed in terms of the opacity, .By first taking the hydrostatic equation and the definition of optical depth, where  is the optical depth.Then noticing that  = , we can rewrite the pressure as, where  ∼ 2/3 is the optical depth at the photosphere.Assuming constant gravity and opacity, the photospheric pressure can then be described as.
Because pressure and optical depth vary exponentially in our atmosphere, and because the opacity usually increases with pressure, the opacity value that matters to determine the photospheric pressure is the opacity at the photosphere,  photo .Substituting this back into equation 7, the radiative timescale then becomes, As  increases with pressure, due to pressure broadening of spectral lines,  rad will decrease with increasing surface gravity.The exact scaling relation between opacity and pressure can be found in Freedman et al. (2014).In conclusion, we expect planets with high metallicity and high gravity to be the less efficient at transporting heat from day to night.

Overall model distribution
The bolometric redistribution factor values for all models in the grid are shown in Figure 5, with those models containing TiO/VO on the right and those without on the left.There is a trend towards higher  -factor values, meaning weaker day-night heat redistribution, at higher equilibrium temperature models in our grid.This is a well documented result that has been previously shown in both non-grey (Parmentier et al. 2021) and semi-grey modelling (Komacek et al. 2017;Perna et al. 2012), and also predicted with energy balance models (Cowan & Agol 2011).The scaling is caused by the strong temperature dependency of the radiative timescale ( 3 from equation 12) compared to the scaling of the wind and wave speeds ( and √  from equations 6 and 5 respectively).Overall, hot planets do not move heat around fast enough compared to the time needed for a parcel of gas to cool down radiatively.This naturally explains why the night-sides of hot Jupiters are almost independent of the equilibrium temperature (Parmentier et al. 2021).
There is also a noticeable kink in the distribution of  -factors for models not containing TiO/VO.This is due to the competition between the radiative, advective and wave timescales.Below ∼ 1600 the photospheric pressure stays roughly constant (see Figure 4 in (Parmentier et al. 2021)), as opacity is increasing with the temperature but the Plank function is shifting to shorter wavelengths.Over this range of temperature,  -factor does not vary strongly with  eq .Over the same range of  eq , the radiative and advective timescales stay very close to each other.This is because the advective timescale is not inversely linear but inversely affine (linear minus constant alpha, see Parmentier et al. (2021) equation 6).As a consequence, when both radiative and advective timescales cross each other they do not scale as  2 but  3 /( eq − ), which shows a much weaker variation with  eq .
When TiO and VO are present, right panel of Figure 5, the heat redistribution is worse, making the  -factor larger.The reason is multi-fold.First, the opacities do not decrease with decreasing wavelength anymore, meaning that, contrary to the case without TiO, as the temperature increases the photospheric pressure decreases, leading to shorter radiative timescales.Second, a portion of the stellar energy is absorbed at very low pressure and re-emitted directly from the stratosphere, rather than transported to the IR-photosphere, hence decreasing the overall rate of day-night heat transport.
The main takeaway here, however, is that at each equilibrium temperature we find that there is a significantly large spread in the heat redistribution based on the other planetary parameters.When gravity, period and metallicity are varied together, they can alter the efficiency of the heat redistribution as much as the equilibrium temperature itself.If we compare the spread in these values to that found in Parmentier et al. (2021), showing the change in redistribution factor between cloudless and cloudy atmosphere models when altering only the equilibrium temperature, we find that the intrinsic properties of the planet are having a comparable effect on the atmosphere to that of night-side clouds.Furthermore, the spread is also equivalent to that found in Komacek et al. (2017) between models with very low drag timescales ( drag = 10 3 ) and infinite drag timescales.
The spread in heat redistribution factor at a given equilibrium temperature in our models, due to the combined variability of rotation period, gravity and metallicity, is well in line with the spread observed by May et al. (2022) who measure a heat redistribution factor varying between 0.6 and 2.6 for planets between 1100-2400K (in the Spitzer bands the  -factor can be below 1).May et al. (2022) also display the heat redistribution to have a very wide range within a narrow variation of equilibrium temperatures (  =0.6-2.4 between 1100-1700K).This is an important point to note, as due to the wide range of effects that different parameters have on the redistribution, there is a large amount of degeneracy between the efficiency of heat redistribution between cases.This greatly complicates the characterisation of hot Jupiter atmospheres, as the intrinsic properties of the planet, chemical constituents of the atmosphere and the presence of night-side clouds all act to effect the redistribution in a wide range of potentially offsetting ways.

Specific trends with parameters
We now look into more detail at the trends with each parameter and discuss whether it is in agreement with simple timescale analysis.To take a closer look at how each planetary parameter is individually affecting heat redistribution, we can plot the  -factor as a function of temperature for a subset of models within our grid (ignoring crossterms between the parameter space for now).This can be seen in Figure 6.
The change in  -factor with atmospheric metallicity is found in the left panel of Figure 6.Here, we find the calculated  -factor to increases with higher atmospheric metallicity.This is caused by higher metallicity increasing the opacity and therefore shifting the photosphere to a lower pressure (as explained in Section 4.2.1).This reduces the radiative timescale at the photosphere, which increasingly outweighs the dynamical timescale in the atmosphere and therefore reduces the day-night heat redistribution efficiency.
In the right panel we see that the  -factor increases with surface gravity.From the re-written radiative timescale equation (see Equation 12), we can see that, as the atmospheric opacity increases with surface gravity, the radiative timescale at the photosphere decreases with surface gravity therefore leading to a decrease in heat redistribution (again see Section 4.2.1).
At a given equilibrium temperature, a longer orbital period is shown to increase the heat redistribution efficiency (decrease the factor) in the middle panel Figure 6.This is caused by the increased strength of the combined rotational and divergent components of the circulation, as both transport heat from the day to the night-side and are both dependent on rotation period (Hammond & Lewis 2021;Lewis & Hammond 2022).As we will see in more details in Section 4.3.1, the opposite is true for the offset, for which it is the relative strength of divergent and rotational components of the circulation that matter.
Comparing the models excluding and including TiO/VO shows us that the presence of these species reduces redistribution across the board.Strong absorption in the near UV/optical range by these chemical in the high altitude day-side atmosphere, where the radiative timescale is short, reduces the efficiency of energy redistribution and significantly increases the temperature compared to the nightside.The average effect of the TiO/VO on the heat redistribution is Δ  = 0.216, but this varies significantly with temperature.

Combined parameter trends: heat redistribution
We now discuss the quantitative effect that each parameter has on the heat redistribution.Because the change of heat transport with, say, period, depends on the other parameters, such as temperature, we propose to look at both the average effect of the parameters and the variation of this average effect within the grid.For this, we measure the difference in heat redistribution between two models when varying only one parameter and keeping the other ones constant: Here  1 is our model with a period normalization factor of  norm = 0.37 and  2 the model with  norm = 1.We then take the statistical mean between all the Δ values for the combinations of  eq , log(M/H) and log(g).This gives us the statistical average for the change we would expect to see in the heat redistribution when transitioning between consecutive grid points for the period, i.e Δ   1,2 .This process is then repeated for each of the other parameters.The standard deviation in the Δ  values is also calculated, i.e    1,2 .The mean and standard deviation for each parameter Δ  is shown in Figure 7.
The size of the bar is the averaged change of the heat redistribution parameter when changing metallicity, pressure or gravity.The size of the error bar is the variation of this change when evaluated at different values of the other parameters.For example, one can see that the impact of period on the heat redistribution is similar for all metallicities and gravities.Hence it has a small variance.On the contrary, the change of heat redistribution due to a change of period depends strongly on the specific metallicity and gravity of the planet, hence the larger variance.
Overall, we see for the case without TiO/VO that gravity, period and metallicity have a quantitatively similar effect on the heat redistribution at a given equilibrium temperature.For the case with TiO/VO (right panel), the variation with gravity has a much smaller effect than the other parameters but both metallicity and period contribute in a similar way.However, atmospheric metallicity is largely unconstrained in comparison to period and gravity, which vary by no more than an order of magnitude within the population (see Figure 1).If the metallicity variation is strong enough across the population, its heat redistribution trend could therefore dominate the observed population trend.

Bolometric Phase Curves
The phase curve of a planet, characterised as the change in brightness temperature observed during a full orbit, contains information about the longitudinal structure of the planet's atmosphere.As the phase curve amplitude and shape are determined by longitudinal inhomogeneities for transiting planets, they contain useful information about the 3D nature of planets.This includes the contrast in brightness temperature between day and night, which is directly linked to the phase curve amplitude, and the asymmetry of the brightness distribution, which is directly obtained by the measurement of the phase curve offset.

Offset
The phase curve offset is a measure of the longitudinal offset of the brightest hemisphere from the substellar point (0 • ).This has a positive value if offset eastward (respective to the direction of planetary rotation) from the substellar point, meaning that the brightest hemisphere would be observed before secondary eclipse.Over the population of planets, phase curve offset is firstly determined by equilibrium temperatures, with offsets varying between 50 • in the colder cases to a few degrees for the hottest planets.Within each temperature bin, however, the variation of the phase curve offset is extremely large.For example, for planets at 1400K it ranges between 50 • for low gravity, low metallicity cases to 10 • for high metallicity, high gravity cases.Additionally, the phase curve offset is predicted to vary in a non-monotonic way with temperature (Komacek & Showman 2016;Zhang et al. 2018).
The bolometric phase curve offsets for a subset of our models are displayed in Figure 8.For the most part the trends we find can be explained using the same logic as that of the heat redistribution.A reduction in the efficiency of heat redistribution naturally leads to lower phase offsets at higher temperatures and the trends in both increasing metallicity and surface gravity can again be explained by the reduced pressure of the photosphere, where a shorter radiative timescale will result in a lower hot-spot shift from the substellar point.The presence of TiO and VO significantly reduce the phase offset in all cases, again due to the reduction in heat redistribution efficiency caused by the strong near UV/optical absorption in the upper atmosphere of the day-side.As mentioned in Section 4.2, the presence of TiO/VO causes a portion of the energy to be re-emitted directly from the stratosphere instead of the IR-photosphere, which leads to a smaller hot-spot offset.
The phase offset appears to peak between equilibrium temperatures of 1400-1600K, with models at lower temperatures having reduced offsets.This behaviour is also found in Komacek et al. (2017) (see their Figure 9).This is due to homogenisation of temperature at these lower equilibrium temperatures, caused by longer radiative timescale, which can be seen in Figure 2. Because of this, their phase offset is determined by small, wavenumber-2 variations of the thermal structure.The affect of rotation period on the phase offset (via the implied planetary rotation rate) is more complex, which we discuss in the next section.

Offset trend with period
As we saw in Section 4.2, larger rotation period leads to more efficient heat redistribution (smaller  -factor).Based on simple energy balance models (Cowan & Agol 2011), we would therefore expect phase curve offsets to increase with increasing period, as was proposed by May et al. (2022) based on Spitzer phase curve observations.However, as shown in Figure 9, this is only true for small enough rotation periods.Indeed, our models systematically show a turning point, where, past a critical rotation period, the phase curve offset starts decreasing with increasing period.This turning point, however, appears to be a function of equilibrium temperature.
We can explain this trend, and the fact that heat redistribution and phase curve offsets have different behaviour, by considering the competition between the divergent and the rotational component of the circulation (as in Hammond & Lewis (2021); Lewis & Hammond (2022)) and their impact on the Doppler shifting of Rossby waves.Whereas both components contribute to the heat redistribution, only the rotational component contributes to the phase curve offset.Therefore, planets with a strong divergent component can have a small phase curve offset despite having a significant heat transport efficiency.The position of the peak due to the rotational circulation is approximated by Equation 12of Hammond & Pierrehumbert (2018), from which we derive the following expression for the longitudinal shape () of the wavenumber-1 rotational circulation: where  is longitude (with  = 0 at the substellar point),  rad is the radiative timescale,   is the Rossby wave frequency (which increases with increasing Ω), and Ū is the jet speed.This equation is valid only in the rotationally dominated regime (e.g. for short rotation periods).It shows the balance between the radiative damping (1/ rad term), that tends to drag the offset towards zero, the Rossby wave (−  term), that tends to drive the offset towards negative values, and the jet speed (+ Ū term), that tends to increase the hot spot offset.For a given equilibrium temperature (and thus a given radiative timescale), a change in rotation period would change both terms in the imaginary part of the denominator.In our grid, as the rotation period increases from  norm = 0.37 to  norm = 1, the coriolis force decreases and   decreases.At the same time, the strength of the wind speed stays roughly constant, but the size of the equatorial jet increases (see Figure 3), so the term + Ū increases.Overall, this means that the eastward shift from the jet becomes increasingly larger than the westward shift from the Rossby waves, leading to an increase of the eastward shift with increasing period.
When the normalised period goes from  norm =1 to 2.23, both the jet speed and the Rossby wave speed reduce significantly.Therefore the divergent circulation (Hammond & Lewis 2021;Lewis & Hammond 2022), which by definition has no offset, becomes more and more dominant compared to the rotational circulation.This leads to an offset that decreases with increasing period.
In summary we find that for a given equilibrium temperature, at short orbital periods the circulation is dominated by the rotational component, which, because the Rossby wave speed decreases with increasing period while the jet becomes wider, leads to an offset that increases with period.Then, for larger rotation periods, the divergent circulation becomes more and more dominant, leading to an offset that decreases with increasing period.
To illustrate this behaviour, the thermal structure at the photosphere for our models at differing periods is plotted in Figure 10.Taking the 1800K row, in the leftmost plot we can clearly see the jet Doppler-shifting the Rossby wave towards the evening terminator, with the divergent circulation present at higher latitudes.In the centre column, the Doppler shift has increased, shifting the overall hot-spot further eastward.In the rightmost column, the magnitude of the rotation has now decreased to the point where the divergent circulation is now pulling the hot-spot back towards 0 • .Figure 9 shows that the maximum phase offset occurs at different orbital periods for different equilibrium temperatures.We can explain this by estimating that the peak phase shift will occur when Ū is larger than   in Equation 14.We estimate the ratio of these two quantities as the thermal Rossby number (Mitchell & Vallis 2010), which at the equator is where Ω is the rotation rate,  is the characteristic length scale (for which we use the planetary radius), and  is the characteristic velocity scale (for which we use the jet speed).
The right panel of Figure 9 displays the phase curve offset plotted against   .For those sets of models which display turning points, the maxima of the fitted curves are all found between values of   ∼ 0.5 − 1, supporting the argument that the rotational phase shift depends on the balance of the jet speed and Coriolis parameter.This trend breaks down for lower equilibrium temperatures.However, as explained in the previous section, these planets have a very uniform brightness temperature due to their long radiative timescale, and so their phase offset is determined by small variations in the thermal structure.
It should be noted that although all of the offset values calculated for our model grid are positive, negative values of the phase curve offset have been previously observed in certain wavelength bands (Esteves et al. 2015;Dang et al. 2018).The mechanisms causing westwards offsets are currently unclear, although MHD has been suggested as a potential reason (Rogers & Komacek 2014;Hindle et al. 2021).

Combined parameter trends: offset
To investigate the combined effect of the parameters we perform the same calculation from Section 4.2.4 for the phase curve offset.The mean and standard deviation for each parameter Δ is shown in Figure 11.The offset, as expected, has a much more complex dependency on the parameters.Increasing both surface gravity and metallicity reduce the offset across the board, with the metallicity having the larger impact of the two.Whereas increasing period first increases the offset then reduces it, as explained in the previous section.However, the standard deviation in the impact of period is large, indicating that cross parameters can drastically affect this relationship.
Phase curve offsets are most sensitive to gravity, metallicity and periods in the 1000-1600K range.The same trends can be seen in the models containing TiO/VO, although notably here the offset is much less significantly affected by the parameters, due to the substellar temperature inversion having a dominating effect on the hotspot position.

Amplitude and offset are set by different mechanisms
The relative amplitude is a measure of the maximum brightness contrast between hemispheres, and so is intrinsically linked to the heat redistribution which we discuss in Section 4.
The less efficient that day-to-night heat redistribution is, the larger the phase curve amplitude is expected to be.Because of this, it follows the same trends as seen in Figure 6 for the  -factor: increasing with temperature (Komacek & Showman 2016;Pluriel et al. 2020), metallicity and surface gravity, but decreasing with orbital period.A figure displaying this behaviour can be found in Appendix B.
Energy balance models based on advection of heat by a jet would lead to the conclusion that phase curve offset and amplitude should correlate with each other, with high winds expected to produce efficient heat redistribution and a large phase shift.However, as seen in Figure 12, our models do not see a systematic correlation between phase curve amplitude and phase curve offset This highlights that phase curve offset and phase curve amplitude are not set by the same component of the atmospheric circulation.The amplitude of the phase curve increases when heat redistribution becomes less efficient.As both the jet and the day-to-night flow con-0.30.4 0.5 0.6 0.7 0.8 0.9 Phase Curve Amplitude tribute to the heat transport, the phase curve amplitude is sensitive to the sum of the rotational and divergent components of the circulation.On the contrary, rotational and divergent parts of the circulation affect the phase curve offset in opposite directions.The rotational part tends to increase the offset whereas the divergent part tends to reduce it.As a consequence, the offset depends on the relative strength of divergent and rotational flow whereas the amplitude depends on the combined strength of divergent and rotational circulation.
This means that the relationship between phase curve offset and amplitude should not necessarily be linear, which we can see in Figure 13 where there are two clear regimes in the amplitude vs offset trend based on the thermal Rossby number (discussed in Section 4.3.1.For planets with short periods (  ≤ 1) the amplitude increases when offset decreases, For those with long periods (  > 1) it is the opposite.We do see that a linear relationship holds for fast rotating planets, the low (  ≤ 1) branch in Figure 13, with offset=−103.5 • ×amplitude+99.6 • .However, this relationship does not hold for larger rotation periods (those with   > 1).
The models deviating from the linear trend are those with  norm = 1 (planets around ∼G0 stars, Table 2) between  eq = 1000−1600 and those with  norm = 2.23 (planets around ∼F6 stars, Table 2) between  eq = 1000 − 2200.However, we expect the linear relationship between amplitude and offset to hold for planets around ∼K2 dwarfs.
There are, however, some processes our models omit which could alter the relationship between phase curve offset and amplitude.Our model does not include H2 dissociation, which can have a nonnegligible impact on the phase curve offset and day-night temperature contrast below 2400K (Tan & Komacek 2019).H2 dissociation reduces the day-night contrast and resulting phase curve amplitude, by transporting additional heat from the day to the night-side (Bell & Cowan 2018;Komacek & Tan 2018;Roth et al. 2021).This can reduce the speed of winds at a constant rotation period and reduce the phase-offset.Including this additional heating source could change the non-monotonic trend found between offset and rotation period to a linearly increasing trend for high temperature models.However, the effects of the H2 dissociation are greatly lessened in the presence of TiO/VO, which are likely to be present in the atmosphere at these high temperatures.Additionally, Beltz et al. (2022); Rogers (2017); Rogers & Komacek (2014) show that MHD can affect hot Jupiter atmospheres from ∼1500K.However, these effects do not become very strong until ∼2000K when there is sufficient coupling between the atmospheric magnetic field and circulation.As the magnetic drag timescale decreases the flow structure can be altered, initially providing pure drag then moving to oscillatory mean flows and eventually stationary westward mean flows.This could cause time variable winds and result in variable, or even westward, phase curve offsets.For the high temperature models, this would likely reduce the phase curve offsets, potentially altering the location of the turning point.Our models also do not include any clouds, which have been shown capable of drastically changing the observable features of hot Jupiter atmospheres Parmentier et al. (2021).The presence of clouds, particularly if mostly contained on the night-side, can increase the phase curve amplitude and decrease phase curve offset.

Temperature Contrast Between the Limbs
During transit spectroscopy one probes the limb of the planet.Limbto-limb variation in temperature can lead to variations in chemistry and cloud coverage.All of these lead to potential bias when interpreting transmission spectra (Line & Parmentier 2016;MacDonald & Madhusudhan 2017;Pluriel et al. 2020) and can have specific signatures in the transit light curve both at low spectral resolution (Line & Parmentier 2016) and at high spectral resolution (Ehrenreich et al. 2020;Wardenier et al. 2021).
Figure 14 shows the temperature contrast between the east and west limbs of our models at the photospheric pressure level, corrected for the transit chord optical depth by the geometric factor √︁ 2  / (Fortney 2005).The temperature is calculated as the average over all latitudes and ±20 • of longitude, in order to take into account the lineof-sight effect and the planet rotation during transit (Wardenier et al. 2022).A positive value here denotes an evening limb that is hotter than the morning limb.We can see that the temperature contrast does become increasingly large, although the upwards trend begins to flatten out at higher equilibrium temperatures.This is due to the combined affects of both the heat redistribution and hotspot offset.In Section 4.2 we see that heat redistribution becomes increasingly less efficient at higher equilibrium temperatures, which will naturally lead to a hotter evening limb and colder morning limb.Then, in Section 4.3 we see that the phase curve offset, which tracks the hotspot offset in the case of a cloudless atmosphere, reduces with equilibrium temperature.This will naturally lead to a more symmetrical structure and thus a smaller east/west contrast.
As the phase offset peaks for models with   ∼ 1, models with   < 1 or > 1 will have lower phase offset and therefore more symmetrical atmospheres.At low equilibrium temperatures ( eq <1600K) the heat redistribution efficiency is strong and most models have   > 1, so the limb temperature contrast is small.As the equilibrium temperature increases, the efficiency of heat redistribution reduces, but the offset also reduces.In our models we find the limb temperature difference starts to increase.Implying that the reduction to heat redistribution with equilibrium temperature must be more impactful on the limb temperatures in this regime than the offset reducing with equilibrium temperature.This behaviour is maintained between 1000-2000k.Above 2000K, the heat redistribution trend with temperature flattens out as it reaches close to its minimum (e.g large  -factor).However, the offset is still reducing with  eq , so the limb temperature difference stops increasing for models with   ∼ 1 and even decreases for those with   < 1.
In the presence of TiO/VO, the heat redistribution is much worse at lower temperatures.The offset is also significantly reduced.From the right panel of Figure 14, we can see that in this case the relationship between the limb temperature contrast appears to flatten at much lower temperatures, which is likely due to the decrease in variation for both the heat redistribution and offset when TiO/VO are present.

COMPARISON WITH OBSERVATIONS
We now look at two different types of population level eclipse observations: the population of planets observed with the Spitzer Space Telescope and the one observed with the Wide Field Camera 3 on board the Hubble Space Telescope and compare them to our population of modelled planets.

Spitzer Band Secondary Eclipse
The band-averaged measurements provide the brightness temperature of the planet day-side in a given band-pass, which can be used to measure a "band averaged" heat redistribution parameter.Observations of exoplanet atmospheres do not generally probe the bolometric quantities described above, but band averaged or spectrally resolved properties of the atmosphere.Because different band-passes probe different atmospheric depth, the band-averaged quantities can be different from the bolometric ones.Thus we will now look at the model predictions for the Spitzer eclipse measurements.
Over its lifetime Spitzer observed 457 eclipses of 122 hot Jupiters in it's 3.6 and 4.5 micron observing modes (Deming et al. 2023).Eclipses from both channels spectra, where CH4 and CO features are expected to be present depending on the planets equilibrium temperature, have been re-analysed in homogeneous ways by several authors (Baxter et al. 2020;Garhart et al. 2020;Deming et al. 2023).The distribution of planet parameters observed matches well with the overall transiting planet population, as seen in Figure 1, making it a good representation for diversity in population level studies.Baxter et al. (2020) find that there is a transition between seeing CO in absorption to CO in emission at T eq = 1600 ± 100K, and are the first to attribute this switch to a transition between non-inverted and inverted temperature profiles.We hereafter use the recent re-analyses of Deming et al. (2023), because it is the most recent and comprehensive one.Deming et al. (2023) identify a statistically significant unexpected and sudden rise in the brightness temperatures of the population, with  the brightness temperature increasing by 291±41 K between equilibrium temperatures of 1714-1818 K. Deming et al. (2023) propose a few explanation to explain this sudden rise in brightness temperature, such as the onset of magnetic drag which would inhibit longitudinal heat redistribution or the rapid dissipation of day-side clouds at these temperatures.
Here we propose a novel mechanism to explain this rise in temperature: that the sudden appearance of TiO and VO changes the ability of the atmosphere to redistribute heat.Indeed, as shown in Figure 5, the presence of TiO and VO in a planetary atmosphere reduces the heat transport and thus increases the heat redistribution factor (Δ  avg = 0.216).We will now look at the specific effect of TiO/VO on the brightness temperature of the Spitzer band-pass.
For this, we calculate the brightness temperature in the Spitzer band-pass in our grid of models by first fitting a blackbody to the spectral region of the Spitzer band-passes, we can then use the temperature of this blackbody as the brightness temperature and compare with the equilibrium temperature of our models.
We show the difference between day-side brightness temperature and equilibrium temperature in Figure 15.This difference would be zero for a planet with full heat redistribution and would be (2.66 1/4 − 1) *  eq for a planet with no heat redistribution.In between these two extremes, the band averaged redistribution value depends on both the day/night heat transport and on the spatially varying chemical composition of the planet, thus sometimes leading to negative values (Dobbs-Dixon & Cowan 2017).Figure 15 shows two tracks for the brightness temperature versus equilibrium temperature, corresponding to models in chemical equilibrium and to models where TiO/VO have been artificially removed from the atmosphere.Both tracks have a brightness temperature that increases relative to the equilibrium temperature as the equilibrium temperature gets larger.This is a direct consequence of the heat redistribution parameter that increases with equilibrium temperature, as shown in Figure 5.
Additionally, we can see that the with and without TiO/VO track are offset by, approximately, 50-300K depending on the equilibrium temperature.This was discussed in Section 4.2 and is due to TiO/VO intercepting and re-radiating part of the stellar flux at very low pressures, where the radiative timescales are very short.Overall, this reduces the heat transport to the planetary night-side and increases the planet brightness temperature.
Figure 15 also compares our population of models with the best fit to the population of planets proposed by Deming et al. (2023), which consist of a rising slope at low equilibrium temperature, a jump around 1714-1818 K and a gentle slope at high equilibrium temperature.Without any adjustment, we can see that the best fit to the observed population of planets falls right inside our population of models, with the low temperature fit corresponding to the no TiO/VO track and the high temperature fit corresponding to the TiO/VO track.We do note, however, that towards very high temperatures the Deming et al. (2023) fit lies slightly below the TiO/VO track.Given that the parameter distribution of the Spitzer observations is representative of the wider hot Jupiter population (see Figure 1 for this comparison), and not weighted towards planets with parameters that would give lower values here, one possible cause is partial depletion of TiO/VO in these atmospheres, which we discuss with more detail in Section 5.2.
We now show in 16 the brightness temperature difference between models with and without TiO/VO, with all other parameters being constant (e.g.size of the "jump" in Figure 15 between equivalent models).We see that the position and amplitude of the jump observed by Deming et al. (2023) (black cross in Figure 16) aligns very well with the model predictions if the population of planets jumps from the no TiO/VO track to the TiO/VO track around 1800K.We therefore conclude that the appearance of TiO/VO between 1714 K and 1818 K can explain the jump in brightness temperature observed in the population of hot Jupiters.This would align well with observational evidence suggesting the appearance of gaseous TiO/VO in atmospheres at around this temperature range.Parmentier et al. (2016), for instance, predict the transition to occur at ≤1900K based on the large apparent albedos for Kepler-76b and HAT-P-7b, and better fitting of Kepler light-curves for planets above this temperature.This transitional temperature is consistent with cold trap models (Beatty et al. 2017;Parmentier et al. 2016Parmentier et al. , 2013) ) and the temperature predicted for deep cold traps to possibly occur from micro-physical cloud models (Powell et al. 2018).
We further note that, although our modelled planets have a large scatter along the temperature trend ( =∼ 104K, which increases to  =∼ 166K when the error bars from the observed Spitzer population are included), this scatter is significantly smaller than the one seen in the population of Spitzer data points ( ∼ 221K).We conclude that either additional mechanisms lead to an additional scatter around the mean value compared to our expectation (e.g. chemical disequilibrium, presence of night-side clouds, atmospheric drag), or that the Spitzer IRAC instruments have additional noise that increases the scatter of the population (Hansen et al. 2014).
We now wonder whether this jump should be seen with other instruments, such as HST/WFC3.We show in Figure 17 the difference between the brightness temperature minus the equilibrium temperature for a single model at 1800K with and without TiO/VO as a function of wavelength.We can see that the size of the jump is highly dependent on the wavelength.This is because different wavelengths probe different depths in the atmosphere.Because TiO/VO changes both the heat redistribution and the vertical structure of the thermal profile, at some wavelengths (or some pressure levels) the reduction of the heat redistribution due to the presence of TiO/VO is compensated by the cooling of the atmosphere below the thermal inversion due to the anti-greenhouse effect (Parmentier & Guillot 2014).This is particularly prevalent in the WFC3 band-pass, where no specific jump in the population would be expected.
We would therefore predict that if a similar population study using HST/WFC3 observations were to be carried out, then the brightness temperature increase would be absent altogether.However, over a wavelength range such as that probed by the JWST/NIRSPEC/G395 instrument, 2.87-5.14, the signature of this transition would be larger.

Spectral Feature Strength
The thermal flux emitted by a planet's day-side can be measured by observing its eclipse behind the star.When the spectroscopic eclipse is measured, the spectra can provide insights into the thermal profile of the planet.In particular, the spectra can be used to determine whether the thermal profile is increasing with decreasing pressure (an inverted atmosphere that produces emission features) or whether it is decreasing with decreasing pressure (a non-inverted atmosphere with absorption features).Therefore, another way to probe the TiO/VO to no TiO/VO transition at the population level is to look at the spectral variation of the water feature within the HST/WFC3 band-pass, where water has a strong molecular feature between 1.35-1.48.When this feature transitions from absorption to emission, it indicates that a temperature inversion is present in the atmosphere.For hotter planets the combination of thermal dissociation and the appearance of H-opacities reduces the spectral feature of water and brings the spectrum close to a blackbody.If the transition from TiO/VO to no TiO/VO happens at cooler temperatures than the onset of dissociation, we would expect a sharp transition between planets with inverted thermal profile and thus emission features (such as Cont et al. (2022); Evans et al. (2016); Coulombe et al. (2023)) to planets with non-inverted thermal profile and thus absorption features (such as Line et al. (2016); Kreidberg et al. (2014)).
The Hubble Space Telescope observed the eclipse of 21 planets with the WFC3 instruments.A spectral index was determined by Mansfield et al. (2021) in order to reduce these points (Initially 19 in Mansfield et al. (2021), with later additions in Mansfield et al. (2022); Fu et al. (2022)) into a single index value.The WFC3 index measured the departure of the 1.4 micron feature (usually due to water) to the underlying continuum.To calculate it we first fit a blackbody to two 'out-of-band' regions between 1.22-1.33and 1.53-1.61,where water absorption in the band-pass is minimal.The temperature of this blackbody is then referred to as the observed day-side temperature, T day .In Mansfield et al. (2021), this temperature is plotted against the water feature strength.However, for ease of comparison to Section 5.1, we use the equilibrium temperature here.The water feature strength in the 1.35-1.48'in-band' region is then calculated by (Mansfield et al. 2021): Based on this definition, positive values correspond to a water feature observed in absorption and negative values to a water feature observed in emission.
In Figure 18 the water feature strength for all models within the grid can be seen when compared to observations from HST/WFC3 seen in Mansfield et al. (2021Mansfield et al. ( , 2022)); Fu et al. (2022).In our models we see that, when TiO/VO are absent from the models, the water feature is seen in absorption (negative index value) and the strength of the water feature decreases with increasing temperature as the dayside thermal profile becomes more isothermal for more irradiated planets.When TiO/VO are present, the water feature is initially seen in absorption below ∼ 1500K, where the TiO/VO begin to condense out of the atmosphere in chemical equilibrium.The feature rapidly swaps to emission above this temperature, with the strength of emission decreasing towards higher temperature models.Both models containing TiO/VO and those without are seen to converge towards the zero line at high day-side temperatures.This occurs as increasing quantities of water are thermally dissociated in the atmosphere at high temperatures and H-opacities start filling the gap in the water bands (Parmentier et al. 2021).
Figure 18 further shows that our models are able to encompass all the observed data-points.Particularly, we see that the spread of the water feature strength index becomes smaller at high equilibrium temperature, where the dissociation of water and the H-opacity reduce the water feature strength in all planets.
Looking in more details, we see that the observations seem to lie in between the TiO/VO and no TiO/VO model tracks for the 1400-2000K equilibrium temperature range.Particularly, we do not see the large emission features that should be prominent for planets in the 1600K-1800K range.Similarly, we do not see the non-inverted water features expected by the no TiO/VO model in that specific range of temperatures.
One possible explanation is that TiO/VO is indeed present in these atmospheres, but at a smaller abundance than in the chemical equilibrium prediction.This could happen due to deep or day/night partial cold trap mechanisms (Spiegel et al. 2009;Parmentier et al. 2016;Beatty et al. 2017;Pelletier et al. 2023), which are omitted from our 3D models, or to original elemental Ti and V depletion.Pelletier et al. (2023), for instance, recently detected VO but no TiO in WASP-76b ( eq = 2228K), and suggest condensation of TiO to be responsible.Such a partial depletion of TiO/VO would naturally bring the models between the TiO/VO and no TiO/VO track in Figure 18, allowing a closer match to the data.
In Figure 18, we also highlight the transition temperature between TiO/VO and no TiO/VO expected from the Spitzer data in Section 5.1.The current HST data does not show a steep difference in water feature strength around this temperature, however, the step might be too small to be measured given the error bars on the data.
We note that the partial reduction of TiO/VO that we argue can explain the HST data set does not necessarily rule out our explanation for the sharp increase in brightness temperature in the Spitzer data set (see Section 5.1).Indeed, the TiO/VO model track in Figure 15 lies above the observations and partial depletion of TiO/VO could easily lower the model track and lead to a better agreement with the observations.To determine if this works quantitatively would require models with partial TiO/VO abundances, which are out of the scope of the current paper.Nonetheless, we believe that the current Spitzer data provides compelling argument for future ground-based optical hi-res observations in order to determine the TiO and VO abundances in these atmospheres.
Another possible explanation for the damped features in the HST bandpass would be that the H-abundance is larger than expected by our models.This could be due to an increase in alkali abundances, which are the prime source of electrons for the H-.As for the partial TiO/VO abundance scenario, we leave this for a future study.

Spitzer Phase Curve offset and amplitude
Spitzer observed 21 phase curves of hot Jupiters in its 3.6  and 4.5  bands.For those planets within the equilibrium temperature bounds of our grid, we take the values of phase curve offset and amplitudes from Parmentier & Crossfield (2018) with a correction for HAT-P-7b found in Parmentier et al. (2021).For phase curve data published after this we use the values summarised in May et al. (2022) for Qatar-1b, Qatar-2b, WASP-34b, WASP-52b and WASP-104B, Beatty et al. (2019) for KELT-1b and Bell et al. (2021) for KELT-16b and MASCARA-1b.For CoRoT-2b, we take the offset  2023), and points not containing TiO/VO above the the transition temperature, have been reduced to emphasise the transition.The blue points are models without TiO/VO, the red points are the models with TiO/VO and the black points with error-bars are the observational data from Mansfield et al. (2021Mansfield et al. ( , 2022)); Fu et al. (2022).The yellow shaded region displays the temperature range of the abrupt increase in emergent flux from Deming et al. (2023) value from Dang et al. (2018), but use the re-analysis of Bell et al. (2021) to calculate an amplitude.We also exclude HD149458b from our data set, as done in Zhang et al. (2018).For WASP-43b we use the most resent re-analysis of May & Stevenson (2020).
Figure 19 displays the phase curve offsets and amplitudes calculated from the model grid compared to the Spitzer/IRAC 4.5 band observations.The 3.6  comparison can be found in Appendix C. The model grid shows the amplitude to increase with the equilibrium temperature and decrease with orbital period.Overall the trend with equilibrium temperature matches the observations quite well, with a few outliers at lower temperatures.These outliers (Qatar-2b, WASP-34b May et al. (2022), CoRoT-2b Bell et al. (2021) and WASP-43b May & Stevenson (2020)), with amplitudes close to unity, are likely due to night-side clouds, which are omitted from our models and would significantly reduce the night-side brightness temperature and increase observed amplitude.In particular, for WASP-43b night-side clouds have been postulated to be the reason for the observed amplitude (Kataria et al. 2015;Venot et al. 2020;Parmentier et al. 2021).One feature that seem to be present in the data and not in the models is the sudden jump to lower amplitudes for planets with orbital periods larger than 1 day.However we do not have an explanation for this potential jump.
As our model grid omits all processes which could cause negative phase curve offsets, we will focus on comparison to only ones with positive values.The trend in offset with equilibrium temperature is remarkably consistent between the models and observations, decreasing for higher equilibrium temperatures.The observational data also matches the model grid well when compared by period.In our models we clearly see the turning point behaviour detailed in Section 4.3.1.However, there is no observed data beyond periods of ∼4 days for which the offset would begin to noticeably decrease if following the model grid trend.
Overall, from this comparison we can see that the spread of values observed for both phase curve amplitude and (positive) offset can be induced by the varying planetary parameters within our model grid.

CONCLUSIONS
In this paper, we have created the largest library (345 simulations) of 3D hot Jupiter models to date.These models were produced using the non-grey SPARC/MITgcm global circulation model (Showman et al. 2009), under the assumption of local chemical equilibrium with local rain-out of condensate materials (Visscher et al. 2006(Visscher et al. , 2010)).Phase curves and spectra were further calculated via the plane-parallel radiative transfer code of (Marley & McKay 1999), at increased spectral resolution.The parameters of equilibrium temperature, orbital period, metallicity, surface gravity were all varied throughout the grid, including cross-variation.In addition, models with the strong photo-absorbing molecules TiO and VO artificially removed from chemical equilibrium were also run.These intrinsic planetary prop- erties were all varied systematically and jointly over the parameters space.
We show that the day-night heat redistribution in hot Jupiter atmospheres is as sensitive to the combined effects of atmospheric metallicity, rotation period and surface gravity as it is to the equilibrium temperature.We therefore expect a large scatter along the general trend of increasing heat redistribution with temperature.As a consequence, a large sample of planets will be needed to measure any trends in the hot Jupiter population by observation, unless one can control some of these parameters through careful target selection.
We find that the relationship between phase curve offset and orbital period is non-monotonic within our models.With an offset that gen-erally increases with rotation period for slowly rotating planets and decreases with rotation period for fast rotating planets.We propose that the turning point is due to the balance between the rotational circulation (i.e. the jet) and the overturning circulation, and can be characterised by the thermal Rossby number of the atmosphere (see Equation 15).This relationship breaks down for low equilibrium temperature models which are in a different circulation regime than the hotter planets.
We further find that phase curve offset and phase curve amplitude do not necessarily track each other.This is because the rotational and divergent part of the atmosphere both advect heat from the day to the night-side of the planet, but have counteracting effect on the hot spot offset.This competition between the rotational and divergent effects create two clear regimes within our models.For those with low periods (  ≤ 1) the amplitude increases when offset decreases, and those with long periods (  > 1) the opposite is found.This produces a scatter in the relationship between the offset and amplitude.
We find that the combined effects of the heat redistribution and phase shift also impacts the temperature contrast between the planetary limbs.As heat redistribution becomes increasingly less efficient at higher temperatures, the evening limb becomes increasingly hotter than the morning limb, increasing the east/west temperature contrast.For the hottest planets, however, the reduced hot spot shift reduces the limb temperature asymmetries.We predict the largest limb-to-limb temperature differences for planets with  eq ≈ 2000.
By comparing our spectra to the re-analysed Spitzer/IRAC secondary eclipse data from Deming et al. (2023), we propose a novel mechanism that explains quantitatively the abrupt increase in the brightness temperature seen in this population: the appearance of TiO/VO in the atmosphere at  eq ≈ 1800 reduces the efficiency of heat transport and creates a sudden brightening of the flux in the Spitzer band-pass.
Further comparing our models to the spectral variation of the water feature strength observed within the HST/WFC3 band-pass (Mansfield et al. 2021), we find that most observation between 1400-2200K are poorly matched by our grid of models.We suggest that this could be due to the partial depletion of TiO/VO due to non local condensation processes.We leave for future work whether an atmosphere partially depleted of TiO/VO can indeed fit both Hubble and Spitzer population level trends.
Finally, by looking at observed phase curve parameters, we find that the intrinsic scatter of our model grid aligns well with the observed spread in the data.We also find that the trend that phase curve amplitude increases with equilibrium temperature is well matched by our models.However, we do not predict the possible sudden jump towards lower phase curve amplitudes for planets with period larger than 1 day.We also find that the offset trends in both equilibrium temperature and period are well matched with our models, although our model is unable to reproduce any of the observed negative offsets.
Overall, this GCM library is a powerful tool that can be used to investigate the the effects of individual parameters, and a combination of parameters, on key observable features.There are many further avenues to explore using these models, including a deeper look into the dynamics of hot Jupiter atmospheres.Not only this but the parameter space spanned by the grid lends it to being used as model data base in the future, from which fast interpolated PT structure and spectra can be calculated.factor has negligible change after the 150 day mark, suggesting that we have in fact reached a pseudo steady-state at the photospheric pressure level in these models.The dashed lines show the value after 900 days and the solid lines the value at the integration time within the grid (150 days for models with equilibrium temperatures above 1400K and 300 days for models below).The points marked with a cross show the redistribution factor of the initial profile.

Figure 1 .
Figure 1.Distribution of known exoplanets (0.3<  <3  ) with equilibrium temperature vs orbital period.Planets that fall within our grid boundaries are color-coded by their gravity.Planets that fall outside our grid boundaries are left grey.The blue dashed lines on the scatter plot represent the minimum and maximum values of constant period multiplicity for our model grid, as defined in Section 3. The outlined solid blue shape shows the area that is covered by our model grid in the equilibrium temperature vs log(period) space.The histogram to the left of the scatter plot displays the distribution of planets in the equilibrium temperature space (blue bars) when compared to the distribution of the measured Spitzer population from Deming et al. (2023) (red bars).The histogram above the scatter plot shows the same population level comparison for the period distribution.

Figure 2 .
Figure2.Latitude vs longitude thermal structure, temperature normalised by equilibrium temperature, for all models within the grid, at the pressure level closest to the 1.4 micron photospheric pressure.Equilibrium temperature is sorted into columns and atmospheric metallicity into rows.The surface gravity increases between the values seen in Table1down each bracket of rows and the normalised orbital period increases in the same manor along columns.Models excluding TiO/VO are separated into the top half, with models including the two species in the bottom half.

PFigure 3 .
Figure3.Zonal-mean zonal wind speed for all models within the grid.Within each subplot, the x-axis shows latitude and the y-axis pressure.Equilibrium temperature is sorted into columns and atmospheric metallicity into rows.The surface gravity increases between the values seen in Table1down each bracket of rows and normalised orbital period increases in the same manor along columns.Models excluding TiO/VO are separated into the top half, with models including the two species in the bottom half.

Figure 5 .
Figure 5. Bolometric redistribution factor ,  -factor, a measure of the day-night atmospheric heat redistribution, for the model grid vs equilibrium temperature.Left Panel: Shows models with TiO/VO removed.Right Panel: Shows models containing TiO/VO.The shape of each point depicts the orbital period, and the color depicts the metallicity of the model.The grey dashed lines are drawn at key values corresponding to specific redistribution regimes.

Figure 7 .
Figure 7. Relative impact of each planetary parameter on heat redistribution.The y-axis of each plot is the delta in the observable feature divided by the delta in the log of each planetary parameter.The x-axis is the Equilibrium temperature.Light blue denotes the change due to metallicity, dark blue the change due to period and red the change due to surface gravity.The error bars display the standard deviation in the impact of each parameter.Left Panel: Models with TiO/VO removed.Right Panel: Models containing TiO/VO.

Figure 9 .Figure 10 .
Figure9.Left Panel: Phase curve offset for all models plotted against log(period).Right Panel: Phase curve offset for all models plotted against log(  ).In both plots, models with a larger offset at P Norm =1 than at P Norm =0.37 or 2.23 are color coded in equilibrium temperature.The black crosses display the maxima of the curve fitted to these models.Any models that do not satisfy this criteria are in grey.Models with TiO and VO are shown in dashed lines and those without in filled lines.

Figure 11 .Figure 12 .
Figure 11.Relative impact of each planetary parameter on phase curve offset.The y-axis of each plot is the delta in the observable feature divided by the delta in the log of each planetary parameter.The x-axis is the Equilibrium temperature.Light blue denotes the change due to metallicity, dark blue the change due to period and red the change due to surface gravity.The error bars display the standard deviation in the impact of each parameter.Left Panel: Models with TiO/VO removed.Right Panel: Models containing TiO/VO

Figure 13 .
Figure 13.Phase curve offset plotted against amplitude for all models within the grid.Colour denotes the thermal Rossby number of the model and shape denotes the normalised period.The red dashed line displays the best fit to the models with value of log(  )≤-0.1.

Figure 14 .Figure 15 .
Figure 14.Temperature difference between the evening and morning limb pressure-temperature profiles at the transit photospheric pressure level plotted against equilibrium temperature.Left Panel: Without TiO/VO.Right Panel: With TiO/VO.

Figure 16 .
Figure16.Difference between the day-side brightness temperature for models containing and not containing TiO/VO.The yellow shaded region displays the temperature range of the abrupt increase in emergent flux fromDeming et al. (2023).The black cross shows the exact value found inDeming et al. (2023).Left Panel: Spitzer 1 (3.6 microns).Right Panel: Spitzer 2 (4.5 microns).

Figure 17 .
Figure 17.This figure shows   TiO/VO −   noTiO/VO for an 1800K model over all wavelengths.The red, light blue and dark blue shaded regions show the HST/WFC3, Spitzer 1 and Spitzer 2 band-passes respectively.Dashed lines of the same colors represent the average value within each band-pass.

Figure 18 .
Figure18.figure shows a comparison between the water feature strength calculated for the 3D GCM grid and observational data.The opacity of points containing TiO/VO with equilibrium temperatures below the transition temperature found inDeming et al. (2023), and points not containing TiO/VO above the the transition temperature, have been reduced to emphasise the transition.The blue points are models without TiO/VO, the red points are the models with TiO/VO and the black points with error-bars are the observational data fromMansfield et al. (2021Mansfield et al. ( , 2022));Fu et al. (2022).The yellow shaded region displays the temperature range of the abrupt increase in emergent flux fromDeming et al. (2023)

Figure 19 .
Figure19.This figure shows a comparison between the phase curve offset and amplitudes in the Spitzer/IRAC 4.5 band calculated for the 3D GCM grid and the observational data.The blue points are models without TiO/VO, the red points are the models with TiO/VO and the black points with error-bars are the observational data.Left Panels: shows the offset (lower) and amplitude (upper) vs equilibrium temperature.Right Panels: shows the offset (lower) and amplitude (upper) vs orbital period.

Figure C1 .
Figure C1.This figure shows a comparison between the phase curve offset and amplitudes in the Spitzer/IRAC 3.6 band calculated for the 3D GCM grid and the observational data.The blue points are models without TiO/VO, the red points are the models with TiO/VO and the black points with error-bars are the observed data.Left Panels: shows the offset (lower) and amplitude (upper) vs equilibrium temperature.Right Panels: shows the offset (lower) and amplitude (upper) vs orbital period.

Figure E1 .
Figure E1.Redistribution factor as a function of model run time for a range of equilibrium temperatures in the model grid.The dashed lines show the value after 900 days and the solid lines the value at the integration time within the grid (150 days for models with equilibrium temperatures above 1400K and 300 days for models below).The points marked with a cross show the redistribution factor of the initial profile.

Table 1 .
Free parameters used in the 3D GCM grid.

Table 2 .
(Hauschildt et al. 1999 stellar models, we use the NextGen stellar spectra(Hauschildt et al. 1999) for these stars.
star [M sun ]  star [R sun ]  star [K] normalised Period Multiplier Absolute Period Range [days] (Coulombe et al. 2023), and VO was recently detected in WASP-76b Calculated bolometric  -factor values for a subset of models within the grid.Each panel displays the variation in  -factor with  eq due to one of: metallicity, orbital period and surface gravity, whilst the other two parameters are kept constant (when constant Log(M/H)=0.0,P Norm =1, log(g)=1.3).Crossterms between parameters are not included.Models excluding and including TiO/VO are displayed with solid and dashed lines respectively.
This timescale can range from a few seconds to thousands of years at deeper pressure levels (Showman & Guillot Left Panel: Variation in  -factor with metallicity.Middle Panel: Variation in  -factor with different normalised orbital period.These tracks corresponds to planets tidally-locked around different stellar types.Right Panel: Variation in  -factor with surface gravity. Phase curve offset for a subset of models within the grid.Each panel displays the variation in offset with  eq due to one of: orbital period, surface gravity and metallicity, whilst the other two parameters are kept constant (when constant Log(M/H)=0.0,P Norm =1, log(g)=1.3).Cross-terms between parameters are not included.Left Panel: Variation in offset whilst changing metallicity.Middle Panel: Variation in offset with different normalised orbital period.These tracks corresponds to planets tidally-locked around different stellar types.Right Panel: Variation in offset whilst changing surface gravity.