Inverse Multview I: Multi-Calibrator inverse phase referencing for Microarcsecond VLBI Astrometry

Very Long Baseline Interferometry (VLBI) astrometry is a well established technique for achieving $\pm10~\mu$as parallax accuracies at frequencies well above 10~GHz. At lower frequencies, uncompensated interferometer delays associated with the ionosphere play the dominant role in limiting the astrometric accuracy. Multiview is a novel VLBI calibration method, which uses observations of multiple quasars to accurately model and remove time-variable, directional-dependent changes to the interferometer delay. Here we extend the Multiview technique by phase referencing data to the target source ("inverse Multiview") and test its performance. Multiple observations with a four-antenna VLBI array operating at 8.3~GHz show single-epoch astrometric accuracies near $20~\mu$as for target-reference quasar separations up to about 7 degrees. This represents an improvement in astrometric accuracy by up to an order of magnitude compared to standard phase referencing.


INTRODUCTION
High accuracy astrometry at radio frequencies provides fundamental information for many fields of astronomy and astrophysics (Reid & Honma 2014). Very Long Baseline Interferometry has provided trigonometric parallaxes with accuracies of 10 µas or better for masers associated with massive young stars throughout the Milky Way and approaching this accuracy for evolved stars (e.g. Reid et al. 2019;VERA Collaboration et al. 2020). In addition, parallaxes for X-ray binaries (e.g. Miller-Jones et al. 2021) and pulsars (e.g. Deller et al. 2019) have been critical to characterizing these sources and, in some cases, using them to test General Relativity. Finally, measurements of extra-galactic proper motions approaching ∼ 1 µas y −1 accuracy have been accomplished (e.g Brunthaler et al. 2005). The most accurate radio frequency astrometry has been achieved at observing frequencies above ∼ 10 GHz where ionospheric delay errors, which scale with observing frequency as ν −2 , are typically small. However, below this frequency, these delay errors become the dominant source of astrometric error, and improving astrometric accuracy at frequencies below ∼ 10 GHz requires new approaches to measure and remove ionospheric delays (Rioja & Dodson 2020) The electron density in the ionosphere varies with a strong diurnal signature above a location on Earth. Since the dominant source of ionization is solar ultraviolet radiation, the highest electron densities track the sub-solar point and this can produce strong gradients in the electron distribution both in Geographic longitude and latitude. Models of the total electron content (TEC) as a function of time and location on the Earth can be used to partially correct for propagation delays through the ionosphere, but these models can have uncertainties of 20% or more (Walker & Chatterjee 1999). The result is that residual phase-delays of ∼ 0.1 nsec can be present in radio interferometer data at 8.3 GHz. Importantly, these phase-delays can have significant and long-lived gradients ("ionospheric wedges") across ∼ 10 • of sky, which can degrade the accuracy of relative position measurements between nearby sources (e.g. Reid et al. 2017).
Multiview is a novel approach for calibration of Very Long Baseline Interferometric (VLBI) observations designed to achieve the highest possible astrometric accuracy, particularly at frequencies below about 10 GHz where ionospheric effects are the dominant source of position errors. The core idea of Multiview is that observations of multiple calibrators, which surround the target on the sky, allow for removal of directional, timevariable residual phase-delays from the calibrators to the target, using 2D spatial and temporal interpolation in the visibility domain. Initial trials of Multiview calibration, called cluster-cluster phase referencing (Rioja et al. 1997(Rioja et al. , 2002, involved simultaneous observations of a target and multiple calibrators by utilising multiple telescopes at a single site. This method showed promise in removing residual ionospheric delays which plague low-frequency astrometry. However, the availability of multiple telescopes at many sites is extremely limited. The next iteration of Multiview used only a single telescope at each site and source switching (Rioja et al. 2017). Observations at 1.6 GHz were conducted using the Very Long Baseline Array (VLBA) and structured with three calibrators (C 1 , C 2 , C 3 ) surrounding a target OH maser (T ) with the observing sequence: We refer to this approach as direct Multiview. A critical requirement is that all sources have to be observed within the atmospheric coherence time. Tests showed that excellent results were achieved when one entire sequence was completed in 5 minutes at 1.6 GHz (Rioja et al. 2017). Reid et al. (2017) investigate fitting and removing a single 'positional' gradient to the measured positions of multiple calibrators after these had been phasereferenced to the target and imaged at each epoch. As this fitting was done in the image domain, we refer to this approach as 'image-based' Multiview (imMV), and it has been shown to improve astrometric accuracy at 6.7 GHz (Sakai et al. 2019;Zhang et al. 2019).
Standard phase-referencing (PR) for VLBI observations involves 'nodding' all telescopes between a calibrator and a target, measuring the phase on the calibrator and transferring it to the target (Alef & Porcas 1986;Beasley & Conway 1995). Inverse phase referencing (iPR) is commonly used for astrometry of astrophysical maser sources when the target is strong and the calibrator may be weak. Here the phase of the target is transferred to the calibrator, and ultimately the measured offset position of the calibrator is used to infer that of the target. In PR/iPR (or imMV) one will often observe a sequence of N calibrators as T, C 1 , T, C 1 , T, . . . T, C N , T, C N , T . . .
(2) thereby allowing either iPR or PR to be used (e.g. if the target is weaker than expected). The main benefit of the iPR technique for Multiview applications is that it only requires that adjacent observations of the target (e.g. T, C1, T ) are spaced by less than the coherence time of the atmosphere. This is especially valuable at frequencies around 7 GHz, where coherence times are typically a factor of 5 shorter than at 1.4 GHz.
In this paper, we introduce and test "inverse Multiview" (iMV), a combination of iPR and Multiview. In order to test iMV we selected strong 'target' quasars surrounded by compact extragalactic radio sources acting as calibrators, and observed them in multiple sessions over 4 months. By using the target as the phasereference, we tested the positional repeatability after application of normal iPR, iMV and imMV, and compared the results. Also, by using calibrators that surrounded the target at increasingly larger separations, we could evaluate how far separated from the target one can use calibrators without violating the assumption of planar ionospheric wedges.

Sources & Observations
We selected three groupings of quasars at right ascension 6, 13 and 19 hours with declination < 0 • from the catalogue (rfc 2019b) of Petrov et al. (2019). All quasars had a catalogued 8.3 GHz unresolved flux density ≥ 100 mJy (Table 1) and where possible, were chosen to have little to no structure. Five or six calibrator quasars were selected to be distributed in a thin 'ring' around the target quasar ( Figure 1) with mean radii of θ sep = 2.8, 6.7 and 7.5 • respectively.
Observations were constructed of 3 types of blocks: fringe-finder blocks (FFBs) for fringe-alignment and clock determination in the correlator, and preliminary electronic delay and phase calibration (i.e., manual phase calibration); geodetic-like blocks (Geoblock) for advanced delay calibrations; and iMV blocks. These blocks were scheduled in the following repeated sequence: Geoblock, FFB, iMV block.
FFBs are comprised of 3-4 scans on strong (≥ 1 Jy) quasars over a combined duration (on-source and slew) of < 15 mins. These quasars are optimally located close to the target to minimise slewing time and at > 45 • elevation at all telescopes.
Geoblocks are short geodesy-style (Heinkelmann 2013) observation periods that consist of 10 − 15 ICRF2 quasars (Fey et al. 1991) with sub-milliarcsecond accurate positions spread in elevation at each telescope site. These 30 min blocks are scheduled every 3 hr and allow for correction of post-correlation residual tropospheric and clock delays at each telescope (Honma et al. 2008;Reid et al. 2009a;Reid & Honma 2014).
The iMV blocks are a modified version of conventional inverse phase-referencing nodding sections (Equation 2) and are placed between the calibration blocks. The iMV blocks form the primary data for the experiment. For a reference target (T ) with N calibrators (C 1...N ), iMV blocks have the following sequence: T, C 1 , T, C 2 , T, . . . , C N , T, C 1 , T, C 2 , T, . . . , C N , T (3) where target scans bracket sequentially different calibrator scans. This allows all individual calibrators scans to be phase-referenced to the target in under the atmospheric coherence time (∼ 4 min at 8.3 GHz).
The iMV blocks were ∼ 150 min duration, with each full cycle (C 1 → C N ) taking on average about 15 mins with individual on-source times of 50 sec, and the remaining time reserved for slewing. Therefore each 21 hr observation contained seven calibration blocks and six iMV blocks, allowing two iMV blocks for each ring cluster spanning seven hours, which kept source elevations above 30 • .
In this way, we conducted four 21 hour observations of the three ring clusters on 2019 February 16, March 17, April 13 and May 4. At these epochs the angular separation from the Sun was 120, 100, 85 and 70 • for the smallest 2.8 • ring; 45, 72, 100 and 120 • for the 6.7 • ring and; 120, 150, 175 and 160 • for the largest 7.5 • ring. The angular separation from the Sun is expected to influence the size of the ionospheric gradients encountered.
Data were recorded at 1024 Mbps in right circular polarization covering frequencies between 8200 and 8456 MHz, with Nyquist sampling and 2-bits per sample. Baseband data were correlated using the DiFX-2 software correlator (Deller et al. 2011) at 0.5 MHz spectral resolution. Raw correlated FITS files can be provided upon request.

Preliminary Calibrations
Data were analyzed in AIPS (Greisen 1990(Greisen , 2003 using standard VLBI tools and with the assistance of the python wrapper software ParselTongue/Obit (Kettenis et al. 2006). Correlated FITS files were loaded into AIPS using the task FITLD and data observed during off-source periods or windstows were flagged with task UVFLG. This generally included only a few seconds at the start of some scans, with the exception of a ∼ 1 hr period at the end of the second epoch where Ceduna 30m was windstowed. Task ACCOR was used to correct amplitude in the cross-correlation spectra arising from Table 1. Target and calibrator positions, separations and flux densities. Columns: Target (1) and calibrator (2) names, correlated positions in right ascension (3) and declination (4), calibrator-target offset in right ascension/East-West (5) and declination/North-South (6), total separation (7), catalogue 8.3 GHz flux density (8), average synthesised image integrated intensity (9). digitizer sampler threshold errors and antenna system temperatures were applied with task ANTAB. Updated Earth Orientation Parameters (Seidelmann 1982) were downloaded 1 and applied with the task CLCOR/EOP, and Total Electron Content (TEC) maps based on global positioning system data were downloaded 2 and applied with task VLBATECR. Data observed in the geoblocks were separated and processed as follows: 1. The task FRING was run on a single FFB scan to fit single-band delay and phase (whilst zeroing the delay-rate); then the task CLCAL was used to apply this solution to all geoblock quasars. This 1 gdc.cddis.eosdis.nasa.gov/vlbi/gsfc/ancillary/solve apriori/ 2 gdc.cddis.eosdis.nasa.gov/gnss/products/ionex/ removes the time-invariable delays between intermediate frequency bands.
2. FRING was then used to fit multi-band delays and rates on all geoblock quasars. These data, along with the antenna elevations, were used to determine tropospheric zenith delays, clock offsets and drift rates at each telescope using an external, DELZN-equivalent program (described in Reid et al. 2009b). These values were applied to the iMV data with the task CLCOR/ATMO.
Data for each ring group in the iMV blocks were processed as follows: 3. A single FFB scan was used to remove the singleband delays and phases (as with step 1 on the geoblocks).   4. The target quasar was fit for phase and rate using the task FRING, and the solutions were applied to both the target quasar and its ring calibrators using the task CLCAL.
5. The calibrators were imaged using the task IMAGR, and their peak emissions were fitted with the task JMFIT.
6. The above process was repeated for all four epochs and the positions of the ring calibrators were shifted to the mean measured offset with the task CLCOR/ANTP. This minimizes possible phase wrapping (see Section 3.1). After the position corrections were applied to the data, a second iteration of the same calibration process was undertaken (from step 3, skipping this step).
7. The data for the ring calibrators were averaged in frequency using the task SPLAT and a fringe fit for phase was performed using the task CALIB. The task TBOUT was used to print out the CALIB phase solutions so that they could be used for iMV fitting.
8. The target quasars were averaged in frequency with task SPLIT.
At VLBI resolutions, many sources are somewhat resolved and this can lead to additional phase variations (structure phase). For sources that are partially resolved, preliminary self calibration may be used (between steps 3 and 4) to correct for the structure phase of each calibrator and target separately. We found that the targets and calibrators selected for these observations were sufficiently compact as to not require this correction.

Phase Wraps
The initial correlation used the catalogue quasar positions and antenna location which have a reported accuracies of 0.3 mas and 1 cm respectively, sufficient to yield residual interferometric phase variations of 120 • over a track (at 8.3 GHz with a 3500 km baseline). With this effect alone, phases should therefore remain in the same −180 • ,180 • wrap over the track. However, the residual tropospheric and ionospheric delays generally introduce phases which can vary many times this amount over a track.
Due to the abundance of relatively strong (> 300 mJy) quasars in our data, we could measure the residual delay present after preliminary calibration (but before phasereferencing). We performed a fringe fit for multi-band delay (with task FRING) on target quasars and FFB quasars, and this revealed that there were residual de-lays in our data that had variations between 3 to 10 cm (Figure 3), equivalent to 300 to 1000 • of phase. After phase referencing to the central target, the ring sources would then be expected to have differential phase errors of up to φ ≈ 1000∆θ • , where ∆θ is the angular separation in radians. For ∆θ ≈ 0.1 radians, we expect these atmospheric phase errors to approach 100 • , and the total amount of phase variation to possibly exceed 180 • . This can lead to some phase-wrapping (Rioja & Dodson 2020), which needs to be corrected for before the iMV fitting process.
In order to correct phase-wrapping, we minimised the absolute phase difference between consecutive scans on the same calibrator quasar, separated by approximately ∆t = 15 min (the average ring cluster duty-cycle), by adding or subtracting 360 • as needed to keep the phasedifference, We used the phase measured at the middle of each track as a reference, since this is when the average antenna elevation is largest and we expect the phase differences due to the unmodeled residual atmospheric delays to be at a minimum. For the data presented here, only a few baselines required correction for phase wrapping, rarely for the 2.8 • ring, and more commonly for the larger separation clusters (1-2 times per track on average over all baselines).

Inverse Multiview Fitting
We modelled the phase-screen over a quasar grouping as a 2D plane or 'wedge' taking the form: where ∆α i cos δ T and ∆δ i are the angular offset from the target position for the i th calibrator (in degees on the sky). The subscript jk indicates an interferometer baseline. The measured phase on the i th calibrator for the jk baseline is given by φ i,jk in degrees of phase, and φ T,jk is the phase at the target position. The parameters A jk , B jk are the phase slopes in units of degrees of phase per degree of offset in the East-West and North-South directions, respectively. Least-squares fitting was used to determine the phase slopes and target source phase in a sliding 15 min window interpolated to target quasar scan times ( Figure 4). As there will be a separate phase screen per antenna and all antennas have been referenced to a common antenna, we only solved Equation 5 for the n − 1 baselines to that reference antenna (in this case n = 4). The reference antenna was always chosen to be either Ceduna 30m or Hobart 26m.
By comparing the residual phases and least-squares fit parameters over time, we were able to estimate if a quasar had an uncorrected phase wrap ambiguity; typically these occurred only once or twice per track. For example, in Figure 4 top panel, only the two red data points near 14:00 UT (which originally were at −120 • ) and the two purple points at 16:30 UT (which originally were at −170 • ) required correction for this effect.
Finally, we flagged data where the mean absolute residuals exceeded 1 radian, as these indicated that the assumption of a planar phase gradient was not sufficient to explain the phase screen (for example, the break in solid black line at 16:15 UT in Figure 4 top panel). In general, only a few minutes of data was lost per track due to this, often near the beginning and/or end and accounting for less than 4% of the total iMV block baseline hours observed.
The target phase solutions, φ T,jk (solid black line, top panel in Figure 4), determined from the iMV process were loaded into AIPS using the task TBIN to the SPLIT target data (from step 8 in Section 2.3) as a solution (SN) table and then applied to the target data when  imaging. Note that the phase-referencing during pre-iMV calibration had been applied to the target quasar (in addition to the ring quasars), leaving it with essentially zero phase. This process effectively transfers all residual phase errors for the target quasar equally to all ring quasars, removing phase shifts due to both atmospheric and position errors. However, the iMV phases, φ T,jk , which are defined at the position of the target quasar, will reflect phase shifts only associated with the target quasar's possible position shift, and the iMV ionospheric phase slopes are essentially side products.
The target quasar was imaged using the task IMAGR and the peak emission in the images were fit with the task JMFIT. Measured target quasar positions are given in Table 2. Since quasars should have essentially zero (sub-microarcsecond) motions over our observations, we used the scatter in the sky positions (x, y) (in µas) of the target quasars over the four epochs as the estimate for a single-epoch positional accuracy. In the East-West direction, the single-epoch positional accuracy (σ x in µas) was estimated with the standard deviation: where x is the mean of x. The uncertainty in the singleepoch position accuracy (SE σ ) was estimated with: (Rao 1973) where N is the number of epochs. The single-epoch position accuracy and uncertainty in the North-South direction (σ y ,SE σy ) were estimated using analogous equations.

Comparison with iPR and imaged-based Multiview
In order to compare iMV with iPR accuracy, we replicated the analysis process that would be undertaken if each calibrator was separately referenced to the target (e.g. Equation 2). In such a case, each calibrator would be separately imaged, have its offset from the phase centre measured at each epoch, and then the (negative of the) offsets would be used over time to infer target apparent motion. The results from each calibrator could then be combined. Therefore, we took the scatter in offset for each calibrator over time (Equation 6), then took the average of these scatters for each ring (black stars and circles in Figure 5).
In the comparison with imaged-based Multiview, the calibrators are first separately referenced to the target, imaged and have offsets measured at each epoch (as with iPR). Next, these offsets are used to construct an 'artificial quasar' at the location of the target. The motion of this artificial quasar should then mirror that of the target. Therefore we followed the "Method-2" from Reid et al. (2017) -fitting a tilted plane (with slopes S x , S y and constants C x , C y ) to each of the x and y offsets of the calibrator quasars at each epoch as a function of their angular separation from the target: where ∆α i cos δ T and ∆δ i are the respective East-West and North-South angular separations for the i th calibrator quasar from the target at position (α T , δ T ) (Table 1). The slopes S x , S y should reflect the average delay gradient at each epoch, and the constants C x , C y give the likely position shift of the targets at each epoch. It is from these position shifts that we determined the imMV scatter (blue triangles in Figure 5).

Thermal Uncertainty
The thermal noise in an astrometric image limits how precisely the position of a feature in said image can be determined. We refer to this accuracy limit as the thermal uncertainty (σ pos,th ), and estimate its magnitude as σ pos,th ≈ 0.5 θ B /SNR (Reid et al. 1988). Here θ B is the synthesized beam (set by the ratio of the observing wavelength λ to the baseline length |B|) and SNR is the signal-to-noise ratio for the feature in the images.
For the final astrometric images of the three targets, the average synthesized beam was θ B = 2.3 × 1.3 mas at position angle 60 • (hereafter using the geometric mean θ B = 1.73 mas), and the average SNR was 175, 110, and 70 for the 2.8, 6.7, and 7.5 • separation targets respectively. This gave a respective estimate for the thermal uncertainty as σ pos,th = 5, 8 and 12 µas for each of the three separations.

RESULTS & DISCUSSION
For iMV the single epoch accuracy is ±20 µas in both coordinates for calibrator separations from the target of at least 6.7 • . As expected, for iPR the accuracies are separation dependent, growing from about ±100 µas at 2.8 • to about ±300 µas at 6.7 • . The accuracies of imMV are between those of iPR and iMV, since imMV is, essentially, an observing-track average of a changing phase screen.
There is one possible outlying measurement in our data in Table 2 at epoch 3 in the 7.5 • ring (J1336-0829) for the East-West (x) direction. Removing this measurement would yield a standard deviation of the remaining three points of 45 ± 22 µas, more in keeping with our other results. One possible reason for having an outlier is an unresolved phase wrap ambiguity, which is expected to be more likely for larger target-calibrator separations. Another potential source of problems for iMV calibration is tropospheric water-vapour fluctuations, which would not be expected to give a planar "phase wedge" over the larger angular separation of some of the rings we observed, and in such cases would not be properly corrected by iMV. Experience with standard phase referencing at 22 GHz suggests significant degradation in image quality if a target is separated from a calibrator by more than about 3 • on the sky. For a given (non-dispersive) delay error, interferometer phases scale linearly with observing frequency, and the 3 • "limit" at 22 GHz scales to about 8 • at our observing frequency of 8.3 GHz. Even when including this possible outlier point, iMV still gives a much better singleepoch uncertainty than either imMV at 250 µas or iPR at 325 µas for the same coordinate and separation.
These tests of the iMV technique show a dramatic improvement in accuracy for VLBI astrometry at 8.3 GHz compared to standard phase-referencing ( Figure 5). For example, Reid et al. (2017) found evidence for singleepoch position errors of ≈ 100 µas per degree of separation between a target and a calibrator at an observing frequency of 6.7 GHz. Scaled to an observing frequency of 8.3 GHz, this error estimate is close to the dot-dashed line in Figure 5 inferred from our iPR results.
Comparing the iMV results with those that would have been obtained from imMV imply that fitting in the visibility domain is preferable to fitting in the image domain, and gives a factor of 2-10 improvement in accuracy ( Figure 5). This can be likely explained by the ability to measure and correct the for 'spatially dynamic' components of the phase-screen that changes on hour-scales (e.g. A jk (t), Figure 4) whereas imMV effectively averages this effect over the whole track and whole array.
The iMV astrometric results are consistent with a single-epoch accuracy near ±20 µas, which is approaching the thermal uncertainty estimates from Section 3.4. Furthermore, the iMV results appear largely separation independent, as would be expected if the accuracy were dominated by thermal processes (Rioja & Dodson 2020). This is a dramatic improvement over standard phase referencing at this frequency, where typical residual ionospheric path-delays of 5 cm would translate to roughly ±100 µas (for a baseline length 3500 km and targetcalibrator separation 2 • ).
With accuracies near 20 µas for iMV at 8.3 GHz using calibration sources with separations up to ≈ 7 • separations, we have matched or exceeded some of the best results obtained at 22 GHz (e.g. Reid et al. 2019;VERA Collaboration et al. 2020) where ionospheric effects are about seven times smaller. As the direct Multiview analysis would be expected to have much the same performance as iMV, we expect that future observations will reveal that it achieves comparable astrometric errors.

CONCLUDING REMARKS
To date, the highest astrometric accuracies have been obtained from inverse phase referencing at 22 GHz using the VLBA, a homogeneous array of 10 antennas with a maximum baseline length of around 8500 km. These benchmark results are obtained only for small (< 2 • ) angular separations between the target and reference source, coupled with accurate "geodetic-block" modeling, where non-dispersive delays, owing mostly to water vapor in the atmosphere, dominate. Residual dispersive delays at 22 GHz, after applying global TEC models, are generally small (∼ 1 cm path delay). Here we have shown that through application of the new technique of inverse Multiview it is possible to obtain comparable astrometric accuracy for small arrays of heterogeneous antennas with shorter baseline lengths, at lower frequencies and in the presence of much larger residual delays. Multiview methods will be central to achieving ultraprecision astrometry on the next generation of instruments, such as ngVLA and SKA (for further reading, see sect. 7.2 in Rioja & Dodson 2020).
The Bar and Spiral Structure Legacy (BeSSeL) Survey and the VERA project have measured parallaxes to ≈ 250 massive young stars which display water or methanol maser emission. This has resulted in a partial map of the spiral structure of the Milky Way. We have begun parallax observations of 6.7 GHz methanol masers associated with massive young stars that cannot be seen with telescopes in the northern hemisphere in order to complete this map. Astrometric results at this frequency using iMV, which we will explore in paper II, have achieved single-epoch accuracies approaching 20 µas, consistent with the results of the tests documented here. This should allow parallax accuracies of ≈ 10 µas, which translate to 10% distance uncertainty at 10 kpc from the Sun.

DATA AND CODE AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author. The programmes and scripts used for data reduction and analysis are available from https://github.com/lucasjord. Table 2. Measured positional offsets for the three target quasars over the four epochs and the determined positional accuracy. Columns: Average target-calibrator separation (1), target name (2), direction corresponding to data (3), measured position over four epochs (4-7), average positional offset from phase centre (8)