© ESO, 2016

1. Introduction

Gaia is an astrometric satellite launched in late 2013 and designed to produce a three-dimensional map of the local part of our Galaxy from which our Galaxy’s composition, formation, and evolution can be reconstructed. The satellite will measure the positions, proper motions, and parallaxes of at least one billion stars in the Milky Way. In addition, Gaia will detect thousands of exoplanets, asteroids, and about half a million distant quasars in the optical spectrum. Since the accuracy is at the micro-arcsecond level, these precise measurements help to improve the optical realisation of the International Celestial Reference System (ICRS) by a few orders of magnitude compared to its current realisation by the Hipparcos and Tycho catalogues (ESA 1997). It will also provide a number of new tests of the general theory of relativity.

The ICRF (Ma et al. 1998) is a quasi-inertial reference frame that was originally defined by the measured positions of 212 extragalactic radio sources derived from ground-based very long baseline interferometry (VLBI) and has its reference point at the barycentre of the solar system. In general relativity there is no true inertial reference frame; the extragalactic sources (quasars) used to define the ICRF are so far away, however, that any net angular motion is almost zero. The ICRF is now the standard reference frame used to define the positions of astronomical objects. It has been adopted by International Astronomical Union in 1998 (Ma et al. 1998). In 2009, the second realisation, ICRF2 (Ma et al. 2009; Fey et al. 2015), was adopted. This included improved models and concepts and was based on 3414 compact astronomical sources. ICRF2 defines the reference frame to an accuracy of 40 μas and includes 295 defining sources uniformly distributed on the sky and selected on the basis of positional stability and the lack of extensive intrinsic source structure.

The Hipparcos and Tycho catalogues currently serve as the corresponding optical realisation of the International Reference System (ICRS), but they will be superseded by the Gaia mission in the coming years (Mignard 2011). This frame must be constructed with the same principles as the ICRS, that is, overall a kinematically non-rotating system with the same orientation as the radio ICRF. It is estimated that Gaia will astrometrically measure some 500 000 quasars, and their repeated measurement over the estimated five year mission will lead to a new kinematically defined inertial reference frame in the optical.

Quasars located at the centre of distant active galaxies are characterised by extremely compact and bright emission. They are at cosmological distances and therefore show negligible transverse motion. However, recent observations (Taris et al. 2011; Porcas 2009; Kovalev et al. 2008) of active galactic nuclei (AGNs) and theoretical studies (Popović et al. 2012) indicate that variability in the accretion disk and dusty torus surrounding the central black hole can cause photocentre shifts of up to the milliarcsec level, therefore it is probable that quasar variability will affect the reference frame as well. This paper investigates the statistical impact of variability on the Gaia reference frame based on simulated astrometric observations.

Another interesting aspect is the acceleration of the solar system that was first pointed out by Bastian (1995). It causes a pattern of proper motions that can be solved for while determining the reference frame orientation and spin parameters. We present an estimate of how well the acceleration vector can be determined based on realistic Gaia simulations that also assess the impact of quasar variability. In addition to this effect, there may be several more subtle effects present in the real data. Many of these will be weak; nevertheless, it is interesting to review and compare them, and if they are deemed strong enough, we can simulate them together with the acceleration of the solar system to assess their impact on the solution and the potential to distinguish the different effects.

The organisation of the paper is as follows: Sect. 2 presents various proper motion effects that could have an impact on the Gaia measurements, and we quantify their respective magnitudes. In Sect. 3 we present the numerical simulations. In Sect. 4 we explain how the Gaia reference frame is determined together with the patterns of proper motion discussed in Sect. 2. In Sect. 5 the simulation results are presented and discussed, and finally in Sect. 6 conclusions are given.

2. Proper motion effects in the Gaia data

Titov (2010) noted that individual apparent proper motions of distant radio sources are generally attributed to the internal structure of AGNs. However, he pointed out that there are a number of other effects that could give rise to systematic apparent motions of quasars, and we briefly consider these here. The various sources of proper motion are summarised in Table 1. We note that the terminology used in the literature describing the different effects varies and can be inconsistent, which is mainly due to the historical progression of the topics.

2.1. Photocentre variability

Perryman et al. (2014) considered the various sources of quasar instability. They concluded that the most important effect is optical photo-centric motion, where the internal structure of the AGN could result in variability of typically less that 60 μas but up to 100 μas in extreme cases (Taris et al. 2011). In the present simulations we considered the extreme case where these effects combine to produce maximum variability distortions of 100 μas. The timescales of quasar variability have been studied by Smith et al. (1993) and Taris et al. (2011) and were found to range from 3–15 years, peaking between 6–9 years. We chose to use 2 and 10 years, which roughly represents the range of timescales to which Gaia is most sensitive, assuming a mission duration of 5 years.

2.2. Acceleration of the solar system

The solar system orbital velocity of ~200 km s^-1 around the Galactic centre results in an aberration effect of about 2.5 arcmin (Perryman et al. 2014) in the direction of motion. The acceleration of the solar system in its Galactic orbit causes this effect to change slowly, which results in a proper motion pattern for all objects on the sky (Bastian 1995; Kovalevsky 2003; Gwinn et al. 1997; Sovers et al. 1998; Klioner 2003; Kopeikin & Makarov 2006; Malkin 2014). This effect has also been referred to as secular aberration drift (Titov 2010). The effect is generally assumed to be towards the Galactic centre where most of the mass is concentrated, but of course the unknown distribution of dark matter may affect the direction. The Galactocentric acceleration can be calculated as a = v²/r. For example, assuming a circular velocity v = 223 km s^-1 and radius r = 8.2 kpc (Ghez et al. 2008), the resulting proper motion pattern is (1)where c is the speed of light and θ is the angular distance between the object and the Galactic centre (Bastian 1995).

The local group and local supercluster also accelerate the barycentre of the solar system similar to the Galactocentric acceleration. We estimated ã for some 33 galaxies of the local group, using a = GM/d² . The seven largest examples are listed in Table 2. The values of ã are very low (<0.1 μas yr^-1), and Gaia will not be able to measure them. We similarly calculated the value of ã for some clusters of the local supercluster, and in this case the values were found to be very low as well and not measurable by Gaia.

If we consider some non-axisymmetric potential for the Milky Way that is due to the bar and the spiral arms, then the acceleration vector will not point exactly towards the Galactic centre. For example, in the potential used by Feng & Bailer-Jones (2013), we find the acceleration vector to be offset by about 2°.

2.3. Cosmological proper motion

The instantaneous velocity of the solar system with respect to the cosmological microwave background (CMB) will cause extragalactic sources to undergo an apparent systematic proper motion. The effect is referred to as cosmological or parallactic proper motion (Kardashev 1986). This effect depends on the redshift z, and fundamental cosmological parameters can in principle be determined from the motion. The velocity (v) of the solar system with respect to the observable Universe produces a dipole pattern in the CMB temperature with ΔT/T = v/c. Observations of the CMB indicate that v = 369 ± 0.9 km s^-1 in the apex direction with Galactic longitude l = 263.99° and latitude b = 48.26° (Hinshaw et al. 2009; Planck Collaboration XXVII 2014). This motion should produce a parallactic shift of all extragalactic objects towards the antapex at the angular rate (2)where β is the angle between the source and apex directions, and(3)is the transverse comoving distance, which for a flat universe equals the line-of-sight comoving distance, illustrated in Fig. 1. (Hogg 1999; Hobson et al. 2006; Weinberg 2008). The quantities Ω_m, Ω_r, Ω_k, and Ω_Λ are the dimensionless energy densities of matter, radiation, curvature, and cosmological constant, respectively. For very low redshift (z ≪ 1), μ ≃ (vH₀/cz)sinβ, or 1–2 μas yr^-1 for z = 0.01 (Kardashev 1986). In principle, this could be within the reach of Gaia, depending on the number of available objects and the precision of the observations.

Both the acceleration of the solar system and the cosmological proper motion give rise to dipole patterns in the proper motions of distant objects. However, the former does not depend on the redshift, while the latter does, which makes it possible to separate the effects.

Fig. 1 Variation of comoving distance with redshift.

2.4. Primordial gravitational waves

Primordial gravitational waves could give rise to systematic proper motions over the sky, composed of second-order transverse vector spherical harmonics (Gwinn et al. 1997). Changes in the space time metric that are due to gravitational waves alter the optical path length, preserve the sources brightness and intensity, but produce oscillations in the apparent position of the source. If the interval of observation is shorter than the period of the gravitational wave, then this will be seen as a systematic proper motion on the sky. Pyne et al. (1996) and Gwinn et al. (1997) have investigated the low frequency observational constraints on gravitational waves in detail, which could arise naturally in inflationary theories of cosmology. It is anticipated that Gaia is not sufficiently sensitive to detect primordial gravitational waves directly, but it should be able to place upper limits on their energy density, which is comparable to that of pulsar timing measurements (Book & Flanagan 2011).

2.5. Transverse redshift drift/cosmic parallax

Kristian & Sachs (1966) and Titov et al. (2011) pointed out that an anisotropic expansion of the Universe would result in a distortion, as a function of redshift, of all distant objects on the celestial sphere in a particular direction. The effect gives a direct measurement of space time curvature, which is similar to a gravitational lens, but in this case is due to the cosmological curvature and not to a single body. The time-dependent components of the distortion would result in patterns of proper motion that could be a function of the redshifts and can be measured in principle. Titov et al. (2011) has shown that the dipole term does not vary significantly and agrees with the predicted estimates of 4–6 μas yr^-1. For the quadrupole anisotropy, which could be interpreted as an anisotropic Hubble expansion or as an indicator of primordial gravitational waves, no detection has been made.

Quercellini et al. (2009, 2012) used the term cosmic parallax for the varying angular separation between any two distant sources, caused by the anisotropic expansion of the Universe. They considered two different scenarios:

• Bianchi 1 models in which the observer is centrally embedded in an intrinsically anisotropic expansion of the early Universe. In this case, the anisotropic stress of dark energy can induce an anisotropic expansion of the Universe at late times that cannot be constrained by the CMB background measurements. Their simplified model, assuming an anisotropy of about 1%, gives a proper motion pattern of about 0.2 μas yr^-1.

• Lemaitre-Tolman-Bondi (LTB) void models in which the observer is off-centre and the Universe is inhomogeneous and isotropic. For the models considered by Quercellini et al. (2012), they derived effects of about 0.02 μas yr^-1.

2.6. Peculiar proper motion

In analogy to the peculiar motion of stars, which is defined as the star velocity relative to the local standard of rest, the peculiar motion of a galaxy is its velocity relative to the Hubble flow. The transverse component of this motion causes a proper motion of the galaxy that we call peculiar proper motion. The typical transverse velocity (v_pec) is of ~600 km s^-1, but in rich galaxy clusters it may be as high as ~1500 km s^-1 (Sparke & Gallagher 2000). The proper motion is then given by (4)where d(z) is the comoving distance Eq. (3)and μ_pec ≃ v_pecH₀/cz for low redshifts. This gives a peculiar proper motion in the range 3–7 μas yr^-1 at z = 0.01. The effect is random and is not expected to give a systematic pattern of proper motion.

2.7. Quasar microlensing

Belokurov & Evans (2002) pointed out that microlensing events can be detected by Gaia both from the brightening of the source star (photometric microlensing) and the positional excursion (astrometric microlensing). The optical depth is given by them as 10^-7 and 2.5 × 10^-5 for photometric and astrometric microlensing, respectively, which corresponds to about 1300 photometric and 25 000 astrometric events during the course of the nominal five-year Gaia mission. The individual images of microlensed sources cannot be resolved by Gaia, but the centroid shift along the trajectory of a source can be measured. In the case of stellar microlensing of quasars, the centroid shift can be tens of μas (Lewis & Ibata 1998).

Sazhin et al. (2011) derived detailed expressions for the apparent proper motions caused by weak microlensing. In Appendix A we show that the statistical effect is related to the optical depth for photometric microlensing. If we consider that a source is a quasar with zero proper motion (μ_Q = 0), then by considering a large number of lensing objects with random proper motions, we find that (5)where (μ_Q′)_RMS is the root mean square (RMS) value of the quasar proper motion, (μ_L)_RMS is the RMS value of the lens proper motion, and τ is the probability that the source Q is within the Einstein radius of some lens, which is simply the optical depth for photometric microlensing.

If we assume that the number of galaxies is 0.05 Mpc^-3, the lens mass is 10¹¹M_⊙, the lens distance is 5 Gpc, and the source distance is 10 Gpc, then the value of τ ~ 0.08 in the extragalactic case. Similarly, in case of Galactic microlensing, if we take the number of sources as 0.1 pc^-3, the mass as 1 M_⊙, the lens distance as 5 kpc, and the source distance as 10 kpc, then τ ~ 10^-6. Clearly, (μ_Q′)_RMS ≪ 1 μas yr^-1 in both cases. This means that even though a distant galaxy acting as a lens may have a tiny motion in the transverse direction with respect to the observer, the source quasar being at rest, the apparent proper motion of the quasar is too weak to be measured by Gaia, and the same applies to Galactic microlensing.

2.8. Other apparent motions of quasars

The peculiar proper motion of quasars is expected to be very small because of their large distances. For example, for a quasar at z = 3, μ_pec = 0.02 μas yr^-1 according to Eqs. (3)and (4). This proper motion is very small and hence the quasars can in this respect be assumed to be stationary. However, radio observations of quasars have shown that a significant number of them have apparent proper motions exceeding 50 μas yr^-1 (MacMillan 2005; Titov et al. 2011). Sazhin et al. (2011) have pointed out several possible causes, including apparent superluminal motions in radio jets, gravitational waves and weak microlensing by stars and dark bodies in our Galaxy. Sazhin et al. (2011) has given a number of examples of apparent motions that are due to microlensing, and the order of proper motions in some cases has been shown to be several tens of μas yr^-1. However, the number of such events is estimated to be very small.

Fig. 2 Markov chain photo-centre variability for a VLBI quasar with timescales, τcor, of 2 (left) and 10 years (right). In each case three examples are given for the same sources.

3. Simulations

3.1. Data sets

To simulate the quasars and galaxies, we took the positions, magnitudes, and redshifts of the objects from the GUMS (Robin et al. 2012) simulated data set. This data set includes a realistic sky, magnitude, and redshift distribution of about 1 million quasars and 38 million unresolved galaxies¹ expected to be seen by Gaia. As a significant fraction of the galaxies will not be observed by Gaia because of their extended structure, we focus on the simpler and most point-like objects down to G = 20 magn by selecting only ellipticals and spirals (Hubble types E2, E-S0, Sa, and Sb with a bulge-to-total flux ratio of 1.0, 0.9, 0.56–0.57, and 0.31–0.32, respectively) in the redshift range 0.001 to 0.03. This results in just over 100 000 objects, which is sufficient for our simulations.

The GUMS data set of nearly 1 million quasars represents an idealised case of those objects that in principle could be used for the frame determination assuming that they can be correctly classified by means of the Gaia observations (Bailer-Jones et al. 2008). We randomly selected half a million quasars from the full data set for our simulations.

As a very conservative alternative, we also considered using the Initial Gaia Quasar List (IGQL; Andrei et al. 2009, 2012), which currently consists of about 1.2 million quasars compiled from various catalogues. It will be used in the early Gaia data processing to simplify the identification of quasars. IGQL contains a snapshot of the best information of optical quasars available just before the launch of Gaia. The quasars in the IGQL are not uniformly distributed on the sky, but there are a number of bands and high-density regions corresponding to the various surveys used to compile it (Véron-Cetty & Véron 2010; Souchay et al. 2012; Andrei et al. 2009, 2012; Shen et al. 2011; Pâris et al. 2014). Notably, the IGQL is significantly lacking sky coverage in the southern hemisphere, and the impact of this non-uniform distribution is discussed in Sect. 5.

For both quasar data sets, a random subset of about 3300 objects were assumed to have accurate VLBI positions that simulate the positions in the ICRF2 catalogue.

Fig. 3 Left: astrometric weight (w) of GUMS galaxy observations (see Hobbs et al. 2009), where , is the total statistical weight of the assumed along-scan standard error of Gaia observations. Below is the number of GUMS galaxies as a function of G magnitude. The two left plots show that the highest astrometric weight is from G = 12 to G = 16, but the number of galaxies peaks around G = 17. Right: redshift distribution of the galaxies from GUMS and the magnitude of the cosmological proper motion effect as a function of redshift.

3.2. Simulating quasar observations

To include the Galactocentric acceleration in these data sets, we simulated a pattern of proper motions of the form (Kopeikin & Makarov 2006) (6)

where μ_α^⋆ = μ_αcosδ and is the acceleration divided by the speed of light, Eq. (1), in ICRS. We assumed acceleration components of μas yr^-1 in ICRS, which corresponds to (4.3,0,0)μas yr^-1 in Galactic coordinates, that is, directed towards the Galactic centre.

To simulate quasar photocentric variability, we used a Markov chain with an exponentially decaying correlation (Pasquato et al. 2011; Chatfield 2013; Doob 1942) with a characteristic correlation timescale τ_cor. This results in random variations in positions at time t_i that are both Gaussian and Markovian and generated by (7)where Δt_i = t_i − t_{i − 1} and the and are sampled from a Gaussian distribution with zero mean and standard deviation (8)σ_var is the standard deviation of the random variations Δα_∗(t) and Δδ(t), which are used to perturb the source positions for each quasar independently. We assumed photocentre variations of σ_var = 100 μas and τ_cor = 2 and 10 years (Sect. 2.1). Two examples of Markov chain photo-centre variability are shown in Fig. 2.

We did not simulate the offset of the optical positions of VLBI quasars from their radio positions. This could amount to tens of milliarcseconds or more (Makarov et al. 2012; Orosz & Frey 2013; Zacharias & Zacharias 2014), which could strongly affect the determination of the orientation of the optical reference frame. However, since we are primarily interested in determining proper motion patterns, only the time variation of the offset matters. We assumed that this variation is much smaller and covered by the Markov chain model.

3.3. Simulating galaxy observations

As explained in Sect. 3.1, we simulated only 100 000 galaxies (GAIA_GALAXY) from GUMS, which were used for simulations involving redshift-dependent proper motion patterns. This number is very conservative considering there are 38 million galaxies in the GUMS catalogue, but we also restricted ourselves to a narrow redshift range (0.001−0.03) where the effect is strongest (see Fig. 3). In this redshift range the angular diameter of a typical galaxy core (assuming a core diameter of 1 kpc) is between 46″ and 1.5″, respectively. In the simulation we assumed that these galaxies can be observed with the same accuracies as point-like sources of the same total magnitude. This assumption is certainly unrealistic, given their extended structures. In reality, Gaia may detect multiple compact bright features in a single galaxy, and the centroiding accuracy for the galaxy will depend on the number and brightness of these features. Since the purpose here is to investigate whether these proper motion patterns are potentially detectable by Gaia, we adopted an optimistic approach pending a more detailed study of the centroiding accuracies for such objects. The simulation results must be interpreted keeping this assumption in mind. Clearly, increasing the redshift range will increase the number of objects with small angular diameter, but the amplitude of the proper motion patterns also decreases at higher redshift.

The cosmological proper motion was calculated according to the assumptions in Sect. 2.3. The proper motion components were then obtained by projecting this velocity on the tangent plane defined by the unit vector in the direction of the CMB antapex. In addition, we calculated random proper motion components of the galaxies by assuming a velocity dispersion of v = 750 km s^-1 and again computed the proper motion components from the line-of-sight comoving distance. These different effects were simulated together with the Galactocentric acceleration using Eq. (6), but with the modified values for the coefficients.

3.4. Tools and other simulation parameters

For our simulations we used AGISLab (Holl et al. 2012), which is a simulation tool used to help develop and test concepts and the corresponding algorithms for the astrometric global iterative solution (AGIS; Lindegren et al. 2012). AGISLab implements much of the functionality of the real Gaia data reduction software, but in addition is able to generate simulated observations using realistic estimates of the observation noise as a function of magnitude. Using this tool, sets of true and noisy astrometric parameters can be generated for a range of different types of sources distributed on the celestial sphere. In addition, attitude parameters for Gaia can be determined in the iterative least-squares solution (AGIS) using splines. For both the source and attitude parameters, different initial systematic and random noise values can be added and are then compared to the true values as an iterative solution proceeds until an acceptable level of convergence is achieved. AGISLab generates the observations based on, for example, the satellite CCD geometry and its orbit. Additionally, the direction to a source is computed using the full Gaia relativity model (Klioner 2003, 2004). A set of observation equations are used to construct the least-squares problem for the astrometric parameters using normal equations (Lindegren et al. 2012), which are then solved using a conjugate gradient algorithm described in Bombrun et al. (2012). To the conventional features of AGISLab we added the option to include Markov chain photo-centre variability as a perturbation to the transit times on the Gaia CCDs. Finally, AGISLab contains a number of utilities to generate statistics and graphical output.

In the simulations presented here, we used AGISLab without scaling down (see Holl et al. 2012) the mission and assumed an attitude-modelling knot interval of 30 s. For the starting noisy astrometric parameters we used 100 mas random errors and 10 mas systematic errors, while for the attitude parameters we assumed a nominal attitude error of 10 μas. For the accuracy of the reference frame we compared the final positions and proper motions with different realisations of their true values, assuming a random variation of 100 μas for VLBI_QUASARS in position and 10 μas yr^-1 in proper motion for all quasars. For the simulations we used an antapex direction in Galactic coordinates of (l = 83.99°, b = −48.26°). We assumed a Galactocentric acceleration of 4.3 μas yr^-1, CMB velocity of 369.0 km s^-1, a Hubble constant of 67.3 km s^-1 Mpc^-1, and a five-year mission.

4. Determining the reference frame and proper motion patterns

The relative measurement principle of Gaia results in astrometric parameters (positions and proper motions) that are determined with up to six degrees of freedom in the orientation ϵ = (ϵ_x,ϵ_y,ϵ_z) and the spin ω = (ω_x,ω_y,ω_z) of the Gaia reference frame relative to the ICRS at an adopted frame rotator epoch, t_fr, taken to be the same as the reference epoch of the astrometric parameters. To express the final astrometric results in a celestial reference frame that closely matches the ICRS, the ϵ and ω parameters must be estimated from several sets of sources in a least-squares solution. The determined parameters can then be used to correct the reference frame to coincide with the ICRS. This frame rotation determination is made subsequent to AGIS using three different sets of sources:

• VLBI_QUASAR (called S_P in Lindegren et al. 2012) – a subset consisting of the optical counterparts of a few thousand radio VLBI objects with known positions and proper motions in the ICRS independent of Gaia, which help to constrain ϵ. They are typically the optical counterparts of extragalactic objects with accurate positions from VLBI. These quasars can also be used to calculate the spin parameter (ω), but being very small in number, their contribution to the determination of spin is small.

• GAIA_QUASAR (S_NR) – a larger subset consisting of hundreds of thousands of quasar-like objects (~10⁵–10⁶) taken from ground based and photometric surveys that help to constrain ω. These quasars do not have accurately known positions in ICRS and therefore cannot be used to calculate ϵ. They are assumed to define a kinematically non-rotating celestial frame.

• ICRS_STAR (S_PM) – a subset of primary sources that have positions and/or proper motions that are accurately determined with respect to the ICRS independent of Gaia. This could include radio stars observed by radio VLBI interferometry, or stars whose absolute proper motions have been determined by some other means.

For the present investigation we did not consider the subset ICRS_STAR set, but added another subset of point-like sources at low redshift (z< 0.03):

• GAIA_GALAXY – a subset of point-like galactic nuclei sources with known redshift that may have measurable proper motion, albeit small. This subset can be used to probe redshift-dependent proper motion patters (see Sect. 2.3).

We determined the orientation and spin parameters together with the acceleration parameters, combined into a single parameter array using a least-squares estimation of the positions and proper motions of subsets GAIA_QUASAR and VLBI_QUASAR in the two frames. For GAIA_QUASAR the expression for the apparent proper motion in the Gaia reference frame is given by Eq. (108) in Lindegren et al. (2012), namely

and the corresponding Eq. (110) for VLBI_QUASAR simplifies to

if the position differences are measured at the reference epoch t_ep. In these equations p = [ −sinα, cosα, 0 ] ′ and q = [ −sinδcosα, −sinδsinα, cosδ ] ′.

After the parameters ϵ, ω and were determined, the GAIA_GALAXY subset of low-redshift galaxies can be used to estimate the redshift-dependent parameters² To determine the velocity of the solar system relative to the CMB, we must assume to know the Hubble constant and then solve in a similar way for , assuming ϵ, ω, and H₀ are known. Conversely, if we assume that the velocity of the solar system relative to the CMB is known, then it is possible to determine H₀ by solving for . We note that it is not necessary to solve for the acceleration term in each case, but the acceleration effect is also present when using low-redshift galaxies, and hence it is useful to include it as a consistency check and to obtain the correlation. For these secondary calculations³, we cannot use the quasars because the effect rapidly decreases with increasing redshift and most quasars would contribute very little to the solution (see Sect. 3.3).

5. Results

To determine how well we recover the reference frame, we ran one hundred different realisations of the various simulation cases described below. The results of a comparison between the results and the true values gives the errors from which the mean and standard deviation values are found.

Table 3 gives a summary of the different simulations. Case A is a reference data set where only Gaia nominal observation noise and the nominal attitude noise were added. Using the same noise assumptions, cases B and C then assessed the impact of adding quasar variability with different characteristic timescales. Finally, case D shows the results of using the ground-based quasar list IGQL (Sect. 3) under nominal conditions similar to case A to determine the impact of using a non-uniform sky distribution.

In each of these cases, a first simulation was made using quasars to accurately determine the reference frame parameters, ϵ and ω, and the acceleration of the solar system, . In two independent secondary simulations we then used the values of ϵ and ω found using the quasars to determine either the instantaneous velocity of the solar system relative to the CMB, or the Hubble constant using low-redshift galaxies as described in Sect. 4. As these two parameters are degenerate, we cannot solve for both at the same time and hence we had two secondary simulations. For the secondary galaxy simulations no variability was assumed, and we also calculated the acceleration of the solar system, , which is less accurate because we now used fewer sources. The assumption of using fewer sources may be pessimistic given that the number of objects detected by Gaia is much higher, but a detailed study of the centroiding accuracy for galaxies would be needed for a more realistic simulation of this effect.

Table 4 shows the results of the reference frame determination based on quasars for cases A to D. For the reference case A the mean value for Galactocentric acceleration, which is determined simultaneously with the reference frame, is as yr^-1 with a standard deviation and mean formal error of ~0.26 μas yr^-1 compared to the simulated input value of 4.3 μas yr^-1. This shows that under optimal conditions, the Galactocentric acceleration can be determined to a few percent. The formal errors (⟨ σ ⟩) typically underestimate the standard deviation of the actual errors (Std) because of the source variability and the attitude noise, which are not included in the formal errors. The mean values are consistent with the standard deviations except in case D discussed below. The general trend is as expected when comparing cases A to D, as each case adds more complexity and the standard deviations and the mean formal errors all increase slightly.

The addition of photo-centre variability in cases B and C with characteristic timescales of 2 and 10 years, respectively, does not significantly degrade the solution for Galactocentric acceleration. The mean value for Galactocentric acceleration is very similar at as yr^-1, but both the standard deviation and the mean formal errors of these runs only increase very modestly and roughly agree with, but are slightly better than, the predictions of Perryman et al. (2014). These differences can be partially explained by the use of different numbers of quasars in the two estimates – 170 000 in theirs and 500 000 in ours.

In addition to using the GUMS quasar data, which are well distributed on the sky apart from the Galactic plane, we also used the IGQL data set in case D. The results for acceleration are poorer (for the IGQL we derive a value of as yr^-1 with significant components directed away from the Galactic centre), but the standard deviations and formal errors do not change considerably from Case A. The reason for different mean value is the non-uniform sky distribution of the data set, particularly in the southern hemisphere, where the data set is rather sparse (see Fig. 4) compared to the GUMS catalogue. This results in the spin parameter (ω) being poorly determined for the reference frame, and this consequently affects the determination of the acceleration vector. By comparing the resulting correlation matrices for GUMS and IGQL, we show in Table 5, highlighted in red, that a weak correlation exists between the components of ω and in the GUMS catalogue, but this is much stronger in the IGQL data set. The IGQL will be used in the early Gaia processing and care should be taken when interpreting the reference frame and Galactocentric acceleration results. However, as Gaia detects more quasars itself, the sparse distribution found in the IGQL catalogue will be filled in and the results presented here based on the GUMS catalogue should become more representative. We note that many of the elements of the correlation matrix are zero, which results purely from our choice of reference frame epoch and the source reference epochs, which are all the same. Choosing different reference times would result in non-zero entries in the correlation matrix, but for comparison purposes, the current choice is more useful.

Fig. 4 All-sky maps showing the median number of sources per pixel in an equatorial Hammer-Aitoff projection. Left: quasar distribution on the sky from the GUMS catalogue. Right: quasar distribution on the sky from the IGQL catalogue. The cyan line denotes the ecliptic.

In Table 6 the results from simulations using low-redshift galaxies from the GUMS catalogue are presented. We first determine the individual reference frame parameters ϵ and ω, using the quasar simulations tabulated in Table 4, which are then fixed input for theses secondary simulations. Two independent secondary simulations were made. First, we determined the instantaneous velocity of the solar system with respect to the CMB together with the acceleration. Second, the Hubble constant was solved together with the acceleration. It is not necessary to make these two estimates of the acceleration based on low-redshift galaxies, but doing so allows the correlation with the parameters of interest to be estimated, and it provides a consistency check with the value determined from the quasar data.

If we compare the results for the acceleration determined from the low-redshift galaxies with that determined from the quasars in Table 4, we see that the actual values do not change significantly. However, the standard deviation and mean formal errors have increased by factors of ~1.5 and ~1.8, respectively, for Case A, which is slightly lower than the expected value of due to the difference in the number of objects used. This discrepancy is probably a result of the lack of variability in the galaxy observations and the somewhat higher brightness of the galaxies compared with the quasars (the number of galaxies peaks at ~17 G mag – see Fig. 3 bottom left). Little information on the centroiding accuracy of point-like galaxies (or components thereof) is currently available, and our values here certainly underestimate the real errors. The differences of the acceleration determined from the two secondary simulations are well within the expected variance, and thus we have only reported those from the velocity case. In addition, for case D, the acceleration results are much better. This is because case D in Table 4 used the IGQL quasar data set, with a non-uniform sky distribution, to estimate the values of ϵ, ω, and while in Table 6 we use the uniformally distributed galaxies to estimate the acceleration, which now gives better results. Clearly, the sky distribution is important.

The results for the instantaneous velocity of the solar system with respect to the CMB are also shown in Table 6. The recovered components of the velocity, magnitude, and direction are given. In the simulations we used an apex direction of (l = 263.99°, b = 48.26°) and a magnitude of |v| = 369 km s^-1. The recovered mean values and direction agree well (within five percent) in all simulations. The magnitude of the standard deviation is just over 100 km s^-1, while the magnitude of the mean formal error is of the order of 170 km s^-1. These errors are relatively large, but we recall that we used only 100 000 objects in the redshift range z< 0.03, while the largest contribution to this effect is dominated by a small number of objects in the range z< 0.01 (see Fig. 3 right). The results given in Table 6 for case D are slightly poorer because the fixed reference frame parameters ϵ and ω were less accurately determined using the IGQL quasars. They therefore also have a weak effect on determining the proper motion patterns in these secondary simulations.

For the Hubble constant the recovered mean values also agree well (to within four percent) in all simulations. Again, the errors are relatively large, the standard deviation is 11−12 km s^-1 Mpc^-1, while the mean formal error is of the order of 18 km s^-1 Mpc^-1. Clearly, detecting more sources in the low-redshift range z< 0.01 would greatly improve these results, but it remains unclear for now how many more objects could be used for this determination and how well these extended objects can be centroided. We here used a very optimistic assumption concerning the centroiding accuracy, which means that measuring these proper motion patterns from galaxies will remain a challenging task.

6. Conclusions

We have assessed how well Gaia can improve the optical realisation of ICRS using sets of quasars to define a kinematically stable non-rotating reference frame with the barycentre of the solar system as its origin. Photo-centric variability of the quasars was simulated using a Markov chain model at the level of 100 μas and does not appear to significantly perturb the results for either the reference frame determination or the proper motion pattern that is due to the acceleration of the solar system.

The initial quasar list (IGQL) does not have an ideal sky distribution to determine the reference frame or proper motion patterns, and care should be taken in early Gaia data processing when interpreting theses results. However, as Gaia detects more quasars itself, the sky distribution will improve and eventually become more uniform and representative of the GUMS catalogue used in these simulations. It should be noted that other surveys exist, such as the WISE project (e.g., Secrest et al. 2015), which could improve the situation in the early Gaia solutions.

We reviewed the various proper motion effects that could affect the Gaia data processing, and we tried to quantify their respective magnitudes. Most of the effects are either not detectable by astrometry or do not seem to be strong enough to be reliably detected by Gaia. However, a cosmological proper motion pattern that is due to the instantaneous velocity of the solar system with respect to the CMB might be just about measurable. In contrast to the acceleration of the solar system, this cosmological proper motion pattern is redshift dependent and only significant at low redshift. Thus it cannot be determined from the more distant quasars, and instead low-redshift point-like galaxies must be used. This source type has not yet been considered in the baseline processing for the Gaia mission.

From this proper motion pattern our velocity relative to the CMB can be determined provided we assume that the Hubble constant is known. Alternatively, if we assume our velocity is know from other missions, such as Planck (Planck Collaboration XXVII 2014), we can estimate the Hubble constant from this pattern of proper motions. We find that both measurements might be just within the reach of Gaia, provided that a suitable set of low-redshift point-like objects can be identified and used in the processing. The onboard detection and data processing (centroiding) of such low-redshift point-like objects will ultimately limit whether this effect can be measured in practice.

Acknowledgments

The authors thank the Swedish National Space Board (SNSB) for financial support, without which this project would not have been possible. We also thank the referee, V. Makarov, for useful comments and suggestions on the manuscript.

References

Appendix A: Statistical proper motion of quasars that are due to the random motion of lensing objects

Consider a point source Q and a lensing object L at a distance D_Q and D_L from an observer O (see Fig. A.1). Let ∠QOL = α and ∠Q′OL = β, where Q′ is the shifted image of source Q that is due to the presence of the lens. Then the lens equation is given by (Weinberg 2008) (A.1)where θ_E is the Einstein radius. We assume that the point source is a quasar so that we can assume it to be at rest. Then the proper motion is μ_Q = 0. If the lens has a low velocity v_L along the direction perpendicular to the line of sight, then the proper motion of lens is given by The negative sign is chosen so that the lens motion and the image shift Q′ are in the same direction. If the lens is moving, then the image Q′ will also be moving, resulting in an apparent proper motion From Eq. (A.1)The roots of this equation are given by and differentiating with respect to time gives Assuming that α ≫ θ_E, and taking the positive sign for an image in one direction only, we obtain (A.2)

Fig. A.1 Lensing parameters.

To calculate the combined statistical effects of many lenses at different distances, we first consider a shell at a certain mean distance from us. All the lenses in this shell are assumed to have the same mass and therefore have the same Einstein radius. Let N be the projected density of the lenses in this shell. The mean-squares from Eq. (A.2)is then where the second equality follows from the statistical independence between μ_L and α. With source Q at the centre, we draw a circular ring of radius α and thickness dα. Then the ring will contain 2παN dα lenses and the expectation value of α^-4 is given by Therefore, (A.3)where is the probability that Q will be inside the Einstein ring of some lens in the shell, also known as the photometric optical depth of the shell. Summing the contributions from all the shells, we see that Eq. (A.3)is valid also for the total effect, if τ is the total optical depth. We finally arrive at (A.4)

All Tables

All Figures

	Fig. 1 Variation of comoving distance with redshift.
In the text

	Fig. 2 Markov chain photo-centre variability for a VLBI quasar with timescales, τcor, of 2 (left) and 10 years (right). In each case three examples are given for the same sources.
In the text

Fig. 3 Left: astrometric weight (w) of GUMS galaxy observations (see Hobbs et al. 2009), where , is the total statistical weight of the assumed along-scan standard error of Gaia observations. Below is the number of GUMS galaxies as a function of G magnitude. The two left plots show that the highest astrometric weight is from G = 12 to G = 16, but the number of galaxies peaks around G = 17. Right: redshift distribution of the galaxies from GUMS and the magnitude of the cosmological proper motion effect as a function of redshift.

In the text

	Fig. 4 All-sky maps showing the median number of sources per pixel in an equatorial Hammer-Aitoff projection. Left: quasar distribution on the sky from the GUMS catalogue. Right: quasar distribution on the sky from the IGQL catalogue. The cyan line denotes the ecliptic.
In the text

	Fig. A.1 Lensing parameters.
In the text