Upgraded modelling for the determination of centre of mass corrections of geodetic SLR satellites: impact on key parameters of the terrestrial reference frame

Journal of Geodesy https://doi.org/10.1007/s00190-019-01315-0 ORIGINAL ARTICLE Upgraded modelling for the determination of centre of mass corrections of geodetic SLR satellites: impact on key parameters of the terrestrial reference frame José Rodríguez1 · Graham Appleby1 · Toshimichi Otsubo2 Received: 22 March 2019 / Accepted: 22 October 2019 © Springer-Verlag GmbH Germany, part of Springer Nature 2019 Abstract We present improved modelling to compute centre of mass (CoM) offsets for spherical geodetic satellites and derive new values for the spacecraft LAGEOS, LAGEOS-2, Etalon-1/2, Starlette, Stella, LARES, and Ajisai for all system configurations of the laser ranging network, past and present. The updated CoM values, i.e. the average distances between the centres of mass and the optical reflection points, reveal that a major contributor to the error budget in satellite laser ranging (SLR) for reference frame realisation has been a limited accuracy of the CoM corrections hitherto employed to model some of the systems of the SLR network. Our modelling strategy takes into account target characteristics, detection hardware, and mode of operation. For the optical modelling, we have employed high precision data from the single-photon, high-repetition rate SLR system at Herstmonceux (UK), performing the relevant computations for all laser frequencies ever used for tracking operations in the SLR network. For stations that do not operate in single-photon mode, we performed Monte Carlo simulations to estimate the CoM offsets, in this way taking into consideration previously neglected aspects of the detection process. The newly derived values reduce by up to 50% the biases estimated for the two biggest targets, Etalon and Ajisai, and eliminate the generally positive trend in range biases previously observed across the network. These CoM values, mostly lower than the previous ones, imply higher estimated station heights, with a consequent increase in the scale of the terrestrial reference frame that is realised by SLR on the basis of LAGEOS and Etalon observations. This increase in scale, of about 0.65 ppb, is confirmed by geodetic solutions reported here using nineteen years of observations from the ILRS network of tracking stations. We also examine the consequences for the determination of the geocentric gravitational constant, GM, that result from using the new CoM corrections. Keywords SLR · ITRF · LAGEOS · Etalon · Starlette · LARES · Ajisai · GM 1 Introduction Satellite laser ranging (SLR) is in principle perhaps the most straightforward of the three satellite-based space geodetic techniques currently contributing to the realisation of the global terrestrial reference frame (SLR, Doppler Orbitography Integrated by Satellite, and Global Navigation Satellite Systems). This is true at several levels, from the simple, elegant concept of the technique to the technology underpinning B José Rodríguez josrod@nerc.ac.uk 1 NERC Space Geodesy Facility, Herstmonceux Castle, Hailsham, East Sussex BN27 1RN, UK 2 Geoscience Laboratory, Hitotsubashi University, 2–1 Naka, Kunitachi, Tokyo 186–8601, Japan its realisation, as well as for the comparatively modest complexity of the analysis required to exploit its measurements. Of course, this directness does not free the SLR technique from the presence of an array of potential pitfalls that may compromise its accuracy. Among the various design limitations, operational issues, and modelling problems that may arise, of paramount importance are an accurate determination of the offset between the centre of mass of the tracked spacecraft and the laser reflection points on their retroreflector arrays. For spherical geodetic satellites (e.g. Starlette, LAGEOS), designed with very low area-to-mass ratios and with a favourable symmetry for modelling non-conservative forces, the centre of mass offset is simply a scalar relating the optical and the mass centres. However, the arrangement of the retroreflectors over the spherical surfaces of these satellites causes a spread of the reflected optical pulses used to per- 123 J. Rodríguez et al. form the observations (Appleby 1992), the correct modelling of which, as well as that of the detection systems employed to measure them, is essential to derive accurate centre of mass values. The computation of the optical responses of spherical satellites has been tackled previously for several geodetic targets (Otsubo and Appleby 2003; Otsubo et al. 2014a), and CoM values derived for different kinds of laser ranging systems. However, the observational data employed at certain stages of those computations, at least in the case of the LAGEOS, Etalon, and Ajisai satellites, were of a previous generation of ranging systems, with the consequent poorer precision and much reduced data volumes afforded by the traditional low-repetition rate lasers. Additionally, CoM modelling for stations operating at multi-photoelectron levels of detection was never satisfactorily solved, given the uncertainties associated with the behaviour of these systems and their more complex characterisation. For these reasons, in conjunction with the ever present need and efforts to identify and mitigate systematic errors in all the space geodetic techniques (e.g. Arnold et al. 2015; Appleby et al. 2016; Capdeville et al. 2016; Tourain et al. 2016; Exertier et al. 2017; Ŝtêpánek and Filler 2018; Bergstrand et al. 2018; Dach et al. 2019; Luceri et al. 2019), it is therefore timely to revisit this topic in an attempt to improve the standards of modelling. 1.1 Ground-based testing of LAGEOS and LAGEOS-2 In principle, it would appear reasonable to assume that the problem of CoM determination could readily be solved by testing the laser retroreflector arrays prior to launch. In reality, the diverse and continuously evolving landscape of tracking systems deployed worldwide (Wilkinson et al. 2018) all but completely invalidates this proposition, as the experience gained with the LAGEOS and LAGEOS-2 satellites demonstrates. As part of these two missions, a great deal of effort was devoted to characterise the optical behaviour of the laser targets before launch. The outcome of these tests is detailed in two extensive reports by NASA (Fitzmaurice et al. 1977; Minott et al. 1993), which include the results of multiple experiments to determine the range correction and the lidar cross-section of the retroreflector arrays. The testing done for LAGEOS was somewhat limited technically relative to the standards that were to come in only a few years time, driven by developments in laser and detector technology. As an example, the accuracy requirements for LAGEOS mentioned in Fitzmaurice et al. (1977) are 2 cm, a long way from the current 1 mm long-term geodetic goal set by GGOS (Plag and Pearlman 2009). The tests carried out in 1977 investigated the effects of laser pulse spreading, range corrections for leading edge and peak detection methods, and the pulseto-pulse variation of the returned signal. The average range correction values obtained were 251 mm for leading edge half maximum, and 249 mm for peak detection, using laser 123 pulses estimated to be 60 ps long. These values would not quite agree with those obtained over a decade later during the pre-launch optical testing of LAGEOS-2 (Minott et al. 1993), where technical limitations and certain neglected sources of systematic errors in the previous testing were mentioned as the source of the disagreement. By 1993, the accuracy requirement for LAGEOS observations was set to 1 cm, still an order of magnitude worse than current demands. A battery of tests and methods were employed to investigate, first and foremost, the range corrections that apply to LAGEOS-2, a satellite designed and built to be as identical as possible to the original LAGEOS. But, for all the care and attention to detail that this report entailed, no definitive values with a practical, real-world use could be provided, as the authors of the report themselves acknowledged and forewarned. In essence, the different methods used to optically test LAGEOS-2 were measuring different things; some techniques worked purely in the optical domain (analysis of far-field diffraction pattern obtained with a continuous wave laser), and others attempted to replicate the measuring technology of a particular kind of SLR system present in the network (micro-channel plate photomultiplier, MCP-PMT, detection with constant fraction discrimination). Results obtained with different methods did not agree with each other, which, given the various limitations and caveats discussed at length in the report for each measuring technique, is unsurprising. For laser pulse lengths similar to those employed for the testing of LAGEOS, the newly determined range corrections were 1–2 mm longer (peak and leading edge detection), but depending on the particular method, the reported corrections ranged from 247 mm (streak camera measurements with 128 ps pulses, centroid) to 251 mm (MCP-PMT with constant fraction discrimination and 130 ps pulses). In conclusion, after painstaking efforts to measure directly the range corrections for LAGEOS, the outcome is the characterisation of the LAGEOS retroreflector array under a multitude of experimental conditions, which although invaluable for those interested in the practical technicalities of the SLR measuring process, it nevertheless does not provide the CoM quantities required for the analysis of laser ranging data. This is not caused by specific limitations of the testing performed for LAGEOS and LAGEOS-2, but simply because even if those were overcome, direct measurement of hardware-specific range corrections would necessitate as many different detection systems as there are deployed on the ground, which, ignoring the problem of systems that do not yet exist at the time of testing, is wholly unrealistic. 1.2 Previous modelling of centre of mass offsets The centre of mass values for spherical geodetic satellites employed currently by users of SLR data have been computed on the basis of semi-empirical modelling of the optical Upgraded modelling for the determination of centre of mass corrections of geodetic SLR… response of the retroreflector arrays (Otsubo and Appleby 2003; Otsubo et al. 2014a). The optical responses can in principle be computed on the basis of theoretical considerations and the physical characteristics of the array, or they can be measured on the ground with the adequate equipment. The solution arrived at in Otsubo and Appleby (2003) was to use laser ranging observations made at the so-called single-photon level, exploiting the fact that the distributions of laser returns obtained with this ranging policy are related to the optical responses of the satellites in a simple mathematical way. This enables the recovery of the optical responses of the arrays, which are subsequently employed to model the reflected pulses obtained with any laser ranging system, and therefore estimate the relevant centre of mass correction. Beyond the limited precision, by current standards, of the data that had been employed to derive the CoM offsets for the two LAGEOS, Etalon and Ajisai, an outstanding issue has been the difficulty of modelling SLR systems operating at the multi-photoelectron level. The simple relation between the incoming distribution of photons reflected off the retroreflector arrays, and the resulting measured distribution of photoelectrons, only holds for systems ranging at the single-photon level or in any case for single-stop systems, where data acquisition and processing take place in the digital domain. For legacy multi-photon systems, the additional delays and signal transformations in the detector and timing electronics, of an analogue nature, are not the same for terrestrial calibration and for satellite data collection, necessitating additional modelling to derive an appropriate CoM value. Previously, only optical delays (resulting from laser pulse spreading by the retroreflector arrays) were explicitly considered, along with a general estimate of the effects caused by additional delays in the detection systems. Given the limited information available on the timing and electronic behaviour of the various devices involved in the detection of retroreflected laser pulses, a system-specific computation taking into account these effects had not been attempted. 2 Computation of satellite optical response The starting point for the present work is identical to that in Otsubo and Appleby (2003) and Otsubo et al. (2014a), namely the computation of the theoretical optical response function of the various targets. The method followed is essentially that demonstrated in these previous reports, with several improvements mentioned in the following as relevant. The input data required here are the physical characteristics, dimensions, and locations of the individual cube-corner retroreflectors. Assumptions and simplifications include not considering the polarisation state of light (only its average), supposing perfect geometry for individual retroreflectors with no dihedral angle offsets, and incoherent behaviour of the retroreflected pulses. First, we compute a 2D reflectivity map of a cube-corner retroreflector by ray tracing, using Fresnel equations and modelling the path of the incoming light as a refraction upon entering the cube corner; reflections at each of the three internal surfaces; and a final refraction as the ray exits the front window of the retroreflector. Depending on whether the posterior sides of the cube corners are coated with a metallic layer or not, the internal reflections are modelled by total internal reflection (if the incident angle is smaller than the critical angle for the specific orientation being considered), or as reflections off metallic surfaces of the specified material (e.g. aluminium for the Etalon retroreflectors and silver for Starlette’s). Relevant details documenting this aspect of the modelling can be found in Arnold (1979). The effective area at every incident angle can be modelled with a simple expression for retroreflectors with a circular front window (Arnold 1979). Shapes other than circular require taking into account the loss of symmetry in the azimuth direction, which in general means that numerical methods are needed to compute the effective area for each incident angle. For the truncated triangular front windows of Ajisai, the same effective area maps computed for Otsubo and Appleby (2003) were used here. For the hexagonal entrances of Etalon, we found that the large number of retroreflectors (2146) makes the deviation from circularity negligible on average, so that no special treatment is required. The optical response function for a particular orientation is computed as the sum of each contributing reflection, taking into account each cube-corner retroreflector visibility, effective area, and reflectivity. Averaging overall possible orientations gives the optical response function for the entire array. The intensity of each reflection is modelled according to the relation Intensity ∝ arean ∗ reflectivity (1) where the exponent n is an adjustable parameter that takes into account, in an average sense, the complex heterogenous intensity distribution of the far-field diffraction pattern that results from nonzero dihedral angle offsets, imperfections in the construction of the cube corners and retroreflector array, thermal effects, and velocity aberration. For an ideal circular retroreflector and an observer in the on-axis direction, the farfield diffraction pattern is the Airy disc, and the value of the parameter n would be 2. Values for n that adequately describe the average intensity in real-world situations will always be lower than this. It must be noted that the absolute lack of knowledge about the precise as-constructed characteristics and state of each individual cube corner renders it impossible to model and predict the retroreflected intensity from a first principles perspective. Thus far, the result of the theoretical computation is a single-parameter family of functions; it is necessary to 123 J. Rodríguez et al. Table 1 Input data for empirical distributions. Full-rate single-photon observations (post-filtering) from Herstmonceux station Observations (M) Pass sections 15 1069 LAGEOS-2 9.5 700 LARES 10 314 Etalon 1.0 518 Starlette 4.5 159 Ajisai 5.3 117 LAGEOS determine the values of the exponent n that best describe the optical behaviour of each satellite, by comparison with empirical distributions collected at the single-photon level. This operation policy is observed by some stations of the ILRS (International Laser Ranging Service) network (Pearlman et al. 2002, 2019), whereby the signal level of the returned laser pulses is kept sufficiently low to ensure, statistically, that no more than a single photon triggers a photodetection in the receiver electronics. This is done to maximise operational consistency by avoiding intensity fluctuations, which can distort the shape of the distributions of detections as well as cause electronic biases in the detection equipment. An outcome of this mode of operation is that the probability distribution of detected photoelectrons is simply the convolution of the satellite optical response function and the temporal response, or noise profile, of the receiver. By modelling or measuring the system noise (e.g. with a terrestrial calibration to a flat target or single cubecorner retroreflector), the problem is reduced to comparing the family of functions arising from expression 1, convolved with the system noise, to the empirical distributions of detections. For this work, we employed the full histogram of a terrestrial calibration as the system noise. We obtained empirical distributions for each satellite using single-photon data from the Herstmonceux station. Since the end of 2014, this station operated almost exclusively with a high repetition laser (1 kHz) with very short pulses (8 ps), actively controlling the return rate to keep the detection level strictly at single photons (a policy followed since 1992). Night-only passes of a certain minimum duration, to avoid possible artefacts during residual flattening, were centred about a reference point in the leading edge of each pass distribution, processed with a Poisson filter to eliminate possible multi-photon events (Rodríguez et al. 2016b) and accumulated to obtain the required empirical distributions of detections. This process was manually supervised to detect and discard any blunders that could distort the resulting distributions. The input data used to build the empirical distributions are detailed in Table 1. In previous works, a limited number of satellite passes (sometimes as few as a single pass) were used to build the 123 empirical distributions on the assumption that high repetition lasers, with two orders of magnitude higher throughput than conventional ones, made up for the economy of input pass data. During this investigation, we found out that this assumption is not valid; there are orientation-specific signature effects that cannot be captured with the limited observation geometry afforded by a limited number of passes, requiring data accumulation over a longer time period to ensure that they are properly averaged and taken into account. Theoretical optical response functions were computed for n = [1, . . . , 2] at steps of 0.05 and convolved with a terrestrial calibration. The search for the parameter n was done by comparing these functions with the empirical distributions, after alignment and area normalisation. As far as possible, the range of the distributions employed included data far in excess of what is used by SLR stations/processing centres to compute normal points (usually an iterative mean with up to 3 σ rejection). This way we minimise variability in the selection of the best function, which is sensitive to the truncation level of the distributions under comparison. The selection of best-fit function was done by computing three different statistics: χ 2 , Kolmogorov–Smirnov, and the Bhattacharyya coefficient, common ways to compare similarity between pairs of distributions. These metrics agreed well with each other in most cases; when they did not, we used the median of the p parameter value indicated by the three methods. The empirical distributions and the best-fit functions, convolved with system noise, for all satellites are shown in Fig. 1. The quality of the fit is near-perfect for LAGEOS, excellent for LAGEOS-2, LARES, Starlette, and Ajisai, and only fair for Etalon. It appears that the current modelling for the Etalon satellites does not reproduce some features visible in the empirical distribution of detections, which show a slight slope change in the tail of the curve. Since we are not applying any weighting scheme, the best-fit function will be a compromise between fitting different parts of the empirical distribution, and the resulting function is that which best reproduces its overall shape. For the purposes of this work, this is satisfactory, as local deviations from the theoretically expected curve do not change greatly the statistics employed for data reduction and normal point formation, i.e. the mean and standard deviation of iteratively clipped distributions. For reference, the standard deviation of the 2.5 σ clipped empirical distribution for Etalon is 43.4 mm, and for the bestfit function 44.4 mm, with also a difference in their averages of 1 mm. This level of agreement exceeds by an order of magnitude the current estimated errors for this satellite for most stations of the network and will not be a practical limiting factor in the accuracy of the CoM values derived here. We note that for Ajisai the length of the distribution employed is truncated at approximately 85 cm from the theoretical centre of the satellite. This truncation is imposed by a feature of Upgraded modelling for the determination of centre of mass corrections of geodetic SLR… Table 2 Best-fit optical response functions (λ = 532 nm), from this work and Otsubo and Appleby (2003), Otsubo et al. (2014a) n parameter This work Previous results LAGEOS 1.25 1.1a LAGEOS-2 1.10 1.1a LARES 1.30 1.0 Etalon 1.40 1.3 Starlette 1.40 1.2 Ajisai 0.75 1.2 a LAGEOS Fig. 1 Optical response functions. Blue: empirical histograms; orange: best-fit functions; black: difference. Top to bottom: LAGEOS, LAGEOS-2, LARES, Starlette, Ajisai, and Etalon. Note the different horizontal scales, representing the distance from the satellite centre the data collected by Herstmonceux station, which presents a strong secondary reflection caused by the internal optics used to limit the return rate. This reflection does not affect in any way the regular data reduction procedure (equivalent to 2.5 σ clipping), as its location is outside the clipped region used for normal point formation. Etalon data do not show this feature because the signal intensity from these satellites satellites treated together in previous studies is very weak and need no limiting to keep detections at the single-photon level. The values for the parameter n corresponding to the bestfit optical response functions are detailed in Table 2. Several items are noteworthy here. First, in some cases the results do not agree with those obtained previously within what was thought to be the uncertainty of the method (approximately ± 0.1 in Otsubo and Appleby 2003). However, this level of uncertainty was probably too optimistic, and certainly that is the case for the results obtained with low-repetition data in 2003 for LAGEOS, Etalon, and Ajisai. The n values for Starlette (and the identical Stella) previously computed (Otsubo et al. 2014a) ranged from 1.1 to 1.3, and a higher upper limit would obtain if greater resolution for the test n values had been used. For LARES, previous computations resulted in n = 1, on the basis of fitting data from a more limited number of passes (Neubert and Otsubo 2012; Otsubo et al. 2014a). For these satellites, as well as for the LAGEOS pair, the limitations and uncertainties in the determination of the best-fit n parameter of the optical response functions can reasonably explain the disparities between previous results and the present ones, which are subject to less variability arising from greater accumulation of input data. On the other hand, the difference in Ajisai results is more significant. We note that the distributions employed to fit the data in previous work were clipped by the reduction algorithm and therefore had a usable range shorter by almost 5 cm compared to this study. The comparison between asymmetric distributions is sensitive to the length of the tail, with a tendency to overfit the peak region if not enough data after the peak of the distribution are included. Another factor is the limited resolution of the histograms afforded by low-repetition data, 6 mm in 2003 (12 Hz rate laser) compared to the current 0.25 mm, which would only coarsely realise the features of the underlying distributions. For testing purposes, we recomputed the best-fit functions for Ajisai using the old, already clipped, data and obtained results in agreement with those previously published (Otsubo and Appleby 2003), within realistic uncertainty. We conclude that with the several lim- 123 J. Rodríguez et al. imately in our normal point return rate estimations), whereas from their slightly different average altitudes the opposite would be expected (−18% approximately according to 1/R 4 space losses). 3 Modelling of detection systems Fig. 2 Difference in average range bias for LAGEOS and LAGEOS2. Error bars represent the 95% confidence intervals of the differences between means. Stations ordered, left to right, by productivity over the solution period (1993–2018). In red, the weighted average of the range bias difference, −1.2 mm (95% CI 0.3) itations and uncertainness mentioned above, some of which were much reduced for this work, the present result for Ajisai is the current best value to represent the optical behaviour of this satellite. For satellites LAGEOS and LAGEOS-2, we have computed the optical response function separately for each retroreflector array, instead of assuming identical behaviour as in previous works. The large number of observations that we accumulated and their high precision allowed us to notice small differences in the shape of the empirical distributions of detections, which imply slightly different optical response functions for the two objects. The consequence is a small difference in the computed CoM corrections of up to 0.8 mm, dependent on detection systems and reduction procedures. The engineering tolerances in the construction of LAGEOS2 clearly exceeded geodetic requirements, e.g. the acceptable range for the requested average equatorial radius of the satellite was 0.127 mm, and the actual measured values for both hemispheres differed from the central point of that range by 12.5 µm (Cogo 1988). Nevertheless, the long-term averages of range bias estimates of these two satellites indicate a small difference in their apparent size, which we estimated as −1.2 mm (95% CI 0.3) (see Fig. 2). This would imply that the apparent size of LAGEOS relative to LAGEOS-2, as seen by the SLR detection systems, is greater by this amount, a result already noted by Ries and Cheng (2014). Of course, this conclusion could only hold under the assumption that the optical responses of these two satellites are identical, which are not. Although the retroreflector arrays are geometrically identical, differences in their optical behaviour can explain the disparity. The returned laser pulses contain photons reflected off the visible section of the satellites (limited by the cubecorner acceptance angles); small differences in the average reflectivity profile of individual cubes will result in a different average optical response, with different location and spread (our findings). We note that we also know that LAGEOS and LAGEOS-2 are not optically identical from examination of the average return rates achieved across the network, where higher return rates are obtained for LAGEOS (+ 10% approx- 123 Having determined the optical response functions for the satellites of interest, we have to model the behaviour of the different kinds of detection systems present in the ILRS network. Previously, in the works upon which the current CoM values are based, three cases had been considered: single-photon ranging, compensated single-photon avalanche diodes (CSPAD), and leading edge photomultiplier (PMT) systems. The distinction between the first two cases does not correspond to any fundamental, qualitative difference between modes of operation, but rather to a concept—single-photon ranging—that differs from other operation policies only in a quantitative way, and that strictly speaking is never perfectly realised. Here, we have modelled all stations operating at single-photon or otherwise using single-stop detectors (SPADS, CSPADS) in an identical way; stations using PMT or micro-channel plate (MCP) detection devices were modelled separately. 3.1 Single-photon operation and single-stop detection Stations restricting the intensity of the laser return so that, at most, only about 10% of the laser pulses result in a successful detection, are said to operate in the single-photon mode. It is trivial to demonstrate, with Poisson statistics, that at this detection level the probability of obtaining multiple photoelectron events is very low (< 0.5%). This ensures that detections are not biased towards the leading edge of the returned pulses, a situation that obtains with increasing probability as return rates rise. The detection probability of individual events is then proportional to the probability distribution of the return pulses, mediated by the system noise profile, a situation we exploited to obtain the optical response functions as detailed in the previous section. The detection probability for single-stop detection systems not operating at single-photon rates can be easily computed for any average rate of detections (Neubert 1994). Of course, at low return rates the computed probability distributions converge to the single-photon case. Thus, we model all these systems identically, using for each satellite actual average return rate data, therefore taking into account deviations from the single-photon regime in a quantitatively correct way, avoiding guesswork about return intensity levels and not relying on the partial, qualitative, and at times inaccurate information relayed through the station logs available from the ILRS. In Upgraded modelling for the determination of centre of mass corrections of geodetic SLR… order to compute the CoM values, we reduced the computed distributions for each system with the relevant algorithms detailed in the site logs for each stations (e.g. iterative rejection at N × σ ). We note that we are ignoring detector time-walk effects in the laser range measurements that inevitably arise when operating at high-intensity levels of return. Many stations employ CSPAD detectors, which greatly compensate for time-walk (Kirchner and Koidl 1999); however, the compensation can never be perfect since the calibration of these devices is only strictly valid for light pulses of a specific shape; whereas, the pulse width of satellite reflections is ultimately given by the different sizes of the retroreflector arrays. 3.2 Multi-photon operation Stations employing PMT or MCP detectors, very common in the previous generation of SLR systems, normally do not decrease the signal intensity to low levels, but instead attempt to obtain as high return rates as possible (ideally without exceeding the linear behaviour of the detection electronics). From a modelling perspective, the problem that arises is that these are analogue devices, whose output signals are then discriminated and timed, processes we can no longer assume to be identical between terrestrial calibrations and satellite ranging and must be taken into account somehow. Here, we have attempted to model these systems on the basis of the relevant key hardware specifications available from the station logs. Monte Carlo simulations are a common tool in fields such as particle physics and medical imaging device development, aiding in the exploration and assessment of new experimental setups involving photodetector devices (e.g. Assié et al. 2004; Liu et al. 2009). We followed the same strategy here: for each station we performed a Monte Carlo simulation of the photodetection and pulse timing, for both terrestrial calibrations and satellite ranging, for the six targets under consideration. The differences between the computed average timing points of the calibration and satellite cases are the CoM corrections sought. We did not attempt to simulate the signal amplification process inside the detectors, or any other aspects for which no information is in any case available, such as gain variations across photodetector area, nonlinear behaviour, or the contribution to the system jitter of discriminator units and cabling. In essence, we simulated the temporal behaviour of ideal detection systems comprised of detector, discriminator, and timer. We computed input light pulses on the basis of Poisson statistics from the average return rates of normal point data, drawing random (Poisson) samples from the convolution of the laser pulses and the optical response functions; each sample, representing a single photoelectron event, was simulated as a Gaussian function whose width corresponds to the 10–90% rise time of the detector under Fig. 3 Simulation of multi-photon data example. Calibration and LAGEOS pass for Matera system (7941). Bottom plot shows an actual pass consideration, convolved in turn with the known transit time spread and precision of the timing device. The sum of the single photoelectron signals represents the signal output of the detector for each simulated pulse, as seen by the timer, which we discriminate with a straightforward leading edge method (0.25 fraction of the peak). Repeating this procedure, a sufficient number of times we obtain the average distributions of detections, which we processed with the relevant normal point reduction algorithm in use by each station. An example of the simulated distributions of detections for both calibration and satellite ranging is shown in Fig. 3. The dispersion of the simulated distributions is in good agreement with the average values reported in the normal point data (about 1.2 and 3.5 mm for calibration and LAGEOS, respectively). The simulated data for the satellite case are slightly skewed compared to the actual data from the satellite pass, but the differences in the averages of the distributions caused by this appear very small (sub-millimetric in the case of LAGEOS satellites). 4 Results 4.1 Model validation The CoM values are computed as the difference between the reference points of the modelled distributions of satellite and calibration data. Clearly, the accuracy of the resulting CoM 123 J. Rodríguez et al. Fig. 4 RMS of expected distributions plotted against RMS of normal point data for all known ILRS system configurations, shown by different ranging policies and detector types. Data split according to retroreflector array size. Orange: small targets (LAGEOS, LARES, and Starlette); blue: big targets (Ajisai and Etalon). Horizontal error bars represent the standard deviation of the empirical NP RMS. Insets: histograms of differences between RMS of modelled and NP data. The number of system configurations plotted in each case is noted values depends on the accuracy of the simulated distributions of detections (in the case of multi-photon stations), and on the accuracy of the distributions directly computed on the basis of available hardware specifications (in the case of single-photon/single-stop systems). Although a direct validation on the basis of the CoM values derived here is impossible to realise, the comparison between the RMS values of the expected distributions and those from the normal point data provides an indirect validation of the models employed. This comparison, for all the system configurations of the network, is shown in Fig. 4 for different operation policies and detector types. In each plot, points along the diagonal correspond to cases where perfect agreement between the modelled and empirical data distributions has been achieved. With respect to ranging policy, in general the level of agreement obtained is superior for single-photon systems, either because the modelling of this type of system is better, because of the inherent observational stability advantages of this operation mode, or because of both. In the case of multi-photon systems, the scatter appears greater than for single-photon ones, although the greater number of system configurations with this ranging policy increases the apparent dispersion of the data. We explored the differences between different kinds 123 of detectors (SPADs and PMT/MCPs) for multi-photon stations, since the operation policy is the same but the underlying models employed in this work are not. Model performance is similar for either kind of detectors; however, in general for stations using SPADS the agreement between simulated and empirical data is slightly worse for multi-photon operation relative to single photon. We note that the scatter plots of Fig. 4 mask the density of data points that lie close to the diagonal, which can be better seen from the histograms of the differences between modelled and empirical data (insets in the figure). These reveal that the agreement between modelled and real distribution RMS for the majority of systems is actually fairly good, and that the overall differences between ranging policies and detectors are rather small. For all system configurations, the RMS of modelled distributions is within 5 mm of the average of the empirical ones in 54% of cases, and 72% are within 10 mm. For the top 25 performers of the network since 1983, providing 75% of the NP data, 59% of the system configurations are modelled within 5 mm of the NP RMS, and 78% within 10 mm. In any case, assuming identical errors in the computation of system calibrations and satellite data, inaccuracies in the dispersion of the modelled distributions translate to much less than one-to-one to inaccuracies in the computed CoM values. For instance, a 5 mm RMS underestimation in the computed distributions of LAGEOS for the Herstmonceux single-photon station would result in 1 mm overestimation of the CoM, and 10 mm would translate to approximately a 1.6 mm overestimation. For multi-photon stations, the range of dispersion values for the simulated distributions is quite restricted by the input data (hardware, operational and post-processing parameters); errors of ∼ 10 mm in the resulting distributions possibly originate from significantly erroneous/incomplete system information. The previous figures do not take into account the dispersion in the quantity we are using to assess the performance of the model, the RMS of the NP observations (horizontal error bars in Fig. 4). This dispersion is in many instances very large due to instability of the distributions of detections or in some cases due to a small number of data points. For each system configuration (and satellite), the NP RMS was computed as the average RMS from all the NP data available. Plotting the difference between model predictions and empirical data against the standard deviation of the NP RMS reveals that the agreement is correlated with this quantity (Fig. 5). For system configurations whose NP data are consistent to within ∼ 5 mm, the models predict accurately the RMS of the NP distributions, with 71% simulated within 5 mm and 89% within 10 mm. Clearly, lack of consistency in the empirical data is responsible for a good share of the model performance at the poor end. Unsurprisingly, predicted distributions tend to agree worse with empirical data for the earlier systems (bottom plot of Fig. 5), reflecting a combination of less capa- Upgraded modelling for the determination of centre of mass corrections of geodetic SLR… Fig. 5 Differences between the RMS of modelled NP distributions and empirical data for all satellites and system configurations, plotted against the standard deviation of the NP RMS (top), and as a time series (bottom) ble hardware, poorer quality control procedures, and possibly a less accurate knowledge of system specifications. In summary, we can indirectly assess the performance of the models that we employed to model the distributions of detections, upon which the computation of CoM values relies. The comparison with NP data indicates that the simulations perform quite well for the majority of combinations of system configurations and satellites, both for single- and multi-photon ranging policies and detector types. Many of the poor predictions arise from actual inconsistencies in the NP data, which greatly limits the attainable accuracy. 4.2 CoM values Following the methods described in the previous section, we computed the CoM values for the six satellites under study for all the system configurations for which we could find information. For this, we relied on the ILRS site logs, current and historic, searched the SLR literature published on the proceedings of engineering and technical workshops and contacted quite a few individuals who are or have in the past been involved with the running of stations with some details missing. The number of stations and unique system configurations identified in this work are greater than those for which CoM values had been computed previously (83 system configurations in 66 stations in the old values, which are the current ILRS standard, for 193 system configurations in 77 stations in the new ones). Many of the differences between the newly identified individual system configurations are minor; they correspond mostly to documented hardware changes in any of the components that can affect the CoM computation. Fig. 6 Distribution of the CoM values currently employed and the newly computed ones, per satellite, for all system configurations (different number for the two datasets). Note the different horizontal scales Although many of these do not give rise to a significantly different CoM value, they are included for completeness and correctness, to avoid setting an arbitrary threshold below which slightly different systems would be treated as one, and to have a traceable correspondence between known site changes and the applied CoM values. A visual comparison of the new CoM values with the old ones is shown in Fig. 6, where we plotted the distributions of values for each satellite as normalised histograms. Although the distributions of the old and new values are not strictly comparable because of the different number of system configurations, it is apparent that the new figures are, on average, lower than the previous ones for all satellites. For LAGEOS and LAGEOS-2, there is a difference of about 4.5 mm in the average of the two sets of values. For the smaller targets, Starlette and LARES, the differences are smaller, at approximately 2 mm in the average. We note that the distributions of old CoM values for these two satellites are strongly bimodal, which makes the modest differences in their averages misleading, as they suggest a greater similarity between both sets than there actually is. This bimodality does not arise as markedly in the new CoM results, possibly because the new models depend on more input parameters, with the consequence that the output is more granular and smoothed out over the range of possible values. The greatest differences are seen for the biggest retroreflector arrays, i.e. satellites Ajisai and Etalon, for which the difference in the averages of the distributions is over 20 mm. A significant factor driving the CoM results towards lower values is the use in this work of actual signal return rates, computed from NP data. Photoelectron signal strength was probably overestimated previously, even for low-Earth orbit satellites, despite their much higher 123 J. Rodríguez et al. link budgets than high-altitude objects. Here, we found that many stations following a multi-photon ranging policy actually only achieve a modest return rate, which when taken into consideration effectively reduces the computed CoM values. Equally relevant, in addition to modelling the optical pulse spreading caused by the finite depth of the retroreflector arrays, as done by Otsubo and Appleby (2003); Otsubo et al. (2014a), we considered the electronic signal pulse spreading caused by the detection hardware, which also reduces the location of the reference points in the range distributions. We note that previously published studies did not provide station-by-station corrections, but ranges of values for different ranging technologies and operational circumstances used only later on to derive CoM offsets for individual system configurations. Here, we have taken into account from the onset the specific characteristics of each station, instead of exploring the CoM corrections resulting from the whole gamut of possible combinations of system hardware, satellite targets, and ranging policies. A strict comparison between the values arising from the present and previous models is therefore not readily available, as the differences emerging from comparing individual station values may have multiple causes (e.g. system configuration details and actual operational conditions). Nevertheless, an inspection of the corrections computed for LAGEOS and Etalon reveals that many of the greatest differences are found for high-energy multi-photon stations, for which a new modelling strategy has been employed here. Noticeably, different CoM values are also obtained for a number of stations employing SPAD detectors with unclear ranging policy, as stated in their site log files (e.g. single-photon to few or multi-photon); this could in some cases indicate model improvement, or perhaps more accurate corrections thanks to the use of return rate data. 5 Range bias analysis with new CoM values The newly computed CoM values present in many cases significant differences at the level of several millimeters. Changes in the CoM offsets cause a 1:1 change in the reduced measured distances to the satellites, which in turn will impact the estimation of coordinates, albeit with a different scale factor. In practice, SLR is most sensitive to the UP component of station coordinates, so CoM changes affect mainly the estimated height, mapping onto it with a slightly higher than 1:1 ratio, depending on observation geometry. In the previous section, we have shown that the general tendency, for all satellites, is for the new CoM values to be lower than those from the previous CoM tables. The implication is that applying the new CoM offsets to correct the observations will, in general, result in greater station heights, i.e. a greater SLR network scale. However, it is now recognised in current modelling of SLR data that it is advantageous to estimate range 123 bias parameters to absorb any possible observational errors and modelling deficiencies Appleby et al. (2016); Luceri et al. (2019). In this analysis context, employing a different set of CoM offsets will therefore impact primarily the estimated range errors for each station, in principle without affecting station coordinates, which should already be free of biases. The behaviour of such biases, estimated with the use of the new CoM tables, and their comparison with those previously obtained, is discussed next. We have computed orbital solutions for all satellites estimating station coordinates, and range biases, applying either the old CoM tables or the new ones presented here. For satellites LAGEOS-1/2 and Etalon-1/2, we used the software SATAN (Sinclair and Appleby 1988) to set up solutions in an analogous fashion to the products delivered by the official ILRS Analysis Centre NSGF: computing 7-day arcs solving for initial conditions, daily EOPs, weekly site coordinates and range biases in combined 4-satellite solutions (2000– 2018). For the low-Earth orbiters, we employed the c5++ code (Otsubo et al. 2016) to compute yearly 6-satellite solutions (LAGEOS-1/2, Ajisai, Starlette, Stella and LARES), solving for initial conditions (LAGEOS-1/2 every 5 days, LEOs every 3 days), 5 × 5 gravity field, station coordinates and range biases estimated annually for the period 2015–2017. In these solutions, the RB for the satellite pairs LAGEOS-1/2 and Starlette/Stella were combined. The RB for LEO satellites obtained with the application of the new CoM tables were obtained indirectly from the estimates by adding the difference between the applied CoM values and the new ones; in the case of LAGEOS and Etalon, this was unnecessary because we had the two sets of solutions readily available. The differences between the absolute values of the RB estimates averaged over the period of analysis are shown in Fig. 7 for the most productive ILRS stations. Range biases obtained with the old CoM were subtracted from those obtained with the new CoM values, i.e. negative bars in the figure represent stations for which the absolute estimated biases are smaller, on average, when employing the new tables. Conversely, positive bars indicate, where present, the magnitude by which the absolute value of the estimated biases is increased by applying the new CoM offsets. The results for the LAGEOS pair are very similar, a consequence of the small difference in the values of the newly computed CoM figures for these satellites. The estimated RB increases for a few stations, although for the majority (about two thirds) it does not change or is smaller with the new CoM values. Etalon and Ajisai, the biggest targets, show a similar behaviour where many stations (all but four for Etalon) see their estimated biases reduced notably by up to over 15 mm and 20 mm, respectively. The validity of the CoM values for Etalon and Ajisai has been questioned in recent years from analysis of averaged range bias time series, which showed the presence of considerable positive biases Upgraded modelling for the determination of centre of mass corrections of geodetic SLR… Fig. 8 Weekly scale parameters between combined LAGEOS-1/2 and Etalon-1/2 solutions computed with new and old CoM values, with and without estimating errors, plus linear fits the new values has a negligible effect on the averaged estimated biases (especially given the big uncertainties of some of the estimates), and a minority for which differences of up to 2.5 and 5 mm, respectively, are seen. These differences do not appear to be in the same direction, i.e. the overall change is approximately neutral in terms of increasing or reducing the biases present in the network. In summary, the application of the new CoM values results in very significant changes in the overall biases estimated for satellites Etalon and Ajisai, less notable but consequential changes for the LAGEOS pair, and a more modest impact on the errors recovered for LARES and Starlette. 5.1 Geometric frame scale Fig. 7 Difference of absolute values of mean range biases per station, estimated with the new and old CoM models. Negative values indicate a reduction in estimated RB when using the new CoM tables. Error bars represent the 95% CI of the difference between means. Note the different vertical scales, and the use of logarithmic scale outside the region delimited by dashed lines for LARES, Starlette, and Ajisai for many stations Dunn and Otsubo (2014); Rodríguez et al. (2016a). The magnitude and sign of these biases, together with the lack of correspondingly large values for the RB estimates for other satellites, made the CoM offsets employed suspect. The effect of applying the new values presented here is to reduce to a large extent the range errors estimated for these satellites. In the cases of LARES and Starlette, there are two groups of stations: those for which the application of In principle, any effects arising from CoM inaccuracies will have already been taken into account in orbital solutions where range errors are estimated (such as those currently under development by the ILRS analysis centres). But given the magnitude and, especially, the general tendency towards lower values in the newly computed CoM offsets, it is relevant to examine here the impact on station coordinates that follows from their application alone. The effect on station heights is best seen at the network level, by comparing the frame scales obtained from solutions where the two different sets of CoM values have been applied, with no RB estimation. The time series of the relative scale parameters for these two solutions are shown in Fig. 8, with linear fits over the computation period. The mean scale difference at the approximate mid-point, 2009, is 0.65 ppb (95% CI 0.08), with the greater scale corresponding to the solutions where the new CoM offsets were applied. Comparing the scale parameters between solutions with range errors estimated shows no difference, as expected. Notably, the scale change caused by using the new CoM tables is of a very similar magnitude (and same direction) to what is obtained if range errors are estimated for all stations of the network regardless of CoM model: 0.72 ppb (SE 0.13) in Appleby et al. (2016). In principle, this would suggest that most of the measurement errors recovered in their estimation are simply model errors in the computation 123 J. Rodríguez et al. Fig. 9 Average RB estimated per station and satellite, obtained with the new CoM model. Error bars represent the 95% CI. Note that for the y axis the scale is logarithmic outside the interval [−15, 15] mm, indicated by the dotted line of CoM corrections. We examine next whether this is true at the individual station level; i.e. to what extent the application of the new CoM offsets and the associated removal of previously observed range errors leads to unbiased observations. 5.2 Remaining biases Shown if Fig. 9 are the averages of the absolute values of estimated range biases for the six satellites under study over 123 the relevant analysis periods. For LAGEOS and LAGEOS-2, the averaged estimated range errors remaining after correcting the observations with the new CoM values are below 5 mm for about 50% of the stations, no greater than 10 mm for ∼ 30% of the stations, and greater than 15 mm for 15% of the stations. Large errors of over 10 mm are unlikely to be caused by CoM mismodelling for these targets, indicating problems with the data of another nature. Although using the new tables the estimated errors are reduced considerably for Etalon and Ajisai (Sect. 5), there remain for many stations unexplained range biases of over 1 cm. For these satellites, we cannot discard that an important contribution towards the estimated errors may be originated by imperfect CoM modelling, although it is reassuring that the distribution of averaged errors is centred around zero. In the case of LARES and Starlette, for stations providing sufficient observations, the remaining estimated biases are of a similar magnitude to those obtained for LAGEOS. Stations with sparse tracking coverage of low-Earth orbiters, the solutions are consequently poor. Still, it is interesting to note that there is a high level of agreement between the averaged errors of the three LEO satellites for stations at the low-productivity end (in general with big confidence intervals in the figure), which may indicate that their root cause is a genuine problem at those stations. For other, more productive, stations there is substantial agreement between the range biases obtained for LARES and Starlette, but not so much for Ajisai. We must note that for LAGEOS-1/2 and Etalon-1/2 these are values averaged over a long period and therefore do not represent current network performance. However, they are useful as a general assessment of the realistic accuracy achieved at the level of individual stations throughout the latest 18 years, which is approximately half of the whole period of SLR history currently employed for global reference frame realisation. A notable outcome of employing the new CoM offsets is that the remaining average biases are more evenly distributed around zero, i.e. the number of stations with positive errors is similar to that with negative ones. This is true in particular for satellites Etalon and Ajisai, for which very large positive biases across the network were obtained with the old CoM corrections. Although to a lesser extent, this also applies to LAGEOS and LAGEOS-2, where predominantly positive biases were also present, driving the scale of the SLR network down towards lower station heights (Appleby et al. 2016). The picture that arises from this analysis is that despite the neat coincidence between the global scale differences obtained with the CoM values presented here, and those resulting from explicit error estimation, the situation for individual stations is more complex. While the global network scale parameter for the two solution types is similar, the individual station heights that determine it are different wherever observational errors are still present. To illustrate Upgraded modelling for the determination of centre of mass corrections of geodetic SLR… Fig. 10 Height difference relative to ITRF2014 of four selected stations, estimated with and without range biases, using the new CoM values this, in Fig. 10, we show the heights relative to ITRF2014 for four selected stations, estimated with and without range biases, using the new CoM values. The first two stations, 7840 (Herstmonceux, UK) and 7841 (Potsdam, Germany), are examples where the range biases remaining after using the new CoM values to correct the observations is small (less than 4 mm long-term average for LAGEOS-1/2 and Etalon-1/2). Thus, the heights obtained either estimating biases or simply using the new CoM values are quite similar, with the benefit of reduced noise levels in the solutions with no RB estimation. The next two stations, 7110 (Monument Peak, US) and 7832 (Ryadh, Saudi Arabia), present biases of increasingly greater magnitude (7110: negative ∼ 8 and ∼ 14 mm for LAGEOS-1/2 and Etalon-1/2, respectively; 7832: positive ∼ 21 and ∼ 17 mm). As a result, systematic differences in their estimated heights appear (modest for 7110, severe for 7832): the new CoM offsets do not and were not expected to free the observations of other possible error sources. The conclusion is that notwithstanding the remarkable agreement between the network scale offsets discussed in Sect. 5.1, the application of the new CoM values alone, without taking into account other range errors, may lead to inaccurate station heights and therefore potentially incorrect vertical site velocities as well. 6 Dynamic scale The most accurate determination of the dynamic scale of the Earth, the terrestrial gravitational constant GM, is pro- vided by laser ranging. The adopted IERS standard is the estimate by Ries et al. (1992), obtained from analysis of five years of LAGEOS observations (1986–1991), with a value of 398,600.4415 ± 0.0008 km3 s−2 , where the time unit is TT second. This estimate was larger than previous ones primarily because the centre of mass correction employed was larger than what had hitherto been applied (240 mm), by a remarkable 11 mm difference. The origin of this difference was a revision of the range corrections derived for LAGEOS prompted by the pre-launch optical testing of LAGEOS-2, where the previous analyses appear not to have used the larger figures derived from the optical testing of LAGEOS. A few other estimates of GM were computed subsequently using more recent range data from satellites LAGEOS and LAGEOS-2 (and in some cases LEO satellites), resulting in general in somewhat higher GM values, e.g. 398,600.4419 km3 s−2 (Dunn et al. 1999); 398,600.44187 km3 s−2 (Smith et al. 2000); 398,600.4417 km3 s−2 (Pavlis 2002); 398,600.44163 km3 s−2 (Ries 2007); 398,600.44184 km3 s−2 (Otsubo et al. 2014b); 398,600.4416 km3 s−2 (Ries and Cheng 2014). As all these values fall within the uncertainty of the current IERS standard, the motivation to revise it has been limited. More recently, from analysis of the range errors computed from observations to the LAGEOS satellites for all stations of the ILRS network, we presented an indirect estimation of GM on the assumption that all range errors could be attributed to an inaccurate apriori value for this constant (Appleby et al. 2016). The conclusion was that, in light of the range errors identified across the network, the GM value resulting from the observations made during the recent 2008–2014 period would be higher (∼ 398,600.4419 km3 s−2 ) than the current standard. However, when only three of the most stable stations of the network, believed to be largely free of bias (Herstmonceux, Matera, and Yarragadee), were used for the computation, the average result (∼ 398,600.44155 km3 s−2 ) was in line with the value from Ries et al. (1992). The extensive changes in the computed centre of mass corrections presented in this work will have consequences for the estimation of GM from observations of the satellites treated here. Although the average CoM difference for LAGEOS relative to the old values is less than 50% of that which prompted the revised GM value in 1992, it still is clearly very significant. For the purpose of testing the impact of the new CoM corrections on the estimated GM constant, we setup combined 14-day arcs LAGEOS-1/2 and Etalon-1/2 solutions (2000–2018) where, except for range biases, which are highly correlated with GM, we solved for the same set of parameters detailed in Sect. 5, plus a global GM parameter for each arc. We computed two identical solutions, correcting the range observations with both sets of CoM values. The average of the estimated GM values obtained with the use of the old CoM corrections was 398,600.44200 km3 s−2 (95% 123 J. Rodríguez et al. have been the result of averaging over a very specific period. Be it as it may, our current estimation, consistent with the combined 4-satellite solutions used for terrestrial reference system realisation, computed over an 18-year-long period of time, and performed with the most up to date models available, agrees with the standard GM value within the much reduced statistical uncertainty of the new estimate. 7 Conclusions Fig. 11 Satellite laser ranging determinations of GM from various authors. The horizontal extent of the traces indicates the period of each analysis. Error bars represent different measures of uncertainty/dispersion depending on the study (SE for Ries92, Dunn99, Ries07 and Ries14; SEM for Pavlis02, and Appleby16; SD for Otsubo14; 95% CI for Rodriguez19). NB more than one solution is available from some of the studies CI 0.00005), in agreement with our previous considerations based on an indirect estimation. The geometric frame scale of this solution was essentially identical to that obtained with range bias estimation, indicating that the GM parameters absorbed errors in the (CoM corrected) range observations, affecting station coordinates in a similar way, on average. The result obtained when the new CoM values are employed is 3,986,400.44146 km3 s−2 (95% CI 0.00005), unchanged from the current standard. These results are plotted in Fig. 11, along with several previous determinations of GM for reference. These results indicate that previous GM estimations suffered from inaccuracies in the available CoM corrections, which were on average slightly larger than those presented here. As a consequence, a higher GM value was required to accommodate the longer ranges, relative to the situation we obtain with the new corrections. This does not seem to explain the value of the current standard (or other results close to it), obtained applying a single centre of mass value of 251 mm for all stations. Results presented in Dunn et al. (1999) indicate that the stability of the annual GM determinations on the basis of LAGEOS before 1990 was poor, with a clear tendency towards lower values compared to the later estimates (Fig. 2 in the reference). It is not known what may have caused this upwards trend and later stabilisation of GM, but a similar behaviour may be assumed in the solutions that provided the standard value in use today, which would 123 The realisation of a global reference frame that satisfies the ever more stringent needs of the scientific and engineering community of users relies currently on the relative strengths of four space geodetic techniques. Of those, and although research on the suitability of other techniques to make a contribution is ongoing (Haines et al. 2015; Männel and Rothacher 2017; Couhert et al. 2018; Kuang et al. 2019), currently satellite laser ranging remains the standard to realise the origin of the frame. In addition, the coordinates of the SLR network of stations define its scale, together with those of VLBI sites. A long-standing disagreement in the network scales of these two techniques, as determined in the official ITRF realisation, has been the driver of several investigations, focused on the need to uncover systematic errors that could compromise their accuracy (Altamimi et al. 2016). For SLR, it is possible to co-estimate errors in the observable along with the parameters of actual interest, at the cost of increased solution noise (Appleby et al. 2016). Unfortunately, errors thusly estimated offer no clues as to their origin, but their successful quantification is valuable to, first and foremost, obtain bias-free solutions, and secondly to direct the attention of the analysis and engineering SLR communities towards outstanding problems in the technique. This kind of analysis prompted us to question the accuracy of the centre of mass corrections for the Etalon satellites, and by extension those for the rest of SLR targets whose CoM values relied on the same models and assumptions. In this paper, we have presented refined models to compute centre of mass offsets for SLR observations to spherical geodetic satellites, applying them to derive values for the main targets tracked by the ILRS, for all stations of the network. The new values represent the most detailed attempt ever to compute these corrections, of critical importance for the accuracy of the technique. A direct validation of the new CoM values is not possible, but we have presented a metric that provides an indication of the goodness of the models employed in the computation. The comparison of the dispersion of the distributions of detections obtained empirically with that arising from the models shows an acceptable level of agreement in the majority of cases. Severe disparities between empirical and simulated data are nonetheless present, the causes of which include inconsistent observa- Upgraded modelling for the determination of centre of mass corrections of geodetic SLR… tions, unreliable or unknown system information, and of course, model deficiencies. The newly derived values are similar to the previously computed ones for the smaller targets examined here, with differences increasing for the bigger retroreflector arrays. Notably, the newly computed CoM offsets result in a significant reduction in the estimated range biases for satellites Etalon-1/2 and Ajisai, which should prove advantageous for the consistency and quality of multi-satellite combinations for the computation of terrestrial reference frame solutions and other geodetic parameters of interest (Sośnica et al. 2014; Bloßfeld et al. 2018). Most relevant for reference frame purposes, our CoM values also show significant differences across the network for the LAGEOS satellites, with consequent changes in the estimated heights for many sites. The net result is a change in the scale of the SLR network of similar magnitude to what is obtained when range errors are estimated. However, as those estimated biases include other (and in principle all) sources of errors, aside from CoM mismodelling, the application alone of the new offsets does not result for many sites in the same station heights. A contribution that our new CoM values offer resides in advancing our knowledge of the SLR error budget structure. Although the models cannot be considered perfect (we have demonstrated that they are not) and the resulting station coordinates are still contaminated by the presence of errors (if these are not accounted for), the fact that we obtain a difference in the network scale simply by applying the new values is highly significant. In order to achieve this, the remaining errors estimated for many stations have changed in magnitude and/or been subject to a sign reversal, which, accepting the superiority of the new CoM models over the previous ones, leave us with a better picture of the lingering issues that may be affecting the performance of the network. In addition to the problem of model improvement for the determination of station coordinates, we have probed the impact of the new corrections on the determination of the terrestrial gravitational constant, as CoM inaccuracies are considered the main limiting factor influencing the value of the dynamic scale. Our findings, from a new GM estimation based on solutions consistent with the official ILRS product, are in very close agreement with the currently employed standard value. Acknowledgements We wish to acknowledge the ILRS for providing the data, and in particular the stations of the network for their continuous tracking efforts. We also acknowledge many current and former ILRS colleagues who provided valuable information about tracking stations. Data Availability Statement The results obtained in this work will be made available to download from the ILRS website in the shape of centre of mass tables for the satellites discussed, as well as example software to interrogate them and retrieve appropriate corrections for the requested epochs and station. References Altamimi Z, Rebischung P, Métivier L, Collilieux X (2016) ITRF2014: a new release of the internationall terrestrial reference frame modeling nonlinear station motions. J Geophys Res Solid Earth 121(8):6109–6131 2016JB013098 Appleby G, Rodríguez J, Altamimi Z (2016) Assessment of the accuracy of global geodetic satellite laser ranging observations and estimated impact on ITRF scale: estimation of systematic errors in LAGEOS observations 1993–2014. J Geod 90(12):1371–1388. https://doi.org/10.1007/s00190-016-0929-2 Appleby GM, (1992) Satellite signatures in SLR observations. In: 8th international workshop on laser ranging instrumentation. Annapolis, USA Arnold D, Meindl M, Beutler G, Dach R, Schaer S, Lutz S, Prange L, Sośnica K, Mervart L, Jäggi A (2015) CODE’s new solar radiation pressure model for GNSS orbit determination. J Geod 89(8):775– 791. https://doi.org/10.1007/s00190-015-0814-4 Arnold DA, (1979) Method of calculating retroreflector-array transfer functions, Technical Report 382, Smithsonian Astrophysical Observatory Assié K, Breton V, Buvat I, Comtat C, Jan S, Krieguer M, Lazaro D, Morel C, Rey M, Santin G, Simon L, Staelens S, Strul D, Vieira JM, Van de Walle R (2004) Monte Carlo simulation in PET and SPECT instrumentation using GATE. Nucl Instrum Methods Phys Res Sect A Accel Spectrom Detect Assoc Equip 527(1):180–189. https://doi.org/10.1016/j.nima.2004.03.117 Bergstrand S, Herbertsson M, Rieck C, Spetz J, Svantesson CG, Haas R (2018) A gravitational telescope deformation model for geodetic VLBI. J Geod. https://doi.org/10.1007/s00190-018-1188-1 Bloßfeld M, Rudenko S, Kehm A, Panafidina N, Müller H, Angermann D, Hugentobler U, Seitz M (2018) Consistent estimation of geodetic parameters from SLR satellite constellation measurements. J Geod 92(9):1003–1021. https://doi.org/10.1007/s00190018-1166-7 Capdeville H, Ŝtêpánek P, Hecker L, Lemoine J-M (2016) Update of the correctivemodel for Jason-1 DORIS data in relation to the South Atlantic anomaly and a corrective model for SPOT-5. Adv Space Res 58(12):2628–2650 Cogo F (1988) Weight discrepancy analysis between LAGEOS-1 and LAGEOS-2 satellites, Technical Report LG-TN-AI-035, Aeritalia Space System Group Couhert A, Mercier F, Moyard J, Biancale R (2018) Systematic error mitigation in DORIS-derived geocenter motion. Geophys Res Solid Earth 123(10):142–10161. https://doi.org/10.1029/ 2018JB015453 Dach R, Suŝnik A, Grahsl A, Villiger A, Schaer S, Arnold D, Prange L, Jäggi A (2019) Improving GLONASS orbit quality by reestimating satellite antenna offsets. Adv Space Res. https://doi. org/10.1016/j.asr.2019.02.031 Dunn P, Otsubo T, (2014) Etalon and Ajisai observations from NASA’s SLR network. In: 19th international workshop on laser ranging. Annapolis, US Dunn P, Torrence M, Kolenkiewicz R, Smith D (1999) Earth scale defined by modern satellite ranging observations. Geophys Res Lett 26(10):1489–1492 Exertier P, Belli A, Lemoine JM (2017) Time biases in laser ranging observations: a concerning issue of space geodesy. Adv Space Res 60(5):948–968. https://doi.org/10.1016/j.asr.2017.05.016 Fitzmaurice MW, Minott PO, Abshire JB, Rowe HE (1977) Prelaunch testing of the laser geodynamics satellite (LAGEOS), Technical Report 1062, NASA Technical Paper Haines BJ, Bar-Sever YE, Bertiger WI, Desai SD, Harvey N, Sibois AE, Weiss JP (2015) Realizing a terrestrial reference frame using the 123 J. Rodríguez et al. global positioning system. J Geophys Res Solid Earth 120:5911– 5939. https://doi.org/10.1002/2015JB012225 Kirchner G, Koidl F (1999) Compensation of SPAD time-walk effects. J Opt A Pure Appl Opt 1:163–167 Kuang Da, Bertiger Willy, Desai Shailen D, Haines Bruce J, Yuan Dah-Ning (2019) Observed geocenter motion from precise orbit determination of grace satellites using gps tracking and accelerometer data. J Geod. https://doi.org/10.1007/s00190-019-01283-5 Liu S, Li H, Zhang Y, Ramirez RA, Baghaei H, An S, Wang C, liu J, Wong. W (2009) Monte Carlo simulation study on the time resolution of a PMT-quadrant-sharing LSO detector block for timeof-flight PET. IEEE Trans Nucl Sci 56(5):2614–2620 Luceri V, Pirri M, Rodríguez J, Appleby G, Pavlis EC, Müller H (2019) Systematic errors in SLR data and their impact on the ILRS products. J Geod (in press) Männel B, Rothacher M (2017) Geocenter variations derived from a combined processing of LEO- and ground-based GPS observations. J Geod 91:933–944. https://doi.org/10.1007/s00190-0170997-y Minott PO, Zagwodzki T, Selden M (1993) Prelaunch optical characterization of the laser geodynamics satellite (LAGEOS 2), Technical Report TP-3400, NASA Technical Paper Neubert R (1994) An analytical model of satellite signature effects. In: Proceedings of the 9th international workshop on laser ranging instrumentation. Canberra, Australia Neubert R, Otsubo T, (2012) The centre of mass correction for LARES and LAGEOS for single photon detection. In: 12th international technical laser workshop. Frascati, Italy Otsubo T, Appleby GM (2003) System-dependent centre of mass corrections for spherical geodetic satellites. J Geophys Res Solid Earth. 108(B4):2201. https://doi.org/10.1029/2002JB002209 Otsubo T, Sherwood R, Appleby G, Neubert R (2014a) Center-of-mass corrections for sub-cm-precision laser-ranging targets: Starlette, Stella and LARES. J Geod 89(4):303–312 Otsubo T, Matsuo K, Appleby G, Sherwood R, Neubert R, (2014b). Scale parameters of the earth sensitive to the optical response of spherical SLR targets. In: AOGS 11th annual meeting. Sapporo, Japan Otsubo T, Matsuo K, Aoyama Y, Yamamoto K, Hobiger T, Kubo-oka T, Sekido M (2016) Effective expansion of satellite laser ranging network to improve global geodetic parameters. Earth Planets Space 68(1):65. https://doi.org/10.1186/s40623-016-0447-8 Pavlis EC (2002) Dynamical determination of origin and scale in the earth system from satellite laser ranging. In: Ádám J, Schwarz KP (eds) Vistas for geodesy in the New Millennium, Vol. 125 of International Association of Geodesy Symposia. Springer, Berlin, pp 36–41 Pearlman MR, Degnan JJ, Bosworth JM (2002) The international laser ranging service. Adv Space Res 30(2):135–143 Pearlman MR, Noll CE, Pavlis EC, Lemoine FG, Combrink L, Degnan JJ, Kirchner G, Schreiber U (2019) The ILRS: approaching 20 years and planning for the future. J Geod. https://doi.org/10.1007/ s00190-019-01241-1 123 Plag H-P, Pearlman M (eds) (2009) Global geodetic observing system: meeting the requirements of a global society on a changing planet in 2020. Springer, Berlin. https://doi.org/10.1007/978-3642-02687 Ries JC, (2007) Satellite laser ranging and the terrestrial reference frame: principal sources of uncertainty in the determination of the scale. In: Geophysical research abstracts, Vol. 9, 10809. Vienna, Austria Ries JC, Eanes RJ, Shum CK, Watkins MM (1992) Progress in the determination of the gravitational coefficient of the Earth. Geophys Res Lett 19(6):529–531 Ries JC, Cheng MK (2014) Satellite laser ranging applications for gravity field determination. In: 19th international workshop on laser ranging. Annapolis, USA. https://cddis.nasa.gov/lw19/docs/ 2014/Papers/3117_Ries_paper.pdf Rodríguez J, Appleby G, Otsubo T, Dunn P (2016a) Accuracy assessment of CoM corrections for Etalon geodetic satellites. In: 20th international workshop on laser ranging. Potsdam, Germany. https://cddis.nasa.gov/lw20/docs/2016/papers/ 15-Rodriguez-paper.pdf Rodríguez J, Appleby G, Otsubo T, Sherwood R, Wilkinson M (2016b) Assessing and enforcing single-photon returns: Poisson filtering. In: 20th international workshop on laser ranging. Potsdam, Germany. https://cddis.nasa.gov/lw20/docs/2016/ papers/45-Poisson_paper.pdf Sarti P, Abbondanza C, Petrov L, Negusini M (2011) Height bias and scale effect induced by antenna gravitational deformations in geodetic VLBI data analysis. J Geod 85(1):1–8. https://doi.org/ 10.1007/s00190-010-0410-6 Sinclair AT, Appleby GM (1988) SATAN: programs for determination and analysis of satellite orbits from SLR data, SLR Technical Note 8. Royal Greenwich Observatory, Cambridge Smith DE, Kolenkiewicz R, Dunn PJ, Torrence MH (2000) Earth scale below a part per billion from satellite laser ranging. In: Schwarz KP (ed) Geodesy beyond 2000, Vol. 121 of International Association of Geodesy Symposia, 3–12. Springer. ISBN 978-3-642-64105-3 Sośnica K, Jäggi A, Thaller D, Beutler G, Dach R (2014) Contribution of Starlette, Stella, and AJISAI to the SLR-derived global reference frame. J Geod 88(8):789–804. https://doi.org/10.1007/ s00190-014-0722-z Ŝtêpánek P, Filler V (2018) Cause of scale inconsistencies in DORIS time series. Studia Geophysica et Geodaetica 62(4):562–585. https://doi.org/10.1007/s11200-018-0406-x Tourain C, Moreaux G, Auriol A, Saunier J (2016) DORIS Starec ground antenna characterization and impact on positioning. Adv Space Res 58(12):2707–2716. https://doi.org/10.1016/j.asr.2016. 05.013 Scientific Applications of DORIS in Space Geodesy Wilkinson M, Schreiber U, Procházka I, Moore C, Degnan J, Kirchner G, Zhongping Z, Dunn P, Shargorodskiy V, Sadovnikov M, Courde C, Kunimori H (2018) The next generation of satellite laser ranging systems. J Geod. https://doi.org/10.1007/s00190-018-1196-1