Titan’s “Average” Ionospheric Structures from Cassini

The decadal observations of Titan’s neutral atmosphere and ionosphere by the Ion and Neutral Mass Spectrometer mass spectrometer and other instruments on board the Cassini spacecraft at the Saturnian system provide a precious data set concerning the large-scale structure of Titan’s upper ionosphere. An attempt is made in this study to generate average 3D ion density distributions for different ion species by using a simple approximation of the solar zenith angle dependence. Both altitude dependence and neutral gas number density dependence of the ion density distributions will be presented. This empirical approach allows a comprehensive visualization of the global properties of Titan’s ionosphere that could be useful as engineering models for future missions to Titan.


Introduction to Titan's Ionosphere
The Voyager 1 flyby of the Saturnian system in 1980 showed that Titan's thick atmosphere with a surface pressure of 1.5 bars is composed of nitrogen (95%) and methane (4%) (Broadfoot et al. 1981;Hanel et al. 1981;Hunten et al. 1984;Tyler 1987). Another important result was the upper limit of 3 × 10 3 cm −3 of the electron density derived by the radio occultation experiment (Lindal et al. 1983). Using the number density profiles of neutral molecules from a 1D photochemical model of Titan's atmosphere by Yung et al. (1984), Ip (1990) produced a simple ionospheric model with HCNH + and C 2 + H 5 as the major ions. More sophisticated theoretical models were subsequently developed by Keller et al. (1992), Gan et al. (1992), and Keller et al. (1998), taking into account the effects of photoionization and photoelectron impact ionization. Bird et al. (1997) reanalyzed the VGR 1 radio occultation data and derived an ionospheric electron peak density of 2400 ± 1100 cm −3 just before the launch of the Cassini spacecraft.
The first result of Cassini measurements of the neutral atmospheric composition and ionospheric composition by the Ion and Neutral Mass Spectrometer (INMS) observations during the T5 encounter was reported by Cravens et al. (2006). The most abundant ion species are + CH 5 , C 2 + H 5 , and HCNH + as predicted by photochemical models. T5 was a nightside pass; the main ionization agent is magnetospheric electrons. In comparison to the primary ions like + CH 5 , heavy ions with long chemical lifetimes could be transported from the dayside to the nightside as well as production via ion-molecule reactions (Cravens et al. 2006(Cravens et al. , 2009Agren et al. 2007;Cui et al. 2009). Richard et al. (2015a) examined in detail the nightside ion production rates under different Saturnian magnetospheric conditions. Titan's dayside ionosphere is generated by photoionization and photoelectron impact ionization. A focus is therefore on the energy spectra of the secondary electron flux (Lavvas et al. 2011;Richard et al. 2011Richard et al. , 2015bMukundan & Bhardwaj 2018a). In addition, impact ionization effects of magnetospheric electrons (Snowden et al. 2013) could be important. Because of the draping of the magnetic flux tubes, magnetospheric electrons tend to precipitate on the wake side atmosphere according to the theoretical calculations of Edberg et al. (2015). The atmospheric heating rates resulting from the collisional impacts of the magnetospheric ions and auroral electrons, respectively, have been estimated to be an order of magnitude smaller than from the solar extreme ultraviolet heating (Snowden & Yelle 2014). Because of the rapid change in the plasma environment in the vicinity of Titan with a timescale of hours to days, the ionospheric number density distributions and structures could vary significantly from encounter to encounter. A good example is the large variability between T118 and T119 as discussed by Edberg et al. (2018). Below about 1300 km altitude, the ionospheric composition and structure can be described by photochemical equilibrium calculation. Above this altitude, dynamical transport effects due to magnetospheric interaction will modify the ionospheric structure with the ions of short chemical lifetimes being less affected (Cravens et al. 2009;Westlake et al. 2012).
The Ion and Neutral Mass Spectrometer (INMS) measurements of the ionospheric composition have been used to constrain theoretical models in detail. The data sets of ionmolecule reactions compiled by Anicich & McEwan (1997) and McEwan & Anicich (2007) provide inputs to the chemical networks of hydrocarbon production and loss processes. A noteworthy study by Vuitton et al. (2006) pointed out the possibility that some ion masses like CH 2 N + H 2 , C 3 H 4 N + , and C 3 H 6 N + could be produced by nitrogen chemistry involving proton exchange. On the other hand, Westlake et al. (2012) proposed that reactions of HCNH + with H 2 , C 2 H 2 , and C 2 H 4 could lead to loss of HCNH + and production of CH 2 N + H 2 , C 3 H 4 N + , and C 3 H 6 N + at the same time. From the inclusion of an additional neutral reaction, namely, N( 2 D) + HCN → CH + N 2 introduced by Hébrard et al. (2013), Dobrijevic et al. (2016) are able to fit the HCNH + number density profile very well. In a model-data comparison study, Cravens et al. (2009) showed that neutral minor species like HCN, C 2 H 2 and C 2 H 4 could be important to the ion composition and large-scale structure of Titan's ionosphere. It is noted that, among these neutral minor molecules, HCN as an important coolant can control the thermosphere temperature through a chemical feedback network (Yelle 1991;Agren et al. 2007;Bell et al. 2010). Another problem in the photochemical model calculations has to do with the error propagation caused by uncertainties in the electron dissociative recombination coefficients. For example, Shebanits et al. (2017) showed how tuning of these coefficients might allow the achievement of better agreement with the observed electron number density profiles of Titan's dayside ionosphere. The existence of different possible chemical path ways means that Titan's ion chemistry is still not fully understood.
As far as ion composition is concerned, Cassini additionally made two very exciting discoveries. The first one has to do with the detection of heavy positive ions with masses up to 350 amu and beyond by INMS and the Cassini Plasma Spectrometer's Ion Beam Spectrometer (CAPS-IBS; Waite et al. 2007;Crary et al. 2009). The second one deals with the finding of negative heavy ions with masses up to 13,800 amu (Coates et al. 2007(Coates et al. , 2009Waite et al. 2007;Wellbrock et al. 2013Wellbrock et al. , 2019. Both positive and negative heavy ions have strong implications on the formation of the dense haze layers investigated by the Cassini Ultraviolet Imaging Spectrograph experiment (Liang et al. 2007;Tomasko et al. 2008;Vinatier et al. 2010;Koskinen et al. 2011). In the case of the positive heavy ions up to the mass of benzene, their existence in Titan's ionosphere was analyzed by Wilson & Atreya (2003 by way of ion-molecule chemistry. Their formation processes were subsequently discussed by De La Haye et al. (2008) and Vuitton et al. (2008). Westlake et al. (2014) compared three different production mechanisms: namely, (1) direct ionization of ambient neutrals, (2) proton exchange of the major ion species ( + CH 5 , + C H 2 5 and HCNH + ) with ambient neutrals, and (3) ionmolecule reactions of the form, C x H y + C 2 H 2 and C x H y + C 2 H 4 , and concluded that the third process should be most promising. The Westlake et al. (2014) chemical models showed that a density peak of the C1, C2, and C3 positive ions should coexist with the ionospheric electron peak at about 1200 km, while the peak altitude of heavier positive ions (C4 to C13) should decrease as a function of the mass. A recent study by Vuitton et al. (2019) compared the Cassini measurements with theoretical model calculations using a chemical network taking into account high-resolution isotopic photoabsorption cross sections of N 2 and other factors. It was shown that the observational data can be generally reproduced by this 1D coupled ion-neutral photochemical kinetic and diffusion model.
The discovery of negative heavy ions by the CAPS-ELS instrument came as a surprise (Coates et al. 2007;Waite et al. 2007). These heavy anions distributed between 950 km and 1400 km can be divided into several mass groups: 12-30, 30-55, 55-95, 95-130, 130-190, 190-625, and 625 up to 13,800 amu (Coates et al. 2009). The number density peak (∼950 cm −3 ) is at about 1021 km altitude (Wellbrock et al. 2013). Statistical study of the CAPS-ELS measurements in many Titan encounters showed that the highest concentration can be found near the terminator and in the vicinity of the winter pole (Coates et al. 2009;Wellbrock et al. 2019). Such a global distribution of the negative heavy ions suggests that they are mainly produced on the dayside and be transported across the terminator, to the nightside, or the regions where solar insolation is weak. This observational effect is in contrast to the previous theory that the negative heavy ions and negatively charged aerosols should be generated mainly in the nightside ionosphere (Whitten et al. 2007).  There are three major source mechanisms of the negative heavy ions. They include (1) dissociative electron attachment, (2) radiative electron attachment, and (3) ion-neutral associative detachment, i.e., H − + HCN → CN − + H 2 and C 3 N + e → C 3 N − , according to Vuitton et al. (2009Vuitton et al. ( , 2019, Dobrijevic et al. (2016) and Mukundan & Bhardwaj (2018b). The negative ions like CN − and C 3 N − are precursor of the aerosols (or tholins) in the lower atmosphere (Waite et al. 2007). Lavvas et al. (2013) described the growth of the negatively charged aerosol particles from the hydrocarbon building blocks by fast absorption of positive ions. Interestingly, it was also shown that, because of the presence of the relatively large amount of negatively charged aerosols, the electron number density would be reduced by a factor of about 1.2 and the positive ion number density by a similar factor of 1.3. Titan's ionosphere is therefore a natural laboratory of highly complicated processes of coupled interactions among neutrals, ion, heavy ions, and aerosol particles of both positive and negative charge.
From the Saturn Orbit Insertion on 2004 July 15 to the final plunge into Saturn's atmosphere on 2017 September 15, the Cassini mission provided long-term measurements of Titan's ionospheric variation during one solar cycle. By comparing the time variability from solar minimum (∼2005-2006) to the solar maximum (∼2012-2015), Edberg et al. (2013) found that the peak electron number density is proportional to the square root of the solar ultraviolet energy flux as expected for photochemical equilibrium. Madanian et al. (2016) found a consistent enhancement of ion peak densities at solar maximum, especially for primary ions like + CH 3 which is a proxy for the photoionization rate of N 2 . On the other hand, the altitude of the ionospheric  density peak tends to decrease, which might indicate the "shrinking" of Titan's neutral atmosphere. Madanian et al. (2016) showed that this effect can be most clearly seen in other primary ions, like + CH 5 and suggested that this can be understood in terms of the smaller abundance of CH 4 at the time near solar maximum because of the larger photodissociation rate of this molecules (Westlake et al. 2014).
For a thin atmosphere like that of the Earth or Mars, the dependence of the peak electron density (n e ) on the solar zenith angle (SZA) is generally proportional to cos 1/2 (SZA) and there will be no solar insolation once across the terminator. Richard et al. (2015b) modeled the production rate profile of + N 2 and + CH 4 as a function of SZA. Edberg et al. (2015) examined the electron density measurements from the Radio and Plasma Wave Science (RPWS) experiments using an empirical equation to approximate the dependence of the peak n e value on SZA: where a and b are coefficients specific to the individual data sets. The best-fit average values obtained form T16 to T71 observations are a = 3198 cm −3 and b = 0.68, respectively. In order to explore the time variations of Titan's ionospheric structure and ion composition systematically, this study will apply the Edberg et al. (2015) formulation to electron and different ion species in a 2D cross-section by combining the INMS and RPWS data. This is a first attempt to produce Titan reference ionospheric models.

Observations
The data used came from the experimental observations by the INMS on the Cassini spacecraft. The INMS has two different observation modes providing neutral and ion density measurements, respectively (Waite et al. 2004;Mandt et al. 2012). In the open-source mode, the incident ions are guided through a quadrupole mass analyzer to select specific mass-tocharge ions. The INMS data consist of a sequence of counts in mass channels from 2 to 99 Dalton. We have compared the total ion number density distributions with the electron number density distributions obtained from the same Titan flybys. In most cases, good agreement can be found for altitude >1200 km. At lower altitude, the ion number densities are generally larger than those of the electrons because of the existence of a significant population of negative ions in this region (Coates et al. 2009;Wellbrock et al. 2019 Table 1. The criterion for selecting an encounter is basically the availability of a complete density distribution covering the altitude from 1000 to 2000 km. Figure 1 is a summary of the Cassini spacecraft trajectories with INMS measurements used in this study projected on the R − ρ coordinate system, where R is the radial distance of the Cassini spacecraft from the center of Titan and ρ = (y 2 + z 2 ) (1/2) is the perpendicular distance of the spacecraft to the Titan-Sun axis in the x-direction. In addition to the ion mass spectrometer . Comparison with spectra of T5, T40, and T48 flyby measurements averaged between 1000 and 1200 km. The gray area shows the dayside (T40 and T48) ion mass spectra and the black line shows the nightside (T5) ion mass spectra.
observations, important information was obtained by the RPWS instrument for electron density distribution (Edberg et al. 2010) and the CAPS instrument for both positive and negative heavy ion measurements (Coates et al. 2010). This figure highlights the fact that because of the flyby nature of the Cassini measurements, each Titan encounter can only provide information on the plasma environment along its trajectory with the information on its global structure missing. As discussed earlier, Titan's ionosphere is quite variable as demonstrated by the RPWS observations of electron density on two successive Titan encounters with very similar trajectories (Edberg et al. 2018). On the other hand, it is interesting to search for systematic three-dimensional changes in terms of diurnal and solar cycle effects.
The Titan flybys took place at a variety of possible combinations. An attempt to isolate the plasma effect of magnetospheric interaction like electron impact ionization and flow scavenging can be made by defining four regions of flyby geometry, namely, ram, wake, Saturn, and anti-Saturn. In principle, the INMS and/or RPWS data can be organized by the values of the solar zenith angle and the altitude according to the statistical binning. If Titan's ionospheric structure is assumed to be axially symmetric along the Sun-Titan line, the cumulative data can be used as inputs to 3D grid models for the electron and different ion species.
Because of the extended structure of Titan's neutral atmosphere, the usual thin atmosphere approximation of the Chapman  layer model is not applicable. The most important difference is that the solar UV photons will be able to penetrate the atmospheric region beyond the nominal terminator where the solar insolation angle SZA > 90°. In other words, the solar UV ionizing flux (F(z, λ)) is determined by the attenuation effect along the light of sight that is an integral function of the neutral density distributions and absorption cross sections of different neutral species. In the calculation of the N 2 photoionzation rate, for the sake of comparison, the vertical distribution of the neutral N 2 number density was obtained by averaging all Cassini data of the N 2 number density distribution. The spectral distribution of the F10.7 solar flux was adapted from Keller et al. (1992), and the photoabsorption cross sections for photodissociation and photoionization given in Huebner et al. (1992) and https://phidrates. space.swri.edu/ have been used in the calculation. Figure 2 illustrates the photoionization rate of the + N 2 ions. It shows the formation of a thick Chapman layer structure as well as an extension of the photoionization effect to regions with SZA as large as 135°. As a consequence, the ion densities determined by ionization equilibrium will no longer be of the form of n i ∼ cos 1/2 (SZA). In fact, because of the impact ionization effect of Titan's magnetospheric charged particles and electrons, the nightside atmosphere is subject to continuous ionization effect though at a substantially lower level. Figure 3 shows the number density distributions of ion mass 15 (i.e., + CH 3 ) obtained in different flybys between the SZA range of 0°(noon) and 180°(midnight) over a sequence of (18) altitudinal intervals from 1000 to 2000 km. It can be seen that the + CH 3 number density profile varies from encounter to encounter with the corresponding values changing by nearly one order of magnitude at the lower ionosphere below 1100 km altitude. The ion number density profiles also have relatively large-amplitude swings above 1500 km altitudes as a consequence of plasma dynamics. The subpanels illustrate the assemblage of the data points in different bins.
The dayside ionosphere is mainly affected by photochemistry. An examination of the mass spectra of different flybys including the dayside and nightside encounters could also reveal the diurnal variations of different ion species. Figure 4 compares the INMS mass spectra of T5 (nightside) and T40 and T48 (both dayside). It can be seen that there is an apparent trend for some light ions with mass between 10 and 40 in the dayside ionosphere to be more abundant relative to their counterparts in the nightside ionosphere, while some heavy ions with mass 70 display the opposite trend. This phenomenon is a global feature of Titan's ionosphere as discussed below.

Model Description
To a first approximation, the Titan ionosphere can be considered to have an axially symmetric structure. In this event, the INMS ion data can be collected into a grid system with binning in altitude and SZA. Figure 5 is a composite plot of a representative set of ion data grids in z and SZA. In each frame the data points are best fit to a statistical curve in the form of n i = a × cos(b × SZA). This means that for each ion species, its 2D ionospheric configuration between 1000 and 2000 km altitude could be described by 18 pairs of (a and b) coefficients. The parameter a determines the maximum density value in each altitude interval, and the parameter b determines the degree of density decrease along SZA.
To make the presentation of the results more comprehensive to the subsequent analysis, a polynomial fit of the parameters a and b are needed to smooth the results as shown in Figure 6. Figure 7 illustrate the patterns of the three major ion species, namely, + CH 5 (mass 17), HCNH + (mass 28) and + C H 2 5 (mass 29), respectively. It is interesting to note that while the profiles of + CH 5 (generated by + CH 4 + CH 4 → + CH 5 + CH 3 ) are generally similar to those of + CH 3 because they both represent the primary ions from direct photoionization, the spatial distributions of HCNH + and + C H 2 5 are quite different. The simple reason is that these protonated ions are produced in a number of steps of ion-molecule reactions at z ∼ 1100-1200 km thus smoothing the SZA dependence. Figure 8 shows the contour plots of the ion number density distributions in the rectangular SZA-Z coordinate system with 75 km × 5°bins at altitude above 1550 km and 50 km × 5°b ins at altitude below 1550 km. Another way to present the results is to reproject them in the cylindrical SZA-R coordinate system, where +x-axis point to the Sun (where SZA = 0) as displayed in Figure 9.
For large-mass ions like C 3 + H 5 (mass 41), HC 3 NH + (mass 52), and C 4 H 3 NH + (mass 66) as shown in Figure 7, two characteristics can be found. The first one is that the altitudes of their density peaks (z max ) are lower the larger the mass. For example, z max decreases from 1100 km to about 1000 km for mass 41 to mass 66, in comparison with the nominal value of z max ∼ 1200 km for the major ions like + CH 5 , HCNH + and + C H 2 5 . The second important feature is that the diurnal variations of the heavy ions are significantly less than those of the primary ions ( + CH 3 and + CH 5 ) and the major ion species. As discussed in Cui et al. (2009), the long chemical lifetimes of the heavy ions could be caused by a combination of ion-molecule reactions and a day-to-night ionospheric transport. The aforementioned variability can also be illustrated by comparing the number density profiles at  different SZA for different ion species according to the best-fit data-driven ionospheric models (see Figure 7 and Tables A1-A3).

Neutral Number Density Dependence
One drawback of the global 3D ionospheric model described above has to do with the fact that Titan's upper atmosphere is known to be highly variable, with the N 2 number density (n(N 2 )) varies by a factor of 3-4 at 1000 km altitude and up to 10 at 1500 km altitude (Westlake et al. 2011;Hsu & Ip 2019) between different Titan encounters. Such large-amplitude variation can happen in solar minimum and solar maximum. Because the peak of the photoionization rate is mostly controlled by the vertical distribution of the neutral gas dominated by N 2 , it is therefore prudent to examine the global ionospheric distribution in terms of n(N 2 ) as was considered in Westlake et al. (2014) and Madanian et al. (2016). One limitation is that only 13 Titan flybys have both ion density and neutral density measurements from the INMS experiment, namely, T5, T18, T26, T32, T39, T40, T48, T50, T57, T59, T65, T71, and T83). Figure 10 shows the total ion density profiles (both inbound and outbound) obtained from the 13 Titan flybys according to n(N 2 ). The variations seen are largely the result of the SZA dependence in addition to the time variability of the solar ionizing flux. The ion density distributions corresponding to different ranges of SZA are illustrated in Figure 11. The data points are best fit to a function defined by Equation (1). It can Figure 11. Total ion density variations as a function of the solar zenith angle (SZA) from 10 5 to 10 10 N 2 number density (cm −3 ). The black line shows a fitting curve of data in the N 2 number density interval. be seen that once n(N 2 ) > 10 7 cm −3 or z > 1400 km, the ion density values will have wide scattering reflecting the dynamical influence of the Saturnian magnetospheric plasma flow.
After going through a similar procedure of statistical smoothing as explained in Figure 6, a set of the (a, b) coefficients of Titan's total ionospheric density distribution with the N 2 number density as the dependable variable can now be achieved (Figure 12). These numerical data will allow us to produce representations of the total ion density distribution in different coordinate systems (see Figure 13). By the same token, the number distributions of various ion species as functions of n(N 2 ) and SZA can also be derived resulting in Figures 14-16. For completeness, the coefficients (a, b) are also given in Tables A4-A6 for future reference.

Discussion and Summary
In this study, we explore the average 3D structure of Titan's extended ionosphere by considering the INMS measurements from different flyby trajectories. Using an empirical expression originally used by Edberg et al. (2015) for the electron density distributions, a set of semi-quantitative ionospheric models for different ion species can be generated. The results are useful to discern the diurnal and latitudinal variations of the ion densities of different chemical characteristics. These include the primary ions like + CH 3 and + CH 5 that reflect mainly the photoionization (and electron impact ionization), the major ions like HCNH + and + C H 2 5 , and the heavy ions produced by ion-molecule reactions. Another way to examine the global ionospheric structure is to base on the neutral gas number density instead of the vertical altitude. This is because the Cassini INMS observations have shown clearly that Titan's upper atmosphere could be highly variable between consecutive Spacecraft encounters (Westlake et al. 2011;Hsu & Ip 2019), the neutral gas number density, say that of N 2 , could therefore be a more reliable parameter (Westlake et al. 2014;Madanian et al. 2016). This approach was applied to a data set of 13 Titan flybys without taking into account the variation of the solar ionizing flux.
As discussed in Waite et al. (2007), the production of tholins and heavy negative ions detected by CAPS (Coates et al. 2007(Coates et al. , 2009Wellbrock et al. 2013Wellbrock et al. , 2019 was initiated by the ionization and photochemical processes at high altitudes (∼>1000 km). In turn, the presence of tholins and large amounts of negative ions below the altitude limit of Cassini could have a significant impact on the thermal structure and density distribution of the neutral atmosphere. There could   therefore be quite intriguing coupling between Titan's neutral atmosphere and ionosphere at different timescales from day-today variation to solar cycle.
Strictly speaking, the formation of Titan's ionosphere must be the combined effect of photoionization, magnetospheric electron impact ionization, pickup ion sputtering, and most importantly the structure of Titan's nitrogen-rich neutral atmosphere. Even during solar quiet time, Titan's rotation and orbital motion would introduce time-dependent phenomena into the global structure of the ionosphere. Because of the limited number of Titan encounters with INMS measurements and the relatively sparse spatial coverage in terms of the orientation of the corotating plasma flow and the Saturn-facing angle, very little is known in this regard. As far as the general features are concerned, a wealth of information probably can still be derived by a joint study of the CAPS and RPWS data from Cassini. However, many of the details of the physical properties of Titan's ionosphere must be deciphered by a Titan orbiter mission in future.
We thank the reviewers for useful comments. This work was supported by grant No. 107-2119-M-008-012 of MOST, Taiwan.

Appendix
In the paper, the coefficients (a, b) are used to describe the best fit to a statistical curve in each ion data grid. Tables A1-A3 list the results of the coefficients (a, b) in the altitude-SZA grid system, and Tables A4-A6 list the results in the N 2 density-SZA grid system. Figure 16. Ion density distributions in the cylindrical SZA-N 2 density coordinate system. The ion species included are the same as in Figure 14.