Laser and Radio Tracking for Planetary Science Missions - A Comparison

At present, tracking data for planetary missions largely consists of radio observables: range-rate range and angular position. Future planetary missions may use Interplanetary Laser Ranging (ILR) as a tracking observable. Two-way ILR will provide range data that are about 2 orders of magnitude more accurate than radio-based range data. ILR does not produce Doppler data, however. In this article, we compare the relative strength of radio Doppler and laser range data for the retrieval of parameters of interest in planetary missions, to clarify and quantify the science case of ILR, with a focus on geodetic observables. We first provide an overview of the near-term attainable quality of ILR, in terms of both the realization of the observable and the models used to process the measurements. Subsequently, we analyze the sensitivity of radio-Doppler and laser-range measurements in representative mission scenarios. We use both an analytical approximation and numerical analyses of the relative sensitivity of ILR and radio Doppler observables for more general cases. We show that mm-precise range normal points are feasible for ILR, but mm-level accuracy and stability is unlikely to be attained, due to a combination of instrumental and model errors. ILR has the potential for superior performance in observing signatures in the data with a characteristic period of greater than 0.33-1.65 hours This indicates that Doppler tracking will typically remain the method of choice for gravity field determination and spacecraft orbit determination in planetary missions. Laser ranging data, however, are shown to have a significant advantage for the retrieval of rotational and tidal characteristics from landers. Similarly, laser ranging data will be superior for the construction of planetary ephemerides and the improvement of solar system tests of gravitation.

We show that mm-precise range normal points are feasible for ILR, but mm-level accuracy and stability in the full analysis chain is unlikely to be attained, due to a combination of instrumental and model errors. We find that ILR has the potential for superior performance in observing signatures in the data with a characteristic period of greater than 0.33-1.65 hours (assuming 2-10 mm uncertainty for range and 10 µm/s at 60 s for Doppler). This indicates that Doppler tracking will typically remain the method of choice for gravity field determination and spacecraft orbit determination in planetary missions. ILR data will be able to supplement the orbiter tracking data used for the estimation of parameters with a once-per-orbit signal. Laser ranging data, however, are shown to have a significant advantage for the retrieval of rotational and tidal characteristics from landers. Similarly, laser ranging data will be superior for the construction of planetary ephemerides and the improvement of solar system tests of gravitation, both for orbiter and lander missions.
Keywords Interplanetary Laser Ranging · Radio Tracking · Planetary Missions

Introduction
For both Earth-orbiting and planetary missions, a variety of tracking observables is available from which the trajectory of the spacecraft can be reconstructed. A tracking data type that will become available for future planetary missions is Interplanetary Laser Ranging (ILR) (Degnan, 2002). The technology for such a system derives strongly from Satellite Laser Ranging (SLR; Pearlman et al., 2002), Lunar Laser Ranging (LLR; Murphy, 2013) as well as Laser Time Transfer (LTT; Exertier et al., 2014).
The key difference between SLR/LLR and ILR is that ILR requires a transponder on both the ground and space segments, as opposed to the purely passive space segment in SLR/LLR (retroreflectors).
Currently, deep-space missions largely rely on the use of radio tracking for their orbit determination and the associated parameter estimation. In particular Doppler data has been the primary data type for this application, (e.g. Konopliv et al., 2011;Iess et al., 2012;Mazarico et al., 2014). Radio-based range data supplement the Doppler data by providing the absolute distance between spacecraft and ground station. The angular position of the spacecraft in the sky can be measured by means of Very Long Baseline Interferometry (VLBI) (Duev et al., 2012). In the so-called PRIDE (Planetary Radio Interferometry and Doppler Experiments) setup, VLBI and Doppler data are obtained concurrently (Duev et al., 2016). Unlike the Doppler data, range and VLBI data are used primarily for the estimation of solar system ephemerides (Fienga et al., 2009;Jones et al., 2015;Dirkx et al., 2017), which provide crucial input for experimental relativity (e.g. Will, 2014). The selection of ILR as a tracking type in future missions is contingent upon its data being able to provide scientific results that are complementary/supplementary to, or competitive with, the results obtained from existing systems, such as those mentioned above.
There has been a number of experimental demonstrations of ILR, both in one-way (Abshire et al., 2006;Neumann et al., 2008;Noda et al., 2017) and two-way  modes. The only operational implementation of ILR to date has been on the Lunar Reconnaissance Orbiter (LRO), using the laser altimeter system (Zuber et al., 2010;Bauer et al., 2016;Mao et al., 2017). None of the demonstrations of ILR have been performed with dedicated hardware, and the attainable measurement accuracy has not yet been pushed to the limit of state-of-the-art. In recent years, there has been substantial development in LTT (Exertier et al., 2014;Samain et al., 2015;Exertier et al., 2017;Prochazka et al., 2017a), which is in many ways similar to a transponder ranging system. Laser communication technology is also maturing for use in planetary missions, as demonstrated by the LLCD demonstrator on the LADEE lunar orbiter. As an ad hoc product, LLCD communications data were used to obtain two-way ranging data with a precision of several cm (Stevens et al., 2016).
Our goal is to provide a comparison of radio-Doppler and laser-range data for planetary missions, and to identify the areas where the addition of ILR data would be beneficial for attaining the scientific objectives of planetary missions. A preliminary set of results is given by Dirkx (2015). Instead of focusing on a single mission concept in detail, we take a broad view and quantitatively analyze the relative strength of the range and Doppler observables for parameters of interest in a variety of planetary missions. As a result, we identify classes of science products/mission scenarios for which ILR will be a competitive design option.
We give an overview of tracking data in Section 2 and discuss the expected error budget of ILR in Section 3. In Section 4, we present the methods we use to compare ILR and Doppler data. We use both an analytical approach, and a numerical covariance analysis based on simulated orbit determination/parameter estimation. The numerical technique serves to identify the level of applicability of the analytical approach, and to provide guidelines on how to apply it in mission design and analysis. In Section 5 we show the results of our comparison of the two data types. In Section 6, we discuss the implications of these results for the use of ILR data , with a focus on geodetic observables. We conclude with a discussion on the overall science case of ILR in Section 7.

Planetary Tracking Data
In this section, we give a general overview of planetary tracking data, including the models for both range and Doppler observables.

Radio Data
The typical precisions of radio tracking data are currently at the level of 0.02-0.1 mm/s at 60 s integration time for range-rate measurements at X-band (e.g., Thornton and Border, 2000;Marty et al., 2009;Konopliv et al., 2011;Iess et al., 2014;Bocanegra-Bahamón et al., 2018) and 0.5-5 m for range measurements (e.g., Thornton and Border, 2000;Folkner et al., 2014). Detailed discussions on sources of both systematic errors and random noise are given by Thornton and Border (2000); Moyer (2000); Asmar et al. (2005); Iess et al. (2014);Molera Calvés et al. (2014). For Doppler measurements, the systematic errors are typically close to negligible (e.g. Iess et al., 2014). In contrast, the level of systematic errors of radio range measurements can be quite large, comparable to the random noise, at the 1 m level.
The radio tracking noise is dominated by the propagation effects in the interplanetary medium, and depend strongly on the solar elongation angle. For Mars Express, Duev et al. (2016) show that PRIDE Doppler noise at X-band is indeed dominated by the these effects, with a median one-way value of ∼ 30µ/s at 10 s integration time. Combining observations at multiple wavelengths allows the removal of the majority of these errors (Reasenberg et al., 1979;Bertotti et al., 2003). The combined X-and Ka-band approach was/is used for Cassini (Kliore et al., 2004) and Juno (Iess et al., 2018), and is anticipated for use on the BepiColombo and JUICE missions, which have tracking data quality requirements of 0.01 mm/s at 60 s integration time and 3 µm/s at 1000 s integration time. Additionally, these missions will use an advanced radio ranging system, allowing two-way range measurements with a predicted accuracy of down to 20 cm.

Laser Ranging -Measurement Concepts
SLR has been used for Earth-orbiting satellites for more than 50 years (Pearlman et al., 2002), and provides twoway range data with sub-cm accuracy. In ILR, the use of retroreflectors is no longer feasible due to the large target distance requiring an active system in both the ground and space segments.
Two main types of active laser ranging systems (sometimes termed transponder laser ranging) can be distinguished for use in planetary missions (Degnan, 2002;Birnbaum et al., 2010, see Fig. 1): -One-way laser ranging. A laser pulse is transmitted from a ground station and detected by a (satellitebased) receiver (or the other way around). An important issue with this method is that the transmission and reception times are measured by different, unsynchronized, clocks (Bauer, 2017). -Two-way asynchronous laser ranging. Both the space and ground segment fire laser pulses towards one another independently. By pairing a range measurement from the up-and down-link, a two-way range measurement is obtained (Birnbaum et al., 2010), Laser transmitter Receiver Laser transmitter Receiver Laser transmitter Receiver which does not suffer from the clock error issue of a one-way system. Two-way asynchronous laser ranging has been the method of choice in most ILR mission proposals, due to its higher data quality and less stringent clock requirements. In this article our main focus is on two-way systems.
Both the one-and two-way observable are created from time tags of the transmission (τ t ) and reception (τ r ), in their local proper time scales 1 . The error-free two-way raw range measurement is then created from the coordinate times t as: Due to uncalibrated errors in the measurements (Section 3) of τ t and τ r , the range measurement quality is degraded. We denote the measured time (transformed to coordinate time) ast, so that the measured ranges becomes: 1 In the context of this article, we use the symbol t to denote either a coordinate time such as Barycentric Coordinate Time (TCB) or a scaled coordinate time such as Dynamical Barycentric Time (TDB). Although TDB is decidedly not a coordinate time, as discussed in detail by Klioner (2008), the distinction is not relevant for the purpose of our discussion. Details of the conversion between coordinate and proper time scales are discussed by Soffel et al. (2003), while the influence on the data analysis is discussed by  where ε s is the lumped range error. Note that this range value includes the effects of atmospheric and relativistic effects.

Observation Models
We will denote the one-way range observation between point A (transmitter) and point B (receiver) by s (1) BA . From their position functions, denoted r A (t) and r B (t), respectively, the one-way range is obtained from: where the formulation is referenced to the reception time t r , here equal to a given t 2 . The term ∆s (1) BA denotes range corrections due to e.g., propagation medium and relativistic effects.
For a one-way range-rate (or Doppler) observable, denoted here asṡ (1) BA , with an integration time denoted by ∆t i , the observable is modelled as (Moyer, 2000): where ∆t i = t 4 − t 2 . Here, two one-way range observables with reception times t 4 and t 2 (and associated transmission times t 3 and t 1 ) are used. Typical values of ∆t are 1-60 s, but may be > 1000 s in certain cases. The two-way range is modelled as the combination of two one-way ranges: where δt B represents the delay between the reception and retransmission of the signal at station B, typically < 1 ms for radio data (Bertone et al., 2018), up to the order of 1 minute for ILR .

Laser ranging data -Error Sources
The error budgets of radio tracking systems are well understood and quantified (Section 2.1). Here we analyze and discuss the various sources of error in ILR measurements/analysis.

Inherent Measurement Uncertainty
The primary measurement uncertainty of ILR is a convolution of two main contributors: the laser pulse profile and the detector impulse response. For a pulse with a perfectly Gaussian temporal profile and single-photon intensity detection energies, the pulse profile introduces a purely Gaussian measurement error (Murphy, 2001;Dirkx et al., 2014a). The one-way range precision limit due to only this effect is 1.3 to 13 mm single-shot root mean square (RMS) for laser pulses with a full-width half maximum (FWHM) of 10 to 100 ps, respectively. Note that the actual single-shot precision may be limited by the detector random error (Section 3.2). In cases where the number of detectable photons is larger than one, additional range biases at the several mm level may be incurred (Dirkx et al., 2014a), if traditional threshold detection is used. Degnan (2017) has proposed the use of centroid detection that could be used to overcome this bias, by combining all incoming detections from a single pulse into a single waveform. In ILR, however, detection energies are expected to be at the single-photon level.
The contribution due to the atmospheric turbulence (Gardner, 1976) is in most cases below the 0.5 mm level (Kral et al., 2005) (see Section 3.3 for discussion of troposphere model errors).
In SLR/LLR, the retroreflector signature causes significant distortion of the temporal pulse shape (Otsubo and Appleby, 2003), which is not the case for ILR. Therefore, this aspect of the data stability could be better for ILR than for SLR, as the incoming pulse energy temporal profile is more predictable in ILR.

Measurement Errors
For ILR, the hardware-derived error sources are similar to those of SLR, which were summarized by Exertier et al. (2006). Although much of the error budget remains close to that given there, we note several recent changes in the coming sections. An important difference between SLR and ILR stems from the fact that in ILR part of the active hardware is on the space segment, and no passive reflectors are used.
For the characterization of the space segment, we rely in part on the development of space-grade detection systems that are currently in operation, such as T2L2 (Exertier et al., 2014), and those that are under development, such as the single-photon ELT . Although the optical components of T2L2 and ELT are very different from those in an ILR system, the stable single-photon detection system has very similar characteristics.
The measurement error due to the detector, often a Silicon Photon-Avalanche Diode (SPAD) is, at 3-6 mm, the largest contributor to hardware-induced error budget of SLR (Exertier et al., 2006). Space-grade detectors showing sub-ps stability have been developed for the ELT project Kodet and Prochazka, 2012). The contribution to the ILR precision of the best photon-counting detectors is typically 3 mm RMS (Prochazka et al., 2017b). Exertier et al. (2006) give values of several mm for the influence of jitter in the event timer. A novel type of event timer, developed and applied by Panek et al. (2010); , provides sub-ps precision and a stability of several fs over a period of minutes to hours. The use of this technology allows the event timer to have an almost negligible contribution to the range error budget.
For ILR, the influence of clock noise is substantially different from SLR/LLR (Degnan, 2002;Dirkx et al., 2015). For two-way systems, the clocks only need to be sufficiently stable over short periods of time (two-way light time). A stability of about 10 −15 over a typical ILR light time of 1000 s will result in 1 ps timing error (0.3 mm one-way range error), and is achievable by H-masers (e.g. Dehant et al., 2017). For the space segment, stability is only required over the time δt, putting sub-mm errors well within the capabilities of present spaceborne systems. A one-way range system requires clocks at both ends of the link to be stable over longer time periods (Bauer et al., 2016), making clock noise a significant issue.
Delays in various components of the electrical and optical system of the detection assembly must be accurately characterized to realize a high-quality range measurement. Kirchner and Koidl (2014) show that ground station calibration consistency on short time scales is at the several ps level averaged over 10 s, comparable to the value of 3 ps given by Prochazka et al. (2012) for ELT. Both are in line with the value of 1 mm given by Exertier et al. (2006). Nevertheless, consistently obtaining mm-level system calibration on the space segment will be challenging.
Data from existing two-way ILR experiments cannot be used to set up a measurement error budget. The only two-way ILR experiment thus far used the non-dedicated hardware on the MESSENGER spacecraft , which is not representative of the state-of-the-art. Two-way links have been demonstrated on laboratory scales. Chen et al. (2013) obtain range measurement errors below the 0.2 mm level (averaged over 1000 measurements). Blazej et al. (2014) have shown time transfer with an accuracy of 3 ps (≈ 1 mm) using two representative ground segment hard-ware packages. These experiments show the capabilities of laboratory-scale experiments with well-controlled hardware, indicating the potential for (sub-) mm range accuracy.

Data Analysis Uncertainties
Even in the case of perfect range measurements (ε s = 0), the science products obtained from the data will not be error-free. Errors in the evaluation of Eq. (4) will degrade the fidelity of the results.
The position of the ground station in the Geocentric Celestial Reference System (GCRS) requires a timedependent position in the International Terrestrial Reference Frame (ITRF), as well as a rotation between the two. Inaccuracies in these models limit the accuracy of the ground station position function at the sub-cm level (Altamimi et al., 2011;Rothacher et al., 2011;Sośnica et al., 2013).
ILR analysis must be done in the Barycentric Celestial Reference System (BCRS). As a result, uncertainties in the a priori Earth ephemeris will enter the error budget of the ground station position model. The ephemeris of the Earth is orders of magnitude less accurate than the expected cm-accuracy of ILR measurements (Fienga et al., 2011). This indicates that the Earth's ephemeris should be among the estimated parameters during ILR data analysis.
For landers on solar system bodies, the general issues in modelling the barycentric state are of a similar nature as for Earth ground stations. Depending on the target body, however, the uncertainty may be limited and accounted for by the addition of a number of estimated parameters. In fact, the signatures of these effects will often be key science objectives of the lander tracking (e.g. Le Maistre et al., 2013;Dirkx et al., 2014b). Fulfilling the modelling requirements may require significant theoretical work. For planetary dynamics, the uncertainty in asteroid masses and orbits is presently the limiting factor in the dynamical models (Fienga et al., 2009). For the dynamics of natural satellite systems, the consistent coupling between translational and rotational dynamics and tidal deformation will be challenging to model at the mm-level (e.g. . For orbiters, uncertainties in non-conservative forces, as well as target body gravity field variations (e.g. Marty et al., 2009), can limit the accuracy to which the dynamics can be modelled. This requires the state estimation to be performed arc-wise (a typical arc length is several days). Dynamical model error can become the dominant source of uncertainty in the estimated parameters, and is a key reason for the true estimation errors often being significantly higher than the formal estimation errors (e.g. Konopliv et al., 2011;Mazarico et al., 2014).
Finally, models for tropospheric correction (as part of ∆s BA ) have an accuracy of 5-8 mm (Exertier et al., 2006). This level can be reduced to the (several) mm level using ray-tracing models (Hulley and Pavlis, 2007). Detailed models for relativistic range corrections have been developed (e.g. Teyssandier and Le Poncin-Lafitte, 2008), so model uncertainties for this contribution of ∆s BA will be negligible if state-of-the-art models are applied.

Total Uncertainty -Summary
In Sections 3.1-3.3, we have presented the main error sources that enter the data realization and analysis chain of ILR. The main sources of measurement error are the detector uncertainty (at 3 mm), and the finite pulse length (3 to 13 mm RMS, for laser pulses 10 to 100 ps FWHM respectively).
As a result, 1.0 to 4.3 mm precision averaged over 10 measurements (for 10 to 100 ps pulse length) may be achieved. Considering the existing detector and timing devices performance (Section 3.2), a limiting precision (but not accuracy) <0.1 mm can be achieved when averaging higher number of individual single photon measurements.
Hardware imperfections, as they are deduced from current SLR and space-based laser transmitter and detector systems, will induce systematic errors at the several mm level, as is the case for SLR.
It is especially the instabilities in the systems that will be an issue for the quality of ILR data (as well as for SLR/LLR). A system bias that is constant for an extended time interval can be mitigated by the estimation of a single long-arc range bias. Instabilities (e.g. random walk behaviour) cannot be removed in this manner without introducing an excessive number of parameters. Therefore, ILR data accuracy will be limited to the level of several mm at best. Considering the current performance of SLR systems, sub-cm accuracy will be feasible. A major design driver for the space segment will be the stable system delay calibration in an (inter)planetary environment. Existing model errors in tropospheric correction and ground station positioning will limit ILR data modelling to the several mm level (Section 3.3). These issues are also crucial in space geodesy, and are a topic of active research.
Dynamical model errors of both the space segment and the Earth will limit the accuracy to which the data can be interpreted. The degree to which these errors will affect the estimation of parameters of interest is strongly dependent on the correlation of the signatures of these parameters with the model errors. This error source is similar for both Doppler and range data, and can prevent a data set from being exploited to its full potential.

Data-Type Comparison Methodology
In this section, we present the methods we use to compare ILR with radio data. We outline our concept for analytical comparison and numerical covariance analysis in Sections 4.1 (based on Dirkx (2015)), and Section 4.2, respectively.

Analytical Approach
The sensitivity of an observable h to a parameter q is quantified by its associated partial derivative ∂h/∂q, (e.g., Montenbruck and Gill, 2000). Since the range observable s and range-rate observableṡ are related through Eq. (5), we have: To quantitatively compare the data types, we define a signal-to-noise (SNR) criterion for an observable h and a parameter q, denoted SNR h;q , which is computed as follows: where σ h is the noise level of the measurement h. Now, we define the following figure of merit to compare the relative sensitivity of range and Doppler observables to a parameter q: where the maximum is taken over the observational period. To first order, we can set Ξ q < 1 as a criterion when ILR would become a feasible alternative to Doppler data for determining a parameter q. The interpretation of numerical values of this criterion is discussed in Section 6.1. We start by using an analytical model for Eq. (8), in which the influence of a parameter q is manifested in the range measurements as a purely sinusoidal signal of amplitude A and angular frequency ω (period denoted as T ), so that: and we obtain the following from Eq. (8): This approximates the situation where all planets are in circular orbits around a static Sun, with the spacecraft in a circular orbit around one of these bodies, and the parameter q imparting an N -cycles-perorbit sinusoidal signal on the data (with N an integer). In the data, the orbital frequencies of the spacecraft, Earth and the target planet will then all be visible 2 . These assumptions are a reasonable approximation for our analysis, as validated in Section 5.2 from numerical results. In Appendix A, we discuss how to extend the method to elliptical orbits.
For ∆t T , we obtain the following from Eq. (12): so that: We apply this limit approximation in Sections 5 and 6 to compare the analytical and numerical approaches.

Covariance Analysis -Numerical Models
To assess the validity of our analytical results, and gain insight into how the results obtained from them should be interpreted, we perform covariance analyses (e.g. Montenbruck and Gill, 2000;Milani and Gronchi, 2010) for a number of representative cases. We generate formal errors for a set of parameters when using only Doppler data, and when using only range data. We denote the formal error of parameter p, using data type h, as p,h . 3 We simulate two scenarios: a Mars lander mission, and a Mars dual-orbiter mission. The orbits of the spacecraft are both low-altitude, nearly polar and nearly circular, with initial condition e = 0.01 for both; a = 3850 km and i = 88 • for one orbiter and i = 92 • and a = 3800 km for the other orbiter (similar to spacecraft such as Mars Odyssey and Mars Reconaissance Orbiter). The lander is placed equatorially. For both 2 For various parameters, the sinusoidal signature will be modulated by a linear trend, increasing the amplitude with time. However, if the observation time is much larger than the period of the signal T , the impact of this linear trend on Eq. (10) will be small.
3 Note that we use the symbols ε (lumped range error) and to represent different physical quantities.
cases, we assume a simplified Mars rotation model, with fixed pole right ascension and declination α M and δ M , and a fixed rotation rate ω M . For the orbiter simulation, we estimate the initial states x i of both orbiters (i = 1, 2) w.r.t. Mars' center of mass over 1-week arcs, with each arc j starting at time t j over a 2-year period (j = 1..104). For both the lander and the orbiter simulations, we estimate Mars' dynamics over a single 2-year arc w.r.t. the barycenter x M (t 0 ) (see Table 1)..
We set up our state transition matrix Φ(t, t 0 ) = ∂x/∂x(t 0 ) (with x = [x M ; x 1 (t 1 ); ...; x 1 (t 104 ); x 2 (t 1 ); ...; x 2 (t 104 )]) in such a way that the coupling terms ∂x i (t)/∂x M (t 0 ) are referenced to the correct arc, so that for a given arc j: with δ ij the Kronecker delta function. By doing so, the estimation of spacecraft state and natural body state can be done concurrently, as the influence of a change of natural body state is directly mapped to a change in (barycentric) spacecraft state. This approach is in contrast to the typical approach of orbit determination and ephemeris generation, where the spacecraft's orbit is first estimated using Doppler data only, and the planetary ephemerides are then estimated using range/VLBI data (without adjusting the spacecraft orbit). In this traditional approach, the direct coupling between the planetary and spacecraft orbits is omitted. When incorporating ILR, however, the laser data will have significant and useful information on the dynamics of both natural and artificial bodies, requiring the coupling to be incorporated into the simulations. Except for the modification in the computation of Φ(t, t 0 ) as in Eq. (15), our covariance analysis follows that outlined by, e.g., Montenbruck and Gill (2000).
The list of parameters we estimate for both scenarios is given in Table 1. The relevance of the parameters in the context of Mars missions is discussed by e.g, Konopliv et al. (2011);Rivoldini et al. (2011)). The notation (C, S) l,m is used as a shorthand for the combination of spherical harmonic coefficients C lm and S lm of Mars. For the Sun, we consider only the degree-two zonal coefficient J 2, . The Love numbers are denoted as k lm , and β denotes the PPN parameter (Will, 2014). Our simulations represent a reduced analysis, not a full mission analysis. Instead, it is geared towards investigating the contribution of the range and Doppler data types to the estimation of various parameters, and comparing the results to those obtained from Section 4.1.
For both data types, we use a daily tracking pass of 2 hour duration, generating one independent measure-  ment every 60 s. The laser ranging data is weighted at 1 cm, and the Doppler data at 0.01 mm/s (see Section 6.1 for more detailed discussion). We do not simulate data for a solar separation angle smaller than 5 • . For our simulations, we use the Tudat software 4 .

Results
Using the methods outlined in Section 4, we give the results of our analytical and numerical analyses.

Analytical Results
Here, we present the results of the analysis method described in Section 4.1. For current radiometric systems, we use accuracies of 1 m in range and 0.04 mm/s in range rate for ∆t i =60 s. For upcoming state-of-theart systems, we assume 0.2 m accuracy in radiometric range, and both 0.01 mm/s at 60 s integration time and 0.002 mm/s at 1000 s integration time (Section 2.1). For ILR, we assume an upper and lower bound on data quality of 2 and 10 mm, respectively, a broad range which we derive from our discussion in Section 3.4. Evaluating Eqs. (11) and (12) to approximate the behaviour of the partial derivatives gives the results for the SNR h,q from Eq. (9), shown in Fig. 2.
The figure shows that the dual-frequency Doppler data at ∆t i =60 s (σ=0.01 mm/s) and the ILR curves cross in the area of 0.33-1.65 hours (for range data precision of 2-10 mm), and therefore Ξ q < 1 for lower T and Ξ q > 1 for larger T . This time interval is a particularly interesting one, as it contains the orbital period of many spacecraft orbiting rocky/icy bodies. For the ∆t i =1000 s Doppler data, this time interval is shifted by a factor of 5 (as the 1000 s data has a precision 5 times better than the 60 s data). For cases where Fig. 3: Histograms of formal error ratio s / ṡ for all arcs of the orbiter position component estimation, generated using the settings described in Section 4.2. Vertical black lines represent the analytical ratios Ξ q at the spacecraft orbital periods.
T < ∆t i , however, a 2π ambiguity arises in the estimation, requiring the use of a smaller ∆t i . Doppler data for orbiters typically have an integration time of 60 s or smaller. Reducing the integration time of the data does result in greater data volume, reducing the formal estimation errors if the noise is not correlated in time.

Numerical Results
In this section, we show the results of the covariance analyses described in Section 4.2. These results are shown in Figs. 3-6. For these results, we use the formal error ratio q,s / q,ṡ as figure of merit, instead of Ξ q from Eq. (14). Formal errors are a more robust indicator of the relative strength of the observables, as it includes the difference in correlations when performing Doppleronly and range-only estimation. In the analytical analysis of Section 5.1, formal errors cannot be obtained, as no estimation is performed. In the remainder of this section, we analyze how well the two figures of merit compare to one another. The analytical approach is a reasonable approximation if Ξ q ≈ q,s / q,ṡ .
The orbiter arc initial position error ratios ri,0,s / ri,0,ṡ are shown in Fig. 3. The mean ratio (≈ 1.0) is slightly larger than the theoretical ratio Ξ q (≈ 0.95, see Fig. 2, with settings from Section 4.2). Deviations from the analytical value are primarily due to the correlations between the various parameters, which are slightly worse for the range-only case than for the Doppler-only case.
The ratio r M ,s / r M ,ṡ (Mars initial state) for the orbiter and lander simulations is shown in Figs. 4a and 4b, respectively. The theoretical ratios for several relevant periods (as taken from Fig. 2) are also indicated. For the lander simulations, the theoretical ratio is close to the value expected from a periodic signal at the Mars and Earth orbital frequencies. For the orbiter simulations, however, the ratios are a factor 5 higher than these analytical values (Fig. 4a). The correlations between x M (t 0 ) and the other parameters are low for both data types. Since the analytical value for r M ,s / r M ,ṡ is much more closely attained in the lander estimation (Fig. 4b), the discrepancy in the case of orbits indicates that the once-per-orbit signature of the spacecraft continues to have a significant effect on the estimation of Mars' initial state when obtained through the use of Eq. (15).
Rotational properties and Love numbers of Mars as estimated from orbiter data both have formal error ratios of around 1.0 (Fig. 5a). This is comparable to the orbiter initial-state estimation ratios (see Fig. 3), and indicates that the primary signature of these parameters is derived from a once-per-orbit signature. The range data perform slightly worse than the analytical result would indicate, due to stronger correlations between the parameters.
The error ratio for the gravity field coefficients C l,m and S l,m is shown in Fig. 6. These results show a clear trend in their formal error ratios, with C l,l and S l,l having a ratio about twice that of the initial state r 0 , increasing to roughly (l + 1) times that of the initial state for C l,0 . Under the assumptions of the model outlined in Section 4.1, this would indicate that these parameters induce a periodic signature with a frequency of twice (for m = l) to (l + 1) times (for m = 0) the orbital frequency.
Since the orbits we have used are close to polar, the time-behaviour of the perturbation due to a given gravity field coefficient is mostly determined by ∇P l,m (sin φ), with P l,m associated Legendre polynomials at degree l and order m, and φ the body-fixed latitude of the spacecraft (e.g. Montenbruck and Gill, 2000). A spherical harmonic coefficient at this degree and order causes an acceleration on the spacecraft that is linearly proportional to P l+1,m+1 , P l+1,m and P l+1,m−1 . For P l,m , the number of zero crossings over a full orbit is 2 for l = m, while it is l for m = 0 and m = 1.
Our results in Fig. 6 show that the influence of a function P l,m may be simplified to a sine function with the same number of zero crossings, for the purposes of our analytical comparison. This allows the analytical criterion in Fig. 2 to be used for approximating the relative contribution of range and Doppler data to gravity field estimation.   The error ratios of the parameters estimated in the lander scenario are shown in Fig. 5b. They clearly fall into two categories: those for which the formal error ratio is close to that predicted for signatures at Mars' rotational period (dotted black line; q,s / q,ṡ ≈ 0.07), and those for Mars' orbital period (full black line; q,s / q,ṡ ≈ 10 −4 ). All of Mars' rotational properties fall into the former category, as expected. The parameters impact-ing Mars' orbit (β and J 2 ) fall into the latter category, but their formal error ratio is about a factor two smaller than expected from our analytical approximation (assuming the signal to be at Mars' orbital period). This difference is due to the correlation with the Mars initial velocity parameters, which is < 0.2 for both β and J 2 in the range-only simulations, and 0.6-0.7 for the Doppler-only simulations. Finally, one of the compo-

Discussion
Having shown the results of our analytical and numerical analyses in Section 5, we now discuss the potential application of ILR in planetary missions.

Observation Uncertainty Comparison
Section 5.1 gave the results of a conceptual criterion for comparing range and Doppler data using an analytical formulation. This was shown in Section 5.2 to approximate the results of the numerical covariance analysis reasonably well. However, both the analytical method and the covariance analysis suffers from common limitations: they ignore any differences in the measurement uncertainty probability distributions, and dynamical model errors are neglected. Doppler tracking is close to being bias-free, and in the case of dual-frequency tracking, it has a noise spectrum that is close to Gaussian (Section 2.1). For laser tracking, the single-shot uncertainty distribution will be defined by the Gaussian pulse profile, convoluted with the detector impulse response in the case of singlephoton detection (Section 3.1). Both the pulse pro-file and the detector response can be characterized to high accuracy (Section 3.2). However, as is the case in SLR/LLR, biases and instabilities at the several mmlevel will likely continue to be an issue (Section 3.4).
The issue of stability will be especially significant for lander tracking, where model uncertainties are expected to be the dominant source of estimation error (Section 3.4). For orbiter tracking, the dynamical model error will be more significant (compared to landers) and the impact of the precise observation noise spectrum less pronounced. Dirkx et al. (2014b) analyzed the impact of unresolved constant ILR biases for the Phobos Laser Ranging (PLR; Turyshev et al., 2010) mission. They used 5 mm Gaussian measurement noise, and 5 mm constant unresolved bias, and obtained results indicating that the biases contribute about an order of magnitude more uncertainty to the estimated parameters than the Gaussian noise. For biases that vary quasi-randomly from pass to pass, this impact will in part average out. However, it will still cause the range-only simulations presented here to be more optimistic than the Doppleronly simulations, especially for the lander case. As such, any results from the criterion in Eq. (14) that indicate similar signatures on Doppler and range data should be interpreted as indicating Doppler data will likely continue to be the better choice of data type. Nevertheless, reducing system and observation biases in laser ranging systems is a continuous and ongoing priority in ILRS activities.

Gravity Fields and Orbit Determination
Radio-range measurements are poorly suited for estimating planetary gravity fields. The Doppler data are used as the primary input data type for current missions (e.g. Marty et al., 2009;Konopliv et al., 2011;Mazarico et al., 2014). The results in Section 5.1 indicate that the use of ILR becomes competitive for signatures with a period of 0.33-1.65 hours for range data precisions of 2-10 mm (corresponding to range data accuracy at approximately sub-mm to several-mm level, see Section 6.1). The analytical analysis was shown in Section 5.2 to provide a good approximation to the full numerical simulations for gravity field estimation (see Fig. 6). This indicates that Doppler data are the superior choice for gravity field estimation, with the possible exception of very low gravity field degrees. For low degrees (2-4), ILR could meaningfully complement the Doppler data, which would require exceptional stability of the range data. Since Doppler data are essentially bias-free, the competitive estimation of low-degree gravity fields will only be achievable by means of ILR if it is similarly (close to) bias-free (Section 6.1).
The comparatively poor sensitivity of ILR data to gravity field coefficients shows that a laser-only tracking system is unlikely to be a suitable design choice for a planetary orbiter. Any uncertainty in the target body gravity field will propagate into an increased error in the orbit determination of the spacecraft (e.g. Mazarico et al., 2012). Even if ILR is used on an orbiter for other applications, it should always be in tandem with a Doppler tracking system, which can be used to measure the high-frequency variations in the spacecraft dynamics, chiefly the influence of the target body's gravity field. One possible exception is in the case of lunar missions. The gravity field of the Moon has been determined to such extreme accuracy  that the impact of its uncertainty on a typical spacecraft's orbit determination will likely be negligible.

Rotational and Tidal Characteristics
For tidal parameters, the signature on orbiter dynamics will show a combination of the spacecraft's orbital period, the rotation rate of the body being orbited and the period of the tidal forcing. The results in Fig. 5a show that it mainly is the once-per-spacecraft-orbit signature that is dominant in determining the formal error ratio s / ṡ .
Similarly, rotational variations such as librations can have a broad range of periods, but significant variations typically do not have a period that is much smaller than the rotational period (e.g., Konopliv et al., 2006;Petit et al., 2010). As with tidal parameters, our results in Fig. 5a indicate that it is again the once-per-spacecraftorbit signature that is dominant in determining s / ṡ . For most tidal and rotational parameters, the analytical approximation (≈0.95) slightly overestimates the strength of the range data, as Fig. 5a shows formal error ratios ≈ 0.9−1.25 (note that a higher ratio indicates a relatively weaker contribution of the range data).
This indicates that orbiter ILR data could be used to estimate tidal and rotational characteristics at a level that is comparable to that from Doppler data. Doing so will require a low level of unresolved systematic error in the range data, at the several mm level, assuming that (non-conservative) force modelling on the spacecraft does not limit the estimation quality. Moreover, as discussed in Section 6.2, the extraction of signals from the dynamics of orbiters will require the inclusion of a Doppler system, to prevent the uncertainty in short-periodic perturbations from degrading the quality of the estimation results. As a result, ILR could be used to supplement Doppler tracking in the determination of rotational and tidal characteristics from orbiter dynamics, but would not be the optimal choice as a dedicated system.
Estimation of rotational properties from lander data is shown in Fig. 5b. Our analytical model closely predicts the s / ṡ ratio, at the once-per-Martian-day frequency. Compared to orbiters, lander missions are more favorable for ILR, as the rotation period of a body (e.g. its day) is typically longer than the orbital period of a spacecraft.
Considering the typical rotational periods of bodies in our solar system, we can confidently state that laser range measurements to landers will be better suited for the estimation of rotational parameters than Doppler measurements will be. Fig. 5b shows a factor 20 improvement of Mars rotational parameters from ILR compared to Doppler data. For some fast-rotating bodies such as Phobos, this factor would be reduced to <10. For bodies with much slower rotational periods, such as Ganymede and Mercury, the formal error from ILR could be >100 and >1,000 times smaller than from Doppler data. However, at these levels of observational accuracy, many models will need to be improved to make full use of the data quality that would be available (Dirkx et al., 2014b).
Nevertheless, our results clearly show the exceptional strength that ILR can have in characterizing rotational motion and tidal deformation in the solar system, especially for landers. Similarly, it will be well suited for the determination of the h 2 and l 2 Love numbers (and possibly k 2 , due to its influence on rotational dynamics (e.g. Williams et al., 2001)). Fig. 4 indicates that ILR is preferred over Doppler data for the determination of planetary ephemerides, which is unsurprising considering the current role of radio range data (Section 1). The relative contribution of range and Doppler data to the estimation of ephemerides is well approximated by the criterion in Eq. (14) for the case of lander missions. For orbiter data, Eq. (14) overestimates the relative contribution of the range data by a factor of up to 5. Clearly, the Doppler data continues to contribute to the determination of planetary ephemerides, by accurately extracting the orbiter dynamics from the observations. As a consequence, ephemeris determination will benefit from the combination of Doppler and ILR data, for reasons discussed in Section 6.2. Fig.  5b shows that the parameters estimated jointly with ephemerides (here only β and J 2, ) will benefit greatly from the use of ILR. These parameter are crucial in relativistic experiments (Will, 2014). For the case of landers, their uncertainty is reasonably approximated by the analytical formulation. However, properly decorrelating these parameters using ILR data may require laser data to multiple targets.

Solar System Ephemerides
Due to the scarcity of ILR data when it will be first implemented, there will be an imbalance of several orders of magnitude between the measurements used for creating solar system ephemerides. For instance, when using ILR data from an Earth-Mars link, there will be mm-level range measurements for Earth and Mars, mlevel (radiometric) range measurements for other solar system bodies at which orbiter/lander tracking data is available and km-level range measurements (radar) to bodies where no data from spacecraft tracking techniques is available. For many small bodies, no range data will be available at all, and ephemeris generation must be performed from astrometric data alone. This effect, as well as dynamical model error (Section 3.3) will degrade the fidelity of the orbit estimation of the bodies between which an ILR link is set up. Quantifying the exact requirements for a mission profile, tracking schedule, estimation settings, etc. for optimally exploiting ILR data in ephemeris generation will require a dedicated study.

Conclusion -The Science Case for ILR
The purpose of the article has been twofold. Firstly, we have given a detailed overview of the sources of uncertainty in both the realization and analysis of ILR data (Section 3). Secondly, we have compared the performance of ILR and radio Doppler data both analytically and numerically, to clarify the science case for ILR (Sections 4-6).
Mm-precision normal points are feasible for ILR. Sub-cm accuracy of the data will be attainable, but reaching the mm-level accuracy is hindered by both measurement instabilities and model uncertainties, similar to SLR/LLR (Section 3.4).
We have derived an analytical approximation of the sensitivity of ILR and radio Doppler data types, which is shown in Fig. 2, under the assumption of sinusoidal signatures on the data. This figure indicates that the signatures with a period of 0.33-1.65 hours can be observed in the ILR and radio Doppler data at a similar signal-to-noise level. The results of our numerical covariance analysis largely validate the analytical approach, allowing Fig. 2 to be used as a conceptual design tool.
However, instabilities in the range data accuracy, which are not directly included in the simulations, will limit the performance of ILR data to the upper bounds of the 0.33-1.65 hour range. That is, ILR will start to be competitive for determining signatures of >1.5-2 hours. For effects with a longer period than several hours, ILR data unambiguously provide a more accurate estimation, while Doppler data will continue to be the optimal data type for short-periodic effects.
ILR's weak sensitivity to short-periodic effects makes it a poor choice for gravity field estimation, with the possible exception of low-degree coefficients. Any orbiter with typical orbit determination requirements will continue to require Doppler data (with the possible exception of lunar missions). For typical mission profiles, ILR and Doppler data have a similar sensitivity to onceper-orbit effects (period of 1.5-2 hours). Therefore, ILR will be valuable in complementing the Doppler data for the estimation of tidal and rotational characteristics from orbiter tracking. As expected, ILR is exceptionally well suited to generating ephemerides, although our analytical approximation overestimates its strength by a factor 5 for the case of orbiter tracking.
In addition, when used for lander tracking, ILR will be an excellent method for the estimation of both tidal and rotational characteristics of the target body. For a Mars lander, ILR data produces estimates that are about a factor 10 more accurate than the estimation from Doppler data. The strength of ILR improves even more for bodies with slower rotation rates.
We have focused on the estimation quality of orbits and geodetic parameters, concluding that the science case for landers is excellent, and it could serve a complementary role for orbiters. Owing to the highly accurate data that ILR will deliver, improvements in the various models entering not only the analysis, but also the interpretation of the estimation results, must be brought to a level where all data can be used to their full potential. Not only will this require significant theoretical effort, but it implies that more accurate knowledge must be obtained of various quantities that cannot be obtained from tracking data alone. Examples of synergistic data are magnetic field, heat flow, geological and seismic measurements, which will be important for the full characterization of a body's interior structure and composition. By combining these data from nextgeneration space missions with ILR, the full set of these measurements can be exploited to their full potential, allowing the study of planetary interiors to be brought to the next level. Fig. 7: Ratio of figures of merit Ξ r0 (e)/Ξ r0 (e = 0) for elliptical orbits and circular orbits as a function of viewing geometry. The color scale denotes this ratio. Note that the periapsis distance r p is set equal in the e = 0 and e = 0 simulations. Doppler data, and to three anonymous reviewers, whose comments improved the clarity, conciseness and completeness of the manuscript. Part of this work was performed in the FP7 ESPaCE project, financially supported by the EC FP7 Grant Agreement 263466.

A Semi-analytical Approach -Eccentric orbits
A key exception for which the assumptions of Section 4.1 will not hold is when a spacecraft orbits with a substantial eccentricity, e.g., Juno (e ≈ 0.95 w.r.t. Jupiter), Messenger (e ≈ 0.7 w.r.t. Mercury), Mars Express (e ≈ 0.57 w.r.t. Mars). To extend the method of Section 4.1 to non-spherical orbits, we continue to use the criterion of Eq. (9), but compute the partial derivatives from the orbits directly using the method of (Moyer, 2000), instead of imposing them to behave sinusoidically. We limit ourselves to taking q as the initial position r 0 of the spacecraft w.r.t. the planet it is orbiting.
We vary the eccentricity from 0 to 1 and assess the influence on the behaviour of Ξ q . For an eccentric orbit, the influence of the geometry of the orbital plane w.r.t. the observation line-of-sight needs to be analyzed. We parameterize this geometric dependency by two angles: the angle between the line-of-sight vector and the spacecraft orbital plane, and the angle between the line-of-sight vector and the spacecraft orbital velocity vector at periapsis.
Using the methodology outlined above, we have generated values of Ξ q (with q the initial position r 0 , and computing the associated Ξ r 0 as the root sum square of the constituent vector components of ∂h/∂q) for elliptical Kepler orbits under a full range of observational geometries w.r.t. the observer line-of-sight.
As a test case, we use a Mars orbiter for our analysis. When generating the values of Ξ r 0 for two cases, one with e = 0 and one with e > 0, we find that using the same spacecraft semi-major axis for both does not lead to insightful results. Instead, when using the same spacecraft periapsis distance r p = a(1 − e) for the two cases, the results for zero and nonzero eccentricities can be related much more intuitively. We show the results of this analysis in Fig. 7, were we plot the ratios of Ξ r 0 (e)/Ξ r 0 (e = 0), using the same r p to compute the two values of Ξ r 0 .
From this figure, it can be seen that even for substantial eccentricities, the ratio Ξ r 0 (e)/Ξ r 0 (e = 0) remains close to 1 for a broad range of eccentricities and observational geometries. In our simulations (e ≤ 0.9) , the results deviate from the analytical ratio by less than 50 % in all cases and less than 10 % in most cases. This indicates that the analytical criterion shown in Fig. 2 continues to be largely applicable even for large eccentricities. This result greatly simplifies the first-order analysis of the added value of an ILR system, as it allows the analytical results for spherical orbits to be used for arbitrary eccentricities.
When using Fig. 2 to analyze the comparative strength of range and Doppler data for elliptical orbits with period ω P , the following scaling should be applied to obtain the value of ω at which Fig. 2 is to be read: ω = (1 − e) −1.5 ω P so that for elliptical orbits the Doppler data continues to be the preferred method over ILR for a larger range of values of ω P .