The raw RSS data of the two orthogonally polarized channels are sequences of complex samples made of in-phase and quadrature components (IQ samples) of the incoming signal sampled with a fixed frequency (16 kHz for this dataset)58. Data processing was performed using a basic averaged periodogram technique, looking for the best trade-off between spatial resolution on the surface, frequency resolution, and SNR. We carried out absolute amplitude calibration of the samples from system noise temperature models alongside an equalization to remove ripples introduced in the noise floor by the open-loop radio science receivers. A Gaussian function was fitted to the calibrated waveforms to measure the spectral broadening of the echoes at their full-width-half-maximum. The total reflected power was obtained in both polarizations by integrating the power spectral density over a frequency band spanning four times the measured signal spectral width45. The relative dielectric constant εr of the surface can be derived from the circular polarization ratio (CPR) of the same sense and opposite sense received powers through the Fresnel voltage reflection coefficients. The power ratio can be reduced to a ratio between the surface reflectivity in the two polarizations36, which depends on the relative dielectric constant of the reflecting surface and the incidence angle of observation (\(\theta\)).

$${CPR}=\frac{{P}_{R}}{{P}_{L}}=\frac{{\left|{R}_{R}\right|}^{2}}{{\left|{R}_{L}\right|}^{2}}$$

(1)

where

$${R}_{R}=\,\frac{{R}_{V}+{R}_{H}}{2\,}$$

(2)

and

$${R}_{L}=\,\frac{{R}_{V}-{R}_{H}}{2}$$

(3)

\({R}_{V}\) and \({R}_{H}\) are the Fresnel voltage reflection coefficients for horizontal and vertical polarizations, and they can be expressed as

$${R}_{H}=\frac{cos \theta -\sqrt{{\varepsilon }_{r}-{sin }^{2}\theta }}{cos \theta+\sqrt{{\varepsilon }_{r}-{sin }^{2}\theta }}$$

(4)

and

$${R}_{V}=\frac{{{{{{\mathrm{\varepsilon }}}}}\; }cos\theta -\sqrt{{\varepsilon }_{r}-{sin }^{2}\theta }}{{{{{{\mathrm{\varepsilon }}}}}}\,cos\theta+\sqrt{{\varepsilon }_{r}-{\sin }^{2}\theta }}$$

(5)

Inverting this whole system of equations results in

$${\varepsilon }_{r}=\left(\frac{{tan }^{2}\theta }{{CPR}}+1\right){sin }^{2}\theta$$

(6)

For a surface of a given relative dielectric constant, there is a unique incidence angle that yields the circular polarization ratio to unity. This is called the Brewster angle and can be easily retrieved as:

$${{{{{{\rm{\theta }}}}}}}_{{{{{{\rm{B}}}}}}}={{{{\mathrm{arctan}}}}}\left(\sqrt{{{{{{{\rm{\varepsilon }}}}}}}_{{{{{{\rm{r}}}}}}}}\right)$$

(7)

Cassini’s bistatic observations were usually performed at specular incidence angles θ ranging 50°−70°, thus slightly larger than the Brewster angle expected for Titan’s liquid surfaces (εr1.5–2 corresponds roughly to θB 50°−55°).

Theoretical studies35 found a relation between the broadening of the received echoes and the large-scale (up to 100 s of λ20, which is on the order of meters at X-band) surface roughness in terms of rms-slope for a surface with Gaussian statistics and s λ, where s is the root mean square (rms) variation in the surface height. Previous analyses of the Cassini RADAR altimeter and VIMS datasets acquired over Titan’s seas didn’t reveal any significant surface roughness at large-scale (10 s of km) but only at cm-scales or shorter. These observations were consistent with an upper limit on s of a few mm16,21,34,37,59, suggesting incredibly smooth surfaces with perhaps the presence of capillary waves at the wavelength-scale (λKu). As a result, the mm-scale roughness of Titan’s seas fails to satisfy the assumptions of previous models when compared to an X-band wavelength of 3.6 cm35. However, no assumption prevents the CPR of the received echoes from being used to estimate the relative dielectric constant of the surface. Once the effective surface dielectric constant is known, the bistatic scattering coefficient measured from either channel can be used to constrain the surface s at the wavelength scale for mainly coherent surface reflections39. While, in general, surface scattering consists of coherent and incoherent components, the flatness of Titan’s seas will cause a substantial suppression of the incoherent component. According to the Rayleigh Criterion s<λ/(8sinγ), with γ the depression angle, reflections from a surface observed at 50°−70° incidence angle can be considered predominantly coherent with s < 1.3 cm at X-band39. Unlike single-polarization monostatic RADAR measurements, where one cannot separate the effects of surface microscale roughness from composition26,37,60, with dual-polarization bistatic RSS data we are able to infer the effective values of both parameters. To investigate the wavelength scale of sea surface roughness, we use the power measured in each of two orthogonally polarized channels. More reliability is accorded to the results obtained using the RCP channel for flyby T101 and T102, due to their higher SNR values acquired at θ > θB. The results that we obtain through the use of the LCP channel give similar (where not identical) outcomes, but with larger errors.

Amplitude calibration

The computation of the true ratio of orthogonally polarized reflected powers of echoes from Titan’s surface to Earth played a pivotal role in our retrieval of correct values of both dielectric constant and small-scale surface roughness. We now discuss briefly the accuracy and possible sources of error in the estimation of the actual level of received power.

The values of the surface dielectric constant were derived from the true ratio of orthogonally polarized reflected powers. The power of Titan’s echoes can be computed from the difference between the total signal power received by the antenna and the noise power. When received on Earth at one of the DSN antennas, the complex samples at the output of the radio science receiver (RSR) are attenuated and/or amplified at different steps of the reception chain to increase the SNR and avoid saturation of the analog-to-digital converter (ADC). The exact amplifying gain, in general different for the two orthogonally polarized receiving channels, is unknown a priori, and needs to be estimated to scale the amplified powers down to their true values, which are on the order of 10−21 Watt (zW).

A typical RSS experiment performed at the closest approach to Titan lasts several hours and is made of three distinct parts: the BSR ingress observation, the radio occultation and the BSR egress observation. Before the ingress the amplifying gain is set, and it remains constant through the entire radio science observation.

The noise power recorded during the observation at the output of the RSR can be expressed as

$${N}_{{RSR}}=k\cdot {T}_{{sys}}\cdot B\cdot {G}_{{amp}}$$

(8)

where k is the Boltzmann constant, Tsys is the system noise temperature, B is the receiving bandwidth, and Gamp is the total amplifying gain through the reception chain. The variability of system noise temperature drives the noise power variations with time as long as the amplifying gain is stable. Since the system noise temperature mainly depends on the weather, the elevation angle of the receiving DSN antenna, and blackbody radiation from Saturn, it is possible to predict its behavior with time. Instabilities in the amplifying gain would appear as unexpected variations of noise power, not consistent with the previously mentioned phenomena. Such unexpected behavior is not present in the X-band datasets for the T101, T102, T106, and T124 flybys. Only the last of these observations features strong bumps of noise power, but they were recognized to be consistent with rainfalls detected at the DSN complex in Canberra during the observation of Punga Mare. Hence, we assumed the amplifying gain to be constant throughout each radio science observation.

Real-time measured profiles of system noise temperature are not available for the entire time of observation, but a calibration routine is scheduled before and after each flyby where a bistatic experiment is executed in order to provide enough information to model the noise history (we call these routines respectively Pre-cal and Post-cal). Before Bistatic Ingress and after Bistatic Egress, when the receiving DSN antenna is pointed at the zenith, the RSR connection is switched between the antenna itself, a well-known ambient load and a noise-generating diode nominally working at 12.5 K. A proper analysis of the noise jumps induced by the switches allows for the determination of the actual noise diode temperatures of the two orthogonally polarized channels before and after the entire closest-approach observation. At specific times (usually before Ingress, after Egress, and in-between the two) three short ‘minicals’ (about 5 minutes) are also scheduled to compute on-the-spot amplifying gain and system noise temperature. Noise diode temperatures serve as inputs to the minicals. For all the flybys we were able to determine diode temperatures TND for the two channels with less than 1.6 K of uncertainty (Table 2). High values of the RCP diode temperature, with respect to the nominal condition, are consistently observed in all four flybys, both during Pre-cal and Post-cal (NASA operators confirmed this issue in a report from flyby T119). This phenomenon should not represent a problem for our calibration algorithm, and although the RCP diode was working at a higher temperature compared to nominal, there is no reason to expect it to compromise the results of the BSR observations.

Table 2 DSS 43 (Canberra) noise diode temperatures

Table 3 shows values of inherent receiver temperatures for the four flybys of interest in this paper. These numbers are intermediate results of our calibration routine since the main output of pre-cal and post-cal analysis are the noise diodes’ temperatures, which are useful for the minicals analysis. We do not discuss all the details of the calibration process here, as we are planning to include them in a technical document on BSR calibration procedures that is currently in preparation. The main result shown in Table 3 is that the RCP channel is 4–6 K noisier than the LCP channel for all four flybys. The consistency between pre-cal and post-cal outputs, observed also in the noise diodes’ temperatures, was expected and strengthens our confidence in the results.

Table 3 DSS 43 (Canberra) inherent receiver temperatures

Phase calibration

Before and after a flyby observation, the Cassini orbiter was scheduled to perform 15 min of free-space baseline, transmitting an unmodulated carrier directly to Earth. As also observed during the Mars Express mission36, the received X-RCP carrier always comes with an apparent X-LCP leakage signal about 24 dB weaker. Such a small effect does not impact the computed values of the relative dielectric constant and thus was left uncorrected.

Equalization

The digitized transfer function of the receiving system is characterized by a spectrum that is not flat at the output of the analog-to-digital converters and produces ripples in the noise floor of the received signal. This may be a relevant issue when processing weak echo components whose amplitude could be on the order of the detected ripples (<2% of the noise pedestal). For completeness, the equalization was applied before computing the noise pedestal and the reflected power coming from Titan’s seas. A noise power spectrum was produced by integrating over tens of minutes before or after the BSR observation to accumulate an equalizing spectrum that was later smoothed by a moving average technique. The square root of this was used to scale the voltage spectra of the signal at the output of the RSR. Different equalizing spectra were derived for the two orthogonally polarized channels, but the same small ripple amplitude was detected.

Model of reflected power

For a proper retrieval of the received power, careful modeling must be carried out independently for the two polarization channels. After choosing the perfectly conducting sphere model and after applying an absolute calibration to the samples, potential power losses through the transmission chain, reception chain, and the medium should be addressed. We are aware that some minor physical phenomena may be missed, but we incorporated as many sources of loss as possible.

For example, we found that reception gain loss due to varying the elevation angle of the Canberra DSN antenna causes a slight change in the amount of received power61. Sky conditions at the DSN station can also change throughout the BSR observation and cause slight power losses that are difficult to model. However, qualitative station operators’ notes about the average weather during BSR observations and weather condition tables are publicly available. Ranging from a clear to a cloudy sky, we estimated that the difference in power loss would be only about 0.05 dB, which is negligible61. We assumed the average clear sky for all the observations and using the formula included in ref. 61 for atmospheric attenuation due to the weather conditions recorded during data reception, we found that the values retrieved are usually <0.5 dB. We applied this correction to the data for the estimation of surface roughness. Note that the T124 observation of Punga was affected by rain at Canberra during reception of the data (≤35 mm/h, ≤75% relative humidity), while during the subsequent observation of the main body of Kraken Mare, acquired about 45 minutes later, the weather was clear (about 0 mm/h and 45% relative humidity). Despite this, we have successfully modeled the increase in noise temperature due to the attenuation from a wet atmosphere (see the comparison between noise level in echoes from Punga and from the other seas in Fig. 3).

Gain losses due to wind loading on the DSN antenna are negligible. We estimated the pointing error of the Cassini antenna with respect to the center of the specular reflection to be <0.1 dB. For this calculation, we used the Ku-band antenna patterns. Because the X-band antenna −3 dB aperture is larger than the Ku-band aperture, we conclude that the pointing error is negligible. We estimated that the pointing accuracy of the Canberra DSN antenna at X-band (0.032 deg) produced losses <0.01 dB and is, therefore, also negligible.

The total power received on the ground from Titan’s specular reflections can be modeled using the basic bistatic radar equation20,35:

$${P}_{R}=\frac{{P}_{T}{G}_{T}}{4\pi {\left|T\right|}^{2}}\sigma \frac{{G}_{R}{\lambda }^{2}}{{\left(4\pi \left|R\right|\right)}^{2}}$$

(9)

where PT is the power transmitted by the spacecraft, GT and GR are, respectively the gains of the transmitting antenna aboard Cassini and the receiving antenna on Earth (Canberra DSS in our case). These quantities were calibrated by exploiting the 15-minute-long free-space baselines scheduled before and after each BSR experiment, when the spacecraft was transmitting an unmodulated signal directly to Earth. Referring again to Eq. (9), T and R are the distances of transmitter and receiver from the center of Titan, and σ is the bistatic radar cross-section of the observation target, which in this case is the rough surface of a planet. Fjeldbo derived a simple model for the radar cross-section of a rough surface under the Kirchhoff approximations (KA)35. If we consider reflectivity to be the ratio between true received power and reflected power expected from a perfectly conducting sphere20,45, then the radar cross-section can be expressed as:

$$\sigma=\frac{4\pi {\left|T\right|}^{2}{{R}_{P}}^{2}cos \theta }{({R}_{P}cos \theta+2\left|T{\prime} \right|)({R}_{P}+2\left|T{\prime} \right|cos \theta )}*{\varGamma }_{r}$$

(10)

where RP is the planet radius, θ is the incidence angle, T’ is the distance between the transmitter (Cassini) and the specular point, and Γr is the reflectivity affected by wavelength-scale surface roughness. The reflectivity, which describes how power is redistributed among two orthogonal circular polarization senses after reflection, depends on incidence angle and relative dielectric constant through the Fresnel reflection coefficients36. For a random isotropic surface, roughness has the same effect on reflectivity in both polarizations39, and is modeled with a scale factor function of s:

$${\varGamma }_{r}=\varGamma*\exp \left\{-4{\left(\frac{2\pi }{\lambda }{{{{{\rm{s}}}}}}cos \theta \right)}^{2}\right\}$$

(11)

where Γ is the reflectivity from a perfectly smooth surface. s can be estimated by inverting the system of formulas reported above for every single received power measurement PR acquired within each area of interest and using the relative dielectric constant εr obtained from the CPR. From these estimates, we can identify a mean value for s and relative 1σ errors.