Titan in Transit: Ultraviolet Stellar Occultation Observations Reveal a Complex Atmospheric Structure

Transit spectroscopy is a key tool for exoplanet atmospheric characterization. However, transit spectrum observations can be limited by aerosol extinction when gas opacities are weak. The ultraviolet wavelength range contains a variety of strong molecular and atomic features, potentially enabling gas species detection even when atmospheric hazes are present. To understand the interplay between aerosol extinction and ultraviolet molecular opacities, we investigate transmission through the atmosphere of Saturn’s moon Titan during an occultation observed with the Ultraviolet Imaging Spectrometer (UVIS) on board NASA’s Cassini orbiter. We analyze the derived ultraviolet transit spectrum of Titan using exoplanet-relevant atmospheric retrieval models that both include and exclude treatments for hazes. Our retrieved gas column densities are consistent with previous studies analyzing UVIS occultation data. Despite the apparent haze impact on the underlying occultation data, our treatments fail to correctly characterize the haze in fits derived from simulated transit observations. This suggests that oversimplified haze parameterizations can hinder detection of atmospheric hazes in transit. Our work indicates that continued characterization of exoplanets in the ultraviolet wavelength regime can provide novel atmospheric constraints even if transit spectra are dominated by haze extinction at longer wavelengths.


Introduction
Understanding the atmospheres of exoplanets provides essential insight into the formation, evolution, and potential habitability of these systems (Seager & Deming 2010;Madhusudhan 2019). Over the last two decades, transit spectroscopy (Seager & Sasselov 2000;Brown 2001;Hubbard et al. 2001) has emerged as the leading technique for characterizing exoplanet atmospheres. Here, atmospheric opacity sources can lead to small variations in the wavelength-dependent dimming of a stellar host during an exoplanet transit event. Despite the subtle nature of this effect, spectroscopic transit observations have yielded detections of atmospheric species in a diversity of exoplanet atmospheres (Charbonneau et al. 2002;Tinetti et al. 2007;Swain et al. 2008;Stevenson et al. 2010;Fraine et al. 2014;Line et al. 2014;Sing et al. 2016;Benneke et al. 2019;Tsiaras et al. 2019).
Of course, not all spectroscopic transit observations have revealed atmospheric features. Especially for lower-mass exoplanets, observations have sometimes revealed flat, featureless transit spectra Kreidberg et al. 2014;de Wit et al. 2016), at least to within measurement uncertainties. Here, the presence of high-altitude aerosols is often used to explain such flat transit spectra. As the transit geometry implies that transit spectra probe long pathlengths along the limb of an exoplanet, even hazes or clouds with small vertical optical depths can appear opaque (Fortney 2005).
In theory, observing transit spectra in wavelength ranges where molecular or atomic opacities are relatively large will probe the upper reaches of exoplanet atmospheres, thereby potentially avoiding the obscuring effects of hazes and providing stronger detections of atmospheric species. This, for example, leads to a key strength of NASA's upcoming James Webb Space Telescope (Gardner et al. 2006), whose spectral coverage overlaps strong molecular rotation-vibration bands in the near-and mid-infrared with relevance to exoplanet transit spectra (Deming et al. 2009;Beichman et al. 2014; Barstow & Irwin 2016;Greene et al. 2016;Batalha & Line 2017). Another key wavelength regime with the potential for strong molecular and atomic opacities-the ultraviolet-was originally suggested as a range with likely high utility (Hubbard et al. 2001) and has recently been exploited to study the atmospheres of WASP-121 b (Sing et al. 2019) and HAT-P-41 b (Wakeford et al. 2020), both building on earlier efforts in the ultraviolet for HD 189733 b by Sing et al. (2011). In these works, Sing et al. (2019) observe strong Fe II and Mg II features for WASP-121 b, while both Wakeford et al. (2020) and Sing et al. (2011) observe sloped transit spectra in the ultraviolet that are consistent with hazes for HAT-P-41 b and HD 189733 b, respectively. In new modeling work, Lothringer et al. (2020) detailed how strong opacities at ultraviolet wavelengths due to metals and metal-bearing species could help probe rainout chemistry in exoplanet atmospheres.
To further explore the interplay between aerosol extinction and absorption due to atmospheric gas species in ultraviolet transit observations, we turn to Titan. For Titan, solar photons, solar wind, galactic cosmic rays, and magnetospheric charged particles drive pervasive atmospheric chemistry resulting in multitudes of higher-order hydrocarbons and the carbon nitrogen aggregate tholins (Podolak et al. 1979;Yung et al. 1984;Toublanc et al. 1995;Lavvas et al. 2008;Vuitton et al. 2009;Lavvas et al. 2011;Carrasco et al. 2018). Additionally, Titan has a seasonally dependent detached haze layer located between 300 and 500 km altitude (or at pressures lower than 10 −5 bar; Lavvas et al. 2009;West et al. 2011West et al. , 2018 that is potentially analogous to the hazes responsible for some featureless exoplanet transit spectra. However, it should be noted that due to drastically different thermal conditions, aerosol haze composition and particle sizes for hot exoplanets may differ from those seen in Titan's atmosphere (Lavvas & Koskinen 2017;Lavvas et al. 2019;Moran et al. 2020). For a recent review of Titan's atmosphere and climate, see Hörst (2017).
Here, we use Titan atmospheric stellar occultation observations from the Ultraviolet Imaging Spectrometer (McClintock et al. 1993;Esposito et al. 2004) aboard NASA's Cassini spacecraft to effectively study a hazy world in transit. These occultation data are converted to exoplanet-like transit spectra following techniques developed by Robinson et al. (2014) and Dalba et al. (2015). Critically, these occultation observations have already been used to derive key atmospheric properties for Titan, including number density profiles for various trace hydrocarbons (Koskinen et al. 2011), thereby helping to confirm aspects of our transit spectral analysis. Moreover, it has been suggested that Titan's hydrocarbon rich atmosphere may be representative of a fairly common class of exoplanets, reinforcing the necessity of understanding the interplay of molecular opacities and haze extinction (Lunine 2010).
Below, we begin by describing our adopted occultation data set and technique for converting this to an ultraviolet transit spectrum of Titan. We then present the details of an atmospheric retrieval model designed to interpret our derived transit spectrum. Following our retrieval analyses, we discuss the impact of Titan's detached haze layer on transit spectra. Finally, we conclude by interpreting our results with respect to the current state of exoplanet observations.

Methods
The following subsections describe our approach to data reduction, modeling, and analysis. First, we briefly describe the underlying occultation data set and how this was transformed into an effective transit spectrum for Titan. Next, we present a simple forward model that we use to fit our ultraviolet transit spectrum of Titan. Finally, we describe our Bayesian approach to atmospheric characterization using our simulated transit spectrum and forward model.

Data Reduction and Transformation
Occultation data were acquired with the far-ultraviolet channel of the Ultraviolet Imaging Spectrometer (UVIS) on board the NASA Cassini orbiter, as detailed in Koskinen et al. (2011). Spectra from the far-ultraviolet channel span 110-190 nm with a spectral resolution of 0.28 nm. Further details regarding the optical specifications and design of the UVIS instrument are described in Esposito et al. (2004). Data used in our analyses are from Cassini flyby T41 I on 23 February 2008, where occultation observations probed 6°S. Of the 12 independent occultation observations presented in Koskinen et al. (2011), the T41 I data set is the best option for both high signal-to-noise-ratio data as well as a strong signature of two detached haze layers. Transmissivity as a function of wavelength and altitude were calculated based on the ratio of the transmitted stellar spectrum along the instrument line of sight and the unocculted reference stellar spectrum (Koskinen et al. 2011). This inherently requires that all spectral variations are atmospheric in nature. It should be noted that variations in host star flux also affect the light transmitted through the atmosphere, which is discussed in more detail in Section 4. Quoted observational uncertainties are derived based on photon counting (Esposito et al. 2004;Koskinen et al. 2011). The altitude-and wavelength-dependent transmission data from Cassini flyby T41 I are shown in Figure 1. Note the two distinct vertically-isolated opacity sources in the atmosphere near 500 km and 700 km. The lower Figure 1. Altitude-and wavelength-dependent transmission data for Titan's atmosphere from Cassini flyby T41 I. An altitude of 0 km corresponds to Titan's average surface radius (2575 km), and darker colors indicate lower transmission. We note that these data are also shown in Koskinen et al. (2011, their Figure 8a). region corresponds to the detached haze layer within Titan's atmosphere (Rages and Pollack 1983;Porco et al. 2005). The upper region can also be attributed to higher-order hydrocarbons, and the structure may be due to atmospheric propagation of gravity waves (Strobel 2006;Koskinen et al. 2011). We note that the chemistry and dynamics of this atmospheric haze are complex and are coupled to aspects of the deep atmosphere to which our study is not sensitive (West et al. 2018).
Although the process of computing transit spectra from altitude-dependent transmissivities inherently smooths over the smaller-scale effects induced by gravity waves, how gravity waves may manifest in transit analyses merits some discussion. Gravity waves transfer momentum and energy within an atmosphere, which can cause variation in thermal profiles and, consequently, density profiles of some species (Strobel 2006;Koskinen et al. 2011). These perturbations can be substantial, particularly for hot, rotationally-locked exoplanets (Watkins and Cho 2010). Local increases in wind speed due to the deposition of energy via gravity waves may result in Doppler shifted spectral features. Doppler shifted spectral features in exoplanet transmission spectra have been previously reported for wind speeds of order 1-10 km s −1 (Snellen et al. 2010;Kempton & Rauscher 2012).
The transit depth spectrum corresponding to the altitudedependent atmospheric transmission data was calculated following methods described in Robinson et al. (2014). Given the wavelength-dependent transmissivity (t λ,i ) on a grid of impact parameters (b i ; taken as the radial distance of closest approach for a ray), the effective area of the Sun that would be blocked by an atmospheric annulus spanning levels i to i + 1 if Titan were observed in transit is given by ) . The wavelength-dependent transit depth is then computed by summing over annuli, with, where R p is the solid surface radius of the planet or, for worlds without a solid surface, a sufficiently deep reference radius that all transmissivity values have reached zero. Corresponding transit depth uncertainties are calculated following standard Gaussian error propagation techniques and are described further in Section 3.1 (Taylor 1997). Equation (2) can be solved for the wavelength-dependent planetary radius ( l R p, 2 ), and subtracting the solid surface radius of Titan from this yields the so-called effective transit altitude (or height), z eff . Unlike in Robinson et al. (2014), losses due to refraction can be ignored here as the occultation observations probe much greater altitudes (i.e., very low pressures and number densities). In all analyses, UVIS data between 185-190 nm are omitted due to higher-order hydrocarbon features for which we do not have sufficient absorption cross section data. Species suspected to be responsible for these unidentified features include toluene and xylene (Koskinen et al. 2011). Figure 2 shows the transit spectrum that results from application of the previously described techniques, depicted here as effective transit altitude to help indicate where (vertically) in the atmosphere the transit spectrum probes.
Despite the detached haze layers present over this altitude regime, the modeled spectrum is rich with molecular features. Several notable features include the broad methane absorption centered near 130 nm, acetylene and diacetylene features between 140 nm and 150 nm, and the ethylene features around 170 nm.
As a check of our UVIS-derived transit spectrum, we generated a transit depth spectrum from the species slant path column number density profiles inferred from our adopted occultation data in Koskinen et al. (2011). Here, the inferred slant column number densities were combined with known gas opacities to compute the transmissivities for Equation (1). Figure 3 shows the transit depth spectra derived from the UVIS data and from the Koskinen et al. (2011) species profiles. The dark and light blue spectra correspond to transit depth spectra derived with and without the tholin haze profile, which helps isolate and quantify the impact of haze on the UVIS-derived transits spectrum. Notably, haze impacts the transit depth spectrum at wavelengths where the altitudes probed are at or below roughly 800 km. This observed influence of haze on the transit spectrum will be used in later sections to help explore the efficacy of cloud parameterizations often used in the exoplanet atmospheric retrieval literature.

Forward Model
We seek to define a forward model that (1) will enable retrievals of key atmospheric parameters when applied to our derived transit depth spectrum and (2) is analogous in complexity to similar models for exoplanets. Following Benneke & Seager (2012) and Robinson et al. (2014), if the number density of an extincting species is distributed exponentially with scale height H in an atmosphere, and if the extinction cross section for this species is pressureindependent, then the optical depth integrated along a slant path for an impact parameter b and for a single species is given by where a subscript "j" indicates the species, N 0,j is the vertical column number density above a sufficiently deep reference radius R 0 , σ λ,j is the absorption cross section for the species, and K n is a modified Bessel function of the second kind. Assuming all species are distributed with the same scale height, summing over all species yields the total slant optical depth, Given the total slant optical depth, and assuming that extinction optical depth is equivalent to absorption optical depth (see Robinson et al. 2017), the transit depth is obtained by integrating over impact parameters, Finally, for numerical implementation it is convenient to define a dimensionless parameter, β = b/H, so that where a grid of β values can be straightforwardly designed to ensure vertical resolution that is (at least) finer than H. It is important to note that our forward model is consistent with common assumptions made within the exoplanet retrieval literature of isothermal atmospheres and constant mixing ratios of each species (e.g., Heng & Kitzmann 2017). Following Koskinen et al. (2011), our gaseous opacity sources emphasize hydrocarbon and nitrile species of substantial number densities in Titan's atmosphere and that have non-negligible absorption cross sections at deep ultraviolet wavelengths. Species used in this study are shown in Table 1 along with wavelength coverage, measurement temperature, and a reference paper for our adopted opacity data. Opacity data were selected on the basis of wavelength range and temperature relevance. Absorption cross sections for the eight Figure 3. The UVIS-derived Titan transit spectrum (black), and models that do (dark blue) and do not (light blue) include the spectral impacts of a high-altitude haze detected in Koskinen et al. (2011). Models assume altitude-dependent chemical compositions from fits in Koskinen et al. (2011). 115-168 173 Ferradaz et al. (2009) species used in this study are shown in Figure 4. To treat opaque atmospheric haze layers within our forward model, we defined an altitude below which the optical depth is inflated to simulate complete extinction. We explore two simple haze treatments in our forward model. Though we exclude a treatment of the haze chemistry or diffusion within the atmosphere, which is known to seasonally alter Titan's ultraviolet and visible spectrum, this approach will determine whether a simple haze structure can be extracted from the ultraviolet and visible transit spectrum of a Titan-like exoplanet without a priori information of the haze chemistry. Our two-parameter haze model specifies the lower boundary of a haze layer z h , measured above R 0 as well as a gray haze vertical optical depth τ h , which is distributed uniformly over one scale height. Slant optical depths are computed following the path distribution approach presented in Robinson (2017). A one-parameter haze model simply assumes the atmosphere is opaque for impact parameters that probe below z h (see, e.g., Bétrémieux & Swain 2017;Kempton et al. 2017), which serves to block all light that would probe atmospheric layers deeper than the haze altitude and, as we show later, is instructive for interpreting results from our two-parameter haze treatment. These simplified haze treatments are typical of many exoplanet atmospheric retrievals (Barstow 2020, 2021).

Markov Chain Monte Carlo Implementation
We retrieve atmospheric parameters by fitting our forward model to our UVIS-derived transit spectrum using a standard Markov chain Monte Carlo (MCMC) approach to Bayesian inference. Specifically, we adopt the widely-used MCMC tool emcee, developed by Foreman- Mackey et al. (2013). Bayesian statistics describe the posterior conditional probability distribution of a set of parameters using the prior probabilities of those same parameters and their likelihood probabilities given the data. The posterior distribution enables inference of physical parameters as well as uncertainties on these parameters.
Briefly, Bayes' Theorem states, where y are the observed data, θ is our set of atmospheric parameters (i.e., our atmospheric state parameters), P(θ) is the prior probability for these parameters, and P(y|θ) is the socalled likelihood function. The value of the likelihood function depends on the data, its associated errors, and the atmospheric state parameters. When evaluating the likelihood function, parameter values that result in poorly fit data are penalized. In our approach, the log of the likelihood function (often called the log likelihood) is given by where y k is the kth spectral data point, σ k is the uncertainty for this data point, and f (θ) is the transit depth (forward) model. We adopt uninformed prior probability distributions that simply set physical limits on the values that each parameter can take. Model free parameters, their descriptions, and corresponding prior probability limits for our analyses are shown in Table 2. Our prior probabilities are largely based on the physical limits of the properties that the parameters represent (e.g., all lengths cannot be smaller than zero). Column number densities are retrieved in log space to ensure full exploration of parameter space over multiple orders of magnitude. We note that the lack of prior knowledge regarding the expected range of exoplanet column number densities would require careful attention to ensure the priors do not unphysically truncate the posterior distributions. In later sections, we use the Bayesian Information Criterion (BIC) to determine if the addition of a haze layer is warranted in our forward model. The BIC has become increasingly common in the astronomy community by providing a robust means for model comparisons (e.g., Littenberg & Cornish 2009;Feng et al. 2016;Sharma 2017), and is defined by where ν is the number of model free parameters, N is the number of spectral data points, and q P y max ( | ) is the maximized log likelihood. The change in the BIC, ΔBIC, between two models is then evaluated to support or reject adding additional parameters (e.g., z h and τ h in our case) to improve the model fits. The model with a lower BIC value is preferred over another. Following Kass and Raftery (1995), ΔBIC values between 0 and 2 suggest that the evidence against the higher BIC model is "not worth more than a bare mention." Values between 2 and 6 suggest positive evidence against the larger BIC model. Values between 6 and 10 suggest strong evidence against the larger BIC model. Values greater than 10 suggest very strong evidence against the larger BIC model. A more detailed description of the BIC and how to compare models can be found in Kass and Raftery (1995).

Results
Clear sky (i.e., haze-free) and hazy models were applied within our retrieval framework to our derived transit spectrum of Titan. As discussed immediately below, retrievals were initially performed assuming only photon counting errors in the observations. Later retrievals-more analogous to an exoplanet case-were performed with inflated error bars. As described in Section 2.2, adopted haze models are simplified and similar to those used in some exoplanet atmospheric retrieval studies. Given that the haze has a notable impact on the transit spectrum (Figure 3), our hazy retrieval models offer an opportunity to test the utility of commonly-applied exoplanet cloud/haze parameterizations.

Error Scaling
We initially applied our clear sky and two-parameter haze retrieval models to the UVIS-derived transit spectrum of Titan with only propagated photon counting errors. Initial occultation data errors are 1σ uncertainties from Koskinen et al. (2011). Using Equations (1) and 2, and following standard Gaussian error propagation techniques described in Taylor (1997), these wavelength-and altitude-dependent transmissivity uncertainties (s l t i, ) were mapped to effective transit altitude uncertainties following: Effective transit altitude uncertainties are then propagated to transit depth uncertainties with where z eff is the effective transit altitude and s z eff is the 1σ effective transit altitude uncertainties. Figure 5 shows fitted spectra with uncertainties for the clear sky case. While the modeled spectra appear qualitatively reasonable, the reduced chi-squared for the clear sky best-fit model was 72.1. Similarly, the reduced chi-squared for the two-parameter haze best-fit model was 72.2. Moreover, the probability this reduced chisquared value would occur by chance is infinitesimal, further indicating an oversimplified model or underestimated data error (Flannery et al. 1992). Here, quantitatively poor fits likely stem from a combination of incomplete opacity data and overlysimplified model assumptions. For the former, and as can be seen in Table 1, the laboratory opacity data are often measured in a temperature regime that is much warmer than Titan's upper atmosphere (which has characteristic temperatures of 150-200 K). Additionally, it is worth noting that both best-fit models (our work and Koskinen et al. 2011) struggle to reproduce the acetylene and ethylene features between 145-155 nm despite having cross section data for both species at Titan-relevant temperatures. This may reflect uncertainty in the cross section measurements, which are estimated to be at 10% or less (Wu et al. 2001. Regarding model assumptions, previous analyses of the UVIS occultation observations reveal that the absorbing gaseous species do not have number density profiles that follow simple exponential decreases with a uniform scale height (Koskinen et al. 2011), as is assumed in our model. A comparison of Figures 3 and 5 demonstrates the error associated with our model assumptions.
To produce better fits (as indicated by the reduced chisquared) and to better mimic noisy exoplanet data, we opted to inflate the error bars on our transit spectrum via a uniform multiplicative scaling factor (as is a relatively common practice; Tremaine et al. 2002;Hogg et al. 2010;Foreman-Mackey et al. 2013;Line et al. 2015). Figure 6 demonstrates how increasing the data uncertainty affects the resulting reduced chi-squared for a best-fit clear sky model. Based on this analysis, we chose a multiplicative scaling factor of 8 to produce reduced chi-squared values reliably close to unity without losing spectral information and potentially overfitting the data. As is discussed later, error scaling does not impact any of the conclusions we draw from subsequent retrieval analyses.

Clear Sky Retrievals
Retrieved atmospheric parameters for the application of our clear sky model to the UVIS Titan transit spectrum with inflated uncertainties are shown in Figure 7. This so-called "corners" plot depicts, along the diagonal, the posterior distributions for all retrieved parameters marginalized over all other parameters. Off-axis plots show two-dimensional posterior distributions where all but two retrieved parameters have been marginalized over. For the one-dimensional marginal distributions, values at the 16th, 50th, and 84th percentile are indicated above each subplot (i.e., the distribution mean and ±1σ), and the reduced chi-squared for the best-fit model was 1.126. All one-dimensional marginal distributions are roughly Gaussian in shape. The reference radius (R 0 ) is anti-correlated with all gas column number densities over a narrow range of parameter space as roughly fixed number densities aloft can be maintained when the reference radius is increased but the column number densities down to this radius are decreased. These anti-correlations with R 0 then lead to correlations between all gas column number densities over narrow ranges of parameter space. Finally, Figure 8 shows the 1σ (68.2%) and 2σ (95.5%) spread in model spectra derived from our retrieval.  . Reduced chi-squared values for the best-fit model from our retrieval framework when increasing the multiplicative error scaling factor on the UVIS-derived transit spectrum. The dotted line represents a reduced chi-squared of unity.

Hazy Retrievals
Marginalized posterior distributions for our two-parameter haze model applied to the UVIS-derived transit spectrum are shown in Figure 9. The resultant spectra are shown in Figure 10. As can be seen in the spectra, fits produced by the two-parameter haze model show little to negligible difference in quality as compared with the clear sky model. More quantitatively, the reduced chi-squared for the best-fit twoparameter haze model is 1.130, indicating that the models are nearly identical in their ability to fit the data.
Correlations between gas column number densities and the reference radius parameter in the two-parameter hazy model analysis are similar to those seen in the clear sky analysis. However, as compared with the clear sky analysis, the twoparameter hazy one-dimensional posterior distributions are generally wider and non-Gaussian (excepting the scale height distribution). The detached haze vertical optical depth (τ h ) is poorly constrained, but is generally found to be optically thick in the horizontal (i.e., slant) direction. A tail in the marginal distribution for the reference radius (R 0 ) to smaller radii correlates with z h to maintain the floor of the transit spectrum near 3100 km, or about 500 km above the solid body radius of Titan, which corresponds to the deepest altitudes probed in our transit spectrum. This, in turn, leads to tails in the marginal distributions for the gas vertical column abundances toward larger values. Lastly, it should be noted that the primary benzene feature (∼180 nm) probes near 500 km in the transit spectrum. For this reason, the distribution is likely impacted by the haze parameterization, which sets the floor of the transit spectrum near 500 km.
Results from our one-parameter haze model can be used to further understand the correlations and wide distributions seen in our two-parameter haze treatment. Marginalized posterior distributions from the application of our one-parameter model are shown in Figure 11 and the resulting spread in model spectra are shown in Figure 12. Here, the reference radius and haze altitude combine to set the floor of the transit spectrum near 500 km above Titan's solid body radius. With the haze layer fixing the floor of the transit spectrum, R 0 can take on a wide range of values and, correspondingly, the column number densities also vary over a wide range of values. These behaviors are also seen in the two-parameter haze model results and are analogous to a degeneracy between the reference radius, pressure, and gas abundances discussed for retrievals applied to transiting exoplanets in Heng & Kitzmann (2017).

Discussion
Titan is a well-studied world in the Solar System, which enables intercomparisons between our retrieval analyses and previous works, as is discussed immediately below. Here, we also explore how the BIC parameter enables model selection among our various haze treatments (including haze-free models), and we discuss implications that stem from our nondetections of Titan's tholin haze in our retrievals. Finally, we describe how our transit spectrum results could apply to future exoplanet observations.

Comparisons with Previous Studies
To compare with previous studies, we convert our vertical column number densities (measured above R 0 ) to slant path column number densities. These two column number densities are related through where N s is the slant path column number density of some species at a tangent altitude of z above the solid body radius. Using Equation (14), we randomly re-sampled our MCMCderived distributions for each species to derive slant column number density distributions at 700 km altitude (selected to best compare to results in Koskinen et al. 2011). With these derived distributions, we directly compare our slant column number densities to previous Cassini/UVIS retrievals. Critically, our retrieved column number densities generally agree with previous studies that used Cassini/UVIS observations, as shown in Table 3. Thus, with sufficient signal-tonoise, transit observations can constrain atmospheric properties with uncertainties comparable to orbital measurements. Despite general agreement between our retrievals, Koskinen et al. (2011), andShemansky et al. (2005), column number densities for ethane (C 2 H 6 ) and benzene (C 6 H 6 ) show more substantial deviation and merit some discussion. Ethane is a known photochemical product within Titan's atmosphere (Lavvas et al. 2008). Unfortunately, there are no apparent distinct ethane features within Titan's transit spectrum. Rather, it provides relatively gray opacity between 120-135 nm, with a slight slope longward. Using photochemical calculations and mixing ratio constraints Koskinen et al. (2011) estimate an upper limit for the column number density of ethane. Our Figure 8. Titan transit depth spectrum with 8×inflated error bars (light gray). Modeled transit spectra from our clear sky retrieval analysis are shown as 1σ (68.2%) and 2σ (95.5%) spreads, dark blue and light blue respectively. Note that the difference between the 1σ and 2σ spreads is most discernible in the 175-185 nm range. retrieved column number density, which is obtained without photochemical modeling, is in agreement to within an order of magnitude of previous UVIS results. Additionally, retrieved column number densities for benzene show some variation. These discrepancies likely stem from simplifying assumptions in our model-namely that atmospheric species are distributed exponentially with scale height. Koskinen et al. (2011) demonstrate that column densities for benzene have additional structure with altitude, especially near 700 km. Generally, all retrieved column number density values across studies are fairly consistent.
Comparisons with previous results do present some issues, however. Potentially most importantly, the results from our transit retrievals can be most sensitive to the abundance of a given species at a range of altitudes that is quite distinct from the 700 km altitude level adopted for comparisons. For example, in the far-ultraviolet (∼130 nm) the transit spectrum probes much greater altitudes than the comparison altitude of 700 km. As methane is responsible for the majority of the opacity in this region, our methane constraint is primarily applicable at altitudes markedly larger than 700 km and our assumption of a constant scale height for number density distributions leads to a discrepancy at 700 km. Figure 13 demonstrates this effect by comparing our idealized slant path column number density profiles with those inferred in Koskinen et al. (2011). Our simplified model profiles generally agree with previous UVIS-derived profiles in altitude regions where our best-fit model is most sensitive to the species number density. Sensitivity was determined based on the contribution function, or ∂z eff /∂N j . Deviations between our retrieved profiles and those of Koskinen et al. (2011) occur where our retrievals have limited or no sensitivity, either due to the extent of altitudes probed by the transit spectrum or due to the strength of the given species' opacity.

Model Selection Considerations
To determine a preference for either our clear sky or hazy models, we calculated the BIC for each model using Equation (11). Recall that differences in the BIC between fitted models can quantitatively indicate if, say, the improvement in a fit through the addition of model parameters is warranted. The ΔBIC value was 6 when comparing the clear sky and one-parameter haze model, with the higher individual BIC corresponding to the hazy model. This evidence suggests that there is no reason to include an additional free parameter to fit for an opaque haze layer. For the clear sky and twoparameter haze model comparison, the ΔBIC value was 14. This indicates strong evidence that supports excluding the twoparameter detached haze layer treatment. These findings remain true for data even without inflated error bars, indicating that error inflation did not result in substantive losses in spectral information.
It is tempting to interpret results from our two-parameter haze treatment in Figure 9 as showing a slight preference for models with a detached haze layer near 500 km above Titan's solid body radius with a slant optical depth near unity (i.e., τ h near 0.1). This would correspond to the haze density maximum near 500 km seen in the altitude-resolved transmission data (Figure 1). Future work might find stronger evidence for the detection of this "detached" haze layer (and, potentially, the haze layer/maximum near 700 km) in our transit spectrum by including wavelength-dependent (i.e., nongray) opacity data for organic hazes and/or exploring haze models that do not assume a uniform distribution of aerosol opacity across a scale height. Regarding the latter point, the transmission data in Figure 1 show a sharper increase in aerosol number density near the base of the haze layer near 500 km (Koskinen et al. 2011). Along these lines, future work might also investigate how atmospheric constraints are degraded as observational uncertainties are artificially inflated to increasingly larger levels than are investigated here.

Relevance to Exoplanet Cloud/Haze Parameterizations
Given the quantifiable impact of haze on the UVIS-derived transit spectrum of Titan (Figure 3), it is striking that none of the hazy retrieval models confidently detected a haze layer (regardless of error scaling). This suggests that common models for parameterizing the distribution of hazes/clouds in exoplanet atmospheric retrieval models could lead to false negative nondetections of aerosols. Furthermore, and as can be seen by comparing Figures 7 and 9, false negative results can also be associated with degradations in the inferred distributions of other atmospheric parameters. In the particular case here, the false negative and degradations likely stem from two distinct over-simplifications. First, hazy retrieval models applied here assume a gray haze opacity. Future work could explore other parameterizations of haze opacity, although we caution against any heavy reliance on Titan's aerosol optical properties as these may not be particularly relevant given conditions in hotter exoplanet atmospheres (Lavvas & Koskinen 2017;Lavvas et al. 2019;Moran et al. 2020). Additional laboratory work addressing the optical properties of relevant haze compositions will be increasingly important given the varied and extreme conditions of the numerous classes of exoplanets (He et al. 2020;Moran et al. 2020).
A second over-simplification stems from our assumption that the haze either truncates the transit spectrum below a fitted altitude (i.e., our one-parameter haze model) or that the fitted haze optical depth is distributed uniformly over one scale height above a fitted altitude (i.e., our two-parameter haze Figure 10. Same as Figure 8, except for retrieved spectra that include a two-parameter haze treatment. model). In reality, results from Koskinen et al. (2011) show that the haze profile is continuous below about 1000 km with regions of local increases in density (e.g., near 500 km and 700 km for the data set adopted here). While a two-parameter haze model is common in the exoplanet retrieval literature (for a review see Barstow 2020Barstow , 2021, Titan's high-altitude haze demonstrates the importance of exploring other haze parameterizations (e.g., an exponential vertical haze distribution, potentially with its own distinct scale height) when applying retrieval analysis to exoplanet observations. The inability of our retrieval models to detect Titan's tholin haze highlights the importance of appropriately balancing simplicity with physicality in treatments of clouds/hazes.

Considerations for Future Exoplanet Observations
While our results demonstrate the rich atmospheric information that can be provided by observations of exoplanet transits at ultraviolet wavelengths, the acquisition of such data comes with additional complications. Perhaps most fundamentally, the overall faintness of Sun-like and cool stars at ultraviolet wavelengths implies that reaching even 10-100 ppm accuracy on transit spectrum observations is difficult. Specifically, shortward of 140 nm, Sun-like stellar spectra exclusively consist of atomic emission lines and minimal continuum severely limiting the SNR (Young et al. 2018). Beyond this, the so-called "transit light source effect" will be most pronounced at short wavelengths (see, e.g., Rackham et al. 2018). Here, Figure 11. Same as Figure 7, except with the addition of a single-parameter haze treatment. occulted heterogeneities on the stellar photosphere-namely spots and faculae-can introduce systematic biases in transit spectra, and the contrast between these heterogeneities and the background photosphere is strong in the ultraviolet regime (Oshagh et al. 2014;Llama & Shkolnik 2015). To consider how our results for Titan might translate to other star-planet systems, we note that the transit features we observe span 6-12 scale heights. For a sub-Neptune with a gravity of 10 m s −2 and an atmospheric temperature of 500 K, this range of scale heights translates to 100-200 km for a water-dominated atmosphere and 1000-2000 km for a H 2 /He-dominated atmosphere. This range of thicknesses (100-2000 km) corresponds to transit depths of 10-100 ppm for a Sun-like host and 200-3000 ppm for a mid-M dwarf (at 0.2 R e ). For comparison, Sing et al. (2019) reach uncertainties of ∼20 parts per thousand in 1 nm bins for WASP-121b with HST/STIS while Wakeford et al. (2020) obtain 100-1000 ppm uncertainties in 10 nm bins for HAT-P-41b with WFC3/UVIS aboard HST-both at longer ultraviolet wavelengths than studied here. Thus, while the quality of atmospheric characterization obtained from Cassini/UVIS-derived transit spectra of Titan may be out of reach for HST, it may be that next-generation ultravioletcapable space telescopes-e.g., the Large UltraViolet-Optical-InfraRed (LUVOIR; Roberge & Moustakas 2018) surveyor or the Habitable Exoplanet Observatory (HabEx; Gaudi et al. 2018)-could better leverage exoplanet transit observations at ultraviolet wavelengths.