© The Authors 2026

1 Introduction

At the turn of the 21st century, the fundamental frame of reference used to define astronomical positions and the orientation of the Earth was created using bright, compact, extragalactic radio sources monitored with very long baseline interferometry (VLBI; Arias et al. 1995). The geodetic VLBI technique is especially suitable for realizing a celestial reference frame because the positions of radio sources on the sky and the positions of the radio telescopes on Earth are solved for simultaneously as a global solution that directly links the celestial reference frame to the terrestrial reference frame, yielding the Earth orientation parameters.

Following the first International Celestial Reference Frame (ICRF1; Ma et al. 1998), the ICRF has undergone a process of gradual expansion and refinement, initially being defined solely using the geodetic 2.3 and 8.4 GHz bands (respectively S/X), but later also defined at 24 GHz (K) and 32 GHz (Ka) with the release of the third and current ICRF (ICRF3; Charlot et al. 2020). In comparison to ICRF1 and ICRF2 (Fey et al. 2015), ICRF3 has fewer distortions in the frame because the defining sources were selected to be distributed more uniformly across the celestial sphere. Additionally, ICRF3 was corrected for Galactic aberration, with Charlot et al. (2020) adopting a value of 5.8 μas yr⁻¹ (see also MacMillan et al. 2019). The rotational stability of the ICRF is typically estimated to be a few tens of μas, with 30 μas estimated for ICRF3 S/X (Charlot et al. 2020).

This error, however, is not due to measurement noise from a finite number of defining sources. As we emphasize here, most of the frame instability is systematic and likely to be mainly caused by effects inherent to the astrophysical properties of the defining sources. One such effect that is commonly discussed in the literature is the frequency-dependent core shift (Kovalev et al. 2008), which causes astrometric offsets of ~0.1–1 mas on average in phase-delay astrometric measurements, and offsets in group delay-based positions such as those that were used to construct the ICRF when the jet is not in equipartition (see, e.g., the discussion in Porcas 2009). The majority of sources that have extended structure at VLBI scales show evidence for core shift variability (Plavin et al. 2019), and the deviations from the core distance-frequency dependence r(ν) ∝ ν⁻¹ equipartition prediction (Blandford & Königl 1979) regularly exceed ~0.3 mas (see Sokolovsky et al. 2011 and corresponding discussion in Plavin et al. 2019). As noted by Petrov & Kovalev (2025), particle flares in the jet, which cause core shift variability, disrupt equipartition. As these flares are a universal feature of radio jets, group delays may not in general be considered impervious to core shift. Cigan et al. (2024) recently showed based on S/X position time series that ICRF3 sources exhibit a typical excess astrometric noise of ~0.1–1 mas (their Fig. 8), and half of the defining sources showed excess noise greater than ~0.2 mas (their Table 1).

The selection of defining sources for a fundamental reference frame such as the ICRF therefore requires care to select the most stable and well-behaved objects. In ICRF1 and ICRF2, this was done by selecting sources with observing sessions that covered at least 2 yr, with cuts on formal VLBI position uncertainties made for ICRF1 and position stability and source structure made for ICRF2. In ICRF3, special emphasis was placed on selecting defining sources based on uniformity of sky distribution. Position stability over several epochs, however, is no guarantee that a source will remain astrometrically well behaved. An extraordinary example of this is 3C48, which, despite a sub-milliarcsecond astrometric position uncertainty in ICRF3 for 18 observing sessions preceding 2017 September 7, exhibited a shift of 57 mas in declination shortly thereafter¹. Astrometric shifts of similar magnitude are known for at least three other sources (Titov et al. 2022), which raises the disturbing prospect that an apparent astrometric stability is of limited benefit when choosing sources for a reference frame. None of the four sources studied by Titov et al. (2022) is defining, but Cigan et al. (2024) recently identified several ICRF3 defining sources with excess astrometric dispersions of ~1 mas (see their Table 2).

It is clear that in order to establish a long-term stable ICRF, it is desirable to determine not just which sources, but which type of sources are astrometrically better behaved across epochs. This is especially important now that the ICRF is multifrequency, with Gaia-CRF3 (Gaia Collaboration 2022) being the official instantiation of the ICRS system at visual (optical) wavelengths². Frequency-dependent source physics, such as the aforementioned core shift, limit efforts to create a consistent multiwavelength reference frame³. Understanding the processes that create frequency-dependent positions, such as synchrotron opacity, contamination by the host galaxy, or offsets between the optically bright accretion disk and the radio core should therefore be at the forefront of reference frame research.

Toward this goal, Secrest (2022) showed that selection on photometric variability has a profound effect on the prevalence of offsets between the Gaia optical and ICRF radio positions, which affect about 10%–15% of the ICRF3 sources overall. While about 20%–25% of the least photometrically variable ICRF sources exhibit optical-radio offsets, this prevalence drops to just a few percent in objects with a photometric variability exceeding ~20% (see Fig. 2 in Secrest 2022). Objects with large fractional photometric variabilities are likely blazars: sources whose relativistic jet is pointed close to our line of sight. Through a comparison of optical color and γ-ray counterpart incidence, Secrest (2022) showed that this is unambiguously the case for ICRF objects. This result was further strengthened by an analysis of geodetic VLBI radio position variability published by Lambert & Secrest (2024), who showed that VLBI position stability is indeed correlated with photometric variability, with blazars such as flat-spectrum radio quasars (FSRQs) and BL Lacs having the most stable positions (see also Lambert et al. 2025).

These results are naturally interpreted as arising from two factors inherent to blazars. First, the direction of the jet near to the line of sight minimizes projected position offsets, so the apparent positions of blazars are largely unaffected by the structural variations of their jets. Second, the small viewing angle strongly boosts the apparent luminosity of the jet at all frequencies, so the contributions of nonjet structures such as the accretion disk are proportionally smaller, which further improves the positional consistency. A key prediction made by Secrest (2022) is that, as long as the direction of the jet remains close to the line of sight, blazars are positionally stable, which makes them ideal objects for a stable reference frame in the long term. We discuss the merits of the assumption of long-term blazar stability in Sect. 4.1 and show that it is likely justified for the purpose of creating a stable reference frame.

In this work, we quantify for the first time the stability of the celestial reference frame as a function of astrophysical source type and show that selection on blazars dramatically improves the stability of the frame. In Sect. 2 we detail our formalism and methods. Our results are presented in Sect. 3, followed by a discussion in Sect. 4. Our main conclusions are listed in Sect. 5.

2 Methods

2.1 Formalism

We adopted the l ≤ 2 vector spherical harmonic (VSH) model to parameterize the spin (l = 0), dipolar glide deformations (l = 1), and l = 2 quadrupolar distortions of the reference frame relative to a fiducial standard (see Mignard & Klioner 2012). The VSH equations listed by Titov & Lambert (2013, see also Charlot et al. 2020) were solved using generalized least squares, $$\hat{\mathit{\epsilon }}={\left({\mathit{A}}^T{\mathbf{\Sigma }}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\mathbf{\Sigma }}^{-1}\mathbf{\Delta },$$(1)

where A is a 2n × 16 array of design matrices for n sources, Σ is the 2n × 2n block diagonal matrix of position covariances (the covariance terms between different sources were set to zero), and Δ is a 2n × 1 column matrix of the position offsets in right ascension and declination. The formal covariance of Eq. (1) is $${\mathbf{\Sigma }}_{\hat{\mathit{\epsilon }}}={\left({\mathit{A}}^T{\mathbf{\Sigma }}^{-1}\mathit{A}\right)}^{-1}.$$(2)

Equation (2) assumes that the data covariances are both accurate and Gaussian. This is not usually true of VLBI data, for which the practice has been to apply a scale factor of 1.5 to the formal errors, followed by the addition in quadrature of some noise floor (in the case of ICRF3 S/X, 30 μas; Charlot et al. 2020). The noise floor was derived using decimation tests in which the data were split by epoch to estimate positional variability. For the purposes of this work, however, these post hoc uniformly applied corrections serve only to dilute the astrometric quality differences between source types, as we show below.

2.2 The bootstrap

In order to quantify the stability of the celestial reference frame as a function of astrophysical source type, we used random sampling with replacement (the bootstrap; Efron 1979). The bootstrap is a means of estimating the sampling distribution of an estimator given a particular set of data, so it therefore allows estimation of the frame stability performance of a source class as a population. Bootstrapping also has another advantage in that it provides a means of estimating the effect of source variability within a class. This is because individual sources are uncorrelated, so the state of one source at the time of a particular observation may be taken as representative of the state of a different source at the time of some future observation that occurs beyond the timescale on which the observed position of a source is correlated with its previously observed position. Timeseries analysis of the VLBI positions of ICRF objects suggests that this timescale is months to years (Cigan et al. 2024). Given that there are physical processes in AGN that cause apparent position variation with time and wavelength, such as variations in synchrotron self-absorption, variations in Doppler factor, flares in the shock structure of the jet, superluminal motion, or accretion disk variability, and given that these physical processes are uncorrelated between sources, so different sources are in different states at different times (e.g. a flare in one source is unrelated to a flare in another source), it follows that comparison of different but randomly selected sets of sources is logically equivalent to comparison of one set of sources at a different epoch. The variation in frame orientation between samples beyond the dispersion expected from measurement or statistical error is directly due to source positional variability. This argument suggests, however, that frame stability is properly tested through random sampling without replacement. Unfortunately, there are not enough sources to create the subsamples for this test. For example, there are 180 blazar-like objects in the sample of 1938 ICRF3-Gaia sources developed by Secrest (2022), so the subsamples created by random sampling without replacement just 10 times would be dominated by statistical error.

To assess the validity of allowing replacement, we drew 180 random sources from the objects studied by Secrest (2022). We generated a large number of simulated epochs and assigned to each epoch a randomly oriented offset drawn from the empirical distribution of statistically significant offsets (having covariance-normalized values >5; see Secrest 2022) observed between the ICRF3 S/X and Gaia DR3 positions of our sample (Sect. 2.3). To reflect the prevalence of significant optical-radio offsets seen in real data, we assigned offsets to 20% of the sources, which varied randomly between epochs. We then permuted the positions of the simulated offsets by their formal measurement covariances from Gaia and ICRF3. For a random subset of epochs, we made bootstrap resamples of the sources and fit each resample with the l ≤ 2 VSH model discussed in the previous section. In Fig. 1, we compare the variance of the VSH components between bootstrap resamples to the true variance across all epochs. The bootstrap resampling of a set of sources at a particular epoch gives an estimate of the multi-epoch stability of the rotation, glide, and quadrupole parameters that is close to truth on average. It therefore provides a means of assessing the frame stability.

2.3 Sample

The 1938 ICRF3-Gaia matches used by Secrest (2022) were selected from an initial match between ICRF3 S/X and Gaia-CRF3 (Gaia Collaboration 2022) to within 100 mas (3142 objects), required to have measured Gaia G-band fractional photometric variability (fractional_variability_g, hereafter F_var; 2426 objects), no evidence of source extent or evidence for close secondary sources (1953 objects), and no evidence for spurious parallaxes or proper motions (1938 objects), which can be further evidence for multisources (Souchay et al. 2022; Makarov & Secrest 2022). As stated previously, 180 have F_var > 0.4, which Secrest (2022) showed are likely pure line-of-sight blazars, 1455 have F_var < 0.2 and have properties consistent with radio-quiet quasars, so they likely do not have jets pointed as close to the line of sight, and 303 have 0.2 < F_var < 0.4, which are likely objects intermediate between quasars and pure blazars in which the jet is significantly beamed, but not pointed directly at the observer.

These three fractional variability classes of sources, which we refer to as blazars, quasars, and intermediate sources, formed the basis of the three samples we used. In order to facilitate a meaningful comparison of the frame stability as a function of source class, however, we required that all three samples have the same size to control for purely statistical error and to be as closely matched on the sky as possible to control for the susceptibility to particular distortions. The high F_var blazar-like category, which has the fewest objects, set the sample size n, and we therefore selected the n closest quasar-like and intermediate sources to the blazars. For sources with S/X, K, and X/Ka positions, n = 180, n = 75, and n = 68, respectively. To facilitate consideration of a reference frame with defining sources selected from bona fide variability-selected blazars, we also produced a matched sample of the ICRF3 defining sources.

Position offsets were taken between Gaia and the ICRF3 source positions. At the epoch of Gaia DR3 (2016.0), Gaia positions of unresolved sources generally have reliable position covariances, with uncertainties likely underestimated by ~6% (see discussion in Secrest 2022). Since the ICRF3 and Gaia data we used differ in epoch by only 0.5 yr, the effect of Galactic aberration, which creates a glide of about 5 μas yr⁻¹, was not considered. Because the post hoc corrections to the VLBI position uncertainties applied to ICRF3 (Charlot et al. 2020, Table 6) only dilute the observable (the presence of systematic position error due to astrophysical effects), we restored the original VLBI position uncertainties.

Finally, we note that while high fractional variability reliably selects for bona fide blazars, this metric is probably substantially incomplete because many genuine blazars did not undergo flaring during the epochs in which Gaia observed them. Consequently, many genuine blazars undoubtedly have F_var values below the cut employed by Secrest (2022) and Lambert & Secrest (2024). Employing other criteria such as cuts made on the qso_variability parameter (Butler & Bloom 2011) provided in Gaia DR3, inclusion of additional data such as multi-epoch photometry from surveys such as the Zwicky Transient Facility (Bellm et al. 2019) or the Near-Earth Object WISE (NEOWISE) mission (Mainzer et al. 2011, 2014), or consideration of γ-ray detection by the Fermi Large Area Telescope (Abdollahi et al. 2022) as a blazar indicator will improve the completeness significantly. The size of the blazar sample we explored should therefore not be construed as the limit set by nature.

Fig. 1 Ratio of the average bootstrap sample variance to the true population variance of the 16 terms in the l ≤ 2 VSH model we employed for simulated multi-epoch data in which position offsets vary with epoch. The position offsets were sampled from the empirical distribution of statistically significant ICRF3-Gaia offsets observed in the sample used by Secrest (2022) and this work. The error bars denote the standard errors of the means.

Fig. 2 Distribution of l = 0 (rotation) and l = 1 (glide) values between ICRF3 S/X and Gaia DR3, separating objects by fractional variability into quasar-like (Fvar < 0.2), intermediate (0.2 < Fvar < 0.4), and blazar-like (Fvar > 0.4) samples. A dramatic improvement in frame stability and rigidity is enabled by selection on blazars, exceeding even the sources used to define the ICRF3 frame.

3 Results

3.1 Effect of photometric variability class

The results of the bootstrap are shown in Figs. 2, 3, and 4 for S/X, K, and X/Ka, respectively. The difference in the frame stability and rigidity between frames defined using low-variability quasar-like sources and high-variability blazar-like sources for the S/X band is striking, with the difference being less but still notable for K and X/Ka. The mean offsets of the bootstrap samples from zero, most notable at S/X, suggest that selection on quasar-like objects also introduces significant frame distortions. The intrinsic dispersions of the bootstrap samples of all 16 VSH parameters are given in Tables 1, 2, and 3 for S/X, K, and X/Ka. The intrinsic dispersion was estimated by subtracting the statistical error from the total observed dispersion in quadrature. We estimated the statistical error numerically by permuting a fiducial position of each source by the formal covariances of the VLBI and Gaia positions to create a simulated version of the source, where the observed VLBI-Gaia offset is entirely attributable to observation error. We used the Gaia positions as the fiducial positions for convenience because they remained the same for the S/X, K, and X/Ka analyses. Between S/X and Gaia, blazars define a frame that is about six times more stable than quasar-like objects on average, and they are over twice as stable as the defining sources themselves. For K, blazars are similarly as stable as the defining sources, but twice as stable as quasars. At X/Ka, blazars again outperform the defining sources and are again twice as stable as the quasars. A reference frame defined using blazars is therefore consistently and significantly more stable than a frame defined using other object types, and it can be expected to significantly improve upon the current set of defining sources for S/X and X/Ka.

Fig. 3 Distribution of l = 0 (rotation) and l = 1 (glide) values between ICRF3 K and Gaia DR3, separating objects by fractional variability into quasar-like (Fvar < 0.2), intermediate (0.2 < Fvar < 0.4), and blazar-like (Fvar < 0.4) samples.

Fig. 4 Distribution of l = 0 (rotation) and l = 1 (glide) values between ICRF3 X/Ka and Gaia DR3, separating objects by fractional variability into quasar-like (Fvar < 0.2), intermediate (0.2 < Fvar < 0.4), and blazar-like (Fvar > 0.4) samples. The large value the D3 term is a known issue and was discussed by Charlot et al. (2020).

3.2 Optimal weighting based on variability

A key prediction of Secrest (2022) was that variability can be used to optimally weight the optical-radio celestial reference frame link. To test this prediction, we binned the full sample of ICRF3-Gaia matches from Secrest (2022) by F_var and fit the VSH model to each bin, returning the reduced chi-squared fit statistic as a function of F_var. As shown in Figs. 5, 6, and 7, the fit statistic shows a clear dependence on F_var, dramatically improving for sources in the intermediate or blazar-like classes (F_var > 0.2). This suggests a means of weighting the covariance matrix in Eqs. (1) and (2).

We found that a remarkably simple scaling can be applied to the covariance matrix Σ. Noting that F_var is dimensionless, we determined a weighting vector $\mathit{w}={\mathit{F}}_{\text{var}}^{-x}$ such that $$\mathbf{\Sigma }\to \mathit{w}\mathbf{\Sigma }={\mathit{F}}_{\mathrm{v}\mathrm{a}\mathrm{r}}^{-x}\mathbf{\Sigma },$$(3)

where F_var is the vector containing the fractional source variabilities. We chose x such that the survival probability of the χ² distribution was close to ~0.5 given the degrees of freedom (d.o.f.) ν. For the original covariances, this probability is zero for all three ICRFs. We find x_S/X = 1.61, x_K = 1.15, and x_X/Ka = 0.87 with d.o.f. ν_S/X = 3860, ν_K = 846, and ν_X/Ka = 756. Because sources selected on K or Ka are a priori more likely to be flat-spectrum blazars, the smaller values of x make sense.

Using the original covariances, the 16-parameter VSH fit gives formal frame rotational uncertainties smaller than 5 μas for S/X and smaller than 10 μas for K and X/Ka, while after applying Eq. (3), the formal frame rotational uncertainties reach 13 μas for S/X, ~ 19 μas for K, and ~18 μas for X/Ka. Thus, weighting the formal astrometric covariance of a source according to its astrophysical properties, namely its photometric variability, yields more realistic frame alignment uncertainties without the need to apply a single post hoc error correction to all sources, as was done by Charlot et al. (2020).

Repeating the bootstrap-based stability analysis on all 1938 objects, we found that weighting the covariance matrix using Eq. (3) resulted in a substantial improvement in frame distortions and stability, comparing S/X to Gaia (Fig. 8), with an average improvement of a factor of 2 in the distortions and instabilities (Table 4). The improvement is less notable in K, although the distortions and instabilities are slightly reduced for most VSH coefficients (Fig. 9 and Table 5), which is consistent with the weaker dependence of the covariance scaling on F_var because objects selected in K are more likely to be blazars. ICRF3 X/Ka continues to be dominated by the distortions noted by Charlot et al. (2020), so there is no significant advantage to scaling the covariance matrices (Fig. 10 and Table 6).

We note that scaling the astrometric covariance might be interpreted as suggesting a variability-dependent, per-source position uncertainty correction factor, but we caution against this interpretation. The modification of the individual source position errors in a way that effectively scales the VSH fit covariance leads to covariance-normalized individual position offsets X that are heavily skewed to low values because many sources with a low fractional variability are nonetheless blazars or sources that otherwise do not have notable position offsets. Any per-source position error correction method should ideally produce values of X following a Rayleigh distribution with σ = 1. The covariance scaling method shown above only accounts for the average behavior of sources as a class. This allows optimal weighting of the ties between celestial reference frames.

Fig. 5 Reduced chi-squared statistic of the 16-parameter VSH model as a function of fractional variability using the ICRF3 S/X and Gaia DR3 positions, showing that the variability-based covariance scaling given by Eq. (3) correctly weights the source positions.

Fig. 6 Reduced chi-squared statistic of the 16-parameter VSH model as a function of fractional variability using the ICRF3 K and Gaia DR3 positions, showing that the variability-based covariance scaling given by Eq. (3) correctly weights the source positions.

Fig. 7 Reduced chi-squared statistic of the 16-parameter VSH model as a function of fractional variability using the ICRF3 X/Ka and Gaia DR3 positions, showing that the variability-based covariance scaling given by Eq. (3) correctly weights the source positions.

Fig. 8 Sixteen transformation parameters between ICRF3 S/X and Gaia DR3 before and after weighting by Fvar using Eq. (3) for 1938 sources (d.o.f. = 3860). The error bars show the 1σ dispersion from bootstrapping.

Fig. 9 Sixteen transformation parameters between ICRF3 K and Gaia DR3 before and after weighting by Fvar using Eq. (3) for 431 sources (d.o.f. = 846). The error bars show the 1σ dispersion from bootstrapping.

Fig. 10 Sixteen transformation parameters between ICRF3 X/Ka and Gaia DR3 before and after weighting by Fvar using Eq. (3) for 386 sources (d.o.f. = 756). The error bars show the 1σ dispersion from bootstrapping. The large distortions in X/Ka are a known issue for ICRF3 (compare with Fig. 19 in Charlot et al. 2020).

4 Discussion

We reiterate that the bootstrap method we employed is statistically equivalent to sampling the astrometric behavior of a single set of objects across multiple epochs after accounting for the statistical error. It is therefore a direct measure of the frame stability over time. Our results strongly support a picture in which the small line of sight to blazars improves their astrometric properties, and we have presented a method of optimally weighting all sources to improve frame distortions and stability, especially the stability of the S/X frame.

Nonetheless, we anticipate two objections to the prospect of creating a blazar-based celestial reference frame. The first objection is that as their photometric variability implies, blazars are dynamic objects that are expected to exhibit phenomena such as superluminal motion or changes in jet orientation over human timescales, a quality that might make a researcher hesitate to give blazars preference in reference frame work. The second objection is that blazars are relatively rare objects, and there might therefore be concern about their availability and the constraints placed by this lack of availability on the creation of a future reference frame. In Sect. 4.1, we review the dynamical processes observed in jets and discuss why blazars appear to be astrometrically stable in spite of these effects; and in Sect. 4.2 we comment on the availability of blazar-like sources for the creation of a celestial reference frame.

4.1 Long-term astrometric stability of blazars

It is important to note that changes in the line-of-sight jet angle of radio AGN have been observed on human timescales. Using 20 years of VLBI observations, Cui et al. (2023) found that the jet in M87 precesses around an angle of ~$1_{\cdot }^{\circ }25$ over the span of ~11 yr. von Fellenberg et al. (2023) similarly reported precession with a period of ~7 yr in M81, with an amplitude of 7°. von Fellenberg et al. (2023) also reported, however, a slow and apparently linear drift in the jet angle of $0_{\cdot }^{\circ }5$ yr⁻¹, which they showed can be interpreted as the true precession of the jet with a long but almost unconstrained period of between 200 and 1800 yr, indicating that the 7 yr sinusoidal variation is nutation. While precession can change the alignment of jets by a few degrees on timescales of decades, jets can reorient on much longer timescales. Babul et al. (2013) presented a model in which the random orientation of a newly formed geometrically thin accretion disk induces a torque that causes the black hole spin axis to change on ~Myr timescales. X-ray data of clusters indicate that intercluster heating is approximately isotropic, suggesting that disparate episodes of thin-disk formation can swivel the black hole spin axis across ~4π sr over 𝒪(100) Myr.

With nutation on 𝒪(10) yr timescales and precession on 𝒪(> 100) yr timescales, the Doppler factors δ of blazars should remain high over human timescales. Given a typical blazar with Doppler and Lorentz factors around δ ~ Γ ~ 10 (e.g., Liodakis et al. 2017; Zhang et al. 2020, also Lister et al. 2019), corresponding, for example, to β = 0.995 and line-of-sight angle θ ~ 6°, the Doppler factor is expected to vary between δ ~ 4 and δ ~ 20 for a projected nutation amplitude of 5°. The flux-beaming factor, which scales as δⁿ, is therefore expected to vary around a mean value of ~90 within a factor of five, given n = 2 (Cohen et al. 2007), and the jet is expected to continue to dominate over the accretion disk at visual wavelengths, as indicated by the redder color of blazars (e.g., Fig. 4 in Secrest 2022). The optical-radio offset is then the projected distance between the optical photo-center along the jet and the radio peak that defines the VLBI position, ostensibly the frequency-dependent position at which the core becomes optically thin. In sources with estimated radio and optical positions along the jet, the optical position is down-stream of the radio position about 80% of the time at distances of 𝒪(10) pc (e.g., Plavin et al. 2019). Nutation with decade-scale periods is therefore unlikely to change the optical positions of jets, which are at comparable light-crossing time distances, so nutation is expected to mainly affect the apparent positions of the radio cores. The recent review of the well-studied jet in M87 by Hada et al. (2024) summarized the deprojected core positions of the jet as a function of frequency (their Fig. 8), showing that the frequency-dependent core is within ~0.1 pc of the black hole. Then, for the typical blazar described above, the optical-radio offset is expected to be about 0.1 mas for a z = 1 source (h = 0.7) and will remain stable to within about 1 μas over the nutation amplitude.

It is likely, however, that the deprojected radio core positions of typical z ~ 1 sources lie at considerably larger distances from the central engine than what is observed in nearby M87, because VLBI observations are sensitive to far more of the extended jet. The preceding discussion demonstrated only that nutation, which occurs over human timescales but is limited to jet structures within the light-crossing time, is likely to be a negligible source of astrometric variation in blazars. Astrometric variation in the jet might therefore be limited to superluminal motion or to pseudo-motion induced by flares in quasi-static shocks.

In the case of superluminal motion, the maximum apparent velocity is β_app. ~ Γ = 10 and corresponds to θ ~ 6°, the jet line-of-sight angle in our example. If the structures corresponding to the VLBI and optical positions are superluminal to this degree, then they are expected to exhibit proper motions up to ~0.4 mas yr⁻¹, comparable to the average total Gaia proper motion uncertainty of our sample. The proper motion

uncertainties in about 20% of the sample are small enough for such large superluminal motions to be detectable, but as outlined in Sect. 2.3, Secrest (2022) removed objects with formally significant proper motions because these objects were shown by Souchay et al. (2022) and Makarov & Secrest (2022) to often be multisources such as dual AGN. We therefore restored these objects to the sample, which also entailed allowing those with Gaia counterparts not in Gaia-CRF3, yielding 16 objects in addition to the original 1938. We calculated the covariance-normalized proper motions χ using Eq. (1) in Makarov & Secrest (2022) and determined the statistical consistency of χ between blazars and nonblazars using the Anderson-Darling test. When nonblazars are defined as objects with F_var < 0.2 and blazars as those with F_var > 0.4, there is no compelling evidence that the distributions of total proper motions |μ| are different (p = 0.08), but mild evidence that the covariance-normalized proper motions χ do differ (p = 0.03). When we relaxed the definition of a blazar to objects with F_var > 0.2, the statistical power improved and yielded a highly significant (p ~ 0) difference in total proper motions, with blazars having a mean ⟨|μ⟩ ~ 0.4 mas⁻¹ and nonblazars having ⟨|μ|⟩ ~ 0.3 mas yr⁻¹, but the difference in covariance-normalized proper motion is insignificant (p = 0.24), which indicates that the apparently higher proper motions observed in the blazar class are an artifact of larger measurement errors. The sum of the average proper motion covariance matrix in the nonblazars is 0.23 mas² yr⁻², and for blazar-like objects with F_var > 0.2 the sum is 0.31 mas² yr⁻², although it is 0.23 mas² yr⁻² when blazars are defined as F_var > 0.4. This explains the difference in the significance tests. When we instead compare the mean proper motions of objects with significant values of χ, we find that for χ > 3, 4, 5 the mean proper total motion of nonblazars is 0.81 ± 0.18, 0.58 ± 0.16, 1.11 ± 0.41 mas yr⁻¹ and for blazars (F_var > 0.2), it is 0.63 ± 0.12, 0.69 ± 0.28, 0.46 ± 0.22 mas yr⁻¹. There is therefore no significant evidence that blazars exhibit systematically larger proper motions than nonblazars for objects with (formally) significant proper motions.

The question remains, however, why Gaia does not appear to be sensitive to superluminal motions in the jets of these objects. The likely explanation is that given the moderate redshifts of ICRF objects, Gaia measures the full jet, which dominates over individual superluminal knots. For example, the jet in M87 has several superluminal knots (e.g., Biretta et al. 1999), but they are all much fainter than the full jet, which would be subsumed within the Gaia point spread function at z = 1. The optical astrometric stability of blazars might therefore not necessarily extend to low redshift where the centroid may correspond to a superluminal component of the jet. This is already true of low redshift sources in general, however, which can induce spurious proper motions in Gaia astrometry due to source extent (Souchay et al. 2022; Makarov & Secrest 2022). Fortunately, it is generally straightforward to remove extended sources, in the case of Gaia by excluding sources with large phot_bp_rp_excess_factor.

While we did not detect significant differences in the apparent proper motions of blazar- and quasar-like ICRF sources, stochastic flares might cause variability-induced motion (e.g., Makarov & Goldin 2016), which in Gaia data manifest as a significant astrometric_excess_noise parameter. In the original sample of 1938 objects, we again found no significant differences between blazars and quasars. Objects with F_var < 0.2, F_var > 0.2, and F_var < 0.4 all have a mean astrometric_excess_noise of ~0.2 mas and a mean astrometric_excess_noise_sig of ~0.3, which is close to the expected value of $\sqrt{1/2\pi }\sim 0.4$ for a normal distribution with negative values set to zero⁴. About 50% of the objects in both classes have astrometric_excess_noise = 0, and an Anderson-Darling test did not indicate inconsistencies in the astrometric excess noise or its significance between the two classes. In VLBI data, flares might induce strong changes in the apparent position angle of blazar jets, which are common and might arise because the jet opening angle is comparable to the viewing angle (e.g., Lister et al. 2013; Kostrichkin et al. 2025). For a typical position angle standard deviation of ~10° and core offset of ~1 mas (Lister et al. 2013), we might expect a VLBI position dispersive component of ~0.2 mas. However, Lambert & Secrest (2024) showed that ICRF objects with the lowest VLBI position variability exhibit both the highest fractional photometric variability and the highest prevalence of γ-ray detection. This indicates that blazars are in fact more astrometrically stable in VLBI data, so any astrometric variability caused by large changes in jet position angle does not appear to be significant.

It is clear from the considerations discussed above that although AGN jets are highly dynamic phenomena that exhibit precession, nutation, superluminal motion, and stochastic flaring, these phenomena likely do not significantly detract from the astrometric stability gained by giving preference to sources with jets aligned close to the line of sight where relativistic Doppler beaming enhances the consistency between the optical and radio positions and any inconsistencies are minimized in projection. These arguments explain the empirical fact, demonstrated in the preceding sections, that selection on blazars yields a dramatic improvement in the stability of the celestial reference frame, an improvement that exceeds the stability of the current set of ICRF defining sources, as well as the finding that weighting by source variability allows for use of all sources that have variability information. This confirms the prediction made by Secrest (2022).

4.2 Availability of sources

As noted in Sect. 2.3, the sample we used is only a subset of the sources available to construct a frame that takes source variability into account. While we focused on a sample of 180 likely blazars from ICRF3 identified by Secrest (2022), incorporation of other multi-epoch photometry and larger parent samples of VLBI sources will significantly increase the number of identified blazars. Moreover, as we showed in Sect. 3.2, knowledge of the photometric variability of sources allows for a simple and effective means to weight the frame alignment, and this flexibility can create a sample of defining sources that is more uniformly distributed across the sky. This is a priority of the current ICRF. With the upcoming Gaia DR4 (2026) and Gaia DR5 (2031) releases, which will provide two and four times as many months of data as DR3, respectively, it is likely that the completeness of blazar catalogs based on Gaia data will improve considerably. Several other ground-based surveys such as ZTF and LSST will further improve completeness, and it is likely that infrared data from the WISE/NEOWISE catalogs, which span over ten years, will be useful as well. Finally, dedicated near-infrared monitoring campaigns of VLBI sources within the Galactic plane, where reddening prevented detection by Gaia, might prove invaluable for the source selection or weighting of a future infrared astrometric mission such as GaiaNIR (Hobbs et al. 2021).

5 Conclusions

Using ICRF3 sources with measured fractional photometric variability F_var in Gaia DR3, we have conducted the first study of the effect of astrophysical source type (quasar-like, intermediate, or blazar-like) on the stability of the celestial reference frame. Using the 16-parameter vector spherical harmonic model employed by Charlot et al. (2020) to characterize ICRF3, we quantified the time-dependent stability of the reference frame by noting that random sampling with replacement (the bootstrap) is statistically equivalent to sampling a particular class of source (e.g., blazars) over multiple epochs, given the independence between sources. Our main conclusions are:

1. Photometric variability strongly predicts stability of the celestial reference frame, with the most variable, blazar-like objects defining a frame that is six times more stable than the most photometrically constant, quasar-like objects, comparing ICRF3 S/X to Gaia. In the higher-frequency K and X/Ka bands, this difference persists, but is smaller, likely because the sources are a priori more likely to be blazars;

2. At S/X, blazar-like objects define a frame that is over twice as stable as even the ICRF3 defining sources. This suggests that selecting defining sources based on photometric variability is superior to a selection based on source compactness or apparent VLBI positional stability, as was the case for ICRF3 and previous ICRFs, although we emphasize that we showed that blazar-like objects have more stable VLBI positions in Lambert & Secrest (2024). A mild improvement of a factor of 1.4 is also seen in X/Ka, although blazars perform indistinguishably from the defining sources in K;

3. We showed that scaling the original VLBI source position covariances by $F_{\text{var}}^{-x}$ provides a remarkably simple means of optimally weighting the tie between celestial reference frames. The value of the exponent x is steepest for S/X and flattest for X/Ka, as expected if blazars are more prevalent among sources selected at K and Ka. The values we determined in Sect. 3.2 leave no unexplained variance in the 16-parameter VSH transformations between frames and substantially improve the distortions and stability of the link between ICRF3 S/X and Gaia.

Our results build upon the discovery by Secrest (2022) that blazars are systematically less likely to exhibit optical-radio position offsets, and on the discovery by Lambert & Secrest (2024) that blazars also exhibit significantly less astrometric variability in multi-epoch VLBI data (see also Sect. 6 in Blinov & Arshinova 2024). In addition to quantifying the effect of blazar source selection and weighting on frame stability, we further explored the physical mechanisms of this stability, which superficially seems counterintuitive given the dynamic nature of blazars. We argued that time-varying phenomena such as jet precession and nutation are unlikely to be important on human timescales, and that the effect of superluminal motion is likely negligible for moderate-redshift objects, which comprise the majority of the ICRF. These arguments provide a theoretical explanation for the empirical fact, demonstrated here and predicted by Secrest (2022), that blazars are inherently more astrometrically stable, and that source variability can be used to optimally tie celestial reference frames.

As discussed by Lambert & Secrest (2024) and emphasized by Lambert et al. (2025), photometric variability monitoring campaigns are much less observationally expensive than VLBI monitoring campaigns, so regular observations by networks of one-meter-class telescopes can therefore serve to augment reference frame work by better distinguishing the astrometrically stable blazars from the unstable quasars. With nearly 22 000 known VLBI sources (Petrov & Kovalev 2025), our results show that determining the photometric variability properties of these sources has the potential to significantly improve the precision and stability of the celestial reference frame.

Acknowledgements

The authors thank the anonymous referee for their thoughtful review that significantly improved the quality of this work. N.J.S. gratefully acknowledges V. Makarov for helpful mathematical discussions and clarifications. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/Gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/Gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work made use of Astropy⁵: a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013; Astropy Collaboration 2018, 2022), as well as TOPCAT (Taylor 2005).

References

All Tables

All Figures

Fig. 1 Ratio of the average bootstrap sample variance to the true population variance of the 16 terms in the l ≤ 2 VSH model we employed for simulated multi-epoch data in which position offsets vary with epoch. The position offsets were sampled from the empirical distribution of statistically significant ICRF3-Gaia offsets observed in the sample used by Secrest (2022) and this work. The error bars denote the standard errors of the means.

In the text

	Fig. 2 Distribution of l = 0 (rotation) and l = 1 (glide) values between ICRF3 S/X and Gaia DR3, separating objects by fractional variability into quasar-like (Fvar < 0.2), intermediate (0.2 < Fvar < 0.4), and blazar-like (Fvar > 0.4) samples. A dramatic improvement in frame stability and rigidity is enabled by selection on blazars, exceeding even the sources used to define the ICRF3 frame.
In the text

	Fig. 3 Distribution of l = 0 (rotation) and l = 1 (glide) values between ICRF3 K and Gaia DR3, separating objects by fractional variability into quasar-like (Fvar < 0.2), intermediate (0.2 < Fvar < 0.4), and blazar-like (Fvar < 0.4) samples.
In the text

	Fig. 4 Distribution of l = 0 (rotation) and l = 1 (glide) values between ICRF3 X/Ka and Gaia DR3, separating objects by fractional variability into quasar-like (Fvar < 0.2), intermediate (0.2 < Fvar < 0.4), and blazar-like (Fvar > 0.4) samples. The large value the D3 term is a known issue and was discussed by Charlot et al. (2020).
In the text

	Fig. 5 Reduced chi-squared statistic of the 16-parameter VSH model as a function of fractional variability using the ICRF3 S/X and Gaia DR3 positions, showing that the variability-based covariance scaling given by Eq. (3) correctly weights the source positions.
In the text

	Fig. 6 Reduced chi-squared statistic of the 16-parameter VSH model as a function of fractional variability using the ICRF3 K and Gaia DR3 positions, showing that the variability-based covariance scaling given by Eq. (3) correctly weights the source positions.
In the text

	Fig. 7 Reduced chi-squared statistic of the 16-parameter VSH model as a function of fractional variability using the ICRF3 X/Ka and Gaia DR3 positions, showing that the variability-based covariance scaling given by Eq. (3) correctly weights the source positions.
In the text

	Fig. 8 Sixteen transformation parameters between ICRF3 S/X and Gaia DR3 before and after weighting by Fvar using Eq. (3) for 1938 sources (d.o.f. = 3860). The error bars show the 1σ dispersion from bootstrapping.
In the text

	Fig. 9 Sixteen transformation parameters between ICRF3 K and Gaia DR3 before and after weighting by Fvar using Eq. (3) for 431 sources (d.o.f. = 846). The error bars show the 1σ dispersion from bootstrapping.
In the text

	Fig. 10 Sixteen transformation parameters between ICRF3 X/Ka and Gaia DR3 before and after weighting by Fvar using Eq. (3) for 386 sources (d.o.f. = 756). The error bars show the 1σ dispersion from bootstrapping. The large distortions in X/Ka are a known issue for ICRF3 (compare with Fig. 19 in Charlot et al. 2020).
In the text