GCM-Motivated Multidimensional Temperature Parametrization Scheme for Phasecurve Retrieval

We present a novel physically motivated, parametrized temperature model for phase-curve retrieval, able to self-consistently assess the variation in thermal structure in multidimensions. To develop this approach, we drew motivation from both full three-dimensional general circulation models and analytic formulations, accounting for the dominant dynamical feature of tidally locked planets, the planetary jet. Our formulation shows notable flexibility. It can generate planetary jets of various characteristics and redistribution efficiencies seen in the literature, including both standard eastward and unusual westward offset hotspots, as well as more exotic configurations for potential future observations. In our modeling scheme we utilize a tractable set of parameters efficient enough to enable future Bayesian analysis and, in addition to the resolved temperature structure, we return physical insights not yet derived from retrievals: the amplitude and the phase offset, and the location and the extent of the equatorial jet.


Introduction
The multidimensional nature of planetary atmospheres presents many challenges to modeling efforts. For decades our understanding of complex exoplanetary atmospheres has depended largely on two space telescopes, Spitzer and Hubble, and several ground-based telescopes (Knutson et al. 2012(Knutson et al. , 2014Kreidberg et al. 2014bKreidberg et al. , 2014aLockwood et al. 2014;Brogi et al. 2018). Although the number of confirmed exoplanets and numerical and analytical approaches to characterize their atmospheres have grown significantly over the years, the limited wavelength coverage and spectral resolution of the telescopes have allowed us only partial insight into their envelopes (Deeg & Belmonte 2018). With some confidence we have been able to speak about global hemispheric averages and one-dimensional (1D) quantities, but the complicated three-dimensional (3D) structure has remained elusive.
With the recent James Webb Space Telescope (JWST), Transiting Exoplanet Survey Satellite (TESS), and CHaracterizing ExOPlanet Satellite (CHEOPS) launches and the expectation of two new missions in the near future, PLAnetary Transits and Oscillations of stars (PLATO) and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL), our prospects of unveiling the multidimensional atmospheric structure are dramatically improving (Stevenson et al. 2016;Edwards et al. 2019). The number of exoplanets available for a thorough atmospheric characterization will quickly jump by an order of magnitude, yielding opportunities to select excellent targets over a wide range of properties (e.g., radius, mass, distance to the host star, stellar type).
Attempts to assess the multidimensional atmospheric structure in retrieval started in 2016, when Feng et al. (2016) explored the biases resulting from 1D assumptions within a Bayesian atmospheric retrieval. They showed that for WASP-43b synthetic HST and Spitzer phase-curve data the methane abundance was constrained to incorrect values under the 1D assumption. In Blecic et al. (2017) we explicitly compared the secondary-eclipse temperature structure produced with a 3D simulation and the retrieved 1D model. We show that the retrieved model did not exist anywhere on the planet, demonstrating in this way that tidally locked planets cannot be characterized in 1D. This idea has been furthered in forward modeling frameworks by several groups. Flowers et al. (2019) and Beltz et al. (2021) calculated very-high-resolution emission spectra from a suite of GCM models for a range of planetary rotation rates, while Caldas et al. (2019) developed a tool to generate transmission spectra from 3D temperature profiles and applied it to GJ 1214b and WASP-121b (Pluriel et al. 2020).
Based on these works, a growing number of groups are incorporating multidimensionality into retrievals. Broadly speaking, these approaches define several representative temperature-pressure profiles and combine their respective contributions using novel geometrical approaches that account for the changing viewing angle throughout the planet's orbit. In this manner, Changeat et al. (2019) defined three distinct profiles (day, night, terminator), while MacDonald et al. (2020) used varying toy models for east and west terminators to demonstrate that retrievals that assume a single temperaturepressure profile give erroneously cold temperatures at the terminators. Feng et al. (2020) extended their 2016 work to allow for viewing from arbitrary orbital phases and applied it to the phase curve of WASP-43b. Taylor et al. (2020) used an analytic formulation and toy models to explore the signal-tonoise ratio required to correctly identify an inhomogeneous temperature structure using two temperature-pressure profiles. Irwin et al. (2020) began the exploration of latitudinal temperature dependence by imposing a cosine term when calculating disk-averaged spectra. Finally, a slightly different approach is taken by the recently initiated work of K. Chubb et al. (2022, in preparation) who derive a 3D temperature model based on the diffusion equation. While this attempts to use a physical phenomenon when going from 1D to 3D, there may be limitations associated with this approach as horizontal advection cannot be fully modeled as a diffusive process.
Though these approaches recognize the multidimensional nature of the atmosphere, none of them start from the intrinsic atmospheric processes and address the fundamental dynamical feature of close-in tidally locked planets, the circumplanetary equatorial jet. Short-period planets (from hot Jupiters to terrestrial planets), presumably tidally locked and highly irradiated, tend to exhibit strong dynamical forcing (e.g., Showman & Guillot 2002;Dobbs-Dixon & Lin 2008). The uneven and intense stellar irradiation of their atmospheres leads to large longitudinal temperature gradients closely tied to supersonic equatorial and mid-latitude jets that advect energy around the planet. Inherently 3D, complicated atmospheric dynamics of these planets leads to significant compositional and temperature variations, most prominently seen with changes in longitude. To adequately constrain a plausible atmospheric structure, one must rigorously compare our understanding of these physical and chemical processes with the observational constraints. Thus, combining our theoretical knowledge gathered from GCMs with a robust statistical, observationally driven retrieval approach, utilizing photometric and spectroscopic phase-curve observations, is our best way forward to assess the complicated 3D atmopsheric structure of exoplanetary objects.
Phase-resolved emission spectrometry (phase-curve observations) carries the most comprehensive information about planetary envelopes, atmospheric dynamics, chemical processes, and radiative energy balance (e.g., Knutson et al. 2007;Cowan & Agol 2008;Stevenson et al. 2014;Venot et al. 2020). These observations allow us to map the entire surface of a planet, as for a tidally locked planet the orbital phase is equivalent to its rotational phase, providing, in addition to the altitudinal information, the longitudinal and latitudinal constraints on the atmospheric composition, thermal structure, and energy distribution. Furthermore, in the thermal infrared, the amplitude of the phase variation determines the day-night temperature contrast and the phase offset reveals the longitude of the planet's hottest point. These dynamical features derive primarily from the equatorial jet redistributing heat from the hot dayside to the cooler nightside.
In this paper we use analytic formulations and GCM models to account for the planetary jet, and fundamentally link all orbital phases in a more physically consistent way. Our parameterization scheme accounts for the energy redistribution around the planet by including both advective and radiative timescales. Here, we present a forward model and we explore the flexibility of our formulation. In a following paper, we will present the implementation of this model in retrieval, utilizing phase-curve observations. Such application will allow us to self-consistently tie together temperature profiles across all longitudes and latitudes, enabling a simultaneous retrieval of spectra at all orbital phases, and returning important constraints on the atmospheric dynamics.
The paper is structured as follows. In Section 2 we briefly outline our methodology. In Section 3 we define our 2D temperature model, list the model parameters, and explore the versatility of our assumptions. In Section 4 we describe additional assumptions for the full 3D representation, and in Section 5 we summarize our results.

Methodology
We use radiative-hydrodynamic (RHD) simulations, employing a double-gray radiative scheme to model both the radiative and the advective contributions to the total energy budget. The dynamical portion of our RHD model is described in Dobbs-Dixon & Agol (2013).
The radiative contribution was chosen to match the analytic formulation of Guillot (2010), while to isolate and parameterize the advective contribution we explicitly use the RHD outputs (see Sections 3.1 and 3.2). Our prescription models all the crucial features seen in the GCM simulations responsible for the advection of energy from the day-to the nightside of the planet. By using an amenable set of parameters, we allow for a large flexibility in our models, covering both solutions seen in the GCM simulations and physical scenarios not yet explained by the theory but observed among exoplanets (e.g., Dang et al. 2018;von Essen et al. 2020). Our temperature parameterization is physically motivated and fully analytic and it allows us to calculate the planetary 3D temperature structure in microseconds, as opposed to days or weeks using standard GCM simulations. A reasonable set of model parameters also permits for an efficient Bayesian analysis, allowing for thousands of models to be explored in a short period of time. Such formalism returns the physically consistent resolved temperature structure, but also, for the first time, the fundamental properties of the jet.

2D Temperature Profile
To develop our 2D (p, long) temperature parameterization scheme we start with a static, plane-parallel, isotropically scattered atmosphere in local thermodynamic equilibrium, assuming the condition for a purely radiative equilibrium: Here, κ ν is the absorption coefficient, J ν is the source function, and B ν is the Planck function. To derive the temperature structure, we follow Guillot (2010) and introduce two Eddington coefficients, f kth and f hth , for the thermal radiation pressure and flux, respectively (Eddington 1916). We decouple the thermal and visible radiation and parameterize the opacities using the mean visible and thermal opacities, defining their ratio as γ = κ υ /κ th . The gas temperature as a function of the optical depth, τ, can then be written as The first term in this equation describes the spherically symmetric interior energy contribution, where T int is the internal temperature of the planet, τ is the infrared optical depth, given as τ = p κ th /g, p is the pressure, g is the planetary gravity, and m = * long lat cos cos ( ) ( ), with "long" being the longitude and "lat" being the latitude. The second and third terms define the angle-dependent irradiative contribution, where T irr is the irradiation temperature, given as with T * , R * , and a being the stellar surface temperature, radius, and semimajor axis, respectively. Equation (2) is easily solved analytically, and we denote this radiative temperature contribution as T rad 4 . However, this purely radiative approach does not account for the important advection of thermal energy by the jet. The planetary jet links temperature profiles across all planetary phases by removing energy from the dayside and adding energy onto the nightside.
At a fundamental level, one can introduce an advective flux by modifying the radiative equilibrium condition to be where q ▿ T represents the horizontal advective flux. Total radiative absorption and emission must now be additionally balanced with thermal energy advected in and out of an area. The vertical temperature profile can then be written as  where T rad 4 carries both internal and irradiation terms (first, second, and third term in Equation (2)) and T adv 4 carries the last two advective terms in Equation (5) (the τ dependence is dropped for clarity). Several other authors (e.g., Guillot 2010;Hansen 2008;Burrows et al. 2008) have formulated similar expressions to Equation (5), but a purely analytic solution to the total temperature structure has so far proven elusive. Thus, to formulate an expression for T adv 4 we have turned to results from our GCM simulations. We note that T rad 4 and T adv 4 are not truly decoupled, as several parameters (γ, f hth , and f kth ) appear in both terms of Equation (6). However, since the model was formulated for phase-curve retrieval, it will reveal the global planetary properties (the consistent 3D temperature structure and jet characteristics) from the coupled, total planetary temperature (T p ) model (not the individual components), making irrelevant (or secondary) the fact that the individual components are in fact coupled.

Analytic Radiative Component
To parameterize individual contributions to T p from Equation (6), we simulate a planet resembling HD 189733b, utilizing our RHD model with a dual-gray radiative scheme. We utilize two sets of GCM simulations: one fully selfconsistent, cloud-free model and the other identical except for the introduction of a very large artificial damping of the velocities. By damping the velocities to near zero, we are effectively shutting off any advective contribution. The resulting solution is purely radiative (and we denote it as T rad ) and agrees quite well with the analytic solution from Guillot (2010; see Figures 1 to 3 and Section 3.2) and the expression given in Equation (2), which is fully described by five parameters: T int , T irr , γ, f hth , and f kth (refer to Guillot 2010 andBlecic et al. 2017for an exhaustive study of the effect of these parameters on the temperature-pressure profiles). We note that in our parameterization the Eddington parameters ( f hth and f kth ) can be effectively excluded from the model. As Guillot (2010) states in their Section 2.4, Eddington coefficients can carry only a few possible values; however, only the solution assuming an isotropic radiation field ( = f hth 1 3 , = f kth 1 2 ) closely matches the theoretical solution at higher pressures. In our subsequent analysis, we adopt these values, but we also leave the user an option to use different values. This leaves us with only three free parameters for the radiative solution (T int , T irr , γ).

Parameterized Advective Component
By subtracting the GCM's radiative solution, T rad 4 , from the full GCM solution, we are able to isolate the advective contribution, T adv 4 , and subsequently model it. In this section, we present our parameterization of the advective component in 2D, as a function of pressure and longitude. In Section 4, we introduce additional parameters for a full 3D representation (as a function of pressure, longitude, and latitude).
To derive our parameterization, we ran three different GCM simulations, a hot, medium, and cold, assuming irradiation temperatures of 2500, 2000, and 1500 K, respectively, and carefully studied all the features seen in the results. Figures 1 to  3 show a comparison between the radiative, advective, and planetary temperatures generated using these GCM simulations and our parametrized model. We also show the temperature differences for each component, calculated by subtracting the GCM results from the most closely matching parametrized model.
For the radiative component (T rad ; Figures 1, 2 and 3, left double panel), in both the damped GCM simulations (left sides) and our parametrized solutions (right sides), as expected, the profiles are getting colder going from the substellar point toward the nightside. On the nightside, the temperature of the radiative profiles remains constant, as there is no contribution from the stellar flux and no mechanism to transfer the energy to the nightside of the planet. We see nice matches between the GCM simulations and the parametrized models in all three cases.
For the advective component (T adv 4 ; Figures 1, 2 and 3, middle double panels), we see that for each case the temperatures are negative on the dayside and positive on the nightside, illustrating the extraction of energy from the dayside and deposition on the nightside. The left side panels, generated by subtracting the strongly damped GCM simulation from the full GCM simulation, clearly show that the efficiency of 2D advection is a strong function of both longitude and pressure. Therefore, to parameterize T adv 4 and allow for sufficient flexibility in the parameters, such that one can model a wide range of jet characteristics, we identify a number of different characteristics evident in the figures. First, in longitude, we  distinguish three segments across the planetary surface: east dayside (seg 1 , 0°-90°), nightside (seg 2 , 90°-270°), and west dayside (seg 3 , 270°-360°). In each of these segments T adv 4 has a different functional dependence in longitude (see next paragraph). Second, to account for the pressure dependence, also clearly seen in the figures, the magnitude of this term within the segments is allowed to vary. In particular, we find that the values of T adv 4 at the interfaces between the segments (three in total, approximately corresponding to 0°, 90°, and 270°longitude) are well modeled as Gaussians in log pressure. At each interface, we introduce parameters for the amplitude (A i , in units of K 4 ), width (σ i , in units of bar), and center (c i , in units of bar) of these Gaussians, which can freely vary. We further note that the cores of equatorial jets found in GCM simulations of tidally locked planets are usually isobaric across the planet, allowing us to define only a single central pressure for all three interfaces (e.g., Rauscher & Menou 2010). This leaves us with a total of seven parameters (A 1 , A 2 , A 3 , σ 1 , σ 2 , σ 3 , c). Utilizing these, the value of T adv 4 at each of the three interfaces, denoted by Γ(p), can be calculated as Here, p is the pressure and the index-i denotes the interfaces: (i) between east and west daysides; (ii) between east dayside and the nightside; and (iii) between the nightside and west dayside. Figure 4 shows an example of the Gaussian T adv 4 profiles for a set of fiducial values.
To parameterize the behavior of T adv 4 across the three regions in longitude (seg 1 , seg 2 , seg 3 ) we again turn to the GCM results and find that the east and west dayside regions are well modeled with the sum of parabolic and error functions, while the nightside is best fit with a linear function ( Figure 5). We further introduce longitudinal offset parameters (f 2 and f 3 , in units of degrees), that can shift the two interfaces (dotted vertical lines in Figure 5, right panel) in any direction, determining in this way where most of the energy transferred from the dayside will be deposited. This induces eastward or westward shifts of the temperature profiles, suggestive of the observed phase offsets (see also Section 3.2.2). Thus, we model each of these segments analytically as The parabolic, linear, and error functions are given as In the above equations, Γ i (p) is defined by Equation (7), p is pressure, and j is longitude. Equations (9) and (10) denote the eastern and western parabolic expressions, Equation (11) gives the linear expression for the nightside change in longitude, while Equations (12) and 13 denote eastern and western error functions. Figure 5 shows a comparison between the GCM advective component for the medium case (Figure 2, panel 2) at a pressure of ∼0.1 bar, and the corresponding parametrized  Figure 6 shows each of the segments in a different color for all pressures (given as hue gradient) and longitudes, for the same set of model parameters as in Figure 2, panel 2, right side. For illustration purposes, we have set the location of the interfaces to be 0°, 90°, and 270°on both Figures 5 and 6, though as shown in the above equations these can vary.
Taken together, we can now calculate a fully parameterized T adv 4 as a function of longitude and pressure. These parameterized quantities are shown in Figures 1, 2, and 3, panels 2, right sides and, as seen, they agree quite well with the GCM results shown in the same figures, panels 2, left sides.
Adding the advective term (T adv 4 ) to the radiative term (T rad 4 ), which is calculated analytically using Equation (2), we get the total 2D temperature structure for each case. Figures 1, 2, and 3, panels 3, show the comparison between the full GCM simulations and the parametrized total planetary temperature model including both radiative and advective components. As we can see, the GCM simulations agree nicely with our bestmatched parametrized planetary temperature solution.
The bottom panels in Figures 1, 2, and 3, showing the temperature differences between the GCM simulations and parametrized models, demonstrate that the differences are between 5% and 10%, on the order of the uncertainties usually seen when retrieving planetary temperatures. Differences in the radiative component are most pronounced at the terminators and at depth. The differences at the terminators are due to the abrupt change in longitude in the analytic formulation that is smoothed out in the highly damped GCM simulations, while the differences in the depth are likely due to the fact that the GCM has not converged at high pressures. For the advective contribution and total planetary temperature, the differences seen in the plots are indicative of wave-like structures known to be present in GCM models, but not included in the analytic formulation.

Model Flexibility
Our analytic temperature parameterization allows us to generate physically consistent T(p, long) profiles throughout the entire planet, subject to a set of three parameters governing the radiative solution, T rad (κ υ , κ th , T int ; see Equation (2)), and nine governing the advective solution, T adv (A 1 , A 2 , A 3 , σ 1 , σ 2 , σ 3 , c, f 2 , f 3 ). The flexibility of our parameterization, also seen in Figures 1-3, allows one to explore a wide range of jet parameters and characteristics, crucial for use in retrieval.  (7)). The Γ 1 parameter describes the transfer of energy between ∼0 and ∼90°, Γ 2 between ∼90°and ∼270°, and Γ 3 between ∼270°and ∼360°longitude.   Figure 7 shows several temperature structures generated using our prescription, demonstrating the effectiveness of our assumptions (model parameters are listed in Table 1). The left column shows temperature profiles at the equator across all longitudes, corresponding to the different jets (T adv 4 ) shown in the middle column, together with the global planetary temperatures across all pressure shown in the right column. As presented, we can decide on the location of the jet, whether the energy from the dayside will be taken from wider or narrower pressure range, and on which longitudes and pressures most of this energy will be deposited. The differences in the radial extent, strength, and direction of the jets are also clearly visible when comparing the three middle figures. This allows us to explore the standard eastward temperature offsets, but also unusual configurations such as a westward offset (e.g., Corot-2b, WASP-33b, as in Dang et al. 2018 andvon Essen et al. 2020, respectively).

Flux, Phase Curves, Spectra
Once we have specified the temperature structure, we can easily generate flux for each wavelength and at each point on the planetary surface. To calculate intensities, we perform radiativetransfer calculations using the open-source PYRATBAY framework for exoplanet atmospheric modeling, spectral synthesis, and Bayesian retrieval (Cubillos & Blecic 2021). In Figures 8  and 9, we use our 2D temperature parameterization scheme and assume hydrogen-dominated atmospheres, accounting for major spectroscopically active species and their opacities (CO, CO 2 , CH 4 , H 2 O, HCN, NH 3 ). Our code can, however, assume any chemical composition and account for any opacity sources from ExoMol 4 and HITRAN/HITEMP 5 , as well as other databases (see Sections 2.3.1-2.3.5 in Cubillos & Blecic 2021).
To calculate flux, we first divide the atmosphere into equally spaced longitude-latitude boxes (visible pixels on Figures 8  and 10), where each box has its own temperature profile, calculated either using our 2D temperature parameterization scheme (Section 3), or the 3D scheme (Section 4). Then, we calculate the observer's projected area (Seager 2010) at each phase ( Figure 10 shows examples of the observer's projected area at phases 0°s and 270°). Finally, using PYRATBAYʼs  Table 1. Left panels: temperature profiles at the equator across all longitudes, with the contribution function envelope (given in dashed gray) calculated at 4.5 μm. Middle panels: T adv 4 at the equator. Right panels: corresponding total planetary temperatures as a function of longitude and pressure. radiative-transfer routines, we calculate the intensity per box, convolve it with the corresponding observer's box area to get the "flux per box", generating in this way the total spatially dependent emission flux across the observable hemisphere. Repeating this routine for each phase allows us to create a phase curve for the desired wavelength. The phase curve then returns three fundamental properties of the jet: the phase offset, amplitude, and pressure extent. Figure 8 shows the spatially dependent emission flux in examples of the eastward and westward hotspots, and the corresponding phase curves with the phase shifts at 170°and 200°, respectively (for observers this corresponds to the peaks that occur at +10°, top panel, and −20°, bottom panel, before and after opposition, respectively). These offsets are produced for a planet observed at 4.5 μm at secondary eclipse. The pressure location of the jet can be recovered from the center parameter, c, while the amplitude and the phase offset can be easily extracted from the phase curve. The flexibility of our approach allows us to generate the offset hotspot at any longitude on the dayside (or even nightside) of the planet.
To generate spectra at each phase across a range of wavelengths, for each wavelength we calculate "flux per box" and add them over the observable atmosphere to get the surface flux. Figure 9 shows the phase-dependent spectra from a planet with and without an advective contribution. The increase in the flux from the nightside due to the jet's influence is clearly seen on the right panel.

3D Planetary Structure
Our ultimate goal is to develop a physically motivated 3D T (p, long, lat) approach for phase-curve retrieval to model the Top row 4.0 × 10 −3 10 −2 550 0 0 −10 12 1.3 × 10 13 10 13 10 −5 10 −  true multidimensional nature of planetary atmospheres. This will allow us to link both the planetary longitudes and latitudes in a physically consistent way. To extend our 2D T(p, long) model to 3D T(p, long, lat), we again utilize the results of GCM simulations. The equatorial jet seen in tidally locked atmospheres has a latitudinal extent closely connected to the Rhines scale (Showman & Guillot 2002), thus dependent on the wind velocity and rotation rate. Fast winds and slow rotation rates (compared to Jupiter) imply wide jets. Simulations indicate that the jet usually extends between ∼±25°latitude around the equator (e.g., Dobbs-Dixon & Agol 2013). To take advantage of the north/south symmetry and mimic this behavior, we assume another Gaussian with a variable width in latitudinal direction (σ lat ) symmetric around the equator, and integrate it into our parameterized T adv 4 model. The radiative part of our T(p, long, lat) model will also be attenuated in latitude by including lat cos( ) in the μ parameter, following Guillot (2010).
Since we are already dividing the planet into longitudinallatitudinal cells, this allows a smooth transition to 3D when calculating flux, spectra, or phase curves. In Figure 11 we show examples of 3D advective temperature structures, for two different latitudinal extents of the jet, while in Figure 12 we show the planetary temperatures across all longitudes and latitudes at several pressure layers. We note that utilizing a Gaussian in latitude is a simplification, and future work should utilize the results of GCMs to improve the latitude temperature formulation in a way similar to our longitudinal parameterization.

Summary
In this paper, we presented a novel physically motivated, parameterized multidimensional model for phase-curve retrieval that accounts for the advection of energy due to the planetary jet. This parameterization, informed by GCM models, accounts for the energy redistribution around the planet by including both radiative and advective processes. Instead of treating each phase independently, we fundamentally link them together and self-consistently redistribute energy throughout the entire planet via the equatorial jet. In addition to the resolved temperature structure, the model also returns essential information about the planetary jet structure, the location and extent of the jet, and its amplitude and phase offset. These recovered jet properties provide us with physical insight into Figure 10. Observers projected area on the planet at orbital phases of 0°(left) when looking directly at the nightside, and 270°(right). The total hemispherical area sums to one. Figure 11. Examples of advective temperature structures at a pressure of 0.02 bar, generated utilizing our 3D T(p, long, lat) model and the jet characteristics shown in the top panels of Figure 7. Left: T adv 4 with a Gaussian width of 30°in the latitudinal direction. Right: same as left panel with a Gaussian width of 55°. As seen, the advection shape, direction, and strength seen in Figure 7 are well preserved in both left and right panels.  Figure 11 and the jet characteristics shown in the top panel of Figure 7.