The Science of Fundamental Catalogs

This review paper discusses the science of astrometric catalogs, their current applications and future prospects for making progress in fundamental astronomy, astrophysics and gravitational physics. We discuss the concept of fundamental catalogs, their practical realizations, and future prospects. Particular attention is paid to the astrophysical implementations of the catalogs such as the measurement of the Oort constants, the secular aberration and parallax, and asteroseismology. We also consider the use of the fundamental catalogs in gravitational physics for testing general theory of relativity and detection of ultra-long gravitational waves of cosmological origin.


I. INTRODUCTION
Fundamental astronomy is currently an integral part of modern gravitational physics and astrophysics, which rely upon two pillars -astrometry and celestial mechanics. The defining task of fundamental astronomy is to build the inertial celestial reference frame, which is used to determine coordinates, velocities, and accelerations of astronomical bodies and to predict their past, present, and future dynamical evolution in the course of time. Reference objects of fundamental astronomy are stars and quasars, which are benchmarks materializing the inertial reference frame on the sky. Quasars are located at large distances from the Solar System so that their trigonometric parallaxes are difficult or even impossible to measure even now. Historically, the inability to detect the annual parallax was used as an argument against the Copernicus heliocentric theory, but it was based on misapprehension of astronomical distance scales. Parallaxes of just a few dozen stars could be measured until the early 20th century, when specialized astrograph telescopes and photographic plates facilitated a breakthrough in precise astrometry. Celestial coordinates and velocities are traditionally determined as angles on the celestial sphere and their time derivatives (proper motions), in contrast to the Cartesian three-dimensional coordinates and velocities. For a long time, astrometric observations of stars and planets were interpreted in the framework of spherical astronomy and the Newtonian celestial mechanics in the Euclidean space.
The situation had changed about 100 years ago with the advent of a new generation of large optical telescopes and the emergent understanding of the true nature of spiral nebulae and the size of the universe, accompanied by the advent of special and general theory of relativity and flourishing of astrophysics [1]. Fundamental astronomy has undergone dramatic technological changes and explored a broader range of electromagnetic spectrum radiation operating from the ground and from space. A cardinal improvement in the precision of astrometric measurements of positions and parallaxes of stars and quasars has been achieved (in some cases approaching the level of ≃10 µas) by making use of the Very Long Baseline Interferometry (VLBI) [2][3][4] and the Gaia space satellite [5].
One should notice that the measurement of time is an integral part of fundamental astronomy as time is one of the four coordinates of spacetime manifold, which is the arena where all the physical and astronomical phenomena take place. Precise time determination is a cornerstone discipline and technology on its own but its synergy with fundamental astronomy has never been more essential than today. Modern technology allows one to manufacture atomic clocks with a relative fractional instability of the order of a few parts in 10 −18 [6,7]. A global network of such most stable and accurate clocks produces a practical realization of the International Atomic Time (TAI) scale at the level of 10 −16 [8], which is the basic time scale of planetary and lunar ephemerides [9], binary star orbital dynamics [10], exoplanetary search [11,12], pulsar timing astronomy [13], etc. nuclei (AGN) are much further away from the Solar System and their true proper motions are believed to be negligibly small. They are much better celestial reference objects to build an inertial reference frame on the sky. However, galaxies have noticeable angular sizes preventing precise measurement of their photometric centers in the optical wavelengths, which reduces the accuracy of the catalog. This problem is much less essential for quasars, which are luminous nuclei of the most active and distant galaxies. Therefore, the International Astronomical Union (IAU) decided to build the International Celestial Reference Frame (ICRF) by making use of quasars as fundamental reference objects. Millions of quasars have been cataloged using various methods in the optical and near-infrared domain, but only several thousand are bright enough in the radio to be observed by the VLBI, which is a self-sufficient and global technique inherently linked to the Earth orientation parameters and precision timing, thus closely approaching the concept of absoluteness.
The first realization of the ICRF1 was completed in 1998 and contained 212 defining quasars [28]. It was adjusted to FK5, in order to preserve continuity with the series of the stellar fundamental catalogs, and delivered an order of magnitude improvement in the realization of the inertial reference frame on the sky. A decade later, the ICRF2 catalog consisting of 295 reference quasars, was released [29]. It was subsequently replaced by the ICRF3 [30] as adopted by the IAU during its XXX-th General Assembly in 2018 [31]. The ICRF3 contains positions of 4536 quasars as measured at frequencies 8.4 and 2.3 GHz (X/S band). Among them 303 sources are identified as defining sources. These are rather uniformly distributed on the sky and determine the axes of the frame. The quasar positions at 8.4/2.3 GHz are supplemented with positions of 824 sources at 24 GHz (K band), and with position of 678 sources at 32/8.4 GHz (Ka/X band). The positions of ICRF3 sources are adjusted to those in the ICRF2 at the epoch J2015.0 with an accuracy of 30 µas. The advanced feature of the ICRF3 is the accounting for the galactocentric acceleration of the Solar System, which causes a secular aberration effect [32] of 5.8 µas/yr in the apparent proper motions of quasars [33].
The attempt to introduce alternative catalogs obtained by different VLBI configurations at higher electromagnetic frequencies is related to the problem of frequency-dependent astrometric position of interferometric cores located along parsec-scaled jets of AGNs [34] which may reach ∼ 1 mas, thus, compromising the formal precision of ICRF coordinates. The effect of frequencydependent astrometric core positions in ICRF3 measured by the phase and group delay interferometric technique has been studied more thoroughly by Porcas [35], who argued that if the core shift, ∆x, is described by the frequency power-law, ∆x ∼ ν β , with the index β = 1, then, there is no difference between the phase and group delay positions of the core and the absolute radio astrometry gives the position of the AGN jet base without appreciable error. Recently, a more comprehensive investigation of this problem has been conducted [36]. It revealed that the offsets between the core positions is frequency-dependent and may reach nearly 0.5 mas at 2 and 8 GHz. It can also change in the course of time with a characteristic variability of the individual core positions of about 0.3 mas.
The practical realization of ICRF3 heavily relies upon the advanced astronomical and geophysical modeling of VLBI observations [30]. The modeling has been gradually improved over the past 30 years resulting in the latest exhaustive version of the VLBI data processing algorithm documented in the Technical Note 36 of the International Earth Rotation and Reference Systems Service (IERS) [37]. It includes a comprehensive description of the earth's rotation parameters along with its general-relativistic theory of orbital motion in the Solar System, propagation of radio waves in the ionosphere and troposphere, and a self-consistent post-Newtonian theory of reference coordinate systems and gravitational Shapiro's time delay [9,38]. New IERS conventions are currently in preparation and are to be released in 2022.

C. Hipparcos
A great boost to the development of the optical fundamental catalog was given by the launch of the Hipparcos astrometric satellite in 1989 [https://www.cosmos.esa.int/web/hipparcos]. After ∼ 3.5 years of the generally successful mission, the main Hipparcos catalog and several ancillary data sets were published in 1997 [39]. It contains 118,218 entries (single and multiple stars) distributed over the entire sky with a mean density of approximately 3 stars per square degree. The Hipparcos catalog was adjusted to the ICRF1 through a combination of radio star positions and by a few additional indirect techniques. Combining the Hipparcos 5-parameter astrometry with older ground-based observations lead to the creation of the sixth fundamental catalog (FK6) published in two parts: FK6(I) in 1999 [40], and FK6(III) 2000 [41]. The FK6 contains several solutions for the defining stars. The classical one is the single-star mode solution (SI mode) with a mean error in proper motion of 590 µas/yr. The other two solutions are called the long-term prediction (LTP) and short-term prediction (STP) modes. They have been introduced in response to the quasi-instantaneous nature of the Hipparcos measurements which could not discriminate between a binary star with orbital period of a few years and a single star. The LTP and STP modes are the most precise solutions for single Hipparcos stars with a typical mean error of proper motion of 930 µas/yr.
The Hipparcos technology revolutionized the field of astrometry and galactic astronomy [42,43]. It significantly extended the reference frame and achieved a much higher precision of standard astrometric solutions for a single star based on five model parameters -right ascension α, declination δ, parallax ϖ, and two components of proper motion in the plane of the sky, µ α and µ δ [17]. Combining the Hipparcos catalog with the sixth observable parameter -the radial velocity of the starmeasured independently by means of stellar spectroscopy, effectively allowed one to start studying kinematics of stars in the six-dimensional phase space. This provides a complete description of the trajectories of stars in the neighborhood of the Solar System that is vitally important for an adequate understanding of the processes of formation and dynamical evolution of the Milky Way [43].

D. Gaia
The Gaia astrometric satellite [44], launched on 19 December 2013 [https://www.cosmos.esa.int/web/gaia/launch], is the next step forward after the Hipparcos. The astrometric objectives of the mission are to collect and build a 3-dimensional catalog of more than 1 billion stars and 500,000 quasars and to discover thousands of new asteroids and Jupiter-size exoplanets. Like the Hipparcos mission, Gaia measures positions, proper motions, and parallaxes of stars, which (along with precise ground-based and Gaia's own spectroscopic measurements of radial velocities) provide the most comprehensive census of about 1% of the total stellar population of our Galaxy. The Gaia catalog is released in several steps. As of today, two releases of the Gaia catalog -DR1 [https://www.cosmos.esa.int/web/gaia/dr1], [45] and DR2 [https://www.cosmos.esa.int/web/gaia/dr2]have been published.
The Gaia DR1 catalog includes about 1.14 billion sources. Its five-parameter astrometric solution is given for about 2 million stars from the intersection with the Tycho-2 catalog [46] which were treated as single stars without taking their radial velocity into account. The Gaia-CRF1 reference frame is aligned with the ICRF2 radio catalog at the 0.1 mas level at epoch J2015.0, and is non-rotating with respect to the ICRF2 to within 0.03 mas/yr [47]. On the other hand, it was found that the Hipparcos reference frame has a residual rotation with respect to the Gaia-DR1 frame at a rate of 0.24 mas/yr [48].
The DR2 extends the CRF to a larger number of sources amounting to 1.69 billion stars with the mean epoch J2015.5. The DR2 five-parameter astrometric solution significantly extends that of the DR1, and includes about 1.33 billion stars. Moreover, the DR2 provides radial velocities for more than 7.2 million stars with stellar magnitudes in the range 4 ≤ m ≤ 13 and an effective temperature in the range 3550 ≤ T ≤ 6900 K. This yields a full six-parameter astrometric characterization for those stars. The reference frame of the DR2 is a densification of the Gaia-CRF1 from 2191 frame-fixing sources to 556869 sources uniformly distributed across almost the entire sky except along the Galactic equator [49]. The densification of the quasar grid improved the determination and removal of the residual global rotation of the celestial frame caused by the observed proper motion of the frame-fixing sources. The densification does not change the orientation of the Gaia-CRF2 coordinate axes which is still defined by a much smaller number of radio-loud quasars with accurate positions in the ICRF3 catalog.
At the time when this review was in preparation, the Gaia consortium has announced a third data release (EDR3) which includes a few catalogs and documentation, available online at the Gaia website [https://www.cosmos.esa.int/web/gaia/ early-data-release-3]. The contents and survey properties of EDR3 are summarized in [50]. The source list of EDR3 is mostly the same as that for Gaia DR2 but it does feature new additions and there are some notable changes. The creation of the source list for Gaia EDR3 includes enhancements that make it more robust with respect to high proper motion stars and the disturbing effects of spurious and partially resolved sources. Briefly, EDR3 includes astrophotometric solution for approximately 1.8 billion sources brighter than m>21. For 1.5 billion of those sources, parallaxes, proper motions, and the G BP -G RP color are also available. EDR3 also includes 7 million radial velocities from Gaia DR2 after removal of a small number of spurious values. EDR3 represents an updated materialization of the celestial reference system -the Gaia-CRF3 -which is based solely on extragalactic sources in the optical bandwidth and ICRF3 positions. For bright stars (G < 13 mag), however, the Hipparcos catalog was used to correct for a residual spin of the proper motion field. The internal consistency of the EDR3 CRF realization may also suffer from the split of the source list into two categories, viz., 5-parameter solutions, and somewhat less reliable 6-parameter solutions.
Gaia EDR3 represents a significant advance over Gaia DR2, with parallax precision increased by 30 percent, proper motion precision increased by a factor of 2, and astrometric systematic errors suppressed by 30-40 percent for the parallaxes and by a factor 2.5 for the proper motions. The photometry also benefits from higher precision with much better homogeneity across color, magnitude, and celestial position of the sources. Consistent definitions of the G, G BP , and G RP passbands are achieved over the entire magnitude and color range with no systematic above the 1 percent level. The Gaia team works on the data release DR3 which will be based on the improved version of Gaia-EDR3.

E. Other Catalogs and Databases
It is worth mentioning that besides the fundamental catalogs there exists a large number of other astronomical catalogs and databases of various astronomical objects grouped together by a common type, morphology, origin, or method of discovery. Wikipedia lists all existing astronomical catalogs on this website [https://en.wikipedia.org/wiki/List_of_astronomical_catalogues]. Among these, the Tycho-2 catalog of the 2.5 million brightest stars [46] is one of the most useful and frequently used in fundamental astronomy. It contains resolved binary systems with separations down to 0.8 arcsec [51] and proper motions with accuracy up to 2.5 mas/yr. It was amended a few years ago with the data taken from Gaia-DR1 resulting in the new 5-parameter Tycho-Gaia astrometric solution (TGAS), reaching a positional accuracy of stars of 0.3 mas in positions of epoch and 1 mas/yr in proper motions [52].
The VizieR database should also be mentioned, which groups in an homogeneous way thousands of astronomical catalogs gathered for decades by the Centre de Données de Strasbourg (CDS) and participating institutes [53]. This web-accessible database has developed into a powerful tool enabling researchers and users to retrieve and combine astronomical information across various disciplines and domains, including fundamental astrometry.

A. The Oort constants
Astrophysical applications of fundamental catalogs are ubiquitous. Exact positions, trigonometric parallaxes, and proper motions of stars from the catalog along with their radial velocities (when available) form the observational basis to study stellar kinematics of the Milky Way, globular clusters and dwarf spheroidal galaxies from the Local Group. Stellar kinematics of the stars near the Sun allows us to measure the important Oort's constants abbreviated as A, B, C, and K in the linearized twodimensional model of rotation of the galactic disc neglecting its thickness and represented by close orbits of stars. Two of the constants, A and B, which characterize the azimuthal shear and vorticity of the velocity field, were introduced by Oort himself [54] for the idealized axisymmetric model of Milky Way's rotation. The two other constants, currently denoted as C and K, characterize the radial shear and divergence of the velocity field and account for a non-axisymmetric component of the distribution of stellar velocities [55,56] 1 . The four constants define the local rotational properties of the kinematic vector field of stellar velocities and give us insight in the mass density profile of the galactic disc within the solar neighborhood (A and B constants) as well as allows us to evaluate the physically important parameters of the Galactic spiral density wave (C and K constants) [57,58].
Measuring the Oort constants has been a high priority task for fundamental astronomy since 1927 when the Oort's theory of Galactic rotation was formulated [54] and the astrometric measurements of proper motions of stars became sufficiently accurate. Nonetheless, before the Hipparcos mission, the numerical estimates of the Oort constants had a large dispersion depending on the origin of data and methods used to compile different catalogs [59][60][61]. One of the reasons is that prior to the Hipparcos breakthrough, the astrometric measurements of stellar proper motions relied on fundamental catalogs like FK4 and FK5, whose reference systems were tied to the orbits of Earth and major planets and, thus, were sensitive to the precession of the vernal equinox with respect to the ecliptic. Satellite measurements of proper motions, on the other hand, are tied, directly or indirectly, to quasars and distant galaxies and are not sensitive to the possible residual rotation of the reference frame of the fundamental catalog. Additional systematic errors in the values of the Oort constants appeared in previous studies due to the axisymmetric idealization of the Milky Way with almost circular orbits of stars, mathematically represented by a linear tensor. With the growing accuracy and larger samples of stars in astrometric catalogs, these systematic errors became noticeable, so that more accurate models of the motion of the Solar System with respect to the Local Standard of Rest (LSR) had to be worked out along with the necessity to account for the vertical component of the systemic motion of stars and the mode mixing effect [62]. A comprehensive study of the Oort constants using astrometric measurements for 304267 main-sequence stars within a typical distance of 230 pc, taken from the Gaia-DR1 TGAS catalog, was conducted by Bovy [63], who obtained for the local velocity field generated by closed orbits: Will the linearized Oort and Ogorodnikov-Milne models remain adequate in the post-Gaia future? Perhaps, a three-dimensional mathematical representation of the Galactic velocity field with the origin at the Galactic center and a larger number of fitting parameters will become necessary at some point. The vector spherical harmonic formalism is still intrinsically two-dimensional, and therefore, not well suited to capture the increasingly available radial dimension of the field. A galactocentric model based on a vector generalization of solid spherical harmonics, or other combinations of orthogonal vector functions may provide inroads in the formal part of the problem. For a well-conditioned, three-dimensional Galactocentric model, the kinematics of stars beyond the Galactic center has to be accurately known. This goal will require pushing the distance limits of the fundamental catalog to "the other half" of the Galaxy. Interstellar extinction of light within the Galactic disk is the main obstacle, and exploring space astrometry in the infrared, as proposed in the Gaia-NIR successor mission [65], will open up new fascinating possibilities in this area of research.

B. The secular aberration
While the Solar System moves with respect to the LSR with the peculiar velocity u 0 , v 0 , w 0 , the LSR itself moves around the center of the Milky Way along a circular path, which is characterized by the local parameters of Galactic rotation: R 0 -the radial distance to the center of the Milky Way, Θ 0 and (dΘ/dR) R 0 -the linear velocity of the circular orbital motion and its radial gradient, and ω 0 = Θ 0 /R 0 . In the axisymmetric approximation, these parameters are directly related to the Oort constants A and B, that is, The IAU recommended value for the LSR rotational velocity Θ 0 = 220 km/s and the distance to the galactic center R 0 = 8.5 kpc. These values do not directly match the measured values of the Oort constants (1), (2) due to the fact that in the solar neighborhood, non-circular streaming motions are important. Indeed, comparison of (3), (4) with (1), (2) shows that |C| ≃ |K| ≃ 0.2|A| and they cannot be ignored in the determinations of the circular velocity of the LSR from the local kinematics of stars. Since the parameters of Galactic rotation are used for evaluating the total mass distribution inside the Milky Way and testing general relativity with binary pulsars [66,67], their measurement being done independently of the Oort constants plays a crucial role [3,68,69].
It was noticed that the centrifugal acceleration a = dΘ 0 /dt of the circular motion of the LSR is expressed in terms of the product a = ω 0 Θ 0 which can be measured from the analysis of a systematic vector field of the apparent proper motions of distant quasars over the entire sky induced by the the effect of the secular aberration of light. The amplitude of the secular aberration effect was determined at 5.8 µas/yr [33]. It is large enough to be taken into account in the ICRF3 catalog in certain high-accuracy applications [70]. The deviation from this nominal value caused by the peculiar motion of the Solar System with respect to the LSR is less than 0.8 µas/yr [32,33,71,72], which is sufficiently negligible.
The measurement of the secular aberration is a good and independent quality test of the Gaia astrometric solution. The recently announced astrometric solution based on Gaia EDR3 yields the direction of the secular acceleration towards the point on the celestial sphere with the coordinates α = 269.1 deg±5.4 deg, δ = −31.6 deg±4.1 deg, with a value corresponding to a proper motion amplitude of 5.05 ± 0.35 µas·yr −1 [73]. These results are in good agreement with the acceleration vector of the Solar System expected from the current models of the Galactic gravitational potential. It is anticipated that future Gaia data releases will provide the estimates of the secular aberration effect with uncertainties substantially below 0.1 µas·yr −1 .
At the current level of precision of astrometric observations, the secular aberration alone cannot be used for independent measurement of the Galactic rotation parameters. This will be achieved in future astrometric missions (e.g., STARE, NEAT, THEIA [74,75] or Gaia-NIR [76], etc.) that will allow position and parallax measurements to better than 1 µas and proper motions to better than 1 µas/yr. Nonetheless, a recent paper by Bovy [77] demonstrates that measurement of the secular aberration effect in combination with relative accelerations obtained from binary pulsar orbital decays [78] allows one to determine all of the parameters describing the dynamics of our local Galactic environment, including the circular velocity of the Sun, and its derivative, the local angular frequency, the Oort constants, and the Sun's motion with respect to the LSR.

C. The secular parallax
The Solar System moves with respect to the center of the Milky Way, which itself traverses with respect to the cosmological reference frame defined by the isotropy of the Cosmic Microwave Background (CMB) [79]. The overall velocity of the motion of the Solar System barycenter with respect to the CMB is directed towards the galactic coordinates (l, b) = (263.99 • ± 0.14 • , 48.26 • ± 0.03 • ) and has a magnitude of v ⊙ = 369 ± 0.9 km/s [80]. Distant cosmological objects (galaxies, quasars) participate in the Hubble expansion and move fast in the radial direction with respect to the Milky Way. Nonetheless, they also have peculiar velocities with respect to the CMB. The peculiar motion of the Milky Way and distant galaxies is caused by the inhomogeneities of the large scale structure of the universe. Several researchers have recently suggested that catalog's analysis of the proper motions of distant cosmological objects can provide meaningful consistency tests of the standard model of cosmology 3 , place independent constraints on the numerical value of the Hubble constant and the linear growth rate of cosmic structures, and be instrumental in the search for the nature of dark matter [84,85].
As an example, let us consider how the measurement of the vector field of galactic proper motions can constrain the Hubble constant H 0 . The proper motion of each galaxy consists of several components. The first one is caused by the velocity of the Solar System (i.e., observer), v v v ⊙ , with respect to the CMB, which is called the secular parallax effect [86]. This effect manifests as a dipole harmonic in the proper motion field of the galactic population in the direction collinear to v v v ⊙ and amounting to 78 d −1 sin β µas/yr Mpc, where β is the angle between the galaxy and v v v ⊙ . The second component of the proper motion is caused by the peculiar velocity of the galaxy, v v v g , with respect to the CMB. The peculiar velocities v v v g of sufficiently close galaxies are not distributed randomly in space since they are caused by the large scale structure inhomogeneities of the local supercluster [84]. The proper motion caused by v v v g has approximately the same magnitude and distance dependence as the secular parallax proper motion. The third component of the proper motion is caused by the secular aberration. The secular aberration does not depend on the distance. Moreover, it can be measured independently by observing solely the proper motion of very distant quasars [32,71,72] and after that subtracted.
The fact that the peculiar velocities of galaxies are strongly correlated with the large scale structure allows us to compute the best statistical estimates of the corresponding proper motions and subtract them from the relative proper motion to get the secular parallax effect for each galaxy [84]. The evaluation of the secular parallax makes it possible to measure the proper distance d to the galaxy independently of the astronomical "distance ladder" based on standard candles like cepheids, etc. At the same time, a spectroscopic measurement of the redshift z of the galaxy yields the quantity z = H 0 d. Making use of the distance d from the secular parallax, independent constraints on the value of the Hubble constant H 0 can be obtained. This method was used by Paine et al. [85] who estimated a first preliminary limit of 3500 µas/yr Mpc on the secular parallax amplitude using proper motions of 232 nearby galaxies from Gaia-DR2. This limit can be improved by a factor of 10 at the end of Gaia mission to set a tighter constraint on the Hubble parameter.

D. Asteroseismology and parallax zero-point
A completely new application of the data from fundamental catalogs emerged from their growing capability to deliver a consistency check of the theory of the internal structure of stars. Conventionally, it was broad-band photometry and trigonometric parallax, ϖ, of stars that have been used to determine their physical characteristics: the effective surface temperature T * , radius R * , and luminosity L * , under the condition that extinction A λ and reddening have been measured as well [64]. The photometric parallax of a star is given by where m λ and BC λ are the star's apparent magnitude and bolometric correction at a given wavelength, M bol⊙ = +4.75 is the Sun's absolute magnitude, and R bb is the radius of the black body having the effective temperature T * . Notice that, generally speaking, R bb is close but not equal to R * . Asteroseismology provides the theoretical relation between the star's radius and the global seismic parameters [87] where ∆ν is the average frequency separation between the modes of stellar oscillations and ν max is the frequency corresponding to the maximum observed oscillation power, both normalized to the solar values of the helioseismological parameters [88]. Making use of (7) in (5) yields an asteroseismic parallax where A ≃ 1 is the scaling factor taking into account the deviation of R bb from R * and normalizing the global asterosesmological parameters ∆ν and ν max to the solar values ∆ν ⊙ and ν max⊙ [87,88]. The astreroseismic parallax ϖ sei can be independently computed from corresponding observations of stars and compared with the trigonometric parallax π of the same stars from a fundamental astrometric catalog like Hipparcos or Gaia-DR2. This way, the scaling factor A can be determined along with the so-called global parallax zero-point π 0 of the astrometric catalog [89][90][91]. The value of A ≃ 1 is used to calibrate the theory of stellar oscillations while the knowledge of π 0 ≃ 0 provides the correct value of trigonometric parallaxes. The parallax zero-point can be also determined independently by making use of distant quasars which have negligible trigonometric parallaxes for z > 0.1, but a difficulty arises with possible dependence of the zero-point on magnitude, which affects a large number of brighter stars.
A comparison of the different types of parallaxes has been recently conducted by Khan et al. [87] who used the Kepler mission data [https://www.nasa.gov/mission_pages/kepler/main/index.html] to compare the results with Gaia-DR2. They find that there is no absolute standard within asteroseismology, because different seismic approaches to the parallax determination problem extensively produce a parallax zero-point, which is fairly different from the Gaia-DR2 π 0 ≃ −29 µas determined by Lindegren et al. [92] by making use of 555934 sources identified as quasars. Similar analysis was done by Zinn et al. [93] who used asteroseismic data of evolved red giant branch stars from the APOKASC-2 catalog to determine the zero-point offset of Gaia DR2 parallaxes. Testing asteroseismology with Gaia DR2 has been also undertaken by Hall et al. [94] who used a sample of 5576 red clump stars from the Kepler mission field.
Broadly speaking, the parallax zero-point π 0 is the systematic parallax error which applies to all objects in a fundamental catalog and, therefore, profoundly affects many important applications, including the fundamental distance scale. For Gaia, it comes from specific thermally-induced variations of the basic angle with time, which are periodic functions of the Sun's aspect angle [91]. The Gaia basic angle monitoring system (BAM) is not sufficient to track this variation at the required level of accuracy, resulting in a complicated instrument calibration scheme. It is not easy to estimate the parallax zero point because the error may be both magnitude and color-dependent through the functional dependencies of the calibration models and intrinsic poor conditioning. Some estimates of the parallax zero-point in EDR3 catalog are presented by Lindegren et al. [95] along with a complete list of references to previous papers and in-depth discussion of the problem. Notably, compared with Gaia-DR2, the astrometric calibration models in EDR3 have been extended, and the spin-related distortion model includes a self-consistent determination of basic-angle variations, improving the global parallax zero point. The global parallax zero point of EDR3 is about -17 µas [96]. The parallax zero point issue will be hopefully resolved, at the level of a few µas, in future data releases that will rely upon further improved calibration models.

A. Testing general relativity
The most precise test of general relativity before the advent of gravitational wave astronomy was achieved with timing of compact binary pulsars [66,97]. The compact binary pulsars have short orbital periods and move fairly fast (v/c ∼ 10 −3 ), which makes them ideal celestial laboratories for measuring a number of various relativistic parameters. Most of the parameters are associated with the relativistic effects in the conservative (1-st and 2-nd) post-Newtonian approximations. Gravitational waves carry out the orbital energy and angular momentum of the source binary system causing the orbital period, P b , to decay. The rate of the decay was computed in general relativity (GR) with different mathematical techniques (see review by Damour [98] for historical details) and is given by [99][100][101][102][103] where m 1 and m 2 are the masses of the pulsar and its companion, e is the eccentricity of the binary system, G is the universal gravitational constant, and c is the fundamental speed of gravity (equal to the speed of light in vacuum). The observed value ofṖ obs b is contaminated by a number of relativistic effects caused by the galactic motion of the binary pulsar with respect to the Solar System [104] where the numerical value ofṖ gal b crucially depends on our knowledge of the stellar kinematics of the disk of the Milky Way. According to Damour and Taylor [104] where Θ 0 is the circular velocity of the LSR, l is the galactic longitude of the binary pulsar, δ = d/R 0 , d is the distance between the binary pulsar and the Sun, R 0 is the distance from the LSR to the center of the Milky Way, A and B are the Oort constants. Testing general relativity in binary pulsars is based on measurement of the five classic (Keplerian) parameters of the orbit along with any two (post-Keplerian) relativistic parameters of the 1-st post-Newtonian approximation of general relativity [105]. Usually, the two post-Keplerian parameters are -the secular drift of periastron <ω >= fω(m 1 , m 2 , e, P b ), and the Dopp-ler+Einstein gravitational delay γ = f γ (m 1 , m 2 , e, P b ) -which depend on the masses of the pulsar and its companion, orbital eccentricity and the orbital period. The eccentricity e and the orbital period P b are known from the measured Keplerian parameters. Hence, measuring <ω > and γ allows to determine the masses m 1 and m 2 . These values can be substituted to the right hand-side of (9) along with P b and e yielding Ṗ b /P b GR pred , which can be compared with the observed value Ṗ b /P b GR .
Unfortunately, the direct comparison is impossible, as one can see from (10), because of the galactic contribution Ṗ b /P b gal to the observed orbital decay rate. The galactic contribution is a complicated function (11) of the kinematic parameters of the Solar System and the binary pulsar, which are available from fundamental catalogs. Current determinations of the kinematic parameters entering the right hand-side of (11) are so reliable that the galactic contribution can be accurately computed and removed from the right hand-side of (10). The remaining deviation of the predicted value of the orbital decay rate from its observed value is i.e., the theory is validated within 0.16% [67,106]. It is remarkable that this type of testing of general relativity can be inverted. As recently emphasized by Chakrabarti et al. [78], assuming that general relativity holds and, therefore, the gravitational-wave orbital decay contribution can be predicted and accurately subtracted, the observed orbital period changes can be converted into measurements of the relative Galactic acceleration at the binary pulsar's and Sun's locations. Chakrabarti et al. [78] used such relative accelerations derived from observations of 14 binary pulsars to constrain the local gravitational field produced by the Galactic potential. They found the Oort limit, that is the total volume mass density in the Galactic mid-plane, to be equal to 0.08 0.05 −0.02 M ⊙ pc −3 . Accounting for the baryonic mass budget derived in [107], Chakrabarti et al. [78] obtained a local dark matter density equals to ρ DM = −0.004 0.05 −0.02 M ⊙ pc −3 , which is lower than that obtained by McKee et al. [107] but consistent within the current standard estimates of its measurement in the solar neighborhood by independent methods [108].

B. Gravitational waves
Gravitational waves are the solutions of Einstein's equations describing the propagation of ripple perturbations of spacetime curvature with the speed of light. In fact, the speed of gravitational waves is set to be equal to the speed of light by the Einstein's postulate that the fundamental speed of gravity c g appearing in the gravity sector of general relativity 4 has the same numerical value as the speed of light c from the electromagnetic/matter sector [111,112]. This postulate has been recently confirmed experimentally by measuring the bending of light by the moving gravitational field of Jupiter [113,114] and in observations of gravitational waves [115], in complete agreement with general relativity.
Gravitational waves were detected in 2015 by the Laser Interferometry Gravitational-wave Observatory (LIGO) [116,117]. This discovery opened a new window to study the nature of black holes, neutron stars and astrophysical processes in the very early universe, which are not accessible for observations by any other means [118]. The gravitational-wave spectrum is distributed approximately over 10 decades of frequency, from the high-frequency normal modes of oscillating neutron stars down to the lowest frequencies produced by various cosmological mechanisms such as inflation, phase transitions, cosmic strings, etc. [119,120]. Plane gravitational waves produce periodic fluctuations in the apparent positions of distant astronomical sources and in the time of propagation of electromagnetic signals from them. These fluctuations have a characteristic amplitude that is proportional to the strain of the gravitational wave [121,122]. Direct astrometric measurements of gravitational waves from a single source are not possible with the current pulsar timing or VLBI technology. However, it was noticed that the stochastic ensemble of these waves produced by various mechanisms can be detected in a roundabout way by studying temporal and angular correlations in the times of arrival (TOA) of radio pulses from the pulsar timing array (PTA) [123][124][125][126][127], or in the redshifts of quasars [128], or in the pattern of their proper motions in an astrometric fundamental catalog [129][130][131][132][133][134][135][136].
The temporal two-point cross-correlation function in the times of arrival of radio signals is known as the Hellings-Downs curve [123,137] ζ HD ab = sin 2 θ ab 2 3 ln sin where h c ( f ) is the characteristic strain of the gravitational wave metric perturbation at frequency f , θ ab is the angle between two sources of radio waves labeled a and b respectively, δ ab is the Kronecker symbol. The magnitude of the strain at frequency f is related to the dimensional energy density of gravitational waves Ω g ( f ) at this frequency by where H 0 is the present value of the Hubble constant. The shape of the Helling-Downs curve (14) provides a means to measure polarization properties of gravitational waves and to detect possible violations of general relativity [137,138]. Measuring the amplitude of the temporal cross-correlation function (13) allows us to evaluate the overall energy density of the ensemble of gravitational waves as well as to get information about their spectrum. The PTA technique is the most sensitive to gravitational waves in the frequency band ranging from 1 nanohertz to a few tens of microhertz, where the mergers of supermassive black hole binaries (SMBHBs) are considered to be the most promising sources of gravitational waves. Preliminary results on the spectrum, energy and other characteristics of gravitational waves emitted by SMBHBs have been recently obtained from analysis of the 11-year data release of monitoring of 45 millisecond pulsars by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [139] and published in [140][141][142].
The angular two-point cross-correlation function has been derived by Pyne et al. [129], Gwinn et al. [130], who focused mainly on the contribution of the second spherical harmonic in the expansion of the proper motion of light sources caused by stochastic gravitational waves over the entire sky. Recently, these calculations have been extended to all spherical harmonics by Book and Flanagan [134] and Klioner [143] and we summarize them below for the reader's convenience. Let the source of light be observed in the direction defined by a unit vector n n n = (n i ) = (cos α cos δ, sin α cos δ, sin δ) where the angle α is the right ascension and δ is the declination. By definition, the proper motion is a time derivative of the unit vector, µ µ µ =ṅ n n which is a vector field on a unit sphere. The angular cross-correlation function of proper motions of distant sources caused by a stochastic gravitational wave ensemble, is where the angular brackets denote the averaging over the ensemble of stochastic gravitational waves, the (geometric) vector indices i and j take values {1, 2, 3}, the labels a and b numerate the sources from the astrometric catalog seen in the directions of unit vectors n n n a = (n i a ) and n n n b = (n i b ) respectively, and µ i a =ṅ i a , µ i b =ṅ i b . The Book-Flanagan cross-correlation function reads [134] ζ BF ab = 1 2 − 6 sin 2 θ ab 2 tan 2 θ ab 2 ln sin θ ab 2 where θ ab is the angle between the two sources, that is, cos θ ab = n n n a · n n n b , and three mutually-orthogonal vectors A A A = sin −1 θ ab (n n n a × n n n b ) , B B B = n n n a × A A A , The vector field of proper motions caused by stochastic gravitational waves can be decomposed with respect to the basis of vector spherical harmonics where Y Y Y E lm (n n n) and Y Y Y B lm (n n n) are the electric-and magnetic-type transverse vector spherical harmonics [144]. The cross-correlation matrix of the angular spectral decomposition is where the labels Q ∈ {E, B} and Q ′ ∈ {E, B}. The coefficients [134] C QQ ′ lml ′ m ′ = d 2 Ω n n n a d 2 Ω n n n b Y * lm (n n n a )Y l ′ m ′ (n n n where Y lm (n n n a ) are the scalar spherical harmonics, the integration is over the unit sphere, ∇ a i and ∇ b i are the normalized 3dimensional gradients with respect to x i a and x i b entering the definition of the unit vectors: n i a = x i a /|x x x a | and n i b = x i b /|x x x b |. Computation by (22) yields [134] where Statistical power of the expectation value of the proper motion induced by the stochastic gravitational waves is where and the numerical coefficients α l characterize the power of the l-th spherical harmonic in the proper motion spectrum. Gwinn et al. [130] calculated α 2 = 5/6. Book and Flanagan [134] provided numerical estimates for next 10 coefficients up to l = 11. Klioner [143] derived a general mathematical expression for all harmonics, The expansion coefficients α l decrease asymptotically as l −5 as the harmonic degree l increases. The formalism of the astrometric detection of ultra-long gravitational waves described above can be used in future astrometric missions having much better precision than Gaia [145]. The current level of astrometric technology does not allow us to set useful physical constraints of the density of primordial gravitational waves in the early universe [122], although such an attempt has been done with radio-quasars [129,130]. The main reasons for the insufficient sensitivity of astrometric detection of gravitational waves are relatively large accidental and systematic errors of the coefficients of vector spherical harmonics representing the observed proper motion field of CRF-defining quasars and AGNs as well as the occurrence of post-fit residuals with perturbations (outliers) beyond the statistical expectation. The latter emerges in comparison analyses of Gaia and ICRF astrometry [146], and may partly be caused by the morphology of the optical counterparts. Users of the radio-optical CRF should also remember that the global spin of the Gaia celestial coordinate system has been technically adjusted to zero on the specially selected sample of reference objects.

V. FUTURE OF FUNDAMENTAL CATALOGS
There are two basic directions in the progressive development of astrometric catalogs establishing the inertial reference frame for various applications in astronomy and fundamental physics. They are concerned with the continuous improvement of astrometric precision in positions, proper motions, and parallaxes, and the increase in the number of catalog's objects. Currently, the astrometric precision of fundamental radio and optical catalogs is about 100 µas [30,31,49]. It was recognized long ago that at this level of precision, the data processing of observations used for the construction of catalogs would require an elaborate general-relativistic model of celestial and terrestrial reference frames [9,147,148] and propagation of light [149,150], along with a substantially accurate description of motion of the Solar System bodies [151][152][153][154] and the observer (Earth rotation parameters) [155][156][157]. These relativistic models have been worked out by a number of research groups around the globe and summarized in the form of standards of the International Earth Rotation and Reference Systems Service (IERS) [37].
The near-future goal of fundamental astronomy is to attain the precision of 1 µas [3,[158][159][160][161] and to progress beyond this threshold as soon as technology permits [162][163][164]. In a more distant future, it may be possible to achieve a sub-microarcsecond astrometric resolution [165]. Data processing of sub-microarcsecond astrometric observations will require taking into account a great deal of relativistic effects both in the motion of the observer and the source of light as well as in the propagation of light rays [145]. Relativistic models of equations of motion of the Solar System will have to be extended, at least, to the level of the 2-nd post-Newtonian approximation and to incorporate a large number of gravitational relativistic multipoles of the Solar System's bodies [166][167][168] while the propagation of light must include various relativistic effects in the time delay and deflection angle caused by the higher-order post-Newtonian terms [169,170] along with the orbital and rotational motion of the light-ray deflecting bodies [171][172][173] and their complicated multipolar structure [174][175][176].
Beyond the sub-microarcsecond threshold, one will see a new population of celestial physical phenomena caused by the presence of primordial gravitational waves from the early universe and different localized astronomical sources like binary stars, space-time topological defects (cosmic strings), moving gravitational lenses, the time variability of gravitational fields of the Solar System and the sources of light, and many others -see the review by Vallenari [15] for more detail. Adequate physical interpretation of these yet undetectable sub-microarcsecond phenomena cannot be achieved within the currently applied models, which will require significant development and extension into the field of relativistic gravitational physics [177] and comprehensive understanding of the astrophysical origin of the astrometric jitter in position, proper motion, and structure of the observed light sources [15].
A practical realization of sub-microarcsecond fundamental catalogs can be achieved with a number of different astronomical techniques in various ranges of the electromagnetic spectrum. In particular, radio astrometry at this level of accuracy adhering to the time-tested concept of absolute measurements requires technical facilities like VLBA [https://science.nrao.edu/ facilities/vlba/docs/manuals/oss], VERA [178] and the Square Kilometer Array (SKA) [3,158,179,180] with baselines ranging from the Earth's radius to the size of the lunar orbit. First successful steps in this direction have been recently taken with the space radio-interferometric mission RadioAstron [181][182][183] which proved to be very productive [http://www.asc.rssi. ru/radioastron/publications/publ.html].
Looking into the future of precision optical astrometry, the prospects seem to be less well defined, because it is not obvious which technological innovation would bring a significant jump in accuracy and productivity compared to the currently operating Gaia mission. The Space Interferometry Mission (SIM) was a bold initiative to move beyond 1 µas in differential astrometry by replacing a traditional telescope with an imaging camera with a Michelson-type optical interferometer [184]. The theoretical astrometric precision of phase referenced interferometers is defined by the ratio λ/B, where λ is the effective wavelength and B is the baseline length, for photon shot noise limited, monochromatic interferometric fringes. The corresponding performance limit of a traditional imaging telescope is proportional to λ/D, with D being the entrance pupil diameter. Thus, interferometry can in principle provide much better astrometric precision by using baselines exceeding the largest feasible collecting areas of space telescopes. SIM, however, was much limited by budgetary constraints and the basic technical design of a single-body space craft.
The astrometric capabilities of optical interferometry could be indefinitely improved by operating a set of siderostats in a formation-flying mode with a separate beam combiner space craft, similar to the now defunct DARWIN concept [185], but without the nulling function for direct imaging. The technical challenges of monitoring and maintaining long baselines between free-floating telescopes within microns are formidable, however. In the domain of proven technology, the Theia initiative [75] promises to surpass the 1 µas threshold for brighter targets, but only in the narrow-angle regime of differential astrometry, which can hardly contribute to fundamental astrometry. The proposed Gaia-NIR project [65], a successor to Gaia in the near infrared, will be truly a global astrometry mission greatly improving the Celestial Reference Frame by obtaining a second epoch of reference object positions, observing deeper into the dusty Galactic belt and beyond, and measuring fainter infrared quasars and AGNs, but only a moderate progress in the single measurement accuracy is anticipated.