Comparing Vienna CRF solutions to Gaia-CRF2

We are using various models and analysis strategies, such as galactic aberration, ray-tracing etc., to create different Vienna celestial reference frame (CRF) solutions. These solutions are then compared against the Gaia reference frame (Gaia-CRF2). This is done using a degree 2 vector spherical harmonics approach. The estimated parameters are used to investigate the impact of the various analysis methods on the differences between Gaia and the Very Long Baseline Interferometry (VLBI) CRF. We ﬁnd that correcting for galactic aberration reduces the difference between the Gaia-CRF2 and the VLBI CRF signiﬁcantly (30µas in D 2 and 13µas in D 3 ). Furthermore, we ﬁnd that using a priori ray-traced tropospheric delays in addition with low absolute constraints on tropospheric gradients reduces the a e20 parameter by 20 µas. Using these analysis strategies we can explain almost all signiﬁcant differences between the Gaia-CRF2 and the VLBI CRF. However, the vector spherical harmonic (VSH) parameter a e20 is still highly signiﬁcant and can not be explained by modeling and analysis choices from the VLBI technique.


Introduction
The rotation about the Galactic center causes an acceleration of the Solar System Barycenter (SSB) towards the center of the Galaxy. The galactic aberration (GA) is the aberration of positions of distant objects resulting from the revolution of the SSB about the Galactic center with a period of 250 million years. Over decades of observing, it imprints an apparent source proper motion of a few µas per year.
Source positions estimated from geodetic Very Long Baseline Interferometry (VLBI) are on an accuracy level where an effect of this magnitude can be calculated and consequently has to be corrected in the analysis. Several groups estimated the GA from VLBI data (see Titov et al. 2011;Xu et al. 2012;Titov and Lambert 2013;Titov and Krásná 2018). The reported values range from 5:2˙0:2 to 6:4˙1:1 µas per year with the center of the Galaxy at 17 h 45 min 40 s in right ascension and 29 ı 00 0 28 00 in declination.
The International VLBI Service for Geodesy and Astrometry (IVS) established a working group with the general purpose of investigating the issues concerning the incorporation of the GA in IVS analysis, see MacMillan et al. (2019). This working group agreed upon a value of 5:8 µas per year which was also applied in the creation of the ICRF3, see Charlot et al. (2020).
One of the largest error sources in geodetic VLBI is the troposphere. It was demonstrated by Mayer et al. (2017) that using a priori ray-traced tropospheric delays in a global CRF solution significantly influences the source coordinates. Further, they showed that constraining tropospheric gradients to their a priori values influences the source coordinates as well.
The second data release from the Gaia satellite, Gaia Data Release 2 (Gaia DR2), includes a celestial reference frame (Gaia-CRF2) of comparable accuracy to VLBI, see Lindegren et al. (2018) and Mignard et al. (2018) for more information. Since Gaia was launched at the end of 2013 the effect of GA on its positions is negligible. Further, Gaia is a satellite in space and, therefore, unaffected by tropospheric disturbance. Additionally, the satellite's rotation and precession (the so called scanning law) are designed to maximize the uniformity of the sky coverage. The scanning law introduces some systematic effects, presented in Lindegren et al. (2018), however, they can most likely be ignored when large scale global systematic effects are concerned.
VLBI suffers from an uneven network distribution, which could result in a global deformation of the frame when the troposphere is not modeled correctly. Further, the unmodeled effects of GA introduce a global systematic deformation of the frame over the years. Since the Gaia-CRF2 is not affected by these effects it provides a perfect independent source for external validation of the effects of GA and the troposphere.

DataandAnalysis
We generated five global geodetic VLBI solutions with different analysis options and compared them to the Gaia-CRF2. The basic solution setup is described in the beginning of this section. Other solutions are based on the basic solution with some changes in modeling and analysis described at the end of this section. All VLBI solutions presented here utilize all geodetic VLBI sessions that were used for the creation of the ICRF3. In total this data set includes about 13 million observations of about 4,500 sources. The VLBI CRF solutions were generated using the Vienna VLBI and Satellite Software (Böhm et al. 2018), which is developed by TU Wien. These solutions follow the IERS Conventions 2010 by Petit and Luzum (2010) for reducing the observations and geophysical modeling. Also, antenna thermal deformation (Nothnagel 2009) and atmospheric pressure loading (Wijaya et al. 2013) were taken into account.
A priori positions for stations (including velocities) and sources are taken from the ITRF2014 (Altamimi et al. 2016) and ICRF2 (see Ma et al. 2009;Fey et al. 2015) respectively. A standard geodetic analysis is performed. This results in an updated ITRF and ICRF as well as EOP time series. Station Reference solution Solution 2 Solution 1 + a priori ray-traced delays Solution 3 Solution 2 + removing of absolute constraints on tropospheric gradients Solution 4 Solution 3 + GA model Solution 5 Solution 4 + error scaling coordinate adjustments were estimated as global parameters with No-Net-Rotation and No-Net-Translation conditions applied to the positions and velocities of a group of 21 stations. Most of the radio source adjustments were estimated as global parameters after a No-Net-Rotation constraint is applied to the positions of the 295 ICRF2 defining sources. The coordinates of the 39 special handling sources (Fey et al. 2015) were, however, estimated once per session. A priori hydrostatic zenith delays were determined from local pressure values and then mapped to the elevation using the Vienna Mapping Functions 1, see Böhm et al. (2006). Additionally, tropospheric gradients from the NASA/GSFC Data Assimilation Office (DAO) model (see MacMillan 1995;MacMillan and Ma 1997) are utilized for all solutions. The wet zenith delays, north and east troposphere gradients, and clock values were estimated every 30 min, 6 h and 1 h, respectively. Tropospheric gradients are constrained to their a priori values. This was realized with absolute and relative constraints of 0.5 mm (after 6 h), which removes unrealistic gradient estimates but affects the declination of estimated sources. Other solutions are based on the same parameterization with slight amendments, see Table 1 for an overview. In the second solution a priori ray-traced tropospheric delays were included. Since the absolute constraints on a priori tropospheric gradients tend to influence source declination we also created the same solution where we removed these constraints. In the fourth solution GA was corrected a priori with the recommended value of 5.8 µas/year with the center of the Galaxy at 17 h 45 min 40 s in right ascension and 29 ı 00 0 28 00 in declination. As a reference epoch we chose 2015.0 since the Gaia positions epoch is close to this epoch. In VieVS the correction was realized by modification of the conventional group delay equation as proposed by Titov et al. (2011). Additionally, we generated a solution where we scaled the errors to more realistic values. This was done using a scaling factor of 1.5 and a noise floor of 30 µas, see Charlot et al. (2020) for more information on these values.
As a reference a subset of sources from the Gaia-CRF2 solution described in Mignard et al. (2018) was used. This subset consists of the positions of 2,820 sources which have an ICRF3 counterpart. A detailed analysis of the differences of the catalogue and an ICRF3 prototype solution can be found in Mignard et al. (2018) and Petrov et al. (2018). Methodology The resulting catalogs are compared to the Gaia-CRF2 using vector spherical harmonics (VSH) as described in Mignard and Klioner (2012). Global features of the differences such as a rotation of the catalogs and the so called glide parameters are reflected in degree 1. Degree 2 describes the quadrupole deformations between the catalogs. The whole transformation reads where R i are the three rotation parameters, D i are the three glide parameters and a m;e lm are the quadrupole parameters of electric (e) and magnetic (m) type.
The data set used here is rather small (about 2,800 sources) and we are only interested in global effects. Therefore, we decided to stop the VSH expansion at degree 2. Outliers were eliminated, see Sect. 3.1, and transformation parameters were estimated using the classical least squares method. Variances and correlations are used to weigh the differences.

Outlier Detection
There are many options to choose from when eliminating outliers. Using the normalized separation was proposed by Mignard et al. (2016). This approach takes the correlation between right ascension and declination into account and introduces an arbitrary cut off of 10 mas (angular separation). However, we expect to have differences in positions due to galactic aberration which are on the level of the formal errors. This is especially critical for sources, which have a long observing history with VLBI. These sources tend to be most affected by GA. Further, their long observing history implies that they have been observed many times and, therefore, have small error bars. The outlier test proposed by Mignard et al. (2016) is based on the assumption that the differences follow a Rayleigh distribution. However, it was found by Petrov and Kovalev (2017) that the differences between Gaia-CRF2 and VLBI CRF deviate from a Rayleigh distribution. For the reasons mentioned above we decided to not use the outlier test proposed by Mignard et al. (2016) but rather use our own method of outlier detection.
Since the number of intersecting sources is rather small the outlier detection can be realized in the parameter space. This was accomplished by removing each source once from the standard solution and calculating the VSH parameters. At the end we have a set of about 2,800 VSH parameters which can be used to calculate the standard deviation of each VSH parameter. Outliers can then be found with a simple three sigma cut off for each parameter. The source is excluded when one of the VSH parameters is outside of this cut off. This means that a source that significantly changes one of the VSH parameter by itself is considered an outlier. With this approach we find that about 7% of sources are flagged as outliers. This list of outliers is then used for every solution presented here. In order to compare the different transformation parameters we decided to stick to one list of good sources. This was realized by doing the outlier elimination on the standard solution and then using this list of outliers for each other solution.
When plotting the detected outliers one can see that the outlier test proposed by Mignard et al. (2016) clearly removes some kind of systematic effect, see right part of Fig. 1. When removing the outliers with our technique the systematic is not that clear but might still be there, see left part of Fig. 1. It is important to note here that other outlier elimination techniques that were tested showed similar results.

Results and Discussion
During our research we created five global VLBI solutions. For each solution one analysis option was changed, see Table 1. However, the previous changes were not revoked but rather applied alongside. The VSH parameters of these five solutions with respect to Gaia-CRF2 are depicted in Fig. 2. The center of the bar represents the VSH estimate while the length of the bar reflects the formal uncertainties of the estimate. Different solutions have different color codes and different names. Solution 1 reflect the standard Vienna CRF solution, as described in Sect. 2. This solution experiences significant deviations from Gaia-CRF2. Solution 2 is similar to solution 1 with the difference that a priori ray-traced tropospheric delays are used in the analysis. Using this model succeeds in reducing the a e 20 parameter. However, at the same time the D 3 parameter is increased. The D 3 parameter can be decreased again by loosening the constraints on the tropospheric gradients, which is reflected by solution 3. We further succeed in reducing the D 2 and D 3 parameter when adding the GA model, see solution 4. Scaling the formal uncertainties (solution 5), which is a mandatory task when creating a VLBI CRF, does not change the parameters significantly. Note that correlations between the VSH parameters are relatively low with a maximum of 0:37 between D 2 and R 1 .
When looking at Fig. 2 the deformations of degree 1 and degree 2 are particularly interesting, since they describe real differences between the frames. The rotation between Gaia-CRF2 and the VLBI CRF is statistically significant, but less interesting because they do not reflect real effects. During the creation of Gaia-CRF2 the frame was rotated onto an ICRF3 prototype solution, therefore, no rotations should be present between the two frames. However, the number of sources used for the rotation differs. The Gaia-CRF2 was rotated onto 2,844 matching ICRF3 sources while we use 2,588 sources after outlier detection. Further, the sources selected as outliers are the ones affecting the parameters most, therefore, a large rotation can be expected. However, the meaning of this rotation is questionable and will not be further discussed. Note that the ICRF3 prototype solution was rotated onto the same 295 ICRF2 defining sources that were used for the solutions discussed here.  Table 1, is depicted in green. In red the standard solution with a priori ray-traced delays is depicted. This is solution 2. The blue solution (solution 3) illustrates the effect of removing absolute constraints on tropospheric gradient estimation. In the magenta solution the GA model was applied, see solution 4. The final solution (solution 5), where errors are scaled, is depicted in black Using the GA model clearly removes systematic effects between the two frames (difference of 30 µas in D 2 and 13 µas in D 3 ). This can be explained by the fact that the two techniques have very different time spans. The first VLBI data used for the creation of the VLBI CRF dates back to 1979, which is a total of 40 years of data. This means the effects of GA had time to accumulate, which in turn effects the source positions calculated with VLBI. The Gaia satellite was only launched in 2013. Therefore, GA did not have enough time to affect source positions calculated from Gaia data. This means that a VLBI CRF corrected for GA (with a reference epoch close to Gaia) should agree better with the Gaia-CRF2 than a VLBI CRF not corrected for GA. This is exactly what we see in Fig. 2.
Using a priori ray-traced tropospheric delays succeeds in reducing the a e 20 parameter, which is the most significant deformation between the frames, by 20 µas. This is particularly important, since no other model or analysis choice affects this parameter. Unfortunately, the result is not that clear because the D 3 parameter is increased by about the same amount at the same time. However, we found that the D 3 parameter, which is directly connected to the source declination, is highly susceptible to models and analysis choices. Reducing the absolute constraints on tropospheric gradient estimation, for example, succeeds in reducing the D 3 parameter by 20 µas.
The results from Fig. 2 can also be plotted on the sphere for easier interpretation. Figure 3 depicts the deformations of the standard Vienna CRF solution, marked as green in Fig. 2. Plot (a) depicts the glide parameters (D 1 , D 2 and D 3 ) with the addition of the galactic plane (red dashed line), the galactic center (black circle) and the ecliptic (black dashed line). One can clearly see the dominant effect of GA with arrows roughly from the galactic center to the anti center. Plot (b) depicts the deformations of degree 2, which are dominated by the a e 20 parameter. Figure 3 can now be compared to Fig. 4 which depicts the deformations of the final best fitting solution, which is marked in black in Fig. 2. One can see that deformations of degree 1 are insignificant. Further, the deformations of degree 2 are decreased but still significant.
For completeness the VSH are also calculated using the outlier detection method proposed by Mignard et al. (2016), see Fig. 5. The error bars and estimates are much smaller, when using this technique. This can be explained by the number of outliers found by the different techniques. With the technique proposed by Mignard et al. (2016) about 500 sources are flagged as outliers while only about 200 are flagged by the outlier elimination technique described in this paper. In this comparison the outliers were detected for each solution separately. Unfortunately, also the effect of the models is slightly different. Using a priori ray-traced tropospheric delays does also affect other deformations of degree 2. The correction of GA does move the parameters by roughly the same amount. However, since the parameters are small to begin with the correction moves the parameter into the negative. Scaling the formal uncertainties seems to affect the parameters when using this outlier elimination technique. This demonstrates that the VSH parameters between the VLBI CRF and Gaia-CRF2 are very susceptible to the outlier elimination technique used.

Conclusions and Outlook
We produced several VLBI CRF solutions with different models and analysis choices. These solutions were then compared to the Gaia-CRF2 using VSH of degree 2. We find three (D 2 , D 3 and a e 20 ) significant parameters between the standard Vienna CRF solution and the Gaia-CRF2. The a e 20 can be reduced by 20 µas by using a priori ray-traced tropospheric delays in the analysis. However, this also increases the D 3 parameter by the same amount. We find that the D 3 parameter, which is directly connected to the source declination, is very susceptible to models and analysis choices. The parameter can be significantly reduced (in this case about 20 µas) when lowering the absolute constraints on tropospheric gradients. Further, we find that using GA succeeds in reducing the D 2 parameter by 30 µas and D 3

Fig. 5
The figure is similar to Fig. 2 with the exception that the outlier elimination technique proposed by Mignard et al. (2016) is used to determine the VSH parameters parameter by 13 µas. We can explain all the deformations of degree 1 between the VLBI CRF and the Gaia-CRF2 with analysis choices made by the VLBI analyst. However, the a e 20 can only be reduced but not fully explained by choices made by the VLBI analyst. Furthermore, we find that using the normalized separation method proposed by Mignard et al. (2016) to detect outliers removes sources that are needed to see the systematic differences between the VLBI and Gaia frame.
A more in depth analysis with more VLBI CRF solutions can be found in Mayer (2019).