Pseudo-2D Modelling of Heat Redistribution Through H$_2$ Thermal Dissociation/Recombination: Consequences for Ultra-Hot Jupiters

Thermal dissociation and recombination of molecular hydrogen, H_2, in the atmospheres of ultra-hot Jupiters (UHJs) has been shown to play an important role in global heat redistribution. This, in turn, significantly impacts their planetary emission, yet only limited investigations on the atmospheric effects have so far been conducted. Here we investigate the heat redistribution caused by this dissociation/recombination reaction, alongside feedback mechanisms between the atmospheric chemistry and radiative transfer, for a planetary and stellar configuration typical of UHJs. To do this, we have developed a time-dependent pseudo-2D model, including a treatment of time-independent equilibrium chemical effects. As a result of the reaction heat redistribution, we find temperature changes of up to $\sim$400 K in the atmosphere. When TiO and VO are additionally considered as opacity sources, these changes in temperature increase to over $\sim$800 K in some areas. This heat redistribution is found to significantly shift the region of peak atmospheric temperature, or hotspot, towards the evening terminator in both cases. The impact of varying the longitudinal wind speed on the reaction heat distribution is also investigated. When excluding TiO/VO, increased wind speeds are shown to increase the impact of the reaction heat redistribution up to a threshold wind speed. When including TiO/VO there is no apparent wind speed threshold, due to thermal stabilisation by these species. We also construct pseudo-2D phase curves from our model, and highlight both significant spectral flux damping and increased phase offset caused by the reaction heat redistribution.


INTRODUCTION
Offering both valuable insight into extreme planetary atmospheric conditions, and an abundance of available data due to detection biases, the hydrogen-dominated atmospheres of highly irradiated ultra-hot Jupiter (UHJ) exoplanets have recently become a focus of both the modelling and observational astrophysics communities. These planets are broadly assumed to be tidally-locked, due to strong gravitational interactions with their host star, causing huge irradiative disparities between their permanently stellar-facing dayside and much colder night-side. Sophisticated multi-dimensional atmospheric models exist for these planets, and high speed variable winds have been predicted through modelling (Guillot & Showman 2002;Mayne et al. 2014) and inferred from observations (Snellen et al. 2010). For UHJs, such as WASP-12b, the day-side equilibrium temperature can reach ≥ 2500 K (Chakrabarty & Sengupta 2019), with night-side temperatures ≥ 1000 K colder (Cowan et al. 2012). This extreme temperature gradient between the two hemispheres drives global circulation in the form of particularly strong longitudinal winds, dominated by an equatorial jet, in some cases reaching speeds of over 5 km s −1 (Komacek et al. 2017). ★ alexander.roth@physics.ox.ac.uk For the range of typical atmospheric pressures, temperatures between 2000 -4000 K cause molecular hydrogen (H 2 ) to thermally dissociate into atomic hydrogen (H). Energy is absorbed during the endothermic dissociation reaction (H 2 → H + H) and released during the exothermic recombination reaction (H + H → H 2 ). Bell & Cowan (2018) have shown that dissociation on the day-side followed by recombination on the much cooler night-side, driven by the large temperature difference between the substellar and antistellar hemispheres, reduce the overall day-night temperature gradient. This heat redistribution effect has been demonstrated to impact photometric measurements for ultra-hot Jupiters, specifically transit photometry, both theoretically and observationally (Komacek & Tan 2018;Tan & Komacek 2019;Mansfield et al. 2020). The early work by Bell & Cowan (2018) is based on a simple atmospheric energy balance model. In this paper we further investigate this process using a more sophisticated atmosphere model that includes equilibrium chemistry, radiative transfer and a parameterised horizontal transport.
A more thorough investigation on the consequences caused by H 2 thermal dissociation/recombination must include the treatment of complex chemical feedback within the atmosphere. One-dimensional (1D) radiative transfer and chemistry models, which were built to investigate planetary temperature structures and the associated chemical composition by solving for radiative-convective equilibrium, are a typical starting point for this. Most of these models make the assumption of local chemical equilibrium (e.g., Mollière et al. 2015;Gandhi & Madhusudhan 2017), which is likely a good assumption for very hot atmospheres. Some models have also incorporated the ability to consistently calculate the temperature profile with nonequilibrium chemistry, including vertical mixing and photochemistry ). These models have since been used to investigate a range of planetary atmospheres, such as hot Neptunes (Agúndez et al. 2014b) and brown dwarfs (Phillips et al. 2020), but have predominantly been used for hot Jupiters (Moses et al. 2013;Goyal et al. 2020).
The use of 1D models for planets like ultra-hot Jupiters, which display such extreme longitudinal atmospheric change, clearly has significant limitations. These models, by design, average incoming stellar flux to a single value, severely limiting their use for any study on effects caused by longitudinal temperature or chemical variations. Three-dimensional (3D) global circulation models can be used to model more intricate planetary dynamics (e.g. Amundsen et al. 2016), and have been applied to the H 2 dissociation/recombination effect (Tan & Komacek 2019). However, due to their increased complexity and computational cost, handling of the expensive chemistry calculations is typically simplified compared with the approach in 1D models.
It therefore stands to reason that, in order to properly investigate the affects of H 2 thermal dissociation/recombination, some fundamental changes and additions must be made to existing models. Both pseudo-2D (Agúndez et al. 2014a) and time-dependent (Iro et al. 2005) models have previously been adapted from pre-existing 1D models and used to investigate hot Jupiters. In this work we further develop an existing 1D/2D atmosphere model (ATMO) to include a time-dependent pseudo-2D capability. We then use this to study the effect of H 2 thermal dissociation/recombination, and its importance within UHJ atmospheres. Among the most important effects previously unaccounted for, and a primary focus of this study, are the numerous interactions between H 2 thermal dissociation/recombination and important atmospheric opacity sources, particularly TiO and VO, which are known to drastically affect the pressure-temperature structure within these atmospheres (Fortney et al. 2008;Parmentier et al. 2018).
A brief review and details of new developments to the ATMO model are presented in Section 2. Our results and discussion on heat redistribution from the thermal dissociation/recombination of hydrogen can be found in Section 3. A conclusion is then drawn on all the findings in Section 4.

Atmosphere model: ATMO
The 1D/2D atmosphere code ATMO (Tremblin et al. 2015;Drummond et al. 2016;Tremblin et al. 2017;Goyal et al. 2018) has been widely used to model substellar atmospheres. Most typically the model is used in a 1D form, solving for the 1D radiative-convective equilibrium profile of an atmosphere. An extension of the model by Tremblin et al. (2017) included the ability to solve for the stationarystate 2D (longitude-altitude) atmosphere. ATMO has been applied to hydrogen-dominated, irradiated exoplanet atmospheres (e.g. Goyal et al. 2018;Drummond et al. 2019), brown dwarf atmospheres (e.g. Tremblin et al. 2016;Phillips et al. 2020) and as a retrieval code for constraining observations (e.g. Evans et al. 2017;Nikolov et al. 2018;Wakeford et al. 2018).  John (1988) In its 1D form, ATMO solves for hydrostatic balance and radiativeconvective equilibrium, with an internal heat flux (we use an intrinsic temperature of int = 100 K) and irradiation at the top of the atmosphere as boundary conditions. The radiative transfer equation is solved in 1D plane-parallel geometry and includes isotropic scattering. We include absorption by CH 4 , H 2 O, CO, CO 2 , NH 3 , Na, K, Li, Rb, Cs, TiO, VO, and Has opacity sources. H 2 -H 2 and H 2 -He collision-induced absorption is also included. All references for the opacity source data used in this study can be found in Table 1. The most up-to-date description of the line-lists and calculation of the opacities can be found in Goyal et al. (2020). The model uses the correlated-approximation with the method of random overlap to calculate the combined opacity of the mixture (Lacis & Oinas 1991;Amundsen et al. 2017), with 32 bands to compute the radiative flux in the radiative-convective equilibrium iterations.
Chemical equilibrium abundances are obtained by minimising the Gibbs free energy, following the method of Gordon & McBride (1994). The thermodynamic properties of all chemical species are expressed in terms of NASA polynomials, using coefficients from McBride et al. (1993McBride et al. ( , 2002. The Gibbs minimisation method allows for the depletion of gas phase species due to condensation, but we only consider a gas-phase composition here. This is likely to be a good assumption for the results presented in this paper, since the high-temperatures of the atmospheres that we focus on in this study mean that condensation is expected to be unimportant. The pressure and temperature dependence of chemical equilibrium abundances creates the requirement for consistency between the radiative-convective and chemistry calculations. As the pressuretemperature ( -) structure evolves, chemical abundances shift to updated values which in turn changes the opacity. This in turn affects the radiative transfer and therefore the -structure. By periodically cycling between the radiative-convective solver and the Gibbs energy minimisation scheme, ATMO can account for this feedback mechanism and determine -and chemical profiles that satisfy radiative-convective equilibrium.
ATMO also includes the option to model time-dependent nonequilibrium chemical effects, such as vertical mixing and photochemistry. However, for this study, we only include time-independent chemical equilibrium processes.

A time-dependent approach to solving for 1D radiative equilibrium
Most models, including previous applications of ATMO, use a timeindependent flux-balancing approach when solving for radiativeconvective equilibrium. We refer to this approach as the 'timeindependent' approach as it involves applying incremental perturbations to the -profile until a profile that satisfies radiativeconvective balance is achieved. An alternative approach is to calculate heating and cooling rates and increment the temperature on each pressure level until the heating rate reduces to zero, to within some tolerance, for all model levels (Iro et al. 2005;Malik et al. 2017). We refer to this second approach as the 'time-dependent' approach, since it involves iterating through time-steps until a steady-state is found.
Within ATMO, we have recently implemented a time-dependent approach to solving for radiative equilibrium, following the same basic method outlined in previous studies (Iro et al. 2005;Malik et al. 2017). The governing energy equation is where is temperature, is time, is the specific heat capacity, is the gas density, and is altitude. We note that Equation Eq. (1) can be equivalently expressed in terms of the pressure gradient, rather than the altitude gradient. The net radiative flux, , is calculated using the radiative transfer scheme in ATMO (Amundsen et al. 2014;Drummond et al. 2016). Convective flux is not presently included, though as UHJ atmospheres are dominated by radiative energy transfer, and for the most part convectively stable, we expect convective flux to be unimportant. The -profile that satisfies radiative equilibrum corresponds to a steady-state solution of Eq. (1) (i.e. / = 0, to within some tolerance).
The newly implemented time-dependent radiative equilibrium solver has been validated by comparing the steady-state radiative equilibrium -profile against a -profile calculated using the well-established, and well-tested, time-independent flux-balancing scheme within ATMO. We note that the latter has previously been compared against other 1D atmosphere models for both brown dwarfs and exoplanets (Baudino et al. 2017;Malik et al. 2019). Details of this test are presented in Appendix A.
The atmosphere is discretised onto a pressure grid that is uniformly spaced in the log 10 ( )-space. The vertical grid is then defined by the maximum and minimum pressure and the number of pressure levels. Altitude increments (Δ ) between model levels are calculated assuming local hydrostatic balance.
An adaptive time-stepping method is used to handle the fact that the deep atmosphere evolves much more slowly than the upper atmosphere, where the radiative timescale is much shorter, as applied in Iro et al. (2005) and Malik et al. (2017). When solving for the radiative equilibrium temperature profile, the timestep is allowed to vary independently in each model level, allowing much larger timesteps for model levels at higher pressures. This allows for convergence to be reached in a reasonable number of iterations in the deep atmosphere, while maintaining a timestep that is short enough to remain numerically stable in the upper atmosphere. Applying a different timestep across model levels is acceptable as long as the steady-state solution is sought, though it means that the profiles at intermediate steps before convergence is reached are not physically consistent.
At the start of a new calculation, the timestep in each level is initialised based on the radiative timescale, where is gravity and is the Stefan-Boltzmann constant. At the top of the atmosphere (for = 1 Pa) the timestep is of order ∼ 10 s, while at the bottom of the atmosphere (for = 10 7 Pa) the timestep is typically of order ∼ 10 8 s, though these values strongly depend on temperature. Following each iteration the timestep in each model level is increased by 10%, except if several criteria are met that are introduced to maintain model stability. If the heating rate in a model level has switched to an opposite sign compared with the previous timestep (indicating an oscillating heating rate in time) then the timestep is reduced by 33%. Additionally, if the heating rate in a model level has the opposite sign to both of its neighbouring model levels (indicating a heating rate that is oscillating in space) then the timestep is reduced by 33%. We note that when considering a timedependent phenomenon (such as a varying instellation, as described later) we use a fixed and uniform timestep for all model levels instead.
We also introduce a limit on the magnitude of the temperature increments within each level. This prevents excessively large adjustments to the temperature, which can be present especially in the early stages of a calculation when starting from an initial condition far from the equilibrium state (e.g. an isothermal profile). We set this temperature increment cap to ±50 K.
When solving for radiative equilibrium it is important to define a set of criterion to signal when convergence has been reached. For the steady-state we test the normalised flux gradient, which approaches zero as the profile approaches radiative equilibrium, following the approach of Malik et al. (2017). The criterion for convergence is that for all model levels, where is some set tolerance. We use a tolerance of = 10 −3 . In this work, the time-dependent radiative equilibrium solver described in this section is used to generate initial profiles for the pseudo-2D model described in the next section.

Pseudo-2D modelling
With time-dependency introduced in ATMO, the next step is to continuously rotate the atmospheric column around the equator, creating a pseudo-2D model. This introduces a simple parameterisation of horizontal advection. We follow the method used in Iro et al. (2005). A time-dependent modulating function, , is applied to the incoming stellar flux, which adjusts incident radiation based on the zenith angle of the host star, so that radiation is at a maximum when the star is directly overhead, and zero at all points on the night-side. is the same time as found in Eq. (1). is the column rotation period (not to be confused with the planetary rotation period). The value of is selected, on a case by case basis, in order to achieve the desired effective wind speed at the planetary equator. To determine if a periodic steady-state has been achieved, each subsequent full model rotation is qualitatively compared to the last, as frequently done in GCMs (e.g. Drummond et al. 2020). Qualitative comparison is often used in this way as the deep atmosphere, where the radiative and dynamical timescales are extremely long, does not reach a steady-state in a realistic time due to computational limitations ). Once a periodic  steady-state has been achieved, the model time, , is used as an analogue to the equatorial longitudinal position of the model. All parameters used in the model (Table 2) come from The Extrasolar Planet Encyclopaedia 1 . For the stellar spectra we use the Kurucz spectra 2 for HD 209458b and the BT-Settl models 3 (Allard et al. 2012) for WASP-12b. An example of a full steady-statestructure and the effects of pseudo model rotation, using parameters for HD 209458b, can be seen in Fig. 1. An atmospheric hotspot is found in the day-side atmosphere, with a temperature inversion corresponding to a TiO/VO absorption region. The temperature inversion shifts to much lower pressures on the night-side. Of the four positional profiles shown in Fig. 1a, the substellar point (0 • ) -profile is closest to the -profile calculated using the 1D time-dependent model. The temperature structure for HD 209458b obtained by our pseudo-2D model is broadly similar to that predicted by 3D GCMs (Parmentier et al. 2013, see their Figure 3), for the equatorial region.
It should be reiterated that this simple parameterisation of horizontal advection follows flow in the Lagrangian frame by rotating the column, but does not simulate any material or heat transport horizontally out of the column bounds. It only represents bulk advection by modulating the incident stellar irradiation on the atmospheric column. Hence this model is pseudo-2D, not a full 2D treatment of the atmosphere. In HJ atmospheres, circulation is dominated by 1 http://exoplanet.eu/ 2 http://kurucz.harvard.edu/stars.html 3 https://phoenix.ens-lyon.fr/Grids/BT-Settl/AGSS2009/ the equatorial zonal jet, so pseudo-2D models are useful for these systems.

Thermal dissociation and recombination of hydrogen
When molecular hydrogen, H 2 , reaches temperatures between ∼2000 -4000 K, depending on pressure, it thermally dissociates into atomic hydrogen, H. Atomic hydrogen can then recombine to reform H 2 when the temperature drops, as the reaction is reversible.
Note that a collision partner (indicated as "M" above) is not needed to calculate the thermodynamics of the thermal dissociation/recombination of hydrogen, but is needed to calculate its kinetics (e.g., the chemical timescales).
Reaching the required pressure-dependent temperatures for H 2 dissociation is not uncommon in the day-side of UHJ atmospheres. However, the night-side temperature is far colder. As a result, H 2 in the day-side atmosphere thermally dissociates before being pushed onto the night-side by the equatorial jet, where it recombines (Bell & Cowan 2018). Due to the endothermic/exothermic nature of these processes, the reactions remove energy from the day-side atmosphere and redistribute it to the night-side. If this energy source/sink is large enough, it affects the temperature structure of the atmosphere. As measurements of transit photometry for UHJs also rely on their relative day-night temperature contrast, any significant changes to the temperature structure from this effect can impact observation.
Hydrogen is by far the most abundant element in UHJ atmospheres, so this dissociation/recombination effect has a large impact on the mean molecular weight within the system, and therefore on its heat capacity. In our pseudo-2D model, we use this to introduce the reaction heat redistribution by adding an additional term in the specific heat capacity (Gordon & McBride 1994), where, for a species : is the number of moles per kilogram, 0 is the molar enthalpy, and 0 , is the molar standard heat capacity. The first term on the right of Eq. (5) is the 'frozen' heat capacity, which is simply an abundance weighted sum of the individual heat capacities of the mixture. The second term on the right is the 'reaction' heat capacity, which accounts for the temperature-dependent chemistry. This reaction term becomes significant when a species with high molar enthalpy changes abundance over a small temperature range, which is the case for thermal dissociation/recombination of hydrogen. For a full derivation of this equation see Appendix B. The partial derivative ln ln is evaluated using a Newton-Raphson iterative method analogous to that used to calculate the equilibrium abundances, as described in Gordon & McBride (1994, see their Section 2.5.1).
The specific heat capacity as a function of pressure and temperature (for a solar elemental composition) is shown in Fig. 2, when only the frozen heat capacity is included ( Fig. 2(a)) and when both the frozen and reaction heat capacities are included ( Fig. 2(b)). The specific heat capacity significantly increases for particular pressure-temperature combinations, corresponding to the dissociation/recombination of H 2 /H. This means the reaction only has a significant impact on atmospheric structure when these pressure-dependent temperature combinations are achieved, and the atmosphere horizontally transitions between an H 2 and H-dominated regime.
This method of introducing heat redistribution from the dissociation/recombination reaction is equivalent to the method used in Bell & Cowan (2018), whereby the heat redistribution is input directly through the governing energy equation. The tracer method used in Tan & Komacek (2019) is also equivalent. Our approach is simply the 'chemical' method, as opposed to statistical physics. However, by altering the heat capacity calculation, our model is more generalised and can be used to recreate the heat of any other chemical transition reactions we may want to study in the future. We choose to focus on H 2 dissociation/recombination here as it is the most impactful for the -range of UHJs. This method has also been used in the latest version of the Met Office Unified Model (UM), creating a field to implement variable heat capacity (Drummond et al. 2018a).

THE REACTION HEAT REDISTRIBUTION IN UHJS
WASP-12b has an unusually high day-side temperature, in places reaching ≥ 3000 K (Swain et al. 2013), and an equilibrium temperature of ∼2500 K (Chakrabarty & Sengupta 2019). This atmosphere represents an interesting case between planetary and stellar conditions. We select WASP-12b as indicative of UHJs, but recognise that in some respects this planet is not typical for this class of objects (Bell et al. 2019). However, as we are interested in the mechanism of high temperature thermal dissociation and recombination, and not specifically characterising this planet in detail, omission of the additional processes relevant to WASP-12b will not hinder our wider conclusions.
In our pseudo-2D model, TiO and VO are initially excluded. The    impact of including these species is explored in Section 3.2. Here we use a timestep of 200 s and a column rotation period of 1.728 × 10 5 s, giving a wind speed of ∼5 km s −1 (4830 m s −1 ), which is scaled proportionally for alternate wind speeds investigated in Section 3.3. This is a reasonable estimate given that the zonal-mean zonal wind speed is found to saturate at this value for temperatures above 1500 K (Komacek et al. 2017).

Excluding TiO/VO
To analyse the effects of the reaction heat redistribution within this atmosphere, the pseudo-2D model is first run without the reaction heat capacity included (Eq. (5)). The initial profile is calculated from the 1D time-dependent model (Section 2.2). The model is then run again, this time including the reaction heat capacity. Fig. 3 shows the fractional change in steady-state heat capacity structure once the reaction heat capacity is applied. The fractional change in steady-state chemical abundances of H/H 2 and corresponding -structure can be seen in Fig. 4 and Fig. 5 respectively. First, as shown in Fig. 3, we see that when the reaction heat capacity is applied the atmospheric heat capacity increases drastically on the day-side, most prominently between pressures of 10 −1 -10 −2 bar. This region corresponds to an area in the atmosphere at which the largest variations in H 2 and H abundance occur, due to the dissociation/recombination affect, and as a result the -profile aligns with the reaction heat capacity's optimal -range (seen in Fig. 2).
This can clearly be seen when inspecting the steady-state chemical distributions found in Fig. 4. As thermal dissociation of H 2 occurs, between the substellar point and evening terminator, it cools the dayside atmosphere. Consequently the temperature is lowered and falls below the optimal reaction -range, thus reducing the amount of dissociation occurring. This means, when including the reaction heat capacity, there is an increased abundance of H 2 at the peak dissociation point, with a decrease at the peak recombination point and throughout the night-side atmosphere. The same logic, in reverse, applies to the H abundance. Recombination heating leads to increased abundance just beyond the evening terminator and throughout the night-side, but a decrease in the day-side atmosphere between the substellar point and evening terminator.
To see the effects of heat redistribution caused by the reaction heat capacity, we can look at a comparison between the steady-state -structures (Fig. 5). Firstly, it is interesting to note that a slight temperature inversion is present in the atmosphere at 10 −1 bar, despite the lack of TiO and VO as opacity sources. This temperature inversion occurs as the radiative timescale above the ∼ 1 bar pressure level is short enough for the day-side atmosphere to respond to stellar heating. At higher pressures, the radiative timescale is much longer than the rotation timescale so the atmosphere does not respond rapidly enough to the radiative heating, and a temperature inversion forms. This inversion is decreased when the reaction heat capacity is included in our model. As the radiative timescale is roughly proportional to the atmospheric heat capacity, increased heat capacity due to chemical transition reactions makes the atmosphere more stable against the radiative heating. In other words, a higher incident energy would be required to maintain the temperature inversion for the same wind speed.
By calculating the relative change in longitudinal temperature, seen in Fig. 6, we observe a clear reduction in the global temperature gradient, as the day-side cools and night-side warms. This is, unsurprisingly, most prominent between 10 −1 -10 −2 bars, where we observe the largest horizontal chemical variation. In this case, the atmosphere between the substellar point and evening terminator is seen to cool by over ∼10% in places, or ∼400 K. The night-side atmosphere is even more drastically affected, heating by up to ∼20% just beyond the evening terminator. This also corresponds to a temperature change of ∼400K, as the percentage change is relative to the original temperature structure. Perhaps the most interesting consequence of these changes in temperature is a significant overall shift in the atmospheric hotspot towards the evening terminator.
In our model, the temperature changes are mostly confined to the upper regions of the atmosphere, at pressures below 1 bar, where strong irradiation can penetrate. It should again be mentioned that our model lacks time-dependent non-equilibrium reactions, as we assume local chemical equilibrium, meaning thermal dissociation/recombination occurs instantly. If a more realistic timescale for the thermal dissociation/recombination were to be considered, the   atmospheric hotspot would likely shift even further from the substellar point, depending on the competition between chemical, radiative and advective timescales.

Effects of TiO and VO
TiO and VO are thought to be the primary drivers for steep temperature inversions observed on the day-side of HJ atmospheres (Désert et al. 2008). Both TiO and VO display strong photo-absorption in the near UV/optical range, causing strong heating in the upper atmosphere, and acting as a soft irradiation barrier to lower altitudes. However, there is an ongoing debate about the role of TiO/VO in atmospheres with particularly large day-night temperature contrasts. Although previous studies have shown that these compounds are likely to be thermally dissociated on the hotter day-side, specifically for WASP-121b which has a similar (albeit slightly cooler)structure than WASP-12b (Parmentier et al. 2018), they may condense out of the night-side atmosphere entirely (Parmentier et al. 2018;Merritt et al. 2020). Furthermore, both TiO and VO are currently very difficult to detect directly (Mikal-Evans et al. 2020), as this would require high resolution observations in the near-UV/optical range. In this section, we investigate the effects of the reaction heat redistribution on atmospheric structure, with TiO and VO included in our model.
The fractional change in steady-state heat capacity structure, now including TiO and VO, can be seen in Fig. 7. From this we see that the heat capacity now drastically increases in a thin band within the day-side atmosphere at pressures lower than 10 −2 bar. This corresponds to the edge of the TiO/VO photo-absorption region. Because of the temperature inversion, the pressure-dependent temperature for the reaction heat capacity is now met predominantly at the peak of TiO/VO absorption. This can clearly be seen in the steady-state distributions of both H and H 2 abundance (Fig. 8), where the vast majority of thermal H 2 dissociation is now confined to the day-side atmosphere above the 10 −2 bar pressure level. However, the dissociation/recombination has an increasingly limited impact on the heat capacity inside the region, although the optimal pressure-dependent temperature is met throughout, as the TiO/VO opacities almost en-tirely dominate the temperature structure. The result is a large increase to the heat capacity in a band around the edges of the TiO/VO absorption region; with changes on the morning edge caused by H 2 dissociation and on the evening edge caused by H recombination. Some H 2 dissociation is also found to occur in the upper atmosphere, above the 10 −4 bar pressure level, as strong absorption by TiO/VO maintains a small temperature inversion throughout the night-side atmosphere and the temperature required for the H 2 to H transition reduces towards lower pressures. This causes some cooling at these altitudes on the night-side.
By looking at the the periodic steady-state temperature structures for the simulations (Fig. 9), we can see the H 2 dissociation therefore reduces the magnitude of the temperature inversion on the morning edge of the hotspot and shifts its peak to lower-pressure levels. Conversely, on the evening side of the hotspot, the thermal H recombination heating contributes to significantly increase the temperature at the edge of the TiO/VO absorption region. As previously seen in the case when TiO/VO have been excluded, this again causes a strong shift in the atmospheric hotspot towards the evening terminator.
The dominance of TiO/VO within the primary dissociation region also results in a significant disparity between the heating and cooling in this case. From Fig. 10 we can see there is a large magnitude increase in the heating peak, from ∼20% (∼400 K) up to over ∼40% (∼800 K), when comparing with the TiO/VO exclusion case. Whereas the cooling remains at a similar value of ∼5 -10% (between ∼200 -400 K). We do, however, now see equal cooling at high altitudes in the night-side atmosphere, as predicted by the increased H abundance seen in Fig. 8.

Effects of wind speed
In order to consider the sensitivity of the reaction heat redistribution to our selected wind speed, the model is re-run a number of times with different column rotation periods, to simulate different dynamical conditions excluding, and then including, the reaction heat capacity. For each column rotation period, the resulting percentage temperature change has then been determined (Fig. 11). This is done both excluding and including TiO/VO. As wind speed is increased from ∼1 km s −1 up to ∼10 km s −1 (R=864000 s and R=86400 s are used, giving speeds of 966 m s −1 and 9660 m s −1 respectively), some key trends become apparent. Other than for the 10 km s −1 TiO/VO exclusion case, which is discussed shortly, the magnitude of both relative cooling and heating from the reaction heat redistribution is seen to increase with increasing wind speed in both cases excluding and including TiO/VO. Heating from H recombination is also consistently affected more than the cooling from H 2 dissociation in the presence of TiO/VO, for the same reason as discussed in Section 3.2 (dominance of TiO/VO opacity in the primary photo absorption region). Protracted heating, further from the primary recombination zone, and low pressure cooling (≤ 10 −4 bar) throughout the night-side atmosphere is additionally increased with higher wind speeds.
These effects can be explained by considering the balance between advective and radiative timescales within the atmosphere. As the wind speed is set by the rotation rate of the column, this rate emulates the advective timescale which decreases as the model wind speed is increased. The atmospheric temperature takes a finite amount of time to respond to heating, and the radiative timescale is independent of rotation, so the radiative heating has diminishing effects at higher wind speeds as the advective timescale becomes comparatively small. The chemical timescale, however, is time-independent regardless of the model wind speed, due to the assumption of chemical equilibrium. This means that, at higher wind speeds, the cooling and heating from the reaction heat redistribution increasingly outweigh the radiative heating and cooling within the atmosphere, and therefore have a progressively larger effect on the -structure.
The exception to this trend is seen when comparing the increase from 5 to 10 km s −1 in the TiO/VO exclusion case. Here the reaction heat redistribution is actually observed to decrease at the higher wind speed. This can be explained by considering the effects that a faster jet has on the atmosphere. The temperature inversion decreases as the dynamical timescale begins to outweigh the radiative timescale and the day-night temperature contrast decreases as the faster jet more thoroughly homogenises the atmosphere. As a consequence, the decreased temperature progressively falls below the optimal temperatures for the H 2 /H thermal transition, effectively capping the associated reaction heat capacity at a value set by the threshold wind speed and minimising the overall impact of the reaction heat redistribution. This threshold effect is not observed when TiO/VO are present, as in this case the temperature inversion is set by the photoabsorption of these species rather than purely by radiative heating, so these chemicals stabilise the day-side atmosphere at the required transition temperature.
Unfortunately, the pseudo-2D radiative solver currently lacks any dynamical or chemical feedback on the wind speed, so it is unclear if the reaction heat redistribution would directly impact the jet, particularly at the equator where it is strongest. Conversely, the wind speed could be the primary moderator within the system, and control the dissociation/recombination fractions. Studies of this effect using a 3D GCM have shown the speed of the equatorial jet to decrease if the rotation period is fixed (Tan & Komacek 2019). However, we show here that chemical feedback is also very important. To fully explore the dynamics, a 3D global circulation model, with time-dependent chemical feedback, would need to be used. This is a promising avenue to explore in the future, but is not possible until the treatment of chemical feedback within higher dimensionality models is much improved.

Consequences for the pseudo-2D phase curve
By computing the emission spectra for individual steady-stateprofiles in the atmosphere at 5 • intervals, using 500-band correlatedcross-sections, we produce pseudo-2D phase curves for our model. A sample of these emission spectra can be seen in Appendix D. These pseudo-2D phase curves, at 3.6 m and 4.5 m for cases both excluding and including the TiO/VO species and reaction heat capacity, are shown in Fig. 12. Given our model assumptions (i.e equilibrium chemistry, zero convection, no dynamical feedback, no clouds or haze) and its 2D nature, meaning only longitudinal contributions to the atmospheric emission are considered, these cannot be considered as true phase curves and should not be taken as completely realistic. They can, however, provide a very reasonable approximation to how H 2 dissociation/recombination may impact observation, especially in the presence of TiO/VO. Phase curve data in Bell et al. (2019) shows roughly a 4000 ppm ratio in day-side to stellar flux at 3.6 m and 4500 ppm at 4.5 m. This aligns well with our model, particularly in the cases including TiO/VO where we calculate ∼3500 ppm and ∼4500 ppm ratios respectively. In both scenarios we explore, the spectral flux amplitude damping we would expect to observe due to the reduction in global temperature gradient is present. This damping is much more significant in the case excluding TiO/VO, ∼400 ppm (∼13.5%) compared to ∼100 ppm (∼3%) at 3.6 m and ∼300 ppm (∼9%) compared to ∼75 ppm (∼2%) at 4.5 m, which is also to be expected as we have already seen the reaction heat redistribution to have a more protracted impact on the atmosphere in this instance. This damping has been studied in a number of previous papers, and has been used to explain some flattening observed in UHJ phase curves (Tan & Komacek 2019). The other notable feature in both pseudo-2D phase curves is a clear phase offset increase when the reaction heat redistribution is considered. This offset increase is very similar in both cases at 3.6 m, with a value of ∼15 • . At 4.5 m it is slightly more significant in the TiO/VO exclusion case, where the offset is increased by ∼15 • , compared to ∼ 10 • when TiO/VO are included. Of course, this will also vary significantly on a case by case basis, and be heavily influenced under more rigorous dynamical treatment. As we have seen, increasing wind speeds push the hotspot, and the resulting phase offset, further from the substellar point as the reaction heat redistribution increasingly outweighs the radiative temperature changes.

CONCLUSIONS
During this study a time-dependent pseudo-2D model for ultra-hot Jupiter exoplanets has been constructed and used, based on the 1D/2D atmosphere model ATMO. This new model is capable of accounting for the heating/cooling effect of any chemical transition reactions, and has been applied to the heat redistribution caused by the thermal dissociation and recombination of hydrogen, under time-independent equilibrium chemical conditions. From the results of this study, the following key conclusions can be drawn: • Building upon the work of Bell & Cowan (2018), displaying the importance of H 2 dissociation/recombination in atmospheres exceeding 2000K, we have shown how the presence of certain species (TiO,VO) changes the impact of any resulting heat redistribution.
• The day-side atmospheric temperature is found to be cooled by up to ∼10% (∼400 K) due to thermal H 2 dissociation in our model. Thermal H recombination heating has an even higher relative impact, reaching over ∼20% (∼400 K) just beyond the evening terminator. With TiO/VO included within the model, cooling at the TiO/VO absorption peak is seen to remain at a similar value of ∼10% (∼400 K), due to the dominance of TiO/VO opacity within the photo-absorption region. However, heating at the evening terminator is found to be much higher at over 40% (∼800 K).
• The reaction heat redistribution causes a shift in the atmospheric hotspot towards the evening terminator, which is found to be more significant when excluding TiO/VO from our model. This consequently results in a positive shift in phase offset for the pseudo-2D phase curves. The increased heat redistribution has also been demon-strated to cause significant spectral flux amplitude damping in the pseudo-2D phase curves. This amplitude damping is far more significant when TiO and VO are not present.
• With TiO and VO excluded, increasing the model wind speed, up to a certain threshold value, has been found to increase the magnitude of atmospheric heating and cooling and the distance of protracted heating and cooling throughout the night-side atmosphere. At a threshold wind speed, the reaction heat redistribution is capped by reducing atmospheric temperatures and above this the magnitude of both heating and cooling decreases as the atmosphere becomes increasingly homogenised. TiO/VO photo-absorption stabilises hotspot temperatures, meaning that in the presence of these species the H 2 /H transition temperatures are achieved regardless of wind speed. Protracted night-side cooling at pressures below 10 −4 bar is also observed in this case. Any increase in the heating and cooling acts to increase the resulting hotspot shift at higher wind speeds (or lower the shift above the threshold wind speed when TiO/VO are not present).
Although it is noted that a more rigorous investigation, one including dynamical feedback, into this effect is required to draw more robust conclusions.
There are of course limitations to our model which could be addressed. By far the most significant omission comes from inherent limitations due to dimensionality, and therefore limited dynamical feedback, with the wind speed being a prescribed model parameter. Because of this, the relative heating and cooling calculated here is likely to be much larger than expected and observed. Extending our model into 3D would drastically reduce both, as the reaction has its strongest effect at the equator. The presence of night-side clouds, which are thought to significantly impact atmospheric dynamics and radiative transfer in hot Jupiters (Helling 2020), are also not included here. These clouds have been shown capable of altering observational phase offset by counterbalancing day-side atmospheric emission, meaning that the offset does not necessarily track the atmospheric hotspot (Parmentier et al. 2020).
Any variations in planetary equilibrium temperature, and resulting day-night temperature gradient, would also alter the reaction heat redistribution. The parameters used in our model, roughly simulating WASP-12b, highlight a particularly extreme case, with a much higher equilibrium temperature than average, so we would expect most UHJ exoplanets to have a decreased global temperature gradient and therefore display diminishing effects from the reaction heat redistribution. Some assumptions, such as chemical equilibrium, would have only a limited effect on the transition reaction. The inclusion of time-dependent non-equilibrium chemical timescales would act to increase the protracted nature of heating and cooling, as in our model any dissociation/recombination happens instantly, although this would only be significant for high wind speeds due to the fast chemical timescales. For a wind speed of ∼5 km s −1 , the H 2 -H interconversion chemical timescales is calculated, on the basis of the present equilibrium abundances, to be shorter than the advection timescale below the 10 −5 bar pressure level at 0 • longitude (for H 2 →H+H) and below the 10 −2 bar pressure level at 180 • longitude (for H+H→H 2 ). Variation in planetary mass, and corresponding surface gravity would also have a significant impact, see Appendix C for a brief discussion on this. Interactions between the highly ionised atmosphere and strong magnetic field have also been shown to impact the phase offset in drastic and unpredictable ways (Rogers & Komacek 2014). These magnetic effects are far beyond the scope of this paper, but do provide a promising further avenue for future improvements in UHJ modelling.
Fortunately, very recent progress within Exeter's exoplanet group has provided breakthroughs when running consistent chemistry with 3D models (Drummond et al. 2018b(Drummond et al. , 2020. This study could therefore soon be repeated, utilising the Met Office's 3D atmosphere model , to more thoroughly investigate feedback between the reaction heat recirculation and atmospheric dynamics, whilst maintaining chemical consistency.

APPENDIX A: TESTING THE TIME-DEPENDENT RADIATIVE EQUILIBRIUM SOLVER
The new time-dependent 1D model is benchmarked against a calculation using the well-establish time-independent method implemented in ATMO for HD209458b (see Fig. A1), to check for consistency. Parameters used can be seen under HD209458b in Table 2. A heat redistribution factor of 0.5 and daytime average zenith angle of 60 • were also used for the time-dependent 1D model. These tests are done both with and without TiO and VO included as opacity sources. These two molecules are strong UV/optical stellar flux absorbers, causing a temperature inversion in the atmosphere of HD209458b Figure A1. Comparison between the steady-state -profiles from the time-independent 1D model (dashed lines) and the time-dependent 1D model (solid lines) for the model parameters used in Section 2.3, with TiO/VO either excluded (red lines) or included (blue lines) as atmospheric opacity sources and many other hot Jupiters (Désert et al. 2008). As expected, the additional molecules put more strain on the radiative solver, causing some minor deviations to appear at higher altitudes. Overall, however, there is a remarkable degree of agreement between the two approaches, given they use separate methods to determine atmospheric -structure. The time-independent radiative-equilibrium solver has previously been compared against other 1D atmosphere models with generally very good agreement (see Drummond et al. 2016;Baudino et al. 2017;Malik et al. 2019).

APPENDIX B: THE HEAT CAPACITY OF A MIXTURE
The heat capacity of a substance at constant pressure is defined as Here is the molar heat capacity with units J mol −1 K −1 . The specific heat capacity, , or the heat capacity per unit mass, has units J kg −1 K −1 . The specific and molar heat capacities are intensive properties of the substance. By contrast, the extensive form of heat capacity , measured in [J K −1 ], is related to the intensive forms via where is the mass of the substance, and is the total number of moles. Following Gordon & McBride (1994), we now derive the specific heat capacity of a mixture at chemical equilibrium. We start with the total enthalpy of a mixture, which is given by where ℎ is the total specific enthalpy of the system in J kg −1 , is the number of moles per kilogram of species in mol kg −1 and 0 is the molar enthalpy of the species in J mol −1 . Now, using the definition of the heat capacity, we find the total specific heat capacity of the mixture: where 0 is the molar standard heat capacity of the species . It is clear that the total specific heat capacity of the mixture is a combination of two parts. The first term on the right in the final equation represents the sum of the molar heat capacities of the individual species (multiplied by the number of moles per kilogram, in this case). This component is referred to as the 'frozen' heat capacity by Gordon & McBride (1994). The second term, called the 'reaction' heat capacity, contains a derivative of the number of moles (per unit mass) with temperature. This term becomes important when the abundance of a species varies rapidly with the temperature, for instance at a pressure-temperature region where the chemistry transitions from a CO-dominated to a CH 4 -dominated atmosphere, or in the region of a condensation curve. The molar heat capacity can be easily obtained by multiplying the specific heat capacity by the mean molecular weight of the substance where is in kg mol −1 .

APPENDIX C: EFFECT OF PLANETARY MASS/SURFACE GRAVITY
In order to briefly investigate the effect of planetary mass, and corresponding surface gravity, on the reaction heat redistribution we ran the model selecting a much lower planetary mass of 0.32 (log( )=2.36). A comparison between the heat redistribution in this scenario and a much higher mass of 1.47 (log( )=3.02) can be seen in Fig. C1. In the lower mass case, the impact of the reaction heat redistribution is significantly increased both with or without TiO and VO. This is because lower surface gravity increases the radiative timescale (Eq. (2)), causing the radiative heating and cooling of the atmosphere to have a decreasing relative effect when compared to the reaction heat redistribution. As the deep atmosphere takes much longer to respond to the stellar irradiation, the primary locations of both dissociation and recombination are also seen to shift towards lower pressure levels. Protracted low pressure night side cooling is reduced in the lower mass case for the same reason. As the deep atmosphere experiences less radiative heating, the high altitude cooling required for energy balance is decreased, see Section 3.3 for more detail.

APPENDIX D: EMISSION SPECTRA
In order to construct the pseudo-2D phase curves seen in Section 3.4, individual -profiles are extracted from the steady-statestructures at 5 • longitudinal intervals. These -profiles are then used in ATMO's atmospheric emission routine, to calculate a series of emission spectra from which the pseudo-2D phase curves are constructed for wavelengths of 3.6 and 4.5 . A sample of these emission spectra, excluding (Fig. D1) and including (Fig. D2) TiO/VO, ranging from 0.3-5 at 30 • intervals in longitude can be seen here. This paper has been typeset from a T E X/L A T E X file prepared by the author. Figure C1. Percentage change in temperature when including the reaction heat capacity with TiO/VO either excluded (top) or included (bottom) as atmospheric opacity sources and with a planetary mass of 0.32 (log( ) = 2.36) (left) or 1.47 (log( ) = 3.02) (right), respectively. Wind speed is set at 5 km s −1 . Figure D1. Thermal emission spectra calculated from -profiles situated every 30 • longitude around the equator, for the models excluding TiO/VO as atmospheric opacity sources. Figure D2. Thermal emission spectra calculated from -profiles situated every 30 • longitude around the equator, for the models including TiO/VO as atmospheric opacity sources.