VLBI Astrometry of Radio Stars to Link Radio and Optical Celestial Reference Frames: Observing Strategies

The Gaia celestial reference frame (Gaia-CRF) will benefit from a close assessment with independent methods, such as Very Long Baseline Interferometry (VLBI) measurements of radio stars at bright magnitudes. However, obtaining full astrometric parameters for each radio star through VLBI measurements demands a significant amount of observation time. This study proposes an efficient observing strategy that acquires double-epoch VLBI positions to measure the positions and proper motions of radio stars at a reduced cost. The solution for CRF link compatible with individual VLBI position measurements is introduced, and the optimized observing epoch scheduling is discussed. Applying this solution to observational data yields results sensitive to sample increase or decrease, yet they remain consistently in line with the literature at the 1-sigma level. This suggests the potential for improvement with a larger sample size. Simulations for adding observations demonstrate the double-epoch strategy reduces CRF link parameter uncertainties by over 30% compared to the five-parameter strategy.


INTRODUCTION
In this era of thriving multi-band astronomy, the establishment and maintenance of a unified celestial reference frame covering all bands is of utmost importance for various research fields such as astronomy, geodesy, deep-space navigation, and more.Currently, the two main bands in which each has established a high-precision celestial reference frame are radio and optical bands.At the radio band, the latest realization of the International Celestial Reference Frame (ICRF; Ma et al. 1998) adopted by the International Astronomical Union (IAU) is ICRF3, based on VLBI measurements of accurate coordinates for 4536 extragalactic radio sources with a coordinate noise floor of 30 as (Charlot et al. 2020).At the optical band, the third data release of Gaia mission (Gaia Collaboration et al. 2016), Gaia DR3 (Gaia Collaboration et al. 2023), contains the astrometric parameters for more than 1.6 million quasar-like sources, defining the Gaia Celestial Reference Frame 3 (Gaia-CRF3) in the optical band with comparable noise floor to ICRF3 (Gaia Collaboration et al. 2022).The conventional axes of the International Celestial Reference System (ICRS) are defined by ICRF, so Gaia-CRF needs to be aligned to ICRF for consistency.It is now aligned with ICRF3 ★ E-mail: zb@shao.ac.cn through common extragalactic objects to an uncertainty of about 7 as (Gaia Collaboration et al. 2022).
Although Gaia aims to provide a globally consistent reference frame for all types of objects, there still exist subtle differences among them, depending on the magnitude, color, position on the celestial sphere, and other factors that produce small shifts of image centroids (Lindegren et al. 2018).There is an indication that the bright CRF ( band magnitude 13) of Gaia DR2 rotates with respect to the quasars by more than 100 as yr −1 (Lindegren et al. 2018;Brandt 2018).This deviation is likely caused by a deficiency in the astrometric instrument calibration model (Lindegren 2020a, Appendix B).For bright ( 13) and faint ( > 13) sources, different "window classes" (WCs) of the calibration model are used.The WCs are schemes for sampling the pixels around a detected source on the charge-coupled device (CCD), and the image centroids derived from different WCs may include systematic errors.The Gaia-CRF3 (Gaia Collaboration et al. 2022) took this rotation into consideration and applied a spin correction = (−16.6,−95.0, +28.3) as yr −1 on three axes [ , , ] derived from the proper motion differences between H (van Leeuwen 2007) and Gaia DR2 (Gaia Collaboration et al. 2018) during its production (Lindegren et al. 2021a).The uncertainty in the link of the H celestial reference frame (HCRF) to ICRF at epoch J1991.25 was ±0.6 mas in each axis (Kovalevsky et al. 1997).Verifying the bright Gaia-CRF through the comparison with HCRF is mainly limited by the position error in HCRF, so this method will be of little use in future Gaia data releases (Lindegren 2020a).
Therefore, other inspections, including both internal and external ones, are needed to realize a consistent CRF at different magnitudes.As an example of internal inspections, Cantat-Gaudin & Brandt (2021) characterized magnitude-dependent systematics in the proper motion of bright sources using resolved binaries and open cluster members spanning above and below = 13, and the systematic spin parameters at magnitude bin 12.75 < < 13 are = (+34.9,+68.9, −2.9) ± (1.7, 2.1, 1.9) as yr −1 .As for external inspections, which are beneficial for the alignment of optical and radio CRFs, objects observable in both bands are needed.Since most quasars are faint in the optical band ( ≥ 14), other kinds of objects are required to estimate the correction parameters of the bright Gaia-CRF, and VLBI astrometry of radio stars is the most direct way (Kovalevsky et al. 1997;Lestrade et al. 1999).To this end, Lindegren (2020a,b) linked bright Gaia-CRF2 to ICRF3 with VLBI astrometric data of the best-fitting 26 radio stars, and obtained a set of link parameters of = (−19, +1304, +553) ± (158, 349, 135) as, = (−68, −51, −14) ± (52, 45, 66) as yr −1 for orientation and spin parameters, respectively.After the release of Gaia-CRF3, because it has improved both the number and precision of sources used to establish the reference frame compared with Gaia-CRF2, the CRF link parameters at the bright end also need to be re-estimated.Bobylev (2022) assembled a sample of 84 radio stars with proper motion measurements, and derived a spin of = (+60, +80, −100) ± (60, 70, 80) as yr −1 between Gaia-CRF3 and ICRF3.Lunz et al. (2023) linked Gaia-CRF3 and ICRF3 using radio stars with improved models of stellar motion and expanded VLBI results, and found that a residual spin of = (+22, +65, −16) ± (24, 26, 24) as yr −1 is still significant for bright sources, suggesting that the proper motion correction based on the H data involved in Gaia DR3 is not effective enough.
The advantage of using VLBI observations of radio stars to determine the link parameters between the bright Gaia-CRF and ICRF is its independence and ability to estimate the orientation difference at a given epoch.It is crucial to conduct the VLBI observations before the end of Gaia mission for an accurate orientation parameter estimation.As a part of our project aims to improve the link accuracy of the bright Gaia-CRF and ICRF by an order of magnitude, in addition to making use of archive data, we have also observed about a dozen radio stars with VLBI and obtained/updated their astrometric parameters (e.g., Chen et al. 2023).However, the number of available radio stars remains limited, and the cost of observation is very high.For typical VLBI astrometric observations, more than four epochs spanning one year are required to solve such many parameters to obtain five-parameter astrometric data, i.e., Right Ascension (R.A.), Declination (Dec.), parallax, and proper motion (in R.A. and Dec.).On the other hand, because of their faintness in the radio band, detecting and locating radio stars requires very long integration times.Therefore, more efficient observing strategies are needed to increase the number of available radio stars.
In this paper, we propose to improve the observing efficiency with a new observing strategy, the double-epoch strategy.The mathematical formulation of a compatible CRF link solution for this strategy is described in Sect. 2. Then we applied this solution to existing observational data in Sect.3. In Sect.4, we generated a simulated dataset, assessed and compared different observing strategies, and estimated the effect of adding different numbers of new radio stars.The method to optimally schedule epochs for the double-epoch observation is de-scribed in Sect. 5. Finally, we summarize and look forward in Sect.6.

SOLUTION
The double-epoch strategy is to observe each radio star for only two epochs, approximately separated by an integer number of years, so that two individual positions can be obtained.The purpose of the integer-year time interval is to cancel out parallactic offsets, so an unbiased proper motion measurement can be obtained, which is essential to the spin parameter estimation.Here, we give a solution on how to incorporate the paired position measurements into CRF link.
The data for a radio star may contain a variety of different forms if new observations (paired position measurements) are added to the historical data.In this case, the solution for five-parameter input data presented in Lindegren (2020a, Sect.2) needs to be slightly modified to be compatible with combined forms of input data, including five-parameter astrometric data and individual positions, etc.The modified solution is referred to as the "compatible solution", and we provide the mathematical formulation of it here in outline.
Suppose there is no translation between the two CRFs, and only a time-dependent rotational offset ( ) exists.Here we assume the angular velocity of the rotation is constant, i.e., ( ) = ( ) + Δ , where Δ is the time difference between VLBI epoch and Gaia reference epoch .Visualizing ICRF and Gaia-CRF as two sets of orthogonal unit vectors [ , , ] and [ ˜ , ˜ , ˜ ] separately, the firstorder coordinate differences between the two CRFs can be expressed as where ( ) = [ ( ), ( ), ( )] ′ denotes the orientation parameters at epoch .The proper motion differences are given by * − ˜ * = + cos sin + sin sin − cos , where * denotes cos , = [ , , ] ′ denotes the spin parameters.
The input data consists of two parts, the Gaia part and the VLBI part.Suppose there are stars in total, the form of the Gaia data for source The parameters to be estimated include the six unknown CRF link parameters = [ ( ), ( ), ( ), , , ] ′ and the corrections to Gaia five-parameter astrometric data (3) These parameters can be estimated using a least-square algorithm.
The computation procedure of the algorithm is briefly described here.
The loss function to minimize is where is the coefficient matrix from Eqs. ( 1)-( 2), is the Jacobian matrix of the propagation model of Gaia data at = 0, and = − is the difference between the observed data and propagated Gaia parameters.
For different forms of VLBI data, different , and are used.In the case of five-parameter data, is a 5 × 1 vector, and is a 5 × 5 matrix.In the case of individual position, the propagation model is based on linear proper motion and parallactic ellipse calculated from ephemeris.The Cartesian coordinates [ ⊙ , ⊙ , ⊙ ] (unit: AU) of the Sun relative to the Earth at each epoch are provided by the Jet Propulsion Laboratory's (JPL) DE440 ephemeris (Park et al. 2021), based on which we have is a 2 × 1 vector, and is a 2 × 2 matrix.Radial velocity, Roemer delay, and other effects are currently not taken into consideration in either case of propagation models, because of their relatively small scale compared to the present measurement accuracy.Other forms of VLBI data are also possible to be included in the solution by modifying the corresponding , and .Minimizing can be achieved by solving the normal equation Substitute Eq. ( 9) into Eq.( 8), ˆ can be first solved from where Then ˆ can be solved from Eq.( 9).The covariance matrix of ˆ is ( =1... ) −1 , and the loss can be calculated from where So the total loss, loss of each star, and the loss of each data item of the star can all be calculated.To assess the goodness of fit, the reduced chi-square can be calculated as The amount of information on ( ) and contributed by each source can be derived from the first and last three diagonal elements of respectively: = trace ( ) and Ω = trace ( ) .
For more detail about the interpretation of the mathematics, see Lindegren (2020a).

APPLICATION TO OBSERVATIONAL DATA
In this section, we applied the compatible solution to observational data as a basis for the following simulation.41 radio stars are included in the final solution in Lindegren (2020a, see Table 1).Lunz et al.
(2023, see Table E.1, E.3) added one-epoch measurements for 32 radio stars to the sample, and applied correlations to the positions of 29 stars in Lindegren (2020a) for homogenization, for example, update their reference source positions to ICRF3 coordinates.55 radio stars in total are included in the solution of Lunz et al. (2023), and we also use this dataset and adopt their correlations here.The uncertainties of the 32 one-epoch measurements we adopted are the inflated uncertainties they provided, taking into account both thermal and systematic errors.We also make use of new measurements of two radio stars, HD 199178 and AR Lac, given in Chen et al. ( 2023).The optical counterparts are collected from Gaia DR3, and we adopted the parallax zero-point correction recipe provided by Lindegren et al. (2021b), which takes -band magnitude, color, and ecliptic latitude into account.
The direct application of the CRF link solution to the dataset will be seriously biased by outliers.To eliminate radio stars with high loss, we applied a similar iteration procedure as Lindegren (2020a); Lunz et al. (2023) used.The solution is first applied to the entire dataset, then the star with the highest / is discarded and the solution is applied again to the remaining stars.The iteration is repeated for = 0, 1, • • • rejected stars, and the evolution of max( / ), 2 red and is shown in Figure 1.In Panel (a), the downward trend of max / flattens out after = 10 (in logarithmic coordinates, same below), while in Panel (b), a similar turning point occurs at = 6.The estimated CRF link parameters shown in Panels (c) and (d) have non-negligible changes as increases, especially the spin parameters.We arbitrarily choose = 16 here.Three subsets of data are used to evaluate the effect of adding new observations to the radio star samples collected by Lindegren (2020a): 41 stars from Lindegren (2020a) only; 55 stars from both Lindegren (2020a) and Lunz et al. (2023); 55 stars from the entire dataset.The same iteration procedure as above is also applied to the first two subsets, and = 11 and 16 are chosen for them respectively.The results are shown in Table 1.
The application to the entire dataset gives a result basically consistent with Cantat-Gaudin & Brandt (2021) and Lunz et al. (2023): Systematic rotation is most significant on the Y-axis, followed by the X-axis, but not on the Z-axis.However, the consistency is not very significant (merely around the 1-level).During the procedure of addition (from subset (1) to ( 2) and ( 3)) and removal (the iteration procedure shown in Figure 1) of stars, the CRF link parameters show significant fluctuations.This can be explained by the small sample size, which is still far from giving a stable CRF link.
The estimated parameters and their uncertainties do not change much with 32 one-epoch measurements added for subset (2) relative to subset (1).As shown in Table 2, the one-epoch observations added in subset (2) only slightly increase the weight of HD 199178 and AR Lac, while the new observations of the two stars given in Chen et al. ( 2023) contribute quite a lot in the solution.We suggest that multiepoch measurements perform better in the CRF link, and this will be verified in Sect.4.2.

SIMULATION FOR FUTURE OBSERVING STRATEGY DESIGN
The unstable CRF link result shown in 3 indicates the need for a larger sample of available radio stars.However, as mentioned in Sect. 1, to obtain five-parameter astrometric data is high-cost.Therefore, in this section, we conducted a simulation to compare the three observing strategies: single-epoch, five-parameter, and double-epoch strategies.For the single-epoch strategy, the positions of stars are only measured once.For the five-parameter strategy, each star is simulated to be observed over six epochs within a whole year, and the astrometric parameters are estimated from the six individual positions.For the double-epoch strategy, the positions of stars are measured over two distinct epochs separated by a full year.To ensure comparable observation costs across the three strategies, we kept the total number of epochs the same.

Simulated dataset
The evaluation of the strategies needs a dataset containing three different forms of simulated data: (1) optical five-parameter astrometric data ] ′ . is the number of radio stars, while is the epoch number of each radio star.Since the goal of the simulation is to help design future observations, the generation of simulated dataset requires a list of radio stars that are observable with VLBI but not yet observed.Here we first collected radio stars (those exhibiting non-thermal radio emission detectable by VLBI, such as RS Canum Venaticorum variables and X-ray binaries) from the Very Large Array Sky Survey (VLASS, Gordon et al. 2021) and Gaia.Specifically, we crossmatched the radio star samples from VLASS with Gaia DR3, and used the Gaia DR3 five-parameter astrometric data at epoch Gaia = 2016.0as the "true" data for dataset generating.However, VLASS only covers the sky area with > −40 • , so we also collected some radio stars distributed in the deep-southern sky from Wendker (2015).The stars already used in Sect. 3 were excluded.To evaluate the impact of sky distribution, the 111 radio stars collected from VLASS were classified into dataset A, while the 35 deep-southern radio stars collected from Wendker (2015) were classified into dataset B. In total, = 146 radio stars are collected.Note that the candidate selection criteria only require that the radio flux density and compactness of the radio star are sufficient for VLBI detection and that there is a Gaia counterpart.So the dataset only reflects the radio star distribution on the celestial sphere, which is shown in Figure 2, and cannot be directly used as a list of candidate sources for future observations.The simulation would be meaningless if no noises were added to the data, since all the input data are perfectly matched in CRF connection.Although the errors of astrometric parameters include both systematic errors and random errors, in this work, only random errors are considered, because systematic errors have different origins and characteristics, making them difficult to describe with simple models in the simulation.With Gaussian noise (using the uncertainties given by Gaia DR3 as standard deviation ) added to the "true" data, the optical five-parameter astrometric data { } are generated.The individual VLBI positions { } are generated using the same propagation model as Eq. ( 7) shows.The number epoch = 6 and time distribution of epochs (evenly distributed between 2023.5 and 2024.5) of all stars are kept the same to preserve the inner consistency of the dataset.Gaussian noise is also added to the VLBI positions and the standard error VLBI = 50 as of the noise is used as position uncertainty.This accuracy has been verified to be achievable through MultiView (Rioja et al. 2017) observation at C-band, e.g., Hyland et al. (2022Hyland et al. ( , 2023)).VLBI five-parameter astrometric data { } were then generated.Here a least-square algorithm is used to estimate the five astrometric parameters and the corresponding covariance matrices from individual position sequences.
All the data generated above are under the same CRF, so we need to rotate the optical data to another CRF.The adopted CRF link parameters = [ ( Gaia ), ( Gaia ), ( Gaia ), , , ] ′ shown in Table 1 are used in the rotation, and the optical five-parameter data Now the rotation is complete and the simulated dataset containing { ˜ }, { }, { } is ready for observing strategy evaluation.

Strategy evaluation
The evaluation procedures for the three strategies are discussed separately below:  (1) For the single-epoch strategy, stars are randomly sampled from { ˜ } and { }.Each star has one VLBI position measurement at either 2023.5 or 2024.5, and corresponding optical five-parameter astrometric data.
(3) For the double-epoch strategy, stars are randomly sampled from { ˜ } and { }.VLBI position measurements at both 2023.5 and 2024.5 are included.
To assess the effect of adding new observation results to existing data, the simulated samples are combined with the "real" data, and the CRF link parameters are then estimated with both five-parameter and double-epoch strategies.
1,000 rounds of sampling were conducted for each strategy, and the mean values of the results are shown in Table 3.The first three lines of Table 3 are at the same observation cost so a comparison can be made between the three strategies.The single-epoch strategy performs significantly worse than the other two strategies, suggesting that multi-epoch measurement is necessary for the CRF link.Between the multi-epoch strategies, the double-epoch strategy is a more efficient choice, i.e., it obtains lower parameter uncertainties at the same cost, which are over 30% lower than that of the fiveparameter strategy.
The comparison between the datasets with and without the deepsouthern radio stars shows that the improvement of the uniformity of sky distribution of radio stars helps to reduce the difference in the uncertainties of the CRF link parameters on the three axes.So it is valuable to add new observations of radio stars distributed in the deep-southern sky.
We compared the results of adding different numbers of new stars to the "real" data with both five-parameter and double-epoch strategies under VLBI = 50 as and VLBI = 100 as respectively, which is shown in Figure 3.The results show that the double-epoch strategy is always a better choice as the added epoch number grows, and the orientation parameters are more sensitive to VLBI position uncertainty relative to the spin parameters.Adding 60 new stars with the double-epoch strategy to the "real" data under VLBI = 50 as obtains 25% and 71% lower uncertainties in orientation and spin parameters respectively, while adding 120 new stars further improves these to 34% and 76%.
The future Gaia data release will improve astrometric parameter accuracy by 40% on the precision of all optical data The three subsets are: (1) 41 stars from Lindegren (2020a), with position corrections from Lunz et al. (2023); (2) 55 stars from both Lindegren (2020a) and Lunz et al. (2023); (3) Same 55 stars as (2) but with new measurements for HD 199178 and AR Lac from Chen et al. ( 2023).Formal uncertainties of the parameters are calculated from the inverse normal matrix.The "Adopted" line is the result of subset (3) but with uncertainties normalized by 2 red , so that the uncertainties can be compared with the simulation results in Table 3.
The sky distribution (Aitoff projection in equatorial coordinates) of the radio stars used in this work.The blue solid diamonds denote the 111 radio stars collected from VLASS, while the purple dots denote the 35 radio stars distributed in the deep-southern sky collected from Wendker (2015).The red hollow squares and the orange hollow triangles indicate the 39 accepted and the 16 rejected "real" radio stars, respectively.The Galactic Equator is plotted in green.and Ω are the weights for the determination of orientation and spin parameters, respectively.
is the sum of the dimensions of all data of a star included in the solution.
products and a factor of almost three for the proper motions (Gaia Collaboration et al. 2021).So we tried to multiply all Gaia DR3 parameter uncertainties by a factor Gaia = 0.5 to see how much this will contribute to CRF link accuracy improvement.The result shows that the contribution is insignificant, indicating that the quantity and quality of VLBI data are more important limiting factors.
Although the simulation is not a perfect reflection of the real situation, the conclusions will help to design future observations of radio stars for the link of the bright Gaia-CRF to ICRF and quantitatively estimate the effect of adding new observations.

Discussion on radio-optical offsets
The radio-optical offsets caused by the physical properties of radio stars are not considered in the simulation, such as differences in radiation regions arising from varying radiation mechanisms.Assuming that the directions and magnitudes of the offsets of different stars are random, a larger sample size will statistically result in averaging out the offsets to zero in the aggregate.
The double-epoch strategy will add about three times more new stars than the five-parameter strategy, therefore in practice the influence of unmodeled offsets will decrease, and the CRF link result will be more robust.In practice, a large proportion of radio stars are discarded due to their excessive residuals, which is inevitable.In the column "dataset", A denotes that only radio stars with > −40 • collected from VLASS are used, A+B denotes that radio stars collected from Wendker (2015, < −40 • ) are added, and A+B+R denotes that radio stars are collected from dataset A+B in addition to "real" data, thereby evaluating the effect of adding new observations.Gaia denotes the multiplier for the uncertainties of Gaia astrometric parameters. 2 red is the reduced chi-square of the least-square estimation.For the first four lines of the table, the total observation costs star × epoch are kept the same to compare the efficiency of the strategies.Formal uncertainties of the parameters are calculated from the inverse normal matrix.Since the uncertainties of the simulated data are consistent with the standard deviation of the Gaussian noise added to the parameters, normalization is not needed for the simulated data, while the covariance matrices of all "real" data are scaled by their reduced chi-square to normalize the result.* star and epoch are the numbers of stars and epochs sampled from the simulated dataset respectively, do not include the "real" data used in the last two lines of this table.
The larger sample size will also allow some of the samples to be discarded, rather than gambling on a few.Observations before the end of Gaia mission will contribute more to the derivation of the orientation parameters, thus eliminating the instrument error and allowing the true radio-optical offsets to be derived, and the physical nature of radio-optical offsets can then be studied.Although the observation of Gaia mission is coming to its end, the closer VLBI epoch is to Gaia reference epoch, the larger contribution it will make.This is another valuable aspect of these observations.

PRACTICAL OBSERVING EPOCH SCHEDULING
In this part, we discuss how to practically schedule the observing epochs for the double-epoch strategy.The strategy is designed not to be affected by parallax, however, the time interval between the two epochs is practically difficult to be exactly an integer number of years, which will bring parallactic displacement.In the compatible CRF link solution, this displacement is corrected by an offset calculated from Gaia parallax .Suppose there exists a bias Δ between and the "true" parallax true , the displacement would not be perfectly corrected.This will affect the independence of the CRF link, and reduce the accuracy of the position and proper motion of the radio star.It is necessary to quantitatively estimate the impact of Δ , and we give the mathematical form for the estimation here.
The parallactic displacement is caused by the change in the position of the Sun relative to the Earth on an annual cycle.Since the error estimation here does not need high accuracy, the Earth's orbit can be approximated to be a circle.So the ecliptic longitude of the Sun at epoch can be expressed as where 0 = 280 • is an approximation of the ecliptic longitude of the Sun at the beginning of a year.The Cartesian coordinate [ , , ] (unit: AU) of the Sun can be calculated as where = 23.4 • is an approximation of the obliquity of the ecliptic.The impact of Δ on position can then be derived: Δ ( ) = Δ (sin ⊙ sin cos − cos ⊙ cos sin − sin ⊙ cos sin sin ) .(21) If the time interval Δ between the two epochs is close to an integer number of years, a linear approximation can be used to approximately estimate the impact of Δ on proper motion.Differentiate Eqs. ( 20)-( 21) by : Δ ( ) ′ = 2 Δ (cos ⊙ sin cos the impact of Δ on proper motion can be expressed as where Δ = |Δ − round(Δ )| is the time difference between Δ and its proximal integer number of years.Assuming Δ ≤ 0.05, which is feasible in actual epoch scheduling, the as yr −1 -scale bias is negligible with current measurement precision ( = 0.02-0.03mas at < 15, Gaia Collaboration et al. (2021)).
For epoch within possible observing time, plotting the trends of the impact of Δ will help to find the best epoch.Two coefficients can be used to quantitatively express the impact of Δ on position and proper motion respectively: Δ −1 in the equation is for cancelling out Δ in Eqs. ( 20)-( 23), therefore is a function of epoch and position ( , ), while is in addition related to Δ /Δ .Two examples are shown in Figure 4, with fixed Δ /Δ = 0.05 and different positions in the sky.For the one far from the Ecliptic, the ranges of variation of and are small, and is always higher than because of the extra coefficient 2 Δ /Δ of ; while for the other one close to the Ecliptic, and change significantly with time, and is sometimes lower than .Therefore, for most cases that stars are far from the Ecliptic, the impact of Δ is not very sensitive to epoch ; While for the stars close to the Ecliptic, it is important to carefully choose the epoch to minimize the impact on position or proper motion depending on the needs: A smaller will enable the star to contribute more to the orientation parameters, while conversely a smaller will enable the star to contribute more to the spin parameters.
Since is also (even more) sensitive to Δ /Δ other than the choice of epoch , it is possible and better to reduce the impact of Δ on proper motion through reducing Δ /Δ : extend the time interval between two epochs Δ to more years, and/or limit the time interval even closer to an integer number of years.
As a conclusion, for the minimization of the impact of Δ , it is advised to choose a pair of epochs with a lower , and extend the time interval between epochs and tightening constraints on observation dates to reduce .

SUMMARY AND FUTURE OUTLOOK
In this study, we proposed a double-epoch observing strategy for radio star VLBI observation, aiming to validate the consistency of Gaia-CRF at the bright end with high efficiency.Firstly, we introduced a CRF link solution that is compatible with various forms of input data, and then applied the solution to currently available data, which result indicates the need for a larger sample size.We conducted a simulation to compare three strategies for additional observations: single-epoch, five-parameter, and double-epoch strategy, and the results show that the double-epoch strategy is the most efficient one among them.Specifically, the parameter uncertainties of the double-epoch strategy are over 30% lower than those of the five-parameter strategy at the same observation cost.
We also compared the datasets with and without the radio stars distributed in the deep-southern sky ( < −40 • ), and the result shows that the improvement of the uniformity of sky distribution of radio stars makes the uncertainties of CRF link parameters on the three axes more even.We found that the bottleneck of the link accuracy is the quantity and quality of VLBI astrometric data, and the future improvement of Gaia data will not contribute much.Therefore, the addition of high-quality VLBI measurements is the key to improving the link accuracy.Under the conditions set in our simulation (epochs close to 2024.0, position accuracy of VLBI = 50 as, etc.), if the double-epoch strategy is followed to observe 60 new radio stars, the uncertainties of the orientation and spin parameters will be reduced by 25% and 71%, respectively.
We discussed how to schedule new observations based on the double-epoch strategy, and gave a recipe to estimate the impact of possible Gaia parallax bias on position and proper motion measurement.We found that the priority of epoch choice should be given to reducing the impact on position in most cases, while reducing the impact on proper motion can be more effectively achieved by extending epoch interval and tightening constraints on observation dates.
The new observation can also be combined with archival data, e.g., a radio star was once observed but no proper motion was derived, or the proper motion derived previously is not accurate enough, adding one or two epochs is an effective way to obtain an accurate proper motion measurement.The archival data with a long time interval between its epoch and the Gaia reference epoch contributes very little to the estimation of orientation parameters, in this case, adding epochs will greatly increase its contribution.It is also possible to add additional epochs after the double-epoch observation is complete to obtain a full five-parameter astrometric measurement.
Significant improvements in the accuracy of the link of the bright Gaia-CRF to ICRF would require a substantial number of new observations.We have proposed observing some radio stars, including both previously observed and new targets, to improve the proper motion accuracy of archive data and increase the number of available radio stars.In addition, the new facilities anticipated in the near future, such as the Next Generation Very Large Array (ngVLA) and the Square Kilometer Array (SKA), will greatly augment our capacity for observing radio stars by enabling the detection of currently unobservable faint targets.With a large number of new observations of radio stars, a more precise CRF link model can be adopted, i.e., taking into account effects such as Galactic aberration, magnitude/color dependence of the CRF link parameters, etc.The double-epoch strategy can also be applied to other frame links, such as the link between the reference frames realized by VLBI astrometry and pulsar timing respectively.

Figure 1 .
Figure 1.The evolution of loss and CRF link parameters with rejected star number grows on the full dataset.(a) The max ( ) in the solution.(b) 2 red of the solution.(c) Estimated orientation parameters.(d) Estimated spin parameters.Error bars are formal uncertainties directly derived from the covariance matrix ( =1... ) −1 .

Figure 3 .
Figure 3.The results of adding different numbers of new stars with the fiveparameter and double-epoch strategy to the "real" data.star × epoch is the total epoch number: for the five-parameter strategy, star = 6; for the doubleepoch strategy, star = 2. ¯ (top panel) and ¯ (bottom panel) are the mean uncertainties on three axes of orientation and spin parameters, respectively.D50, D100, F50, F100 denote double-epoch strategy with VLBI = 50 as, double-epoch strategy with VLBI = 100 as, five-parameter strategy with VLBI = 50 as, and five-parameter strategy with VLBI = 100 as, respectively.

Figure 4 .
Figure 4.The impact of Δ on position and proper motion as function of epoch .Due to the symmetry of the parallactic ellipse, an epoch range of only half a year is shown in the figure.Δ /Δ is fixed to 0.05.Top panel: far from the Ecliptic, = 16 ℎ , = −60 • , and Ecliptic latitude ≈ −38 • ; Bottom panel: close to the Ecliptic, = 3 ℎ , = +15 • , and ≈ −2 • .The red solid line and blue dashed line denote and respectively.

Table 1 .
Results of the application to Gaia DR3 and VLBI observational data

Table 2 .
Solution statistics for HD 199178 and AR Lac

Table 3 .
Results of the simulation evaluating the observing strategies