Measurement of the solar system acceleration using the Earth scale factor

We propose an alternative method to detect the secular aberration drift induced by the solar system acceleration due to the attraction to the Galaxy centre. This method is free of the individual radio source proper motion caused by intrinsic structure variation. We developed a procedure to estimate the scale factor directly from very long baseline interferometry (VLBI) data analysis in a source-wise mode within a global solution. The scale factor is estimated for each reference radio source individually as a function of astrometric coordinates (right ascension and declination). This approach splits the systematic dipole effect and uncorrelated motions on the level of observational parameters. We processed VLBI observations from 1979.7 to 2016.5 to obtain the scale factor estimates for more than 4,000 reference radio sources. We show that the estimates highlight a dipole systematics aligned with the direction to the centre of the Galaxy. With this method we obtained a Galactocentric acceleration vector with an amplitude of 5.2 $\pm$ 0.2 \mu as/yr and direction $\alpha_G = 281\deg \pm 3\deg$ and $\delta_G = -35\deg \pm 3\deg$.


Introduction
The solar system is rotating around the Galactic centre with a period of approximately 200 million years. The solar system barycentre shows, therefore, rotational acceleration in the direction of the Galactic centre of 5-6 microarcseconds of arc per year (µas/yr) and this effect, the secular aberration drift (SAD), results in a systematic proper motion of reference radio sources around the sky.
The theoretical prediction of the dipole systematics (Fanselow (1983); Bastian (1995); Eubanks et al. (1995); Gwinn et al. (1997); Sovers et al. (1998); Mignard (2002); Kovalevsky (2003); Kopeikin & Makarov (2006)) was recently confirmed after the analysis of a global set of very long baseline interferometry (VLBI) data since 1979 ; Xu et al. (2012Xu et al. ( , 2013; Titov & Lambert (2013; MacMillan (2014)). The magnitude of the dipole systematics is consistent, whereas the estimate of the dipole direction varies mostly owing to the strong influence of the variability of the relativistic jets (Table 1). The individual apparent proper motion of the reference radio sources may reach 0.5 millisecond of arc/yr (i.e. 100 times of the dipole effect). Therefore, the exclusion of astrometrically unstable radio sources from the analysis is a highly labour-intensive process that yields marginal results with a formal error of 1 µas/yr (e.g. Titov & Lambert (2013)). In the meantime, Mignard & Klioner (2012) showed that the Gaia mission will obtain the parameters of the Galactocentric acceleration with a formal error of 0.2 µas/yr, i.e. five times better. Therefore, we are seeking methods to improve the VLBI data results.
Theoretically, the Galactocentric acceleration vector points to the centre of the Galaxy, which is supposed to be the location of the compact radio source Sagittarius A* with the estimated coordinates α G = 267 • for right ascension and δ G = −29 • for declination having the absolute position error of 0 . 12 (Reid & Brunthaler 2004). A comprehensive instruction to calculate the acceleration amplitude using the International Astronomical Unionrecommended Galactic centre distance (R gal = 8.5 kpc) with the Galaxy rotation velocity at R gal (V gal = 220 km/s) and assuming the circular type of motion is given, for example in Kovalevsky (2003) or Titov et al. (2011). This results in a value A = |a| = V 2 gal R gal = 1.85 · 10 −13 km/s 2 , which can be converted to the acceleration vector a with the magnitude of 3.8 µas/yr. There are dozens of papers devoted to revision of the parameters R gal and V gal based on different modern observation data and analysis strategies. For instance, Reid et al. (2009) provided new meanings for R gal = 8.4 ± 0.6 kpc and V gal = 254 ± 16 km/s. The corresponding acceleration amplitude is 2.49 · 10 −13 km/s 2 and the amplitude A of the dipole proper motion equals to 5.1 µas/yr. In one of the most recent publications, de Grijs & Bono (2016) recommended a downward adjustment of the R gal to 8.3 kpc and an upward revision of V gal in the range from 232 km/s to 266 km/s. The corresponding amplitude of the Galactocentric acceleration vector lies therefore within the range (2.10 ÷ 2.76) · 10 −13 km/s 2 followed by the expected amplitude of the dipole proper motion between 4.3 and 5.7 µas/yr. Titov (2011) noted that the conventional equation of the geodetic VLBI group delay (Petit & Luzum 2010, chap. 11) could be altered to accommodate the Galactocentric acceleration. This modification provides a way to estimate the SAD either as a dipole systematics in proper motions or as a systematic effect in the VLBI scale. Section 2 recalls the corresponding Article number, page 1 of 7 arXiv:1802.05347v1 [astro-ph.GA] 14 Feb 2018 A&A proofs: manuscript no. Titov_Krasna mathematical equations and demonstrates the advantage of the latter option. In Section 3, we present the results obtained by applying the new method to analysis of VLBI data during 1979.7-2016.5 together with the SAD estimates from global VLBI adjustments.

Basic equations
The fundamental geodetic VLBI observable, the group delay τ group , is the first derivative of the phase with respect to angular frequency of the incoming signal. A group of eight channels within the frequency band from 8.1 to 8.9 GHz (known as Xband) measures the phases simultaneously to obtain the group delay by the least-squares adjustment. The conventional group delay model approximating the observed VLBI data is given by Petit & Luzum (2010, chap. 11) as where b is the vector of baseline b = r 2 − r 1 , s is the barycentric (solar system barycentre) unit vector of radio source, V ⊕ is the barycentric velocity of the geocentre, w 2 is the geocentric velocity of the second station, c is the speed of light, G is the gravitational constant, M is the mass of the Sun, and R ⊕ is the geocentric distance to the Sun.
The term 2GM c 2 R ⊕ is related to the general relativity effect. The impact of w 2 is small and may be ignored for sake of simplicity. After these alterations equation (1) is given by Equation (2) now comprises only the effects of special relativity of orders V ⊕ c and V 2 ⊕ c 2 . Leaving only the terms of V ⊕ c and using the Taylor series expansion for (1 + s·V ⊕ c ) −1 , the equation (2) comes down to the main aberration effect where the second term represents the correction for the annual aberration Titov (2011) showed that the Galactocentric acceleration vector a can be added to the conventional equation (1) by replacing the barycentric velocity V ⊕ with the sum V ⊕ + a∆t, where ∆t is the time since the reference epoch 2000.0. As a result the corresponding correction for the aberration effect is Since the proper motions affect the source positions in the form s = s + ∆s = s + µ∆t, the vector of proper motion should be represented as a sum of two components µ = µ 1 + µ 2 = s(s · µ) + µ 2 , where µ 2 is the component along the acceleration vector a, and µ 1 has the direction perpendicular to vectors a and s. Then, the additional delay induced by random-like apparent proper motion of a radio source is given by It should be noted that even if the projection of ∆s from equation (4) on the vector s is zero, i.e. s · ∆s = 0, the time delay in equations (5) and (6) describes a projection of ∆s on the baseline vector b. Therefore, in a general case, the delays in equations (5) and (6) as well as both components of the equation (5), τ aberr 1 and τ aberr 2 , are non-zero values.
From the sum of equations (5) and (6), it can be concluded that the first term, proportional to the geometric delay τ geom = − b·s c , is sensitive to the component µ 1 alone. This implies that the estimate of the Galactocentric acceleration with the first term of equation (7) is likely to be more accurate than from analysis of individual proper motion of reference radio sources.
If the model (equation (1)) described the observations perfectly, then the scale factor (F) determined by is equal to unity for all observations, i.e. F ≡ 1 (e.g. MacMillan (2017)). However, if the proper motions are not estimated, then the first term in equation (7) as a part of the group delay (equation (1)) is considered as a contribution to the scale factor. After division in equation (8), the scale factor is given by and manifests itself as a variable parameter depending on the Galactocentric acceleration, radio source position, and the time since a reference epoch. Consequently, the partial derivative of Table 1. Estimates of the SAD parameters from available publications if only a dipole component is estimated. The components of the acceleration vector are given with a 1 , a 2 , a 3 , the amplitude with A, and α G , δ G are the equatorial coordinates of the vector direction.
the group delay w.r.t. the scale factor as a time-dependent effect is given by or it can be expressed as a time-independent effect in the following way:

Analysis of VLBI data
In this paper we analyse almost all VLBI sessions provided by the International VLBI Service for Geodesy and Astrometry (IVS; (Schuh & Behrend 2012)) together with the Very Long Baseline Array (VLBA) Calibrator Survey (VCS) observing sessions VCS1-6 (Beasley et al. 2002;Fomalont et al. 2003;Petrov et al. 2005Petrov et al. , 2006Petrov et al. , 2008Kovalev et al. 2007 (Böhm et al. 2012) and processing of our standard reference solution (P1) followed the International Earth Rotation and Reference Systems Service (IERS) Conventions 2010 (Petit & Luzum 2010) and the technique specific VLBI analysis recommendations. Further details of the parametrisation and analysis set-up for our standard VLBI solution are given in Krásná et al. (2015).

Scale factor
The scale factor is usually not estimated in the routine procedure of geodetic VLBI data analysis. However, there is a possibility to estimate the scale factor either as a single parameter (in a session-wise or global adjustment) or as a specific parameter for each radio source individually within the global solution, similar to the radio source coordinates.

Scale factor as a single parameter
Using equation (10) we obtained the session-wise estimates of the scale factor from the standard solution (Figure 1), where the station positions and velocities were fixed to the a priori values in ITRF2014 . The Earth orientation parameters and positions of the radio sources (no-net-rotation on the defining ICRF2 (Fey et al. 2009(Fey et al. , 2015 radio sources) were estimated in the usual way on a session-wise basis. The green line in the figure shows the smoothed estimates using the local least-squares regression. The annual periodicity visible in the scale factor estimates is predominantly due to the omitted hydrology loading corrections in the station coordinates, as shown by MacMillan (2017), and the weighted mean over the estimated scale corrections of 0.7 part-per-billion (ppb) is in accordance with the definition of the ITRF2014 scale, which is obtained by the arithmetic average of the scale of satellite laser ranging (SLR) and VLBI solutions. The scale discrepancy between these two combined solutions building the ITRF2014 is 1.37±0.10 ppb at epoch 2010.0 . There is an ongoing discussion about the reason for the systematic difference in scale as determined from the two high accurate space geodetic techniques. Recently, Appleby et al. (2016) pointed out that there seems to be a range bias in the case of SLR, which would be responsible for about the half of the obtained scale discrepancy. Figure 2 shows the respective formal errors of the scale factor estimates presented in Figure 1. A rapid decrease of the formal errors can be seen between the early VLBI years and 1995. The most accurate estimates of the scale factor belong to the RDV (Research and Development VLBA) sessions and these are highlighted with the magenta dots.

Scale factor as a specific global parameter for each radio source
In this section we focus on the globally computed scale factor correction, that is from a common adjustment of the VLBI sessions, where we estimated this correction as a specific parameter for each radio source individually (except the 39 special handling sources). We ran a global solution with a standard parametrisation (P1), where corrections to the terrestrial reference frame (i.e. position and linear velocity of VLBI radio telescopes), celestial reference frame (i.e. position of radio sources), and the individual scale factor for each radio source are determined. Using equation (10) for the partial derivative we obtained the scale factor corrections plotted in the upper part of Figure 3 as a function of right ascension (left-hand side) and as a function of declination (right-hand side). As a next step, the conventional equation for the group delay model (Petit & Luzum 2010, chap. 11) is extended by the Galactocentric acceleration vector a following Titov (2011). The SAD effect is reduced from the measurements a priori using the following values: A = 5 µas/yr, α G = 267 • and δ G = −29 • . We refer to this parametrisation as to P2. Applying the P2 we ran a second global solution. The difference in ∆F (lower part of Figure 3) between the standard solution and the solution with reduced SAD effect shows a clear systematic effect arising from the omitted SAD effect in the standard parametrisation. The maximum difference in the right ascension (left-hand side) and declination (right-hand side) coincides consequently with the direction to the Galactic centre and reaches about 0.2 ppb, whereas the weighted mean of the difference is −0.016 ppb.
To dispose of the time dependency of the scale factor we express the partial derivative of the group delay with equation (11). Another two global solutions (with P1 and P2 parametrisation) were computed. The resulting annual drift of the scale factor corrections ∆F/∆t from P1 is plotted as a function of equatorial coordinates in the upper part of Figure 4 and the lower part shows the difference between the P1 and P2 parametrisation. Comparing the Figures 3 and 4 it is obvious that the scale factor corrections, expressed as a time-independent effect in the [ppb/yr] units, are free of large outliers. The annual variation of the scale factor due to the omitted SAD effect reaches about 0.02 ppb/yr, which reflects the mean observing time of about 10 years over all included radio sources.
Since we estimated the source-wise scale corrections together with the terrestrial (TRF) and celestial (CRF) reference frame within the global adjustment, a constant bias in the scale difference propagates also to the estimated terrestrial reference frame. Comparison of the TRF from the parametrisation P1 minus TRF from P2 in terms of the 14-parameters weighted Helmert transformation gives a negligible difference of −0.011± 0.001 ppb (−0.07 ± 0.01 mm) for the scale. It means that omitting the Galactocentric acceleration effect in the VLBI analysis does not have any significant effect on the global scale of the estimated TRF.

Galactocentric acceleration vector from source-wise scale factor
The individual scale factor corrections ∆F/∆t could be fitted by the following model: F = a 1 cos α cos δ + a 2 sin α cos δ + a 3 sin δ (12) using the components of the acceleration vector a 1 , a 2 , a 3 where A is the amplitude of the Galactocentric acceleration and α G , δ G are the equatorial coordinates of the vector direction. The three parameters determining the dipole (a 1 , a 2 , a 3 ) are to be estimated using the standard least-squares method. Then the amplitude A and the coordinates α G , δ G are given by α G = arctan a 2 a 1 ; δ G = arctan a 3 a 2 1 + a 2 2 .
(15) Table 2 shows the estimated dipole components from fitting the annual variation of the scale factor for different cut-off limits for the number of observations N of each radio source in order to mitigate the effect of rarely observed radio sources and to verify the stability of the solution with respect to N. For a wide range of N (from 50 to 20 000), where the compromise between the number of observations and the number of sources is fulfilled, the amplitude of the dipole effect varies between 4.8 and 5.3 µas/yr, and the best formal error is 0.2 µas/yr, i.e. much less than from the analysis of proper motions Titov & Lambert 2013. One can conclude that this approach is less sensitive to the peculiar behaviour of the active galactic nuclei. The estimated direction in declination (from −28 • to −35 • ) gets close to the value estimated by Reid & Brunthaler (2004) within its formal error, whereas there is a offset of 14 • between the estimated direction in right ascension (∼ 281 • ) and the value 267 • reported by Reid & Brunthaler (2004). The bias can be explained by the apparent proper motion of astrometrically unstable radio sources. This instability is caused by internal motion of matter in relativistic jet.
As shown in Figure 1 the time series of the scale factor estimates contains an annual signal originating predominantly from the omitted hydrology loading corrections. Therefore, we computed another solution based on P1, where we applied the hydrology loading displacement time series provided by the NASA Goddard Space Flight Center VLBI group (Eriksson & MacMillan 2014) on the station coordinates a priori. We ran again a global solution where the annual drift of the individual scale factor corrections was estimated in the same way as in the solutions mentioned above. From fitting the estimates belonging to radio sources with more than 50 observations, the following dipole components were obtained: A = 5.1 ± 0.2 µas/yr, α G = 279 • ± 3 • and δ G = −35 • ± 3 • . This yields a difference of 0.1 µas/yr in A and 2 • in α G with respect to the solution based on parametrisation P1, which lies within the formal errors of the estimates and shows that the seasonal hydrology variation in the station coordinates does not have any significant effect on the estimated dipole components.

Galactocentric acceleration vector as a global parameter
The dipole components can be also estimated as global parameters in a VLBI solution in the form of the Galactocentric acceleration vector. The parametrisation P2 allowed us to build the partial derivative of the group delay with respect to the three components (a 1 , a 2 , a 3 ) of a. We ran a global solution in which corrections to the terrestrial reference frame (i.e. position and linear velocity of VLBI radio telescopes), celestial reference  (2014)). The direction in declination of the vector coincides with the value estimated by Reid & Brunthaler (2004) if only the large network sessions were included in the solution, which implies that this procedure is sensitive to the inclusion of weak networks.

Conclusions
We introduce a new method to detect the secular aberration drift induced by the Galactocentric acceleration from the geodetic VLBI measurements. It is based on fitting the scale factor corrections estimated for each source individually within a global solution. In Titov & Lambert (2013) we used the individual proper motion of reference radio sources to visualise the dipole effect caused by the Galactocentric acceleration. However, because of contamination by large apparent motions induced by relativistic jet, the proper motions have large random noise almost hiding the dipole systematics. However, if the dipole components in proper motion are not estimated, the systematic effect comes to the scale effect that, in contrast to the proper motion, is free of the relativistic jet problems. From fitting the individual scale factor corrections of sources with more than 50 observations during 1979.7 -2016.5 we obtained a Galactocentric acceleration vector with an amplitude of 5.2 ± 0.2 µas/yr and direction α G = 281 • ±3 • and δ G = −35 • ±3 • . The SAD was also estimated directly within a global adjustment of the VLBI data. The Galactocentric acceleration vector determined from the selected large network IVS sessions after 1993 (A = 5.4 ± 0.4 µas/yr, α G = 273 • ± 4 • , δ G = −27 • ± 8 • ) is closer to the value reported by recent papers (e.g. Reid & Brunthaler (2004) and de Grijs & Bono (2016)) than the estimate from the entire VLBI history. The formal error of the amplitude of the acceleration vector estimated by means of the scale factor is about five times better than the estimate from the proper motions (Titov & Lambert 2016) and comparable to the accuracy expected to be obtained with the Gaia mission (Mignard & Klioner 2012). This demonstrates the advantage of the new approach with respect to the traditional procedure.
Our result also points out the potential impact of neglecting the Galactocentric acceleration on the scale factor and, finally, on the geodetic positions of the radio telescopes. The scale factor becomes dependent on the positions of the reference radio sources, i.e. larger in the direction of the Galactic anti-centre, and smaller in the direction of the Galactic centre. In a global sense, the quality of the geodetic results is worsening and, moreover, the loss of quality would increase linearly with time. The results presented in this paper were computed with the VieVS software and verified with the OCCAM software package.