© ESO, 2018

1. Introduction

The Planetary Radio Interferometry and Doppler Experiment (PRIDE) technique exploits the radio (re-)transmitting capabilities of spacecraft from the most modern space science missions (Duev et al. 2012). A very high sensitivity of Earth-based radio telescopes involved in astronomical and geodetic Very Long Baseline Interferometry (VLBI) observations and an outstanding signal stability of the radio systems allow PRIDE to conduct precise tracking of planetary spacecraft. The data from individual telescopes are processed both separately and jointly to provide Doppler and VLBI observables, respectively. Although the main product of the PRIDE technique is the VLBI observables (Duev et al. 2016, Paper I of this series), the accurate examination of the changes in phase of the radio signal propagating from the spacecraft to each of the ground radio telescopes on Earth make the open-loop Doppler observables derived from each telescope very useful for different fields of planetary research.

Dedicated deep space tracking systems, for example, NASA’s Deep Space Network (DSN) and ESA’s tracking station network (Estrack), provide data to determine the precise state vector of a spacecraft, based on the spacecraft signal detected at the ground-based receivers. To this end, the tracking systems can provide a variety of radiometric data (e.g., Doppler, range and interferometry data) under different operational schemes (Thornton & Border 2003). The type of tracking data needed for a particular spacecraft depends on the mission stage of the spacecraft and for which means these data will be used. However, this does not imply that several tracking data types cannot be used for the same purpose. In fact, the use of different precise and reliable tracking techniques not only enables a more challenging navigation performance, but could enhance various scientific experiment carried out during the mission (Martin-Mur et al. 2006; Mazarico et al. 2012; Iess et al. 2014).

The Doppler effect due to the relative motion of the radio elements can be retrieved from the signal received at the ground station in different ways. One way is to measure the changes in light travel time of the received spacecraft signal with a closed-loop mechanism (Kinman 2003). In this tracking scheme, the received signal at the station is mixed with a local oscillator signal. Once the carrier frequency is found at the receiving station, a numerically controlled oscillator is set at the same value of the detected frequency and the carrier loop is closed (Tausworthe 1966; Gupta 1975). The bandwidth of the loop is gradually reduced to a preset operational value using its feedback mechanism. Once the “phase-lock” is acquired, the resulting Doppler shifted beat frequency is input into a Doppler cycle counter. The cycle counter measures the total phase change of the Doppler beat over a count interval, thus yielding the change in range over the count interval. The output, i.e., the Doppler cycle count, consists of an integer number from the Doppler counter itself and a fractional term from a Doppler resolver and is used to reconstruct the received spacecraft frequencies, also known as sky frequencies (Morabito & Asmar 1995; Moyer 2005). The precision at which these measurements can be obtained, is limited by the way the time is tagged (i.e., the quality of the timing standards) and the signal-to-noise ratio (S/N) of the measurements. Because of the mechanism used, the data derived is commonly known as Doppler closed-loop data.

The straightforwardness of this technique and the real-time availability of the data make closed-loop Doppler tracking the preferred tracking scheme when performing navigation and telemetry measurements with the DSN and Estrack networks. However, for radio science applications this is not necessarily the case. The term “radio science” includes all the scientific information that can be derived from the interaction of the spacecraft signal with planetary bodies and interplanetary media as it propagates from the spacecraft to Earth (Tyler et al. 1989; Howard et al. 1992; Kliore et al. 2004; Pätzold et al. 2004; Häusler et al. 2006; Iess et al. 2009). In some scenarios, for instance planetary atmospheric occultation (Jenkins et al. 1994; Tellmann et al. 2009, 2013) and ring occultation (Marouf et al. 1986), the received signal can present abrupt changes in frequency and amplitude, yielding a loss-of-lock in a closed-loop tracking scheme. For such cases, an open-loop receiver is preferable. In this case, no real-time signal detection mechanism is present, but instead the frequency spectrum of the detected signal is downconverted, digitized, and recorded with a sufficiently wide bandwidth to be able to capture the high dynamics of the signal (Kwok 2010). The data processing is performed at a later stage with a digital phase-lock loop (PLL), which simulates the real-time PLL-controlled system used in the closed-loop receivers, and a fast Fourier transform (FFT) that estimates the frequency and amplitude of the received signal. The difference resides in the ability of the digital PLL of starting new locking processes once the system is considered out of lock, and the direct estimation of the frequency of the carrier tone at each sampling time. This mechanism allows an observer to directly reconstruct the sky frequency of parts of the detected signal that would be otherwise considered lost. For the post-processing of the open-loop Doppler, although it relies on the same main detection methods (PLL and FFT), there are various spectral analyses approaches that can be used (Lipa & Tyler 1979; Tortora et al. 2002; Paik & Asmar 2011; Jian et al. 2009).

We present our approach for deriving Doppler open-loop data with the PRIDE technique, using a set of radio telescopes from the European VLBI Network (EVN) and the Very Large Baseline Array (VLBA). Although these telescopes are typically used for observations of natural cosmic radio sources, ranging from nearby stars to distant quasars, we demonstrated in the past that our approach, based on precision wideband spectral analysis, is capable of tracking planetary spacecraft signals (Witasse et al. 2006; Duev et al. 2012; Molera Calvés et al. 2014; Duev et al. 2016). Since their conception, the equipment and data acquisition software of the DSN and VLBI networks have been developed in close collaboration between the two scientific communities. For this reason, the characteristics and capabilities of the VLBI network receivers and the DSN/Estrack open-loop receivers, also known as radio science receivers, are very similar. The post-processing techniques, however, may differ even between radio science teams using the same network because once the data are recorded the tracking center delivers these data to the science teams, who use their own software for data processing and analysis.

For these reasons, the two goals of this paper are as follows. First, we seek to present our processing technique to derive the open-loop Doppler data, and provide a clear formulation of the observed and computed Doppler observables and a noise budget of the derived tracking data. In this way we analyze the quality of open-loop Doppler data derived with VLBI telescopes through the PRIDE technique and compare it to the standards of the closed-loop Doppler data provided by the DSNs. Second, since PRIDE uses another network of ground stations, this technique allows for the possibility of acquiring precise Doppler open-loop data independently from the tracking networks of the space agencies. For instance, one potential application is to use PRIDE, with the EVN and VLBA networks, to track spacecraft for navigation and telemetry purposes only when there are closed-loop tracking passes scheduled by the tracking facilities of the corresponding agency. In this way, PRIDE could allow radio science activities to be conducted in parallel by performing shadow tracking on the spacecraft signal. These goals are addressed in this paper in the framework of the PRIDE tracking of the ESA Mars Express (MEX) spacecraft during its flyby of Phobos in December 2013.

The MEX orbiter was launched in June 2, 2003 and has been orbiting the red planet since December 2003 in a highly elliptical polar orbit, with 86° inclination, periaerion of ~ 300 km, apoaerion of ~ 10 100 km, and an orbital period of 6.7 h. Owing to its highly valuable science return the mission has been extended six times beyond its nominal mission duration (Chicarro et al. 2004). The MEX telemetry, tracking, and command (TT&C) subsystem operates in a two-way mode, receives the transmitted uplink signal in X band (7.1 GHz), and provides coherent dual-frequency carrier downlinks at X band (8.4 GHz) and S band (2.3 GHz) via the spacecraft’s 1.8 m high gain antenna (HGA) for all radio science operations of the Mars Express Radio Science Experiment (MaRS) team (Pätzold et al. 2004). On the ground MaRS activities are supported by the 35 m ESA Estrack New Norcia (NNO) station and the 70 m NASA DSN stations, which are all equipped with hydrogen masers as part of the frequency and timing systems. On December 29, 2013, MEX performed a Phobos fly-by at a distance of ~ 45 km from its surface. Under the European Satellite Partnership for Computing Ephemerides (ESPaCE) consortium an opportunity was offered to track the spacecraft with the PRIDE technique using VLBI stations alongside the customary Estrack and DSN stations. The tracking session lasted for 25 h around the flyby event, using 31 VLBI stations around the world, which are also equipped with hydrogen masers as frequency standards, observing at X-band (channel starting at 8412 MHz, recording bandwidth of 16 MHz) in a three-way mode. The PRIDE setup for this particular tracking experiment is described in detail in Duev et al. (2016; Paper I).

The paper is organized as follows. In Sect. 2, the processing pipeline to extract the Doppler detections from the raw open-loop data is explained along with the formulation of the observed and computed values of the instantaneous Doppler observables. In Sect. 3, the observations and the open-loop Doppler detections obtained from the ESA MEX spacecraft in December 2013 are discussed. The quality of the PRIDE Doppler detections is assessed by comparing the Doppler noise obtained by the multiple VLBI stations involved in the experiment with the noise of the Estrack and DSN stations during the same tracking session. The main contributing noise sources are quantitatively discussed. Section 4 summarizes the results and discusses how the findings can improve the planning and enhance the science return of future radio science experiments with PRIDE.

2. PRIDE Doppler observables

In the nominal MEX gravimetry experiments, the orbit perturbations caused by the gravitational fields of Mars and, in this particular case, of Phobos were determined via precise two-way radio Doppler tracking of MEX with dual-frequency downlink during pericenter passes with the Estrack and DSN stations (Hagermann & Pätzold 2009). However, for the MEX Phobos flyby on December 29, 2013, the PRIDE joined the tracking effort in a three-way mode to assess the performance of the technique; in the three-way mode, the signal received was re-transmitted by the spacecraft with a network of radio telescopes none of which is the initial transmitting ground station, also known as shadow tracking.

2.1. Observed values of the Doppler observables

The transmitting/receiving systems at DSN and Estrack stations used for spacecraft radiometric tracking can operate in a closed- or open-loop manner. In these networks, the primary receiver is the closed-loop receiver, which uses a mechanism to phase lock onto the received carrier signal. In this setup, the receiver passband is continuously aligned to the peak of the carrier tone and its bandwidth is gradually narrowed, allowing the retrieval of real-time tracking data and telemetry (Kinman 2003). The open-loop receivers, on the other hand, do not have such a feedback mechanism (Kwok 2010), hence the bandwidth of the receiver passband is predefined and remains fixed during each observation. For this reason, the carrier signal filtering and tracking is performed at a later stage using the accurately timed detection of the signal recorded at the ground station. The radio telescopes used in PRIDE only operate in the open-loop mode. At each station, the received signals are amplified, heterodyned to the baseband, digitized, time-tagged, and recorded onto disks via the standard VLBI data acquisition systems with Mark5 A/B or FlexBuff recording systems (Lindqvist & Szomoru 2014). For the data processing, the disks can be shipped or the data are transferred directly via high-speed networks to the VLBI data processing center at the Joint Institute for VLBI ERIC (JIVE) in the Netherlands.

The Doppler detections are extracted from the raw open-loop data via the PRIDE spacecraft tracking software, consisting of three packages SWSpec, SCtracker, and dPLL¹ (Molera Calvés et al. 2014). In the nominal PRIDE setup we observe two sources, the spacecraft signal and natural radio sources, which are used as calibrators. A large number of the natural sources observed with radio telescopes emit broadband electromagnetic radiation spanning many gigahertz in the frequency domain, however the signal is typically weak. It is therefore desirable to use as wide a frequency band as possible in order to detect the signal. The open-loop receiver systems of the VLBI stations are typically set up to record 4, 8, 16, or 32 frequency channels with 4, 8, 16, or 32 MHz bandwidth per sub-band. However, the spacecraft signal spectrum takes up only a fraction of the sub-band (see Fig. 1a). For this reason, the first processing step is to extract the narrowband containing the spacecraft signal carrier and/or tones present in the spectrum. The SWSpec extracts the data from the channel where the spacecraft signal is located, and subsequently performs a window-overlapped add (WOLA) direct Fourier transform (DFT), followed by a time integration over the obtained spectra. The result is an initial estimate of the spacecraft carrier tone along the observation scan (Fig. 1b). As shown in Fig. 1b, the detected carrier tone has a moving phase throughout the scan, which is caused by the change in relative velocity between the spacecraft and the receiver. The goal is to extract the Doppler shift, first by fitting the changing frequency of the carrier tone by a n-order polynomial, and then using the fit to stop the moving phase of the tone. The latter step is performed with the SCtracker software, which subsequently allows the tracking, filtering, and extraction of the tone in a narrowband. Figure 1c shows the narrowband output signal of the SCtracker. At this point, the spacecraft signal is in a band of a few kHz bandwidth, in contrast to the initial 4–32 MHz bandwidth sub-band. The final step is conducted by the digital phase-locked-loop (dPLL), which performs high-precision reiterations of the previous steps – time-integration of the overlapped spectra, phase polynomial fitting, and phase-stopping correction – on the narrowband signal. After the phase-stopping correction, the power spectrum is accumulated for a selectable averaging interval. Using a frequency window around the tone, the maximum value of the accumulated spectrum is determined. The corresponding frequency of the peak of the spectrum is stored, using as the time tag the middle of the averaging interval. This procedure is conducted throughout the whole range of spectra. The output of the dPLL is the filtered down-converted signal (Fig. 1d) and the final residual phase in the stopped band with respect to the initial phase polynomial fit. The bandwidth of the output detections is typically about 20 Hz with a frequency spectral resolution of ~ 2 mHz (Molera Calvés 2012).

The PRIDE post-processing pipeline allows us to determine the instantaneous Doppler shift of the recorded tracking data, which is different from the integrated Doppler observables that are derived from closed-loop tracking data. For the purposes of orbit determination and the estimation of physical parameters of a celestial body using the Doppler data, it is important that this difference is taken into account when defining the observed and computed values of the Doppler observable. In the closed-loop case, the Doppler observables are derived by computing the change in the accumulated Doppler cycle counts from the spacecraft carrier phase measurements over a time interval at the receiver, as explained in detail in Sect. 13.3 of Moyer (2005). The corresponding modeled values are obtained by taking the difference in range at the beginning and end of the time interval. In the open-loop case, performed by PRIDE as explained in the previous paragraph, the observed values of the instantaneous Doppler observable (for one-way or three-way mode) are derived directly from an estimate of the carrier tone frequency of the spacecraft spectrum. Therefore, the observables are simply retrieved by adding the base frequency f_base (which for the experiment analyzed in this paper was 8412 MHz, as shown in Fig. 1a) of the channel containing the spacecraft signal and the obtained time averaged tone frequencies f_tone at each sampled time t_i, ${\mathit{f}}_{\mathit{R}}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{base}}\mathrm{+}{\mathit{f}}_{\mathrm{tone}}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathit{,}$(1)where f_R is the received frequency.

The integration time is defined by the number of FFT points used at the dPLL on the ~2 kHz bandwidth signal (Fig. 1c). For the gravity field determination experiments, the desired integration time is ~10 s, hence 20 000 FFT points are used in dPLL. The uncertainty of each tone frequency estimate is derived from the final residual phase of the dPLL output.

2.2. Computed values of the Doppler observables

To process the Doppler data obtained with PRIDE, a model is required that provides the instantaneous Doppler shift f_R/f_T, where f_R and f_T denote the observed frequency of the received and transmitted electromagnetic signal, respectively. Fundamentally, this frequency ratio is obtained from $\begin{array}{ccc}\frac{{\mathit{f}}_{\mathit{R}}}{{\mathit{f}}_{\mathit{T}}}\mathrm{=}\frac{\mathrm{d}{\mathit{\tau }}_{\mathit{T}}}{\mathrm{d}{\mathit{\tau }}_{\mathit{R}}}\mathrm{=}{\left(\frac{\mathrm{d}\mathit{\tau }}{\mathrm{d}\mathit{t}}\right)}_{\mathit{T}}\frac{\mathrm{d}{\mathit{t}}_{\mathit{T}}}{\mathrm{d}{\mathit{t}}_{\mathit{R}}}{\left(\frac{\mathrm{d}\mathit{t}}{\mathrm{d}\mathit{\tau }}\right)}_{\mathit{R}}\mathit{,}& & \end{array}$(2)where τ and t denote proper time of the observer and coordinate time, respectively. The R and T subscripts denote properties of receiver and transmitter. The coordinate times of transmission and reception are related via the light-time equation $\begin{array}{ccc}{\mathit{t}}_{\mathit{R}}\mathrm{-}{\mathit{t}}_{\mathit{T}}\mathrm{=}\frac{\mathrm{1}}{\mathit{c}}\left|{{x}}_{\mathit{R}}\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathrm{\right)}\mathrm{-}{{x}}_{\mathit{t}}\mathrm{\left(}{\mathit{t}}_{\mathit{T}}\mathrm{\right)}\right|\mathrm{+}\mathrm{\Delta }\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathit{T}}\mathrm{\right)}\mathit{,}& & \end{array}$(3)where x_R(t) and x_T(t) are the barycentric positions of the receiver and transmitter and Δ(t_R,t_T) the relativistic correction to the light travel times.

Fig. 1 Example of Doppler data processing pipeline using observations of MEX during Phobos flyby by Hartebeesthoek (see Sect. 3). Panel a shows the typical resulting average power spectrum of a scan after running the SWSpec software. In the 16 MHz pre-defined sub-band starting at 8412 MHz, the spacecraft signal is found in the spectrum. A narrowband containing the moving phase of the spacecraft carrier/tone is selected (in this case of 50 kHz bandwidth) to model the Doppler shift using an n-order polynomial frequency fit. Panel b shows a zoom of the spectrum inside the selected narrowband window to perform the fit. Here the moving phase of the carrier tone is visible along the duration of the scan. After the fit is performed, SCtracker applies the polynomial coefficients after converting the sample to the baseband sample to stop the moving phase of the tone. In this way, SCtracker extracts an initial fit of the Doppler shift. Panel c shows the output of the SCtracker, which is a phase-stopped filtered out signal in a 2 kHz narrowband in baseband. Finally, the dPLL performs high-precision iterations of the time integration of the overlapped spectra, phase polynomial fitting, conversion to baseband, and phase-stopping corrections, using narrowband windows around the carrier tone. The iterations stop when the window bandwidth reaches 20 Hz, as shown in panel d, allowing the extraction of the frequency and phase residuals of the spacecraft carrier tone with a 2 mHz frequency spectral resolution.

The main complication in obtaining an explicit expression from Eq. (2) is to compute the terms dΔ(t_R,t_T)/dt_R and dΔ(t_R,t_T)/dt_T. To expand these equations, we use the formalism of Kopeikin & Schäfer (1999), where it is assumed that:

• The metric g_αβ can be expanded to post-Minkowskian order, so that g_αβ(x,t) = η_αβ + h_αβ(x,t), where η_αβ is the Minkowski metric and the metric perturbation h_αβ = O(G).

• All bodies with gravity fields that perturb the null geodesic of the electromagnetic signal can be modeled as point masses.

• All bodies with gravity fields that perturb the null geodesic of the electromagnetic signal have a constant barycentric velocity over the relevant time interval of a single measurement.

Under these assumptions, Δ(t_R,t_T) reduces to the following, neglecting second order terms in v/c: $\begin{array}{ccc}\mathrm{\Delta }\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathit{T}}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{2}\mathit{G}}{{\mathit{c}}^{\mathrm{3}}}\sum _{\mathit{a}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{m}}_{\mathit{a}}(\ln \frac{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{s}}_{\mathit{R}}\mathrm{\right)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{s}}_{\mathit{R}}\mathrm{\right)}}{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}{\mathit{t}}_{\mathit{T}}\mathit{,}{\mathit{s}}_{\mathit{T}}\mathrm{\right)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}{\mathit{t}}_{\mathit{T}}\mathit{,}{\mathit{s}}_{\mathit{T}}\mathrm{\right)}}\\ & & \mathrm{-}\left({k}\mathrm{·}\frac{{{v}}_{\mathit{a}}\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{s}}_{\mathit{R}}\mathrm{\right)}}{\mathit{c}}\right)\ln \mathrm{\left(}{\mathit{r}}_{\mathit{a}}\mathrm{\right(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{s}}_{\mathit{R}}\mathrm{)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{s}}_{\mathit{R}}\mathrm{\left)}\mathrm{\right)}\\ & & \mathrm{+}\left({k}\mathrm{·}\frac{{{v}}_{\mathit{a}}\mathrm{\left(}{\mathit{t}}_{\mathit{T}}\mathit{,}{\mathit{s}}_{\mathit{T}}\mathrm{\right)}}{\mathit{c}}\right)\ln \mathrm{\left(}{\mathit{r}}_{\mathit{a}}\mathrm{\right(}{\mathit{t}}_{\mathit{T}}\mathit{,}{\mathit{s}}_{\mathit{T}}\mathrm{)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{(}{\mathit{t}}_{\mathit{T}}\mathit{,}{\mathit{s}}_{\mathit{T}}\mathrm{\left)}\mathrm{\right)})\end{array}$(4)the static case of which is also known as the Shapiro effect (Shapiro et al. 1971). In the above, the parameter s denotes the retarded time of body a w.r.t. either the signal transmission or signal reception for T and R subscripts, respectively. This time parameter is obtained from the light-time equation of Eq. (3), only now by considering the perturbing body a as the transmitting body, so that $\begin{array}{ccc}& & \\ & & \end{array}$Physically, these times represent the time at which the gravitational signal from body a must be evaluated for its effect on the transmitter at t_T and the receiver at t_R to be modeled, implicitly assuming c to be the speed of gravity. Although we omit any a sub/superscript of the times s, we stress that these times are different for each perturbing body a.

Under these above assumptions, as shown by Kopeikin & Schäfer (1999), Eq. (2) can be written as ${\mathit{f}}_{\mathit{R}}\mathrm{=}{\mathit{f}}_{\mathit{T}}\frac{\mathrm{1}\mathrm{-}{k}\mathrm{·}{{v}}_{\mathit{R}}\mathit{/}\mathit{c}}{\mathrm{1}\mathrm{-}{k}\mathrm{·}{{v}}_{\mathit{T}}\mathit{/}\mathit{c}}\mathit{R}\mathrm{\left(}{{v}}_{\mathit{R}}\mathit{,}{{v}}_{\mathit{T}}\mathit{,}{{t}}_{\mathit{R}}\mathit{,}{{t}}_{\mathit{T}}\mathrm{\right)}\mathit{,}$(7)where v_R and v_T are the barycentric velocity vectors of the receiving station at reception time t_R and at transmission time t_T, respectively. The term R denotes a set of (special and general) relativistic corrections. The unit vector k is the direction along which the radio wave propagates at past null infinity (i.e., when following the signal back along the null geodesic to t → −∞), which can be expressed as ${k}\mathrm{=}\mathrm{-}{K}\mathrm{-}{\beta }\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{s}}_{\mathit{R}}\mathrm{\right)}\mathrm{+}{\beta }\mathrm{\left(}{\mathit{t}}_{\mathit{T}}\mathit{,}{\mathit{s}}_{\mathit{T}}\mathrm{\right)}\mathit{,}$(8)where K is the geometric direction of the propagation of the electromagnetic wave in a flat space-time and β(t_R,s_R) and β(t_T,s_T) are the relativistic corrections as a function of the states of body. a at the retarded times of reception and transmission of the electromagnetic signal, respectively. These vectors are defined as follows, $\begin{array}{ccc}{K}& \mathrm{=}& \\ {\mathit{\beta }}^{\mathit{i}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{2}\mathit{G}}{\mathrm{|}{{x}}_{\mathit{T}}\mathrm{-}{{x}}_{\mathit{R}}\mathrm{|}{\mathit{c}}^{\mathrm{2}}}\\ & & \mathrm{\times }\sum _{\mathit{a}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{m}}_{\mathit{a}}\left[\frac{\mathrm{1}\mathrm{-}{k}\mathrm{·}{{v}}_{\mathit{a}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}\mathit{c}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{a}}^{\mathrm{2}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}{\mathit{c}}^{\mathrm{2}}}}\frac{{\mathit{r}}_{\mathit{a}}^{\mathit{i}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{\mathit{k}}^{\mathit{i}}\mathrm{(}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{(}\mathit{t,s}\mathrm{\left)}\mathrm{\right)}}{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}}\right]\\ & & \mathrm{-}\frac{\mathrm{4}\mathit{G}}{\mathrm{|}{{x}}_{\mathit{T}}\mathrm{-}{{x}}_{\mathit{R}}\mathrm{|}{\mathit{c}}^{\mathrm{2}}}\sum _{\mathit{a}\mathrm{=}\mathrm{1}}^{\mathit{N}}[\frac{{\mathit{m}}_{\mathit{a}}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{a}}^{\mathrm{2}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}{\mathit{c}}^{\mathrm{2}}}}\\ & & \mathrm{\times }\mathrm{[}{\mathit{v}}_{\mathit{a}}^{\mathit{i}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}\mathit{c}\mathrm{-}{\mathit{k}}^{\mathit{i}}\mathrm{(}{k}\mathrm{·}{{v}}_{\mathit{a}}\mathrm{(}\mathit{s}\mathrm{)}\mathit{/}\mathit{c}{\mathrm{)}}^{\mathrm{]}}\ln \left({\mathit{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\right)]\mathit{.}\\ & & \end{array}$The relativistic term R in Eq. (7) can be decomposed as follows: $\mathit{R}\mathrm{\left(}{\mathit{v}}_{\mathit{R}}\mathit{,}{\mathit{v}}_{\mathit{T}}\mathit{,}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathit{T}}\mathrm{\right)}\mathrm{=}{\left[\frac{\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{v}}_{\mathit{T}}\mathit{/}\mathit{c}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{v}}_{\mathit{R}}\mathit{/}\mathit{c}{\mathrm{)}}^{\mathrm{2}}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}{\left[\frac{\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathit{T}}\mathrm{\right)}}{\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathrm{\right)}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{b}\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathrm{\right)}}{\mathit{b}\mathrm{\left(}{\mathit{t}}_{\mathit{T}}\mathrm{\right)}}\mathit{,}$(11)where the first term accounts for the special relativistic Doppler shift, the second term accounts for the general relativistic corrections due to the dτ/ dt terms in Eq. (2), and the final term is (along with the terms β given above) a result of expanding dΔ/dt when inserting Eq. (3) into the middle term on the right-hand side of Eq. (2). The terms a and b are given by $\begin{array}{ccc}\mathit{a}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{+}\frac{\mathrm{2}\mathit{G}}{{\mathit{c}}^{\mathrm{2}}}\sum _{\mathit{a}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{{\mathit{m}}_{\mathit{a}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{a}}^{\mathrm{2}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}{\mathit{c}}^{\mathrm{2}}}}{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{{v}}_{\mathit{a}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathit{/}\mathit{c}}\\ & & \mathrm{-}\frac{\mathrm{4}\mathit{G}}{{\mathit{c}}^{\mathrm{2}}\mathrm{-}{\mathit{v}}^{\mathrm{2}}}\sum _{\mathit{a}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{{\mathit{m}}_{\mathit{a}}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{a}}\mathrm{(}\mathit{s}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{c}}^{\mathrm{2}}}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{v}\mathrm{(}\mathit{t}\mathrm{)}\mathrm{·}{{v}}_{\mathit{a}}\mathrm{(}\mathit{s}\mathrm{)}\mathit{/}{\mathit{c}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{{v}}_{\mathit{a}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathit{/}\mathit{c}}\\ \mathit{b}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{+}\frac{\mathrm{2}\mathit{G}}{{\mathit{c}}^{\mathrm{2}}}\sum _{\mathit{a}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{{\mathit{m}}_{\mathit{a}}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{a}}^{\mathrm{2}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}{\mathit{c}}^{\mathrm{2}}}}\frac{\mathrm{1}\mathrm{-}{k}\mathrm{·}{{v}}_{\mathit{a}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}\mathit{c}}{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{{v}}_{\mathit{a}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathit{/}\mathit{c}}\\ & & \mathrm{\times }[\frac{\mathrm{(}\mathrm{1}\mathrm{-}{k}\mathrm{·}{{v}}_{\mathit{a}}\mathrm{(}\mathit{s}\mathrm{)}\mathit{/}\mathit{c}\mathrm{)}\mathrm{(}{k}\mathrm{\times }{v}\mathrm{(}\mathit{t}\mathrm{)}\mathit{/}\mathit{c}\mathrm{)}\mathrm{·}\mathrm{(}{k}\mathrm{\times }{{r}}_{\mathit{a}}\mathrm{(}\mathit{t,s}\mathrm{\left)}\mathrm{\right)}}{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}}\\ & & \mathrm{-}\frac{\mathrm{(}{k}\mathrm{\times }{{v}}_{\mathit{a}}\mathrm{(}\mathit{s}\mathrm{)}\mathit{/}\mathit{c}\mathrm{)}\mathrm{·}\mathrm{(}{k}\mathrm{\times }{{r}}_{\mathit{a}}\mathrm{(}\mathit{t,s}\mathrm{\left)}\mathrm{\right)}}{{\mathit{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}\mathrm{-}{k}\mathrm{·}{{r}}_{\mathit{a}}\mathrm{\left(}\mathit{t,s}\mathrm{\right)}}\mathrm{+}{k}\mathrm{·}{{v}}_{\mathit{a}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{/}\mathit{c}]\mathit{.}\end{array}$Evaluating this algorithm for the one-way Doppler case requires 2N + 1 solutions of light-time equations, once for Eq. (3) and N for both Eqs. (5) and (6), with N the number of bodies perturbing the path of the signal. These equations are implicit and must be solved iteratively. For Eq. (3), we must compute t_T from a given t_R. We initialize t_T(1) = t_R and iterate to find t_{T(n + 1)} from t_T(n) using the Newton-Raphson method as follows: ${\mathit{t}}_{\mathit{T}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathit{t}}_{\mathit{T}\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{-}\frac{{\mathit{t}}_{\mathit{R}}\mathrm{-}{\mathit{t}}_{\mathit{T}\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{|}{{x}}_{\mathit{R}}\mathrm{-}{{x}}_{\mathit{T}}\mathrm{|}}{\mathit{c}}\mathrm{-}{\mathrm{\Delta }}_{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}{\mathit{t}}_{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathit{T}\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\right)}}{\frac{{{x}}_{\mathit{R}}\mathrm{-}{{x}}_{\mathit{T}}}{\mathrm{|}{{x}}_{\mathit{R}}\mathrm{-}{{x}}_{\mathit{T}}\mathrm{|}}\frac{{{v}}_{\mathit{T}}\mathrm{\left(}{\mathit{t}}_{\mathit{T}\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\right)}}{\mathit{c}}\mathrm{-}\mathrm{1}}\mathrm{·}$(14)The iterative procedure converges when | t_{T(n + 1)}−t_T(n) | ≤ ϵ for some predefined small ϵ. For the evaluation of the term Δ_n, the unit vector k is also updated iteratively using Eqs. (8)–(10), initially assuming k₍₁₎ = −K.

After the convergence of t_T, the values for k, r_a(t_T,s_T), v_a(t_T), x(t_T), and v(t_T) are computed and the values for the general relativistic corrections a(t_T), a(t_R), b(t_T), and b(t_R) are determined. Finally, the instantaneous one-way Doppler frequency at reception time t_R can be determined from Eq. (7).

The explicit formula for the two- and three-way observables for the propagation of the radio signal emitted from the ground station on Earth, with position x_T and velocity v_T at t_T, then received and transponded back to Earth by a spacecraft with position x_S and velocity v_S at t_S, where the superscripts “+” and “–” denote received on uplink and transponded on downlink, and finally received at a ground station on Earth, with position x_R and velocity v_R at t_R, is written as ${\mathit{f}}_{\mathit{R}}\mathrm{=}{\mathit{f}}_{\mathit{T}}\mathit{M}\left(\frac{\mathrm{1}\mathrm{-}{{k}}^{\mathrm{+}}\mathrm{·}{{v}}_{\mathit{S}}^{\mathrm{+}}}{\mathrm{1}\mathrm{-}{{k}}^{\mathrm{+}}\mathrm{·}{{v}}_{\mathit{T}}}\mathit{R}\mathrm{\left(}{\mathit{v}}_{\mathit{T}}\mathit{,}{\mathit{v}}_{\mathit{S}}^{\mathrm{+}}\mathit{,}{\mathit{t}}_{\mathit{T}}\mathit{,}{\mathit{t}}_{\mathit{S}}^{\mathrm{+}}\mathrm{\right)}\right)\frac{\mathrm{1}\mathrm{-}{{k}}^{\mathrm{-}}\mathrm{·}{{v}}_{\mathit{R}}}{\mathrm{1}\mathrm{-}{{k}}^{\mathrm{-}}\mathrm{·}{{v}}_{\mathit{S}}^{\mathrm{-}}}\mathit{R}\mathrm{\left(}{\mathit{v}}_{\mathit{S}}^{\mathrm{-}}\mathit{,}{\mathit{v}}_{\mathit{R}}\mathit{,}{\mathit{t}}_{\mathit{S}}^{\mathrm{-}}\mathit{,}{\mathit{t}}_{\mathit{R}}\mathrm{\right)}\mathit{,}$(15)where again all the positions and velocities are given with respect to the solar system barycenter.

The solution to Eq. (15) (the three-way Doppler predictions) is found in a similar manner as for Eq. (7) (the one-way predictions) by first solving the expression for the uplink (in parenthesis in Eq. (15)), this time estimating the signal reception time at the spacecraft ${\mathit{t}}_{\mathit{S}}^{\mathrm{+}}$ and subsequently determining all the uplink parameters corresponding to ${\mathit{t}}_{\mathit{S}}^{\mathrm{+}}$ with the same iteration procedure as explained above. Once the expression in the parenthesis is solved, the downlink part is found as for Eq. (15), estimating the signal transmission time at the spacecraft ${\mathit{t}}_{\mathit{S}}^{\mathrm{-}}$ and determining the corresponding downlink parameters. In Eq. (15), f_T is the frequency transmitted by the ground station at time t_T and M is the corresponding spacecraft turnaround ratio. The terms R⁺ and R⁻ are the relativistic corrections on uplink and downlink, respectively.

In a practical implementation of this algorithm, it is important to explicitly recompute all state vectors at each step, as any expansion quickly becomes inaccurate on typical deep-space craft light travel times.

3. MEX Phobos Flyby: GR035 experiment

The global VLBI experiment was designed to track MEX 14 h prior to and 11 h after its closest-ever Phobos fly-by at approximately 7:21 (UTC) on December 29, 2013. At the time of the experiment, Mars was at a distance of ~ 1.4 AU from the Earth with a solar elongation of ~ 87°. During the 25 h, MEX was tracked by the Estrack New Norcia (NNO; western Australia) station, DSN DSS-63 (Robledo, Spain) and DSS-14 (Goldstone, California, USA) stations, and 31 VLBI radio telescopes around the world. The latter were organized through the global VLBI experiment GR035. The experimental setup of GR035 was presented in Duev et al. (2016). During the first nine hours, NNO was the transmission station, followed by eight hours of tracking with DSS-63, and finally eight hours with DSS-14. The distribution of the telescopes over the duration of the experiment is presented in Figs. 2, 3, and 7 in Duev et al. (2016). The spacecraft operated in the two-way mode with an X-band uplink (7.1 GHz) and dual simultaneous S/X-band downlink (2.3/8.4 GHz) transponded by the HGA pointing toward the Earth. The Estrack and DSN stations produced two-way S- and X-band Doppler closed-loop data products. From the 31 VLBI stations involved in the experiment, only the detections of 25 stations were used for the analysis presented in this paper, as listed in Table 1, owing to different technical problems with the remaining six stations during the observation.

We formed the Doppler residuals by differencing the Doppler detections, obtained as explained in Sect. 2.1, with the predicted Doppler, derived as explained in Sect. 2.2 via the latest MEX navigation post-fit orbit of December 2013 provided by the European Space Operations Centre (ESOC)². In order to correctly determine the Doppler noise of the residuals, we flagged some data out. In particular, the data obtained during occultations and the flyby event were discarded. Additionally, we discarded a number of scans that presented systematic outliers, for instance, when the transmission mode at the ground station changed. We found that the Doppler residuals obtained with the VLBI stations are in agreement with the DSN and Estrack residuals. Figure 2 shows as an example the frequency residuals found with the 25 m antenna of the Very Long Baseline Array (VLBA) at Kitt Peak (designation Kp) and the residuals of the 70 m DSS-63 and DSS-14 antennas. In this case, the median value of the difference between the fit of Kp and the fit of DSS-63 and DSS-14, respectively, remains below 1 mHz for an integration time of 10 s. For other VLBI stations, the median (after flagging) of the Doppler residuals was found to be approximately 2 mHz.

In order to determine the quality of the PRIDE Doppler detections of the different VLBI stations involved in this experiment, first we had to identify the different sources that contribute to the overall noise of the Doppler residuals. The signal received at the ground stations have random errors introduced by the instrumentation on board the spacecraft and at the receiving system, and the random errors introduced by the propagation of the signal through the different media along the line of sight of the ground station. Additionally, systematic errors can also be introduced, for instance, when calibrating the signal or in the models used to generate the predicted Doppler signature. The calibration of the Doppler observables in relation to the signal delays induced by the ionosphere of the Earth are performed using the total vertical electron content (vTEC) maps available from the International GNSS Service (IGS) on a daily basis with a two-hour temporal resolution on a global grid (Feltens & Schaer 1998). The calibration of the tropospheric signal delays is applied via the Vienna Mapping Functions VMF1 (Boehm et al. 2006) or ray-tracing through the Numerical Weather Models (NWM) (Duev et al. 2011), depending on the antenna elevation. The systematic errors due to the model and orbit used to derive the Doppler residuals and the residual systematic noise resulting from the ionospheric and tropospheric delay calibration, are not characterized in this paper.

The random errors introduced by the instrumentation are analyzed in Sect. 3.1, and the random errors introduced by the propagation of the signal through the interplanetary media are treated in Sect. 3.2. Finally, in Sect. 3.3 the summary of the noise budget is given.

Fig. 2 Comparison of the Doppler residuals obtained with VLBA-Kp (in black), DSS-63 (in blue), and DSS-14 (in red). The median value of the difference between the fit of VLBA-Kp and the fit of DSS-63 and DSS-14, respectively, remains below 1 mHz, for an integration time of 10 s.

3.1. Instrumental noise

The noise budget of the two-way Estrack/DSN Doppler detections and the three-way PRIDE Doppler detections include instrumental noises introduced at the transmitting ground station (electronics, frequency standard, and antenna mechanical noise) and at the spacecraft (electronics), which are common to both observables (Asmar et al. 2005; Iess et al. 2014). Hence, regarding instrumental noises, the difference between these noise budgets resides at the receiving stations: the thermal noise, induced by the ground station receiver system and the limited received downlink power, the frequency and timing systems’ noise, and the antenna mechanical noise. In this paper, only the first two sources of noise are treated since the antenna mechanical noise of the VLBI stations has not yet been characterized for the time intervals relevant to this study.

The thermal noise of the ground station is characterized by the RMS of the random fluctuations of the total system power at the ground station. The one-sided phase noise spectral density S_φ of the received signal gives the relative noise power to the carrier tone, contained in a 1 Hz bandwidth chosen to be centered at a frequency with a large offset Δf from the carrier frequency f_carrier (Vig et al. 1999), ${\mathit{S}}_{\mathit{\phi }}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{sideband}}\mathrm{(}{\mathit{f}}_{\mathrm{carrier}}\mathrm{+}\mathrm{\Delta }\mathit{f}\mathrm{)}}{{\mathit{P}}_{\mathrm{carrier}}}\mathit{,}$(16)where P_carrier is the power of the carrier tone and P_sideband is the power of the 1 Hz bandwidth band.

The S/N is then approximated by 1 /S_φ. As explained in (Barnes et al. 1971; Rutman & Walls 1991), the Allan deviation of white phase noise can be estimated by ${\mathit{\sigma }}_{\mathit{y}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{\approx }\frac{\sqrt{\mathrm{3}\mathrm{B}{\mathit{S}}_{\mathit{\phi }}}}{\mathrm{2}\mathit{\pi }{\mathit{f}}_{\mathrm{0}}\mathit{\tau }}\mathrm{·}$(17)Using Eqs. (17) and (16), and since $\mathit{S}\mathit{/}\mathit{N}\mathrm{\left(}\mathit{\tau ,B}\mathrm{\right)}\mathrm{=}\mathit{S}\mathit{/}\mathit{N}\mathrm{(}\mathrm{1}\hspace{0.17em}\mathrm{s}\mathit{,}\mathrm{1}\hspace{0.17em}\mathrm{Hz}\mathrm{)}\sqrt{\frac{\mathit{\tau }}{\mathit{B}}}$, the Allan deviation for the S/N detections of the different telescopes were determined, as shown in Table 1. During the Phobos flyby science operations conducted with MEX, the two-way closed-loop Doppler data was obtained using a carrier loop bandwidth of 30 Hz at the ground stations with 10 s integration time (Hagermann & Pätzold 2009). For the VLBI telescopes, the three-way open-loop Doppler data was initially recorded in a 16 MHz wide band and then processed with the software described in Sect. 2.1, for a final phase detection of 20 Hz bandwidth, with 10 s integration time. During the MEX orbits, for which NNO and DSS-14 were the transmitting stations, the telemetry was being transmitted except during MEX’s nominal observation phase around the pericenter passage. However, during the orbit where DSS-63 was the transmitting station, in which the Phobos flyby occurred, the telemetry was turned off throughout the whole orbit. For the Doppler error budget, only the time slots planned as radio science passes are taken into account, since the average S/N calculated with Eq. (16) drops when the telemetry is on. In Table 1 we also give the values of the sensitivity of each telescope, which for single-dish radio telescopes is defined as the system equivalent flux density (SEFD)³. This value is useful when comparing the performance between ground stations, since it comprises information about the total noise of the system and the collecting area of the antenna. Also when planning an experiment, it is important to know the nominal SEFD of a station at given frequency, since it can used to compute the expected S/N of a detection.

The 65 m dish Tianma⁴ has the highest S/N detections with a downlink power at reception 15 dB higher than for the smallest stations, i.e., the 12 m Yarragadee, Katherine, Hobart, and Warkworth, as expected because of their smaller collecting area, at elevations higher than 30°. At τ = 10 s the corresponding Allan deviation is of 4.6 × 10^-15. Table 1 uses the average S/N values over the whole coverage of each telescope. However, because of the large variation of the elevation of the antennas during the several hours-long tracking sessions, there are periods of time for which smaller antennas achieve similar S/N levels as larger antennas. Figure 3 shows such an example in which for 3 h the 15 m Hartbeesthoek achieves better S/N levels than the 25 m Urumuqi and similar S/N levels than the 25 m Onsala because of a more favorable antenna elevation (at elevations <20°, several noise contributions at the antenna rapidly increase, such as the atmospheric and spillover noise).

Regarding the noise induced by the frequency and timing systems, the Estrack, DSN, VLBA, and EVN stations are all equipped with hydrogen masers frequency standards, which provide a stability better than 4 < 10^-14 at τ = 10 s (<10^-15 at τ = 1000 s). Hence the noise contributions related to the frequency standard should be on the same order of magnitude for the different networks.

Fig. 3 Signal-to-noise ratio and elevation angle comparison between 25 m Ur, 25 m On, and 15 m Ht radio telescopes from 4:00–7:30 UTC on December 29, 2013 (TX: DSS-63). Left panel shows that although Ht (in pink) has a lower collecting area, it achieves higher S/N levels than Ur (in blue) and similar S/N levels than On (in green). Right panel shows the elevation angle of each antenna during the same time period from which the correlation with the S/N levels of each station is evident.

3.2. Medium propagation noise

The precision of the Doppler detections is also affected by the noise introduced by the propagation of the radio signal through the interplanetary medium, ionosphere, and troposphere. The effects of ionospheric and interplanetary scintillation can be studied using the differenced phases of the signals received in S band, φ_s, and X band, φ_x, ${\mathit{\phi }}_{\mathrm{\Delta }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{\phi }}_{\mathit{s}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{11}}{\mathit{\phi }}_{\mathit{x}}$ (Levy 1977). By subtracting the phases, the contribution of the dispersive plasma scintillations can be isolated. The Allan variance of the differenced phases is related to its two-sided phase power spectrum S_φ(f) (Barnes et al. 1971; Armstrong et al. 1979) by ${\mathit{\sigma }}_{\mathit{y}}^{\mathrm{2}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{S}}_{\mathit{\phi }}\mathrm{\left(}\mathit{f}\mathrm{\right)}\frac{{\mathit{f}}^{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{0}}^{\mathrm{2}}}\frac{{\sin }^{\mathrm{4}}\mathrm{\left(}\mathit{\pi \tau f}\mathrm{\right)}}{\mathrm{(}\mathit{\pi \tau f}{\mathrm{)}}^{\mathrm{2}}}\mathrm{d}\mathit{f}\mathit{.}$(18)As explained in Armstrong et al. (1979), when the phase spectrum can be approximated to S_φ(f) = Af^{− m}, Eq. (18) can be rewritten as ${\mathit{\sigma }}_{\mathit{y}}^{\mathrm{2}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathit{A}{\mathit{\tau }}^{\mathit{m}}}{{\mathit{\pi }}^{\mathrm{2}}{\mathit{\nu }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\tau }}^{\mathrm{3}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\sin {\mathrm{\left(}\mathit{\pi z}\mathrm{\right)}}^{\mathrm{4}}}{{\mathit{z}}^{\mathit{m}}}\mathrm{d}\mathit{z}\mathit{.}$(19)A first-order approximation of the phase spectrum on a logarithmic scale is performed as explained in (Molera Calvés et al. 2014), from which the slope m and the constant A are determined using only the Doppler-mode⁵ observations, where the length of the scan is typically >10 min. For instance, considering the Doppler-mode observations of Hartebeesthoek 15 m antenna, the plasma phase scintillation noise can be characterized. Figures 4 and 5 show the spectral power density of the phase fluctuations of MEX signal in S and X band, respectively. The lower and upper limits of the scintillation band are determined by visual inspection, taking into account the cutoff frequency defined to perform the polynomial fit and amount of fluctuation due to the receiver system noise. The phase scintillation indices obtained with the S- and X-band signal, 0.070 rad and 0.073 rad, correspond to the results for Mars-to-Earth total electron content (TEC) along the line of sight found by Molera Calvés et al. (in prep.), where the dependence of the interplanetary phase scintillation on elongation was studied using various MEX observations. When comparing Figs. 4 and 5, it is noticeable that in the spectral power, the noise band in the S band is much higher than in the X band. This is due to the higher thermal noise of the receiver and larger presence of RFI in this band.

Fig. 4 Spectral power density of MEX signal in S band. The scintillation band extends from 8 mHz to 0.1 Hz, obtaining a value for the slope of –2.471, which is coherent with the spectral index values found by (Woo & Armstrong 1979). The mean phase scintillation index of the signal received from 3h56m to 09h12m (UTC) on 2013-12-29, is 0.070 rad, at an elongation of ~ 87° and distance of ~ 1.4 AU.

Fig. 5 Spectral power density of MEX signal in X band. The scintillation band extends from 8 mHz to 0.45 Hz, obtaining a value for the slope of –2.469. The mean phase scintillation index of the signal received from 3h56m to 09h12m (UTC) on 2013-12-29 is 0.073 rad at an elongation of ~ 87° and distance of ~ 1.4 AU.

Fig. 6 Spectral power density of the differential phases φΔ. The scintillation band extends from 8 mHz to 0.1 Hz, obtaining a value for the slope of –2.372. The mean phase scintillation index of the signal received from 3h56m to 09h12m (UTC) on 2013-12-29 is 0.069 rad at an elongation of ~ 87° and distance of ~ 1.4 AU.

Figure 6 shows the spectral power density of φ_Δ for Ht. The slope found for the scintillation band that extends from 8 mHz to 0.1 Hz is –2.372 with a mean scintillation index of 0.069 rad. These results are in agreement with Molera Calvés et al. (2014). As the phase power spectrum can be described in the form S_φ(f) = Af^{− m}, following Eq. (19) the Allan variance of the plasma phase scintillation is 2.46 × 10^-15 at τ = 1000 s at ~ 87° elongation (4.44 × 10^-14 at τ = 10 s).

More information regarding the origin of the phase fluctuations can be derived by analyzing the spatial statistics of the phase scintillation in multiple stations during the same tracking session. If the φ_Δ phase data of a few pairs of widely spaced stations are cross-correlated, as suggested in Armstrong (1998), it could be determined whether the main contributor to the phase fluctuations is the interplanetary medium or the local impact of the ionosphere at each station. Unfortunately, in the GR035 experiment this analysis could not be performed, since although the three stations operating in Doppler mode are widely spaced (South Africa, Finland, and China), the differential phase φ_Δ could not be successfully retrieved because of the high RFI on the S band of Sh and technical problems with the X-band receiver of Mh (for this reason the values for Mh are not shown in Table 1).

3.3. Noise budget for the Doppler detections of GR035

Table 2 summarizes the Allan deviations found for the different noise sources described in Sects. 3.1 and 3.2. The Allan deviations from the thermal noise of the ground station for the VLBI stations vary from σ_y(τ) = 0.46−1.60 × 10^-14 at τ = 10 s. Despite the differences of the thermal noise σ_y between the VLBI stations and the DSN and NNO stations⁶, the thermal noise of the stations does not dominate the error budget of the observations in this experiment, as shown in Table 2. Because of the long tracking sessions of the antennas during this experiment, the antenna elevations have a higher impact in the S/N of the detections than the collecting area of the antennas. This issue is usually ignored in shorter tracking sessions, since only stations with elevations ~>20 deg are selected to participate in an observation.

The plasma scintillation noise was estimated for Ht, which was one of the stations observing both in S and X band. The plasma scintillation noise is more dominant for Ht than its thermal noise (σ_y = 4.44 × 10^-14 against σ_y = 1.0 × 10^-14 for τ = 10 s). Owing to problems with the receivers, this analysis could not be performed for the other two stations receiving the dual-band link. In future experiments, the contributions of the ionosphere and interplanetary medium could be discerned from one another by correlating the power spectra of the differential phases between every pair of stations.

Armstrong et al. (2008) reported that for the DSN stations, when the propagation noises are properly calibrated, the antenna mechanical noise was the leading noise of their noise budget. Regarding the VLBI stations, Sarti et al. (2009, 2011) have reported one-way path delay variations caused by antenna mechanical noise, however these were computed for VLBI geodetic and astrometric studies, for which the delay stability is evaluated in annual timescales, which is much larger than the integration times relevant for the study at hand. Nonetheless, owing to their size (Armstrong 2016), the expected mechanical noise of the VLBI antennas (except for Tm65) will be considerably less than the 70 m DSN antennas. In fact, simultaneous observations between PRIDE and DSN stations could help improve the sensitivity of the 70 m DSN antennas. Following the approach presented in Armstrong et al. (2008) in future experiments, stations of the global VLBI network close to the deep space tracking complexes could be used to remove the antenna mechanical noise of the larger antennas during simultaneous two- and three-way Doppler passes, for instance, the 2-m VLBA-Ov close the DSS-14, the 14 m Ys close to the DSS-63, the 12 m Ye telescope close to NNO, and the 12 m Atacama Pathfinder Experiment (APEX) telescope close to Estrack’s Malargüe station.

4. Conclusions

With the PRIDE setup, Doppler tracking of the spacecraft carrier signal with several Earth-based radio telescopes is performed, subsequently correlating the signals coming from the different telescopes in a VLBI-style. Although the main output of this technique are VLBI observables, we demonstrated that the residual frequencies obtained from the open-loop Doppler observables-which are inherently derived in the data processing pipeline to retrieve the VLBI observables-are comparable to those obtained with the closed-loop Doppler data from NNO, DSS-63, and DSS-14 stations (see Fig. 2). Figure 2 shows the best case found, where the median value of the residuals fit achieved with VLBI station Kp remains within 1 mHz of the residuals fit obtained with DSS-63 and DSS-14. The median of the Doppler residuals for all the detections with the VLBI stations was found to be ~2 mHz.

The fact that this experiment involved long tracking sessions makes the variability of the elevation angle of the antennas a factor in the characterization of the noise that cannot be ignored. At elevations <20°, noise contributions due to larger tropospheric path delays and larger spillover noise have a larger impact on the total system temperature compared to those of the receiver temperature. For this reason, there are cases for which antennas with smaller collecting areas reach similar S/N levels as larger antenna dishes, as shown in Fig. 3, due to a more favorable antenna elevation. The derived Allan deviations due to thermal noise at the VLBI stations vary between σ_y(τ) = 0.46−1.60 × 10^-14 at τ = 10 s. For this particular experiment, at the DSN stations the expected σ_y(τ) from thermal noise was 6.5 × 10^-15 at τ = 10 s. Although only four of the VLBI stations have comparable Allan deviations (Table 1) to those of the DSN stations, the thermal noise is not the most dominant contribution to the overall noise budget of this experiment.

Although they were not included in this particular experiment (other than the 65 m Tianma station), PRIDE has access through the EVN to multiple radio telescopes that are similiar in size or larger than the DSN antennas; these include the 64 m Sardinia, 100 m Effelsberg, and 305 m Arecibo, which can be scheduled for radio science experiments. The use of these large antennas can result in an advantage when conducting experiments with limited S/N, such as radio occultation experiments of planets and moons with thick atmospheres.

Open-loop Doppler data, such as those collected with PRIDE experiments, present advantages for certain radio science applications compared to closed-loop data. However, closed-loop Doppler tracking is routinely performed in the framework of navigation tracking and does not require post-processing to retrieve the Doppler observables. Although the Estrack/DSN complexes have the capability of simultaneously gathering closed-loop and open-loop Doppler data, this is not an operational mode required for navigation nor telemetry passes, which generally operate in closed-loop mode only. In this sense, PRIDE Doppler data could complement the closed-loop tracking data and enhance the science return of tracking passes that are not initially designed for radio science experiments.

Acknowledgments

The EVN is a joint facility of European, Chinese, South African, and other radio astronomy institutes funded by their national research councils. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Australia Telescope Compact Array is part of the Australia Telescope National Facility, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. T. Bocanegra Bahamon acknowledges the NWO–ShAO agreement on collaboration in VLBI. G. Cimó acknowledges the EC FP7 project ESPaCE (grant agreement 263466). P. Rosenblatt is financially supported by the Belgian PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office. We express gratitude to M. Pätzold (MEX MaRS PI) and B. Häusler for coordination of MaRS and PRIDE tracking during the MEX/Phobos flyby and a number of valuable comments on the manuscript of the current paper. Mars Express is a mission of the European Space Agency. The MEX a priori orbit, Estrack and DSN tracking stations transmission frequencies, and the cyclogram of events were supplied by the Mars Express project. The authors would like to thank the personnel of the participating stations. R.M. Campbell, A. Keimpema, P. Boven (JIVE), O. Witasse (ESA/ESTEC), and D. Titov (ESA/ESTEC) provided important support to various components of the project. The authors are grateful to the anonymous referee for useful comments and suggestions.

References

All Tables

All Figures

Fig. 1 Example of Doppler data processing pipeline using observations of MEX during Phobos flyby by Hartebeesthoek (see Sect. 3). Panel a shows the typical resulting average power spectrum of a scan after running the SWSpec software. In the 16 MHz pre-defined sub-band starting at 8412 MHz, the spacecraft signal is found in the spectrum. A narrowband containing the moving phase of the spacecraft carrier/tone is selected (in this case of 50 kHz bandwidth) to model the Doppler shift using an n-order polynomial frequency fit. Panel b shows a zoom of the spectrum inside the selected narrowband window to perform the fit. Here the moving phase of the carrier tone is visible along the duration of the scan. After the fit is performed, SCtracker applies the polynomial coefficients after converting the sample to the baseband sample to stop the moving phase of the tone. In this way, SCtracker extracts an initial fit of the Doppler shift. Panel c shows the output of the SCtracker, which is a phase-stopped filtered out signal in a 2 kHz narrowband in baseband. Finally, the dPLL performs high-precision iterations of the time integration of the overlapped spectra, phase polynomial fitting, conversion to baseband, and phase-stopping corrections, using narrowband windows around the carrier tone. The iterations stop when the window bandwidth reaches 20 Hz, as shown in panel d, allowing the extraction of the frequency and phase residuals of the spacecraft carrier tone with a 2 mHz frequency spectral resolution.

In the text

	Fig. 2 Comparison of the Doppler residuals obtained with VLBA-Kp (in black), DSS-63 (in blue), and DSS-14 (in red). The median value of the difference between the fit of VLBA-Kp and the fit of DSS-63 and DSS-14, respectively, remains below 1 mHz, for an integration time of 10 s.
In the text

Fig. 3 Signal-to-noise ratio and elevation angle comparison between 25 m Ur, 25 m On, and 15 m Ht radio telescopes from 4:00–7:30 UTC on December 29, 2013 (TX: DSS-63). Left panel shows that although Ht (in pink) has a lower collecting area, it achieves higher S/N levels than Ur (in blue) and similar S/N levels than On (in green). Right panel shows the elevation angle of each antenna during the same time period from which the correlation with the S/N levels of each station is evident.

In the text

	Fig. 4 Spectral power density of MEX signal in S band. The scintillation band extends from 8 mHz to 0.1 Hz, obtaining a value for the slope of –2.471, which is coherent with the spectral index values found by (Woo & Armstrong 1979). The mean phase scintillation index of the signal received from 3h56m to 09h12m (UTC) on 2013-12-29, is 0.070 rad, at an elongation of ~ 87° and distance of ~ 1.4 AU.
In the text

	Fig. 5 Spectral power density of MEX signal in X band. The scintillation band extends from 8 mHz to 0.45 Hz, obtaining a value for the slope of –2.469. The mean phase scintillation index of the signal received from 3h56m to 09h12m (UTC) on 2013-12-29 is 0.073 rad at an elongation of ~ 87° and distance of ~ 1.4 AU.
In the text

	Fig. 6 Spectral power density of the differential phases φΔ. The scintillation band extends from 8 mHz to 0.1 Hz, obtaining a value for the slope of –2.372. The mean phase scintillation index of the signal received from 3h56m to 09h12m (UTC) on 2013-12-29 is 0.069 rad at an elongation of ~ 87° and distance of ~ 1.4 AU.
In the text