Observation of 2nd Schumann eigenmode on Titan’s surface

. This work presents the results obtained from an updated data analysis of the observations of extremely low frequency (ELF) electromagnetic waves performed with the HASI-PWA (Huygens Atmospheric Structure and Permit-tivity, Wave and Altimetry) instrumentation after Huygens Probe landing on Titan’s surface in January 2005. The most signiﬁcant signals observed at around 36 Hz throughout the descent in the atmosphere have been extensively analyzed for several years, and subsequently interpreted as the signature


Introduction
One of the main objectives of the PWA-ELF experiment was to measure the ELF spectral density of natural electromagnetic waves received by an electric dipole during the descent of the Huygens Probe throughout Titan's atmosphere (Grard et al., 1995).A clearly identified signal at around 36 Hz was received for about 2 h 25 min of descent and furthermore interpreted as the second eigenmode of an atypical Schumann resonance (Béghin et al., 2007(Béghin et al., , 2009;;Simões et al., 2007).The main differences between Earth's and Titan's Schumann resonances (hence referred to as SR) are summarized in Table 1 according to most recent studies (Béghin et al., 2012).Other planets such as Venus, Jupiter and Saturn, had been proposed well before the Huygens-Cassini mission by Nickolaenko and Rabinowicz (1982) as potential candidates to experience lightning triggered Earth-like SRs.However, Titan revealed itself as a unique case, because the only source of available energy could presumably result from the Saturn's magnetosphere interaction with Titan's atmosphere (Béghin et al., 2007), a view furthermore supported by the established absence of any lightning in Titan's atmosphere (Fischer and Gurnett, 2011).The present study is a follow-up to the work of Béghin et al. (2012), in which the persistency of the 36 Hz signal is reported for about 4 to 6 s after the landing of the Huygens Probe on the surface of Titan (Lebreton et al., 2005).The signal seemed to disappear or to be hidden in the instrumental noise once the gondola was definitely at rest.Nevertheless, a possibility to detect it appeared after a comparative study performed on the statistical distributions Published by Copernicus Publications on behalf of the European Geosciences Union.
of the ELF signal amplitudes recorded during the three following sequences: i. the instrument check-out #10 performed during the interplanetary cruise, when the electric dipole antenna was stowed and shielded under the thermal cover (Béghin et al., 2009); ii. the first part of the descent, between 140-110 km altitude, under the large parachute and before its jettison at Mission Time (MT) 900 s, when the antenna booms were supposedly not totally deployed (Hamelin et al., 2007); and iii. the first 32 min on the surface, after the Huygens landing, up to the pre-programmed switch-off of the HASI-PWA instrument.
The analysis of the three sequences has been performed on the data files (Ref.PWA-Data-Base, 2013) of the DFT applied onboard to the waveform of ELF signals collected by the double-probe boom-antenna (Grard et al., 1995).We then compare the bin at around 36 Hz (hence referred to as bin 36) with the other bins of the spectra.We demonstrate in Sects. 2 and 3 that the peculiarities of the amplitude spectral density (ASD) distribution of bin 36, recorded on the surface, clearly reveal the presence of a coherent natural signal that is otherwise barely visible in the averaged spectral density distribution.These peculiarities are ascribed to the data processing methods applied onboard, namely the logarithmic compression and the discrete quantification of the DFT spectral line amplitudes.Two different numerical simulations (Sect.4) of the onboard data processing loop are proposed not only to support the SR detection on Titan's surface, but also to estimate its amplitude from signal-to-noise ratio (hence referred to as SNR).

Amplitude spectral distribution of PWA-ELF data
The mathematical treatments of the wave field data were designed at least ten years before the Huygens landing so as to be performed automatically onboard (Hofe, 2005).The ELF power spectral density (PSD) was computed by applying a DFT after a 16-bit analog-digital conversion of two consecutive waveform samples of 333 ms duration each.The square root of the composite DFT modulus (i.e., the ASD) was logarithmically compressed, and the lower byte (8 bits) was transmitted without the phases to Earth by the telemetry system (hence referred to as TM) via the Cassini orbiter (Lebreton et al., 2005).Therefore, this study is constrained by the limited performances of such usual techniques of wavedata processing of space experiments; nevertheless, we take the best advantage of the proper peculiarities of the PWA-ELF instrument, as briefly described here below.The experiment has been operated in two pre-programmed modes, labeled 131 and 132, dedicated respectively to the high and low altitude ranges (Jernej and Falkner, 2004).In mode 131, between 140 and 60 km, each spectrum consists of 32 frequency bins, centred on frequencies evenly distributed between 3 and 99 Hz with a 3 Hz bandwidth.In this mode, bin 36 is actually centred at 36 Hz and covers the frequency range of 34.5 to 37.5 Hz.In mode 132, below 60 km, the computed spectra consist of 16 bins, centred on frequencies evenly distributed between 6 and 96 Hz with a 6 Hz resolution.More precisely, the second mode is derived from the first one by adding the power density contained in two adjacent bins of the first mode, so that the central frequencies in mode 132 are shifted downwards by 1.5 Hz with respect to those of the even bins in mode 131.For instance, bin 36 covers a bandwidth of 6 Hz, from 31.5 to 37.5 Hz within 3 dB amplitude range (Hofe, 2005).We are then unable to resolve any peculiar frequency within that range.Since the aim of this work is rather to identify the presence of a natural signal in bin 36, for the purpose of simplicity we will ignore the 1.5 Hz leftward shift of the frequency scale that concerns the data in mode 132.Hence, a nominal frequency of 36 Hz will be arbitrarily given to bin 36, and similarly to the other bins.
Unfortunately, due to the regrettable loss of one of the two TM links (channel A), one half of the original bins were lost, so that only the eight harmonics of bin12 in mode 132 were recovered from the surface data set.
The data are decompressed on the ground according to the transfer function (Jernej and Falkner, 2004), which is valid for both modes, where V ADC is the peak amplitude of an equivalent sine-wave signal, with the given frequency injected at the input of the analog digital converter (ADC) and I TM is an integer that represents the lower 8 bits of the log compressed ASD of the average of the two consecutive temporal samples.A portion of the transfer function is shown in Fig. 1.Since we consider composite signals, the DFT power density of which is assumed to fill up the entire bin frequency bandwidth, the root mean square (RMS) value remains V ADC / √ 2 as for sine-waves.The ASDs denoted respectively 1 and 2 in either mode, and the amplitude spectral density of the electric field component received by the antenna in Vm −1 Hz −1/2 unit, are given respectively by the following set of equations,    SR component is assumed to lie (Béghin et al., 2012).Since the precise position of both electric sensors with respect to the local ground is still under investigation, we shall keep in mind for the moment the ratio between the antenna voltages before and after touch down.
It is worth emphasizing here that owing to log compression and the 8-bit quantification, the dynamic range of V ADC is far from linear (Fig. 1).This is the main point that we take advantage of in this work for extracting the weak SR signal from the noise.Even in case of initial normal (Gaussian) distribution, such a non linear process implies that the output decompressed values should exhibit a notably different amplitude distribution.A bias of the amplitude distribution is indeed visible on almost all ELF-PWA spectral data during the descent (Béghin et al., 2009).
In order to prove that this bias is purely experimental and due to the peculiarities of the V ADC versus I TM transfer function, we first consider in Fig. 2 the statistical characteristics of the instrument noise measured during the Cruise checkout # 10, labelled above as sequence (i) and performed in the complete absence of natural signal.The bias of a finite series of discrete samples can be defined in different ways (e.g., Ghahramani, 2000), but we shall consider here that a distribution is biased as long as there is a significant difference between the numbers of samples distributed on both sides of the mean value µ of the series.Note that the mean value becomes the expectation whenever all samples have the same likelihood, as for instance a noise with a normal amplitude distribution.We shall consider also the quartiles Q 1 , Q 2 (or median) and Q 3 , defined as discrete values of samples splitting respectively the lowest 25, 50 and 75 % fractions of a series.The bias of the ELF-PWA data is clearly revealed in Fig. 2 (left panel) by the dissymmetry between the Q 1 and Q 3 amplitude bars, that implies a shift between the mean µ and the median Q 2 .We say that the bias is positive whenever the probability in the vicinity of Q 1 is larger than that in the vicinity of Q 3 .We observe also that Q 2 is shifted leftwards with respect to µ.The amplitude vs. time plotted in Fig. 2 (right panel) confirms indeed that most samples lower than Q 2 are focused in the vicinity of Q 1 .This comes from a substantial excess of moderate and weak amplitude samples which is required to balance the weight of higher values, as a result of the non-linearity of the transfer function (Fig. 1).
The statistical properties of the data collected during the Cruise checkout #10 will serve as a reference for comparison with situations when weak natural signals are present, such as sequence (ii) during the early phase of the descent and sequence (iii) after landing.These two sequences contain respectively 246 and 418 samples of 16 and 8 frequency spectra each.The plots in Fig. 3 represent the frequency distributions of the mean V ADC RMS amplitudes (µ f ), plus/minus one standard deviation (σ f ) for these two sequences.One can check that the values of µ f in both modes, 131 and 132, comply with Eq. (2) for a wideband noise.The V ADC RMS values in mode 132 (right panel) are indeed twice those of mode 131 (left panel), except for bin 36, when there is an obvious contribution from the SR.Moreover, the noise frequency distribution during the Cruise checkout has the same shape in all frequency bands, but bin 36, throughout the Huygens descent, except below 20 km, when an additional contribution is possibly caused by the passage through an atmospheric haze layer (Béghin et al., 2007).We shall therefore consider this noise spectrum shape as an intrinsic characteristic of the instrument.In addition to the definitions of the mean µ and of the three quartiles, we now recall those of other symbols that will be considered below in order to avoid any confusion with the terminology sometimes used in the literature (e.g., Ghahramani, 2000).For large number of samples (N > 100) at a given frequency (f ), the variance (s 2 f ), and the standard variation (σ f ), are given by where X i is the amplitude of an individual sample of index i, and µ N the mean amplitude value of the series.X i represents either the V ADC voltage amplitude or its RMS value.
When the number of frequency samples n f is less than or equal to 16, such as in the plots of the mean values versus the frequency (Fig. 3), we use the linear flicker regression (LFR) method to fit the most probable analytical function representing the spectral distribution of noise in semi-conductor devices (e.g., Marshall Leach Jr., 1994, and references therein).Although the classes of different kinds of flicker noises are "as ubiquitous as they are mysterious" after Milotti, (1995), the PSD analytical shape exhibits usually a 1/f β dependence, with β growing from 0 to 2 with increasing frequency.A transition between β = 1.5 and 2 occurs for frequencies such as f > 1/2π τ , where τ is the characteristic time constant of temporal samples.Moreover, note that in our case of ASD data, the shape would be 1/f α with α = β/2, so that the LFR function is expressed by where A and α are the two coefficients deduced from the mean square regression analysis of a two columns matrix made of n f frequency bins and their associated mean amplitudes x i .The quantity δ is the mean absolute deviation between the actual values and the mean flicker noise µ f .
During sequence (ii), under the large parachute (Fig. 3, left panel), the presence of the signal in bin 36, associated with its side band contribution in bin 30, is well visible above the mean-fit flicker noise with a maximum value of 24 mV RMS.On the contrary, during sequence (iii) on the surface (Fig. 3, right panel), the presence of the signal is barely discerned.The LFR coefficients for the flicker noise are found to be A = 0.66 V Hz α and α = 0.87, with δ = 0.9 mV.Reporting these coefficients in Eq. ( 4) leads to a mean fitted value µ f = 29.2mV RMS (i.e., I TM = 84) at around 36 Hz (Fig. 3, right panel, blue solid line), whereas the actual mean among the 418 data samples of bin 36 is µ 36 = 30.4mV (still ).An increase of 1.2 mV, slightly more than one δ above the mean LFR instrumental noise, remains however a marginal evidence for the presence of a natural signal in bin 36.We will nevertheless see in the next Sections that this value of 30.4 mV lies just below the jump from I TM = 84 to 85, which is a major indicator for the presence of a signal.
We have plotted in Fig. 3 (right panel) the smallest values of each series of 418 samples for all frequency bins of the surface sequence.The LFR fit profile of the smallest values versus the frequency (continuous red plot) is derived from all bins, but the 36 Hz one, in order to avoid biasing the statistics with the possible contribution of an additional signal to the noise.One may indeed reasonably assume that the flicker noise contribution at around 36 Hz should not change significantly the global LFR profile derived from other bins.The actual smallest amplitude of bin 36, in terms of telemetry step is I TM = 72 (i.e., 11.95 mV RMS), whereas the LFR fit in the absence of any signal yields I TM = 69 (9.63 mV).Such increase of 2.32 mV RMS above the LFR profile (Fig. 3, right panel, solid red line), although more significant than that of the mean noise level, does not allow us yet to rule out completely the signature of a purely random event.We now investigate the bias more precisely denoted by the skew of the bin amplitude distribution.

Skew versus SNR
The skew during sequence (ii), before MT 900 s (Fig. 4, left panel), is similar to that observed during the entire descent, and also during the checkout sequence (Fig. 2) used as a reference.Namely, the median Q 2 is shifted leftwards with respect to the mean value µ, and the amplitude distribution bar in the vicinity of the quartile Q 1 is higher than that around Q 3 , which means a positive skew.However, the situation is opposite during sequence (iii), on the surface (Fig. 4, right panel), where Q 2 is shifted rightwards and the higher bar is approximately centred on Q 3 .In addition to singularities of bin 36 discussed in the previous Section, it should be emphasized that this skew reversal is a remarkable feature that occurs only in the surface data and requires further investigation.The conventional definition of skewness for a biased distribution is given by (5) Applying Eq. 5 to the ELF data throughout the descent yields mean values of G 1 of about 0.45, always positive at all frequencies including bin 36 (Fig. 9c in Béghin et al., 2009), whereas the same quantity lies in the range 0.16-0.424during the surface sequence, with an intermediate value of 0.244 for bin 36.A similar situation holds during the cruise checkout sequence # 10, where the values are distributed between 0.25 and 0.58, with G 1 = 0.36 for bin 36.We therefore consider that the conventional definition of skewness does not explain the singularity of bin 36 on Titan's surface.On the other hand, the shift between Q 2 (a discrete value associated to an integer I TM ) and µ (a continuous variable) is not a satisfying quantitative evaluation of the skew.Indeed, whatever N odd or even, the skew with respect to Q 2 should depend on the arbitrary choice whether the individual samples of a bin series are considered strictly equal or smaller than Q 2 .
We rather prefer an alternative definition for the skew applied to the mean value µ which may lie in the vicinity of -but rarely strictly equals-an I TM step.Whenever the I TM digit increases or decreases by one unit, the quantity |X iµ f | exhibits a sudden step in accordance with Eq. ( 1).As a consequence, the probability distribution of samples lying in the vicinity of µ f depends on its relative position with respect to the two bracketing I TM integers.We introduce then the skewed standard deviation σ f , an estimator for each frequency bin, defined as where S f is denoted normalized skew for any bin series with the frequency f , and K µf is the number of samples of the series with an amplitude larger than or equal to µ f .For a normal distribution, K µf ≡ N/2 if N is even, and K µf ≡ (N+1)/2 if N is odd, so that S f = 1 in either case.The distribution is said to be skewed as long as there is an abnormal excess or deficiency of samples above µ f .According to the definition given in the previous Sections, the bias is positive when the skew S f is larger than 1, as observed for most of ELF-PWA data, either due to a pure instrumental noise (e.g., cruise check-out sequence), or when SNR is larger than 1 (Fig. 4 left panel).It will be seen in the next Section that the reverse situation (S 36 < 1) is encountered, when the signal amplitude of bin 36 is slightly larger than the mean instrumental noise, but less than the step I TM = 85.A sudden jump through the value S 36 ≡ 1 occurs indeed whenever µ 36 is nearly equal to -but strictly less than -the amplitude corresponding to the step I TM 85.Since the value of the skewed standard deviation is thought to be due to the peculiarities of the PWA instrument, namely, log compression and 8-bits TM transmission processes, there is no physical justification to  assume that its own frequency distribution obeys the flicker law.We will then apply here the general linear polynomial regression LPR function where q is an integer usually smaller than 3 and the coefficients a and b q are deduced from the least square regression analysis of a two-column matrix, as done above with the LFR fit.
Though the mean amplitude µ 36 of bin 36 during sequence (ii) (before MT 900 s) is well visible above the noise level (Fig. 3, left panel), the plot of the skewed standard deviation σ 36 (Fig. 5 left panel) exhibits an even larger jump above the LPR profile.The situation is reversed for the surface data because S 36 is smaller than 1, then σ 36 < σ 36 .The best LPR fit applied to the value of σ f of all bins, but bin 36, which optimizes the mean absolute deviation δ = 0.061 mV, is obtained with the following parameters: q = 3, a = −1.522mV, b 1 = 510.98mV Hz, b 2 = −5681.1 mV Hz 2 , b 3 = 33 450.8 mV Hz 3 .Introducing these values in Eq. ( 7) for f = 36 Hz, the fit yields σ 36 = 9 mV RMS, whereas the experimental value derived from Eq. ( 6) is equal to 8.02 mV RMS: about 16 δ below the LPR noise measurements in the absence of signal (Fig. 5 right panel, dotted line).
Summarizing the above survey, we retain that the probable presence of a natural signal on Titan's surface should be identified by three indicators observed only with bin 36 Hz, that are respectively by order of significance: (i) a reversal of skew while the major part of the sample distribution lies in the vicinity of a step of the transfer function, (ii) an excess of smallest values compared to the distribution of other bins, (iii) a weak although noticeable increase of the nominal average amplitude.We therefore propose to reproduce such a behavior and to confirm the fact that the involved mechanism is due to the peculiarities of the onboard data processing and TM transmission.Two different numerical simulations of the entire loop were performed, starting from the ADC input, to the DFT process up to the ground data decommutation.

Numerical simulations
The common purpose of the two simulations is to assess the probability of presence of a natural signal in bin 36, and eventually, to estimate the order of magnitude of SNR.In each simulation we are referring to the same experimental data, with the parameters reported in Table 2.The procedure is basically the same as performed on board Huygens for the mode 132 and we make the following assumptions: i. in the absence of signal, the instrumental noise spectrum obeys the LFR function (Eq.4), with the parameters derived in Sect. 2 (A = 0.66 V Hz α and α = 0.87); ii. the mean spectral characteristics of both signal and noise are assumed to be stationary during the ground sequence; iii. the amplitude of each V ADC sample (either Volt, mV or RMS, whenever applicable) is the product of its amplitude spectral density 2 by the square root of the bin resolution (Eqs. 1 and 2); iv. the composite RMS amplitude A 3 of both waveforms (noise A 1 plus signal A 2 ) at the DFT output reads where ϕ 12 is the differential phase between signal and noise components within the DFT complex plane;  vii. the direct transfer function (Eq. 1) is then applied to every I TM sample of each set; and viii.the global results of two successive series of 100 runs each (simulation S1) or 10 000 runs (simulation S2), are reported below, with the relevant parameters summarized in Table 2.
Note that due to the averaging process of two timeindependent individual samples (step v), each combined output I TM value is different from that of the original input values.This is the reason why, after completion of the simulation processes, the values "out" of µ and σ of the background noise, in the absence of signal, may differ from the values "in" (Table 2).Such feature accounts for the fact that we are not allowed considering the transfer function Eq. ( 1) as a biunivocal relation applied through the entire ELF-PWA loop.

Simulation S1
A first phase of 100 runs has been performed by scanning a wide range of SNR from 0 to 0.6.Three different files of noise amplitudes (A 1 ), 418 samples each, have been computed as images of the actual surface noise data bins 12, 24 and 48, after applying the following normalized flicker noise coefficient to the amplitude of each sample where A 1f is the noise amplitude of the nth sample of the relevant data file at f = 12, 24 or 48 Hz, and α = 0.87 (Sect.2).
We have checked that in the absence of signal (A 2 = 0), the simulation yields the theoretical composite value of LFR fit: µ 36 ∼ 41.3 mV (29.2 mV RMS).In the following, since we consider composite signals assumed to fill in the entire 6 Hz bandwidth of mode 132, it is more convenient to make use of RMS voltages, unless stated otherwise.This first phase comprised a total number of 100 runs using consecutively the noise image files at 12, 24 and 48 Hz, with 418 samples each.The mean signal amplitude and standard deviation varied in steps between two successive runs so as to scan the SNR from 0 to 60 % given by the ratio A 2 /A 1 (Fig. 6, second panel from top).One can see that this range covers the composite mean values µ distributed between 29.2 and 32.9 mV, which corresponds to I TM comprised between 84 and 86 (Fig. 6, lower panel).The plot of the indicator (i) i.e., the normalized skew S 36 versus µ (same panel) is a clear demonstration of the bias effect induced by data discretization.This effect occurs when the mean amplitude of the sample series lies in the vicinity of a step of the transfer function, which is actually I TM = 85 for bin 36 on the surface.The simulation retrieves the predicted jump from S 36 < 1 to > 1 at around µ = 30.45mV RMS (i.e., 43.06 mV V ADC ) which corresponds to a jump from I TM 84 to 85 according to Eq. ( 1) and Fig. 1, yielding the amplitude dissymmetry observed between Q 1 and Q 3 bars (Fig. 4, right panel).The surface sequence contains a significant although limited number (418) of experimental data samples, from which we got the values of µ 36 = 30.38mV and σ 36 = 8.7 mV (Table 2).The simulations aimed to retrieve values of the same order of magnitude as actually measured, with some margin of uncertainty.In order to estimate this margin we consider that the mean statistical parameters of both signal and noise were staying stationary during the ground sequence, as usually assumed with terrestrial SRs processing.In such case, the deviations of these parameters should not undergo significant change for a smaller number of samples.We checked indeed that, irrespective of the length of series considered, either the first or the last 200 samples of the surface data, the experimental uncertainties are such that µ = 30.3± 0.2 mV and σ = 8.65 ± 0.25 mV (Table 2), and the skew S 36 = 0.918 ± 0.004.Special attention is paid to the highlighted areas in Fig. 6 just below the step I TM 85 where we retrieve satisfactorily the expected values for µ, σ and S 36 .The self-consistency of these parameters within the spread range of expected values enables us to find out that the most probable range of SNR lies between 24 and 36 % (Fig. 6, middle panel) yielding an expected signal amplitude of 9.05 mV (Table 2) for SNR about 0.3 (Table 2), during the 32 min on Titan's surface.
A second phase of 100 runs was performed by scanning the most likely range of SNR between 0.22 and 0.36: corresponding to experimental values of µ 36 bounded between 30.1 and 30.45 mV (blue filled circle and yellow line in Fig. 6, middle and lower panels respectively).The aim was to estimate the probability that the values predicted by simulation may fit the indicators identified in Sects. 2 and 3. We wanted first to check the anomalous distribution (Fig. 3, right panel) of the smallest I TM value of bin 36 with respect to the other bins i.e., the indicator (ii).We saw in Sect. 2 that the smallest amplitude of bin 36 among the 418 samples of data corresponds to I TM = 72, whereas the predicted LFR fit without signal should be I TM = 69.We have plotted in Fig. 7 the probability distribution of the lowest I TM value obtained from the 2 × 50 runs of 418 samples each, with noise images of bins 12 and 24.The shape of the distribution is obviously far from that of a pure random process.We find a 54 % chance that the lowest I TM sample be larger than 69 (i.e., presence of signal), with a maximum probability at I TM 71 very close to I TM 72 observed with the same number of samples (N = 418) in PWA surface data.Moreover, when we split the surface sequence into four consecutive sets of 100 samples each, the lowest values is again I TM 72 for sets # 1, 2, 4 and I TM 74 for set # 3, which suggests that this parameter might depend on the number of data samples and of their ranking as well.Nevertheless, on the basis of the probability distribution of this indicator alone, we may claim from the plot in Fig. 7 that there is a 54 % chance that the 2nd SR  eigenmode was observed on Titan's surface from the Huygens Probe touch down up to 32 min later.
The second phase of 100 runs allowed also to check another aspect of the indicator (i) which occurs whenever the skew S f is smaller than 1 (Fig. 6, lower panel): when the integer number K µ is larger than 209 for N = 418 (Eq.6).We have plotted in Fig. 8 (left panel) the distribution of K µ as a function of SNR in the range from 0.22 to 0.36.All points are lying well above K µ = 209 because these runs were intentionally selected within the experimental data range of µ from 30.1 up to 30.45 mV (Table 2, columns 4 & 5).It appears that the experimental value K µ = 227 (Fig. 8 left panel), which yields S f ∼ 0.92 (Fig. 6, lower panel, yellow line) might have been obtained for any value of SNR lying in this range.One may however anticipate that for a larger number of runs, the highest probability would rather lie at around K µ ∼ 222 (i.e., S f ∼ 0.94) instead of 0.92, which confirms that such indicators depend upon the number of data samples.But, the probability for K µ = 209 (no signal) is extremely low (1 % in Fig. 8, left panel).
The last result derived from the second series of 100 runs concerns the indicator (iii) namely, the distribution of the mean composite value µ versus SNR (Fig. 8, right panel).
The central position (yellow line) of the expected experimental value (µ ∼ 30.3 mV) corresponds to about a 50 % probability for the SNR lying between 0.24 and 0.36, which encompasses the highlighted range in Fig. 6 and besides emphasizes the self-consistence between the three indicators identified in Sects. 2 and 3.

Simulation S2
The main purpose of this simulation is to evaluate to which extent we may be confident about the statistical parameters deduced from a limited number of experimental data samples, namely N = 418 in our case.By extending the simulation up to 10 000 runs, we consider different approaches compared to simulation S1 for setting up both signal and noise files.Nevertheless, the procedure follows basically the steps listed in Sect.4, with some differences with respect to the simulation S1:  i. the noise files with 418 samples each are computed independently from each other by using the statistical parameters (µ N and σ N ) of a random distribution obeying the LFR fit function (Eq.4) with coefficients A and α derived in Sect.2; ii. each composite sample of any series of 418 is derived from Eq. ( 8) with a constant signal amplitude A 2 corresponding to a given value of SNR with a random phase ϕ 12 where the noise level is the expectation of the global distribution plotted in Fig. 9 (left panel); iii. the only variable quantity for each series of 418 samples is the value of SNR which is scanned from 0 to 0.6 in a large number of runs (10 000); and iv. for each SNR value, two successive sets of 418 composite samples are averaged in pairs, log-compressed and TM converted, as above.
Each file of 418 samples with different background noise amplitudes, introduced in 10 000 different couples, are computed according to a normal random distribution obeying the LFR function (Eq.4) for f = 36 Hz, with the expectation µ N = 30 mV (RMS) and σ N = 8.7 mV (RMS) at the input of the loop.Because Eq. ( 1) is applied successively in reverse and direct ways, according to the ending remark in Sect.4, the output values are slightly smaller (i.e., µ N = 29 and σ N = 8.3 mV).Therefore, both sets of values are denoted "in" and "out" respectively in Table 2.Note that the same effect holds for the simulation S1.From the distributions of µ N and σ N samples plotted in Fig. 9 one deduces that the predicted experimental noise (values "out") on Titan's surface in the absence of natural signal should be characterized by the following parameters µ N = 29 ± 1.5 mV RMS and σ N = 8.3 ± 1 mV RMS.(10) After performing steps (ii) to (iv), the results are summarized by statistical plots of the distributions of the mean amplitude of composite signal (µ) and of its standard deviation (σ )in Figs. 10 and 11, respectively.These distributions allow us to confirm that the most likely range of SNR lies between about 0.2 and 0.35, in good agreement with that derived from simulation S1 (Figs. 6 and 8), although both approaches are using different background noise figures.The new information, however, is that the most likely SNR value should be either 0.35 according to the value of the signal composite, or 0.2 according to the standard deviation.We interpret such discrepancy as due to much larger uncertainties on the standard deviation σ n of background noise than on µ n .At the input of simulation S2, we assume σ n = 8.7 mV RMS, with an overall dispersion from 7.5 to 9.75 mV (Fig. 9 right panel and Table 2).As a consequence, the standard deviation of the output composite noise exhibits a total spread of 8-9.3 mV (Fig. 11 and Table 2), whereas the presence of the signal reduces significantly this spread, thanks to the addition of many values around the composite mean.Therefore, we should rather trust the distribution of the composite mean plotted in Fig. 10, and consider that the most likely value of SNR lies at about 30 %, which reconciles this result with the estimate derived from simulation S1.

Discussion and conclusion
Although both simulations are using different hypotheses for the statistical properties of the noise and different number of runs (200 versus 10 000), we must note that they both lead to similar findings, namely: first of all, the option absence of natural signal in bin 36 (SNR = 0) is strictly inconsistent with the summary charts plotted in Figs. 6, 10 and 11.Second, from the self-consistence of specific indicators and statistical peculiarities of that bin, such as mean and standard deviation, skew and smallest samples of its amplitude distribution, we conclude that there is more than a 50 % chance that the estimated value of SNR lies between 0.24 and 0.36, Geosci.Instrum.Method.Data Syst., 2, 237-248, 2013 www.geosci-instrum-method-data-syst.net/2/237/2013/ with a most likely value of 0.3 (i.e., a signal of about 9.05 mV RMS within 6 Hz bandwidth at the input of the ADC).Introducing in Eq. ( 2) the equivalent amplitude V ADC = 12.7 mV, the nominal effective length of the antenna (l eff = 1.6 m) and G(36) = 13.5, that is, 22.6 dB volt at 34.5 Hz (Jernej and Falkner, 2004), we obtain an induced electric field strength of about 0.12 mV m −1 Hz −1/2 .As a useful comparison, this is 10 times less than the last measurements performed during the last phase of the descent and 4 to 6 seconds after touch down (Béghin et al., 2012).However, since the main component of the conventional SR modes is known to be vertical, the estimate of the actual strength of the incident wave-field vector is most questionable, as it depends on the Huygens motion and tilt of the boom with respect to the local vertical, and whether the sensors are free or in direct electric contact with a lossy dielectric ground (Grard et al., 2006;Béghin et al., 2012).
According to a recent work (Schröber et al., 2012), for about 3 s after touch down, Huygens gondola successively bounced back out of the hole impact, slid and wobbled back and forth five times, after which it commenced a 30-40 cm long slide on a flat surface for about 2 s.During an additional 5 s, the slide motion progressively slowed down until Huygens stayed definitely at rest.Because of the ELF-PWA electronic saturation visible on the first two spectra transmitted 2 and 4 s after impact (MT 0 = 8870 s), the first validated data are received at MT = 8878.375s.Since there is a processingbuffering delay of about 2 to 4 s between the acquisition time of the electric field and the MT dating, we may only assume that the first available measurement after impact was performed during the first 10 s after touch down.On the other hand, we are essentially concerned with the value of the tilt of Huygens Y axis which coincides with the nominal alignment of the PWA boom antenna (Grard et al., 1995).We shall then refer to the measurements performed by the Y-tilt sensor of the Surface Science Package (SSP).After correction of a permanent minus 8 • offset of this sensor (Leese et al., 2012), we shall consider an average tilt of about 2 • between the nominal attitude of our antenna and the local surface during the first 10 s after impact, instead of 10 • initially reported.A tilt of 2 • is consistent with the fact that this value yields the same strength for induced electric field (∼ 1.2 mV m −1 Hz −1/2 ) as that observed several minutes before impact, as well as a few seconds after touch down (Figs. 2 and 4, in Béghin et al., 2012).
Consequently, our estimate of a signal amplitude ten times smaller during the whole surface sequence after the first 10 s should reasonably lead to a tilt also ten times smaller.However, this is not consistent with the observations of other SSP instruments leading to claim that Huygens was definitely stabilized about 10 s after impact (Schröber et al., 2012).For instance, the value measured by the Y-tilt sensor seems to indicate still a permanent tilt of 2 • after offset correction.On the other hand, the mutual impedance (MI) device, a component of the HASI-PWA instrument, designed to measure the ground conductivity and using partly the same sensors as the ELF dipole antenna (Grard et al., 2006), experienced a regular decay in sensitivity from 10 s after impact before reaching an average stable value (i.e., a behavior compatible with that of the ELF bin 36 data).Such coincidence, still under investigation, could perhaps be caused by some motion of at least one of the boom-antenna without any perceptible influence on the Y-tilt sensor.Alternatively, according to further studies in progress, a slow change in conductivity of the near surface, due to the presence of the massive gondola warmer than the dusty surface sediment (Schröber et al., 2012), might explain the simultaneous change of the MI measurements and of the pattern of the local SR wave field.
We have, nevertheless, reached the main objective of the present work which was to identify and quantify the indicators which reveal the presence on Titan's surface of the 2nd SR eigenmode observed during the Huygens descent throughout the atmosphere.We found that the statistical characteristics of the bin 36 spectrum data are inconsistent with the normal flicker noise pattern of the other bins, which implies that we must reject the option SNR = 0 for that bin.Therefore, we estimate that there is more than a 50 % chance that a signal within the frequency range 34.5 ± 3 Hz be permanently present on the surface, with a SNR of 0.3 (i.e., a mean electric field induced in the dipole antenna of about 0.12 mV m −1 Hz −1/2 ).This value is 10 times weaker than those observed several minutes before Huygens impact and 4 to 6 s after landing.Further work is foreseen, using still progressing investigations on both Titan's surface global characteristics and a presumed slow evolution of the local environment for 32 min after Huygens landing, in order to explain how and why both ELF and MI signals have changed during the first 10 to 15 s after the impact.

Fig. 2 .
Fig. 2. Amplitude distribution (left) and waveform (right) of the first hundred V ADC samples of bin 48 of the Cruise checkout # 10.The bar levels represent the fractional amount of samples distributed among 6 classes about 6.26 mV RMS wide each, distributed between the min and max values of the series.The solid red line is the theoretical shape of a normal distribution with expectation equal to µ.

Fig. 3 .
Fig. 3. Frequency distribution of two ELF data sequences with the SNR of bin 36 respectively larger (left panel), and smaller than 1 (right panel).Crosses are mean values and standard deviations.Solid blue lines are LFR fits, excluding bin 36 (see text).Solid red line is the LFR fit of smallest values of the surface sequence excluding bin 36.Dashed lines emphasize the gap between measurements and LFR fits.

Fig. 5 .
Fig. 5. Frequency variation of the normalized skewed standard deviation σ f for the same sequences as in Figs. 3 and 4.

Fig. 6 .
Fig. 6.Global summary chart of simulation S1.Highlighted yellow areas cover the bin 36 experimental ranges of mean amplitude µ, standard deviation σ and normalized skew S, during the ground sequence.The blue disc (middle panel) demarcates the most probable area for SNR.Green, red and blue crosses correspond to simulation runs based upon surface noise bins 12, 24 and 48, respectively.

Fig. 7 .
Fig. 7. Distribution of the smallest value of bin 36 derived from simulation S1; the dark-brown bars correspond to the absence of signal (see text).

Fig. 9 .
Fig. 9. Distribution of mean and standard deviation of LFR noise files derived from 10 000 pairs of 418 samples used in simulation S2, assuming a normal distribution with expected values of µ N and σ N equal to 30 and 8.7 mV RMS, respectively.

Fig. 10 .
Fig. 10.Distribution of mean value µ of the composite bin 36 versus SNR derived from 10 000 runs of simulation S2.The most likely value (SNR ∼ 0.35) lies at the intersection between the green curve and the bright blue line.

Fig. 11 .
Fig. 11.Distribution of the standard deviation σ of the composite bin 36 versus SNR.The most likely value (SNR ∼ 0.225) lies at the intersection between the green curve and the bright blue line as in Fig. 10.Black lines at 7.3 and 9.3 mV are min and max values, respectively, of noise amplitude in the absence of signal.