Measurements of Titan’s Stratospheric Winds during the 2009 Equinox with the eSMA

Saturn’s moon Titan possesses stratospheric zonal winds that places it among a sparse class of planetary bodies known to have superrotation in their atmospheres. Few measurements have been made of these speeds in the upper stratosphere, leaving their seasonal variations still not well understood. We examined observations made with the extended Submillimeter Array in 2009 March (L s = 355°) and 2010 February (L s = 5°), shortly before and after Titan's northern spring equinox. Cassini observations and atmospheric models find equinoctial periods to be especially dynamic. Zonal wind calculations, derived from the Doppler frequency shift of CH3CN near 349.4 GHz, yielded speeds of 128 ± 27 m s−1 in 2009 and 209 ± 48 m s−1 in 2010. We estimated the measured emission to originate from vertical altitudes of 336−88+112 km, equivalent to pressures of 3.8−3.4+19.2 Pa, commensurate with Titan’s upper stratosphere/lower mesosphere. This suggests a possible increase in zonal speeds during this period. The results are then compared to those from previous Cassini-inferred and direct-interferometric observations of winds, as well as general circulation model simulations, to form a more complete picture of the seasonal cycle of stratospheric zonal winds.


INTRODUCTION
Corresponding author: Siobhan Light slight@umd.eduTitan's atmospheric chemistry surpasses in complexity any other known atmosphere in our solar system (Hörst 2017).The atmosphere contains a wealth of organic molecules including hydrocarbons and nitriles (e.g., Wilson & Atreya 2004).Over the years, these molecules have been observed to vary both spatially and temporally.Various studies have demonstrated this observationally (see, e.g., Teanby et al. (2019); Coustenis et al. (2020); Mathé et al. (2020); Vinatier et al. (2020)) by comparing limb composition measurements taken by Cassini/Composite Infrared Spectrometer (CIRS) of various molecules throughout Titan's latitudes over the course of the mission from 2004 to 2017.
One agent responsible for the fluctuations in Titan's atmosphere over time is the seasonally changing solar insolation.The length of a Titan year is the same as Saturn's at 29.5 Earth yr, and Saturn has a similar axial tilt to Earth of about 26.7°.The changing direction of solar forcing drives a seasonally variable atmosphere on Saturn, and likewise Titan which has a similar tilt to the ecliptic.Seasonal changes were observed with the Cassini mission, but since Cassini was not able to witness a full Titan year or sample all atmospheric levels, more observations are needed to construct a complete picture of the seasonal states of the atmosphere (Teanby et al. 2019).
When it comes to Titan's general atmospheric circulation, the moon generally has a single poleto-pole Hadley-like cell, but near the equinoxes, the circulation reverses, resulting in the temporary creation of two equator-to-pole cells (Hourdin et al. 1995;Newman et al. 2011;Teanby et al. 2012;Lombardo & Lora 2023a).Titan's last equinox (2009)(2010) occurred on 2009 August 11, and increasing trace gas abundances were observed over Titan's south pole following this event (Teanby et al. 2012).Titan's equinoctial periods appear to be a high activity point for Titan's meteorology.For example, Rodriguez et al. (2018) reported Cassini observational evidence for active dust storms on Titan during the 2009 and 2010 time period.
Titan's stratosphere and mesosphere (collectively known as the middle atmosphere) exhibit yearround, strong prograde (west-to-east) zonal winds.These winds are in balance with meridional temperature gradients via the gradient wind relation (Sharkey et al. 2021) and are driven by adiabatic cooling above the polar regions by thermal emission.
One molecule detected numerous times in Titan's atmosphere due to its prominent spectral character in submillimeter astronomical observations has been acetonitrile (CH 3 CN).This organic nitrile was first observed with the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope (Bézard et al. 1993).Following this, Marten et al. (2002) used IRAM observations of several CH 3 CN transitions to retrieve a disk-averaged vertical profile of the molecule up to 500 km and also demonstrated that CH 3 CN's distinct lines made it a good candidate molecule to probe the middle atmosphere of Titan.
In recent years, the Atacama Large Millimeter/submillimeter Array (ALMA) has ushered in a new era of more accurate, direct-wind and molecular-abundance measurements on Titan.These new measurements have increased understanding of the latitudinal variability of acetonitrile (e.g., Cordiner et al. 2019;Thelen et al. 2019), but its temporal variability over long periods of time remains to be constrained.Furthermore, several years ago (∼L S =82°), an unexpectedly intense wind jet incongruent with previous theoretical models was detected with ALMA in the thermosphere (up to 1000 km), though winds were detected all the way from the stratosphere through the mesosphere and up to the thermosphere using different molecules and transitions (Lellouch et al. 2019).Subsequent work investigated the mechanisms behind the jet, and it highlights how Titan's atmospheric dynamics are highly temporally variable and altitude dependent (Cordiner et al. 2020).With that in mind, work remains to be done to discern the long-term seasonal variability of Titan's winds.Doppler measurements of Titan's winds have been more frequently measured in the ALMA era, but Titan's long seasonal cycle means that observations made before ALMA remain valuable for addressing this question.In this study, spatially resolved interferometric observations made by the extended Submillimeter Array (eSMA) taken around the time of Titan's last equinox were analyzed in order to see (1) how CH 3 CN emissions changed in the short period preceding and following Titan's equinox and (2) the winds on Titan at this time.Both of these atmospheric variables will be examined in the context of previous work to better understand how they have changed on Titan over time.

OBSERVATIONS
Dedicated interferometric observations of Titan with the eSMA occurred both before and after the Titan equinox, on 2009 March 23 (L s =355°) and 2010 February 12 (L s =5°), respectively (Gurwell et al. 2009;Light et al. 2021).The Submillimeter Array (SMA), located near the summit of Maunakea, is an array of eight 6 m diameter antennas operating at frequencies between 200 and 420 GHz.The eSMA was a joint project to combine the SMA with two other nearby submillimeter facilities into a larger interferometer, the 15 m diameter James Clerk Maxwell Telescope (JCMT) and the 10.4 m CalTech Submillimeter Observatory (CSO; Bottinelli et al. (2008); see Figure 1).
During the 2009 observations, the eSMA was composed of all eight SMA antennas, the JCMT and the CSO.However, in 2010, only 7 SMA antennas were available, along with the JCMT, providing less overall collecting area and a shorter maximum baseline (see Table 1).For both sets of observations, the weather was very good.Additional observational details are outlined in Table 1.In both years, the eSMA was tuned to capture emission from acetonitrile, specifically the CH 3 CN J = 19-18 rotational band, near 349.4GHz.The strongest transitions (k = 0, 1, 2, 3) were measured with ∼ 0.203 MHz resolution within a ∼91 MHz wide band (448 channels), along with another 1600 MHz wide window measured at 3.25 MHz resolution which provided a continuum signal around the CH 3 CN transitions.
In 2009, there was no real-time Doppler tracking of the changing velocity difference between the observatory and Titan, and alignment of spectra (obtained at 30 s intervals) was performed with offline processing.However, in 2010, real-time Doppler tracking was performed by adjusting the local oscillator (LO) frequency, and offline corrections were not needed.

Data Calibration
The interferometric data were calibrated using the Millimeter Interferometer Reduction (MIR) package, a complete suite of tools for calibration of eSMA data, for use in IDL1 .The calibration steps were mostly straightforward and standard; we highlight a few important aspects here.
Initial calibration of complex (amplitude and phase) gain corrections with time were measured using periodic observations of a nearby point source, typically millimeter-wave strong blazars; in 2009 the source was J1058+0133 (4C +01.28), and in 2010 J1229+0203 (3C 273).The complex spectral passband response was measured in both years using 3C 273.The signal-to-noise ratio (S/N) of the data (further detailed in the results section) was found to be high enough that self-calibration of complex gains over time could be completed using the established Titan continuum.In doing so, short-term atmospheric instability was removed from the results, further enhancing their quality.Moreover, the flux density scale was set using the Titan continuum, which is accurate to 3-5% in the submillimeter (Butler 2012).
The calibrated complex interferometric data ("visibility data") of Titan for a broad-frequency continuum band and in each spectral channel were then translated into the UVFITS format.Some further adjustments to the visibility data were applied: 1.A Doppler correction to each time step of the 2009 visibility data was performed to remove the smearing effects of Earth rotation, Titan's orbital motion, and other velocity variation with time, aligning the spectral correlator coverage at each time step with the Titan frame of reference.
2. The 2010 data set was adjusted such that the apparent diameter of Titan (0.8156") was scaled to the 2009 data (0.8422"), for ease of comparison between the years.This was accomplished by scaling the (u,v) coordinates of each visibility by the ratio of the 2009/2010 diameters, and scaling the amplitude of each visibility by the square of the ratio of the 2009/2010 diameters.This maintains the same 2010 synthesized beam relative to Titan's diameter as well as providing a consistent flux scale and S/N, but directly comparable to the 2009 data.
3. The north pole position angle of Titan for each year was about 356 degrees (measured east of north), amounting to a slight rotation of Titan north relative to J2000 sky north.The north pole was rotated modestly (through rotation of the observed (u,v) positions) such that Titan north and sky north are aligned.
4. For ease of imaging, the Titan continuum data for each visibility point was subtracted from the spectral line data, creating a continuum-subtracted visibility data set; these data represent the spectrum of CH 3 CN relative to the continuum (this manifests entirely as emission in these observations).
The visibility data were then imported into the Astronomical Image Processing System (AIPS).Images of Titan in the continuum and in each spectral channel (continuum-subtracted) were performed using a modified "clean" (deconvolution) process, using the task IMAGR.The pixel size was set at 0.02".The spectral line and continuum images were both reconvolved to produce images with shapes reflecting the naturally weighted beam ("natural" beam): a 0.384" × 0.196" ellipsoid with a position angle (PA) of 48.78 • for 2009 and 0.382" × 0.228" with a PA of 53.83 • for 2010).Via a sigma-clipping function, the mean RMS noise on the integrated flux maps for 2009 and 2010 for the natural beam were measured to be 0.0191 Jy bm −1 and 0.0187 Jy bm −1 , respectively.The final results are shown in the continuum-subtracted maps of Titan shown in Figure 2.

RESULTS
Spectra were extracted from the continuum-subtracted images at frequencies near 349.4GHz and the results of this over the disk of Titan are shown in Figure 3.The center of the peak emissions was compared to values for the J = 19-18 rotational band transitions of CH 3 CN listed in the JPL catalog (Pickett et al. 1998) and Cologne Database for Molecular Spectroscopy (CDMS) (Müller et al. 2001(Müller et al. , 2005;;Endres et al. 2016).The two databases effectively agreed on frequency values and associated errors, with the largest error on these values equal to 300 Hz, making it negligible compared to the channel resolution of about 0.203 MHz (∼175 m s -1 ).
Typical thermal profiles for Titan's atmosphere (Lellouch et al. 2019) indicate that thermal broadening dominates the line core shape (Cordiner et al. 2020).As in Moreno et al. (2005), Cordiner et al. (2020), and Lellouch et al. (2019), pressure-broadened line wings also were observed in the line emission.After fitting various line shapes to the spectral peaks, a Moffat distribution (Moffat 1969) was found to more accurately represent the complex line emission from Titan than either Gaussian, Lorentzian, or Voigt line shapes, and was consequently used for determining the Doppler shifts of the observed line shapes, as shown in Figure 3.The equation for the Moffat distribution is with four parameters: amplitude (A), center (µ), width (σ), and exponent (β).The Moffat distribution was able to model both the line core and wings of the CH 3 CN more precisely than other profile types.
To measure the Doppler shift and determine zonal winds, the calculations were performed using seven rectangular regions taken from across the disk of Titan, highlighted in Figure 2. The reason large areas needed to have their emissions integrated together to resolve zonal winds is because the interferometric observations described here, with natural beamwidths approximated as elliptical Gaussians of the FWHM ∼0.38" × 0.20", only coarsely resolve the disk of Titan (∼0.8" diameter, see Table 1).Half of the integrated flux of a Gaussian beam comes from outside the central FWHM, meaning that emission from well outside the central region of the beam contributes to the image at a given pixel.This is especially relevant for these CH 3 CN observations, where the emission (even at the coarse resolution) is strongly concentrated at the limb of Titan as seen in the integrated beam maps (Figure 2).This is a result of submillimeter lines being generally optically thin (e.g.Cordiner et al. 2020) and presenting the strongest thermal emission where path lengths in the atmosphere are largest (i.e., at the limbs).The combination of relatively large beam size and highly concentrated emission means that careful consideration of the beam shape is important.For example, spectra from the disk center show apparent emission, but this is largely limb emission captured by the coarse beam (mostly along the long axis of the beam), not emission from the disk center.The rotational model described below takes this effect into account.
Spectra across latitudes in each of the rectangular regions were then summed together and adjusted by a geometric scale factor determined by where a and b are the beam major and minor axis respectively and ∆l is pixel size.After including continuum error, the spectra were fit at each of the four transitions using the LMFIT package for Python (Newville et al. 2016).The resulting fits for the easternmost region are shown in Figure 3b for the 2009 data.The shift between limbs within the single J = 19-18 rotational band is expected to be essentially the same since they probe the same altitude range (see Moreno et al. (2005)).Calculations performed in Subsection 3.1 also support this choice.Consequently, due to the challenges of interpreting data in this low S/N regime, sensitivity is enhanced by calculating a single frequency shift by finding a weighted average of the individual transition shifts.These frequencies were then converted to velocities via the Doppler shift formula.
To calculate uncertainty on each regional velocity, the same fitting method was applied to 20,000 sets of synthetic Gaussian data based on the spectral resolution, noise, and shape of the original spectral data.The subsequent synthetic data fitting results were filtered to exclude velocities that were outside 4σ so that extreme synthetic data results did not affect the final error calculations.To determine how wide a spectral region to fit over, the fitting methodology and subsequent error determination was tested on regions of 1.015 MHz to 3.453 MHz, limited on the lower end by the spectral FWHM and thermal broadening and on the higher end by the proximity of the 19 1 − 18 1 and 19 0 − 18 0 transitions.After doing so, the root-mean-square errors (RMSE) of the predicted shifts in the synthetic data sets was used to estimate the error on the velocity measurements.The lowest errors and most consistent velocity values were associated with using ± 1.625 MHz (equivalent to ± 8 channels of resolution).Figure 4 displays velocity values from each of the seven rectangular regions, represented for 2009 as blue squares and 2010 as red squares, and their associated error bars.
As mentioned above, the large beam size means that the measured velocity will be offset relative to the true equatorial east-west zonal wind.To account for this, a simple model was made based on a fully zonal wind approximating solid body rotation.Zonal wind speed scales as a function of latitude as V z cos(ϕ), where V z is the equatorial maximum zonal wind speed and ϕ is latitude.Approximation of a roughly solid body zonal wind model is appropriate here, considering the relatively large synthesized beam, that the results of Lellouch et al. (2019) and Cordiner et al. (2020) found a strong zonal wind system at lower latitudes, and that for a viewing geometry with a sub-Earth latitude of zero, 81.8% of the apparent disk is from ±45 deg ("low") latitudes.This distribution was used to generate a smeared velocity map, allowing beam effects to be accounted for when calculating zonal velocity.However, this relied on the assumption that smearing the calculated zonal wind radial velocity by weighting the measured emission and the beam yields the same result.The model also captures how the relative velocity of CH 3 CN emission is averaged within each of the seven rectangular regions (centered on Titan) are affected by the distribution of emission.As the emission maps in Figure 2 show, the model also assumes emission is weaker towards the south pole and stronger towards the north pole.The velocity values in each rectangular region are hence adjusted to reflect the relative velocity computed by the model.Subsequently, we utilize a curve fitting approach to model Titan's solid body rotation, shown in Figure 4.The fitting algorithm uses scipy.optimizemodule from the SciPy library (Virtanen et al. 2020) to minimize using the Nelder-Mead algorithm (Gao & Han 2012) the weighted sum of squared differences differences between the observed velocities and those predicted by the velocity curve.To accommodate potential Doppler tracking errors, we allow for flexibility in the "zero" velocity, rather than fixing it.Thus, the model we were optimize fitting was where V is the predicted velocity, Z denotes the global zonal wind value, r represents the relative velocity of the zonal wind at a given distance in the solid body rotation model, and ω signifies the velocity offset at the center of Titan.Essentially, we are fitting the one-dimensional rectangular velocity data to a solid body rotation curve that acts as a proxy for global zonal winds.Furthermore, we subtract the rotation speed of Titan at the approximate stratospheric altitude of the data, as it is not relevant to the observed zonal winds.Directly fitting the data yields an estimate of approximately 138 m s -1 and 249 m s -1 for the global zonal winds in 2009 and 2010, respectively.However, uncertainties exist in both the distances and velocities associated with each data point, compounded by the limited number of data points available for fitting.To address this, we generated two million sets of velocities and distances, assuming a random normal distribution associated with each point's average distance and measured velocity.Subsequently, these synthetic datasets were fitted repeatedly using the same method, allowing us to derive the average and standard deviation of the synthetic data model parameters.These statistics were then used to calculate the global zonal winds and associated uncertainties.
By this method, the zonal wind speed for 2009 was calculated to be 128 ± 27 m s -1 with an offset of -0.4 m s -1 ± 16 m s -1 .In 2010, the zonal wind speed increased to 209 ± 48 m s -1 with a larger offset of ∼61 m s -1 ± 27 m s -1 .We found that the wind strengthened by a factor of 1.6 ± 0.5, with a 3.2 σ confidence.Overall, the weighted average zonal wind speed for this 2009-2010 time period is 169 ± 27 m s -1 .

Altitude Calculation
Using the Non-linear optimal Estimator for MultivariatE Spectral analySIS (NEMESIS) radiative transfer code (Irwin et al. 2008) and similar methods to Thelen et al. (2019), we calculated relevant functional derivatives from 150 to 1200 kilometers in altitude for a spectrum extracted with the "natural" 2009 beam shape at the Eastern limb.CH 3 CN spectral line parameters were incorporated from the Cologne Database for Molecular Spectroscopy (CDMS) (Müller et al. 2001(Müller et al. , 2005;;Endres et al. 2016) and the HITRAN catalog (Gordon et al. 2022).N 2 -broadening parameters were adopted from Dudaryonok et al. (2015).In turn, this enabled us to estimate the altitude and pressure levels to which the subsequent wind measurements were primarily derived from.Figure 5 shows the normalized contribution from an altitude of ∼150 to ∼1000 kilometers since contributions below or above these altitudes are negligible.The functional derivatives (and thus the extracted spectra) are not sensitive to altitudes below 150 km for reasons outlined in Thelen et al. (2019) based on calculations made from 3 a priori profiles (Marten et al. 2002;Dobrijevic & Loison 2018;Loison et al. 2015) demonstrating that CH 3 CN emission at these frequencies is not sensitive to gas abundances at these low altitudes.
Taking the average of the peak contributions at the four emission line peaks, and using the same weights used to later calculate the Doppler shift, the emission appears to have originated approximately at vertical altitudes of 336 +112 −88 kilometers corresponding to pressures of 3.8 +19.2 −3.4 Pa.These altitudes correspond to Titan's upper stratosphere/lower mesosphere and are similar to previous studies involving CH 3 CN (e.g.Moreno et al. 2005;Lellouch et al. 2019;Thelen et al. 2019;Cordiner et al. 2019)).Although there are secondary peaks in the upper atmosphere, the contribution is orders of magnitudes less than the stratosphere/mesosphere region.Despite seasonal changes in latitudinal concentration, since the molecule has a long photochemical lifetime in the middle atmosphere (Wilson & Atreya 2004;Loison et al. 2015;Vuitton et al. 2019) vertical abundances likely do not strongly vary over the course of a Earth year.In particular, according to Vuitton et al. (2019), CH 3 CN has a photochemical lifetime of about 91 years at 300 kilometers.From 2012 to 2015, there were also no major differences between the northern and southern hemispheres with regard to upper stratosphere CH 3 CN abundances (Thelen et al. 2019).Further work in Cordiner et al. (2019) concluded the molecule appears to be relatively well-mixed at these stratospheric altitudes compared to shorter-lived nitriles and hydrocarbons.

DISCUSSION
Cassini observations around this time period (Teanby et al. 2010) observed trace gas abundance to be enriched at the north pole and depleted at the south pole.The emission maps in Figure 2 are consistent with these observations and show that little significant change occurred in CH 3 CN emissions between 2009 and 2010.Acetonitrile remained concentrated primarily in northern latitudes during these observations, similar to the 2013 ALMA observations described in Thelen et al. (2019), supporting the idea that confinement by the northern winter polar vortex persisted after the equinox (Teanby et al. 2019) and due to the relatively long photochemical lifetime of CH 3 CN in Titan's stratosphere.Close to the second set of observations reported here, Cassini CIRS in January 2010 observed polar warming at high altitudes soon after equinox but no large increases in abundances until approximately a year later (Teanby et al. 2012).There appears to be a slightly higher concentration of acetonitrile in the east of Titan prior to the equinox than after, although noise in the data makes this an uncertain conclusion.Initial analysis of the 2009 data found an increasing CH 3 CN abundance with altitude (Gurwell et al. 2009).
Table 2 shows how the zonal wind measurements in this study compare to others made with similar Doppler shift methodology at comparable stratospheric altitudes (Hubbard et al. 1993;Kostiuk et al. 2001;Sicardy et al. 2006;Kostiuk et al. 2005;Moreno et al. 2005;Lellouch et al. 2019;Cordiner et al. 2020).Consistent with all of these previous observations, the 2009 and 2010 eSMA measurements show a prograde zonal wind flow and support the conclusion of superrotation of Titan's middle atmosphere.
To contextualize Titan's seasons, its northern summer solstice (L s =90°) occurred in February 1988 and May 2017, northern fall equinox (L s =180°) in November 1995 (next in 2025), northern winter solstice (L s =270°) in November 2002, and northern spring equinox (L s =0°) in August 2009.The first direct wind measurements were derived from data taken during a 1989 occultation and noise precluded measurements being known more precisely than to an order of magnitude of around ∼100 m s -1 (Hubbard et al. 1993).Kostiuk et al. (2001) combined observations from August 1993, October 1995, and September 1996 to derive speeds of 210±150 m s -1 peaking at around 210 kilometers in altitude (notably lower than the level probed in this study).Kostiuk et al. (2005) and Sicardy et al. (2006) also probed lower in the stratosphere, deriving winds of 190±90 and ∼100-150 m s -1 , respectively.
Our measurements from 2009 and 2010 suggest a strengthening of Titan's stratospheric zonal winds over its equinox.The strength of the 2009 zonal winds are comparable to those observed half a season prior during northern winter in 2003/2004 (Moreno et al. 2005), and are significantly weaker than those measured shortly before summer solstice in 2016 (Lellouch et al. 2019;Cordiner et al. 2020).The measurements before the northern spring equinox from the eSMA do appear to be lower than those taken at the summer solstice in 2017 (Cordiner et al. 2020) while the measurements taken after the equinox in 2010 are more similar in value.The results from 2009 generally suggest that the equinoctial period may be one of the slowest points in the upper stratospheric jet found to date.Cordiner et al. (2020) suggested that the fastest zonal jets may occur in the years leading up to Titan's solstices and slowing until equinoxes.The results presented here remain consistent with this idea and further suggest that the time period before Titan's northern spring equinox may be a low point in the speed of stratospheric zonal jets.Moreover, the overall 2010 zonal wind speed result appears to be closer in magnitude to those derived shortly before the 2016 summer solstice as reported in Lellouch et al. (2019) and Cordiner et al. (2020).This suggests a rapid increase in the speed of the upper stratospheric jet shortly after the northern equinox.However, more measurements must be made across Titan's seasonal cycle and at different altitudes (particularly the mesosphere and thermosphere) to definitively conclude the seasonal variability of Titan's zonal jets.For example, it would take additional data points in the time period between these 2010 observations and the subsequent 2016 observations to say whether the stratospheric jet maintained its speed through Titan's spring, or whether its speed varied over the course of northern spring.
General circulation models (GCMs) are a useful tool in simulating and diagnosing large scale climate systems in planetary atmospheres.Over the years, several general circulation models (GCMs) of Titan's atmosphere have been developed (Hourdin et al. 1995;Newman et al. 2011;Lebonnois et al. 2012;Lora et al. 2015).Generally, GCMs suggest that zonal winds weaken during winter until hemispheric reversal (Newman et al. 2011;Lora et al. 2015), an idea the low-speed measurement in 2009 would support.Around equinox, upwelling at lower latitudes, closer to the region examined in this study, is suggested by multiple models to bring eastward angular momentum to the upper stratosphere (Newman et al. 2011;Lombardo & Lora 2023a).During equinox, Lombardo & Lora (2023a) observed that the dominant zonal jet shifts from the previously winter hemisphere to the hemisphere entering fall.The 2010 results would support their suggestion that zonal winds reach a high value following the equinox.However, we do see a general increase in disk averaged zonal winds from just before to just after equinox.We compared our results to those of TitanWRF (Newman et al. 2011) and the Titan Atmospheric Model (TAM) (Lora et al. 2015) to attempt to constrain the timeline of the changes in the middle atmosphere zonal winds near equinox.
Recently, Shultis et al. (2022) used TitanWRF simulations from Newman et al. (2011) to constrain the dynamics of Titan's winter stratosphere, a time period the results presented here complement.This model was found to be in general agreement with CIRS observations at pressures above 1 Pa, enabling direct comparisons to be made.Overall, the model tended to have higher simulated temperatures relative to observations but captured general observed temperature and zonal wind trends.In general, the model simulated a minimum in mid-latitude stratospheric zonal winds between winter solstice and spring equinox, a result potentially complementing the low speed observed with the eSMA in 2009.In Figure 6a, the simulated zonal winds of the TitanWRF GCM at the time of the 2009 observations are shown.In this time period shortly before the equinox, this GCM simulates at comparable altitudes a peak zonal speed of ∼150 m s -1 in the northern hemisphere (∼17.5°) and equatorial zonal wind speed of about ∼145 m s -1 .Little appears to directly change by the time of the 2010 observations, with the GCM simulations of this time period shown in Figure 6b.With regards to the seasonal trends of zonal winds, Shultis et al. (2022) simulates for both hemispheres around 10 Pa the winds weaken and move towards low latitudes during fall and reach a minimum around mid-winter.To help understand why this happens, this complements how CIRS observed that the winter pole for the upper stratosphere/lower mesosphere was warmer than lower latitudes likely caused by adiabatic heating of the descending air of meridional circulation (Achterberg et al. 2008(Achterberg et al. , 2011;;Teanby et al. 2019;Vinatier et al. 2020;Sharkey et al. 2021;Shultis et al. 2022;Achterberg 2023).Following this mid-winter weakening, TitanWRF suggests that the zonal winds strengthen slightly and shift poleward just before each hemisphere's respective spring equinox.Following spring equinox, zonal winds are generally thought to increase, corresponding to upper stratosphere angular momentum reaching a maximum (Newman et al. 2011).This is due to cooling from less adiabatic compression/warming from less descending air, which in turn led to a steeper meridional temperature gradient, driving stronger winds.Overall, these results support the increase in zonal winds after equinox, but suggest it may occur more rapidly than the TitanWRF simulation appears to predict.Relative to the southern hemisphere, winter in the northern hemisphere tends to drive stronger general circulation, with asymmetries between the two driven by differences in insolation caused by Titan's orbital eccentricity.As a result, we would expect increases in comparable low-latitude zonal wind speeds following the upcoming 2025 equinox, albeit not as large.
TAM was recently updated to include seasonally varying radiative species (Lombardo & Lora 2023a), which tend to cool the polar latitudes due to increased longwave emission from molecular enrichment.The polar cooling also leads to slightly increased stratospheric zonal wind speeds.Like TitanWRF, TAM simulates a degree of seasonal symmetry between hemispheres.Figure 6c and 6d show TAM zonal wind simulations for the 2009 and 2010 observation times, respectively.Relative to TitanWRF, TAM simulated less change across this time period than TitanWRF, but TAM suggests zonal winds are in the process of strengthening at the time of the observations, peaking part-way through the northern spring (after the 2010 observations).Within the equatorial region examined in this study (roughly latitudes of 30°S-30°N), TAM simulations show that zonal winds appear to peak shortly after equinoxes and reach minimums shortly after solstices.This seasonal variation is more extreme at the upper stratosphere/lower mesosphere level (∼3 Pa) than lower in the stratosphere (∼50 Pa) (Lombardo & Lora 2023a).The eSMA observations also suggest the winds are strengthening, but at a more rapid pace than predicted by TAM.In connection with temperature, this variable appears to generally peak around northern winter solstice followed by an uneven decrease until northern summer solstice.Around northern spring equinox, there appears to be a small temperature dip (Lombardo & Lora 2023b).
One way to see how the GCMs simulate changes in zonal wind over this study's time period is to calculate from their model outputs a disk averaged (considering area via cosine weighting) zonal profile.Calculating the cosine weighted disk average of the zonal winds modeled by TitanWRF as shown in Figures 6a and 6b, corresponding to the timing of the 2009 and 2010 observations, respectively, reveals that TitanWRF simulated a 10% decrease in the overall zonal winds between these two time periods, unlike what the observations suggest.Meanwhile, taking the cosine weighted disk averages of Figures 6c and 6d finds TAM simulated a roughly 8% increase in the zonal winds between the two observation periods.The general increase in zonal winds suggested by the eSMA results during the observed time period would therefore be consistent with TAM simulations, but appears to be of a larger magnitude.
During the same time period as the measurements made here, Cassini observations saw zonal winds in the northern hemisphere gradually decrease as northern winter progressed towards summer, while in the southern hemisphere zonal winds began to increase rapidly after equinox (Sharkey et al. 2021).Overall, higher zonal wind values around L s = 0°would have been driven by higher values in the northern hemisphere.The limited spatial resolution of the eSMA data makes differentiating between the northern and southern hemisphere difficult.The higher zonal wind speeds around equinox includes changes in the wind speeds from both the northern and southern hemispheres.This means the strengthening we observed shortly after equinox likely includes strengthening in the southern hemisphere and weakening in the northern hemisphere, similar to the southern mid-latitude strengthening and slight weakening observed in (Sharkey et al. 2021).However, the thermal wind measurements they reported left winds unconstrained equatorward of about 20 degrees in either hemisphere at the altitudes we are sensitive to, so we are unable to make comparisons as to the finer low latitude meridonal structure of the zonal winds around equinox.

CONCLUSION
New zonal wind measurements of Titan's upper stratosphere shortly before and after the moon's northern spring equinox were made using eSMA observations and Doppler shift techniques.The results appear consistent with the idea that zonal wind speeds decrease during winter on Titan followed by a rapid increase in speed following the equinox and may help inform future GCMs to better understand Titan's atmosphere.Although the spatial resolution of these eSMA measurements are much more coarse than those offered by ALMA in recent years, they nonetheless offer insight into this dynamical time period in Titan's atmosphere.Continued measurements of Titan's atmosphere with ALMA and other sub-millimeter observatories are still needed to form a more complete picture of the seasonal variation in the moon's circulation.For example, knowledge of interannual variability is currently limited to models, so it can not be concluded if the seasonal patterns described here remain constant over future Titan years.By understanding the dynamics of Titan's atmospheric superrotation and circulation, insights have previously been made into the atmosphere of Venus and may also help to understand those of exoplanets as Doppler observations begin to be made of hot Jupiters (Imamura et al. 2020).Improvements in knowledge of Titan's middle atmosphere will complement atmospheric modeling that may support the upcoming Dragonfly mission design (Lorenz et al. 2021;Lorenz 2021).

ACKNOWLEDGMENTS
The material is based upon work supported by NASA under award number 80GSFC21M0002.SLL and CAN received funding for this work from the NASA Solar System Observations (SSO) Program and the NASA Astrobiology Institute.The development of the eSMA was facilitated by grant 614.061.416from the Netherlands Organization for Scientific Research, NWO, and NSF grant AST0540882 to the Caltech Submillimeter Observatory, along with the dedicated work of observatory staff from the JCMT, CSO, and the SMA.The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.We recognize that Maunakea is a culturally important site for the indigenous Hawaiian people; we are privileged to study the cosmos from its summit.Thank you to Erin Flowers insight into TitanWRF.Thank you to Martin Cordiner for help with the code to plot the acetonitrile emissions and fruitful discussions about previous work related to Titan's zonal winds.

Figure 1 .
Figure 1.(a) Photograph of the observatories involved with the eSMA, left to right: SMA, JCMT, CSO.(Photo Credit: Jonathan Weintroub) (b) Map view; small green circles indicate the location of SMA antennas, yellow triangle represents JCMT location, and the pink diamond represents CSO location.(Photo Credit: Google Earth)

Figure 3 .
Figure 3. (a) Moving average of integrated flux over a radius 0.60" centered on Titan, smoothing together three channels at a time.Each channel has a resolution of 0.203125 MHz.The spectrum from 2009 is in purple and 2010 is in black.Individual transition lines are noted with dashed lines since Doppler shift measurements were made relative to these CH 3 CN J = 19 k -18 k (with k=3,2,1,0) transitions.(b) 2009 spectra from the east limb of Titan, with fits overlapped in dashed lines (c) Fit of the 19 3 -18 3 component of the 2009 eastern limb spectra, with a red line showing the center of the fit; the blue, dotted line shows the rest frequency of the transition, and the turquoise, dashed line shows the model fit (d) Similar to (c) but for the western limb

Figure 4 .
Figure 4. Comparison of the measured Doppler shifted velocities from seven rectangular regions spanning all latitudes on Titan, along with their associated errors.Blue squares represent 2009 data, while red circles depict 2010 data.Dashed lines, optimized to align with solid body rotation acting as a proxy for zonal wind, are fit for each respective year.Shaded regions indicate the 95% confidence interval (±1.96σ) on the velocity fits.The SMA beam size corresponds to the east-west projection of the major axis of the 2009 "natural" beam shape.

Figure 5 .
Figure 5. Normalized functional derivatives calculated using NEMESIS of the spectral radiance with respect to 2009 CH 3 CN (J = 19 -18, ν = 0) abundance, pressure/altitude, and frequency.Absolute contour values are shown in side bar (e.g., the white contour shade corresponds to the highest contribution, ∼1.000).

Figure 6 .
Figure 6.(a) Simulations from Newman et al. (2011) TitanWRF for the 2009 zonal wind velocities (b) Simulations from Newman et al. (2011) TitanWRF for the 2010 zonal wind velocities (c) Simulations from Lombardo & Lora (2023a) TAM for the 2009 zonal wind velocities (d) Simulations from Lombardo & Lora (2023a) TAM for the 2010 zonal wind velocities.Dashed horizontal lines represent the average peak contribution pressures for the eSMA wind measurements, and dotted horizontal lines represent the lower and upper bounds based on the associated error bars, indicating the range of uncertainty around the average peak contribution pressures.

Table 1 .
Observational Details Observation Time a D b Lat.c PWV d θ e Apparent Size f Beam Size g PA h B i A j Distance to Titan c Subobserver latitude: https://ssd.jpl.nasa.gov/horizons.cgid Precipitable water va- por e Angular separation between Saturn and Titan, as seen by eSMA f Solid body diameter (5150 km) g Natural beam shape h Position angle (PA) of natural beam i Maximum baseline j Physical collecting area

Table 2 .
Stratospheric Wind Speed Measurements * Order of magnitude