1 Introduction

Since the first demonstration of Satellite Laser Ranging (SLR) in 1964, its accuracy has increased by several orders of magnitude thanks to progress in laser technology, increased computational power of modern computers and the development of more refined models, approaching the current accuracy of few millimeters (Pearlman et al. 2019). Despite its conceptual simplicity, laser ranging technique has implications in several fields from fundamental science, technology and everyday life. Improvement in laser sources (Coddington et al. 2009), and consequently in laser ranging performances, can be beneficial for fundamental physics test ranging from Lense–Thirring effect (Ciufolini et al. 1998, 2019b; Ciufolini and Pavlis 2004; Ciufolini 2007), weak equivalence principle (Ciufolini et al. 2019a), gravitomagnetic effect (Nordtvedt 1988), non-Newtonian gravity (Lucchesi and Peron 2010), dark matter (Yang et al. 2016) and satellite quantum communications (Vallone et al. 2015, 2016; Vedovato et al. 2017; Agnesi et al. 2018; Calderaro et al. 2018). SLR is also a fundamental tool for geodesy and geodynamics, allowing precise length-of-day (LOD) determination (Langley et al. 1981), earth shape and gravitational field measurements (Kaula 1983; Tapley et al. 1985; Cheng et al. 1997; Maier et al. 2012), Earth tidal measurements (Smith et al. 1973; Cheng et al. 1992) and development of the terrestrial reference frame (TRF) (Seeber 2003). Moreover, other fundamental applications of SLR that are the determination of the geocenter coordinates, the determination of the polar motion, the determination of the SLR station coordinates and the determination of the product between the gravitational constant and the Earth mass Pearlman et al. (2019).

The high accuracy of SLR is necessary also for supporting several mission, such as altimeter calibrations (Exertier et al. 2001; Luthcke et al. 2003; Rim et al. 1999) and GNSS positioning (Urschl et al. 2007; Sośnica et al. 2015). Furthermore, by exploiting precise SLR measurements, in the near future it will be possible to perform time and frequency transfer through satellite reaching picosecond level (Samain et al. 2015; Kucharski et al. 2019).

To face the new challenges of modern Earth sciences, such as the precise monitoring of the sea level, the global geodetic observing system (GGOS) sets the accuracy requirements of the terrestrial reference frame at a level of 1 mm with a stability of 0.1 mm/year (Plag and Pearlman 2014). This requirement directly implies an increased precision of the SLR data, since it is one of the fundamental techniques for TRF realization. While the precision of a single measurement of the satellite range is limited by several fundamental and technical reasons, such as the instrumental jitter, retroreflector shape, laser pulse width and atmospheric dispersion, it is possible to increase the accuracy of the satellite orbit determination with an higher number of ranging measurements, i.e., by increasing the system repetition rate. For this reason, many stations of the international laser ranging service network (ILRS) moved from traditional low repetition rate (10 Hz) to high repetition rate regime, in the range of few kHz. The Graz SLR station, for instance, has experimented rates up to 10 kHz, with short laser pulses (Kirchner et al. 2011). Even if the single-shot precision is generally lower for high repetition rate systems, the increase statistics compensates for this factor resulting in a high accuracy obit determination (ILRS 2020). Moreover, the introduction of kHz lasers for SLR led to an increase of data rate and enabled the attitude and spin rate measurement for several satellites and space debris, showing how an increase of the repetition rate can augment the potential of SLR technique (Kucharski et al. 2008, 2010, 2013, 2017; Calderaro et al. 2018; Steindorfer et al. 2019). Finally, the use of high repetition rate and single photon detection allowed the determination of the attitude of the Ajisai satellite, thus removing the satellite signature effect, which is the main contribution of range uncertainty (Kucharski et al. 2015). In Stuttgart, an even higher repetition rate of 100 kHz has been tested (Hampf et al. 2019). However, the long laser pulse duration of 10 ns strongly affected the single-shot precision, preventing spin and attitude determination.

2 Setup description

SLR is a technique that exploits the echoes of short laser pulses reflected by satellites to precisely measure their distance. For this purpose, MLRO uses a 10 Hz train of laser pulses at 532 nm with \(\approx 100\) mJ of energy per pulse and \(\approx 50\) ps pulse width. These pulses are generated using a mode-locking Nd:YVO\(_4\) laser oscillator (ML-Laser), operating at 1064 nm with 100 MHz repetition rate and paced by an atomic clock. One pulse every \(10^7\) is selected with a pulse picker regenerative amplified (PPRA). The 10 Hz pulse train then passes through two single-pass amplifiers (SPA) and is up-converted via a second-harmonic-generation (SHG) stage to the working wavelength of 532 nm. A fast photodiode detects the amplified pulse and generates the start signal. The 10 Hz pulses are sent to satellites equipped with corner-cube retroreflectors (CCRs), using the diffraction-limited Cassegrain telescope of MLRO with a diameter of 1.5 m. After reflection by CCRs, the SLR laser pulses are collected by the MLRO telescope and are detected using a Hamamatsu R5916U-50 fast analog micro-channel plate detector (PMT) placed after a polarizing beam splitter (PBS) used to separate the transmitted beam from the received one. These detection events are the stop signals of the SLR technique. A dedicated time-to-digital converter (TDC) with picosecond accuracy records the start and stop signals. The single-shot measurement of the satellite distance is then estimated from the time difference of these two signals, i.e., the round-trip time.

To demonstrate the possibility of SLR at higher repetition rates, we complemented the 10 Hz MLRO system with an 100 kHz system that worked alongside. The laser pulses were generated by the Ekspla Atlantic 60 high power industrial picosecond laser. The pulses generated by this laser exhibited a 1064 nm wavelength and \(\approx 10\) ps pulse width. The 100 kHz pulse train was then converted to 532 nm wavelength via a SHG stage directly attached to the laser encasing. After the wavelength conversion, each laser pulses had an energy of \(\approx 400\,\upmu \hbox {J}\).

Beam shaping optics were then placed after the Ekspla laser, allowing us to match the beam width and divergence of the 100 kHz pulse train with those of the 10 Hz train. Then, the 100 kHz pulse train was directed towards the Coudé path using dielectric mirrors, where, before entering, it was combined with the 10 Hz laser train using a 50:50 beam splitter (BS). As a consequence, both the transmitted and received intensities are halved. It is worth mentioning that in a system dedicated to 100 kHz SLR (i.e., without 10 Hz as reference) this \(\times 1/4\) factor in signal reduction would not be present. Behind one dielectric mirror, a Thorlabs DET025AL/M silicon photodiode (2 GHz bandwidth, 150 ps rise time) was placed. This photodiode generated the start signal of the 100 kHz laser exploiting a weak transmission from the mirror. To avoid overlap of outgoing and incoming laser pulses we operated in “burst mode”: for each 100 ms interval, the laser has been firing for 25 ms (i.e., longer than the round trip time for all LEO satellites). Then, the laser was turned off and the system switched to detection phase, enabling the detector in free-running (i.e., without further gating). For LAGEOS satellites, the firing interval was extended to 50 ms. With this system, it is possible to remove the noise from back-scattering; however, a maximum of \(50\%\) of duty cycle can be obtained.

The 100 kHz pulse train, together with the 10 Hz standard SLR train, was then directed to the targeted satellites by the MLRO telescope. As in the standard SLR procedure, the pulse trains are then reflected by the orbiting terminals and collected by the MLRO telescope. The receiving apparatus of the 100 kHz beam is comprised of a 50:50 BS to separate the outgoing and incoming beams, a 3 nm FWHM spectral filter with transmission band centered at 532 nm, a focusing lens and a silicon single photon avalanche detector (SPAD), PD-200-CTX provided by Micro-Photon-Devices Srl (Giudice et al. 2007), with \(\approx 50\%\) quantum efficiency, \(\approx 400\) Hz dark count rate and 40 ps of jitter, full width at half maximum (FWHM). These detection events are the stop signals of our 100 kHz SLR system. The start and stop signals are recorded by the quTAG TDC from qutools GmbH with 1 ps resolution and 10 ps RMS interchannel jitter. Also in this case, the 50:50 BS halves both the transmitted and received intensities.

Fig. 1
figure 1

Experimental setup of the 100 kHz SLR system working alongside the standard 10 Hz system. For a detailed explanation, see Sect. 2

Full size image

The 100 kHz system was tested on several ground target CCRs placed at a distance from the telescope ranging from 45 m to 190 m. These acquisitions have been repeated periodically both to calibrate the constant delays introduced by electronics, cables and optical path and to check the stability of the system. As the single CCRs reflects the signal without introducing dispersion, this system can be used to measure the minimum jitter introduced by all the elements in the acquisition chain: laser pulse duration, start and stop detector jitter, TDC interchannel jitter. The expected time response of the detection system is Gaussian peak followed by an exponential tail (Agnesi et al. 2019):

$$\begin{aligned} f(t)=A \mathrm {e}^{-\frac{\left( t-t_{0}\right) ^{2}}{2 \sigma ^{2}}} \varTheta \left( t_{1}-t\right) +B \mathrm {e}^{-\frac{t-t_{1}}{\tau }} \varTheta \left( t-t_{1}\right) \end{aligned}$$
(1)

where \( \varTheta \) is the Heaviside function, \(t_0\) is the peak position of the Gaussian part and \(t_1\) is the time of the crossing between Gaussian and exponential behavior and \(B=A\mathrm {e}^{-\frac{\left( t_{1}-t_{0}\right) ^{2}}{2 \sigma ^{2}}}\) to ensure continuity. This pulse shape is due to the electronic response of the single photon detector, which in our system is the main cause of timing uncertainty. During calibrations we observed that the value of the parameters \(\sigma \) and \(\tau \) have a dependence with the pulse detection probability, when the latter exceeds 0.2. However, in standard SLR operation the mean number of photon per pulse is often below 0.2, due to the high attenuation of the transmission channel. In this case, the measured value of the Gaussian standard deviation is \(\sigma =20.8\pm 4.7\) ps and the exponential constant is \(\tau =196\pm 17\) ps. The measured Gaussian standard deviation matches the theoretical one \(\sigma _{\mathrm {theory}}\), obtained from the detector jitter \(\sigma _{\mathrm {det}}=17\) ps, the electronic jitter \(\sigma _{\mathrm {elec}}=10\) ps and the pulse length \(\sigma _{\mathrm {laser}}=3.8\) ps, giving \(\sigma _{\mathrm {theory}}= \sqrt{\sigma _{\mathrm {det}}^2+\sigma _{\mathrm {elec}}^2+\sigma _{\mathrm {laser}}^2}=20.1\) ps.

3 Data analysis and results

The data from our 100 kHz SLR system were acquired from January 22 to January 26 2018, after a few days of setup optimization and characterization. The data acquisition was limited to night-time only operation, due to the lack of a narrow spectral filter to reduce daylight spurious counts. However, daylight single photon SLR is a well-established technique and there is no intrinsic limitation related to high rate systems in nighttime, while for daytime operation a detector gating could be required to avoid saturation of the data acquisition system.

During the four nights of data acquisition, 12 satellites were tracked, including several LEO and the two LAGEOS, producing 726 normal points overall. It is worth underlining that all the satellites tracked and acquired by MLRO also showed returns with the 100 kHz setup, demonstrating the reliability of the technique. For the data analysis we proceeded according to standard ILRS procedure, as explained in Sinclair (1997). As an example, we show the complete analysis of a single pass of the Jason-3 satellite, whose tracking started on January 24 at 2:34 UTC and lasted slightly less than 5 min. The maximum elevation was of \({44}^\circ \) corresponding to a minimum round-trip time of about 11 ms.

As a first step we computed the prediction residuals as the difference between the actual detection time stamp and the expected time of arrival, computed from the predicted orbits available in the ILRS database. As shown in Fig. 2, the signal-to-noise ratio is excellent, thus allowing a simple and clear distinction of the echoes from the uncorrelated background. The filtering algorithm first performs a two-dimensional histogram with a 1 second spacing along the elapsed time (x-axis of Fig. 2) and 1 ns spacing along the prediction residual (y-axis of Fig. 2). For each time interval (x-axis), the highest values of this histogram are then selected, eliminating only the bins with less than 10 returns. The selected bins are then fit with a 7-order polynomial. Finally, the returns are identified as those for which the prediction residuals lie in a \([-\,500, 3000]\) ps interval with respect to the polynomial value. The noise contribution, which can be estimated as the average bin counts outside the selection interval, ranges from 0.1 to 0.4 counts per second of elapsed time per nanosecond of prediction residual and is much lower than the threshold of 10 returns per bin. The noise corresponds to a free-running dark count rate of \(\simeq 30\,\hbox {kHz}\), mostly due to laboratory illumination.

Fig. 2
figure 2

Prediction residual of a Jason-3 pass before (top) and after (bottom) data filtering

Full size image

After the removal of the noise signal, it is possible to estimate the average detection probability, i.e., the number of returns divided by the number of fired pulses, which ranges from \(0.15\%\) for LAGEOS-1 satellite to \(24\%\) for STARLETTE Fig. 3. It is worth noticing that for most passes the detection probability is less than \(10\%\), thus maintaining the system within single photon regime. On the basis of the LAGEOS detection probability, which ranges from 0.15 to \(1.05\%\), it is possible to estimate the efficiency expected for higher orbit satellites, such as GLONASS, Galileo and Etalon. As the pointing error for LAGEOS and higher orbit satellites is similar, the detection rate, DR, is proportional to the following parameters (Degnan 1993):

$$\begin{aligned} \hbox {DR} \propto \frac{\sigma }{R^4} \end{aligned}$$
(2)

where R is the satellite distance and \(\sigma \) is the satellite cross section, which has been taken from (Arnold 2002; Navarro-Reyes 2014). According to the variation of LAGEOS return rates, as well as range excursion of higher orbit satellites, it is possible to estimate a minimum and maximum detection efficiency. By multiplying it for the transmission rate and the \(50\%\) duty cycle, due to burst mode operation, it is then possible to obtain the overall detection rate, which is reported in Table 1. We also report the estimated detection rate for an optimized system, which removes the 50:50 BS used for merging the 100 kHz line with the 10 Hz system, thus giving a factor \(\times 4\) in return rate. Moreover, by replacing the 50:50 BS for separation of transmission and reception with a polarization beam splitter followed by a quarter wave plate, it is possible to gain a factor \(\times 2\) for depolarizing target and \(\times 4\) for polarization maintaining targets. We assumed depolarization for all targets, as a lower bound on expected rate.

Table 1 Estimation of return rate for several satellites in GNSS orbit
Full size table
Fig. 3
figure 3

Detection efficiency for all satellite passes acquired

Full size image

A polynomial fit of 15th degree is used to remove the trend of the prediction residuals, allowing the determination of the fit residuals defined as the difference of the prediction residuals and the polynomial fit. To remove large outliers, the data outside 2.5 standard deviation are edited and the fit is recalculated. This procedure has been repeated three times, giving a value that quickly converges. The overall dataset is then divided in intervals defined by the standard normal point interval (SNPI), which for Jason-3 corresponds to 15 seconds. In Fig. 4 we show the distribution of the fit residuals within a single SNPI. We used the detector response function to fit the data, obtaining a good agreement and a value of the Gaussian \(\sigma =32\) ps, slightly larger to the one obtained during calibration.

Fig. 4
figure 4

Histogram of the fit residuals within a standard normal point interval for a Jason-3 pass. Data are fitted according to the response of the detector

Full size image

The same analysis has been repeated for all the satellites acquired, computing the mean full-rate RMS for each pass. The results are shown in Fig. 5. The average RMS is limited to 50 ps for all satellites, with the exception of Ajisai and Beacon-C. This fact can be explained by the large target depth of these two satellites, related to their specific shape (Degnan 1993).

Fig. 5
figure 5

Two-way full-rate RMS averaged over a single pass for all satellites acquired

Full size image
Fig. 6
figure 6

Average number of echoes per SNPI for each satellite pass

Full size image

The overall quality of a laser ranging system is measured with the normal point (NP) RMS. According to the ILRS procedure (Sinclair 1997), for each SNPI a NP is generated by taking the experimental range measurement nearest to the mean epoch of the interval, subtracting the corresponding fit residual and finally adding the average of the fit residuals of the interval. The number of returns per SNPI is shown in Fig. 6. To check the validity of the system, we have been analyzing the passes of the two geodetic satellites LAGEOS with the NASA software GEODYN (Putney 1977), both with SNPI of 120 s. The acquisitions of the satellites LAGEOS 1 and 2 from all the ILRS ground stations (including the 100 kHz data) have been analyzed for reconstructing the satellite orbits in the week from 23rd to 30th of January. After the two orbits have been reconstructed on the basis of all data available, it has been possible to measure the difference of each NP from the reconstructed orbit, thus obtaining the orbit residuals. These orbit residuals, reported in Fig. 7, show how the precision of the new system is well within the performance of the ILRS stations, in particular for the LAGEOS-2 satellite, where we find a RMS of 1.3 mm for the pass acquired.

Fig. 7
figure 7

One-way orbit residuals for the two LAGEOS satellites in the week from 23rd to 30th of January 2018 for all ILRS stations. The MLRO station data, both with 10 Hz and 100 KHz setup, are highlighted

Full size image

For LAGEOS-1, the RMS of the orbit residuals is at the cm level, a value comparable to other stations. This is due to both range bias and time bias error, which in the global analysis are not optimized.

To get rid of time and range biases, we performed a linear fit on the NP residuals for each satellite pass. The residuals of this linear fit, shown in Fig. 8, allow us to calculate the NP-RMS, which is 2.0 mm, 1.3 mm and 1.2 mm respectively for the three passes analyzed. These values are close to the MLRO standard performances with 10 Hz laser, which is between 1.2 mm and 1.6 mm, depending on the orbital model used for the analysis (ILRS analysis centers 2019). While the results obtained perform well compared to the other ILRS stations, the NP-RMS is much higher than the expected standard error, namely the full-rate RMS divided by the square root of returns per normal point. The standard error is displayed in Fig. 8 as the blue bar for each NP.

Fig. 8
figure 8

One-way post-fit orbit residuals for the three passes of the LAGEOS satellites. The red shaded area corresponds to the NP-RMS while the blue bars are the theoretical errors, computed from the full-rate RMS divided by the square root of the number of returns

Full size image

4 Discussion

The results obtained during the test with the high repetition rate laser at MLRO have shown a very good single-shot performance for most satellites, with an RMS slightly higher than the one observed on calibration targets. The high number of returns detected in a single SNPI should bring an improved precision of the NP. However, we observe performance at the same level with respect to current stations precision, which does not seem to reflect the advantage of high repetition rate in SLR precision. In order to explain this behavior, we analyzed the stability of the system with ground target calibration runs, which were performed periodically. Similarly to satellite data, we calculated the average round-trip over 120-s intervals divided by two, as well as the RMS of these mean values, the analogous of the one-way NP-RMS. The results have shown in some cases an instability over long period, resulting in a maximum 0.3 mm NP-RMS, a value much higher than the one expected from the standard error, see Fig. 9 top. This value represents the worst case, while in all the other calibrations the NP-RMS was limited to \(< 0.15\) mm, and it reached 0.03 mm in some runs, see Fig. 9 bottom. In general, the instability of the system sets a minimum precision of the LAGEOS NP smaller than 0.3 mm and cannot explain the 1-2 mm NP-RMS measured on data.

Fig. 9
figure 9

Worst-case (top) and best-case (bottom) ground target analysis of system stability. The red shaded area corresponds to the one-way NP-RMS, while the blue bars are the theoretical errors, computed from the full-rate RMS divided by the square root of the number of returns

Full size image

Another explanation for the mm-level NP-RMS could be given by the distance of the satellite center of mass with respect to the optical back-reflected pulse shape. This effect, known as center of mass correction, is generally treated by assuming that the satellite is spinning and therefore within a SNPI the residuals will average any bias related to the satellite position. This turns out to be false for the two LAGEOS satellites, whose spin rate is much slower than the typical integration period. To confirm this fact, we analyzed the fit residuals for 8 consecutive SNPI of the LAGEOS-2 pass of January 25th. As shown in Fig. 10, the fit residuals within each SNPI show a characteristic shape given by the satellite geometry. The change in peak position and relative intensities directly reflects the rotation of the LAGEOS-2 satellite. In this situation, the center of mass correction is expected to change in each SNPI, and applying its mean value to all the pass will introduce a systematic error in the orbit determination. Due to the increased statistics, the use of 100 kHz laser could allow the attitude determination of the LAGEOS satellites, as done in Kucharski et al. (2015), and increase the attitude determination accuracy for GNSS satellites (Steindorfer et al. 2019).

Fig. 10
figure 10

Distribution of the fit residuals for 8 consecutive SNPI of the LAGEOS 2 pass of January 25th. It is possible to observe the change in retroreflected light distribution due to the rotation of the satellite

Full size image

Finally, another possible explanation for the limited NP precision of the high repetition rate system is the uncertainty of the time delay during propagation in atmosphere. This phenomenon, which is due to the variation of the refractive index in the troposphere, has been studied both for satellite and lunar laser ranging, demonstrating in both cases its importance for mm-level accuracy measures with a shot-to-shot variation of \(\simeq 0.5\,\hbox {mm}\) and long terms effects in the order of 1–2 mm (Currie and Prochazka 2014; Hulley and Pavlis 2007). Atmospheric effects can be divided in fast variations of the refraction index due to turbulence and in slow effects, related to the change in azimuth and elevation. While the first effect is averaged over a SNPI and therefore impacts mainly the full-rate RMS, the latter shows a slow dependence, which influences the NP-RMS. This slow effect is considered by SLR data analysis for what concerns the elevation dependency, while the azimuth dependence on refraction model used in the analysis is neglected. To have a better understanding of the potential impact of the repetition rate on the overall NP precision, we use an alternative method for determining the NP-RMS, which removes any slow drift in the system caused by instabilities, satellite signature effects and turbulence effects. It consists in estimating the RMS of the mean value of the fit residuals within each SNPI for a complete satellite pass. To do so, we first computed the polynomial fit residuals from full rate data. We then divided them according to SNPI and calculated the mean of the fit residuals within each SNPI. Finally, we calculated the RMS of these mean values. The results are shown in Fig. 11. The RMS values obtained with this method are much smaller than the ones obtained from the orbit reconstruction and might be better suited for estimating the intrinsic station performance, as they are not influenced by systematic errors. The one-way RMS obtained with the new procedure is lower than 1 mm for most satellites and reaches 0.19 mm, 0.14 mm and 0.07 mm for the three LAGEOS passes acquired, close to the theoretical limit of \(\frac{\hbox {RMS}_{\mathrm{FR}}}{\sqrt{N}}\), where \(\hbox {RMS}_{\mathrm{FR}}\) is the full-rate RMS, reported in Fig. 5, and N is the total number of return per interval.

The results shown in Fig. 11 demonstrate how the use of 100 kHz laser for SLR can improve the system precision by an order of magnitude with respect to the present limit, if only statistical errors are considered. Regarding the system accuracy, which is limited by systematic errors, we note that the use of multi-wavelength high repetition rate SLR, a more frequent calibration and the study of the satellite signature could disentangle the effects coming from system drift, satellite signature and atmospheric effect. This could help to better investigate each of the effects separately for reaching sub-mm accuracy satellite laser ranging.

5 Conclusions

In this work we demonstrated the first SLR operation with 100 kHz repetition rate and few tens of picosecond single-shot jitter. After a few days of alignment and optimization, the system has been fully operational for 4 nights, showing the reliability of the apparatus. The new SLR apparatus has shown a single-shot precision lower than 50 ps for most satellites, a value slightly higher than the one obtained during calibration on ground targets. We observed return rates ranging from 0.15 to 24 kHz, thus at least one order of magnitude higher than the current ILRS kHz stations. Despite the good single-shot accuracy and the high statistics, the NP-RMS calculated from orbit reconstruction is limited to the mm level for both LAGEOS satellites. This limit can be ascribed to several error sources, such as atmospheric effects, satellite signature and system instabilities. To estimate the intrinsic system precision, given by the statistical error, we used an alternative method for estimating the NP-RMS, based on polynomial fit of full rate data. This method showed a sub-mm precision for most satellites, and an error of \( \sim 100 \,\upmu \hbox {m}\) for the LAGEOS satellites. Finally, by analyzing the returns from the two LAGEOS satellites, we observed a clear peak structure, due to the geometrical disposition of the CCRs. The use of high repetition rate laser could provide enough statistics to determine the instantaneous orientation of the satellite, thus eliminating the error on the center of mass bias and allowing the measurement of spin rate also for satellites with extremely long rotation periods, such as the two LAGEOS satellites.

Fig. 11
figure 11

One-way mean RMS of the residuals of the NP with respect to the polynomial fit performed on full-rate data

Full size image