TRIDENT: A Rapid 3D Radiative-transfer Model for Exoplanet Transmission Spectra

A cylindrical coordinate system is a natural choice for radiative transfer, given the symmetry of the problem (Fortney et al. 2010). Consider a Cartesian coordinate system (x, y, z), with x encoding the distance from the day–night boundary (increasing toward the observer), y aligned with the direction of orbital motion, and z pointing to the planetary north pole. From the perspective of a distant observer, y–z defines the terminator plane, located at x = 0. We define a polar coordinate system (b, ϕ) in the terminator plane (Figure 2, left panel), with b the radial coordinate and ϕ the azimuthal angle with respect to the north pole. The advantage of this coordinate system is that a ray has constant values of b and ϕ in the geometric limit (i.e., rays travel in straight lines when refraction and scattering is negligible ² ). The only changing coordinate for a ray is then x. The trajectory of a ray is thus readily specified by a cylindrical coordinate system (b, ϕ, x).

Atmospheric properties, meanwhile, are more naturally described by a spherical coordinate system. While spherical coordinate systems are often defined about the north pole (i.e., longitude and latitude), the day–night irradiation along the x-axis leads us to a different choice. We consider a coordinate system defined by three quantities: the radial distance from the planet's center, r; the azimuthal angle in the terminator plane, ϕ (identical to the polar angle used in radiative transfer); and the local zenith angle with respect to the terminator plane, θ. The coordinate system (r, θ, ϕ) is closely related to a spherical coordinate system about the x-axis, differing only from our choice of θ being with respect to the terminator plane instead of the x-axis. We stress that θ is defined with respect to the terminator plane vertical for a specific ray, $z^{\prime}$, not the north pole axis, z (see Figure 2). The axis $z^{\prime}$ is aligned with the impact-parameter vector. Our (r, θ, ϕ) coordinate system has the advantage of reducing the geometrical complexity of 3D radiative transfer to a series of 2D problems. For a given ray at (b, ϕ, x), the corresponding atmospheric coordinates (r, θ, ϕ) are simply

$$\begin{eqnarray}&&r=\sqrt{{b}^{2}+{x}^{2}}\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&\theta ={\tan }^{-1}(x/b),\end{eqnarray} \tag{ 8 }$$

with ϕ the same in both systems. For a given ϕ, one then traces a series of rays, with different impact parameters, through a 2D day–night slice. The final 3D solution is then just a sum of 2D radiative-transfer solutions.

We can now write the general transmission spectrum expression (Equation (1)) in our cylindrical coordinate system. In what follows, we assume a uniform stellar photosphere (_λ,het = 1) and negligible nightside emission (ψ_λ,night = 1) to focus on 3D effects arising from atmospheric transmission. We consider only full transits (see Section 6.2 for a discussion of ingress/egress and grazing transits), such that ${A}_{{\rm{p}}\ \,(\mathrm{overlap})}=\pi {R}_{{\rm{p}},\,\mathrm{top}}^{2}$. Given these assumptions, a transmission spectrum can be written as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda }=\displaystyle \frac{{R}_{{\rm{p}},\,\mathrm{top}}^{2}-\tfrac{1}{\pi }{\int }_{-\pi }^{\pi }{\int }_{0}^{{R}_{{\rm{p}},\,\mathrm{top}}}{\delta }_{\mathrm{ray}* }(b,\phi )\,{{ \mathcal T }}_{\lambda }(b,\phi )\,b\,{db}\,d\phi }{{R}_{* }^{2}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{ \mathcal T }}_{\lambda }(b,\phi )={e}^{-{\tau }_{\lambda ,\mathrm{path}}(b,\phi )}\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{path}}(b,\phi )={\int }_{0}^{\infty }{\kappa }_{\lambda }(s)\,{ds}.\end{eqnarray} \tag{ 11 }$$

In the final expression, the path optical depth is expressed as an integral of the extinction coefficient, κ_λ (see Section 3.3), over the path, s, taken by a ray through the atmosphere. This functional form assumes no scattering into each beam and negligible forward scattering, with all scattering opacity treated as a loss equivalent to absorption. ³ If one further assumes refraction to be negligible, the path optical depth is identical to the slant optical depth:

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{slant}}(b,\phi )={\int }_{-\infty }^{\infty }{\kappa }_{\lambda }(b,\phi ,x)\,{dx}.\end{eqnarray} \tag{ 12 }$$

The evaluation of Equation (12) determines how an atmosphere shapes a transmission spectrum. We describe the computation of this integral in the following section.

Finally, we demonstrate that Equation (9) is equivalent to a more familiar expression under two further simplifying assumptions. First, if one neglects ray deflection due to refraction or scattering, δ_ray* = 1 for all impact parameters and azimuthal angles. Second, if one assumes the planetary atmosphere has uniform properties (temperature, composition, etc.) at all longitudes and latitudes, varying only in the radial direction (i.e., the 1D atmosphere approximation), then the azimuthal integral in Equation (9) evaluates to 2π. With these approximations, a 1D transmission spectrum is given by

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda }\approx \displaystyle \frac{{R}_{{\rm{p}},\,\mathrm{top}}^{2}-2{\int }_{0}^{{R}_{{\rm{p}},\,\mathrm{top}}}b\,{e}^{-{\tau }_{\lambda ,\mathrm{slant}}(b)}\,{db}}{{R}_{* }^{2}}.\end{eqnarray} \tag{ 13 }$$

Expressions similar to Equation (13) are widely seen in the transmission spectra literature (e.g., Brown 2001; Tinetti et al. 2012), forming the basis for many commonly used radiative-transfer models. In this paper, we instead employ the more general expression given by Equation (9). We now show how to evaluate the optical depth for 3D, nonuniform exoplanet atmospheres.

The computation of path optical depths becomes nontrivial in 3D due to the need to interrelate two coordinate systems. In other words, the atmospheric properties (and hence the extinction coefficient) are defined in the spherical coordinate system (r, θ, ϕ), while the path optical depth is evaluated in the cylindrical coordinate system (b, ϕ, x) (see Figure 2). Here, we first illustrate how to solve this problem for 1D atmospheres (with properties varying only with r), following an approach introduced by Robinson (2017), before generalizing the method for 3D atmospheres.

In 1D, the path optical depth is given by

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{path}}(b)={\int }_{0}^{\infty }{\kappa }_{\lambda }(s)\,{ds}={\int }_{{r}_{\min }(b)}^{\infty }{\kappa }_{\lambda }(r)\,{{ \mathcal P }}_{1{\rm{D}}}(b,r)\,{dr},\end{eqnarray} \tag{ 14 }$$

where ${r}_{\min }$ is the minimum radius ⁴ probed by a ray with impact parameter b, while ${{ \mathcal P }}_{1{\rm{D}}}(b,r)$ is the 1D path distribution (Robinson 2017). The 1D path distribution is defined such that the distance traveled by a ray with impact parameter b passing through an atmospheric layer extending from r to r + dr is

$$\begin{eqnarray}&&{ds}={{ \mathcal P }}_{1{\rm{D}}}(b,r)\,{dr}.\end{eqnarray} \tag{ 15 }$$

Recognizing that the product κ_λ (r) dr in Equation (14) is just the differential vertical optical depth across a layer between r and r + dr, we can thus write

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{path}}(b)=\int {{ \mathcal P }}_{1{\rm{D}}}(b,r)\,d{\tau }_{\lambda ,\mathrm{vert}},\end{eqnarray} \tag{ 16 }$$

where the integral is evaluated between the lowest and highest layers in the atmosphere, and we leave the r dependence explicit to denote the layer corresponding to each differential vertical optical depth. We thus see that the path distribution acts to map the vertical optical depth in an atmosphere to the path optical depth. For a continuous atmosphere, the 1D path distribution can be analytically calculated in the geometric limit. From the geometry in Figure 2, we have

$$\begin{eqnarray}&&{{ \mathcal P }}_{1{\rm{D}}}(b,r)=\displaystyle \frac{{ds}}{{dr}}=2\displaystyle \frac{{dx}}{{dr}}=\displaystyle \frac{2\,r}{\sqrt{{r}^{2}-{b}^{2}}},\end{eqnarray} \tag{ 17 }$$

where the factor of 2 arises from the equal distance traveled on either side of the terminator plane in 1D models.

To numerically compute transmission spectra, atmospheres must be discretized. Consider an atmosphere comprised of distinct layers, r_l , impinged by rays spanning a grid of impact parameters, b_i . The path optical depth for the ith impact parameter is then given by

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{path},i}=\sum _{l=1}^{{N}_{\mathrm{layer}}}{{ \mathcal P }}_{1{\rm{D}},{il}}\,{\rm{\Delta }}{\tau }_{\lambda ,\mathrm{vert},l},\end{eqnarray} \tag{ 18 }$$

where Δτ_λ,vert,l = κ_λ,l Δr_l and the 1D path distribution becomes a matrix, ${{\boldsymbol{ \mathcal Q }}}_{1{\rm{D}}}$. In the geometric limit, the elements of the 1D path distribution matrix are given by (see Appendix B)

$$\begin{eqnarray}{{ \mathcal P }}_{1{\rm{D}},{il}}=\left\{\begin{array}{cl}0, & \unicode{x02460} \\ \displaystyle \frac{2}{{\rm{\Delta }}{r}_{l}}\left(\Space{0ex}{2.25ex}{0ex}\sqrt{{r}_{\mathrm{up},l}^{2}-{b}_{i}^{2}}\Space{0ex}{2.25ex}{0ex}\right), & \unicode{x02461} \\ \displaystyle \frac{2}{{\rm{\Delta }}{r}_{l}}\left(\Space{0ex}{2.25ex}{0ex}\sqrt{{r}_{\mathrm{up},l}^{2}-{b}_{i}^{2}}-\sqrt{{r}_{\mathrm{low},l}^{2}-{b}_{i}^{2}}\Space{0ex}{2.25ex}{0ex}\right), & \unicode{x02462} \end{array}\right.\end{eqnarray} \tag{ 19 }$$

where r_up,l and r_low,l are the upper and lower boundaries of the layer centered on r_l . The conditions of applicability in Equation (19) are

$$\begin{eqnarray}&&\begin{array}{l}\unicode{x02460} :\ {r}_{\mathrm{up},l}\leqslant {b}_{i}\\ \unicode{x02461} :\ {r}_{\mathrm{low},l}\lt {b}_{i}\lt {r}_{\mathrm{up},l}\\ \unicode{x02462} :\ {r}_{\mathrm{low},l}\geqslant {b}_{i}.\end{array}\end{eqnarray} \tag{ 20 }$$

Equation (18) can then be conveniently recast in vector notation to compute a vector of path optical depths:

$$\begin{eqnarray}&&{{\boldsymbol{\tau }}}_{\lambda ,\mathrm{path}}={{\boldsymbol{ \mathcal Q }}}_{1{\rm{D}}}\cdot {\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,\mathrm{vert}}.\end{eqnarray} \tag{ 21 }$$

As noted in Robinson (2017), the main advantage of computing path optical depths in this manner is that the path distribution matrix is wavelength independent. ⁵ The path distribution matrix then needs only be computed once for a given atmosphere, significantly reducing the computational burden of transmission spectra calculations. Given the need for quick spectral computations for atmospheric retrievals (which typically require ≳10⁵ model spectra), the path distribution framework is an excellent foundation for retrieval forward models.

For 3D atmospheres, the path distribution technique must be generalized. The changing temperature across the terminator (see Figure 2) causes the ray path length between two altitudes to become a function of θ, breaking the symmetry between the dayside and nightside. We thus express the path optical depth as

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{path}}(b,\phi )={\int }_{-\tfrac{\pi }{2}}^{\tfrac{\pi }{2}}{\int }_{{r}_{\min }(b)}^{r(b,\theta )}{\kappa }_{\lambda }(r,\theta ,\phi )\,{ \mathcal P }(b,\phi ,\theta ,r)\,{dr}\,d\theta ,\end{eqnarray} \tag{ 22 }$$

where r(b, θ) is the radial distance of a ray ⁶ with impact parameter b at zenith angle θ. Here, we introduce the 3D generalization of the path distribution, ${ \mathcal P }(b,\phi ,\theta ,r)$. We define the 3D path distribution such that the distance traveled by a ray (with impact parameter b and azimuthal angle ϕ) passing through an atmospheric region extending from θ to θ + d θ and a layer extending from r to r + dr is

$$\begin{eqnarray}&&{ds}={ \mathcal P }(b,\phi ,\theta ,r)\,{dr}\,d\theta .\end{eqnarray} \tag{ 23 }$$

With this definition, the 3D path optical depth is

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{path}}(b,\phi )={\int }_{-\tfrac{\pi }{2}}^{\tfrac{\pi }{2}}\int { \mathcal P }(b,\phi ,\theta ,r)\,d{\tau }_{\lambda ,\mathrm{vert}}\,d\theta .\end{eqnarray} \tag{ 24 }$$

Consider now a 3D discretized atmosphere. We divide the atmosphere into columns specified by their azimuthal angle, ϕ_j , and zenith angle, θ_k . Since the temperature profile can be different in each column (see Section 3.2), every column has a distinct radial array, r_jkl , with "l" denoting the layer in each column. This atmosphere is illuminated by a grid of impact parameters, b_i . For such a discretized atmosphere, Equation (24) becomes

$$\begin{eqnarray}&&{\tau }_{\lambda ,\mathrm{path},{ij}}=\sum _{k=1}^{{N}_{\mathrm{zenith}}}\sum _{l=1}^{{N}_{\mathrm{layer}}}{{ \mathcal P }}_{{ijkl}}\,{\rm{\Delta }}{\tau }_{\lambda ,\mathrm{vert},{jkl}}\,{\rm{\Delta }}{\theta }_{k},\end{eqnarray} \tag{ 25 }$$

where Δτ_λ,vert,jkl = κ_λ,jkl Δr_jkl . The path distribution is now a 4D tensor, ${\boldsymbol{ \mathcal Q }}$, with elements given, in the geometric limit, by (see Appendix B)

$$\begin{eqnarray}{{ \mathcal P }}_{{ijkl}}=\left\{\begin{array}{cl}0, & \unicode{x02460} \\ \displaystyle \frac{1}{{\rm{\Delta }}{r}_{{jkl}}{\rm{\Delta }}{\theta }_{k}}\left(\Space{0ex}{2.25ex}{0ex}\sqrt{{r}_{\mathrm{up},{jkl}}^{2}-{b}_{i}^{2}}-\sqrt{{r}_{\min ,{ijk}}^{2}-{b}_{i}^{2}}\Space{0ex}{2.25ex}{0ex}\right), & \unicode{x02461} \\ \displaystyle \frac{1}{{\rm{\Delta }}{r}_{{jkl}}{\rm{\Delta }}{\theta }_{k}}\left(\Space{0ex}{2.25ex}{0ex}\sqrt{{r}_{\max ,{ijk}}^{2}-{b}_{i}^{2}}-\sqrt{{r}_{\min ,{ijk}}^{2}-{b}_{i}^{2}}\Space{0ex}{2.25ex}{0ex}\right), & \unicode{x02462} \\ \displaystyle \frac{1}{{\rm{\Delta }}{r}_{{jkl}}{\rm{\Delta }}{\theta }_{k}}\left(\Space{0ex}{2.25ex}{0ex}\sqrt{{r}_{\mathrm{up},{jkl}}^{2}-{b}_{i}^{2}}-\sqrt{{r}_{\mathrm{low},{jkl}}^{2}-{b}_{i}^{2}}\Space{0ex}{2.25ex}{0ex}\right), & \unicode{x02463} \\ \displaystyle \frac{1}{{\rm{\Delta }}{r}_{{jkl}}{\rm{\Delta }}{\theta }_{k}}\left(\Space{0ex}{2.25ex}{0ex}\sqrt{{r}_{\max ,{ijk}}^{2}-{b}_{i}^{2}}-\sqrt{{r}_{\mathrm{low},{jkl}}^{2}-{b}_{i}^{2}}\Space{0ex}{2.25ex}{0ex}\right), & \unicode{x02464}\\ 0, & \unicode{x02465}\end{array}\right.\end{eqnarray} \tag{ 26 }$$

where the conditions of applicability are

$$\begin{eqnarray}\begin{array}{rcl}\unicode{x02460} :\ {r}_{\mathrm{low},{jkl}} & \leqslant & {r}_{\min ,{ijk}}\\ & & {r}_{\mathrm{up},{jkl}}\leqslant {r}_{\min ,{ijk}}\\ \unicode{x02461} :\ {r}_{\mathrm{low},{jkl}} & \leqslant & {r}_{\min ,{ijk}}\\ & & {r}_{\min ,{ijk}}\lt {r}_{\mathrm{up},{jkl}}\lt {r}_{\max ,{ijk}}\\ \unicode{x02462} :\ {r}_{\mathrm{low},{jkl}} & \leqslant & {r}_{\min ,{ijk}}\\ & & {r}_{\mathrm{up},{jkl}}\geqslant {r}_{\max ,{ijk}}\\ \unicode{x02463} :\ {r}_{\min ,{ijk}} & \lt & {r}_{\mathrm{low},{jkl}}\lt {r}_{\max ,{ijk}}\\ & & {r}_{\min ,{ijk}}\lt {r}_{\mathrm{up},{jkl}}\lt {r}_{\max ,{ijk}}\\ \unicode{x02464}:\ {r}_{\min ,{ijk}} & \lt & {r}_{\mathrm{low},{jkl}}\lt {r}_{\max ,{ijk}}\\ & & {r}_{\mathrm{up},{jkl}}\geqslant {r}_{\max ,{ijk}}\\ \unicode{x02465}:\ {r}_{\mathrm{low},{jkl}} & \geqslant & {r}_{\max ,{ijk}}\\ & & {r}_{\mathrm{up},{jkl}}\geqslant {r}_{\max ,{ijk}}.\end{array}\end{eqnarray} \tag{ 27 }$$

In the expressions above, r_up,jkl and r_low,jkl are the upper and lower boundaries of the lth layer in column jk, while ${r}_{\min ,{ijk}}$ and ${r}_{\max ,{ijk}}$ are the minimum and maximum radial coordinates encountered by a ray with impact parameter b_i in column jk. They are given by

$$\begin{eqnarray}&&{r}_{\min ,{ijk}}=\min \left({R}_{{\rm{p}},\,\mathrm{top},{jk}},\displaystyle \frac{{b}_{i}}{\cos {\theta }_{\min ,k}}\right),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{r}_{\max ,{ijk}}=\min \left({R}_{{\rm{p}},\,\mathrm{top},{jk}},\displaystyle \frac{{b}_{i}}{\cos {\theta }_{\max ,k}}\right),\end{eqnarray} \tag{ 29 }$$

where R_{p, top,jk} is the radial coordinate at the top of column jk, ${\theta }_{\min ,k}$ is the minimum zenith angle of the column, and ${\theta }_{\max ,k}$ is the maximum zenith angle of the column. The conditions in Equation (27) account for all possible ways in which a ray can traverse zenith slices before leaving the atmosphere (see Appendix B).

We efficiently numerically compute 3D path optical depths via matrix operations. By defining an ancillary tensor $\tilde{{\boldsymbol{ \mathcal Q }}}$, with elements given by ${{ \mathcal P }}_{{ijkl}}{\rm{\Delta }}{\theta }_{k}$, we can then succinctly write Equation (25) as

$$\begin{eqnarray}&&{{\boldsymbol{\tau }}}_{\lambda ,\mathrm{path},j}={\tilde{{\boldsymbol{ \mathcal Q }}}}_{j}\,:\,{\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,\mathrm{vert},j},\end{eqnarray} \tag{ 30 }$$

where ":" denotes the double-dot product. In other words, the path distribution tensor for sector "j" is contracted with the vertical optical depth matrix for sector "j" over their common indices "lk". The result is a vector of size N_b containing the slant optical depths for all the impact parameters in sector "j". Computing this double-dot product for all azimuthal sectors yields the path optical depth matrix, τ _λ,path (dimensions N_b × N_sectors). This equation is the 3D generalization of Equation (21). We stress that the functional form in Equation (30) holds for models both with and without refraction. However, the analytical expressions we derive for the path distribution tensor (Equation (26)) are specific to the geometric limit. ⁷ With the path optical depth determined, one can readily compute 3D transmission spectra.

TRIDENT computes transmission spectra in a discretized representation of the geometry in Figure 2. We divide the atmosphere into N_sector azimuthal sectors in the terminator plane, N_zenith zenith slices from the dayside to nightside, and N_layer layers in each column. We assume the day–night transition region has a zenith angular width of β, with atmospheric properties linearly varying over the transition region (see Section 3.2 and Figure 3). Similarly, we assume a morning–evening transition region with an azimuthal angular width of α. Outside the transition regions (i.e., ∣θ∣ > β/2 and ∣ϕ∣ > α/2) the atmosphere varies only with altitude. By convention, the dayside is the first zenith slice (k = 1) and the nightside is the last (k = N_zenith); the evening terminator is the first azimuthal sector (j = 1) and the morning terminator is the last (j = N_sector). The day–night and morning–evening transition regions thus span N_zenith − 2 zenith slices and N_sector − 2 azimuthal sectors, respectively, with uniform angular spacing over the transitions. We further assume north–south symmetry for the background atmosphere (a common assumption, especially in cases with zero obliquity), such that sectors need only be defined for the northern hemisphere.

Zoom In Zoom Out Reset image size

TRIDENT computes transmission spectra by solving the discretized representation of Equation (9). For a collection of rays (b, ϕ), with dimensions N_b × N_sector, the transmission spectrum is given by

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda }=\displaystyle \frac{{R}_{{\rm{p}},\,\mathrm{top}}^{2}-\tfrac{1}{\pi }\sum _{j=1}^{{N}_{\mathrm{sector}}}\sum _{i=1}^{{N}_{b}}{\delta }_{\mathrm{ray}* ,{ij}}\,{{ \mathcal T }}_{\lambda ,{ij}}\,{b}_{i}\,{\rm{\Delta }}{b}_{i}\,{\rm{\Delta }}{\phi }_{j}}{{R}_{* }^{2}},\end{eqnarray} \tag{ 31 }$$

where we take R_{p, top} to be the maximum radial extent of the dayside. ⁸ We take the impact-parameter array to coincide with the lower-layer boundaries of the column with greatest radial extent (i.e., b_i =${\max }_{{jk}}{r}_{\mathrm{low},{jkl}}$). Defining the projected area matrix ⁹ for the atmosphere, ΔA _atm, via ΔA_atm,ij = b_i Δb_i Δϕ_j δ_ray*,ij we finally have

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda }=\displaystyle \frac{{R}_{{\rm{p}},\,\mathrm{top}}^{2}-\tfrac{1}{\pi }\left({{\boldsymbol{ \mathcal T }}}_{\lambda }:{{\boldsymbol{\Delta A}}}_{\mathrm{atm}}\right)}{{R}_{* }^{2}},\end{eqnarray} \tag{ 32 }$$

where the transmission matrix is

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal T }}}_{\lambda }={e}^{-{{\boldsymbol{\tau }}}_{\lambda ,\mathrm{path}}}.\end{eqnarray} \tag{ 33 }$$

Although deceptively simple, Equation (32) provides the practical basis for rapid 3D radiative transfer through nonuniform exoplanet atmospheres. In the geometric limit, which we adopt in the following sections, the optical depth matrix is given by Equation (30). With our radiative-transfer framework established, we proceed to describe the construction of 3D model atmospheres.

Model atmosphere temperature structures are generally constructed by one of two approaches. Forward-modeling studies typically calculate pressure–temperature (P-T) profiles self-consistently with various physical principles (e.g., radiative–convective equilibrium) via an iterative procedure (e.g., Sudarsky et al. 2003; Mollière et al. 2015). Retrieval studies, however, routinely assume an isotherm or a parametric P-T profile to rapidly explore a wide parameter space of potential temperature structures (e.g., Madhusudhan & Seager 2009; Line et al. 2013). Here, we introduce a temperature-structure parameterization appropriate for 3D retrievals of exoplanet transmission spectra.

Multidimensional retrievals require a temperature-field prescription, T(P, ϕ, θ), which must satisfy several requirements. First, it must be sufficiently flexible to capture realistic nonuniform temperature structures, potentially varying on the order of ∼1000 K between different regions of hot Jupiter atmospheres. Second, it must account for the vertical dependence of temperature with pressure, since transmission spectra are sensitive to temperature gradients (Rocchetto et al. 2016; Heng & Kitzmann 2017). Third, the temperature profiles in each column should converge in the deep atmosphere where external irradiation has no influence. Finally, the temperature field should employ the minimal set of free parameters required to capture the information content of 3D transmission spectra. This latter requirement renders common 1D P-T profiles (e.g., Madhusudhan & Seager 2009; Line et al. 2013) suboptimal ¹⁰ for multidimensional retrievals. Instead, we propose a new parameterization that represents a simple construction appropriate for 3D transmission spectra retrievals.

We propose a four-parameter temperature field. Consider a slice through the terminator plane, with the top-of-atmosphere temperature varying linearly between a morning terminator temperature, T_morn, and an evening terminator temperature, T_even, over a transition region with an angular width α. This can be expressed via

$$\begin{eqnarray}{T}_{\mathrm{term}}(\phi )=\left\{\begin{array}{cl}{T}_{\mathrm{even}}, & \unicode{x02460} \\ {\overline{T}}_{\mathrm{term}}-\left(\displaystyle \frac{\phi }{\alpha /2}\right)\displaystyle \frac{{\rm{\Delta }}{T}_{\mathrm{term}}}{2}, & \unicode{x02461} \\ {T}_{\mathrm{morn}}, & \unicode{x02462} \end{array}\right.\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray*}&&\begin{array}{l}\unicode{x02460} :\ \phi \leqslant -\alpha /2\\ \unicode{x02461} :\ -\alpha /2\lt \phi \lt \alpha /2\\ \unicode{x02462} :\ \phi \geqslant \alpha /2,\end{array}\end{eqnarray*}$$

where ${\overline{T}}_{\mathrm{term}}\equiv \tfrac{1}{2}({T}_{\mathrm{even}}+{T}_{\mathrm{morn}})$ and ΔT_term ≡ T_even − T_morn. Once the terminator temperature, T_term(ϕ), is determined for a given azimuthal angle, we can account for the day–night variation. Assuming a linear transition between the dayside and nightside temperatures over a region of angular width β, we have

$$\begin{eqnarray}{T}_{\mathrm{high}}(\phi ,\theta )=\left\{\begin{array}{cl}{T}_{\mathrm{day}}(\phi ), & \unicode{x02460} \\ {\overline{T}}_{\mathrm{term}}(\phi )-\left(\displaystyle \frac{\theta }{\beta /2}\right)\displaystyle \frac{{\rm{\Delta }}{T}_{\mathrm{DN}}}{2}, & \unicode{x02461} \\ {T}_{\mathrm{night}}(\phi ), & \unicode{x02462} \end{array}\right.\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray*}&&\begin{array}{l}\unicode{x02460} :\ \theta \leqslant -\beta /2\\ \unicode{x02461} :\ -\beta /2\lt \theta \lt \beta /2\\ \unicode{x02462} :\ \theta \geqslant \beta /2,\end{array}\end{eqnarray*}$$

where we assume the day–night temperature contrast, ${\rm{\Delta }}{T}_{\mathrm{DN}}\equiv {T}_{\mathrm{day}}-{T}_{\mathrm{night}}$, is constant for all ϕ. Therefore, the dayside and nightside temperatures are automatically specified by ${T}_{\mathrm{day}}(\phi )={T}_{\mathrm{term}}(\phi )+\tfrac{1}{2}{\rm{\Delta }}{T}_{\mathrm{DN}}$ and ${T}_{\mathrm{night}}(\phi )={T}_{\mathrm{term}}(\phi )-\tfrac{1}{2}{\rm{\Delta }}{T}_{\mathrm{DN}}$. We note that if morning–evening gradients are neglected (ΔT_term = 0), Equation (35) is equivalent to the day–night gradient prescription studied by Caldas et al. (2019). Finally, each column has a vertical temperature gradient, assumed linear in log-pressure, running from T_high(ϕ, θ) to T_deep:

$$\begin{eqnarray}T(P,\phi ,\theta )=\left\{\begin{array}{cl}{T}_{\mathrm{high}}(\phi ,\theta ), & \unicode{x02460} \\ {T}_{\mathrm{high}}(\phi ,\theta )+\left[\displaystyle \frac{{T}_{\mathrm{deep}}-{T}_{\mathrm{high}}(\phi ,\theta )}{\mathrm{log}({P}_{\mathrm{deep}}/{P}_{\mathrm{high}})}\right]\mathrm{log}\left(\displaystyle \frac{P}{{P}_{\mathrm{high}}}\right), & \unicode{x02461} \\ {T}_{\mathrm{deep}}, & \unicode{x02462} \end{array}\right.\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray*}&&\begin{array}{l}\unicode{x02460} :\ P\leqslant {P}_{\mathrm{high}}\\ \unicode{x02461} :\ {P}_{\mathrm{high}}\lt P\lt {P}_{\mathrm{deep}}\\ \unicode{x02462} :\ P\geqslant {P}_{\mathrm{deep}},\end{array}\end{eqnarray*}$$

where P_high denotes where higher altitudes are fixed to T_high(ϕ, θ), while P_deep specifies where lower altitudes are at T_deep (common for all columns). We set P_high = 10⁻⁵ bar and P_deep = 10 bar, serving as anchors outside the pressure range probed by low-resolution transmission spectra. Finally, TRIDENT slightly Gaussian smooths each profile (with a standard deviation of three layers), avoiding a discontinuous temperature gradient at P_high and P_deep (not applying this smoothing negligibly alters transmission spectra). Our temperature-field prescription has the strength of being fully specified by just four quantities: ${\overline{T}}_{\mathrm{term}}$, ΔT_term, ${\rm{\Delta }}{T}_{\mathrm{DN}}$, and T_deep. This simplicity avoids overparameterization of the temperature field in multidimensional retrievals.

We visualize our temperature-field prescription for a typical ultra-hot Jupiter in Figure 4. This demonstration shows P-T profiles for a model with ${\overline{T}}_{\mathrm{term}}=1900$ K, ΔT_term = 600 K, ${\rm{\Delta }}{T}_{\mathrm{DN}}=1400$ K, and T_deep = 2500 K. The left panel shows the range of profiles in the terminator plane, while the right panel shows the range of profiles sampled by a ray crossing the day–night transition. We also show a 2D representation for the same planet in Figure 3, illustrating how this temperature field translates into local-altitude grid differences for each column (see Section 3.2.3). Finally, we provide a 3D representation of the same temperature field in Figure 5.

Zoom In Zoom Out Reset image size

We encode the composition of an atmosphere via the volume-mixing ratios (X_q ≡ n_q /n_tot) for a set of chemical species (q = 1, 2, ..., N_species). The total number density is determined, assuming the ideal gas law, solely from the temperature field:

$$\begin{eqnarray}&&{n}_{\mathrm{tot}}(P,\phi ,\theta )=\displaystyle \frac{P}{{k}_{B}T(P,\phi ,\theta )},\end{eqnarray} \tag{ 37 }$$

where k_B is the Boltzmann constant. In general, each chemical species can have its own mixing-ratio field, using the same prescription as Section 3.2.1, by replacing the temperature parameters with log-mixing ratio parameters in Equations (34), (35), and (36) (i.e., ${\overline{\mathrm{log}X}}_{{\rm{q}},\mathrm{term}}$, ${\rm{\Delta }}\mathrm{log}{X}_{{\rm{q}},\mathrm{term}}$, ${\rm{\Delta }}\mathrm{log}{X}_{{\rm{q}},\mathrm{DN}}$, and $\mathrm{log}{X}_{{\rm{q}},\mathrm{deep}}$). If the abundance of a given species does not vary in the terminator plane or across the day–night transition, ${\rm{\Delta }}\mathrm{log}{X}_{{\rm{q}},\mathrm{term}}=0$ or ${\rm{\Delta }}\mathrm{log}{X}_{{\rm{q}},\mathrm{DN}}=0$, respectively. Similarly, isocompositional vertical profiles can be obtained by setting $\mathrm{log}{X}_{{\rm{q}},\mathrm{deep}}=\mathrm{log}{X}_{{{\rm{H}}}_{2}{\rm{O}},\mathrm{high}}(\theta ,\phi )$. Eliminating redundant parameters readily allows a 3D, 2D, or 1D mixing-ratio prescription, allowing the adaptation of model complexity in light of data quality. In our convention, ${\rm{\Delta }}\mathrm{log}{X}_{{\rm{q}},\mathrm{term}}\gt 0$ implies a higher abundance on the evening terminator, while ${\rm{\Delta }}\mathrm{log}{X}_{{\rm{q}},\mathrm{DN}}\gt 0$ implies a higher dayside abundance.

TRIDENT currently supports over 50 chemical species, covering a panoply of molecules, atoms, and ions expected in exoplanet atmospheres (see Appendix C). The main atmospheric constituents for planets ranging from ultra-hot Jupiters to temperate terrestrial planets are included (see Section 3.3). When initializing a model atmosphere, we assign one gas as the bulk constituent, with its mixing ratio in each layer and column determined by the requirement that ∑_q X_q = 1. For giant planets, we treat H₂ and He as a single component with a fixed abundance ratio of ${X}_{\mathrm{He}}/{X}_{{{\rm{H}}}_{2}}=0.17$.

The radial grid of an atmosphere is readily computed from the temperature and mixing-ratio fields. Under the assumptions of hydrostatic equilibrium and inverse-square gravity, the distance from the planetary center to a layer at pressure P in each column is given by

$$\begin{eqnarray}&&r(P,\phi ,\theta )={\left(\displaystyle \frac{1}{{R}_{{\rm{p}},\,\mathrm{deep}}}+{\int }_{{P}_{\mathrm{deep}}}^{P}\displaystyle \frac{{k}_{B}\,T(P,\phi ,\theta )}{G\,{M}_{{\rm{p}}}\,\mu (P,\phi ,\theta )}\,\displaystyle \frac{{dP}}{P}\right)}^{-1},\end{eqnarray} \tag{ 38 }$$

where G is the gravitational constant, M_p is the planetary mass, and μ(P, ϕ, θ) = ∑_q m_q X_q (P, ϕ, θ) is the mean molecular mass in a given layer and column. For multidimensional models, one no longer has the freedom to chose between a reference pressure and a reference radius (e.g., Heng & Kitzmann 2017; Welbanks & Madhusudhan 2019). Choosing a reference pressure in the upper atmosphere (e.g., 1 mbar) would necessitate specifying a different reference radius for every column, due to the nonspherical atmosphere (see Figures 3 and 5). The converse is not true if one specifies a radius at an altitude where the pressure is the same in each column. We therefore assume in Equation (38) that the deep interior of the atmosphere (here, 10 bar) is homogeneous, such that each column shares a common spherical surface defined by r(P_deep, ϕ, θ) = R_{p, deep}. Once the radial grid for a model is determined, one can readily compute the elements of the path distribution tensor (Equation (26)) required for 3D radiative transfer.

The extinction coefficient due to absorption by the chemical species in an atmosphere is given by

$$\begin{eqnarray}&&{\kappa }_{\lambda ,\mathrm{chem}}(P,\phi ,\theta )=\sum _{q=1}^{{N}_{\mathrm{species}}}{n}_{q}(P,\phi ,\theta )\,{\sigma }_{\lambda ,\mathrm{abs},q}(P,T),\end{eqnarray} \tag{ 40 }$$

where n_q = X_q n_tot is the number density of species q and σ_λ,abs,q is the corresponding absorption cross section (in m²). TRIDENT uses a database of precomputed cross sections shared with the POSEIDON retrieval code (MacDonald & Madhusudhan 2017). We calculated pressure- and temperature-broadened cross sections from various line list databases, ExoMol (Tennyson et al. 2020), HITRAN (Gordon et al. 2022), and VALD3 (Ryabchikova et al. 2015), ranging from 100–3500 K and 10⁻⁶–100 bar. We treat each line as a Voigt profile ¹¹ from the line core up to a wing cutoff (at the minimum of 500 Voigt half-width half maxima or 30 cm⁻¹; see Hedges & Madhusudhan 2016), accounting for H₂ and He broadening where data is available. Our cross sections are calculated at a spectral resolution of Δν = 0.01 cm⁻¹, making them suitable for high-resolution models, and generally range from 0.4–50 μm (we extend down to 0.2 μm for some atoms and ions with strong near-UV opacities). Further details on our cross-section computations can be found in MacDonald (2019), where we further find excellent agreement with EXOCROSS (Yurchenko et al. 2018a). We also include bound-free absorption (presently for H⁻) in κ_λ,chem. We show representative cross sections from our database in Figure 6 (top panels); see Appendix C for a full summary of our line list sources.

Zoom In Zoom Out Reset image size

We model Rayleigh scattering as a pure loss term ¹² with extinction expressed as

$$\begin{eqnarray}&&{\kappa }_{\lambda ,\mathrm{Rayleigh}}(P,\phi ,\theta )=\sum _{q=1}^{{N}_{\mathrm{species}}}{n}_{q}(P,\phi ,\theta )\,{\sigma }_{\lambda ,\mathrm{Rayleigh},q},\end{eqnarray} \tag{ 41 }$$

where the Rayleigh cross section can be written as either

$$\begin{eqnarray}&&{\sigma }_{\lambda ,\mathrm{Rayleigh},q}=\displaystyle \frac{24{\pi }^{3}}{{n}_{\mathrm{ref},q}^{2}\,{\lambda }^{4}}{\left(\displaystyle \frac{{\eta }_{q}^{2}(\lambda )-1}{{\eta }_{q}^{2}(\lambda )+2}\right)}^{2}\,{F}_{\mathrm{King},q}\,(\lambda )\end{eqnarray} \tag{ 42 }$$

or, equivalently,

$$\begin{eqnarray}&&{\sigma }_{\lambda ,\mathrm{Rayleigh},q}=\displaystyle \frac{128}{3}{\pi }^{5}\,{\bar{\alpha }}_{q}^{2}(\lambda )\,{\lambda }^{-4}\,{F}_{\mathrm{King},q}\,(\lambda ).\end{eqnarray} \tag{ 43 }$$

In the expressions above η_q is the refractive index of species q (measured at reference number density n_ref,q ), F_King,q is the King correction factor (accounting for nonspherical charge distributions), and ${\bar{\alpha }}_{q}$ is the polarizability. We employ various fitting functions for the refractive indices and King correction factors (e.g., Hohm 1994; Sneep & Ubachs 2005). Where such data are not available, we use static polarizabilities from the CRC handbook (Haynes 2014) and assume F_King,q = 1.

Continuum absorption can also arise from interactions between pairs of chemical species. The most common pair-process opacity in exoplanet atmospheres is collision-induced absorption (CIA), though free–free absorption (e.g., from H⁻) is also a pair process (e.g., see Sharp & Burrows 2007). Unlike single-species absorption or scattering, pair processes are proportional to ${n}_{\mathrm{tot}}^{2}$ with the following functional form:

$$\begin{eqnarray}&&{\kappa }_{\lambda ,\mathrm{pair}}(P,\phi ,\theta )=\sum _{\mathrm{pairs}}{n}_{{q}_{1}}(P,\phi ,\theta )\,{n}_{{q}_{2}}(P,\phi ,\theta )\,{\sigma }_{\lambda ,\mathrm{pair},q}\,(T),\end{eqnarray} \tag{ 44 }$$

where ${n}_{{q}_{1}}$ and ${n}_{{q}_{2}}$ are the number densities of the two species constituting the pair (e.g., for H₂-He CIA they are ${n}_{{{\rm{H}}}_{2}}$ and n_He) and σ_λ,pair,q is the binary absorption cross section (in m⁵) for the interaction. We use temperature-dependent CIA data from the HITRAN CIA database (Karman et al. 2019) and free–free absorption data for H⁻ from John (1988). Where data is not available for a given temperature, we assume the cross section is equal to the nearest temperature with tabulated data. We further assume the cross section is zero at wavelengths outside the data range. We show the range of pair-process cross sections in our opacity database, relevant for both gas giant and terrestrial planets, in Figure 6 (bottom panels).

Aerosols contribute to atmospheric extinction via both absorption and scattering. In the limit where forward scattering is neglected ¹³ we can generally write the aerosol extinction as

$$\begin{eqnarray}&&{\kappa }_{\lambda ,\mathrm{aerosol}}(P,\phi ,\theta )=\sum _{q=1}^{{N}_{\mathrm{aerosol}}}{n}_{q}(P,\phi ,\theta )\,{\overline{\sigma }}_{\lambda ,\mathrm{ext},q},\end{eqnarray} \tag{ 45 }$$

where the mean extinction cross section is a weighted average over a (normalized) particle size distribution:

$$\begin{eqnarray}&&{\overline{\sigma }}_{\lambda ,\mathrm{ext},q}={\int }_{0}^{\infty }{\sigma }_{\lambda ,\mathrm{ext},q}(a)\,f(a)\,{da}.\end{eqnarray} \tag{ 46 }$$

The aerosol cross section for a specific particle size, σ_λ,ext,q (a), can be analytically determined via treatments such as Mie theory. Many studies have already considered the spectral impact of Mie scattering (e.g., Wakeford & Sing 2015; Pinhas & Madhusudhan 2017; Benneke et al. 2019), so we focus here on how spatial variations in cloud opacity impact transmission spectra.

A day–night temperature gradient strengthens the amplitude of transmission spectra absorption features (Figure 10, top panels). The increased band strengths arise from longer ray paths through the dayside (see Figure 3), which effectively increases the equivalent scale height of the atmosphere. This agrees with the analytic derivation by Caldas et al. (2019), who showed that day–night temperature gradients induce a higher slant optical depth than a 1D atmosphere with the same terminator temperature. For a given temperature difference, increasing the opening angle (β) counteracts the increased absorption from the dayside. This effect stems from the finite angular range a ray can sample about the terminator plane. Since the maximum zenith angle a ray samples upon entering the dayside is ${\theta }_{\max }=2{\cos }^{-1}(b/{R}_{{\rm{p}},\,\mathrm{top}})$, we estimate that ${\theta }_{\max }\approx 67^\circ$ for our HD 209458b model (assuming the deepest ray samples down to 0.1 bar). Therefore, as $\beta \to {\theta }_{\max }$ the temperature gradient becomes sufficiently shallow along the line of sight that the 1D limit is recovered (Figure 10, top right).

Zoom In Zoom Out Reset image size

The detectability of day–night temperature gradients from transmission spectra has mixed prospects. Though the residuals from even moderate temperature gradients (${\rm{\Delta }}{T}_{\mathrm{DN}}\gtrsim 600$ K) can exceed 100 ppm, the detection of a 2D day–night temperature difference relies on the inability of a 1D model to compensate for the residuals. Since the residuals for a day–night temperature gradient manifest as a scale height increase, a 1D model can simply increase the planet temperature to compensate. Indeed, Caldas et al. (2019) showed that a 1D retrieval of a spectrum with such a day–night gradient infers a biased temperature warmer than the true terminator temperature. We show in Section 5.2 that this holds only when one assumes uniform abundances, with the presence of a composition gradient substantially improving the detection prospects for day–night temperature gradients.

A morning–evening temperature gradient can strengthen or weaken absorption features (Figure 10, bottom panels). This effect is distinct from day–night temperature gradients, albeit yielding residuals ∼10× less pronounced. The small morning–evening residuals arise from rays sampling both the morning and evening terminators (via the ϕ integral in Equation (9)). The greater effective area of the warmer evening terminator therefore largely cancels with the lower effective area of the colder morning terminator. This agrees with an analytic result from MacDonald et al. (2020), where we showed that morning–evening temperature gradients negligibly alter 2D transmission spectra from their 1D counterpart. However, the superposition of the morning and evening terminators is not exactly identical to a 1D spectrum with the average terminator temperature. The residuals are caused by the different cross-section temperature dependencies of each chemical species in the atmosphere. For example, the residuals in H₂O bands have a different magnitude (and sign) than the 4.3 μm CO₂ feature (see Figure 10). Increasing the opening angle (α) reduces the residuals, since a greater fraction of the azimuthal integral covers sectors with a similar temperature to the average terminator temperature.

A morning–evening temperature gradient may prove detectable in high-precision (∼20 ppm) transmission spectra. A morning–evening temperature gradient distorts the shape of absorption bands (Figure 10). For example, the wings of the 3 μm H₂O band strengthen more than the corresponding absorption peak. Since this residual structure is distinct for different molecular bands and for different molecules, a multidimensional retrieval of a sufficiently precise data set, covering a wide wavelength range, may provide a better fit than a 1D model. While the magnitude of morning–evening temperature gradient residuals is comparable to the ∼20 ppm JWST noise floor (Greene et al. 2016), we show in the next section that the addition of a composition gradient significantly amplifies the residuals.

Pure composition gradients.Day–night composition gradients produce three key differences compared to our 1D reference model (Figure 11, top left):

• 1.  
  A transit depth increase in bands of the species exhibiting the composition gradient (here, H₂O).
• 2.  
  Weaker H₂O bands (e.g., 1 μm) are strengthened more than stronger bands (e.g., 3 μm).
• 3.  
  Large composition gradients lower the short-wavelength continuum transit depth.

Zoom In Zoom Out Reset image size

To understand the first effect, consider a model with 100 × less H₂O on the dayside than on the nightside (${\rm{\Delta }}\mathrm{log}{X}_{{{\rm{H}}}_{2}{\rm{O}},\mathrm{DN}}=-2$). Even though a stellar ray traverses the same distance through the dayside and nightside (for pure composition gradients), the nightside dominates the slant optical depth. Consequently, all H₂O absorption bands are stronger than the 1D reference model with a uniform H₂O abundance equal to that at the terminator plane (10× less than the nightside).

The second effect, unequal band strengthening, comes from geometry: wavelengths with low H₂O opacity allow the transmission of rays with deeper impact parameters. Since low-impact-parameter rays traverse a greater distance in the dayside/nightside before encountering the transition region (see Figure 3, right panel), the high-H₂O hemisphere has a more pronounced effect than for rays only sampling the upper atmosphere.

The third effect is an indirect consequence of a high-H₂O abundance gradient. Despite these optical wavelengths having negligible H₂O opacity, the H₂O abundance dichotomy slightly alters the mean molecular weight of the dayside and nightside. The smaller scale height in the high-H₂O hemisphere then weakens the alkali wings and H₂ Rayleigh opacities.

Finally, we note that models with more H₂O on the dayside (positive gradient) produce identical spectra to those with the same H₂O enhancement on the nightside (negative gradient). This symmetry holds only for pure composition gradients, since the dayside and nightside—with the same physical size and temperature in this case—can be freely interchanged without altering the slant optical depths through the atmosphere.

Composition & temperature gradients.Combined, day–night composition and temperature gradients change the relative amplitudes of absorption bands (Figure 11, top right). Compared to pure composition gradients, a temperature gradient breaks the symmetry between models with positive or negative H₂O gradients. Models with a higher dayside H₂O abundance have significantly amplified H₂O bands due to the extended ray paths through the hot dayside. Conversely, models with a higher nightside H₂O abundance have much weaker H₂O bands due to the shortened ray paths through the cool nightside. In concert, the effective scale height increase (from temperature gradients) and band strengthening (from composition gradients) alters the relative band strengths of the species with a day–night abundance gradient.

One can detect day–night composition and temperature gradients via the distortion of band strengths compared to 1D models. Atmospheres with lower dayside abundances can have lower band amplitudes at wavelengths with strong absorption—alongside higher-band amplitudes at wavelengths with weak absorption—than a 1D model with the same abundance as the terminator plane (e.g., compare the 3 μm and 1 μm H₂O bands in Figure 11, top right). Alternatively, atmospheres with a higher dayside abundance display strongly amplified absorption features. The residuals from day–night composition and temperature gradients are often greater than 100 ppm, while exhibiting a structure different from the scale height enhancement seen for pure temperature gradients. Consequently, 1D retrievals attempting to fit a spectrum with day–night composition gradients experience significant biases in retrieved abundances, temperatures, and planetary radii (Pluriel et al. 2022). However, the same phenomenon causing these biases may offer the solution to reliably identify day–night gradients with multidimensional retrievals. Our results suggest that transmission spectra data sets covering multiple absorption bands of species expected to exhibit composition gradients is the key to detect day–night gradients.

Pure composition gradients. Morning–evening composition gradients produce three differences compared to our 1D reference model (Figure 11, bottom left):

• 1.  
  The bands of the species exhibiting the composition gradient (here, H₂O) weaken at wavelengths where it dominates the opacity.
• 2.  
  Conversely, H₂O bands strengthen at wavelengths where another overlapping species dominates the opacity (e.g., in the K line wings).
• 3.  
  Composition gradients lower the short-wavelength continuum transit depth.

The first effect arises from the superposition of the two terminators (and transition region) not being equivalent to a 1D atmosphere with the average terminator H₂O abundance. Specifically, the weighted sum of regions with lower H₂O opacity and regions with higher H₂O opacity produces a weaker H₂O band than our 1D reference atmosphere (with the same log-mixing ratio of H₂O as the average of the morning and evening terminators). We analytically demonstrated this effect in MacDonald et al. (2020), where we showed that the equivalent 1D temperature of a 2D atmosphere with a morning–evening composition gradient is colder than the average terminator temperature. We note that our derivation relied on the assumption that only the species exhibiting the composition gradient dominates the opacity at a given wavelength. Such an assumption holds for strong H₂O bands, hence explaining the negative residuals in most of the infrared seen in Figure 11 (e.g., for H₂O at 3 μm). Consequently, a 1D retrieval would be biased by retrieving a lower temperature than the true terminator average temperature (MacDonald et al. 2020).

The second effect occurs at wavelengths where absorption features from multiple species overlap. Such an opacity overlap amplifies weak H₂O features relative to the 1D case, resulting in positive instead of negative residuals. We see this effect most clearly in the K line wings (near 0.7 μm and 0.9 μm) and the 4.3 μm CO₂ band. Effects (i) and (ii) have opposing influences on the absorption bands of the species displaying the composition gradient, underscoring the importance of observing both weak and strong molecular bands to detect morning–evening composition gradients.

The third effect is caused by mean molecular weight differences (as with day–night composition gradients). For morning–evening composition gradients, the terminator with the higher mean molecular weight (higher H₂O abundance) has a lower scale height and hence presents a decreased effective area to incoming rays. The terminator with the lower H₂O abundance has a relatively unchanged mean molecular weight (since it is dominated by H₂ and He), so presents a similar effective area. The net effect of the terminator superposition is thus a lower continuum transit depth.

Finally, we see symmetry between models with more H₂O on one terminator and the mirror model with the H₂O gradient reversed. This results from the morning and evening terminators being freely interchangeable for pure composition gradients, since stellar rays sample the full terminator annulus during a transit.

Composition & temperature gradients.Combined, morning–evening composition and temperature gradients modulate the strength of absorption bands and alter the peak-to-wing shape (Figure 11, bottom right). As with day–night combined gradients, the symmetry between positive and negative H₂O gradients breaks with the addition of an interterminator temperature gradient. The first-order effect of a positive H₂O gradient is stronger H₂O bands, arising from the warmer evening terminator presenting both a greater effective area and a higher H₂O opacity. Conversely, a negative H₂O gradient significantly weakens H₂O bands due to the colder morning terminator (with more H₂O) presenting a lower effective area than the warmer evening terminator (with less H₂O). The residuals for morning–evening composition and temperature gradients have a unique signature: they are anticorrelated with the shape of H₂O bands. For positive H₂O gradients, the wings of absorption features strengthen more than the band peak. For negative H₂O gradients, the wings weaken less than the band peak. Consequently, the wing-to-peak slope of absorption bands becomes shallower as the magnitude of morning–evening composition gradients increases.

One can detect morning–evening composition and temperature gradients via the sensitivity of peak-to-wing band shapes to composition gradients. We stress that morning–evening composition gradients are far more detectable than pure temperature gradients. Comparing Figures 10 and 11, we see that the residuals for combined morning–evening composition and temperature gradients often exceed 100 ppm or are 3–5× more prominent than for pure temperature gradients. Given the magnitude and complex structure of the residuals, multidimensional retrievals are required to correctly interpret transmission spectra shaped by morning–evening gradients (MacDonald et al. 2020; Espinoza & Jones 2021). Our results indicate two requirements to detect morning–evening gradients: (i) data with sufficient spectral resolution to resolve band wing shapes (R ∼ 100 is sufficient); and (ii) data with a wide wavelength baseline covering both strong and weak absorption bands for key chemical species. We caution that patchy clouds can also alter band wing shapes in the infrared (see Line & Parmentier 2016; MacDonald & Madhusudhan 2017; and Section 5.3), so uniquely identifying morning–evening composition gradients requires both visible wavelength and infrared observations. Finally, the distinct ways in which day–night and morning–evening gradients influence multidimensional transmission spectra strongly implies that a 3D retrieval code can simultaneously detect gradients around and across exoplanet terminators.

The cloud-top pressure of an Iceberg cloud demarcates the atmospheric layers where an additional continuum opacity source is present. Low-pressure clouds elevate the impact parameter for which the slant optical depth equals unity, resulting in a higher transit depth at all wavelengths (Figure 12, top left). Iceberg clouds do not act as a hard surface, unlike the classical 1D cloud deck, since the nonuniform terminator cloud coverage allows rays with certain azimuthal angles to sample clear atmospheric regions (see Figure 7).

The cloud extinction controls the overall cloud transparency. Our reference extinction, κ_cld = 10⁻⁵ m⁻¹, corresponds to an optically thick cloud at all wavelengths. Partially transparent clouds exist over a narrow range from κ_cld ∼ 10⁻⁷–10⁻⁸ m⁻¹, with the cloud presenting negligible opacity for κ_cld ≤10⁻⁹ m⁻¹ (Figure 12, top right). The long slant paths in transmission geometry cause even moderate cloud extinction to rapidly behave as an opaque cloud (Fortney 2005), explaining the limited range of κ_cld affording partial transparency. The cloud extinction strongly affects spectral regions probing deep layers in the atmosphere, such as the alkali wings and Rayleigh slope in the optical, with little impact on strong absorption bands forming mostly above the cloud-top pressure (e.g., the 3 μm H₂O band).

The cloud zenith angle dependence reveals how far from the terminator plane rays are sensitive to cloud opacity. Clouds commencing in the dayside (θ_0,cld < 0°) are encountered by rays at an earlier stage, resulting in higher transit depths and shallower spectral features (Figure 12, middle left). We see no change in the spectrum between a cloud with θ_0,cld = −10° and θ_0,cld = −20°. Conversely, a cloud extending 10° into the nightside (θ_0,cld = 10°) produces almost the same spectrum as a clear atmosphere. Therefore, our 3D atmosphere has an effective cloud horizon constrained to $| {\theta }_{0,{\rm{cld}}}=\lesssim {10}^{\circ }$. The sensitivity to slightly larger nightside zenith angles is due to the smaller nightside scale height (and hence a shorter slant path traversed per unit zenith angle). From Figure 12, we see that the cloud zenith angle has a similar influence on 3D transmission spectra to the cloud extinction. Therefore, we expect these two parameters to be somewhat degenerate in multidimensional retrievals.

Sensitivity to the cloud azimuthal starting angle is a unique feature of multidimensional atmospheres. A 1D background atmosphere with a 2D patchy cloud (e.g., Line & Parmentier 2016) is insensitive to ϕ_0,cld, since one can azimuthally rotate the cloud's starting location and obtain an identical transmission spectrum. However, a nonuniform background temperature or composition breaks this symmetry. Consequently, multidimensional transmission spectra are weakly sensitive to ϕ_0,cld (Figure 12, middle right). We see that the visible wavelength continuum changes as the cloud rotates, due to the cloud obscuring sectors with different temperatures as ϕ_0,cld varies. Further, absorption bands forming below the cloud top are modulated for species exhibiting a morning–evening gradient. For example, the 1 μm H₂O band is strong when ϕ_0,cld = 40°, since 40°/180° = 22% of the morning half-annulus (containing a higher H₂O abundance than the evening half-annulus) is unobscured by clouds. Conversely, the same H₂O band is weaker when ϕ_0,cld = −40°, since only (40 − 0.1∗180)°/180° = 12% of the morning half-annulus is unobscured. Despite the lower sensitivity to ϕ_0,cld than for the other Iceberg cloud parameters, our results demonstrate a key point: morning–evening temperature and composition gradients break azimuthal symmetry, allowing one to infer the region of the terminator annulus containing 3D clouds.

The terminator cloud coverage controls the fraction of terminator containing cloud opacity. Higher cloud fractions increase the transit depth, decrease the amplitudes of spectral features, and alter absorption band wing shapes (Figure 12, bottom). When f_cld = 1, one recovers the limit of a uniform cloud deck which emulates a hard surface. When f_cld = 0, one obtains a clear atmosphere. We also note that the f_cld and ϕ_0,cld parameters interface in several ways, such as (i) f_cld = 1 means ϕ_0,cld has no effect, and (ii) f_cld = 0.5 yields the same spectrum for ϕ_0,cld and −ϕ_0,cld (due to north–south symmetry of the background atmosphere). Given the strong impact of patchy clouds on transmission spectra—and the unique way f_cld alters spectral features—evidence of patchy clouds has already emerged from Hubble spectra of hot Jupiters (MacDonald & Madhusudhan 2017; Pinhas et al. 2019; Barstow 2020). Our analysis builds on these inferences by demonstrating that multidimensional transmission spectra encode not only the terminator cloud fraction, but also the spatial location of a 3D cloud.

Ingress and egress spectroscopy is one key application of TRIDENT. When a planet partially transits its host star, the fraction of the planet and its atmosphere occulting the stellar disk changes with time. Consequently, transmission spectra are time dependent during partial transits. All transiting planets undergo partial transits during ingress and egress, producing distinct time-averaged spectra for the ingress and egress if the morning and evening terminators differ in their temperature, composition, or aerosol properties (e.g., Fortney et al. 2010; Kempton et al. 2017). We show an ingress and egress spectrum for a typical hot Jupiter in Figure 13 (top panel), computed with TRIDENT using Equations (2) and (A23). With sufficient time resolution and data precision, one could in principle extract multiple spectra as a function of time from ingress/egress observations. Fitting such time-resolved transmission spectra would open the field of transmission mapping, analogous to secondary eclipse mapping (e.g., de Wit et al. 2012; Majeau et al. 2012; Challener & Rauscher 2022). Transmission mapping is a rich area for future theoretical and observational studies.

Zoom In Zoom Out Reset image size

Grazing transit spectroscopy is another application of TRIDENT. A planet with a grazing transit (satisfying b_p R_* > R_* − R_p) has an ingress/egress lasting the entire transit. To date, over 20 confirmed exoplanets have grazing transits (Davis et al. 2020). We show an example grazing transit spectrum for the hot Jupiter TOI-1130c (Huang et al. 2020), assuming a solar-composition atmosphere at 650 K, in Figure 13 (middle panel). For grazing transits, TRIDENT computes spectra at a series of time steps, with the time-averaged spectrum shown in Figure 13. An extreme grazing transit is WD 1856 + 534b (Vanderburg et al. 2020), a planet transiting a white dwarf with R_p/R_* ≈ 8, for which we show a model, assuming a Jupiter-like composition, in Figure 13 (bottom panel). Recently, we demonstrated the principles of TRIDENT's grazing transit modeling in an analysis of Gemini/GMOS observations of WD 1856 + 534b (Xu et al. 2021). With the success of the Transiting Exoplanet Survey Satellite (TESS), the number of planets with grazing transits is steadily increasing. Since TESS planets orbit nearby bright stars, TRIDENT offers a new capability to model grazing transit spectra for the interpretation of JWST observations.

TRIDENT can conduct line-by-line radiative transfer, enabling high-spectral-resolution model applications. We demonstrated this capability by showing a line-by-line transmission spectrum of a hot Jupiter in Figure 8. Since our cross-section database (Appendix C) is precomputed at R ∼ 10⁶ (see Section 3.3), TRIDENT can therefore generate template spectra for high-resolution cross-correlation analyses (e.g., Snellen et al. 2010; Hoeijmakers et al. 2018; Sedaghati et al. 2021). Future improvements to TRIDENT can focus on physical effects that manifest in high-resolution transmission spectra, including planetary rotation and winds (e.g., Brogi et al. 2016; Flowers et al. 2019; Seidel et al. 2020). Another promising avenue is integrating TRIDENT within a high-resolution retrieval framework (e.g., Brogi & Line 2019; Gibson et al. 2020; Pelletier et al. 2021).

Transmission spectra of terrestrial exoplanets is one of the most significant advances enabled by the JWST. While we have focused in this study on giant planet spectra, TRIDENT also contains a growing database of cross sections applicable for radiative transfer through temperate terrestrial exoplanet atmospheres (see Table 1). Our cross-section database has already been used in the literature to simulate transmission spectra of an Earth-like planet transiting a white dwarf (Kaltenegger 2020) and TRAPPIST-1e (Lin et al. 2021). Future extensions to TRIDENT can add a refractive ray-tracing algorithm, which can be important to accurately simulate transmission spectra of terrestrial atmospheres (e.g., Bétrémieux & Kaltenegger 2014; Misra et al. 2014). The inclusion of refraction requires one to numerically compute the path distribution tensor (Robinson 2017), but is otherwise covered by our theoretical framework.