© ESO, 2012

1. Introduction

The VLBI spacecraft tracking technique has been successfully applied to a number of deep-space missions. The technique imposes minimal requirements on the mission’s on-board instrumentation and can operate with almost any radio signal from a spacecraft, provided it is sufficiently powerful and phase-stable. Among VLBI-tracked missions one can mention the VEGA Venus atmosphere balloons (Preston et al. 1986), the Huygens Probe during its descent to the surface of Saturn’s moon Titan (Lebreton et al. 2005), VLBI tracking with the European VLBI Network (EVN) antennas of the controlled impact of ESA’s Smart-1 probe on the surface of the Moon (Pogrebenko et al. 2006), VLBI tracking of NASA’s Mars Exploration Rover B spacecraft during its final cruise phase (Lanyi et al. 2007), VLBI tracking of the Cassini spacecraft at Saturn (Jones et al. 2011), and the latest results achieved with ESA’s Venus EXpress (VEX) VLBI spacecraft observations and the ESA Mars EXpress (MEX) Phobos-flyby (Molera Calvés et al. 2010) monitored by the EVN radio telescopes.

The Planetary Radio Interferometry and Doppler Experiment (PRIDE) is an initiative by the Joint Institute for VLBI in Europe (JIVE). Designed as a multi-purpose, multi-disciplinary enhancement of planetary missions science return, PRIDE is able to provide ultra-precise estimates of spacecraft state vectors based on the phase-referenced VLBI tracking and radial Doppler measurements. These can be used for a variety of scientific applications including planetary science (measurements of tidal deformations of planetary moons, atmosphere dynamics and climatology as well as seismology, tectonics, internal structure and composition of planetary bodies), ultra-precise celestial mechanics of planetary systems, gravimetry and fundamental physics (e.g. tests of general relativity or study of the anomalous accelerations of deep-space probes, Turyshev & Toth 2010). Spacecraft Doppler tracking is known to be the only possible way of detecting low-frequency (10^-5−1 Hz) gravitational waves (Kopeikin & Schaefer 1999). Studies of interplanetary scintillations (Molera Calvés et al., in prep.) are among the potentially rewarding by-products of PRIDE. PRIDE can also provide diagnostics of deep-space missions (“health check”) and direct-to-Earth delivery of a limited amount of critical data (e.g. Friedman et al. 2008).

PRIDE can be conducted with literally any radio-emitting spacecraft. Current PRIDE “customers” include ESA’s VEX and MEX spacecraft and the Russian Space Agency’s Space VLBI mission RadioAstron (Kardashev 2009). PRIDE has been considered in a number of design studies, such as EVE (Wilson et al. 2011), JUICE/Laplace (Blanc et al. 2009; Christophe et al. 2009), TandEM (Coustenis et al. 2009), Kronos (Marty et al. 2009), and several others.

In this paper, we describe the principles and pipeline of the data acquisition, processing and analysis within the scope of PRIDE and present the results of the first large-scale observational EVN campaign conducted in 2009–2011 with ESA’s VEX spacecraft.

2. Observations and data processing pipeline

In phase-referencing VLBI, the relative position of a specific source is determined by observing a known source nearby (i.e. at a short angular distance on the celestial sphere), called a phase-reference calibrator, and applying calibration values of its fringe phase to the fringe phases of a target source. The most accurate measurements can be achieved when both the target and calibrator sources are within the same primary beam of a VLBI antenna or no more than a few degrees apart (Ros et al. 1999). VLBI phase referencing for spacecraft tracking involves alternated observations of the radio signal from the spacecraft and a nearby calibrator, most commonly a quasar. The best switching time between the calibrator and target sources for out-of-beam phase-referencing ranges from tens of seconds to several minutes. In the out-of-beam regime, the ideal scan duration on each source depends primarily on the efficiency of the radio telescopes, their system temperature, and the slewing speed. In the observations described below, the scan duration was set to 2 − 5 min.

Spacecraft VLBI-tracking data processing and analysis tools are based on two major software engines: the ultra-high spectral resolution and phase detection software package SWSpec/SCTrack/phase-lock-loop, developed at the Metsähovi Radio Observatory in collaboration with JIVE¹, and the EVN Software Correlator at JIVE SFXC (Keimpema et al. 2011).

Successful detection, processing, and analysis of spacecraft VLBI-tracking data requires usage of a near-field delay model and computation of near-field analogues of baseline uv-projections. To address these tasks, a software package VINT (VLBI in the near-filed toolkit) was developed at JIVE and is described in detail below.

2.1. Observations

Preparation for observations, scheduling. In standard VLBI observations, coordinates of a source are considered to be constant during the observing run for all participating telescopes. Usually, the source position in the geocentric equatorial inertial coordinate system J2000.0 is used. For near-field VLBI observations of spacecraft, the target can move rapidly across the primary beam of a telescope. Therefore, the position of a spacecraft is calculated at different epochs and is then used for pointing the antenna at those particular epochs. In the “ultra” near-field case, i.e. when the target is within (or very close to) the synthesised aperture of a VLBI array², so-called “cross-eyed” scheduling should be used. This means that the spacecraft J2000.0 positions are different in the schedules for different VLBI stations.

Raw data. The down-converted signals observed by radio telescopes are recorded with the standard VLBI recording systems Mark5 A/B (Whitney 1999), or the PC-EVN (Mujunen & Ritakari 2004). Both systems were designed exclusively to record formatted data of astronomic and geodetic VLBI observations, allowing high data-rate ranging up to 2 Gbps.

Fig. 1 Data processing pipeline of phase-referencing VLBI observations of spacecraft.

Figure 1 represents the data processing pipeline that is described in detail below.

2.2. Broadband correlation of reference-source signal and S/C data band

Initial processing. Once at JIVE, the raw data are correlated with the JIVE software correlator SFXC (Keimpema et al. 2011), which has been developed at JIVE and is currently used as a production correlator for some EVN observations. The SFXC is based on the original design developed for VLBI tracking of the Huygens Probe (Pogrebenko et al. 2008) and is capable of supporting both the far-field and near-field delay models and processing data in various “flavours” of the Mk5 format.

When correlated with the SFXC, the fringe finder³ data are used for the clock search and bandpass correction. Also, a two-way Doppler phase correction (obtained with the SCTrack software, see next subsection) is applied to the spacecraft data to avoid frequency smearing.

The correlation results are saved in a measurement set format (Kemball & Wieringa 2000). After this, the amplitudes and phases of the cross-correlation function (CCF) are extracted for further processing.

Group delay extraction. The complex CCF corrected for clock/clock-rate offsets are smoothed in the frequency domain. After this, the phase is extracted for each averaged spectrum and fitted by a linear function. A group delay for the centre of each averaging interval is derived as the slope of the function τ_gr = dφ_fit/df, where f denotes frequency.

2.3. Narrow-band correlation of S/C signal

SWSpec, SCTrack and phase-lock loop. Narrow band correlation of the spacecraft signal is performed with the SWSpec/SCTrack software package described in Molera Calvés et al. (2011).

The high-resolution spectrometer software (SWspec) allows the initial detection of the carrier tone of the spacecraft and temporal evolution of its frequency over the entire scan. In each SWspec pass, one of the raw data channels is extracted from the input data file. The software performs an accurate window-overlapped discrete Fourier transform (DFT) and time integration of the spectrum. The generated time-integrated spectra are processed to determine the moving phase of the S/C tone frequencies along the scan. The series of frequency detections of the carrier line are fitted to a sixth-order phase-stopping polynomial.

The core of the spacecraft tone tracking software (SCtracker) performs a phase polynomial correction to compensate for the Doppler shift, single and/or multi-tone tracking, signal filtering for each tone and, finally, phase detection of the carrier line with respect to the local H-maser clock.

The final post-analysis of the narrow-bands is performed with the phase-lock loop (PLL) software. The PLL runs high-precision iterations of the phase polynomial fit on the filtered complex narrow-band signals. The implementation of the PLL allows filtering down the tones as often as desired to achieve mHz accuracy around the spacecraft carrier line or other tones in the band. The final residual phase in a specific stopped band is determined with respect to the set of subsequent frequency/phase polynomials initially applied.

Reduction of phases to a common phase centre. Topocentric phase detections for each station are reduced to the common phase centre (usually the geocentre) by applying the geocentric signal delay, τ_gc, for the near-field case (detailed description of the model used is given in the next subsection): ${\mathit{\phi }}_{\mathrm{gc}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{\phi }}_{\mathrm{tc}}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{\tau }}_{\mathrm{gc}}\mathrm{)}\mathit{.}$(1)These geocentric phases are used as Doppler phase corrections for the broad-band correlation with SFXC. Finally, for each baseline i-j, differential phases φ_ij = φ_j,gc − φ_i,gc are formed.

2.4. Resolution of the 2π-ambiguity

A straightforward correlation of the wide-band spacecraft signal allows deriving the group delay τ_gr = dφ/df, while the narrow-band correlation of the S/C carrier tone allows deriving the phase delay τ_ph = φ/f₀, where f₀ is the observational frequency. These are two independent measurables, which should give the same phase solution. The rms error of the group delay is σ(τ_gr) = 1/(ΔF × S/N), with ΔF being the bandwidth of, e.g. 8 or 16 MHz for the calibrator and about 1 MHz for the S/C but with an uneven distribution of signal-to-noise ratio (S/N) over the signal band. The rms error of the phase delay for a narrow-band signal is given by $\mathit{\sigma }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{ph}}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{(}{\mathit{f}}_{\mathrm{0}}\mathrm{\times }\sqrt{\mathit{S}\mathit{/}\mathit{N}}\mathrm{)}$, where the S/N level for the spectral line is very high, so that even the $\sqrt{\mathit{S}\mathit{/}\mathit{N}}$ may be larger than the S/N for the broad band case. In addition to this, f₀ ≫ ΔF, so that the rms error of group delay is significantly larger than that of the phase delay. This in turn results in a different scatter of the phase solution resulting from group and phase delays.

The group delay, although not that accurate, is a “true delay”, while the phase delay is derived from the measured phase and may be subject to a 2π ambiguity. Since phases are 2π uncertain, the phase delay may have a bias of several cycles of  ~ 120 ps for observations at X-band. Checking the phase delay versus group delay consistency by computing the uncertainty of the fit of the phase delay to the group delay and shifting the former by n·2π (where n is an integer number) to obtain the best fit to the latter allows one to resolve the 2π ambiguity.

2.5. Delay model

The VLBI delay model is formulated in the barycentric celestial reference system (BCRS). This means that the delay is computed in the time scale of the BCRS and is subsequently transformed to the time-scale used for timing the measurements on Earth. Thus several time-scale and station coordinate transformations are required.

The models described below are coded in Matlab and aggregated within the VINT software package.

Geocentric terrestrial to celestial frame transformation. As a first step, a transformation of station coordinates from the terrestrial reference system (realised by the international terrestrial reference frame, ITRF) to the geocentric celestial reference system (GCRS) is conducted according to the latest recommendations of the international Earth rotation service IERS 2010 conventions (Petit & Luzum 2010) within the framework of the IAU 2000/2006 resolutions.

Station positions and velocities in the ITRF at the epoch t₀ (e.g. 1997.0 or 2000.0, depending on the realisation used), are reduced to the epoch of observation, taking the tectonic plate motion into account. A simple linear model is used: ${{r}}_{\mathrm{ITRF}}\mathrm{=}{{r}}_{\mathrm{0}}\mathrm{+}{\dot{{r}}}_{\mathrm{0}}\mathrm{·}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(2)At the next step, the station position, velocity, and acceleration in the GCRS are computed as follows: $\begin{array}{ccc}{{r}}_{\mathrm{GCRS}}\mathrm{=}& & \mathit{T}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}{{r}}_{\mathrm{ITRF}}\mathit{,}\\ & & \mathrm{(}{{r}}_{\mathrm{GCRS}}\mathrm{\to }{\dot{{r}}}_{\mathrm{GCRS}}\mathit{,}{\begin{array}{}\mathrm{¨}\\ {r}\end{array}}_{\mathrm{GCRS}}\mathrm{;}\mathit{T}\mathrm{(}\mathit{t}\mathrm{)}\mathrm{\to }\mathit{Ṫ}\mathrm{(}\mathit{t}\mathrm{\left)}\mathit{,}\begin{array}{}\mathrm{¨}\\ \mathit{T}\end{array}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$(3)where the transformation matrix T(t) = Q(t)·R(t)·W(t) with Q(t), R(t) and W(t) being 3 × 3 matrices representing the motion of the celestial pole in the GCRS, rotation of the Earth around the axis associated with the pole and polar motion, correspondingly, and the derivatives $\mathit{Ṫ}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}\begin{array}{}{}_{\mathrm{¨}}\\ \mathit{T}\end{array}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ are calculated numerically.

Displacements of station reference points. Station positions in the GCRS should be subsequently corrected for various geophysical effects. These include displacements caused by the solid Earth tides arising from the direct effect of the external tide-generating potential, ocean tidal loading, diurnal and semidiurnal atmospheric pressure loading, and the centrifugal perturbations caused by the pole tide-driven Earth rotation variations. These effects are taken into account following the IERS conventions 2010 (Petit & Luzum 2010).

Geocentric to barycentric transformation. Finally, the Lorentz transformation needs to be applied to the station GCRS position r_gc, velocity ${\dot{{r}}}_{\mathrm{gc}}$ and acceleration ${\begin{array}{}{\mathrm{¨}}_{}\\ {r}\end{array}}_{\mathrm{gc}}$ to transform it in the solar system barycentric reference frame BCRS: ${{r}}_{\mathrm{bc}}\mathrm{=}\left(\mathrm{1}\mathrm{-}{\mathit{L}}_{\mathrm{C}}\mathrm{-}\frac{\mathit{\gamma }{\mathit{U}}_{\mathrm{E}}}{{\mathit{c}}^{\mathrm{2}}}\right)\mathrm{·}{{r}}_{\mathrm{gc}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}^{\mathrm{2}}}\mathrm{(}{{V}}_{\mathrm{E}}\mathrm{·}{{r}}_{\mathrm{gc}}\mathrm{)}\mathrm{·}{{V}}_{\mathrm{E}}\mathrm{+}{{R}}_{\mathrm{E}}$(4)$\begin{array}{ccc}{\dot{{r}}}_{\mathrm{bc}}\mathrm{=}& & \left(\mathrm{1}\mathrm{-}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}{\mathit{U}}_{\mathrm{E}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{V}}_{\mathrm{E}}^{\mathrm{2}}}{\mathrm{2}{\mathit{c}}^{\mathrm{2}}}\mathrm{-}\frac{{{V}}_{\mathrm{E}}\mathrm{·}{\dot{{r}}}_{\mathrm{gc}}}{{\mathit{c}}^{\mathrm{2}}}\right)\mathrm{·}{\dot{{r}}}_{\mathrm{gc}}\\ & & \mathrm{+}{{V}}_{\mathrm{E}}\mathrm{·}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}^{\mathrm{2}}}\mathrm{(}{{V}}_{\mathrm{E}}\mathrm{·}{\dot{{r}}}_{\mathrm{gc}}\mathrm{)}\right)\end{array}$(5)$\begin{array}{ccc}{\begin{array}{}\mathrm{¨}\\ {r}\end{array}}_{\mathrm{bc}}\mathrm{=}& & \left(\mathrm{1}\mathrm{+}{\mathit{L}}_{\mathrm{C}}\mathrm{-}\frac{\mathrm{(}\mathrm{2}\mathrm{+}\mathit{\gamma }\mathrm{)}{\mathit{U}}_{\mathrm{E}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{V}}_{\mathrm{E}}^{\mathrm{2}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathrm{2}{V}}_{\mathrm{E}}\mathrm{·}{\dot{{r}}}_{\mathrm{gc}}}{{\mathit{c}}^{\mathrm{2}}}\right)\mathrm{·}{\begin{array}{}\mathrm{¨}\\ {r}\end{array}}_{\mathrm{gc}}\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}^{\mathrm{2}}}\mathrm{(}{{V}}_{\mathrm{E}}\mathrm{·}{\begin{array}{}\mathrm{¨}\\ {r}\end{array}}_{\mathrm{gc}}\mathrm{)}\mathrm{·}\mathrm{(}{{V}}_{\mathrm{E}}\mathrm{+}\mathrm{2}{\dot{{r}}}_{\mathrm{gc}}\mathrm{)}\mathrm{+}{{A}}_{\mathrm{E}}\mathit{,}\end{array}$(6)where L_C = 1.48082686741 × 10^-8,   ${\mathit{U}}_{\mathrm{E}}\mathrm{=}\sum _{\mathit{j}\ne \mathit{i}}\frac{\mathit{G}{\mathit{M}}_{\mathit{j}}}{{\mathit{r}}_{\mathrm{E}\mathit{j}}}$ is the Newtonian potential of all solar system bodies excluding the Earth, evaluated at the geocenter, γ is the PPN parameter (equal to 1 in general relativity), R_E, V_E and A_E are the position, velocity and acceleration of the Earth in the BCRS, respectively, and c is the speed of light in a vacuum.

Time-scale transformations. The time-scale used for measurement timing at stations is the coordinated universal time (UTC), which differs from the physically realised international atomic time scale, TAI, by a number of leap-seconds (currently, 34 s). The latter, in turn, differs by 32.184 s from the terrestrial time (TT), which is the theoretical time-scale for clocks at sea-level. The TT is a scaled version of the geocentric coordinate time (TCG, the time coordinate of the IAU space-time metric GCRS), which eliminates the combined effect on the terrestrial clock of the gravitational potential from the Earth and the observatory’s diurnal speed in TCG (Petit & Luzum 2010): $\mathrm{TCG}\mathrm{=}\mathit{TT}\mathrm{+}{\mathit{L}}_{\mathrm{G}}\mathrm{·}\mathrm{(}\mathit{J}{\mathit{D}}_{\mathrm{TT}}\mathrm{-}\mathit{T}{\mathit{T}}_{\mathrm{0}}\mathrm{)}\mathit{,}$(7)where TT₀ is TT at 1977 January 1.0 TAI, JD_TT is TT, both expressed as Julian date, and L_G = 6.969290134 × 10^-10.

The time-scale used in the ephemerides of planetary spacecraft, as well as that of solar system bodies, is the barycentric dynamical time (TDB), a scaled version of the barycentric coordinate time (TCB)  −  the time coordinate of the IAU space-time metric BCRS. The TDB stays close to TT on the average by suppressing a drift in TCB due to the combined effect of the terrestrial observer’s orbital speed and the gravitational potential from the Sun and planets by applying a linear transformation (Petit & Luzum 2010) $\mathrm{TDB}\mathrm{=}\mathrm{TCB}\mathrm{-}{\mathit{L}}_{\mathrm{B}}\mathrm{·}\mathrm{(}\mathit{J}{\mathit{D}}_{\mathrm{TCB}}\mathrm{-}{\mathit{T}}_{\mathrm{0}}\mathrm{)}\mathrm{·}\mathrm{86400}\mathrm{+}\mathit{TD}{\mathit{B}}_{\mathrm{0}}\mathit{,}$(8)where JD_TCB is TCB expressed as Julian date, T₀ = 2443144.5003725, L_B = 1.550519768 × 10^-8 and TDB₀ = 6.55 × 10^-5 s.

The difference between TCG and TCB involves a transformation accounting for the orbital speed of the geocenter and the gravitational potential from the Sun and planets with the accuracy of O(c^-4) given by $\mathrm{TCB}\mathrm{-}\mathrm{TCG}\mathrm{=}{\mathit{c}}^{-2}\left({\mathrm{\int }}_{{\mathrm{t}}_{\mathrm{0}}}^{\mathit{t}}\left[\frac{{\mathit{V}}_{\mathrm{E}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}{\mathit{U}}_{\mathrm{E}}\mathrm{\left(}{{R}}_{\mathrm{E}}\mathrm{\right)}\right]\mathrm{d}\mathit{t}\mathrm{+}{{V}}_{\mathrm{E}}\mathrm{·}\mathrm{(}{R}\mathrm{-}{{R}}_{\mathrm{E}}\mathrm{)}\right)\mathit{,}$(9)where R_E denotes the BCRS position of the Earth and R is the BCRS position of the observing station.

Computation of the difference TCB  −  TCG is performed on the basis of the Fairhead-Bretagnon model (Fairhead & Bretagnon 1990) in its full form, which yields on accuracy of a few nanoseconds. The full time–scale transformation chain is as follows: UTC → TAI → TT → TCG → TCB → TDB.

VLBI delay model for extragalactic sources. For calculating the VLBI signal delay for the extragalactic sources, we use the so-called consensus model recommended by the IERS (Petit & Luzum 2010), which is based on the plane-wave approximation when the distance to a source is assumed to be infinite. In the BCRS, the vacuum delay equation is ${\mathit{T}}_{\mathrm{2}}\mathrm{-}{\mathit{T}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{{K}}{\mathit{c}}\mathrm{·}\mathrm{\left(}{{R}}_{\mathrm{2}}\mathrm{\right(}{\mathit{T}}_{\mathrm{2}}\mathrm{)}\mathrm{-}{{R}}_{\mathrm{1}}\mathrm{(}{\mathit{T}}_{\mathrm{1}}\mathrm{\left)}\mathrm{\right)}\mathrm{+}\mathrm{\Delta }{\mathit{T}}_{\mathrm{grav}}\mathit{.}$(10)The total barycentric vacuum delay is then converted into the geocentric one by applying the Lorentz transformations and is given by $\begin{array}{ccc}{\mathit{t}}_{\mathrm{2}}\mathrm{-}{\mathit{t}}_{\mathrm{1}}\mathrm{=}& & (\mathrm{1}\mathrm{+}\frac{{K}\mathrm{·}\mathrm{(}{{V}}_{\mathrm{E}}\mathrm{+}{\dot{{r}}}_{\mathrm{2}\mathit{,}\mathrm{gc}}\mathrm{)}}{\mathit{c}}{)}^{-1}\mathrm{·}(\mathrm{\Delta }{\mathit{T}}_{\mathrm{grav}}\mathrm{-}\frac{{K}\mathrm{·}{b}}{\mathit{c}}\\ & & \mathrm{·}[\mathrm{1}\mathrm{-}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}{\mathit{U}}_{\mathrm{E}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{-}\frac{{{V}}_{\mathrm{E}}^{\mathrm{2}}}{\mathrm{2}{\mathit{c}}^{\mathrm{2}}}\mathrm{-}\frac{{{V}}_{\mathrm{E}}\mathrm{·}{\dot{{r}}}_{\mathrm{2}\mathit{,}\mathrm{gc}}}{{\mathit{c}}^{\mathrm{2}}}]\\ & & \mathrm{-}\frac{{{V}}_{\mathrm{E}}\mathrm{·}{b}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{·}\left(\mathrm{1}\mathrm{+}\frac{{K}\mathrm{·}{{V}}_{\mathrm{E}}}{\mathrm{2}\mathit{c}}\right))\mathit{.}\end{array}$(11)Here K is the unit vector from the barycentre to the source in the absence of gravitational or aberrational bending, b = r₂(t₁) − r₁(t₁) is the GCRS baseline vector at the time of signal arrival at the first station. The general relativistic delay, ΔT_grav, is given for the J_th gravitating body by $\mathrm{\Delta }{\mathit{T}}_{\mathrm{grav}\mathit{,J}}\mathrm{=}\mathrm{2}\frac{\mathit{G}{\mathit{M}}_{\mathit{J}}}{{\mathit{c}}^{\mathrm{3}}}\ln \frac{{\mathit{R}}_{\mathrm{1}\mathit{J}}\mathrm{+}{K}\mathrm{·}{{R}}_{\mathrm{1}\mathit{J}}}{{\mathit{R}}_{\mathrm{2}\mathit{J}}\mathrm{+}{K}\mathrm{·}{{R}}_{\mathrm{2}\mathit{J}}}$(12)with R_1J = R_1J(t₁) = R₁(t₁) − R_J(t_1J) and ${{R}}_{\mathrm{2}\mathit{J}}\mathrm{=}{{R}}_{\mathrm{2}\mathit{J}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}{{R}}_{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}\frac{{{V}}_{\mathrm{E}}}{\mathit{c}}\mathrm{(}{K}\mathrm{·}{b}\mathrm{)}\mathrm{-}{{R}}_{\mathit{J}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}\mathit{J}}\mathrm{\right)}$, where R₁ and R₂ are the barycentric radius-vectors of the first and second receivers, t₁ is the TCG time of arrival of the radio signal at the first VLBI receiver, t_1J  −  the moment of closest approach of the photon to the gravitating body  J in TCG. The total gravitational delay is the sum over all gravitating bodies including the Earth: ΔT_grav =  ∑ JΔT_grav,J .

Fig. 2 Geometry of VLBI observations of spacecraft in the barycentric celestial reference frame.

Delay model for VLBI observations of spacecraft. For spacecraft VLBI observations, a delay model should be able to account for the “near-field” effects caused by finiteness of the source’s distance. Indeed, the diffracted radio wave is considered to be in the near field if the distance to the source from the aperture R ≲ D²/λ (Born & Wolf 2002), where D is the characteristic size of the aperture and λ is the wavelength. For a VLBI array with a synthesised aperture of  ~1000 km observing at X-band (λ ~ 3 cm) this yields R ≲ 100 AU. Note that near-field delays can differ from “classic” ones by microseconds.

The signal delay in the TDB-frame is defined as the difference between two light travel times LT₁ and LT₂ from the spacecraft to the first and the second station of the baseline (Fig. 2). Given the reception time T₁ in TDB at the first station, the transmission time T₀ is calculated by iteratively solving the light-time equation for the “down-leg” 0-1: ${\mathit{T}}_{\mathrm{1}}\mathrm{-}{\mathit{T}}_{\mathrm{0}}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{01}}}{\mathit{c}}\mathrm{+}\mathit{RL}{\mathit{T}}_{\mathrm{01}}\mathit{,}$(13)where RLT₀₁ is the relativistic term, which accounts for the effects of both special and general relativity (Moyer 2003): $\begin{array}{ccc}\mathit{RL}{\mathit{T}}_{\mathrm{01}}\mathrm{=}& & \frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}\mathrm{·}\mathit{G}{\mathit{M}}_{\mathrm{S}}}{{\mathit{c}}^{\mathrm{3}}}\mathrm{·}\ln \frac{{\mathit{R}}_{\mathrm{0}}^{\mathrm{S}}\mathrm{+}{\mathit{R}}_{\mathrm{1}}^{\mathrm{S}}\mathrm{+}{\mathit{R}}_{\mathrm{01}}^{\mathrm{S}}\mathrm{+}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}\mathrm{·}\mathit{G}{\mathit{M}}_{\mathrm{S}}}{{\mathit{c}}^{\mathrm{2}}}}{{\mathit{R}}_{\mathrm{0}}^{\mathrm{S}}\mathrm{+}{\mathit{R}}_{\mathrm{1}}^{\mathrm{S}}\mathrm{-}{\mathit{R}}_{\mathrm{01}}^{\mathrm{S}}\mathrm{+}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}\mathrm{·}\mathit{G}{\mathit{M}}_{\mathrm{S}}}{{\mathit{c}}^{\mathrm{2}}}}\\ & & \mathrm{+}\sum _{\mathit{B}\mathrm{=}\mathrm{1}}^{\mathrm{10}}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}\mathrm{·}\mathit{G}{\mathit{M}}_{\mathit{B}}}{{\mathit{c}}^{\mathrm{3}}}\mathrm{·}\ln \frac{{\mathit{R}}_{\mathrm{0}}^{\mathit{B}}\mathrm{+}{\mathit{R}}_{\mathrm{1}}^{\mathit{B}}\mathrm{+}{\mathit{R}}_{\mathrm{01}}^{\mathit{B}}}{{\mathit{R}}_{\mathrm{0}}^{\mathit{B}}\mathrm{+}{\mathit{R}}_{\mathrm{1}}^{\mathit{B}}\mathrm{-}{\mathit{R}}_{\mathrm{01}}^{\mathit{B}}}\mathrm{·}\end{array}$(14)Here S denotes the Sun and B runs over all solar system planets and the Moon, R₀ is the spacecraft barycentric position⁴, $\begin{array}{ccc}{{R}}_{\mathrm{01}}& \mathrm{=}& {{R}}_{\mathrm{1}}\mathrm{\left(}{\mathit{T}}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}{{R}}_{\mathrm{0}}\mathrm{\left(}{\mathit{T}}_{\mathrm{0}}\mathrm{\right)}\\ {{R}}_{\mathit{i}}^{\mathit{\alpha }}& \mathrm{=}& {{R}}_{\mathit{i}}\mathrm{\left(}{\mathit{T}}_{\mathit{i}}\mathrm{\right)}\mathrm{-}{{R}}_{\mathit{\alpha }}\mathrm{\left(}{\mathit{T}}_{\mathit{i}}\mathrm{\right)}\mathit{,}\mathit{i}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{1}\mathrm{;}\mathit{\alpha }\mathrm{=}\mathit{S,B}\\ {{R}}_{\mathrm{01}}^{\mathit{\alpha }}& \mathrm{=}& {{R}}_{\mathrm{1}}^{\mathit{\alpha }}\mathrm{\left(}{\mathit{T}}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}{{R}}_{\mathrm{0}}^{\mathit{\alpha }}\mathrm{\left(}{\mathit{T}}_{\mathrm{0}}\mathrm{\right)}\mathit{,}\mathit{\alpha }\mathrm{=}\mathit{S,B}\mathit{.}\end{array}$At each iteration step the correction to the T₀ estimate is given by $\mathrm{\Delta }{\mathit{T}}_{\mathrm{0}}\mathrm{=}\frac{{\mathit{T}}_{\mathrm{1}}\mathrm{-}{\mathit{T}}_{\mathrm{0}}\mathrm{-}{\mathit{R}}_{\mathrm{01}}\mathit{/}\mathit{c}\mathrm{-}\mathit{RL}{\mathit{T}}_{\mathrm{01}}}{\mathrm{1}\mathrm{-}\mathit{ṗ}\mathrm{01}\mathit{/}\mathit{c}}\mathit{,}$(18)where $\mathit{ṗ}\mathrm{01}\mathrm{=}\frac{{{R}}_{\mathrm{01}}}{{\mathit{R}}_{\mathrm{01}}}\mathrm{·}\dot{{R}}\mathrm{0}\mathrm{\left(}{\mathit{T}}_{\mathrm{0}}\mathrm{\right)}\mathit{.}$(19)Given the solution T₀ of Eq. (13), the light-time equation for the down-leg 0 − 2 is solved for reception time T₂ in the same manner as described above.

The difference T₂ − T₁ represents the VLBI delay in the barycentric TDB-frame, which is then transformed into the geocentric TT-frame: $\begin{array}{ccc}{\mathit{t}}_{\mathrm{2}}\mathrm{-}{\mathit{t}}_{\mathrm{1}}& \mathrm{=}& \left(\frac{{\mathit{T}}_{\mathrm{2}}\mathrm{-}{\mathit{T}}_{\mathrm{1}}}{\mathrm{1}\mathrm{-}{\mathit{L}}_{\mathit{C}}}\mathrm{·}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}^{\mathrm{2}}}\left(\frac{{{V}}_{\mathrm{E}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}{\mathit{U}}_{\mathrm{E}}\right)\right]\mathrm{-}\frac{{{V}}_{\mathrm{E}}\mathrm{·}{b}}{{\mathit{c}}^{\mathrm{2}}}\right)\\ & & \mathrm{·}{\left(\mathrm{1}\mathrm{+}\frac{{{V}}_{\mathrm{E}}\mathrm{·}\dot{{r}}\mathrm{2}\mathit{,}\mathrm{gc}}{{\mathit{c}}^{\mathrm{2}}}\right)}^{-1}\mathrm{·}\end{array}$(20)The performance of the geometric delay model described here was checked by comparing it with one of the most accurate up-to-date analytical models by Sekido & Fukushima (2006), which proved to be consistent at the picosecond level.

Propagation and other additional impacts on the signal delay. Random fluctuations of a radio wave phase caused by propagation effects are dominated by the interplanetary scintillations in the solar wind. The power spectrum of these fluctuations is well-represented by the Kolmogorov spectrum (Kolmogorov 1991). The scintillation phase is given by ${\mathit{\phi }}_{\mathrm{sc}}\mathrm{=}{\left[{\mathrm{\int }}_{{\mathrm{f}}_{\mathrm{0}}}^{{\mathit{f}}_{\max }}\mathit{D}\mathrm{\left(}{\mathit{f}}_{\mathrm{0}}\mathrm{\right)}\mathrm{·}{\left(\frac{\mathit{f}}{{\mathit{f}}_{\mathrm{0}}}\right)}^{\mathrm{slope}}\mathrm{d}\mathit{f}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(21)where D(f₀) is the spectral power density, the lower limit of integration f₀ = 1/τ with τ being the length of the phase-referencing nodding cycle, the upper limit f_max is determined by the system noise level and slope is the linear slope of the spectrum in logarithmic scale. For example, with a nodding cycle of τ = 180 s, f₀ = 1/τ ≈ 5.5 mHz and the scintillation phase is proportional to the marked area in Fig. 3. Obviously, the shorter the nodding cycle, the better the target phase calibration.

The Earth’s ionosphere and troposphere also induce an additional delay to the radio signal propagation. To mitigate the effect of the former medium, we use the vertical total electron content (vTEC) maps provided by the international GNSS service (IGS) on a daily basis with a two-hour temporal resolution on a global grid (Feltens & Schaer 1998). The IGS uses the single thin ionospheric layer model, therefore the vTEC should be properly mapped onto the direction of the target source to obtain the slant total electron content (TEC): $\mathrm{TEC}\mathrm{=}\frac{\mathrm{vTEC}}{\cos {\mathit{z}}^{\mathrm{\prime }}}\mathit{,}$(22)where z′ is the zenith distance of the target as seen from the model’s ionospheric layer: ${\mathit{z}}^{\mathrm{\prime }}\mathrm{=}\arcsin \frac{\mathit{R}}{\mathit{R}\mathrm{+}\mathit{H}}\mathrm{·}\sin \mathit{z}\mathit{,}$(23)with R being the mean radius of the Earth, H the height of the ionospheric layer and z the apparent vacuum zenith distance of the target.

The ionospheric delay for each station can then be computed as ${\mathit{\tau }}_{\mathrm{iono}}\mathrm{=}\frac{\mathrm{5.308018}\mathrm{·}\mathrm{TEC}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\mathrm{·}{\mathit{f}}^{\mathrm{2}}}\mathit{,}$(24)where f is the observational frequency (e.g. 8.4 GHz for observations in X-band).

Fig. 3 Kolmogorov spectrum of phase scintillations.

There are two options for mitigating the tropospheric effect on the signal propagation delay. They are based either on “ready-to-use” empirical models or ray-tracing through the numerical weather models (NWM). The former provides the tropospheric delays in the direction of the local zenith and a mapping function to project the delay onto the direction to the source. We use the Vienna mapping functions VMF1 (Boehm et al. 2006), recommended in the IERS conventions 2010. The precision of the empirical models is sufficient for many cases, whereas we employ the latter option when an extreme accuracy is required, e.g. if the target sources are very low above the horizon and/or the angular distance between the spacecraft and the phase-referencing source is large. The latter was the case in the VEX tracking experiment EM081c, which we describe in the next section. We use the ray-tracing algorithm from Hobiger et al. (2008) with some modifications adopted from Duev et al. (2011). The meteorological data from the European centre for medium-range weather forecasts (ECMWF) are used.

Additionally, contributions to the delay due to the antenna axis offsets and smaller size effects due to thermal (Nothnagel 2009) and gravitational (Sarti et al. 2009; Clark & Thomsen 1988) deformations of telescopes are also taken into account in our delay model.

2.6. Estimation of corrections to the S/C state vector

The main task of VLBI tracking of spacecraft is estimating the state vector as a function of time. This task has some similarities with traditional astrometric applications of the VLBI technique. However, because a spacecraft is in the near-field of the VLBI-synthesised aperture, several modifications of traditional approaches and formalisms of far-field VLBI are necessary.

Analogues of uv-coverage for near-field VLBI. The uv-projections of baselines used in the standard far-field VLBI are basically the degenerated (in the limit of an infinite distance to the source) partial derivatives of the predicted geocentric delays τ_ij by the source position lateral to the line of sight: $\begin{array}{ccc}{\mathit{u}}_{\mathit{ij}}\mathrm{=}{\mathit{\partial }}_{\mathit{\varphi }}{\mathit{\tau }}_{\mathit{ij}}\mathrm{·}\mathit{c,}{\mathit{v}}_{\mathit{ij}}\mathrm{=}{\mathit{\partial }}_{\mathit{\theta }}{\mathit{\tau }}_{\mathit{ij}}\mathrm{·}\mathit{c,}& & \end{array}$(25)where c is the speed of light in a vacuum and ϕ and θ are the geocentric spherical angular coordinates of the source.

For near-field VLBI observations at the epoch t the uv-projections of the baselines i-j ($\mathit{i}\mathrm{=}\overline{\mathrm{1}\mathit{,N}\mathrm{-}\mathrm{1}}$, $\mathit{j}\mathrm{=}\overline{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,N}}$, where N is the number of participating stations), because of the finiteness of the distance to the source, should be replaced with the corresponding elements of the Jacobian J: $\mathit{J}{|}_{\mathit{t}}\mathrm{=}\left(\begin{array}{c}\\ \frac{\mathit{\partial }{\mathit{\tau }}_{\mathrm{12}}}{\mathit{\partial \varphi }}& \frac{\mathit{\partial }{\mathit{\tau }}_{\mathrm{12}}}{\mathit{\partial \theta }}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \frac{\mathit{\partial }{\mathit{\tau }}_{\mathrm{1}\mathit{N}}}{\mathit{\partial \varphi }}& \frac{\mathit{\partial }{\mathit{\tau }}_{\mathrm{1}\mathit{N}}}{\mathit{\partial \theta }}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \frac{\mathit{\partial }{\mathit{\tau }}_{\mathit{N}\mathrm{-}\mathrm{1}\mathit{,N}}}{\mathit{\partial \varphi }}& \frac{\mathit{\partial }{\mathit{\tau }}_{\mathit{N}\mathrm{-}\mathrm{1}\mathit{,N}}}{\mathit{\partial \theta }}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \frac{\mathit{\partial }\mathrm{(}{\mathit{\tau }}_{\mathrm{1}}\mathrm{-}{\mathit{\tau }}_{\mathrm{2}}\mathrm{)}}{\mathit{\partial \varphi }}& \frac{\mathit{\partial }\mathrm{(}{\mathit{\tau }}_{\mathrm{1}}\mathrm{-}{\mathit{\tau }}_{\mathrm{2}}\mathrm{)}}{\mathit{\partial \theta }}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \frac{\mathit{\partial }\mathrm{(}{\mathit{\tau }}_{\mathrm{1}}\mathrm{-}{\mathit{\tau }}_{\mathit{N}}\mathrm{)}}{\mathit{\partial \varphi }}& \frac{\mathit{\partial }\mathrm{(}{\mathit{\tau }}_{\mathrm{1}}\mathrm{-}{\mathit{\tau }}_{\mathit{N}}\mathrm{)}}{\mathit{\partial \theta }}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \frac{\mathit{\partial }\mathrm{(}{\mathit{\tau }}_{\mathit{N}\mathrm{-}\mathrm{1}}\mathrm{-}{\mathit{\tau }}_{\mathit{N}}\mathrm{)}}{\mathit{\partial \varphi }}& \frac{\mathit{\partial }\mathrm{(}{\mathit{\tau }}_{\mathit{N}\mathrm{-}\mathrm{1}}\mathrm{-}{\mathit{\tau }}_{\mathit{N}}\mathrm{)}}{\mathit{\partial \theta }}\end{array}\right)\mathit{,}$(26)where τ_i is the geocentric delay for the station i at t and the partials in the Jacobian are computed numerically. In other words, the traditional uv-formalism is merely a special case (a target source at an infinitely large distance) of a more general formalism, the Jacobians.

Spacecraft imaging and state vector estimation. The (“dirty”) image of the spacecraft at the epoch t is recovered by performing the Fourier transform of the sampled “visibility” function (see Thompson et al. 2001) ${\mathit{I}}_{\mathit{t}}\mathrm{\left(}\mathit{l,m}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{\Re }(\mathit{S}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{·}\mathit{V}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{·}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{\pi }\mathrm{i}\mathrm{·}\mathrm{(}\mathit{ul}\mathrm{+}\mathit{vm}\mathrm{)}})\mathrm{d}\mathit{u}\mathrm{d}\mathit{v}{|}_{\mathit{t}}\mathit{,}$(27)with the sampling function being the sum over all baselines of two-dimensional delta functions: $\mathit{S}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}{|}_{\mathit{t}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}\mathrm{-}\mathrm{1}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{N}}\mathit{\delta }\mathrm{(}\mathit{u}\mathrm{-}{\mathit{u}}_{\mathit{ij}}\mathit{,v}\mathrm{-}{\mathit{v}}_{\mathit{ij}}\mathrm{)}{|}_{\mathit{t}}$(28)and the normalised weighted ’visibility’ function given by $\mathit{V}\mathrm{\left(}{\mathit{u}}_{\mathit{ij}}\mathit{,}{\mathit{v}}_{\mathit{ij}}\mathrm{\right)}{|}_{\mathit{t}}\mathrm{=}{\mathit{w}}_{\mathit{ij}}\mathrm{·}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{\pi }\mathrm{i}\mathrm{·}{\mathit{\phi }}_{\mathit{ij}}}{|}_{\mathit{t}}\mathit{,}$(29)where w_ij is the weight coefficient for the baseline i-j at time t.

For a given time range t = [t_s,t_e] , the resulting image is reconstructed as a sum over all epochs within the time range: $\begin{array}{ccc}\mathit{I}\mathrm{\left(}\mathit{l,m}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{t}}{\mathit{I}}_{\mathit{t}}\mathrm{\left(}\mathit{l,m}\mathrm{\right)}\\ & \mathrm{=}& \sum _{\mathit{t}\mathrm{=}{\mathit{t}}_{\mathit{s}}}^{{\mathit{t}}_{\mathit{e}}}(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}\mathrm{-}\mathrm{1}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{N}}{\mathit{w}}_{\mathit{ij}}\mathrm{·}\mathrm{\Re }({\mathrm{e}}^{\mathrm{i}\mathrm{·}\mathrm{2}\mathit{\pi }\mathrm{·}\mathrm{(}{\mathit{u}}_{\mathit{ij}}\mathit{l}\mathrm{+}{\mathit{v}}_{\mathit{ij}}\mathit{m}\mathrm{-}{\mathit{\phi }}_{\mathit{ij}}\mathrm{)}})){|}_{\mathit{t}}\mathit{.}\end{array}$(30)The measurement equation for the epoch t is $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\phi }\end{array}{|}_{\mathit{t}}\mathrm{=}({\mathit{J}}_{\mathit{ij}}\mathrm{·}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\alpha }\end{array}){|}_{\mathit{t}}\mathit{,}$(31)where $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\phi }\end{array}$ is a [N_b × 1]  vector of differential phases on all baselines (N_b = N·(N − 1)/2 is the number of baselines) and $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\alpha }\end{array}$ is the vector of corrections to the a priori geocentric angular position of the spacecraft: $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\phi }\end{array}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{\phi }}_{\mathrm{12}}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{\phi }}_{\mathrm{1}\mathit{N}}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{\phi }}_{\mathit{N}\mathrm{-}\mathrm{1}\mathit{,N}}\end{array}\right)\mathit{,}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\alpha }\end{array}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{\Delta }\mathit{\varphi }\\ \mathrm{\Delta }\mathit{\theta }\end{array}\right)\mathit{.}$(32)Finally, the astrometric solution of Eq. (31) in the least-squares sense for each epoch t is given by $\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\alpha }\end{array}{|}_{\mathit{t}}\mathrm{=}(\mathrm{(}{\mathit{J}}^{\mathrm{T}}\mathrm{·}\mathit{J}{\mathrm{)}}^{-1}\mathrm{·}{\mathit{J}}^{\mathrm{T}}\mathrm{·}\begin{array}{}\mathrm{-}0\mathrm{\to }\\ \mathrm{\Delta }\mathit{\phi }\end{array}){|}_{\mathit{t}}\mathit{.}$(33)The main purpose for the introduction of the Jacobian formalism is to use it for the estimation of the corrections to the a priori position of a near-field VLBI target according to Eqs. (31) − (33). Imaging, although it can be used for the same purpose yielding the same result, is merely a “by-product” of such a formalism, and is used in our paper mostly for the purposes of illustration. In addition to this, the closer the target is to an observer, the more significant the difference between the Jacobians and “normal” uv-projections becomes.

3. VEX VLBI observations and results

ESA’s Venus Express Spacecraft has been observed with the European VLBI Network (EVN) within the framework of the PRIDE initiative under the EVN observing project EM081.

From four requested observing epochs we used only three. The first two observing runs, on 2010 August 23 (2 h long, 3 stations) and on 2010 September 20 (3 h long, 2 stations) were mainly devoted to tuning-up and verification of the technique, namely developing and testing the S/C observation-related scheduling, data acquisition, transfer and processing methods with the near-field delay model.

Here we report the results obtained with the third observational run, which has been the most intensive experiment so far in terms of the number of participating telescopes: the EVN stations Onsala (On, Sweden), Wettzell (Wz, Germany), Medicina (Mc, Italy), Matera (Ma, Italy), Yebes (Ys, Spain), Metsähovi (Mh, Finland), Svetloe (Sv, Russia), Zelenchukskaya (Zc, Russia), Hartebeesthoek (Hh, South Africa), and the NRAO Very Long Baseline Array (VLBA) station St. Croix (Sc, USA).

The observations were conducted on 2011 March 28, from 08:45 to 11:30 UTC. The source J2225-0457 (also known as 3C 446) was used as a fringe finder and for the clock search. The source J2211-1328 with J2000.0 coordinates of ${\mathrm{22}}^{\mathrm{h}}{\mathrm{11}}^{\mathrm{m}}\mathrm{24}\stackrel{\mathrm{s}}{.}\mathrm{0994590}\mathit{,}\mathrm{-}{\mathrm{13}}^{\mathrm{◦}}\mathrm{28}\prime \mathrm{09.723950}″$ (Petrov 2011) was used as a phase-reference calibrator. The estimated coordinate error for this source is 0.11 mas in Right Ascension and 0.12 mas in Declination (Petrov 2011). All involved EVN radio telescopes were configured to observe the spacecraft signal in a phase-referencing mode, switching between the S/C and the phase reference source with a nodding cycle of 4 min, including a 20 second gap for re-pointing and antenna calibration. At On, Mc, Ma, Mh, Hh, we recorded 79 scans per station, including 2 for the fringe finder, 39 for the calibrator source and 38 for VEX. The stations Ys and Sc joined in starting from the 15th scan owing to technical problems and limited visibility, respectively. The eastern-most stations Sv and Zc stopped observing according to schedule owing to the target sources’ setting. During the observing run, the VEX spacecraft had mean coordinates 22^h14^m13^s, − 11°41′22″. At the time of the observations, Venus was at a solar elongation of 36 degrees and at a distance from Earth of 1.24 AU. Figure 4 shows the orbital position of VEX around Venus during the observations as seen from the geocentre. Figure 5 shows the pseudo-uv-coverage (Jacobians) of the EM081c VLBI array.

Fig. 4 VEX orbit around Venus on 2011 March 28, as seen from the geocentre. Red boxes represent the position of VEX at the epochs when imaging was performed.

Fig. 5 Near-field Jacobians for the VEX pointings in kilometers, 2011 March 28, the EVN project EM081. The longest east-west baseline – Ys-Zc, the longest north-south baseline – Mh-Ys.

It should be pointed out that the observational setup for the experiment was not ideal – target declinations ranging from  −11° to  −13° are unfavourable for the EVN array. In addition to this, the angular separation of  ~2.5° between Venus and the calibrator was at the edge of the isoplanatic patch⁵ at low elevations caused by a low declination of the target, which lead to a relatively large unmodelled tropospheric/ionospheric phase difference between the calibrator and VEX.

Fig. 6 Frequency Doppler shift over the path of the ESA tracking station Cebreros – VEX – Metsähovi, Medicina, Yebes (top left), Cebreros – VEX – Metsähovi, Medicina, Yebes – geocentre (top right). The differential frequencies in the geocentre with a standard deviation of  ~5 mHz at 2.5 s sampling are shown at the bottom. 2010 August 23, EVN project EM081.

Fig. 7 Interferometric fringes of the VEX signal in the spectral domain, 1024 frequency points, baseline Onsala-Zelenchukskaya, 2011 March 28, EVN project EM081, using the near-field delay model and two-way Doppler phase correction to avoid frequency smearing when correlating at the SFXC at JIVE. Top: amplitude in relative units; bottom: phase in radians. The horizontal axis – frequency in video band (MHz) relative to 8.411990 GHz. Most of the power of the spacecraft signal is concentrated within a  ~1 MHz-wide band around the carrier line, although the S/N for the carrier is considerably higher than for the rest of the band.

Fig. 8 a), b) Fringe finder (3C 446, 2.5 Jy source at an angular distance of 15° north from the target) and a calibrator (J2211-1328, 0.17 Jy source at an angular distance of 2.5° south from the target) self-calibrated images obtained using the EVN array excluding Hh, c), d) self-calibrated images of the same sources obtained with the whole array. Horizontal axis – deviation from the nominal source’s Right Ascension in mas, vertical axis – deviation from the nominal source’s Declination in mas. Dashed white ellipses at the bottom-right represent the synthesised beam size and shape.

Fig. 9 Radio images of VEX obtained with the EVN array. The horizontal axis represents the deviation from the nominal source’s Right Ascension in mas, vertical axis – deviation from the nominal source’s Declination in mas. a) 09h05m TDB, stations On, Wz, Ma, Mc, Mh, Sv, Zc; b) 09h30m TDB, Ys joined in; c) 09h55m TDB, Sv left; d) 10h20m TDB, Mh left; e) 10h45m TDB, Zc left. No cleaning applied.

The data were acquired using the standard VLBI data acquisition system based on the Mark5 A/B recording equipment. We used a base frequency of 8411.99 MHz, four 16-MHz-wide channels with 2-bit Nyquist sampling, resulting on aggregate data rate of 256 Mbps per station. The data files were either transferred over the high-capacity network connection after the experiment to Metsähovi and JIVE or recorded on disks and then shipped to the processing centre at JIVE.

The narrow-band data processing and carrier frequency/phase extraction were conducted with the SWSpec and SCtracker software at Metsähovi. The Doppler frequency detections from all stations reduced to the geocentre using the near-field delay model (as ) differ on average only by several mHz or even less, proving the consistency between the actual and model delays. As an example, the topocentric frequency detections measured for the stations Mh, Mc and Ys, detections reduced to the geocentre and differential geocentric frequencies for the experiment conducted on 2010 August 23 are shown in Fig. 6.

The broadband correlation was conducted at JIVE using the SFXC. The fringe finder data were used to perform a clock search. The VEX cross-correlation spectrum is smeared due to the intrinsic change of the frequency emitted by the S/C, as it retransmits the signal in the two-way Doppler regime. This change may reach several kHz over a time span of one hour, which is shown in Fig. 6. To avoid this smearing when correlated with the SFXC, a two-way Doppler correction (see formula (1) and the paragraph below it) was applied to the S/C data. Figure 7 represents the amplitude and phase of the CCF of the VEX signal (amplitude and phase) on the baseline On-Mh. Most of the spectral power is concentrated within a  ~1 MHz region around the carrier line, allowing us to resolve the 2π ambiguity between the group delay detected by the broad-band correlation with the SFXC and phase delay detected using the SCtracker software. Figure 8a, b shows the self-calibrated and cleaned images of the fringe finder and the calibrator obtained with the Difmap (Shepherd 1995) and AIPS (Greisen 2003) software using the EVN array excluding Hh.

Unfortunately, an anomalous clock rate offset at Hh did not allow us to obtain its absolute phases. Although good for self-calibrated imaging (see Fig. 8), the Hh data were not used for absolute positioning during the EM081c observational run on 2011 March 28. The Hh LO converter was fixed several days later, based on our analysis of the LO behaviour made with ultra-high spectral resolution SCTrack software package.

Fig. 10 (Top) Measured coordinate offsets with respect to the a priori trajectory, supplied by ESOC. The X-offset is shown in blue, Y-offset  −  in red. X and Y are approximately RA and Dec. The measurement sampling  −30 s, 1 km corresponds to  ≈1.1 mas at the distance of 1.24 AU. (Bottom) Error scatter plots; the light blue ellipses show the a priori estimates of the trajectory accuracy, red ellipses – our measurement’s 3σ scatter for a full time-range (left) and best-calibrated time range (right).

The correlation of the calibrator source for each baseline provides the residual delays. For coherent phase-referencing, the residual delays must be lower than the radio wave period at X-band, which is about 0.120 ns. The measured residual delays and phases were applied to the detected residual phases of the target. With these we performed imaging of the VEX spacecraft and derived the deviations of the lateral coordinates of the spacecraft with respect to a priori ones. Figure 9 shows the imaging results obtained for VEX with the EVN array without any cleaning. The averaging time for each image is 20 min. Figure 10 shows the estimated deviations of the a priori lateral position of the VEX spacecraft. The integration time for each point is 30 s. The coordinate offsets were estimated using Eq. (33). However, for a consistency check, we made S/C images for each short interval of 30 s with a “super”-resolution and measured the peak position on each of them. The results yielded a consistency of the two approaches at the level of 25 microarcsec. The orbit accuracy estimate is at a 3 sigma level of 200−300 m across the track and 500 − 600 m along the track.

The north/south VEX offsets of  ~1 mas from the expected orbital position seen in Fig. 10, although nearly within the error budget, may be caused by systematic phase errors between the target and calibrator due to their low elevation. More experiments with different relative positions of the S/C and calibrator(s) are needed to prove or disprove this claim.

4. Discussion and conclusions

The EVN experiment EM081 enabled us to show that astrometric spacecraft positioning with a very high accuracy is achievable using the PRIDE approach. A relatively large angular distance between VEX and the calibrator source (2.5°), their low elevation and a limited accuracy of the troposphere/ionosphere models did not allow us to include the data from St. Croix (Sc) in the astrometric solution without a baseline fringe phase 2π-bifurcation. Figure 8c, d shows the self-calibrated images of the fringe finder and the calibrator with Sc and Hh added. As seen from this figure, more than a three-fold resolution improvement (at the level of up to 20 mas, 0.1 nrad) is possible. We adopt this value as the aim for our future experiments. To achieve this increase in resolution, we will have to resolve the 2π ambiguities for the longer baselines. For this purpose, closer (in-beam in the ideal case), preferably multiple calibrators should be used in the experimental setup for phase referencing. In addition to this, better troposphere and ionosphere models would be very important, especially for the low elevations of the sources.

The results of the current study in terms of achieved accuracy are compatible with those of the other recent well-established VLBI spacecraft tracking projects, such as the VLBA tracking of the Cassini spacecraft (Jones et al. 2011) and NASA’s Mars Exploration Rover B spacecraft during its final cruise phase (Lanyi et al. 2007).

The PRIDE approach described in this paper proves its applicability to virtually any deep-space mission almost anywhere in the solar system. We note that the approach is not exclusively oriented on planetary science missions, but other deep-space and near-Earth targets. This has been demonstrated by VLBI observations of GLONASS satellites (Tornatore et al. 2010, 2011) aiming to establish a better link between the terrestrial and celectial reference frames.

Acknowledgments

We would like to express our sincere gratitude to T. Morley (ESA ESOC) for persistent support, M. Pérez-Ayúcar (ESA ESAC) for robust assistance and the anonymous referee for useful and constructive suggestions. 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 authors would like to thank the personnel of the participating stations: P. de Vicente (Yebes), J. Quick (Hartebeesthoek), G. Kronschnabl (Wettzell), R. Haas (Onsala), A. Orlatti (Medicina), G. Colucci (Matera), A. Finkelstein, M. Kharinov, A. Mikhailov (Svetloe and Zelenchukskaya of the KVAZAR network) and the staff of the NRAO VLBA. R. Campbell and A. Keimpema of JIVE provided important support to various components of the project. We also thank the Royal Netherlands Meteorological Institute (KNMI) for granting access to the ECMWF meteorological data. D. A. Duev acknowledges the EC FP7 project EuroPlaNet (grant agreement 228319). G. Cimó acknowledges the EC FP7 project ESPaCE (grant agreement 263466). T. Bocanegra Bahamon acknowledges the NWO-ShAO agreement on collaboration in VLBI.

References

All Figures

	Fig. 1 Data processing pipeline of phase-referencing VLBI observations of spacecraft.
In the text

	Fig. 2 Geometry of VLBI observations of spacecraft in the barycentric celestial reference frame.
In the text

	Fig. 3 Kolmogorov spectrum of phase scintillations.
In the text

	Fig. 4 VEX orbit around Venus on 2011 March 28, as seen from the geocentre. Red boxes represent the position of VEX at the epochs when imaging was performed.
In the text

	Fig. 5 Near-field Jacobians for the VEX pointings in kilometers, 2011 March 28, the EVN project EM081. The longest east-west baseline – Ys-Zc, the longest north-south baseline – Mh-Ys.
In the text

	Fig. 6 Frequency Doppler shift over the path of the ESA tracking station Cebreros – VEX – Metsähovi, Medicina, Yebes (top left), Cebreros – VEX – Metsähovi, Medicina, Yebes – geocentre (top right). The differential frequencies in the geocentre with a standard deviation of  ~5 mHz at 2.5 s sampling are shown at the bottom. 2010 August 23, EVN project EM081.
In the text

Fig. 7 Interferometric fringes of the VEX signal in the spectral domain, 1024 frequency points, baseline Onsala-Zelenchukskaya, 2011 March 28, EVN project EM081, using the near-field delay model and two-way Doppler phase correction to avoid frequency smearing when correlating at the SFXC at JIVE. Top: amplitude in relative units; bottom: phase in radians. The horizontal axis – frequency in video band (MHz) relative to 8.411990 GHz. Most of the power of the spacecraft signal is concentrated within a  ~1 MHz-wide band around the carrier line, although the S/N for the carrier is considerably higher than for the rest of the band.

In the text

Fig. 8 a), b) Fringe finder (3C 446, 2.5 Jy source at an angular distance of 15° north from the target) and a calibrator (J2211-1328, 0.17 Jy source at an angular distance of 2.5° south from the target) self-calibrated images obtained using the EVN array excluding Hh, c), d) self-calibrated images of the same sources obtained with the whole array. Horizontal axis – deviation from the nominal source’s Right Ascension in mas, vertical axis – deviation from the nominal source’s Declination in mas. Dashed white ellipses at the bottom-right represent the synthesised beam size and shape.

In the text

	Fig. 9 Radio images of VEX obtained with the EVN array. The horizontal axis represents the deviation from the nominal source’s Right Ascension in mas, vertical axis – deviation from the nominal source’s Declination in mas. a) 09h05m TDB, stations On, Wz, Ma, Mc, Mh, Sv, Zc; b) 09h30m TDB, Ys joined in; c) 09h55m TDB, Sv left; d) 10h20m TDB, Mh left; e) 10h45m TDB, Zc left. No cleaning applied.
In the text

Fig. 10 (Top) Measured coordinate offsets with respect to the a priori trajectory, supplied by ESOC. The X-offset is shown in blue, Y-offset  −  in red. X and Y are approximately RA and Dec. The measurement sampling  −30 s, 1 km corresponds to  ≈1.1 mas at the distance of 1.24 AU. (Bottom) Error scatter plots; the light blue ellipses show the a priori estimates of the trajectory accuracy, red ellipses – our measurement’s 3σ scatter for a full time-range (left) and best-calibrated time range (right).

In the text