{'fraction_non_alphanumeric': 0.02868175044031974, 'fraction_numerical': 0.01972632434629454, 'mean_word_length': 4.834400442652755, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Aura-3D: A Three-dimensional Atmospheric Retrieval Framework for Exoplanet Transmission Spectra

Atmospheric retrievals of exoplanet transmission spectra allow constraints on the composition and structure of the day-night terminator region. Such retrievals in the past have typically assumed one-dimensional temperature structures which were adequate to explain extant observations. However, the increasing data quality expected from exoplanet spectroscopy with JWST motivates considerations of multidimensional atmospheric retrievals. We present AURA-3D, a three-dimensional atmospheric retrieval framework for exoplanet transmission spectra. AURA-3D includes a forward model that enables rapid computation of transmission spectra in 3D geometry for a given atmospheric structure and can, therefore, be used for atmospheric retrievals as well as for computing spectra from General Circulation Models (GCMs). In order to efficiently explore the space of possible 3D temperature structures in retrievals, we develop a parametric 3D pressure-temperature profile which can accurately represent azimuthally-averaged temperature structures of a range of hot Jupiter GCMs. We apply our retrieval framework to simulated JWST observations of hot Jupiter transmission spectra, obtaining accurate estimates of the day-night temperature variation across the terminator as well as the abundances of chemical species. We demonstrate an example of a model hot Jupiter transmission spectrum for which a traditional 1D retrieval of JWST-quality data returns biased abundance estimates, whereas a retrieval including a day-night temperature gradient can accurately retrieve the true abundances. Our forward model also has the capability to include inhomogeneous chemistry as well as variable clouds/hazes. This new retrieval framework opens the field to detailed multidimensional atmospheric characterisation using transmission spectra of exoplanets in the JWST era.


INTRODUCTION
The study of exoplanet atmospheres has seen immense progress in recent years, with spectroscopic observations enabling constraints on a range of physical and chemical processes that occur in exoplanets . One of the most successful observational techniques used to characterise exoplanet atmospheres has been transmission spectroscopy, which measures the wavelengthdependent decrease in flux from the host star as the planet transits in front of it. This approach has led to detections and abundance constraints of chemical species mnixon@ast.cam.ac.uk nmadhu@ast.cam.ac.uk in numerous exoplanets (e.g., Charbonneau et al. 2002;Snellen et al. 2010;Deming et al. 2013;Mandell et al. 2013;Kreidberg et al. 2014a;Madhusudhan et al. 2014;Wyttenbach et al. 2015;Wakeford et al. 2018;Pinhas et al. 2019) as well as providing constraints on the properties of clouds and hazes (e.g., Pont et al. 2008;Bean et al. 2010;Pont et al. 2013;Kreidberg et al. 2014b;Sing et al. 2016; Barstow et al. 2017;Pinhas et al. 2019;Benneke et al. 2019).
One of the most common methods for deriving atmospheric properties from spectroscopic data of exoplanets is atmospheric retrieval (Madhusudhan & Seager 2009). An atmospheric retrieval algorithm involves combining an atmospheric forward model with a statistical sampling algorithm in order to determine estimates of the atmospheric properties that best explain Figure 1. Schematic of the Aura-3D retrieval framework. The algorithm combines an atmospheric forward model in 3D geometry with a Bayesian sampling algorithm to conduct parameter estimation. The forward model is modular in structure, with numerous features that can be incorporated as needed, such as clouds/hazes, refraction, scattering and stellar heterogeneity.
In order to create an accurate and flexible atmospheric model for use in retrievals, a large number of parameters are required to capture all possible compositions, thermal structures, and other properties such as the presence of clouds or hazes. Since statistical sampling algorithms often require a large number of model evaluations, increasing as the number of free parameters increases (Fortney et al. 2021), it is necessary to make some simplifying assumptions to reduce computation time. Common assumptions include the use of a 1D temperature profile and uniform chemical abundances in the region of the atmosphere probed by the observations (Madhusudhan 2018). It is important to re-examine the assumptions made in these algorithms as we enter a new era in the study of exoplanet atmospheres, in which the quality and volume of observational data will improve immensely thanks to next-generation facilities such as the James Webb Space Telescope (JWST; Greene et al. 2016), the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL; Tinetti et al. 2018) and the Extremely Large Telescopes (ELTs; Snellen et al. 2015).
Despite the successes of the retrieval approach for analysing data from existing observatories such as the Hubble Space Telescope (HST), several studies have highlighted the potential for 1D retrieval codes to lead to biased results due to variations in temperature structure and chemical abundances with latitude and longitude, particularly when considering high-quality observations of very hot Jupiters from upcoming facilities such as JWST. For example, assumptions of a 1D atmosphere have been shown to bias retrievals of emission spectra (e.g., Blecic et al. 2017;Taylor et al. 2020) and can lead to unrealistic model fits to phase curves (e.g., Irwin et al. 2020;Feng et al. 2020).
This work focuses on retrievals of transmission spectra. Transmission spectroscopy probes the day-night terminator region of the atmosphere, across which inhomogeneities may be expected (e.g., Fortney et al. 2010). These effects should be particularly prevalent in hot and ultra-hot Jupiters, since day-night temperature contrasts are expected to increase as planetary equilibrium temperature increases (Cowan & Agol 2011;Komacek & Showman 2016). Previous studies have already indicated that 1D retrievals can lead to biases when interpreting transmission spectra of hot giant planets. Caldas et al. (2019) developed a transmission spectrum model that could incorporate a 3D atmospheric structure, which they used to investigate the effect of daynight temperature gradients on model spectra. They used the 1D retrieval code TauREx (Waldmann et al. 2015) to analyse synthetic JWST spectra generated with their 3D model for which the temperature could vary significantly between the day-and nightsides. The retrievals returned terminator temperatures that were biased towards the dayside temperature, and chemical abundances which were not consistent with input values. Pluriel et al. (2022) also found day-night induced biases for a range of synthetic JWST-like spectra of hot and ultra-hot Jupiters, and Pluriel et al. (2020) used the same 3D model to explore the potential for biases caused by day-night chemical heterogeneities. They focused on the ultra-hot Jupiter WASP-121b, for which the H 2 O abundance is expected to be much higher on the nightside than the dayside due to thermal dissociation (Parmentier et al. 2018). A 1D retrieval analysis showed that inferred chemical abundances could be significantly different to the input values. MacDonald et al. (2020) also discussed the potential for biased temperatures and abundances from 1D retrievals (but cf. Welbanks & Madhusudhan 2022). A generalised retrieval framework should be able to address the above issues in order to accurately infer atmospheric properties regardless of the presence of thermal inhomogeneities. This is the goal of the present study.
Unlike the simplified atmospheric models used in retrieval algorithms, General Circulation Models (GCMs) can be used to capture some of the more complex processes at work in planetary atmospheres. GCMs are detailed 3D models that simulate atmospheric dynamics, some also incorporating chemical processes along with radiative transport (Showman et al. 2020). Originally developed to simulate the Earth's atmosphere, a wide range of GCMs have been adapted for application to exoplanet atmospheres (Showman & Guillot 2002;Cooper & Showman 2005;Lewis et al. 2010;Rauscher & Menou 2010;Thrastarson & Cho 2010;Wordsworth et al. 2011;Polichtchouk & Cho 2012;Dobbs-Dixon & Agol 2013;Kataria et al. 2013;Amundsen et al. 2016;Mendonça et al. 2016;Way et al. 2017;Deitrick et al. 2020). These models take various different approaches to computing atmospheric dynamics. In general, GCMs aim to solve the primitive equations of atmospheric dynamics in 3D across the entire planet, which are simplifications of the Navier-Stokes equations that assume hydrostatic equilibrium, a shallow atmosphere compared to the planetary radius, and constant gravity with height. However, some models instead opt to solve the full Navier-Stokes equations (e.g., Dobbs-Dixon & Agol 2013;Mayne et al. 2014;Amundsen et al. 2016).
GCMs have predicted a number of important features of hot Jupiter atmospheres. This includes the finding that hot Jupiters can show large day-night temperature contrasts of up to serveral hundred K (e.g., Showman & Guillot 2002;Cooper & Showman 2005;Amundsen et al. 2016). GCMs of hot Jupiters have also predicted equatorial superrotation that can lead to an eastward shift of a planet's dayside hotspot away from the substellar point (e.g., Showman & Guillot 2002;Cooper & Showman 2005;Rauscher & Menou 2010;Kataria et al. 2013). These predictions have subsequently been confirmed through comparison with observations of infrared phase curves (e.g., Knutson et al. 2007;Komacek et al. 2017;Stevenson et al. 2017).
Model transmission spectra can be generated using 3D temperature structures. Fortney et al. (2010) used precomputed temperature profiles from SPARC/MITgcm simulations of hot Jupiters (Showman et al. 2009) to examine the differences between 1D and 3D model atmospheres. A number of more recent works have also developed multidimensional forward models for transmission spectra. Caldas et al. (2019) used a 3D radiative transfer model to investigate biases caused by day-night temperature gradients using a simple parametric temperature profile. Falco et al. (2022) presented transmission spectrum models in 1D, 2D and 3D, and MacDonald & Lewis (2022) described parametric prescriptions for 3D atmospheric models.
Although GCMs are too computationally expensive to incorporate into retrieval algorithms directly, a number of studies have aimed to bridge the gap between 1D and 3D models in the context of retrievals of transmission spectra. Lacy & Burrows (2020) conducted chemical equilibrium retrievals assuming separate dayside and nightside temperature profiles, showing that in some cases it is possible to constrain day-and nightside temperatures of hot Jupiters from their transmission spectra. Espinoza & Jones (2021) demonstrated how JWST could be used to acquire separate transmission spectra for each of a planet's limbs, and also expanded the retrieval framework CHIMERA ) to enable retrievals of a single transmission spectrum with separate temperature profiles for each limb, under the assumption of chemical equilibrium. Welbanks & Madhusudhan (2022) also presented a 1+1D retrieval framework in which the final transmission spectrum is calculated by a linear combination of two separate spectra, one representing the morning terminator and another representing the evening terminator.
In this work we present Aura-3D, a 3D atmospheric retrieval framework for transmission spectra of exoplanets. In Section 2, we describe our algorithm in detail. We present a forward model that incorporates a 3D temperature structure and which can be used to efficiently generate transmission spectra from the output of a GCM, similarly to several past studies Caldas et al. 2019). We develop a parametric pressure-temperature (P -T ) profile that can match the azimuthally-averaged structure of a GCM while being very fast to compute, making it suitable for atmospheric retrieval. In Section 3.1 we compare our parametric temperature profile to a number of established GCMs. We explore the effects of a day-night temperature gradient on resulting transmission spectra in Section 3.2 by quantifying the difference in transit depths between 3D models and their 1D-averaged counterparts. In Section 3.3 we discuss the possible effects of chemical inhomogeneity on transmission spectra. In Section 3.4 we demonstrate the capability of Aura-3D to carry out retrievals with a multidimensional P -T profile by conducting a retrieval on a synthetic hot Jupiter spectrum generated using our new parametric profile. In Section 3.5 we retrieve a spectrum generated using a GCM temperature structure, comparing results from retrievals both with and without a day-night temperature gradient. We summarise our findings and discuss avenues for the future in Section 4.

METHODS
The retrieval framework developed in this study is shown in Figure 1. Similarly to previous retrieval algorithms, this framework combines an atmospheric forward model with a Bayesian parameter estimation scheme. The key differences between our methods and those of previous retrieval algorithms concern the forward model. Our forward model enables the computation of a transmission spectrum whose temperature structure can vary in three dimensions throughout the atmosphere, and can therefore take the temperature structure calculated by a GCM as its input. However, since the calculation of a temperature structure using a GCM is very computationally expensive, we also describe a number of parameterisations that allow for rapid computation of a multidimensional temperature profile, and which can be incorporated directly into a retrieval. Our retrieval framework also incorporates several components from the existing retrieval code Aura and its subsequent implementations (Welbanks & Madhusudhan 2019;Nixon & Madhusudhan 2020) as well as Aurora (Welbanks & Madhusudhan 2021). These features are discussed later in this section. We note that Aurora can also be applied to H-poor atmospheres, which are not considered in the present work.

Transit Geometry
We define a spherical polar coordinate system (r, θ, φ) with an origin at the centre of the planet (see Figure  2). The coordinates (R p , 0, 0) correspond to the substellar point. The zenith angle θ varies between 0 and π, reaching π/2 at the day-night terminator and π and the antistellar point. The azimuthal angle φ varies between 0 and 2π, increasing in the direction of the trailing hemisphere of the planet. This geometry, based on the star-planet-observer axis, differs from a traditional latitude-longitude grid based on the rotation axis which is often used when describing the geometry of a planet. We adopt this approach since it allows for a straightforward description of day-night variations in atmospheric properties, since the zenith angle θ moves from the substellar to antistellar point. This geometry has been used in previous studies involving transmission spectroscopy (e.g., Fortney et al. 2010;Caldas et al. 2019). Figure 3 depicts a ray of light travelling from the host star to an observer, passing through a portion of the atmosphere of a transiting planet. A single ray is assumed to have a fixed impact parameter b and a fixed azimuthal angle φ as it moves through the atmosphere. The ray travels along the path s, defined so that s = 0 when the ray is directly above the planet's terminator (θ = 90 • ). As it travels along this path, the zenith angle θ increases (moving from the dayside to the nightside of the planet) and the distance between the ray and the centre of the planet (denoted r) also varies. The values of θ and r are given by and cos θ = s r . (2)

Radiative Transfer
The following expression describes the transit depth ∆ λ for the generalised atmosphere of a planet with radius R p transiting a star with radius R * : where Figure 2. Three-dimensional co-ordinate system adopted for our forward model. The zenith angle θ is defined to be 0 at the substellar point and π (180 • ) at the antistellar point. The azimuthal angle φ varies between 0 and 2π (360 • ), increasing in the direction of the trailing limb.
In this equation, the atmospheric height is denoted by H, and τ (b, φ) represents the optical depth at impact parameter b and azimuthal angle φ. The total transit depth ∆ λ is found by integrating the transit depth at each value of φ over all azimuthal angles. A crucial difference between the 3D transmission spectrum model and a 1D approximation is that the attenuation coefficient at a given wavelength, µ λ , depends on r, θ and φ, as opposed to only depending on r. The expression for the optical depth along a ray path s at a given (b, φ) is therefore Previous retrieval codes have included forward models with φ-dependence to capture properties such as inhomogeneous cloud cover (Line & Parmentier 2016;Mac-Donald & Madhusudhan 2017;Welbanks & Madhusudhan 2021). Aura-3D also incorporates φ-dependence to model inhomogeneous clouds (see Section 2.5.1). By including a temperature profile that depends on θ (see 2.3.3), our retrieval algorithm is therefore capable of constraining atmospheric properties in three dimensions.
Our model incorporates absorption from a wide range of chemical species following the methods described in Gandhi & Madhusudhan (2017, using crosssections derived from line list data from a number of different sources. A Voigt function is used to apply temperature and pressure broadening. We assume a H/He- The attenuation coefficient of a given chemical species i can be expressed as where ρ i is the mass density of species i, κ i is the opacity of species i, n i is the number density of species i and σ i is the absorption cross-section of species i. The number density n i of a given species can be related to its volume mixing ratio, X i = n i /n tot , where n tot is the total number density. The volume mixing ratio of each chemical species, apart from H 2 and He, is a free parameter in the model. The volume mixing ratios of H 2 and He are given by where the value of 0.17 is derived from a solar composition (Asplund et al. 2009).
In past implementations of Aura, cross-sections of each species are stored on a three-dimensional grid of λ, P and T . In order to determine the total attenuation coefficient at a given height r in the atmosphere, the cross-sections of each chemical species are interpolated using the values of P and T found at r using the equation of hydrostatic equilibrium and the prescribed P -T profile. This yields a two-dimensional array µ i (λ, r) for each species. In the present work, since each combination of values of φ and θ may have a different P -T profile, the cross-sections must be interpolated for each (φ, θ) pair, yielding a four-dimensional array µ i (λ, r, θ, φ). The attenuation coefficient for all chemical species acting as absorbers in the model is obtained by summing the attenuation coefficients for each individual species.
We also include collision-induced absorption (CIA) due to H 2 -H 2 and H 2 -He (Richard et al. 2012). The attenuation coefficient due to CIA is given by the expression µ CIA = X H2 n 2 tot [X H2 σ H2−H 2 (λ, T ) + X He σ H2−He (λ, T )] (9) where σ H2−H 2 and σ H2−He are the H 2 -H 2 and H 2 -He cross-sections respectively. We note that unlike the cross-sections for chemical species described above, which have units of m 2 , the CIA cross-sections have units of m 5 .

A Multidimensional Parametric Temperature Profile
Our radiative transfer model is capable of taking any 3D temperature structure as its input. The model can therefore be used for post-processing GCMs to create transmission spectra. However, given that a single GCM typically takes at least several days to compute, such a model is not appropriate for atmospheric retrieval, which can require 10 6 model computations in order to analyse a single spectrum (Fortney et al. 2021). The temperature profile must therefore be simplified to enable rapid computation. The simplest approach is to consider an entirely isothermal temperature profile when carrying out retrievals (e.g. Waldmann et al. 2015;Zhang et al. 2019). The analytic temperature profile described in Guillot (2010) has been incorporated into retrieval algorithms (e.g., Benneke & Seager 2012;Line et al. 2013), and the parametric temperature profile described in Madhusudhan & Seager (2009) has also been employed in numerous retrieval frameworks (e.g., Pinhas et al. 2018;Blecic et al. 2022). While the latter two approaches result in a temperature profile that varies with height in the atmosphere, these profiles do not vary with longitude or latitude.
Aura-3D incorporates several different parametric prescriptions for P -T profiles that include temperature variations in multiple dimensions. This enables fast, flexible model computation while also allowing for varying levels of complexity in the temperature structure. We present three possible prescriptions: one in which temperature varies with height in the atmosphere but not with θ, one in which temperature varies with θ but not with height, and one in which temperature varies with height as well as with θ and φ.

Temperature varying with height only
In the case where the temperature profile does not vary with θ or φ, we adopt the 1D parameteric P -T profile from Madhusudhan & Seager (2009). The atmosphere is divided into three layers, defined by P 1,2,3 . The deepest layer is isothermal while the upper layers have thermal gradients controlled by α 1,2 . The full profile is defined as: where (P 0 , T 0 ) defines the conditions at the top of the atmosphere. This can be recast to yield expressions for T (P ): This P -T profile was developed to be capable of emulating observed temperature profiles of solar system planets as well as profiles from self-consistent exoplanet atmosphere models. This flexibility, along with a rapid computation time, makes it an ideal prescription for retrievals. Our new profile presented in section 2.3.3 applies the same philosophy to 3D models in order to extend our retrieval framework to incorporate multidimensional effects.

Temperature varying with θ only
This model provides a simple means of incorporating day-night temperature variations across the terminator region. It is described by four main parameters. The extent of the terminator region over which the temperature varies is fixed by an angle β, similarly to Caldas et al. (2019). The other three parameters are T term , the temperature in the middle of the terminator region (θ = π/2), T day , the temperature at the dayside end of the terminator (θ = π/2 − β/2), and T night , the temperature at the nightside end of the terminator (θ = π/2 + β/2). The temperature at intermediate values of θ is found by linearly interpolating between the two appropriate temperatures. Outside of the fixed terminator region, the temperature remains uniform at either T day or T night . This approach differs from the formalism described in Lacy & Burrows (2020) as it uses three different temperatures rather than just two (T day and T night ). This allows for the temperature gradient between the terminator and the dayside to be different to the gradient between the terminator and the nightside, which has been shown to be possible from GCMs, such as the model of HD 209458b presented in Fortney et al. (2010).

Temperature varying in three dimensions
Our new temperature parameterisation generalises the 1D parametric P -T profile described in Section 2.3.1 to allow for temperature variations in 3D. In order to achieve this, we compute three separate P -T profiles located at three different values of θ. We calculate one profile at θ = π/2, the exact centre of the terminator, as well as profiles at the edges of the terminator which are closest to the day-and nightsides of the planet. We label the dayside profile T d (P ), the nightside profile T n (P ), and the terminator profile T t (P ). The terminator edges are located at θ = π/2±β/2, where β is defined as in the previous section. At a given (P, θ) where θ is intermediate between the two terminator edges, the temperature is obtained by linear interpolation in θ between the two appropriate profiles: where G(φ) is a generic function that can be used to describe the φ-dependence of the temperature profile.
For the purposes of this study we are mostly interested in exploring day-night temperature contrasts, and so we fix G(φ) = 1, however this can easily be modified to include φ-dependence using a formalism similar to the patchy cloud approach described in Section 2.5.1.
Assuming that the P -T profile is continuous throughout, it can be described using 6 free parameters. This suggests that the above parameterisation should require a total of 18 parameters. However, temperature variations with longitude/latitude are typically negligible in the deeper atmosphere (e.g., Showman et al. 2008), and so we can assign common values of P 3 and T 3 for each of the profiles, reducing the total number of free parameters to 14. A demonstration of the efficacy of this parameterisation is presented in Section 3.1.

Inhomogeneous chemistry
Chemical inhomogeneities are likely to persist in a wide range of planets, and will therefore be important to consider when analysing upcoming data. Of particular importance will be day-night variations in chemical abundances. These variations can occur in ultra-hot Jupiters, planets with equilibrium temperatures in excess of 2000 K on which molecules such as H 2 and H 2 O can be thermally dissociated on the dayside (Lothringer et al. 2018;Parmentier et al. 2018;Venot et al. 2020). They may also be present on warm to hot Neptunes, since their atmospheric composition varies strongly as a function of temperature and bulk atmospheric properties such as metallicity and C/O ratio (Moses et al. 2013). High-metallicity atmospheres for these planets are known to exhibit strong day-night temperature contrasts (Lewis et al. 2010), which could lead to a scenario in which the atmospheric composition could vary strongly between the day-and nightsides of the planet.
Our model can incorporate chemical compositions which vary in three dimensions. This extension requires the expansion of the number density of each species n i , as defined in equation 6 from a single number to a three-dimensional array n i (r, θ, φ). For the purpose of this study, we consider a simplified case in which two abundances are specified for each chemical species: a 'dayside' abundance (θ < π/2) and a 'nightside' abundance (θ > π/2). The effect of day-night chemical inhomogeneities on transmission spectra are explored in Section 3.3. Considerations for retrievals with inhomogeneous chemical abundances are reserved for a future study (see section 4.3 for further discussion).

Additional physical effects
Our model incorporates a wide range of other atmospheric properties. These are discussed in more detail in previous works (e.g., Pinhas et al. 2018;Welbanks & Madhusudhan 2021) and are briefly described here.

Clouds and hazes
Our model includes a prescription for inhomogeneous cloud cover adapted from Welbanks & Madhusudhan (2021), which divides the atmosphere into (up to) four azimuthal slices. These slices can have one of four properties: (1) cloud/haze free, (2) covered by a grey cloud deck, (3) covered by hazes, or (4) covered by a grey cloud deck with hazes above the deck. The grey cloud deck is essentially opaque (Fortney 2005) and is therefore implemented by setting the optical depth to infinity for where a is the Rayleigh enhancement factor, γ is the scattering slope, and σ 0 = 5.31 × 10 −31 m 2 is the crosssection due to H 2 Rayleigh scattering at a reference wavelength λ 0 = 3.5 × 10 −7 m (Dalgarno & Williams 1962). The attenuation coefficient due to hazes, µ λ,haze is therefore given by the following expression: where X H2 is the H 2 abundance and n tot (P, T ) is the total number density. The above prescription can easily be simplified to yield the cloud model from MacDonald & Madhusudhan (2017) by considering only slices with properties (1) and (4) (1) and (2).

Stellar heterogeneity
In order to account for the effect of star spots and faculae on transmission spectra, we follow the approach of Pinhas et al. (2018) and incorporate the treatment of stellar heterogeneity from Rackham et al. (2017) into our atmospheric forward model. In cases where stellar heterogeneity is considered, the expression for the observed transit depth, ∆ λ,obs , is where ∆ λ is defined in Equation 3 and E het is a perturbative term describing stellar contamination. E het is defined as where S u and S o are the average spectral energy distributions of the unocculted and occulted regions of the stellar surface, and f het is the areal fraction of the projected stellar disk that is covered with cool spots and/or hot faculae. The spectral components of the star are computed by interpolating from the PHOENIX grid of stellar atmospheric models (Husser et al. 2013) for stellar effective temperatures exceeding 2300 K. For stars cooler than 2300 K we instead interpolate between the DRIFT-PHOENIX model grid (Witte et al. 2011), which enables consideration of stars with temperatures down to 1000 K.

Forward scattering and refraction
Our forward model incorporates the analytic prescriptions for forward scattering and refraction presented by Robinson et al. (2017). The optical depth is typically calculated by integrating along dτ λ = µ λ ds (see equation 5). When forward scattering is included, the optical depth is modified to become where f scat is the correction factor for forward scattering andω is the forward scattering albedo. Refraction is implemented by calculating P max , the pressure at which the effect of refraction is sufficient to cause a ray of light coming from one end of the planet to bend off the far limit of the star (Robinson et al. 2017). The optical depth at pressures greater than P max is set to infinity. The effects of forward scattering and refraction on transmission spectra were investigated by Robinson et al. (2017) for Jupiter-like atmospheres and by Welbanks & Madhusudhan (2021) for mini-Neptunes. In both cases the effects were shown to be small compared to effects such as Rayleigh scattering and collisioninduced absorption.

Synthetic Data and Statistical Inference
For this study we generate synthetic JWST observations and transmission spectra using PANDEXO (Batalha et al. 2017). For all simulated observations, we assume a noise floor of 5 ppm, and a saturation limit of 80% full well. For atmospheric retrieval, we follow the binning strategy of Pinhas et al. (2018) in order to compare the high-resolution forward model with the synthetic data.
Parameter estimation is carried out using PyMultinest (Buchner et al. 2014), a Python implementation of the Nested Sampling algorithm (Skilling 2004;Feroz & Hobson 2008;Feroz et al. 2009). For each parameter θ we aim to find the posterior probability distribution p(θ|d) given the data d being analysed: where L = p(d|θ) is the likelihood, p(θ) is the prior and Z = p(d) is the Bayesian evidence, which is not required for parameter estimation but can be used for model comparison. We assume independently distributed Gaussian errors for each of the spectral data points, meaning the likelihood is defined as whereȳ i and σ i are the mean and standard deviation of the observed data point i, andŷ i is the value of the model data point i. The Nested Sampling algorithm samples the prior space in order to calculate the evidence Z, and in the process of doing so computes the likelihood L which can be used to estimate the posterior distribution.

Model validation
In order to verify that our 3D forward model is accurate, we reproduce effects reported in previous works that have incorporated a 3D transmission spectrum model. We consider the set of model transmission spectra of the hot Jupiter HD 209458b presented in Figure  15 of Caldas et al. (2019). We adopt the parametric P -T profile used in that model, which consists of a twolayer vertical profile, with a constant temperature in the deep atmosphere and a variable temperature from dayto nightside in the upper atmosphere. We produce a set of models with the same temperature structures and transition angles (β) as in Caldas et al. (2019), as shown in Figure 4. Our models show the same dependencies on temperature and transition angle as the previous work, with higher transit depths found for lower values of β.
While we have adopted the temperature profile described in Caldas et al. (2019) for the purposes of validating our model, we note that this prescription is simplified compared to our parametric P -T profile described in Section 2.3.3. Unlike our parametric profile, their parameterisation is not intended to be capable of matching temperature structures from GCMs (see Section 3.1).

RESULTS
In this section we demonstrate that our 3D parametric P -T profile can match azimuthally-averaged temperature structures from established GCMs. We subsequently explore the extent to which day-night temperature gradients and chemical inhomogeneity can affect resulting transmission spectra. We also perform 3D retrievals on a number of synthetic datasets, demonstrating Aura-3D's capabilities to constrain temperature profiles which vary across the terminator and to overcome biases in retrieved abundances that can arise from the assumption of a 1D temperature profile.

Comparison of parametric profile with GCM temperature profiles
We demonstrate that our 14-parameter prescription for the temperature profile is able to approximate a range of azimuthally-averaged temperature structures of several GCMs. Figure 5 shows temperature profiles from the terminator regions of several published 3D simulations of hot Jupiter atmospheres. We consider the models of HD 209458b and HD 189733b presented in Fortney et al. (2010), the model of HD 189733b from Blecic et al. (2017), and the Deep Hot Jupiter benchmark model from Deitrick et al. (2020). For each model we show the azimuthally-averaged temperature profiles across the terminator region of the atmosphere. Since only dayside profiles are presented in Blecic et al. (2017) we only show the temperature profiles at θ = 70 • , 80 • and 90 • for this case. To fit these models we assume a value of β = 40 • . For each model, we therefore fit the temperature profiles at θ = 70 • , 90 • and 110 • with our 14-parameter prescription, with T 3 and P 3 shared between the profiles, using a Trust Region Reflective least-squares algorithm (Sorenson 1982). The full temperature structure is subsequently calculated by linear interpolation between the nearest of those three profiles. Figure 5 includes profiles at θ = 80 • and 100 • in order to demonstrate the accuracy of the interpolated profiles. The full distribution of day-night temperature profiles for the models presented in Fortney et al. (2010) are shown in Figure 6.
Our parametric P -T prescription is able to fit all of the profiles shown in Figure 5 to a good degree of ac-  . For these models, we adopt the parametric temperature profile described in Caldas et al. (2019). We obtain the same effects as described in their paper: namely, that the transit depth decreases as β increases, and the divergence between models increases at larger wavelengths. Lower panels: Difference in ppm between spectra with the non-uniform temperature profiles described above and uniform spectra with averaged temperature profiles. Figure 5. Comparison between azimuthally-averaged P -T profiles from the terminator regions of a range of GCMs Blecic et al. 2017;Deitrick et al. 2020) and parametric profiles as described in Section 2.3. The GCM profiles are shown as solid lines, and parametric fits are shown as dashed lines. In each case the profiles at θ = 70 • , 90 • and 110 • are fit using a modified form of the parameterisation from Madhusudhan & Seager (2009), with the parameters P3 and T3 shared between each profile. The profiles at θ = 80 • and 100 • are calculated by linear interpolation between the nearest two profiles. For the Blecic et al. (2017) model, we only show profiles for θ ≤ 90 • since nightside profiles were not presented in that paper. curacy. The model of HD 209458b from ) is a good demonstration of why interpolation between three temperature profiles is preferred to two. In this case, the difference between the temperature profile at θ = 70 • and the profile at θ = 90 • is much greater than the difference between the profiles at θ = 90 • and θ = 110 • . If the parametric model only interpolated between the profiles at θ = 70 • and θ = 110 • , then the interpolated θ = 90 • would be much hotter than the profile from the GCM and would yield a worse fit.
For the Deep Hot Jupiter case, we did not fit our profile to the deep atmosphere (>10 bar) since this is below the region of the atmosphere probed in transmission, hence the deviation between the GCM and parametric profiles at high pressure. For atmospheric retrieval, the goal is to quickly generate a large number of model spectra which can be compared to observed data. Current retrievals require the generation of 10 6 forward models (Fortney et al. 2021). Generating a temperature structure from a GCM typically takes several days, rendering them unfeasible for retrievals. In contrast, a forward model using our parametric profile is computed in ∼0.6 seconds on a single core. Our parameterisation therefore enables the efficient exploration of a wide range of multidimensional temperature structures.
In order to confirm that our parametric P -T profile can be used to generate transmission spectra closely matching those generated with a full 3D temperature structure, we produce a model spectrum from 0.5-5.5µm using the Deep Hot Jupiter temperature structure from Deitrick et al. (2020), assuming a H/He-rich atmosphere with H 2 O as the sole absorbing species. We compare this to a spectrum generated using our parametric fit shown in Figure 5. The mean difference between the transit depths in the two spectra is 12.1 ppm, which is likely below the noise floor of most JWST instruments (Greene et al. 2016). In contrast, comparison between the full 3D model and a 1D best-fit to the globally-averaged temperature profile yields a difference of 65.5 ppm, indicating that our parameterisation provides a better match to the 3D spectrum than a 1D model. The differences between the spectra are more pronounced at longer wavelengths (see Figure 7).

Effects of a day-night temperature gradient on atmospheric spectra
We use our parametric temperature structure to explore how a day-night temperature gradient can lead to differences in resulting transmission spectra from 1D averaged models. The benefit of using the parametric model to investigate this is that a large number of spectra can be generated quickly for comparision, without the need to run a separate GCM for each temperature structure. In order to explore this effect we generate a grid of model hot Jupiter transmission spectra. We take the planetary and stellar parameters of the hot Jupiter HD 209458b (Stassun et al. 2017), and create models with a variety of temperature structures assuming solar H 2 O abundance. In the middle of the terminator region (θ = 90 • ), we vary the temperature at the top of the atmosphere in steps of 100 K from 1000 K to 1800 K, covering the range of typical hot Jupiters. We also vary the temperature contrast, ∆T , in steps of 5 K, from 0 K to 700 K. This temperature contrast is applied to both the day-and nightsides, so that T d = T t + ∆T and T n = T t − ∆T . We assume β = 40 • , meaning the temperature varies between θ = 70 • and θ = 110 • . The other parameters describing the temperature structure are held constant: the common temperature in the deepest layer of the atmosphere, set by T 3 , remains at T t + 800 K in all cases. The values of log P 1 , log P 2 and log P 3 are fixed to -0.9, -1.0 and 1.4 respectively. This means that none of the temperature profiles have a thermal inversion, since P 1 > P 2 .
We generate transmission spectra over a wavelength range of 0.5-15µm. This covers the wavelength range of several JWST instruments, including NIRSPEC FSS & BOTS, NIRISS SOSS, NIRCAM GRISM and MIRI LRS. Each spectrum is computed at a moderately high resolution (R = 5000). For each combination of T t and ∆T , we calculate the mean difference with the 1D spectrum and present our full results in Figure 7. An example with T t = 1500 K, ∆T = 500 K is shown in Figure  8. These tables show that for a wide range of temperature contrasts, the difference between a spectrum with a temperature gradient and one without is non-trivial. In general, the difference is more pronounced at higher temperatures and at higher temperature contrasts. For all models except T t = 1000 K, the difference between the spectra is >50 ppm for a temperature contrast of 500 K or more. This level of precision would be expected for many JWST transmission spectra of hot Jupiters, depending on factors such as observing time and the brightness of the host star. However, the difference between the spectra is not uniform with wavelength, with less of a noticeable difference at short wavelengths (see Figure 8. This suggests that for some instruments, such as NIRISS SOSS, the temperature gradient will have less of an impact on the observed spectrum, while the effect will be more important to consider when observing with instruments such as MIRI MRS and LRS. It is also important to note that while some HST Wide Field Camera 3 (WFC3) transmission spectra have also attained a level of precision of 50 ppm or better, these models indicate that WFC3 spectra will not be substantially affected by day-night temperature gradients since the difference between the spectra is small at the wavelengths probed by WFC3 (1.1-1.7 µm).

Effect of chemical inhomogeneity on transmission spectra
Since Aura-3D is capable of modelling atmospheres with inhomoegeneous chemistry, we explore the possibility for chemical variability between the day-and nightsides of a planet to affect resulting transmission spectra. Past studies have demonstrated that thermal dissociaton of molecules such as H 2 O on the daysides of ultra-hot Jupiters can substantially affect transmission spectra (Lothringer et al. 2018;Parmentier et al. 2018;Pluriel et al. 2020Pluriel et al. , 2022. In this section we consider a different chemical transition, namely the transition from a CH 4dominated to CO-dominated atmosphere. In chemical equilibrium, the abundances of key atmospheric species in a H 2 -rich atmosphere can vary strongly with temperature. At temperatures below ∼1200 K, CH 4 is expected to be a highly abundant molecule, whereas at higher temperatures CO is expected to dominate (Lodders & Fegley 2002;Moses et al. 2011;Madhusudhan 2012). If a planet's day-and nightsides are above and below this temperature respectively, this could lead to substantially different chemical abundances on either side of the planet.
In order to investigate this effect, we generate models of the hot Neptune HAT-P-11b (Bakos et al. 2010). Given this planet's equilibrium temperature of 878 K, it is possible that the day-and nightsides of the planet are above and below the transition temperature between CH 4 -dominated and CO-dominated chemistry, making the planet a useful test case to examine the effects of chemical inhomogeneity on transmission spectra. We consider two different sets of equilibrium chemical abundances: a CO-dominated regime, which we may expect to see on the dayside of the planet, and a CH 4dominated regime which may be expected on the nightside. The abundances for each regime are shown in Table 1. We choose a model P -T profile with a similar structure to the profiles used in Section 3.2. We take T t = 878 K, the equilibrium temperature of the planet, with ∆T = 500 K. The values of log P 1 , log P 2 and log P 3 are again fixed to -0.9, -1.0 and 1.4 respectively. For the chemistry we consider four different scenarios: (1) dayside abundances across the whole terminator region, (2) nightside abundances across the whole terminator, (3) inhomogeneous chemistry with dayside abundances for θ < 90 • and nightside abundances for θ > 90 • , and (4) averaged abundances across the whole terminator (the mean of the dayside and nightside abundances). We generate models from 0.5-15 µm, which covers the wavelength ranges of several JWST instruments, including NIRSpec, NIRISS NIRCam, and MIRI LRS, and also covers many of the prominent spectral features of the molecules considered (Madhusudhan 2019).
The resulting spectra are shown in Figure 9. We see that the model with nightside chemistry across the entire terminator has a higher transit depth overall than the model with dayside chemistry across the entire terminator. This is a result of the high abundances of H 2 O, CH 4 and NH 3 expected in chemical equilibrium at low temperatures, each of which has prominent absorption features at the wavelengths considered here. The model with averaged chemistry also differs substantially from the model with chemical abundances varying between the day-and nightsides of the planet across the wavelength range considered. The model with averaged chemistry has larger transit depths, since this model has large amounts of H 2 O, CH 4 and NH 3 on the dayside as well as the nightside, and the dayside has a larger scale height since it is at a higher temperature. In the model with variable chemistry, these molecules (particularly CH 4 and NH 3 ) are mostly present on the cooler side of the planet, leading to smaller features overall.
Whether such chemical inhomogeneities can be detected with JWST will need to be determined in future studies. For the remainder of this work we will focus on the possibility of constraining temperature inhomogeneities using JWST-quality observations.

Retrieval of a multidimensional temperature structure
Here we demonstrate the capability of our new retrieval framework to recover a multidimensional temperature structure. We create a parametric temperature Figure 6. Parametric fits to the 3D P -T profiles presented in Fortney et al. (2010) showing a slice at φ = 0 (through the north-south polar plane). As in Figure 5, the profiles at θ = 70 • , 90 • and 110 • follow our parametric prescription while the profiles all other values of θ are calculated by interpolation. For θ < 70 • and θ > 110 • the temperature remains constant, set by the temperatures at θ = 70 • and θ = 110 • respectively. Note that the inner portion of the planet (below Rp) is not to scale with the atmosphere. profile using the prescription described in Section 2.3.3. We use this temperature structure to produce a synthetic transmission spectrum, and carry out a retrieval of this spectrum using the same parameterisation. Our goal is to determine the extent to which the retrieval can constrain the input multi-dimensional temperature structure using JWST-quality data.
The input parameters used to generate the temperature structure are shown in the table embedded in Figure 14. The model interpolates between temperature profiles defined for the dayside (θ ≤ 70 • ), nightside (θ ≥ 110 • ) and the middle of the terminator (θ = 90 • ). We consider a cloud-free, solar composition atmosphere with uniform abundances of H 2 O, CH 4 , NH 3 , HCN, CO and CO 2 . Input chemical abundances are also shown in the table embedded in Figure 14.
To generate our synthetic data set we assume planetary and stellar properties that are similar to the hot Jupiter HAT-P-1b (Bakos et al. 2010). We choose this planet as an intermediate case that is representative of the majority of hot Jupiters that may be observed with JWST. HAT-P-1b has similar bulk properties to HD 209458b, and we find day-night temperature gradients to produce similar differences in the transmission spectrum to those shown in Figure 8. The forward model is first computed at a moderately high resolution (R = 3000, 0.5−5.5µm). Synthetic observations are then generated using PANDEXO (Batalha et al. 2017), assuming a single transit observed with NIRISS SOSS and another single transit observed with NIRSpec G395H. This provides a broad wavelength coverage with JWST that contains multiple spectral features of the chemical species  included in the forward model. The Nested Sampling retrieval uses 2000 live points.
The input data and retrieved spectrum are shown in Figure 10. We find excellent agreement between the input and retrieved spectra across the full wavelength range considered. Marginalised posterior distributions for each model parameter, as well as the retrieved median values with associated 1σ uncertainties, are shown in Figure 14. The estimates of most free parameters are consistent with input values to within 1σ uncertainty, and all retrieved parameters are consistent with input values to within 2σ.
Our input and retrieved temperature profiles are shown in Figure 11. Each of the dayside, nightside and terminator input temperature profiles almost entirely remain within the retrieved 1σ confidence intervals of the retrieved profiles, and lie completely within the 2σ confidence intervals. The uncertainty estimates for the temperature profiles in the model range from ∼100-300 K. We also compare 2D representations of the input and median retrieved temperature structures in Figure  0 100 . Differences between simulated hot Jupiter spectra with a 1D thermally averaged temperature profile and spectra with different day-night temperature gradients. For a range of terminator temperatures Tt (the temperature at the top of the atmosphere at θ = 90 • ) we compute spectra with temperature contrasts ∆T = 0-700 K (∆T = T d − Tt = Tt − Tn). In each case, the spectra are calculated at R = 5000 and a wavelength range of 1-15 µm and the mean difference between the two transit depths is calculated. The difference in transit depth is given in parts per million. The dashed red line at 50 ppm is an approximate indication of the expected precision that can readily be achieved by JWST observations of hot Jupiters.
12. We see that the temperature structures are very similar, highlighting the retrieval's ability to recover a detailed multidimensional temperature profile using high-quality data.
The retrieved abundances of all chemical species are shown in Figure 14. The abundances of H 2 O, CH 4 , CO and CO 2 are all constrained to within 0.5 dex and are consistent with input values to 1σ. We obtain upper limits on the abundances of NH 3 and HCN which are also consistent with input values, noting that these species have low abundances in the input model (log X NH3 = −6.3, log X HCN = −8.0).
This result demonstrates that using JWST observations of a nominal hot Jupiter, it is possible to obtain accurate constraints on a multidimensional temperature profile across the terminator region with reasonable precision. Furthermore, these observations can be used to obtain very precise abundance estimates for a number of chemical species.
As an additional validation test for our retrieval pipeline, we also conduct a retrieval of a synthetic spectrum for a planet with no day-night temperature gradient (i.e. a 1D P -T profile) using our 3D framework. The retrieved dayside, nightside and terminator temperature profiles are shown in Figure 13. We find that there is considerable overlap between each of the retrieved profiles, as would be expected for a planet with no temperature contrast. In this case, the lack of strong distinction between the three profiles suggests that this  Figure 8. Comparison of two theoretical transmission spectra for a HD 209458b-like planet (top right). The orange spectrum uses a parametric temperature structure with a day-night temperature contrast, shown by the blue, green and red temperature profiles in the top left panel, whereas the blue spectrum uses a 1D averaged structure shown by the dashed grey temperature profile. The bottom panel shows the difference between the two spectra as a function of wavelength. We find that the spectra differ by 78 ppm on average, enough that high-quality JWST observations of this planet would require day-night temperature variations to be considered in retrievals. spectrum is better suited to a 1D retrieval (see Section 4.1).

Retrieval of a synthetic spectrum from a GCM
In this section we present retrievals of a simulated JWST transmission spectrum of a hot Jupiter generated with a 3D temperature structure. We conduct retrievals using two different modelling paradigms: one using a 1D parametric P -T profile, and one using a temperature profile that is allowed to vary between the day-and nightsides of the planet.
The temperature structure of our input forward model is adapted from the open-source General Circulation Model THOR (Mendonça et al. 2016;Deitrick et al. 2020). We use THOR to reproduce the "Deep Hot Jupiter" scenario described in Deitrick et al. (2020). In order to create a temperature structure suitable for our forward model, we subsequently extrapolate the output temperature profiles from the GCM to lower pressures. In order to do this, we fit a parabolic P -T curve to each temperature profile in the longitude/latitude grid of outputs. The resulting 3D P -T structure is shown in Figure  16.
Having adapted the GCM temperature structure to be appropriate for our retrieval case study, we now use it to create a model JWST transmission spectrum. Our forward model in this instance is cloud-free and assumes chemical abundances of a 1×solar composition H 2 -dominated atmosphere (see Gandhi & Madhusudhan 2017), including H 2 O, CH 4 , NH 3 , HCN, CO and CO 2 .

HAT-P-11b
Dayside chemistry only Nightside chemistry only Variable day-night chemistry Averaged chemistry Figure 9. Synthetic spectra of the hot Neptune HAT-P-11b. Four different chemical compositions are considered: (1) 'dayside' chemistry everywhere, (2) 'nightside' chemistry everywhere, (3) variable chemistry between the day-and nightsides, and (4) averaged chemistry everywhere (see Table 1). Model 2 generally has a higher transit depth than model 1 due to higher abundances of H2O, CH4 and NH3 which are all prominent at these wavelengths. Model 3 generally lies between models 1 and 2, whereas model 4 lies closer to model 2, since the 'nightside' molecules are also found on the dayside of the planet, where the scale height is larger due to the higher temperatures.
We assume the same planetary and stellar properties for the planet HAT-P-1b as in Section 3.4.
We carry out two retrievals of our simulated JWST spectrum. The first such retrieval uses the 1D parametric P -T prescription as described in Madhusudhan & Seager (2009), while the second uses a temperature profile that varies across the terminator, but not with height (see Section 2.3.2). We fix β to 40 • , which is representative of the extent of the terminator region which affects a transmission spectrum ). We do not allow cases in which the dayside temperature is colder than the terminator temperature, or where the nightside temperature is hotter than the terminator.
A comparison of the two retrieval results is shown in Figure 15. We find that the retrieval assuming a 1D parametric P -T profile leads to some inaccuracies in the retrieved abundances using JWST-quality data. None of the retrieved abundances found using the 1D model are consistent with the input values to within 1σ. The only abundance retrieved accurately to ≤2σ is that of HCN, while the input abundance of NH 3 is within 3σ of the retrieved posterior. Note that these abundances are also the least well-constrained of the species considered. All other abundances fall more than 3σ away from the retrieved posteriors.
This finding is similar to that of Caldas et al. (2019), who found that when applying a 1D retrieval framework to a high-quality spectrum generated with a 3D input model, the true abundance would often be outside the error bars of the retrieval. They also find that the magnitude and direction of the error depends on the specific case; in this particular example, the retrieved abundances of most species are lower than the true values.
In contrast, for the retrieval in which the temperature is allowed to vary across the terminator, each of the retrieved chemical abundances is consistent within 1σ of the true values. Furthermore, the uncertainties on the parameter estimates in this case are generally larger than in the 1D retrieval. This suggests that the additional flexibility introduced by incorporation of a day-night temperature gradient allows the model to ex-

HAT-P-1b
Simulated data Binned best-fit Retrieved model Figure 10. Retrieval of a synthetic JWST spectrum of the hot Jupiter HAT-P-1b. The red points represent the simulated data, and the purple line shows the best-fit model spectrum, with 1σ and 2σ confidence intervals shown by the dark and light shaded regions. The input and retrieved models both incorporate a multidimensional parametric P -T profile. We find good agreement between the input data and the retrieved spectrum.  Figure 11. Input and retrieved multidimensional parametric temperature profiles of a hot Jupiter. The solid lines show the input dayside, nightside, and terminator temperature profiles used to generate the synthetic data. The dashed lines correspond to the median retrieved dayside, nightside and terminator temperature profiles respectively, with shaded regions corresponding to 1σ and 2σ confidence intervals. The retrieved temperature profiles are consistent with the input profiles, with almost all of the input temperature profiles contained within the retrieved 1σ confidence intervals.
plore a wider range of models that can explain the data at hand. This demonstrates that our prescription is able to mitigate against biases that can arise when only a 1D temperature profile is allowed in a retrieval of JWSTquality data.
We also carry out a retrieval of a simulated HST WFC3 spectrum of the same planet for comparison, using a 1D atmospheric model. In this case, the only well-constrained abundance is that of H 2 O, with a precision of 0.8 dex. The retrieved abundance of H 2 O is consistent with the input value, albeit with a larger uncertainty. This suggests that the constraints obtained from retrievals of HST data are not sufficiently precise to necessitate a multidimensional retrieval algorithm.
The retrieved temperature profiles acquired from each retrieval are shown on the right-hand panel of Figure 16. The retrieved 1D temperature profile exhibits a strong temperature gradient with pressure that is not found at the terminator of the input model. For the model including a day-night temperature gradient, the retrieved temperatures are broadly consistent with the dayside, nightside and terminator temperatures in the upper atmosphere. We note that both of these prescriptions are heavily simplified in comparison to a full 3D temperature structure. Therefore, it is difficult to gain complete information about the temperature structure of the planet's atmosphere from these retrievals alone, beyond the fact that a strong day-night temperature gradient appears to exist. However, this example demonstrates importance of having the capability to include a daynight temperature gradient in retrievals in order to en- Figure 12. Two-dimensional representations of the input and median retrieved parametric P -T profiles for HAT-P-1b. The temperature structure is computed by interpolating between the temperature profiles shown in Figure 11. Note that the inner portion of the planet (below Rp) is not to scale with the atmosphere.  Figure 13. Input and retrieved multidimensional parametric temperature profiles of a hot Jupiter with no day-night temperature contrast. The solid line shows the input temperature profile used to generate the synthetic data. The dashed lines correspond to the median retrieved dayside, nightside and terminator temperature profiles respectively, with shaded regions corresponding to 1σ and 2σ confidence intervals. The retrieved temperature profiles overlap significantly with each other and are broadly consistent with the input profile, suggesting that a multidimensional retrieval approach is not required in this case. Table 2. Description of priors for retrievals of simulated JWST transmission spectra of HAT-P-1b generated using a GCM. The temperature profile for the 1D retrieval uses the first six parameters (T0, α1, α2, P1, P2, P3) while the retrieval with a day-night temperature gradient uses the parameters T0, T 0,d , T0,n.  We introduce Aura-3D, a 3D modeling and retrieval framework for transmission spectra of exoplanetary atmospheres. The high data quality expected from upcoming JWST observations of exoplanet transmission spectra motivated the development of this new framework. Our new forward model enables calculation of transmission spectra with any 3D temperature structure, including GCM outputs, as well as computation of model spectra using parametric temperature profiles that include day-night inhomogeneities. Here we summarize the key functionalities of Aura-3D retrieval framework.

Parameter Lower Bound Upper Bound Prior
• The framework includes a 3D parametric temperature profile for use in atmospheric retrieval that can fit a wide range of GCM temperature structures. This parameterisation enables the computa- log X CO 2 Figure 15. Atmospheric retrieval of a synthetic hot Jupiter transmission spectrum generated using a 3D temperature structure. Top: Observations (blue) and retrieved model spectra (green and orange) for the two different model considerations. Shaded regions represent 1σ and 2σ confidence intervals. Bottom: Posterior distributions for the retrieved volume mixing ratios of each molecular species in the model. Input (equilibrium solar) values are shown by solid black lines. The model which does not include a day-night temperature contrast finds abundances that are not consistent with the input values. The model that does include a day-night temperature contrast is capable of inferring accurate abundances. The grey posterior distributions represent a retrieval of the HST WFC3 spectrum of the same model planet for comparison, using a 1D forward model. tion of transmission spectrum models with highly flexible, physically realistic temperature structures incorporating thermal inhomogeneities.
• The radiative transfer in 3D geometry is computationally efficient for retrievals, with computation time for a single model spectrum in 3D geometry 1s. We have demonstrated retrievals using a 3D geometry requiring over 10 6 model evaluations.
This is easily scalable to higher order depending on the requirement.
• When applied to synthetic JWST data of hot Jupiters, Aura-3D is able to accurately constrain separate dayside, nightside and terminator temperature profiles probed in transmission geometry, yielding meaningful information on day-night temperature contrasts. Right: Azimuthally-averaged temperature profiles from the forward model (solid lines), along with retrieved temperature profiles from the two retrievals. The retrieved 1D P -T profile is shown as a purple dotted line, with shaded 1σ and 2σ confidence intervals. The retrieved terminator (green), dayside (red) and nightside (light blue) temperatures from the retrieval incorporating day-night temperature variations are shown as dashed lines, also with shaded 1σ and 2σ confidence intervals. • Aura-3D can be used to accurately retrieve chemical abundances from transmission spectra whose temperature structures are calculated with 3D GCMs, even with JWST-quality data of hot Jupiters where 1D retrievals may lead to biased estimates.
We use Aura-3D to investigate constraints that can be placed on hot giant exoplanets using JWST-quality spectra, finding the following results: • For hot Jupiters with photospheric terminator temperatures 1100 K and temperature contrasts 500 K, the day-night temperature gradient can cause differences in transmission spectra that could be detectable with JWST.
• We demonstrate that it is possible to constrain a multidimensional temperature profile across the day-night terminator of a nominal hot Jupiter to a precision of ∼100-300 K using nominal JWSTquality data.
• We demonstrate that for JWST spectra of some hot Jupiters, a 1D retrieval can lead to biased abundance estimates, in agreement with previous studies. However, we find that in the case presented in this work, the introduction of a multidimensional temperature profile can overcome these biases and retrieve accurate abundances.
• For atmospheres where the terminator may transition between CO/CH 4 dominated composition, models with inhomogeneous composition may cause an observable difference to transmission spectra compared to those with globally-averaged compositions.

Applicability of 1D vs 3D retrievals
The developments presented in this work allow for retrievals with a much more complex forward model than has previously been possible. Our new framework will be extremely useful for analysing JWST-quality spectra of hot Jupiters with inhomogeneous terminators. While it is tempting to apply the most complex forward model available to all retrievals, care needs to be taken in ascertaining whether a 3D framework is necessary to explain a given observed spectrum. Furthermore, simply starting from the most complex possible model may not be computationally feasible. Although we have shown that retrievals with a 3D geometry can be carried out with high computational efficiency, the combination of 3D modelling alongside inhomogeneous clouds/hazes, stellar heterogeneity, large numbers of chemical species and other possible features such as inhomogeneous chemistry would lead to incredibly computationally expensive retrievals with an extremely high number of free parameters. It is therefore prudent to asses a priori what degree of complexity is required for a given spectrum, finding an optimal parameter set for analysing the data at hand.
While we have demonstrated that thermal inhomogeneities will lead to substantial deviations between 1D and 3D models for hot giant planets with large daynight temperature contrasts, we also find that for cooler planets the effect of thermal inhomogeneities is smaller and, therefore, may not be detectable even with JWSTquality data. When analysing future observations, several properties of the planet should be taken into account when considering the best P -T profile prescription, including the planetary radius, gravity and equilibrium temperature, all of which affect the amplitude of spectral features and therefore the degree to which effects such as day-night temperature gradients will be detectable. The extent of the terminator region being probed in transmission is also important to consider. Hotter planets tend to have larger opening angles, allowing for a greater variation in atmospheric conditions across the region probed in transmission (Wardenier et al. 2022). When considering retrievals of real observed data it will be important to determine whether the opening angle is sufficiently large to justify the inclusion of a 3D model.
As well as the planet being targeted, it is important to consider the quality of the observed data when determining the appropriate retrieval paradigm. This study has largely focused on high-quality JWST observations with a broad infrared wavelength coverage. We have found that deviations between 1D and 3D models are less significant when considering HST-quality data. Furthermore, other studies have demonstrated that the assumption of a 1D temperature profile will generally not lead to biases when analysing HST transmission spectra (Welbanks & Madhusudhan 2022) of most hot Jupiters. For future observations, the optimal retrieval paradigm will need to be chosen taking into account the wavelength range being probed as well as the resolution and precision of the spectrum.

Complementary observations using emission spectra and phase curves
In an ideal scenario, in order to obtain complete information about the structure and properties of an exoplanet's atmosphere, we would acquire its transmission and emission spectra as well as its phase curve. This combination of observations would allow for detailed characterisation of a planet's terminator region through transmission spectroscopy and of its dayside through emission spectroscopy, as well as providing insight into spatially-resolved properties through phase curve analysis. However, undertaking such an extensive observing campaign for a single planet would require very large amounts of telescope time, and if the planet does not exhibit strong spatial inhomogeneities, then time-consuming phase curve observations are unlikely to yield additional insights into the nature of its atmosphere.
The results from this work indicate that it should be possible to detect thermal inhomogeneities in certain exoplanet atmospheres using only transmission spectroscopy. We therefore suggest that an optimal observing strategy would be to first observe a planet's transmission spectrum to look for hints of day-night variability, possibly complementing this observation with an emission spectrum to constrain dayside properties and break certain degeneracies (see Section 4.3). Using these observations it will be possible to ascertain the degree to which spatial variations in temperature and chemistry permeate the planet's atmospheric structure, thus determining whether observations of the full phase-curve will lead to further meaningful results.

Future developments
The retrieval framework presented here enables 3D atmospheric characterisation using transmission spectra, a crucial step towards robust analysis of JWST data. However, there are additional atmospheric properties that will be important to consider when interpreting upcoming observations. For example, while inhomogeneous cloud/haze coverage has been considered in past works (e.g., Line & Parmentier 2016;MacDonald & Madhusudhan 2017;Welbanks & Madhusudhan 2021), the combined effect of patchy clouds/hazes and inhomogeneous temperature profiles is yet to be fully explored in a retrieval context.
The parametric P -T profile incorporated into Aura-3D has been sufficiently complex to analyse the syn-thetic JWST spectra considered in this study. Aura-3D has the capability to model thermal variations with both θ and φ, however only variations with θ have been explored to carry out the retrievals presented in this work. However, for certain targets it is possible that an even higher signal-to-noise will be achieved, either due to observations of atmospheres with a very large scale height or via observing multiple transits of a given target. In these cases it is possible that more complexity, including a φ-dependent temperature profile, will be required to accurately retrieve certain parameters. Our current framework has a built-in φ-dependence which has already been used to model patchy clouds (Pinhas et al. 2019), and has the functionality to include a parametric temperature profile with variations in both φ and θ.
It is also important to assess the impact of inhomogeneous chemical abundances on retrievals of transmission spectra. This topic has been explored for ultrahot Jupiters ) but may also be relevant for cooler planets, as discussed in this work. While Aura-3D is capable of modelling planets with different day-and nightside chemical compositions, further work is needed to ascertain the feasibility of retrieving inhomogenous chemical abundances directly from a transmission spectrum. A retrieval allowing for independent dayside and nightside chemical abundances for each species would be subject to very strong degeneracies, since low dayside abundances could be compensated for by high nightside abundances and vice-versa. It would therefore be important to limit the allowed parameter space to physically realistic dayside and nightside abundances motivated by detailed self-consistent models. As discussed above, combining transmission spectra with observations of other regions of a planet's atmosphere, such as emission spectra, could also be useful in breaking this degeneracy.
Ultimately, we should strive to incorporate the most sophisticated forward models possible into retrieval algorithms, so that they can be invoked if necessary to explain observations. Any additional increases in model complexity must be accompanied by concurrent efforts to make these models efficient for retrievals, as we have done in the present work. Efforts are underway to improve the computational efficiency of forward models and retrievals through means such as GPU implementations (e.g., Malik et al. 2017) and Machine Learning (e.g., Márquez-Neila et al. 2018;Nixon & Madhusudhan 2020). These avenues for improvement may enable further steps towards reducing other simplifying assumptions that still remain in place in current retrieval frameworks.
The developments presented in this work represent an important step forwards in our ability to identify and constrain multidimensional effects in exoplanet atmospheres. This opens the field to comprehensive atmospheric characterisation of exoplanets with large spectral features that will be accessible with transmission spectroscopy in the JWST era. Our retrieval framework will be an important tool as we work to uncover exciting new insights into the nature of exoplanet atmospheres using upcoming observations of unprecedented detail.