Single-band VLBI Absolute Astrometry

The ionospheric path delay impacts single-band, very long baseline interferometry (VLBI) group delays, which limits their applicability for absolute astrometry. I consider two important cases: when observations are made simultaneously in two bands, but delays in only one band are available for a subset of observations; and when observations are made in one-band design. I developed optimal procedures of data analysis for both cases using Global Navigation Satellite System (GNSS) ionosphere maps, provided a stochastic model that describes ionospheric errors, and evaluated their impact on source position estimates. I demonstrate that the stochastic model is accurate at a level of 15%. I found that using GNSS ionospheric maps as is introduces serious biases in estimates of declination and I developed a procedure that almost eliminates them. I found serendipitously that GNSS ionospheric maps have multiplicative errors and have to be scaled by 0.85 in order to mitigate the declination bias. A similar scale factor was found in comparison of the vertical total electron content from satellite altimetry against GNSS ionospheric maps. I favor interpretation of this scaling factor as a manifestation of the inadequacy of the thin-shell model of the ionosphere. I showed that we are able to model the ionospheric path delay to the extent that no noticeable systematic errors emerge and we are able to assess adequately the contribution of the ionosphere-driven random errors on source positions. This makes single-band absolute astrometry a viable option that can be used for source position determination.


INTRODUCTION
Method of VLBI absolute astrometry involves observations of many active galactic nuclei roughly uniformly distributed over the sky.Data analysis of group delays derived from these observations is usually performed in the accumulative mode, i.e. all VLBI absolute astrometry and geodesy observations are processed in a single least square solution for estimation of source coordinates, Earth orientation parameters, stations positions, and nuisance parameters, such as atmospheric path delay in zenith direction and clock function.The errors in source position adjustments are mainly due to the ther-mal noise and the inaccuracy of modeling path delay in the atmosphere.They vary in a wide range from 0.05 to 50 mas depending on source flux density, network geometry, and the number of observables, with 1 mas error being typical.The grid of sources which positions are determined from dedicated absolute astrometry observing campaigns forms the basis that makes possible differential astrometry that involves observations of pairs of targets and calibrators with known positions.This distinction is often blurred in literature for brevity.We should be aware of that differential astrometry, despite being capable to determine very precise position differences, in principle cannot provide positions more precise than positions of the calibrators determined with absolute astrometry and it inherits its random and systematic errors.

Petrov
Since the contribution of the ionosphere to group delay is reciprocal to the square of the observing frequency, usually absolute astrometry experiments are performed at two widely separated bands, 2.3/8.4 or 4.3/7.6GHz (Petrov 2021).Processing of dual-band data allows us to eliminate effectively the contribution of the ionosphere.However, there are situations when absolute astrometry observations are available only at one band.
We know that compared with dual-band observations, single band absolute astrometry observations are affected by the contribution of the path delay in the ionosphere.Two questions arise: 1) which data analysis strategy does provide source coordinate estimates with the lowest uncertainties and 2) how to account for the contribution of errors in the ionosphere modeling to reported source coordinate uncertainties?The goal of this study is to provide answers to these questions.The dualband group delays are considered the ground truth free from the impact of the ionosphere in the framework of this study.I took several trial datasets of dual-band VLBI observations from twenty-four hour observing sessions and used them as a testbed.I dropped existing observations at the second band during data analysis and compared results of single-band data against the reference solution that used both bands.I should note that the analysis presented in this article is specific to source positions determined with absolute astrometry and is not applicable to position differences determined with differential astrometry that has its own error model.

VLBI DATASET USED FOR TRIALS
The primary dataset used for analysis is 263 twenty four hour experiments at the Very Large Baseline Array (VLBA) network since April 1998 through May 2021.All these data are publicly available at US National Radio Astronomy Observatory (NRAO) archive1 .The dataset consists of observing sessions under regular VLBI geodesy program RDV (Petrov et al. 2009), astrometric VCS-II program (Gordon et al. 2016), its follow-ups VCS-III and VCS-IV, and geodetic CONT17 campaign (Behrend et al. 2020).The motivation of this choice is to have a long history of observations, a homogeneous network, and both short and long baselines.VLBI absolute astrometry programs at declinations above −40 • are almost exclusively run with VLBA.Therefore, conclusions made from processing trial runs at the VLBA can be propagated directly to the past and future astrometry programs at that network.

Dual-band observations
The impact of the ionosphere dispersiveness on fringe phase is reciprocal to frequency in the first approximation.Therefore, fringe phase in channel i in the presence of the ionosphere becomes where τ p and τ g are phase and group delays, f i is the frequency of the ith spectral channel, f 0 is the reference frequency and where N v -electron density, e -charge of an electron, m e -mass of an electron, o -permittivity of free space, and c -velocity of light in vacuum.Integration is carried along the line of sight.Having substituted values of constants (Klobuchar 1996) and expressing the total electron contents along the line of sight N v ds in 1 • 10 16 electrons/m 2 (so-called TEC units or TECU), we arrive to α = 1.345 • 10 9 sec/TECU.Phase and group delay are computed using fringe phases φ i with weights w i using least squares.The result can be expressed analytically after some algebra: where τ if is the ionosphere-free group delay, TEC is N v ds expressed in TEC units, and f e is the effective ionospheric frequency that depends on weights of spectral channels w i .Typically, the effective ionospheric frequency is within several percents of the central frequency of the observing band.
The best way to mitigate the impact of the ionosphere on group delay is to observe simultaneously at two or more widely separated frequency bands.Then the following linear combination of two group delays at the upper and lower bands, τ u and τ l , respectively is ionosphere free: Here f u and f l are effective ionospheric frequencies at the upper and lower bands respectively.
The residual contribution of the ionosphere in dualband combinations is caused by systematic errors, namely a) higher order terms in the expansion of the dispersiveness on frequency (Hawarey et al. 2005); b) the contribution of frequency-dependent source structure, and c) the dispersiveness in the signal chain.These contributions affect group delay at a level of several picoseconds and they are considered insignificant with respect to other systematic errors.

The use of TEC maps for computation of ionospheric path delay
TEC maps (Hernández-Pajares et al. 2009) also known as Global Ionospheric Model (GIM) derived from analysis of Global Navigation Satellite System (GNSS) data are used for reduction of single band observations.In particular, I used CODE TEC time series (Schaer 1999) available at ftp://ftp.aiub.unibe.ch/CODEsince January 01 1995 with a spatial resolution of 5 • × 2.5 • .Time resolution varied.It was 24 h since 01 January 1995 through March 27, 1998; 2 h since March 28, 1998 through October 18, 2014; and 1 h after that date.The ionosphere is considered as a thin shell at the constant height H i of 450 km above the mean Earth's radius.The ionospheric contribution is expressed via TEC as where M (e) is the so-called thin shell ionospheric mapping function and R ⊕ is the mean Earth's radius and e gc is the geocentric elevation angle with respect to the radius vector between the geocenter and the station.Computation of TEC value at a given moment of time is reduced to computation of the position of the ionosphere piercing point at a given azimuth and elevation and interpolation of TEC at the piercing point at that latitude, longitude, and time.TEC maps from GNSS is a coarse model of the ionosphere.Errors of τ i computed according to expression 5 are greater than the residual ionosphere contribution of ionosphere-free linear combinations of group delays.Therefore, dual-band delay observables are preferred when they are available.However, there are two cases when they are not available: a) dual-band observing sessions with some source detected only at one band; b) single-band observing sessions.In these two cases we resort to computation of the ionospheric contribution to path delays using GNSS TEC maps and evaluation of uncertainties of these contributions.

Ionospheric contribution in dual-band observing
sessions when a source is detected at one band only The simplest way to deal with a mixture of dual-and single-band data is to process experiments three times: 1) using dual-band data of those observations that detected a source in both bands, 2) using low band data, and 3) using upper band data with applying ionospheric path delay computed from GNSS TEC maps.However, typically only a fraction 2 to 20% of observations is detected at only one band; the rest of observations are detected at both bands.Therefore, we can use available dual-band observations at a given observing session to improve the TEC model.
I represent ionospheric path delay at stations j, k as where b j (t) = n i B 0 i (t)b ij is a delay bias expanded over the B-spline basis of the 0th degree, a j (t) = n i B 3 i (t)a ij is the TEC bias expanded over the B-spline basis of the 3rd degree.φ, λ are coordinates of the ionosphere piercing point that depend on positions of observing stations as well as azimuths and elevations of observed sources.
The clock bias occurs due to path delay in VLBI hardware that is different at different bands.This bias is constant for most of the experiments, however occasionally breaks may happen at some stations.Epochs of these breaks coincide with the epochs of breaks in the clock function.Expansion over the B-spline basis of the 0th degree accounts for these breaks.(B-spline of the 0th degree is 1 within the range of knots [i, i+1] and 0 otherwise.)I estimated parameters a j for all the stations and b j for all the stations except the one taken as a reference using all available dual-band observations of a given experiment using least squares with weights  where σ(τ u ) and σ(τ u ) are group delay uncertainties and y is the error floor, 12 ps, introduced to avoid observations with very high signal to noise ratios to dominate the solution.
The time span of the B-spline knot sequence for TEC bias in my solutions was 15 minutes.I applied constraints on the value of the B-spline coefficients and on the first and the second derivatives with the reciprocal weights 5•10 −10 s, 4•10 −14 , and 2•10 −18 s −1 respectively.These constraints were introduced to ensure the continuity of biases and to prevent a singularity in rare cases when too few available observations at a given station could be used for bias estimation at a given spline segment.Figure 1 illustrates estimates of the ionospheric bias.
The resulting total electron contents model TEC j (t)+ a j (t) is more precise than the a priori TEC j (t) taken from GNSS maps because it uses additional information.Using estimates of a j and b j spline coefficients, I compute τ i (t) and its uncertainty according to the law of error propagation using the full variance-covariance matrix of spline estimate coefficients.In order to evaluate the realism of these errors, I processed the trial dataset and computed τ i (t) using the estimates of clock and TEC biases and compared them with the ionospheric contribution derived from dual-band observations.I removed clock biases from VLBI dual-band ionospheric contributions τ vi , formed the differences τ i − τ vi , and then divided them by σ(τ i ) derived from the variancecovariance matrix of a j , b j .I generated the normalized histogram from the dataset of 4,343,782 differences and computed the first two moments of the empirical distribution shown in Figure 2. The fitting parameters of the first and second moments of the distribution are 0.003 and 0.889 respectively.Two factors cause a deviation of the second moment from 1.0 in the opposite way: a) TEC variations not accounted by the parametric model; b) statistical dependence of the estimates of a j , b j , and VLBI path delay.After scaling the variancecovariance matrix by square of 0.889, the distribution of the normalized residuals becomes close to Gaussian.The closeness of the empirical distribution to the normal distribution provides us a confidence that the extra noise introduced by the mismodeled ionospheric path delay after applying clock and TEC biases is properly accounted for.
The closeness of the distribution of normalized differences to Gaussian with σ = 1 is encouraging, but it does not guarantee that the residual errors of the sum of TEC from GNSS and TEC bias adjustments cause no systematic error in estimates of source positions.To characterize the impact of residual errors of the ionospheric contribution on source position, I ran solution XIA that had the following differences with respect to the reference solution: 1) it used X-band group delays; 2) data reduction for the ionosphere accounted for both a priori TEC from GNSS maps and the ionosphere bias adjustment (expression in the parentheses 8); 3) errors of the ionospheric biases σ = a ij α/f 2 e were added in quadrature to reciprocal weights of observables.It should be stressed that parameterization, editing, and delay uncertainties were exactly the same as in the reference dual-band solution.In all trial solutions I varied only input observables, ionosphere-specific reduction model and added noise to the delay uncertainties.Since I analyzed only the differences in solutions, the contribution of other model deficiencies, such as path delay in the neutral atmosphere is canceled when I formed the differences.
I ran also solution SIA that differed from XIA by using S-band group delays and the S-band effective ionospheric frequency f e .For control, I ran solutions XIN and SIN that used the ionosphere-free combinations of group delays, the same weights as XIA and SIA solutions respectively, and the zero mean Gaussian noise with σ = a ij α/f 2 e was added to each observable.The differences in XIA and SIA solutions with respect to the reference solution provide us a measure of the impact of residual ionospheric errors on source position.The differences in XIN and SIN solutions characterize the impact of the residual ionospheric errors if they were Gaussian and totally uncorrelated.
Figure 3 shows the distribution of differences in declinations of source position estimates from XIA and XIN solutions.There is no noticeable deviation from the Gaussian shape.Table 1 lists first two moments of the distribution of position differences in both right ascension and declination.The second moment from SIN solution is close to 1.0, while the second moment from XIN solution is 0.56.The errors of the ionospheric contribution at S-band dominate the error budget.These errors are 14 times less at X-band and are only a fraction of overall group delay errors.The mean biases in right ascension and declination are negligible.The second moment of position estimates from XIA and SIA solutions is 15% higher than moments from positions from XIN and SIN solutions respectively.This increase occurs due to non-randomness of residual ionospheric errors and can be viewed as a measure of unaccounted systematic errors in source positions due to the ionosphere.Analysis of these trial solutions demonstrates that we are able to predict the impact of ionospheric errors on source position with an accuracy of 15% when TEC bias is estimated using dual-band observations.

Ionospheric contribution in single-band observing sessions
TEC biases cannot be computed when an entire session is observed only at one band.Therefore, we have to resort to deriving a regression model to provide estimates of these errors.In the past, Petrov et al. (2019) derived a regression against the so-called global TEC: the integral of TEC over the entire Earth following ideas of Afraimovich et al. (2008) and Krásná & Petrov (2021) derived a regression against the root mean square (rms) of the total ionospheric path delay from GNSS TEC.In  1.The first and second moments of position differences of trial solutions XIA, SIA, XIN, and SIN with respect to the position estimates derived from the reference solution.
The a priori TEC from GNSS and bias adjustment from dualband solutions were applied in XIA and SIA solutions (the 2nd column).The zero-mean Gaussian noise with σ equal to the uncertainty in path delay from the TEC bias adjustment was added in XIN and SIN solutions (the 3rd column).
this study I use the second approach with slight modifications.Following general results of the turbulence theory (see Tatarskii 1971), we can expect that fluctuations at scales x are related to fluctuations at scales y via a power law.I processed the same dataset of 263 twenty-four hour VLBI experiments that was used in the previous section and computed residual ionospheric path delay for each observation as where τ gi is the vertical ionospheric path delay from GNSS TEC maps, τ vi is the vertical ionospheric path delay from VLBI, b j − b k are contributions of the clock bias, and M is the averaged ionosphere mapping function between stations 1 and 2 of a given baseline: M = (M (e 1 )+M (e 2 ))/2.0.The clock biases are routinely adjusted during analysis of VLBI observations and therefore, their contribution on VLBI results, such as source positions, is entirely eliminated.Subtracting them in expression 10, I eliminate their impact on statistics as well.I used only twenty-four hour VLBI experiments for deriving statistics because the ionospheric path delay strongly depends on Solar time, especially at low latitudes, and statistics derived at shorter time intervals may not be representative.
Figure 4 shows the dependence of the rms of residuals τ r on the rms of total ionospheric path delay from GNSS TEC maps τ gi .Each point on the plot corresponds to the rms for a given baseline and a given observing session.I confirm the early result of Krásná & Petrov (2021) but here I used a much larger dataset.The result reported in Krásná & Petrov (2021) was slightly affected by an error in computation of ionospheric group delays for a case when some data are flagged because of radio interference.This error has been fixed and the affected

Normalized residual
Normalized residual Figure 3.The distribution of normalized differences in declination from the trial VLBI solutions with respect to the reference dual-band solution.Left: solution with GNSS TEC maps + adjustments of TEC biases using dual-band VLBI group delays (XIA).Right: solution with added Gaussian noise with σ equal to errors in path delay that correspond to TEC bias adjustments uncertainties (XIN).Blue thick lines show fitted Gaussian distributions with the first and second moments listed in Table 1.
Table 2. Coefficients of the B-spline expansion of the dependence of the rms of residual ionospheric path delay derived from GNSS TEC maps on the rms of the total ionospheric path delay from GNSS TEC maps at 8 GHz.experiments have been reprocessed from scratches.This dependence can be coarsely described as a square root of the rms of the total ionospheric path delay.For a better approximation I sought a regression in the form of an expansion over B-splines of the 3rd degree.The spline coefficients computed using least squares are listed in Table 2.
Using that regression, I developed the following algorithm for computation of errors of the ionospheric path delay from GNSS TEC maps.First, coordinates of K points uniformly distributed over the sphere are computed using a random number generator.Then for each baseline and each time epoch azimuth and elevation angles of that point are computed at both stations of the baseline, and if the elevations above the horizon are greater than 5 • at both stations, that point is selected for further computations.If not, the next point is drawn.Then total ionospheric path delay τ i (A 1 , e 1 , A 2 , e 2 ) is computed using GNSS TEC maps.It is worth mentioning here that unlike to troposphere path delay, τ i (A, e) = τ i (A, π/2)•M (e), since path delay depends on positions of the ionosphere piercing points.It is not sufficient to compute the ionospheric path delay in zenith direction and then map it via M (e): latitude and longitude of the piercing point can be as far as 1000 km from the station.Following this approach, we sample piercing points uniformly distributed within the mutual visibility zone.The process is repeated for 1440 time epochs that cover the time interval of VLBI experiment under consideration with a step of 1 minute.Then for each baseline σ(τ gt ) is computed over a time series of 1440 τ i values.Finally, the estimate of the rms of residual ionospheric path delay derived from GNSS TEC maps is computed from this regression via the rms of the total ionospheric path delay as Baseline-dependent datasets are considered independent for this computation: the mutual visibility at all the stations of the network at a given moment of time is not enforced.For several baselines longer than 96% Earth diameter, this algorithm has a poor performance for selecting points above 5 • .Therefore a minor modification is made for such an extreme case: the elevation angle is fixed to 5 • , mutual visibility is not enforced, and azimuths are selected randomly within a range of [0, 2π] independently for both stations.
In order to evaluate the validity of this regression model of residual ionospheric path delay errors, I computed τ r from dual-band observations, σ gt following the algorithm described above for 263 twenty-four experiments, and then computed the histograms of normalized residuals τ r /σ rr .The histogram is presented in Figure 5.The first two moments of the distribution are −0.083 and 1.214 respectively.Since regression σ rr was found using least squares, the number of observations with σ(τ gr ) less and greater than σ rr for given σ(τ gt ) is approximately equal: the thick blue line cuts the cloud of green points in Figure 5 almost by half.However, the variance of the contribution of those points with σ(τ r ) > σ rr overweights the contribution of those points with σ(τ r ) < σ rr because variance quadratically depends on τ r .This causes a positive bias.After multiplying σ rr by 1.214, the distribution of normalized residuals becomes almost Gaussian.
One can notice that M 2 (e 1 ) + M 2 (e 2 ) in expression 11 is not the same as M = (M (e 1 ) + M (e 2 ))/2 used for computation of the regression.I found that using M instead of M 2 (e 1 ) + M 2 (e 2 ) decreases the second moment from to 1.214 to 1.196, which is negligible, The distributions shown in Figures 2 and 5 are computed for the entire dataset of 4.3 million path delays and they represent the general population over the interval of 23 years.Statistics for an individual observing session may differ.In order to evaluate the scatter of the statistics, I computed the time series of second moments of the distribution of normalized residuals of ionospheric path delays and their uncertainties with and without TEC biases adjusted for each observing session separately.I divided the normalized residuals by scaling factors of 0.889 and 1.196 respectively.I computed the distribution of second moment estimates and showed it in Figure 6.The scatter of the second moments is small when TEC biases are adjusted.That means this statistics is robust.When TEC is not adjusted, the scatter  is significantly larger, but even in that case 90% of the second moment estimates deviate from 1.0 by no more than 30%.This provides us a measure of uncertainties in computation of ionospheric path delay errors in individual single-band observing sessions.

The impact of the residual ionospheric errors in source position in a case of single-band observing sessions
Error analysis presented in the previous section characterizes our ability to predict the first and second moments of the distribution, but it does not guarantee that residual errors due to the ionosphere cause no systematic errors in source positions.I ran trial solutions XIT and SIT that used X-and S-band group delays respectively, applied ionospheric path delays derived from GNSS TEC maps, and inflated reciprocal weights of observables by adding in quadrature errors in ionospheric path delays from TEC maps.
Analysis of differences in source positions from trial XIT and SIT solutions with respect to the reference dual-band solution revealed no peculiarities in right ascensions, but revealed significant systematic errors in declinations (see Figure 7).The pattern of systematic errors in declination from the S-band solution is similar but greater by a factor of f 2 x /f 2 s ≈ 14.We cannot consider a solution with such errors as satisfactory.This was unexpected because prior work of Sekido et al. (2003) I made a number of trial solutions.The following leads turned out productive: 1) to modify mapping function and 2) to scale ionospheric path delays from TEC. Schaer (1999) suggested the following modification of the ionospheric mapping function:   The declination bias is reduced when the parameter ionosphere height in the mapping function is increased.I interpret this result as a deficiency of GNSS TEC maps and I associate the origin of this deficiency with an oversimplification of the mapping function that was used for derivation of the TEC maps from processing of GNSS observations.
Close analysis of these figures reveals that the declination bias has three components: 1) a constant; 2) a linear increase in the declination bias with a decrease of declination; 3) a feature δ(x sin δ + y cos δ) with the minimum at approximately −12 • , where δ is declination.All three components depend on parameters of the mapping function ∆H, α and on scaling factor k. Selecting optimal ∆H, α, k, we can reduce declination biases.I selected ∆H = 56.7 km, α = 0.9782, and k = 0.85.This choice makes the weighted mean declination bias over all sources 0.013 mas.There is an element of subjectivity in this specific choice, since there exist another combinations of ∆H, α, k that provide the zero mean bias.Wielgosz et al. (2021) showed that the mean rms of the difference between GNSS TEC maps and GNSS slant TEC was reduced from 2 to 6% when this modified mapping function was used.That choice of ∆H and α led me to a selection of the specific scaling factor k.
Although scaling and modification of the ionospheric mapping function results in a substantial reduction of the declination bias, the remaining bias is still worrisome.To mitigate the bias even further, I introduce the ad hoc correction for the ionospheric bias in the data reduction model.to observables.Here f e is the effective ionospheric path delay of a given observation.Figure 12 shows the effect of applying the de-bias correction.The bias has gone.I ran two solutions using X-band only data (XIB) and S-band only data (SIB), applied both the a priori ionospheric path delay derived from GNSS TEC maps using the modified mapping function with parameters ∆H = 56.7km,α = 0.9782, k = 0.85 and the declination bias correction.The reciprocal weights of observables were adjusted by adding in quadrature the errors of residual ionospheric errors σ rr modeled according to regression expression 11.The first and the second moments of the normalized differences in sources positions are presented in Table 3.The second moments are less than 1.0 in the X-band solution.This indicates that the ionospheric contribution is not the dominating error source in these solutions.
Unfortunately, it does not look possible to fix the deficiency of the GNSS TEC maps without re-processing GNSS observations, which is well beyond the scope of this article, but nevertheless, it is still feasible to mitigate the impact of the imperfection of GNSS TEC maps on source position estimates.Table 3.The first and second moments of position differences of trial solutions XIB and SIB with respect to the source position estimates derived from the reference solution.The a priori ionospheric contribution was computed from GNSS TEC maps using the modified mapping function with parameters ∆H = 56.7km,α = 0.9782, k = 0.85, the declination bias correction was applied, but no reduction for TEC bias adjustment has been applied.
observation, and dual-band satellite altimeters (Sekido et al. 2003;Hobiger et al. 2006;Dettmering et al. 2011;Hernández-Pajares et al. 2017;Cokrlic et al. 2018;Li et al. 2018Li et al. , 2019;;Xiang & Gao 2019;Wielgosz et al. 2021;Zhao et al. 2021;Motlaghzadeh et al. 2022).A message these publications convey is there is a reasonable agreement between GNSS TEC maps and other techniques and there are no majors problems.Therefore, large systematic errors driven by mismodeling of the ionospheric contribution derived from GNSS maps came as a surprise.Do results presented in this study contradict to prior publications?Comparison in the past was often performed using very short continuous dataset of VLBI observations, from 5 to 15 days (Hobiger et al. 2005;Dettmering et al. 2011;Etemadfard et al. 2021;Motlaghzadeh et al. 2022).The level of the agreement was characterized in terms of additive errors of GNSS TEC model.Unfortunately, when data from a short time interval are analyzed, a distinction between additive and multiplicative errors becomes blurry.Characterizing the differences in terms of biases and rms of their scatter did not turn out productive in revealing systematic errors.Moreover, an un-conscious bias of focusing a study on assessment of an agreement rather than investigation of disagreements, that were noticed and just briefly reported, diverted attention from investigating the differences in depth.
However, reading between lines of published papers, we can find pieces of evidences supporting findings in this work.The distribution of differences in VLBI vertical TEC with respect TEC derived from GNSS TEC maps presented in Figure 6 in Hobiger et al. (2006) shows a very significant skew.The distribution has a negative mean and a much greater left tail than the right tail.The VLBI vertical TEC from that study appeared less than the TEC from GNSS TEC maps, in agreement with what I have found.The authors did not investigate that stark deviation of the distribution of over one million differences from the Gaussian shape, only noting that the bias is less than 3 TECU.Considering the errors are normally distributed and additive, and considering that the derivation of TEC from VLBI does not introduce serious errors, one can expect to arrive to the Gaussian distribution of residuals.For instance, the distribution of differences in Figure 5 indeed does not show any measurable deviation from the Gaussian shape.However, if errors are multiplicative, the residuals will be non-Gaussian and their distribution will be skewed.
Satellite altimetry using Jason satellites (Vaze et al. 2010, and references therein) provides an independent way for assessment of a level of the disagreement between direct vertical TEC measurements and GNSS TEC maps.Li et al. (2018) showed that comparisons of the differences revealed significant systematic biases that depend on geomagnetic latitude.However, an attempt to characterize these additive biases in terms of the rms was not very productive.Liu et al. (2018) presented spatial distribution of the differences (Figures 6-7).That distribution strikingly reminds the average dis-tribution of TEC itself, suggesting the differences are multiplicative.A more recent paper of Dettmering & Schwatke (2022) revealed a highly significant scaling factor between GNSS TEC maps and four altimeter missions.The scaling factors varied from 0.809 to 0.919 which is very close to what has been found in my analysis of astrometric observations of extragalactic radio sources.Li et al. (2019) performed a comparison of TEC maps with Jason satellite altimetry and with ionospheric radio occultations (see the overview of this technique in Bonafoni et al. 2019).They characterized the differences as a superposition of an additive bias and a multiplicative scale factor.The scale factors defined as TEC Jason/GIM and TEC COSMIC /TEC GIM vary with time and latitude and stay in a range of 0.7-0.9 at night and 0.9-1.0 during daytime.Since Jason orbits have altitude 1,350 km and GNSS orbits have altitude about 20,000 km, Jason altimetry misses the contribution from the upper layers of the ionosphere and the plasmosphere.Yizengaw et al. (2008) analyzed TEC between Jason and GNSS satellites and found that the share of plasmosphere to TEC is on average 15% and may reach 60% when TEC is low.However, such large share of the electron density at altitudes above 1,350 km is inconsistent with the thin shell model at the altitude of 450 or 506.7 km that assumes no electron density above that height at all.
It is essential to note that analysis of Jason and COS-MIC data provides estimates of vertical TEC, while analysis if GNSS observations provides slant TEC that is converted to vertical TEC using the mapping function.The dependence of the declination biases on the mapping function found in this study suggests that a simple thin model may not be adequate.Schaer (1999) discussed the dependence of the effective ionosphere height on solar zenith angle.Xiang & Gao (2019) studied it in more details.The height of the peak electron density has annual and diurnal variations.The latter variations have an amplitude of about 100 km, being lower at daytime.They showed that the instantaneous mapping function that accounts for the height of the electron content maximum achieves 8% reduction of mapping errors.It should be mentioned that if the effective height of the ionosphere is changed, the latitude and longitude of the ionosphere piercing point for an observation with a given elevation and azimuth will change as well, which will provide an additional change in path delay.
These works strengthen argumentation in favor of that the thin shell model used for derivation of TEC maps since mid 1990s is oversimplified.The effective height of the ionosphere changes with time and latitude and therefore, a realistic mapping function should also vary with time and latitude.Omission of this complexity results in a deficiency of GNSS TEC maps that manifests in multiplicative (scale) and additive (bias) errors that varies with time and latitude.The non-linear dependence of the declination bias with declination in a form of δ(x sin δ + y cos δ) may be caused by the dependence of GNSS TEC scaling factor with latitude reported by Li et al. (2019).Radio waves from Southern Hemisphere sources observed at the array located in the Northern Hemisphere propagate through regions in the ionosphere with systematically lower latitudes than from Northern Hemisphere sources.Since a fixed mapping function is used for both computation of TEC maps and computation of slant ionospheric path delay from these maps, no modification of that fixed mapping function for data reduction of astronomical data is able to account for this kind of complexity, but as it was shown earlier, it is still possible to mitigate it.The remaining bias can be eliminated by applying the empirical de-bias correction.

Omitted refinements of the ionosphere contribution
Etemadfard et al. ( 2021) processed 60 days of VLBI data and estimated not only TEC for each site using B-splines of the 1st degree, but also TEC derivatives over longitude and latitude that were considered constant over 24 hour periods.They claim that estimating TEC partial derivatives decreases discrepancies of TEC from VLBI to GNSS TEC maps by 36%.Inspired by this results, I introduced estimation of latitude and longitude TEC gradients in a form of B-spline in addition to TEC estimation using dual-band data.I considered the weighted rms (wrms) of the differences between the parameterized model of TEC adjustment to vertical TEC from dual-band data as a metric of an improvement.I did not find a reduction of the rms greater then several percents and abandoned this approach.I should note this result does not disprove findings of Etemadfard et al. ( 2021) since they used a different metric.
Expression for the ionospheric contribution 1 is an approximation.Hawarey et al. (2005) considered the impact of the higher order of expansion on group delay, namely proportional to f −3 .They found that the maximum contribution of the 3rd degree term varied from 3 to 9 ps at 8 GHz depending on the baseline.This contribution is approximately one order of magnitude less than the contribution of residual ionospheric errors after taking into account GNSS TEC maps.It should be noted that the 3rd degree term affects both singleband and ionosphere-free linear combinations of group delays at two bands and therefore, cannot be retrieved in analysis of the differences with respect to the dual-band reference solution.
4.3.The impact of remaining errors in modeling of the ionosphere on source positions As it was shown before, our ability to model the ionospheric contribution using single-band observations is limited.After applying the declination de-bias data reduction, the declination bias is virtually eliminated.The residual ionospheric contribution causes additional random errors.In order to investigate its impact, I ran four trial solutions.I used ionosphere-free linear combinations of dual-band data and added to them the ionospheric contribution derived from VLBI data scaled to the specific frequency of the trial solution.Then I applied the data reduction for the ionosphere to this dataset of modified observables using GNSS TEC maps as if I processed single-band observations and estimated source positions.The differences in source positions from these solutions with respect to the reference solutions are interpreted as an impact of the residual ionospheric errors at different frequencies.The declination dependence of the differences in right ascensions and declinations for four solutions for frequencies 2.3 GHz (S-band), 4.3 GHz (C-band), 8.4 GHz (X-band), and 23.7 GHz (K-band) is shown in Figure 13.
These plots help to quantify additional errors in source positions that would arise if the dataset of 263 twentyfour hour experiments used for this study would have been observed at one band only.We see that estimation of TEC biases reduces the ionosphere-driven errors by a factor of 2-4.Still, even after estimation of TEC biases, the declination dependence of the ionosphere-driven errors remains.
It is known that unaccounted source structure affects source position estimates.In general, the contribution of source structure at higher frequencies is less because jets become optically thin.Therefore, there was an expectation that observing at high frequencies, such as 22-24 GHz (K-band) one may obtain more precise source positions.So far, observational evidence do not support that prediction (Lanyi et al. 2010;Charlot et al. 2020).A detailed comparison of K-band absolute astrometry versus dual-band astrometry of Karbon & Nothnagel Declination (deg) (2019) did not reveal an improvement.Therefore, it is instructive to see what is the contribution of the ionospheric errors at K-band after applying data reduction from GNSS TEC maps. Figure 14 shows additional errors in right ascensions and declinations due to residual contribution of the ionosphere.We see that additional errors in right ascension are about 0.1 mas.Declination errors of Northern Hemisphere sources are at the same level.However, these errors grow with a decrease of declination for sources in the Southern Hemisphere approximately linearly and reach 0.3 mas at declination −40 • .Therefore, the unmodeled ionospheric contribution sets the error floor in position accuracy.This error floor is not accounted in the ICRF3 catalogue, and the unaccounted ionospheric contribution is 300% greater than the noise floor at K-band adopted by Charlot et al. (2020) according to their Table 6.The potential of Kband astrometry cannot be utilized unless the accuracy of ionosphere modeling will be substantially improved.

The rms of residual atmospheric errors
I investigated how the rms of the residual ionospheric contribution for 45 VLBA baselines from the dateset of 263 twenty-four observing sessions varied with time.I computed three statistics for each observing session: and re-scaled them to 8 GHz.Here τ v and τ g are ionospheric path delays from VLBI and TEC maps respectively, a i is the adjusted TEC bias, and b i is clock bias.The statistics are shown in Figure 15.In order to improve readability, the time series were smoothed using Gaussian kernel with parameter σ = 1 year.The first statistics characterizes the impact of the total ionosphere on group delay.The second statistics characterizes the rms of the impact of the residual ionospheric errors on group delay after applying the a priori GNSS TEC model.The final

Time in years
Figure 15.The rms of the mean residual ionospheric contribution at VLBA baselines at 8 GHz for three cases: 1) no a priori ionospheric contribution is applied (upper red curve); 2) the a priori ionospheric path delay computed using GNSS TEC maps is applied (middle blue curve); and 3) the a priori ionospheric path delay computed using GNSS TEC maps are applied and TEC biases are adjusted (lower green curve). 1 TECU causes group delay 21 ps at 8 GHz.
statistics characterizes the impact of residual errors after applying the a priori GNSS TEC and adjusting TEC biases using dual-band VLBI data.

Processing geodetic data from southern and R1/R4 VLBI experiments
It would be instructive to explore to which extent the numerical values found in processing data at VLBA network are representative to experiments conducted at other VLBI networks and to which extent they are specific to VLBA.I ran two solutions XIS and XIR.Solution XIS used 63,984 X-band group delays from 36 VLBI experiments in 2016-2019 for absolute astrometry programs at core Southern Hemisphere stations hart15m, Petrov hartrao, hobart12, hobart26, kath12m, tid-bin64, wark12m, yarra12m with participation of other International VLBI Service for geodesy (IVS) stations.Solution XIR used 6,600,959 X-band group delays from 2154 geodetic VLBI experiments from IVS programs R1 and R4 (Lambert & Gontier 2006) for 2002-2022.Data from R1, R4, and Southern Hemisphere VLBI experiments are available at the NASA Crustal Dynamic Data Information System2 .In total, 36 globally distributed stations participated in experiments under R1 and R4 programs.These experiments were designed for Earth orientation parameter determination.I computed also reference solutions RIS and RIR that used ionosphere-free combinations of dualband group delays.Solution parameterization and editing were identical for pairs of XIS/RIS and RIS/RIR solutions.
I compared source positions from XIS/RIS and XIR/RIR solutions.In both cases declination systematic errors are evident.I found that scaling that makes the mean declination bias of XIS solution close to zero is k=0.78,somewhat smaller than that from solution XIB based on data from the VLBA network.Moreover, estimates of declinations from XIS solution require a different de-bias correction than declinations from XIB solution.
Figure 16 shows differences in declinations from position estimates of 638 sources with formal uncertainties < 0.5 mas from XIR solution with respect to the reference solution RIR based on data from R1/R4 networks.Comparing solutions made with different mapping function parameters k, I found that the mean declination bias vanishes when k = 0.75.The shape of the bias and its magnitude is noticeably different than the bias from declinations derived using group delays from observations at the VLBA network (see Figure 7).
The magnitude of the peak-to-peak declination bias from the solution that uses data from the R1/R4 network is 0.62 mas, while the magnitude of the declination bias from the solution that uses data from the VLBA network is 0.85 mas.The regional VLBA network stretches over 90 • in longitudes, but only within [18 • , 48 • ] in latitude, while the global R1/R4 network stretches over all longitudes and over [−43 • , +79 • ] in latitude.Observations of Southern Hemisphere sources with VLBA are possible in the southern sectors only, and the further south a source, the lower elevations it can be observed.At the same time, a given source can be observed at different azimuths and elevations at a global network.
Figure 17 shows the rms of the mean residual ionospheric contribution at R1/R4 baselines.Comparing it with a similar plot for the VLBA network (Figure 15), we see that the wrms of the ionospheric contribution at the VLBA network is a factor of 2-3 lower.This can be explained by both differences in total ionospheric path delay sensed by stations at two networks, and by the amount of the ionospheric path delay that is absorbed in estimates of parameters, such as source declination, causing biases.If a source is observed only at a given azimuth and elevation, most of the ionospheric path delay will be absorbed by estimated parameters reducing residuals.The wider spread in azimuths and elevations a given source is observed, the smaller share of the ionospheric path delay that will be absorbed by estimated parameters and propagate to residuals.
We can estimate the scaling factor directly from single-band VLBI data using least squares.I ran twelve solutions using data from the VLBA, southern, and R1/R4 networks separately and using data from all three networks altogether and applying three different flavors of the ionospheric mapping function, ∆H=0, α = 1.0; ∆H=56.7,α = 0.9782; and ∆H=150, α = 0.9782.Table 4 displays the estimates.The differences between the estimates of the scaling factors derived from observations at the R1/R4 and Southern Hemisphere networks are statistically insignificant, while the scaling factor estimates from the VLBA network are systematically lower, and these differences with respect to the R1/R4 network are statistically significant.It is also remarkable that scaling factors that minimize postfit residuals differ from those that eliminate the mean declination bias.This implies that scaling does not eliminate systematic errors due to deficiency of GNSS TEC maps but only alleviates their impact.

A combined use of dual-band and single-band data
Negligible remaining biases and the realistic assessment of errors due to the residual contribution allows us to use a mixture of dual-band and single-band group delays in a single least square solution.Source positions derived from single-band data are in general less precise than those derived from dualband data because an additional factor affects the uncertainties of group delays.These additional errors can be computed and accounted for in deriving weights of observables.An increase in uncertainties of group delay observables propagates to uncertainties of source positions.
Combined analysis of a heterogeneous datasets provides a significant advantage because all the data are used.We can fuse dual-band and single-band data in one dataset and use it for estimation of source positions.A fused dataset consists of observables of three types: 1. Dual-band ionosphere-free linear combinations of group delays.
2. Single-band delays of dual-band experiments detected at one band only.Data reduction accounts for the ionospheric path delay computed from the sum of GNSS TEC and TEC bias adjustments.
The uncertainty of the bias adjustment is added in quadrature to the group delay uncertainty for reciprocal weights of such observables.
3. Group delays from single-band experiments.Data reduction accounts for the ionospheric path delay computed from GNSS TEC and applies the de-bias correction.The uncertainty of the ionospheric model computed from the regression model is added in quadrature to the group delay uncertainty used for computation of weights of such observables.
Solving for source coordinates using a fused dataset provides estimates with minimum uncertainties and minimum correlations between right ascensions and declinations since it uses all available data provided the frequency-dependent position biases due to source structure are negligible.This condition may be violated for some strong sources with significant structure.Analysis of position offsets between single band and dual-band solutions will help to identify such sources (see, for example, Petrov et al. 2011), but these cases are encountered infrequently.

Future work
An improvement in modeling of ionospheric path delay when TEC biases are adjusted is quite impressive.However, this requires using additional information about the state of the ionosphere.Specifically, VLBI ionospheric path delays from dual-band data were used for computation of these biases.This information is not available in processing of single-band observations.Ros et al. (2000) considered the use of dual-band GPS observations from collocated receivers for analysis of radio astronomy observations.They have shown that the GPS determination of TEC from ground receivers alone without the use of GNSS TEC maps can be successfully applied to the astrometric analysis of VLBI observations.There are on-going efforts to install advanced GNSS receivers within a hundred meters of each VLBA antenna.The use of geometry-free pseudo-ranges at 1.2 and 1.5 GHz in a similar way as I used VLBI Petrov ionosphere-free group delays for adjustments of TEC biases promises a similar level of improvement.

CONCLUSIONS
Single-band group delays are affected by the contribution of the ionosphere.This contribution is noticeable at frequencies below 30 GHz and becomes a dominating error source at frequencies below 5-8 GHz.Compared with astrometric solutions based on the use of the ionosphere-free linear combinations of dual-band observables, source position estimates derived from a singleband solution are affected by additional random and systematic errors caused by mismodeling the contribution of the ionosphere to path delay.I explored two approaches to modeling the ionospheric path delay using GNSS TEC maps and assessed residual ionospheric errors.
The findings can be summarized as follows: 1.In a case when an experiment was recorded at two widely separated bands, but a fraction of sources were detected at one band only, estimation of the TEC bias in a form of an expansion over the Bspline basis using dual-band data and then applying that bias to GNSS TEC maps provides an unbiased estimates of sources positions.The stochastic model that describes residual errors of TEC bias adjustment predicts an increase of positions errors with an accuracy of 15%.No remaining systematic errors were found.This approach provides source positions with the lowest uncertainties with respect to other approaches.
2. In a case of single-band observations, path delays are computed using GNSS TEC maps.The thin spherical shell model of the ionosphere with the constant height 450 km above the mean Earth's radius causes a strong systematic bias in declinations that reaches 1 mas at 8 GHz an 12 mas at 2.3 GHz.This bias can be virtually eliminated when a) the modified ionospheric mapping function with parameters ∆H =56.7 km, α = 0.9782, k = 0.85 is used; and b) the empirical de-bias correction is applied.
4. I determined the scaling factor of TEC from GNSS maps by processing from 11,457,749 VLBI observations at different networks for 25 years: 0.834 ± 0.001.This estimate is in a good agreement with a totally independent comparison of vertical TEC determined from satellite altimetry against GNSS TEC maps.Since VLBI is sensitive to a delay incurred in the total ionosphere, interpretation of a scaling factor of TEC from altimetry and radio occultation observations as the contribution of upper layers of the ionosphere at altitudes above Jason orbit, i.e. 1,350 km, suggested by Liu et al. (2018) and Dettmering & Schwatke (2022, 2023) is not consistent with presented results.Although the plasmosphere above 1,350 km may contribute to discrepancies between Jason and GNSS TEC, systematic errors in GNSS TEC due to the adopted model of vertical distribution of electron density in a form of a thin shell at height 450 or 506.7 km is another major factor.The contribution of the plasmosphere alone cannot explain the scaling factor.
5. I have found that the scaling factor of GNSS TEC maps that provides zero mean declination bias depends on the used ionospheric mapping function.Therefore, I surmise that the established deficiency of GNSS TEC maps is caused by oversimplification of the ionospheric mapping function used for their derivation that considers the ionosphere as a thin spherical shell.The electron contents in the real ionosphere varies not only in latitude and longitude, but also in height.Diurnal variations of the effective ionosphere height at a level of 100 km, i.e. over 20%, are large enough to cause errors of the magnitude that was found in analysis of VLBI data.
6.The impact of the ionosphere on path delay depends on the Solar cycle.Modeling ionospheric path with GNSS TEC maps reduces the residuals at 8 GHz by a factor of 2 during the Solar maximum and only by 10% during the solar minimum.Estimation of the TEC bias reduces ionospheric errors further by a factor of 2 regardless of the phase of the Solar maximum.
7. The impact of the ionosphere on source position errors can be modeled with an accuracy of 15%.It remains noticeable even at frequencies as high as 22-24 GHz (K-band).In particular, the ionospheric errors even after applying data reduction based on GNSS TEC maps with the modified mapping function and the de-bias correction exceed 0.1 mas.Declination errors of Southern Hemisphere sources observed with VLBA are in a range of 0.1 to 0.3 mas.An assertion that K-band astrometry is able to provide results more precise than 0.1 mas is not true at the current state of our ability to model ionospheric path delay.Considering on-going efforts to install advanced GNSS receivers in the close vicinity of VLBA stations and other radio telescopes, the situation may change in the future.
This study lays the foundation of the single-band absolute astrometry.Dual-band astrometric observations still provide the best accuracy.The use of singe-band data with the procedure of data reduction and weighting described above allows us to get unbiased positions with known added errors.
6. ACKNOWLEDGMENTS This work was done with datasets RDV, RV, CN, UF001, UG003, BL122, BL166, BP133, BP138, GC073, BC204, BG219, and V17 collected with VLBA instrument of the NRAO and available at https://archive.nrao.edu/archive.The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.The author acknowledges use of the VLBA under the USNO's time allocation for some datasets.This work made use of the Swinburne University of Technology software correlator, developed as part of the Australian Major National Research Facilities Programme and operated under license.

Figure 1 .
Figure 1.Adjustment to the ionosphere path delay bias at 8.4 GHz with respect to the path delay derived from GNSS TEC maps at mk-vlba station from processing dual-band observations on April 22, 2015.The bias estimate corresponds to the averaged positions of the ionosphere piercing points within a given time interval.1 TECU causes group delay 21 ps at 8.4 GHz.

Figure 2 .
Figure 2. Empirical distribution of the normalized differences of the ionosphere path delay computed from the GNSS TEC maps adjusted for clock and TEC biases (green dots).The normal distribution with σ = 1 is shown as a reference (solid blue line).

Figure 4 .
Figure 4. Dependence of the rms of residual ionospheric path delay derived from GNSS TEC maps σgr on the rms of the total ionospheric path delay from these maps σgt.No adjustment to TEC has been applied.Path delay is computed for the reference frequency 8 GHz.The blue smooth line shows the regression model in a form of a B-spline that fits the data.

Figure 5 .
Figure 5.The distribution of the normalized differences of the ionospheric path delay computed from the GNSS TEC maps against VLBI ionospheric path delay with clock biases subtracted (green dots).The normal distribution with σ = 1 (solid blue line) is shown as a reference.

Figure 6 .
Figure 6.The distribution of second moment estimates of the normalized differences of ionospheric path delays derived from VLBI dual-band observations and GNSS TEC maps among individual observing sessions.The narrow green curve shows the statistics of the normalized residuals with TEC biases adjusted and the wide blue curve show the statistics of normalized residuals without TEC adjustment.

Figure 7 .
Figure7.Differences in declinations from the X-band solution XIT using data at the VLBA network with the data reduction for the ionosphere using ionospheric path delay from GNSS TEC maps applied with respect to the dual-band reference solution.The thick blue line shows the differences smoothed with the Gaussian filter.
;Hobiger et al. (2006);Dettmering et al. (2011);Motlaghzadeh et al. (2022) claimed a good agreement between TEC derived from VLBI and GNSS.And indeed, the plot of ionospheric contributions in zenith direction from VLBI after removal clock biases against the ionospheric contributions from GNSS (Figure8) shows no peculiarities and fits the straight line τ vi = −3.9ps + 1.06 • τ gi .Although the residuals of this dependence look random, they still cause systematic errors in source positions.

Figure 8 .
Figure 8. Dependence of VLBI ionospheric group delay at 8 GHz against the ionospheric group delay from GNSS.The blue straight line is the least square fit of this dependence.

Figure 10 .Figure 11 .
Figure 10.Smoothed declination biases from four solution with respect to the dual-band reference solution using Xband observables and the ionospheric mapping function with parameters ∆H = 56.7 km and α = 0.9782 with different scaling factors k = 0.7, 0.8, 0.9, and1.0.
I smoothed the biases D(δ) with the Gaussian filter with σ = 8 • -these smoothed biases are shown in Figures 9-11 -and expanded them over the basis of B-splines of the 3rd degree with 13 equidistant knots in the range of [−45 • , 90 • ].The expansion coefficients are presented in the appendix.I added correction

Figure 13 .
Figure13.The rms of the differences in right ascensions (upper row) and declinations (lower row) when single band group delays are used with respect to the positions from the dual-band reference solution.A priori ionospheric path delay using GNSS TEC with the modified mapping function (∆H = 56.7 km, α = 0.9782, k = 0.85) was used adjusting TEC biases from dual-band observations (right) and without adjusting biases but with applying the de-bias correction from expression 13 (left).The upper blue band shows the differences for S-band, next red line shows the differences for C-band, next green line shows the differences for X-band, and the bottom purple line shows the differences for K-band.

Figure 14 .
Figure 14.An increase in source position errors in right ascension (Left) and declination (Right) due to the residual ionospheric contribution at K-band (22 GHz) after applying data reduction based on GNSS TEC maps with the modified mapping function (∆H = 56.7 km, α = 0.9782, k = 0.85)

Figure 16 .
Figure16.Differences in declinations from the X-band solution XIR with respect to the dual reference solution RIR using data at the R1/R4 VLBI, similar to Figure7.Default TEC GNSS mapping function (∆H = 0, α = 1.0, k = 1.00) was used that comparison.

Figure 17 .
Figure 17.Similar as Figure 15 but derived from analysis of observations at R1/R4 VLBI network.Statistics computed for one year intervals.
pleasure to thank Yuri Y. Kovalev, Urs Hugentobler, Robert Heinkelmann, Minghui Xu, and Denise Dettmering for discussions of presented results.