A broadband VLBI system using transportable stations for geodesy and metrology: an alternative approach to the VGOS concept

Abstract

We have developed a broadband VLBI (very long baseline interferometry) system inspired by the concept of the VLBI Global Observing System (VGOS). The new broadband VLBI system was implemented in the Kashima 34 m antenna and in two transportable stations utilizing 2.4 m diameter antennas. The transportable stations have been developed as a tool for intercontinental frequency comparison but are equally useful for geodesy. To enable practical use of such small VLBI stations in intercontinental VLBI, we have developed the procedure of node-hub style VLBI: In joint observation with a large, high sensitivity ‘hub’ antenna, the closure delay relation provides a virtual delay observable between ‘node’ stations. This overcomes the limited sensitivity of the small diameter node antennas, while error sources associated with large diameter antennas, such as gravitational deformation and delay changes in necessarily long signal cables, are eliminated. We show that this scheme does not result in an increased sensitivity to radio source structure if one side of the baseline triangle is kept short. We have performed VLBI experiments utilizing this approach over both short range and intercontinental distance. This article describes the system components, signal processing procedure, experiment, and results in terms of baseline repeatability. Our measurements reveal signatures of structure effects in the correlation amplitude of several of the observed radio sources. We present a model of the frequency-dependent source size for 1928+738 derived from correlation amplitude data observed in four frequency bands.

1 Introduction

Historically, very long baseline interferometry (VLBI) has utilized the S-band (2.3 GHz) and X-band (8.4 GHz), which have been assigned for deep space communication and for geodetic VLBI observations. The International VLBI Service for Geodesy and Astrometry (IVS) is now promoting the VLBI Global Observing System (VGOS) (Petrachenko et al. 2012; Niell et al. 2018) as the next generation geodetic VLBI system, which is targeting 1 mm geodetic precision by observing four 1 GHz bands spread over the 2–14 GHz frequency range. The concept of VGOS was developed to improve the estimation of the atmospheric delay contribution by fast sampling using quickly slewing antennas. To enable measurement of the delay observable with sufficient accuracy, broadband observation, high rate data acquisition, and dual polarization observation were included in the specification. Since the precision of group delay measurements is inversely proportional to the effective bandwidth (Rogers 1970) [see later Eq. (6)], a ten times wider frequency range enables a one order of magnitude improvement in delay precision with the same signal-to-noise ratio (SNR). The SNR itself improves with bandwidth B as

$$\begin{aligned} \mathrm{SNR} \propto D_1\sqrt{\frac{\eta _1}{T_{\mathrm{sys,1}}}} D_2 \sqrt{\frac{\eta _2}{T_{\mathrm{sys,2}}}} S \sqrt{T B}, \end{aligned}$$

(1)

where $D_{{i}}$, $\eta _{{i}}$, and $T_{\mathrm{sys,}i}$ are diameter, aperture efficiency, and system temperature of stations ($i = 1, 2$). S is the flux of the radio source, and T is the integration time. Broadband observation thus provides a dual improvement for geodetic VLBI: Sensitivity improves with high rate data acquisition, and delay precision improves with large effective bandwidth. These features open the possibility of utilizing transportable small diameter VLBI stations for geodesy and distant frequency transfer.

The idea of using transportable VLBI stations has been investigated since the 1970s (e.g., Clark 1979; Davidson and Trask 1985). Mobile VLBI stations with 9 m, 5 m, and 3.5 m diameter antennas were deployed in the CDP (Crustal Dynamics Project) for geodetic measurements in North America and Europe (Ryan et al. 1993). The Geodetic Survey Institute of Japan (GSI) has used several sizes (2.4 m, 3.8 m, and 5 m) of mobile stations for geodetic VLBI measurements in the period from 1986 to 1995 (Takashima and Ishihara 2008). The Transportable Integrated Geodetic Observatory (TIGO) (Hase et al. 2000) is a fundamental station equipped with a 6 m offset-parabola antenna for VLBI in addition to SLR (Satellite Laser Ranging), GNSS (Global Navigation Satellite System), and DORIS (Doppler Orbitography and Radiopositioning Integrated by Satellite) terminals. It has contributed to IVS observation in the southern hemisphere at Concepción in Chile from 2002 to 2014. As discussed by Plank et al. (2015), the inhomogeneous network and overall shortage of VLBI stations in the southern hemisphere are important issues for the geodetic VLBI community. Transportable VLBI stations may help improve the global distribution of geodetic sites.

In addition to geodesy, transportable VLBI can also enable precise frequency comparison over intercontinental distances. Driven by the rapid development of high-accuracy optical frequency standards, a re-definition of the SI second as the unit of time is being discussed in the metrological community (Riehle 2015). Confirming the performance of frequency standards at different locations requires accurate methods for frequency comparison. Frequency transfer by optical fiber-link reaches uncertainties below 10$^{-18}$ (Lopez et al. 2012; Predehl et al. 2012; Calonico et al. 2014), but no trans-oceanic links have been realized so far. Advanced methods of two-way satellite time and frequency transfer (TWSTFT) using the signal’s carrier phase have the potential to reach instabilities on the order of $10^{-17}$ (Fujieda et al. 2016). The two-way exchange of signals allows for a direct calibration of atmospheric excess path length, but the availability of suitable communication satellites and the need to license radio transmissions limit its applicability. Space geodetic techniques are not subject to these limitations.

Because of their ubiquitous character, Precise Point Positioning (PPP) methods using the carrier phase of Global Navigation Satellite Systems (GNSS) can be used as control techniques in this work. Integer-ambiguity techniques (IPPP) (Petit et al. 2015) can also reach instabilities of order $10^{-17}$, but need continuous measurements over longer averaging durations.

VLBI application for time and frequency transfer has been investigated since the 1970s (e.g., Counselman et al. 1977; Saburi 1978; Clark et al. 1979), and an uncertainty of $1.5\times 10^{-15}$ for a time period of 1 day has been reported in an earlier study (Rieck et al. 2012). In a minimal implementation of such a VLBI observation, a pair of radio telescopes performs a differential phase measurement referenced to stable frequency standards at each end. The baseline vector and clock difference between the two stations are obtained by analyzing the observables in terms of signal arrival time difference from distant radio sources. The first derivative of the clock offset with regard to time, termed the clock rate, is also determined during VLBI data analysis.

For atomic clocks used as time references in VLBI stations, which have a relative frequency instability well below $10^{-12}$, the VLBI-determined clock rate is equivalent to the difference in the coordinate rate of the two clocks (see Appendix A). From this comparison of the VLBI station reference clocks, high accuracy frequency standards [such as cold atom clocks or ion clocks (e.g., Ludlow et al. 2015)] operated near the VLBI stations can be compared remotely.

Stimulated by the VGOS concept, we developed the broadband VLBI system named ‘GALA-V’ (Sekido et al. 2016) with the aim to enable intercontinental precise frequency comparison. Transportable broadband VLBI stations can be flexibly installed at desired locations, such as the metrology institutes operating next-generation frequency standards. The 34 m antenna located at Kashima Space Technology Center provides the high sensitivity used to boost the SNR of observations [Eq. (1)]. We introduce an algorithm for this node-hub style (NHS) VLBI, using transportable stations with a high-sensitivity antenna, in Sect. 2.8.

Fig. 1

Kashima 34 m antenna (left), and the NINJA feed in its receiver room (right)

Table 1 Receiver parameters and slew rate performance for Kas34, MBL1, and MBL2

Fig. 2

Block diagram of the receiving signal chain for the data acquisition of Kas34. Dual (V & H) linear polarization signals are transferred from the receiver room of the antenna to the observation room through one single mode fiber (SMF) using a wavelength division multiplexer (Mux) and demultiplexer (DeMux). For efficient data acquisition, equalizers (EQ) are inserted to compensate for power attenuation in the signal path. A notch filter at 6.55 GHz avoids potential saturation by RFI. For the transportable VLBI stations, only the path of (V) linear polarization is implemented. E/O electric–optical conversion; O/E optical/electric conversion, HPF high-pass filter, BPF band-pass filter

This article first describes the observation system, the survey of the radio frequency interference (RFI) environment along with resulting mitigation measures, and the signal processing in Sect. 2. VLBI observations conducted over short and intercontinental baselines are reported in Sect. 3. Section 4 discusses the results of our experiments, including radio source structure observed in correlation amplitude. Section 5 gives concluding remarks.

2 The broadband VLBI observation system ‘GALA-V’

2.1 Broadband receiver

The broadband feed, which is the antenna component converting incoming electromagnetic waves to electric current, is one of the key devices of the receiving system. Almost all VGOS stations use a feed based on the Quadruple-Ridged Flared Horn (Akgiray et al. 2013) developed by the California Institute of Technology, which has a wide beam size around $120^\circ$. This makes it difficult to use in existing high-gain Cassegrain-focus radio telescopes, where the physical distance from focal point to the sub-reflector is not easy to change and the pre-set beam size is usually narrow (e.g., $34^\circ$ for the Kashima 34 m antenna). To use existing Cassegrain antennas for broadband observations, NICT (the National Institute of Information and Communications Technology) developed the NINJA feed system (Ujihara et al. 2018). We installed this on the Kashima 34 m antenna (Kas34) and two 2.4 m diameter antennas (hereafter referred to as MBL1 and MBL2). Our receiver employs a room-temperature low noise amplifier (LNA, Low Noise Factory LNF-LNR4_14B) rather than a cryogenic LNA to save cost and time for development, installation, and maintenance. The NINJA feed allows simultaneous observation of V and H linear polarization. Photographs of the Kashima 34 m antenna and the broadband NINJA feed are displayed in Fig. 1.

Signals of both polarizations are acquired in Kas34, while the 2.4 m antennas only record the V-polarization signal. Although it is preferable to acquire a complete set of polarization data to achieve higher SNR, the restriction to single polarization for only the transportable stations is a compromise between practical considerations and minimizing the loss due to polarization misalignment. Besides better cost–performance, it allows for a quick turnaround of experiments where data recorded at a remote site have to be transferred to a correlation center (e.g., Kashima Space Technology Center). Although performed over high-speed research networks (e.g., Internet2, GEANT, GARR, and JGN), the present data transfer rate is only about half of the data acquisition rate for a single polarization. Acquiring dual polarizations would result in a twofold increase in turnaround time. Thus, we opted to employ single polarization observation at the remote sites and synthesize the correlation data to form an emulated polarization-aligned dataset (Sect. 2.5).

2.2 Broadband VLBI stations

2.2.1 Kashima 34 m antenna

The Kashima 34 m antenna is a standard VLBI station and was built in 1988. Thus, it does not have the fast slew rate of a VGOS antenna system, which has 12 $^{\circ }$/s (Petrachenko et al. 2012). Table 1 gives an overview of antenna parameters. Broadband dual polarization observation is enabled by the NINJA feed. The two polarization signals covering the frequency range of 3–14 GHz are transferred to the observation room by an RF signal transmission system (LTA-40-M, PD-30-M; Optilab Co. Ltd.) via optical fiber. The signal path from the feed to the sampler is displayed in Fig. 2. A notch filter (6R1-6538-100S11, Reactel Co., Ltd.) was installed to avoid saturation by RFI at 6.55 GHz. The system equivalent flux density (SEFD) represents the noise level limiting the sensitivity of the antenna. The SEFD of the Kashima 34 m antenna (Fig. 3, Table 1) shows noise below the requirement of the VGOS specification for frequencies below 7.5 GHz, and only slightly worse in the high frequency region. Due to the almost identical sensitivities, a standard VGOS station is equally suitable to take the place of the Kashima 34 m antenna in node-hub style VLBI observation (see Sect. 2.8).

2.2.2 Transportable 2.4 m antennas

The two small antennas (MBL1 and MBL2) of our project were originally built for the project MARBLE (Ishii et al. 2010) with an S/X receiver system. The MARBLE antenna was designed with transportability in mind. It can easily be assembled or disassembled without heavy machinery. Keeping these easy-installation features, the antennas have been upgraded by changing the main reflector size from 1.6 to 2.4 m diameter. Station coordinates of radio telescopes and antenna performance parameters are listed in Tables  1 and 2 . The designations Kas34, MBL1, and MBL2 refer to the antennas throughout this article. To refer to a VLBI station, we use the conventional name KASHIM34 for the Kashima 34 m station and MARBLE2 for MBL2 installed at NICT (Koganei). The station names MARBLE1 and MARBLE1M refer to the installations of MBL1 at NMIJ (National Metrology Institute of Japan, Tsukuba) and INAF (Istituto Nazionale di Astrofisica, Medicina), respectively.

Fig. 3

Frequency-dependent SEFD of the Kashima 34 m antenna with broadband NINJA feed system for V and H polarizations

Table 2 Station names and nominal coordinates

The signal chain from the receiver system to the data acquisition system (DAS) is almost the same as for the Kashima 34 m antenna (Fig. 2) except that only one polarization is implemented in the transportable stations. The radio frequency (RF) devices in the signal path including filter, amplifier, and electro–optical converter (E/O) are installed in a temperature controlled ‘receiver box’ on the small antennas, visible in Figs. 10 and 14. After E/O conversion, the received radio frequency signal is transmitted from the receiver box to the observation room via a single-mode optical fiber cable.

2.3 RF-direct sampling

The VGOS specification proposes acquiring four radio frequency bands with 1 GHz bandwidth allocated in the 2–14 GHz frequency range (Petrachenko et al. 2012). The NASA proof-of-concept (PoC) system (Niell et al. 2018) has realized this with a flexible up-down frequency converter that selects the observing band for conversion to intermediate frequency (IF) with 1 GHz width. Data acquisition is performed by sampling the IF signal and extracting multiple channels of 32 MHz width by digital filtering and frequency conversion implemented in the sampler.

Our DAS achieves the same goal in a different way. The received RF signal is amplified and separated into lower (0–8192 MHz) and upper ( > 8192 MHz ) signals by a power divider and filters (see Fig. 2). Each of these signals is then directly digitized at 3-bit quantization by a high-speed sampler K6/GALAS (Takefuji et al. 2012; Sekido 2015) (product by Elecs Industry¹) with 16,384 MHz sampling rate. The sampler then extracts the four desired signals with 1024 MHz bandwidth per polarization by digital filters implemented in an FPGA (field-programmable gate array). The location of the frequency bands is flexibly selected by 1 kHz steps. We refer to this as RF-direct sampling (RFDS), and it was initially tested for simultaneous data acquisition of S/X band with the 11 m antenna at Kashima and the Tsukuba 32 m antenna (Takefuji et al. 2012).

Reducing data acquisition to four channels per polarization is a unique characteristic of our observing system. The smaller number of signal paths reduces relative phase variations caused by temperature or physical stress. The overall reduction of analog components likewise improves phase stability. The RFDS technique brings large benefits for phase calibration, polarization synthesis, and wideband bandwidth synthesis in the signal processing (see Sects. 2.5, 2.6 and 2.7). The reduced number of channels and wider frequency band require fewer steps of digital filtering. We estimate that FPGA resources required to extract four 1024 MHz wide channels from data sampled at 16,384 MHz are reduced by factors of 6.4, 3.9, 2.3, and 3.1, respectively, for lookup table, flip-flop, RAM, and DSP compared to 64 channels of 32 MHz width from data sampled at 2048 MHz as in the case of the VGOS PoC system.

The K6/GALAS sampler has four 10GBASE-SR interfaces, and it is capable of providing a total data rate of 16,384 Mbps via eight data streams of VDIF/VTP/UDP/IP Ethernet packets. Each stream conveys 2048 Msps of 1-bit quantized data for each observing band. Two off-the-shelf Linux-OS personal computers (PCs) serve as data recording disk servers for 16,384 Mbps of dual polarization data at KASHIM34 station. The 8192 Mbps data rate for a single polarization requires only one PC at the transportable stations. Each PC is equipped with two 10GBASE-SR network interfaces and a RAID disk system.

2.4 Radio frequency interference (RFI) and mitigation

Broadband observations are more likely to be affected by incoming RFI than traditional narrow band observations. Such RFI caused by artificial radio emission is usually several orders of magnitude stronger than natural radio sources. Thus, it does not only cause loss of signal at the frequency of RFI but can render the entire radio telescope useless by saturating the receiver system.

In 2012, preceding the design of the receiver system, we performed RFI surveys at Kashima, Koganei, and Tsukuba, where initial installation and short-baseline test experiments had been planned (Sect. 3.1). These surveys employed a transportable receiver system composed of a double ridged broadband horn antenna (Schwarzbeck BBHA-9120D) and room temperature LNA (B&Z Technologies BZP118UD1, $F = 0.1$–18 GHz). Strong RFI at 0.8 GHz and 2.6 GHz made a high-pass filter (HPF: RLC F-100 ${F}_{\mathrm{c}}$ = 3.0 GHz) necessary to avoid saturation of the LNA. The system temperature of the receiving system was around 500 K for the frequency range of 2–18 GHz. The measurements were made by max-hold recording over 30 s with a spectrum analyzer (Rohde & Schwarz FSV-30) and directing the BBHA-9120D antenna (beam size: 60$^{\circ }$ for H-plane, 90$^{\circ }$ for E-plane) toward each of the four cardinal directions (elevation angle: 0$^{\circ }$, azimuthal angle: 0$^{\circ }$, 90$^{\circ }$, 180$^{\circ }$, 270$^{\circ }$). Figure 4 shows the spectra recorded at Kashima and Koganei. A significant amount of RFI is present at NICT’s Koganei campus, which is located inside a Tokyo residential area where radio waves from cellular base stations and Wi-Fi transmitters are abundant. Experiments with telecommunications technology further contribute.

Fig. 4

Radio frequency power spectra obtained by maxhold measurement for 30 s during the RFI survey at Kashima (top) and Koganei (bottom) performed in 2012

Based on the survey results, we integrated a filter bank behind the LNA in the signal path of the MBL2 antenna to avoid saturation of the receiver system. RFI signals below 3.2 GHz are attenuated by a high-pass filter (${F}_{\mathrm{c}}$ = 3.0 GHz) in addition to the inherent frequency characteristics of the feed system. The pass-band frequencies of the filters were carefully selected as a compromise between rejecting strong RFI and keeping the sensitive frequency range as wide as possible for observations. Figure 5 shows the results of pointing the MBL2 antenna toward the sun and the empty sky. The measured ratio of detected power confirms that the receiver system has sensitivity for all pass-band frequencies without saturation.

Fig. 5

Signal power received from the sky (blue) and the sun (red) for MBL2 in Koganei (Tokyo)

2.5 Polarization synthesis

The orientation of linear polarization differs by the parallactic angle between observing stations separated by long distance, leading to a reduced cross-correlation amplitude. Dual polarization data acquisition to obtain full Stokes parameters has been demonstrated (Martí-Vidal et al. 2016) to be the best way for recovery of polarization alignment. We follow a different approach in avoiding the loss of correlation due to parallactic angle. Let us use H and V for data of each polarization acquired at the high sensitivity antenna and $V'$ for the V-polarization obtained at the small transportable antenna. Correlation processing is performed for each combination $HV'$ and $VV'$ by using the GICO3 software correlator (Kimura 2008) as the first step. Then, correlation data of the aligned polarization $\widehat{{VV}'}$ are composed by a linear combination of $HV'$ and $VV'$ data as

$$\begin{aligned} \widehat{{VV}'} = {VV}' \cos (\delta \chi (t)) - {HV}' \sin (\delta \chi (t))e^{j(\omega \tau _0+ \phi _0)}, \end{aligned}$$

(2)

where $\delta \chi (t)$ is the parallactic angle difference between the observation sites as a function of time. $\tau _0$ and $\phi _0$ account for the delay and phase difference between the two sets of correlation coefficients. They are caused by the difference of signal paths between V and H polarizations at the hub antenna (Kas34). These parameters are determined to maximize the SNR of the synthesized correlation products $\widehat{{VV}'}$. Figure 6 illustrates the spatial configuration of the two correlations $HV'$ and $VV'$ and the parallactic angle difference $\delta \chi (t)$. The phase difference $\phi _0$ between the two polarizations is assumed to remain constant throughout the session thanks to the RF-direct sampling technique described in Sect. 2.3. Polarization synthesis is only applied for the intercontinental Medicina–Kashima baseline (Sect. 3.2). For the shorter baselines Kashima–Koganei and Kashima–Tsukuba, the parallactic angle is almost identical between stations, and $VV'$ can be regarded as polarization-aligned correlation data. We note that Eq. (2) works only for unpolarized radio sources.

Fig. 6

Emulated polarization-aligned cross-correlation data $\widehat{{VV}'}$ are computed by linear combination of the cross-correlations for $HV'$ and $VV'$, taking into account the parallactic angle difference $\delta \chi (t)$

2.6 Phase calibration with a distant radio source

The phase and amplitude of the signal received by the antenna are affected by propagating through microwave components such as antenna feed, LNA, filters, mixers, and the sampler. Thus, VLBI systems employ compensation techniques to recover the signal for accurate delay measurements. In the conventional S/X VLBI system as well as the proof-of-concept VGOS system, a phase calibration (Pcal) signal, based on a comb-tone pulse train and injected from the frontend of the antenna, is used as a reference and is recorded together with the observed signal. The Pcal phase recovered from recorded data then calibrates the phase–frequency response of the signal transmission line (Whitney et al. 1976).

At present, no sufficiently stable Pcal system is available for our broadband observation system. We therefore use a natural radio source signal for phase calibration (hereafter referred to as ‘RSPcal’) instead of a local Pcal signal. A reference radio source with sufficient SNR is selected and observed as a ‘reference scan.’ A ‘scan’ refers to the VLBI observation of an individual radio source over a certain length of time. The phase spectrum of the correlation coefficient (cross-spectrum) of the reference scan is then stored as phase calibration data of that baseline. The procedure itself is very similar to the manual Pcal, which applies constant phase offsets to the correlation coefficients of each channel so that the fine delay-resolution function (Whitney et al. 1976) is synthesized coherently. While our RSPcal method may be considered a manual Pcal in the broadest sense, there is a distinction: Where manual Pcal is typically applied per-station, RSPcal is baseline based. It thus contains the full differential phase response between the two stations, including not only the signal transmission line, but also the geometrical propagation delay, tropospheric delay, and ionospheric dispersive delay to that reference radio source at the epoch of the reference scan. By assuming that the instrumental phase response is invariant with time, the phase response of the transmission line can be calibrated with this data.

For the broadband NINJA feed, a frequency-dependent phase characteristic may result from the design using a composition of multi-mode horn and dielectric radio lens to form a narrow beam over a wide frequency range. Thus, the capability to calibrate the phase characteristics including the feed system is an advantage of RSPcal, as standard Pcal does not cover the phase response of the feed. While the reliance on an instrumental phase response that is constant over time and antenna motion is a potential drawback of the RSPcal, the RFDS technique mitigates this: (1) Digitization of data without frequency conversion avoids using analog devices and the associated phase changes; (2) high-speed sampling and digital filtering enable acquisition of multiple bands with identical signal path. These features reduce relative phase variation with respect to frequency. We also paid specific attention to the temperature control of the environment of the sampler and the signal cables. As described before, RF components are installed inside a temperature-controlled box. We took care to protect the signal cables from mechanical stress and have not detected any evidence of characteristics changes by antenna motion so far. For the experiments at Medicina campus, the long optical fiber cable between the antenna and observation room was placed in an underground trench. The delay error caused by thermal expansion of the optical fiber is discussed in Sect. 4.1.

While we typically perform RSPcal for each observation session, the response of the system is very stable. Figure 13 in Sect. 3.1 shows results of successfully tracking the clock difference UTC(NMIJ)–UTC(NICT) over a month, obtained with no further RSPcal procedure after the initial calibration. This requires that components in the signal path remain unchanged and that the sampler’s 1PPS synchronization is preserved.

2.7 Extraction of ionospheric TEC and broadband group delay

Group delay is defined as a linear phase gradient over frequency. In VLBI, that corresponds to the phase slope of the correlation coefficients in the frequency domain. This cross-spectrum phase of the raw observation data is affected by not only non-dispersive delay, but also dispersive ionospheric delay. The contribution of the differential ionospheric excess phase delay between two VLBI stations is expressed as

$$\begin{aligned} \delta \phi _{\mathrm{ion}} = -A \cdot {\mathrm {TEC}}/f, \end{aligned}$$

(3)

where ${\mathrm {TEC}}$ is the difference of electron column density in the two line of sights from the observing stations to the common source, f is the observing radio frequency, and $A = e^2/4 \pi \varepsilon _0 m_e c$ is a physical constant composed of elementary charge (e), permittivity in vacuum ($\varepsilon _0$), electron mass ($m_e$), and speed of light (c). The conventional signal processing of S/X bands (Clark et al. 1985) calibrates differential ionospheric excess delay between propagation paths by a linear combination of group delays obtained for S-band and X-band. This is valid when the fractional bandwidth is narrow (12% for current S/X geodetic VLBI). However, in case of broadband VLBI with 3–14 GHz frequency range, ionospheric excess phase can no longer be approximated by a simple linear phase gradient, and the full dispersive phase delay has to be modeled using Eq. (3). A wideband bandwidth synthesis (WBWS) algorithm (Kondo and Takefuji 2016, 2019) developed for our project simultaneously estimates relative ionospheric total electron content (${\mathrm {TEC}}$) and broadband group delay. We give a brief review of the mathematical description of the WBWS process and RSPcal calibration in Appendix B. By assuming invariance of instrumental phase, the group delay and ionospheric TEC are derived as relative quantities with respect to the reference scan.

$$\begin{aligned}&{\Delta } \tau ^i = \tau ^i - \tau ^{{\mathrm {ref}}} \nonumber \\&= \tau ^i_{\mathrm{geo}} + \tau ^i_{\mathrm{atm}} + \cdots - \tau ^{{\mathrm {ref}}} \end{aligned}$$

(4)

$$\begin{aligned} {\Delta } {\mathrm {TEC}}^i= & {} {\mathrm {TEC}}^i - {\mathrm {TEC}}^{{\mathrm {ref}}} \end{aligned}$$

(5)

Equation (4) states that the observable differs from normal delay only by a constant offset. Thus, it can be processed with standard VLBI analysis software with proper treatment of the offset. One has to note that the obtained ${\Delta } {\mathrm {TEC}}^i$ is a relative value in Eq. (5).

For the group delay, precision is inversely proportional to the effective bandwidth EBW and SNR (e.g., Rogers 1970; Clark et al. 1985):

$$\begin{aligned} \sigma _\tau = \frac{1}{2 \pi \cdot \mathrm{SNR} \cdot {\mathrm {EBW}} } \end{aligned}$$

(6)

EBW is defined by the root-mean-square (RMS) separation of the observing frequencies from the mean as ${\mathrm {EBW}} = \sqrt{ \frac{1}{N} \sum _i^N {( f_i - {\bar{f}})^2 }}$, where ${\bar{f}} = 1/N \sum _{i}^{N} f_i$. A test experiment of broadband VLBI between the Kashima 34 m and the Ishioka 13 m antenna, a VGOS station operated by the Geospatial Information Authority of Japan (GSI), demonstrated sub-picosecond delay precision (Sekido 2015).

The delay-resolution function (DRF) of the bandwidth synthesis is determined by Fourier transformation of the correlation coefficients allocated on the frequency axis. Thus, the choice of the radio frequency bands is an important factor in the group delay determination. The DRF consists of a periodic sequence of peaks in the time domain with a sinc-function envelope that has a temporal width reciprocal to the frequency width of each single channel. The per-channel bandwidth of 1024 MHz in our broadband system provides an envelope as narrow as 1 ns. The structure resulting from the multi-channel synthesis has a period that can be determined by the reciprocal of the greatest common divisor (GCD) of the intervals of the frequency array. A periodicity shorter than the realistic range of signal arrival times introduces ambiguity. It is a remarkable property of broadband measurements that channels can be arranged to provide unambiguous identification of the group delay (see Fig. 7). See Sects. 3.1 and 3.2 for the frequency arrays used for short baseline sessions.

Fig. 7

Delay-resolution function for the array of frequency bands centered on 6.0, 8.5, 10.4, and 13.3 GHz with 1024 MHz width each. This array was used for the intercontinental broadband VLBI between Medicina–Kashima–Koganei discussed in Sect. 3.2. The plot focuses on the positive side of the symmetric function to provide greater detail. The inset shows a magnification of the behavior near 0 delay

2.8 Node-hub style VLBI observation

When applying VLBI for long-distance frequency comparisons, small ‘node’ stations are convenient for installation near the frequency standards, but it may then be impossible to directly obtain the delay observable due to limited SNR. Here, we show how joint measurements with a high sensitivity ‘hub’ station overcome this problem by evaluating the closure relation for the full measurement network. We have named this ‘Node Hub style’ (NHS) VLBI.

Fig. 8

Closure delay defined by the difference in arrival time of an identical wavefront from a radio source at all three stations. A plane wavefront represents the idealized case of a point-like radio source at infinite distance

This observation scheme was previously unattractive because low SNR resulted in a large delay error. Now that a large effective bandwidth allows obtaining sufficient delay precision with minimum SNR, NHS VLBI becomes an appealing option. Besides flexibility and low cost, a small VLBI station also avoids gravitational antenna deformation and can help reduce temperature-dependent signal transmission cable length change. To obtain a delay observable in the NHS VLBI scheme, we consider the closure delay relation schematically depicted in Fig. 8. Let us suppose three stations are observing a common radio source, and an identical plane wavefront of a radio signal arrives at a hub station (R) and two node stations (A and B) at epochs of $t_{\mathrm{R}}$, $t_{\mathrm{A}}$, and $t_{\mathrm{B}}$, respectively. The propagation part of VLBI delay observable $\tau _{\mathrm{XY}}$ is defined as the time of arrival at station Y minus the time of arrival at station X, where the latter serves as reference epoch of this delay data (Petit and Luzum 2010). This corresponds to the propagation delay as $\tau _{\mathrm{XY}}^{\mathrm{prop}} = \tau _{\mathrm{XY}}^{\mathrm{geo}}+\tau _{\mathrm{XY}}^{\mathrm{atm}}$, the sum of the signal propagation delay due to the geometry of the XY baseline and the propagation in the atmosphere including the ionosphere. They are represented as follows:

$$\begin{aligned} \begin{array}{rcl} \tau _{\mathrm{RA}}^{\mathrm{prop}}(t_{\mathrm{R}}) &{}= &{} t_{\mathrm{A}} - t_{\mathrm{R}} \\ \tau _{\mathrm{AB}}^{\mathrm{prop}}(t_{\mathrm{A}}) &{}= &{} t_{\mathrm{B}} - t_{\mathrm{A}} \\ \tau _{\mathrm{BR}}^{\mathrm{prop}}(t_{\mathrm{B}}) &{}= &{} t_{\mathrm{R}} - t_{\mathrm{B}} = -\tau _{\mathrm{RB}}^{\mathrm{prop}} (t_{\mathrm{R}}) , \end{array} \end{aligned}$$

(7)

where the delay $\tau _{\mathrm{XY}}$ is defined at the reference epoch $t_{\mathrm{X}}$ of signal arrival at station X. Summing the three equations gives

$$\begin{aligned} \tau _{\mathrm{AB}}^{\mathrm{prop}} (t_{\mathrm{A}}) = \tau _{\mathrm{RB}}^{\mathrm{prop}} (t_{\mathrm{R}}) - \tau _{\mathrm{RA}}^{\mathrm{prop}} (t_{\mathrm{R}}). \end{aligned}$$

(8)

Note that this gives the delay at epoch $t_{\mathrm{A}}$, which is shifted from $t_{\mathrm{R}}$ by a fractional second as given by $\tau _{\mathrm{RA}}^{\mathrm{prop}}(t_{\mathrm{R}})$. We expand the delay of the AB baseline in a Taylor series around epoch $t_{\mathrm{R}}$ as

$$\begin{aligned} \tau _{\mathrm{AB}}^{\mathrm{prop}}(t_{\mathrm{A}}) = \tau _{\mathrm{AB}}^{\mathrm{prop}}(t_{\mathrm{R}}) + \sum _{i=1}^\infty \frac{{\mathrm {d}}^i }{{\mathrm {d}} t^i} \tau _{\mathrm{AB}} \cdot (\tau _{\mathrm{RA}}^{\mathrm{prop}}) ^i . \end{aligned}$$

(9)

Using Eqs.(8) and (9), we then find

$$\begin{aligned} \tau _{\mathrm{AB}}^{\mathrm{prop}}(t_{\mathrm{R}})&= \tau _{\mathrm{RB}}^{\mathrm{prop}}(t_{\mathrm{R}}) - \tau _{\mathrm{RA}}^{\mathrm{prop}}(t_{\mathrm{R}}) \nonumber \\&\quad - \frac{ {\mathrm {d}}}{{\mathrm {d}} t} \tau _{\mathrm{AB}} \cdot \tau _{\mathrm{RA}}^{\mathrm{prop}} - \frac{{\mathrm { d}}^2}{{\mathrm { d}} t^2} \tau _{\mathrm{AB}} \cdot \big (\tau _{\mathrm{RA}}^{\mathrm{prop}}\big )^2 - \cdots \end{aligned}$$

(10)

Note that a theoretical a priori delay value $\tau _{\mathrm{RA}}^{\mathrm{prop}}$ and its derivative must be used here instead of the observed delay $\tau _{\mathrm{RA}}^{\mathrm{obs}}$, which may contain a clock offset of a few $\upmu$s. A good model such as Calc (Gordon et al. 2017) provides a priori delays with an uncertainty of less than $1 \times 10^{-8}$ s. The magnitude of the second-order term of Eq. (10) is less than $3 \times 10^{-14}$ s on any ground-based baseline, and higher-order correction terms can be safely neglected.

As is conventional in VLBI, we choose epoch $t_{\mathrm{R}}$ to represent an integer second. Hereafter, the time argument ($t_{\mathrm{R}}$) is omitted for simplicity, and all the delay values are assumed to be evaluated at $t_{\mathrm{R}}$. The delay observable of the RA baseline is composed of the contribution of the propagation delay including geometry and atmosphere ($\tau _{\mathrm{RA}}^{\mathrm{prop}}$), station-dependent delays such as clock and instrumental delays ($\tau _{\mathrm{R}}^{\mathrm{stn}}$, $\tau _{\mathrm{A}}^\mathrm{stn}$), and the delay caused by radio source structure ($\tau _{\mathrm{RA}}^{\mathrm{str}}$). Then,

(11)

The radio source structure effect of the group delay is the derivative of the visibility phase with respect to the observing radio frequency (Charlot 1990). The difference of delays between the RA and corresponding RB baselines as given by Eq. (11) is then

$$\begin{aligned} \tau _{\mathrm{RB}}^{\mathrm{obs}} - \tau _{\mathrm{RA}}^{\mathrm{obs}} = \tau _{\mathrm{RB}}^{\mathrm{prop}} -\tau _{\mathrm{RA}}^{\mathrm{prop}} + \tau _{\mathrm{B}}^{\mathrm{stn}} - \tau _{\mathrm{A}}^{\mathrm{stn}} +\tau _{\mathrm{RB}}^{\mathrm{str}} -\tau _{\mathrm{RA}}^{\mathrm{str}}. \end{aligned}$$

(12)

By eliminating ($\tau _{\mathrm{RA}}^{\mathrm{prop}} , \tau _{\mathrm{RB}}^{\mathrm{prop}}$ ) with Eqs. (10) and (12), we define a virtual delay observable $\tau _{\mathrm{AB}}^\mathrm{obs*}$ as

(13)

(14)

The delay observable of the AB baseline corresponding to the form of Eq. (11) is

$$\begin{aligned} \tau _{\mathrm{AB}}^{\mathrm{obs}} = \tau _{\mathrm{AB}}^{\mathrm{prop}} + \tau _{\mathrm{B}}^{\mathrm{stn}} -\tau _{\mathrm{A}}^{\mathrm{stn}} + \tau _{\mathrm{AB}}^{\mathrm{str}}. \end{aligned}$$

(15)

The difference between the virtual observable $\tau _{\mathrm{AB}}^\mathrm{obs*}$ and the true observable of the AB baseline ($\tau _{\mathrm{AB}}^{\mathrm{obs}}$) is

$$\begin{aligned} \tau _{\mathrm{AB}}^\mathrm{obs *} - \tau _{\mathrm{AB}}^{\mathrm{obs}} = \tau _{\mathrm{RB}}^{\mathrm{str}} -\tau _{\mathrm{RA}}^{\mathrm{str}} -\tau _{\mathrm{AB}}^{\mathrm{str}}. \end{aligned}$$

(16)

Only the radio source structure effect differentiates the expected true VLBI delay from the virtual delay observable of the AB baseline. This virtual delay observable $\tau _{\mathrm{AB}}^\mathrm{obs*}$ computed by Eq. (13) is the basis of the NHS VLBI scheme.

When the network of the three stations is composed of two long baselines and one short baseline, the radio source structure effect is often negligible for the short baseline ($\tau _{\mathrm{RB}}^{\mathrm{str}} \cong 0$). While it still has significant magnitude on the long baselines, these contributions largely cancel because of their opposite directions: $\tau _{\mathrm{RA}}^{\mathrm{str}} \cong -\tau _{\mathrm{AB}}^{\mathrm{str}}$. The difference between $\tau _{\mathrm{AB}}^\mathrm{obs*}$ and $\tau _{\mathrm{AB}}^{\mathrm{obs}}$ then becomes small. Placing the hub station close to one of the node stations also ensures a maximal mutually visible sky for the three stations. This is the case in our experiment over intercontinental baselines, where RA is 8810 km, and AB is 8780 km, but RB is only 109 km.

An advantage of the NHS scheme is that the delay contribution of the hub station expressed as $\tau _{\mathrm{R}}^{\mathrm{stn}}$ is eliminated in the process of linear combination from Eqs. (11) to (12). Systematic changes of delay induced in signal cables and antenna structure by temperature or mechanical stress are thus canceled.

Another approach to the NHS scheme is placing the hub station at an intermediate place between the two nodes. As the structure effect is a nonlinear function of the baseline vector, the combined structure effects on the two node-hub baselines of intermediate lengths might be of smaller magnitude than that of the long baseline between the nodes. The potential of the NHS VLBI scheme in reducing structure effects needs further investigation.

Although the virtual delay precision is affected by the errors of both observables ($\tau _{\mathrm{RB}}^{\mathrm{obs}}$, $\tau _{\mathrm{RA}}^{\mathrm{obs}}$), this does not lead to a significant degradation in precision since the overall error is typically dominated by the long baseline, where correlated flux is lower. Figure 9 compares the error of the delay observable over long and short baselines [obtained with Eq. (6)] to the error of the virtual delay observable. Over all measurements, the average delay errors obtained as the standard deviation of the post-fit residuals were 3.1 ps and 5.2 ps for short (Kashima–Koganei) and long (Medicina–Kashima) baselines, and 6.4 ps for the synthesized Medicina–Koganei baseline, respectively. The integration time of this data was 30 s or 60 s as described in Sect. 3.2.1. Obtaining VLBI observations with delay errors of a few ps between sensitivity-limited 2.4 m diameter antennas with such short observation time is a notable achievement.

Fig. 9

Histograms of group delay error magnitude for long baseline Kas34-MBL1, short baseline Kas34-MBL2 and NHS virtual delay of MBL1-MBL2 for the session on 25 Dec. 2018. The root-mean-square values for this session are shown by vertical lines and in good agreement with the overall averages listed in the text

2.9 VLBI data analysis

The VLBI data analysis software package Calc/Solve (Gordon et al. 2017) developed by NASA/GSFC (National Aeronautics and Space Administration/Goddard Space Flight Center) was used for analyzing the virtual delay observables between the 2.4 m antenna pair. Atmospheric delay is one of the most significant contributions to microwave delay measurements in space geodesy and needs to be corrected as precisely as possible. In the case of a single intercontinental VLBI baseline, the azimuth angle of the antennas was biased in a particular direction by the mutually visible sources. This made it difficult to estimate anisotropic properties of the atmosphere, which was parameterized by atmospheric gradients. We therefore used the dataset of the Vienna Mapping Function (VMF3) (Landskron and Böhm 2018; re3data.org 2019) to compute the a priori atmospheric delays in the form of an elevation-dependent mapping function:

$$\begin{aligned} \tau _{\mathrm{atm}}&= \tau ^{\mathrm{z}}_{\mathrm{h}} \cdot mf _{\mathrm{h}}(\epsilon ) + \tau ^{\mathrm{z}}_{\mathrm{w}} \cdot mf _{\mathrm{w}} (\epsilon ) \nonumber \\&\quad + mf _{\mathrm{g}} (\epsilon ) \cdot [G_{\mathrm{n}}\cdot \cos (\theta ) + G_{\mathrm{e}}\cdot \sin (\theta ) ] , \end{aligned}$$

(17)

where $\theta$ is the azimuth angle measured clockwise from north, $\epsilon$ is the elevation angle from horizon to zenith, $\tau ^{\mathrm{z}}_{\mathrm{h}}$ and $\tau ^{\mathrm{z}}_{\mathrm{w}}$ are hydrostatic dry and wet delay components, $G_{\mathrm{n}}$, $G_{\mathrm{e}}$ are atmospheric gradient coefficients, and $mf _{\mathrm{h}}(\epsilon )$, $mf _{\mathrm{w}}(\epsilon )$, $mf _{\mathrm{g}}(\epsilon )$ are mapping functions for dry, wet, and gradient components, respectively. Mapping functions are defined as the ratio of delay in the slanted direction relative to zenith direction. The gradient coefficients and good overall accuracy of the VMF3 are generated by ray tracing through the numerical weather model of the ECMWF (European Centre for Medium-Range Weather Forecasts). An intercontinental VLBI baseline requires more observations of sources at low elevation angle, making the atmospheric gradient component particularly important. The coefficients of the VMF3 dataset are provided at 6-h intervals, and we applied linear interpolation for evaluation at an arbitrary epoch. The virtual delays of the baseline between the two 2.4 m antennas were stored in a Mark-III VLBI database and analyzed with the Calc/Solve software suite according to the following steps:

• 1. The geometrical a priori delay model was computed by the Calc software and stored in the database.

• 2. To include the VMF3 data in the analysis, we externally calculated an a priori atmospheric delay for the line of sight from each antenna to the observed radio source as described by Eq. (17): The hydrostatic zenith delay $\tau ^{\mathrm{z}}_{\mathrm{h}}$ was computed from ground pressure data and station location using the model of Saastamoinen (1972). All other terms of the atmospheric model, such as the wet zenith delay $\tau _{\mathrm{w}}^{\mathrm{z}}$, the mapping functions $mf _{\mathrm{h}} (\epsilon )$ and $mf _{\mathrm{w}} (\epsilon )$, and the gradient coefficients $G_{\mathrm{n}}$ and $G_{\mathrm{e}}$, were given by the dataset of the VMF3. All the components (dry, wet, gradient) were summed up and stored in the Mark-III database.

• 3. The Solve software subtracted geometrical and atmospheric delays from the observed delay data for each scan and performed least-squares estimation. The delays of the VMF3 (step 2) replace the built-in atmospheric models of Solve here. The delays after subtraction are then interpreted in terms of station coordinates, clock parameters, and atmospheric model. (a) The remaining atmospheric components are estimated as zenith delay parameters ($\delta \tau _{\mathrm{wh}}^\mathrm{z}$) via piecewise linear functions with 60 min interval for both stations. (b) Regarding station coordinates, the MBL2 antenna at Koganei was fixed for both NMIJ–NICT and INAF–NICT sessions (Sects. 3.1 and 3.2), because it was attached to a solid foundation on the roof of NICT’s building in Koganei. The location of MBL1 remained open to estimation by VLBI, as the antenna did not have a solid attachment at NMIJ. In the INAF–NICT sessions, station coordinates were determined by VLBI after the installation in Medicina. (c) Clock parameter estimation differed slightly between sessions of NMIJ–NICT and INAF–NICT. The local references used for the first experiment were the UTC(NMIJ) and UTC(NICT) signals, which were occasionally steered by step-wise frequency corrections. MARBLE2 provided the master clock and the clock of MARBLE1 was estimated by a piecewise linear function with 60 min interval. For the INAF–NICT sessions, stable free-running hydrogen masers were available at both stations and served as flywheels to maintain the frequencies of the optical frequency standards. Then, a single set of clock offset and rate for the MARBLE2 station was estimated to compare the frequency standards.

• 4. We repeated a process of flagging outliers, re-weighting of observed data, least-square analysis until the reduced $\chi ^2$ converged to unity.

Unlike in conventional geodetic VLBI analysis, ambiguity editing was not necessary because the broadband group delay was unambiguous (see Sect. 2.7).

We confirmed that the a priori data for atmospheric delay provided by the VMF3 appropriately predicted the atmospheric delay by comparison with a GPS PPP solution (see Sect. 3.2). The small magnitude of the atmospheric parameters obtained from the Solve analysis provided additional evidence for the accuracy of the VMF3 a priori delay (see Sect. 4.1).

3 Node-hub VLBI experiments

We conducted two stages of VLBI experiments. Performance was initially tested over short baselines in Japan. In the next stage, intercontinental experiments were conducted for the frequency comparison between two optical lattice clocks reported in Pizzocaro et al. (2021).

3.1 Short baseline experiments (Kashima, Koganei, Tsukuba)

3.1.1 Observations

NMIJ and NICT both generate Coordinated Universal Time (UTC) signals as UTC(NMIJ) and UTC(NICT). Comparisons between them are regularly made available by the Bureau International des Poids et Mesures (BIPM), based on time transfer by GNSS. This makes the NICT–NMIJ baseline a good testbed for frequency comparison. The transportable broadband antenna MBL1 and accompanying broadband DAS were installed at NMIJ in Tsukuba. The MBL2 antenna was installed at NICT Headquarters in Koganei. Both antennas were placed on the roof of the buildings where the respective UTC signals are generated, and incoming signals were transferred from the receiver box of the antenna to the DAS installed in air-conditioned rooms inside the building via a single mode optical fiber cable. The MBL2 antenna was rigidly mounted, while no such fixed attachment was possible for MBL1. Pictures of the two radio telescopes are displayed in Fig. 10.

Fig. 10

(left) MBL1 at NMIJ (Tsukuba) free-standing on protective rubber pads. (right) MBL2 at NICT (Koganei). The temperature controlled ‘receiver box’ including filter, amplifier, and E/O device is visible behind the main reflector (see main text)

VLBI experiments were conducted in this configuration in the period of 2017–2018 as listed in Table 3. The center frequencies of the 1024 MHz wide acquisition channels were 4.1, 6.0, 10.4, and 13.3 GHz. The duration of individual scans was in the range of 10–40 s, and 33.6 scans/h were recorded on average. Sessions were separated by gaps of several days. Data observed at NMIJ were manually transferred to Kashima for correlation processing as physical storage. The data recorded at NICT were directly accessed from Kashima at the time of correlation processing via Network File System (NFS) protocol over the 10 Gbps network provided by JGN (High Speed R&D Network Testbed).

Table 3 VLBI sessions using MARBLE1 and MARBLE2 located at Tsukuba and Koganei

To better achieve the goal of characterizing clock behavior, the session length was not limited to 24 h as usual in geodesy but extended to as long as several days. The measurement instability is characterized based on the post-fit residuals for the individual scans, calculated as the deviation of the observed group delay from the model, as shown for an exemplary dataset in Fig. 11. The scatter of post-fit residuals results not only from the limited delay precision due to SNR, but also from the model fitting of propagation delays, as well as radio source structure effects of the observed sources. The residual plot visually shows no significant time dependence and a mostly statistical distribution characterized by a weighted root-mean-square (WRMS) value of 9.7 ps, roughly equivalent in function to the standard deviation for equally weighted data. The WRMS values over all scans are calculated for each session and included in the table.

Fig. 11

Post-fit delay residuals of the session on 19 Jan. 2018 (MJD 58137)

The initial VLBI sessions to evaluate geodetic performance started in April 2017. Antenna troubles in the summer of 2017 then interrupted the experiments until a series of VLBI sessions for clock comparison started from December 2017. Whereas the sampler clock had previously been reinitialized each session, a continuous run referenced to the local hydrogen maser was maintained after initialization on 18 December 2017, to allow long-term measurements of broadband group delay. After correlation processing by the GICO3 software correlator, post-processing and data analysis were performed as described from Sects. 2.5 to 2.9.

3.1.2 Results

Geodetic results Figure 12 compares obtained baseline lengths to data from two stations of the Japanese GPS observation network GEONET (Hatanaka et al. 2001), located in Tsukuba (0627) and Koganei (3019). The GPS baseline is not identical to the VLBI baseline, but sufficiently close to represent regional crustal deformation. The GPS-based baseline solution (Hatanaka 2003; Nakagawa et al. 2009) data were downloaded from the GSI database² and are displayed in the plot. While the baseline length according to GPS data appear to shorten at a rate of − 5.3 mm/year, the VLBI data show a rate of − 9.5 mm/year. If we consider that the rates observed by VLBI and GPS are consistent over the period from April to June 2017, then the discrepancy in the latter half corresponds to a shift of MBL1 by 4.7 mm east and 3.1 mm north. Such a hypothetical shift (indicated in the plot) might be explained by the installation of the antenna. To protect the roof of the building, this was placed on rubber pads, without rigidly mounting it to the concrete floor (Fig. 3). While the discrepancies do not allow for a conclusive evaluation of baseline length repeatability, we can investigate the stability over periods of several weeks by subtracting an overall linear trend from the data. To account for the two degrees of freedom of the linear fit to the total of eight data points, we estimate the standard deviation representing baseline repeatability by correcting the RMS deviation with a factor of $\sqrt{8/6}$. We then find a repeatability of 1.4 mm for the full span, and 0.7 mm and 1.7 mm for the first and second half, respectively. A 1 mm baseline repeatability approximately coincides with the formal error over the 8 VLBI sessions, corresponding to the 1-sigma error of the least-squares analysis by Solve.

Fig. 12

Time series of the baseline length measured by VLBI sessions in Table 3. Filled circles ‘ ’ indicate NHS VLBI results for the MARBLE1 - MARBLE2 baseline with error bars showing 1-sigma uncertainty of the solution. GPS-based baseline lengths are shown by ‘+’ marks, vertically shifted to match the VLBI results. Open circles ‘ ’ represent a hypothetical shifted antenna position as given in the text

Comparison between UTC (NMIJ) and UTC (NICT) The sampler clocks of the VLBI node stations were referenced to the local time scales UTC(NMIJ) and UTC(NICT). As a result, the VLBI clock difference obtained in the Calc/Solve analysis directly represents the difference between the time scales except for a fixed offset. Hereafter, clock difference data from VLBI measurements are presented as the sum of the estimated clock model and delay residuals of each scan. The precise, unambiguous measurement of the absolute group delay by broadband VLBI then allows for monitoring the long-term behavior of the time scales, even if recording capacity limits individual VLBI sessions to a few days’ duration.

Figure 13 presents the observed time scale difference between UTC(NICT) and UTC(NMIJ) for the period between 18 and 29 December. In these measurements, the result of an initial RSPcal was applied to all measurements. In addition to the VLBI measurements, results from two types of GPS data analysis, using Precise Point Positioning (PPP) and an improved integer ambiguity solution of PPP (IPPP) (Petit et al. 2015), are displayed with shifts for visibility. The clock differences obtained by VLBI and the two GPS-based solutions show identical long-term behavior as displayed in the upper panel, demonstrating consistent determination of absolute group delay across multiple VLBI measurements. The lower panel of Fig. 13 presents the relative differences between solutions. This double-difference must be constant if the two paths of comparison differ only by a constant timing offset. We evaluate the residual slope in terms of the relative clock rate. The rate difference of 2.5 $\times$ $10^{-17}$ s/s observed in VLBI-IPPP suggests that both IPPP and VLBI provide accuracy at this level, whereas PPP exhibits an error on the order of $10^{-16}$ s/s. We attribute the larger scatter (characterized by an Allan deviation of $10^{-13}$ at 300 s interval for VLBI-IPPP) in the VLBI results to the coupling of atmospheric and clock models. This represents the difficulty of estimating atmospheric parameters on a short baseline, where all the stations observe radio sources with almost the same elevation angles. The results nevertheless show that the NHS VLBI measurements enable accurate frequency evaluation and long-term tracking of clock behavior and provide a valuable method for time scale comparison.

Fig. 13

Time scale comparison of UTC(NICT)–UTC(NMIJ) with broadband VLBI, GPS(PPP), and GPS(IPPP). The lower panel shows the pair-wise double difference. Since a constant time offset is of no concern, data in both plots are offset for visibility

3.2 Intercontinental baseline experiments (Kashima, Medicina, Koganei)

3.2.1 Observations

The broadband VLBI antenna MBL1 was moved from NMIJ at Tsukuba to INAF at Medicina in August 2018 (Fig. 14). The transportable design of the antenna, careful preparation of electricity, optical fiber, and antenna foundation allowed for assembly and start-up of the station to be completed within a week. A fringe test experiment with the Ishioka 13 m antenna over an 8000 km baseline was successfully performed only three days later. The MBL1 antenna was installed about 500 m away from the Medicina 32 m VLBI station to avoid blocking sky coverage. The observed RF signal was converted to optical within the temperature controlled receiver box and then transferred to the Medicina VLBI station over approximately 600 m of optical fiber placed in an underground trench. The O/E converter and the K6/GALAS sampler for RFDS data acquisition were installed in the observation room of the 32 m telescope. The purpose of the installation was to compare optical atomic frequency standards between NICT and INRiM (Istituto Nazionale di Ricerca Metrologica, Torino), connected to INAF via optical fiber link (Calonico et al. 2014). Several hydrogen maser frequency references were used as flywheels to bridge interruptions in operation of clocks and links. More details on the frequency link are described in a separate publication (Pizzocaro et al. 2021).

Fig. 14

Broadband 2.4 m antenna MBL1 installed at the campus of Medicina Radio Astronomical Station of INAF in 2018

VLBI sessions were conducted from October 2018 to July 2019 and are listed in Table 4. The observation schedule files were generated using the Sked program (Gipson 2010), which determines the necessary scan length for a specified SNR by evaluating source flux and SEFD of the observation stations. We used slightly higher (i.e., worse) SEFD values to describe station performance conservatively. As a result, the observed SNR exceeded what was necessary for fringe detection with the 2.4 m antennas in the majority of observations. The effective bandwidth of our frequency array (now operated with bands centered at 6.0, 8.5, 10.4, and 13.3 GHz with 1024 MHz width each) was 2.685 GHz, which corresponds to a theoretical limit of 6 ps of delay uncertainty [Eq. (6)] per scan even for an SNR as low as 10. A scan with SNR = 15 can be divided into 2 scans with SNR = 10, and we applied such subdivision targeting a final SNR of around 10. The ‘Subdivided #Scans’ in the table represent this. This procedure increases the number of observations that contribute to statistics and provides a finer temporal resolution of probing the atmospheric conditions, which is also an important design goal of VGOS. After subdivision, the typical scan duration of 30 s or 60 s corresponds to about 50 subdivided #Scans per hour.

A series of sessions during October 2018–February 2019 focused on comparing optical atomic frequency standards to a fractional uncertainty of only a few parts in $10^{16}$, to which the VLBI link contributed less than $10^{-16}$ (Pizzocaro et al. 2021). Since the duration of the individual VLBI sessions was limited by data recording disk capacity, experiments from 30 May to 18 July 2019 explored the feasibility of intermittent operation. The session on 30 May tested a duty cycle of 6 h operation followed by a 12 h break, the sessions of 12 June, 18 and 31 July implemented 8 h operation intervals separated by 12 h breaks. However, we found the stochastic behavior of the hydrogen maser references over long sessions to be inadequately described by the current first order clock rate model. This affects the results through parameter coupling, and we have not yet been able to demonstrate a quantifiable advantage of this approach. An approximately twice larger WRMS on 30 May 2019 indicates that modeling of delay behavior for the intermittent observation in this period was not successful.

Table 4 VLBI sessions between MARBLE1M and MARBLE2 with the NHS VLBI scheme

3.2.2 Results

Observed zenith delay Figure 15 displays the atmospheric total zenith delay from VLBI observation, a priori delay by VMF3, and GPS IPPP and GPS PPP solutions for the stations at Medicina and Koganei. All results show similar values and overall trend, confirming the prediction of the VMF3 dataset. VLBI results are shown as the sum of a priori (VMF3 zenith) model input and the estimated residual ($\delta \tau _{\mathrm{wh}}^{\mathrm{z}}$). Part of the discrepancy between the zenith delay estimated by GPS and VLBI simply represents the anisotropic atmosphere: The azimuthal directions for MARBLE1M and MARBLE2 stations were in the range of $-\,40^{\circ }$ to $+\,50^{\circ }$, and $-\,50^{\circ }$ to $+\,130^{\circ }$, respectively. The mean zenith delay parameter displayed here then absorbs a contribution from the anisotropic components. The estimated instability from atmospheric delay prediction is further discussed in Sect. 4.1.

Fig. 15

Comparison of atmospheric zenith delays for Medicina (top) and Koganei (bottom). Results obtained by VLBI observations are indicated by filled circles ‘ ’ and filled triangles . Predictions by the VMF3 dataset are indicated by open circles ‘ ’ and open triangles in the same order. The values obtained by GPS PPP and IPPP solutions are superimposed with solid and broken lines, respectively

Geodetic results Baseline lengths determined by the sessions in Table 4 are displayed in the upper panel of Fig. 16. Both standard VLBI solutions for the MARBLE1M-KASHIM34 baseline and NHS VLBI solutions for the 2.4 m antenna pair are plotted. The baseline lengths between the station MEDI (at Medicina) and KGNI (at Koganei) of the IGS (International GNSS Service) based on PPP solutions are also included. The plot shows a consistent trend in baseline length for all cases. The WRMS values for the deviation from a linear trend are 15.0 mm, 17.5 mm for MARBLE1M-MARBLE2, MARBLE1M-KASHIM34 by VLBI, and 13.4 mm for MEDI-KGNI by GPS, respectively. Although it may appear surprising to find a better stability for the synthesized MARBLE1M-MARBLE2 baseline, it is consistent with the expected partial cancellation of systematic effects resulting from the large size of the Kashima 34 m antenna. To provide context for the baseline length repeatability (BLR), we re-evaluated IVS data in the R1 and R4 series, obtained in 2011–2013. We chose Wettzell as reference station for clock and station coordinates. One hundred sessions including 201 baselines and 25 stations were analyzed with the Calc/Solve software in the same configuration as used for MARBLE1M-KASHIM34 data, except that atmospheric gradient estimation was re-enabled. Baselines were excluded when less than 10 data points were available, and all the data reported by Tsukuba station before MJD 55,800 were removed due to nonlinear post-seismic drift following the 2011 Tohoku earthquake. In addition, baselines of obvious outliers related to the TIGO station affected by the Chile earthquake in 2010 were excluded. Then, BLR as WRMS deviation from a linear trend was determined for 78 baselines and is displayed in the lower panel of Fig. 16. The BLR obtained by NHS VLBI between our 2.4 m antennas compares well to the results of the standard VLBI sessions. In relation to the baseline length, we observe a fractional repeatability of 1.7 parts-per-billion. This is comparable to the IVS-R1 and IVS-R4 sessions using S/X antennas in the 11–32 m range. The IVS-R1 and R4 sessions take advantage of network observations with five to eleven stations. This provides better sky coverage and a better estimation of atmospheric contributions. Multi-baseline analysis also provides additional constraints on the station coordinates. To estimate the effects on BLR, we reevaluated the baselines Wettzell–Tsukuba (Wz–Ts) and Wettzell–Kokee (Wz–Kk) under different analysis conditions. We compared the network solution to a single baseline solution, each calculated with or without atmospheric gradient estimation. The results of this case-study are summarized in Table 5. We find that the repeatability of single baseline measurements is severely degraded for Wz–Kk, likely because the nearly antipodal station locations for this 10,360 km long baseline limit the number of mutually visible sources. For Wz–Ts, where the baseline length of 8440 km is comparable to Medicina–Koganei (8780 km), a similar effect is visible, but less pronounced. Here, the network measurements show a repeatability of 1.2 ppb, a 20% reduction in scatter compared to the repeatability of 1.5 ppb observed over the single baseline. The analysis suggests that for this baseline length, a large part of the improvement in network measurements results from the improved directional sampling of the atmospheric conditions. If we disable the estimation of the atmospheric gradient coefficients that represent this in the internal model, network measurements no longer result in a dramatically improved repeatability, while there is no significant impact on the repeatability of the single baseline result. We expect that NHS VLBI will similarly benefit from network observation in terms of atmospheric modeling and baseline repeatability as seen for the Wz–Ts measurements.

Fig. 16

Upper panel shows baseline lengths (as deviation from the value specified in the figure) measured by VLBI and GPS PPP solutions of similarly located MEDI and KGNI stations. Lower panel shows BLR for MBL1–MBL2 (‘ ’) and our re-evaluation of IVS-R1 and R4 sessions during 2011–2013 (‘ ’). The solid line represents a regression function of BLR for IVS-R1 and R4 sessions in a preceding study (Teke et al. 2008)

Table 5 Baseline length repeatability (BLR) calculated as weighted RMS for Wz–Ts and Wz–Kk baselines (Wz: Wettzell, Ts: Tsukuba, and Kk: Kokee)

4 Discussion

4.1 Uncertainty of the group delay observable

The group delay observable is affected by various factors, including instrumental delay and source structure effects. In this section, we discuss such contributions in terms of how much they are expected to affect the result of a single scan over the course of an observation session. We will refer to this as the group delay instability in the following, conceptually representing a RMS-deviation from the mean. For application to geodesy and clock frequency comparisons, constant delay-offsets as might be introduced by, e.g., signal distribution or electronic propagation delays do not affect the results. The parameters in the following sections represent the clock comparisons between Koganei and Torino (Pizzocaro et al. 2021) over the period between 14 October 2018 and 14 February 2019. Table 6 gives an overview of the contributions discussed here, with an overall group delay instability of 27.5 ps under favorable conditions.

Table 6 Group delay instability for measurement by NHS VLBI

Sensitivity-dependent instability The broadband group delay precision is related to the bandwidth and SNR as described by Eq. (6). The EBW of our observation frequency array is 2.685 GHz. An observation with SNR = 10 thus yields 5.9 ps of delay precision for a single scan. The average of the delay precision over the sessions is 6.4 ps, which represents the per-scan instability.

Instrumental Delay When we consider the entire signal path from the dish to digital data, temperature and antenna motion may introduce time-varying delay contributions in multiple ways. The GALA-V system uses signal transmission by optical fiber to obtain lower sensitivity to ambient temperature change and mechanical stress. Berenz et al. (2009) have measured the response of signal transmission to temperature variation and mechanical stress from moving radio telescopes of 100 m diameter at Bonn and 32 m diameter at Medicina. Their results (specifically Fig. 7) show the response to the mechanical stress change by motion in azimuth and elevation to be much smaller than the response to ambient temperature, and we adopt a contribution of 0.5 ps to the group delay instability for this. For the larger thermal contribution, we compare a hypothetical group delay measurement at an arbitrary time within the session to an identical measurement at the start. When considering the 600 m long connecting fiber at Medicina, a temperature change of 5 K results in a 7.6 ps change in delay for silica glass with a typical thermal expansion coefficient of $5\times 10^{-7}$/K and index of refraction $n = 1.47$. The temperature uncertainty of 5 K represents the additional insulation provided by placing the fiber in an underground trench. For the Koganei station, we assume a larger temperature uncertainty of 15 K, but the shorter fiber length of 50 m nevertheless reduces the contributed group delay instability to 1.9 ps.

Similarly, we adopt a group delay instability of 10 ps for the combined contribution of the acquisition systems. As these are all operated in air-conditioned environments, we expect thermal effect here to be uncorrelated with the effects on the fibers. At short time scales, we have characterized the sampling jitter of each K6/GALAS system as 0.2 ps by measuring the frequency-dependent carrier phase fluctuation (Takefuji and Ujihara 2013).

Ionospheric Delay The determination of the observed group delay by WBWS corrects for the ionospheric delay contribution based on estimates of the relative ${\mathrm {TEC}}$ from the cross-correlation data. Kondo and Takefuji (2016) evaluated the error of ${\mathrm {TEC}}$ estimation as about 0.1 TECU (1 TECU = $10^{16}$ electrons per m$^2$) for large SNR ($\ge 20$) and a baseline within 100 km. In the case of the intercontinental baseline of KASHIM34-MARBLE1M, the SNR was limited due to reduced correlated flux. Later work (Kondo and Takefuji 2019) extended the algorithms used in WBWS to allow determination of ${\mathrm {TEC}}$ also at low SNR. We expect similar performance as demonstrated for the earlier algorithms, although at this time the only experimental confirmation available is a comparison to GNSS-based global ionosphere maps, which restrict the accuracy to approximately 2 TECU due to their uncertainty and limited spatial and temporal resolution. The uncertainty of TEC measurement by broadband VLBI should be well below the uncertainty of 1 TECU or less in conventional S/X-band VLBI (Sekido et al. 2003), and we take this value as a pessimistic estimate. For a known error in ${\mathrm {TEC}}$, the effect on the observed group delay can be found by considering the correlation of the parameters in the WBWS. We re-evaluated the data of 15 and 25 January 2019 with externally applied bias up to ± 3 TECU to find a linear dependence of 17.2 ps/TECU (Fig. 17). Differentiation of Eq. (3) with regard to frequency yields

$$\begin{aligned} \delta \tau = A \frac{ {\mathrm {TEC}} }{2 \pi f^2} = 1.34\times 10^{-9} \mathrm{s } \cdot \left( \frac{ {\mathrm {TEC}} }{ {\mathrm {TECU}}} \right) \left( \frac{f}{{\mathrm {GHz}}}\right) ^{-2} \end{aligned}$$

(18)

indicating an expected dependence on the selected frequency array. We found 24.5 ps/TECU for a frequency array of (4.912, 6.000, 8.412, 10.712 GHz), and 20.4 ps/TECU for the array of (5.300, 6.242, 8.412, 10.712 GHz) used in our other experiments. Since the simplified calculation of Eq. (18) does not precisely reproduce the observed dependence, we take the experimentally determined coefficient of 17.2 ps/TECU to estimate the effect of the ionospheric delay determination. For the expected ${\mathrm {TEC}}$ uncertainty of 0.1 TECU, this yields 1.7 ps group delay instability across the scans obtained in a measurement session, while the pessimistic value of 1 TECU corresponds to 17.2 ps.

Fig. 17

Correlation between delay and ionospheric TEC. Using WBWS, changes of group delays were computed by shifting the TEC value up to ± 3 TECU from the estimated nominal value

Fig. 18

Means of residual atmospheric zenith delay $\delta \tau _{\mathrm{wh}}^{\mathrm{z}}$ for MARBLE1M ( ), MARBLE2 ( ) during each session, and average over all results ($\blacksquare$). The error bars are root-mean-square deviation from the mean

Atmospheric Delay The VMF3 dataset (re3data.org 2019) provides wet zenith delay and atmospheric gradient parameters every 6 h. Derived from the numerical weather model of the ECMWF, these parameters provide one of the best approximations of real propagation delay. When applied as part of the Solve analysis, we expect the least-squares fit to find only minimal corrections in terms of the zenith atmospheric parameters. Here, we specifically look at the estimated residual atmospheric zenith delay parameter ($\delta \tau _{\mathrm{wh}}^{\mathrm{z}}$) (see Fig. 18), for which we find a weighted mean of only − 3.3 ps over all sessions. This demonstrates the accuracy of the VMF3 dataset on longer time scales, although we find standard deviations of up to 50 ps over the individual residuals obtained in separate sessions. While we consider this a useful measure of the prediction error of VMF3, composed of estimation error and deviation from the interpolation model, the large deviation does not manifest in the observed instabilities shown as WRMS in Table 4. This suggests that the least-squares estimate of the residual zenith delay ($\delta \tau _{\mathrm{wh}}^{\mathrm{z}}$), which is updated at an hourly rate, provides additional suppression of the resulting delay instability. We look at the variations between the subdivided scans (Sect. 3.2.1) to investigate atmospheric instabilities over shorter times. For these subdivisions of continuous scans of the same radio source, we expect delays from ionospheric effects, instrument or source structure to be effectively constant, leaving only the contributions from limited SNR and atmospheric changes. We suppose the variance between successive delay residuals, separated by 90 s or less, to be composed of a sensitivity-limited $\sigma _{\mathrm{sens}}$ and an atmospheric component represented as follows:

$$\begin{aligned} \sigma _\tau ^2 (\alpha )= & {} \langle (\delta \tau _{{\mathrm {i}}} - \delta \tau _{\mathrm {i+1}})^2 \rangle /2 \end{aligned}$$

(19)

$$\begin{aligned}= & {} \sigma _{\mathrm{sens}}^2 + (c_{{\mathrm {atm}}} \cdot \alpha )^2, \text {with} \end{aligned}$$

(20)

$$\begin{aligned} \alpha= & {} \sin ^{-1}(\epsilon _1) + \sin ^{-1}(\epsilon _2), \end{aligned}$$

(21)

where $\langle \cdot \rangle$ denotes average and $\delta \tau _{\mathrm{i}}$ are delay residuals of the subdivided scan. The coefficient $c_{{\mathrm {atm}}}$ quantifies an atmospheric instability proportional to the elevation parameter $\alpha$ defined by Eq. (21), approximating the contributions of mapping functions for two stations with different elevation angles $\epsilon _1$ and $\epsilon _2$. The individual differences are binned according to $\alpha$, and a RMS average is calculated. We take this as an estimate of the 1-sigma uncertainty characterizing the instability of a single residual. The data are displayed in Fig. 19, fitted with the model of Eq. (20). We obtained $\sigma _{\mathrm{sens}}$ = 6.5 ps and $c_{{\mathrm {atm}}} =$ 1.76 ps. Combined with a quadratic mean of elevation parameter $\alpha$ = 4.64 over all sessions, this yields 4.64 $\cdot$ 1.76 ps = 8.2 ps for the expected per-scan instability of atmospheric delay. The fitted $\sigma _{\mathrm{sens}}$ carries a significant uncertainty due to the extrapolation to $\alpha = 0$ (no atmosphere) and is consistent with the expected delay precision of 6.4 ps.

Fig. 19

Short-term atmospheric contribution to the delay residual. Symbols show the per-scan instability $\sigma _\tau$ as defined in Eq. (19) binned according to the intervals marked on the horizontal axis. The analysis represents data obtained from October 2018 to February 2019. Outlier detection eliminated less than 1% of the overall data. The source 2022+542, with an atypically large RMS deviation, was found to dominate the results at low $\alpha$ and excluded from the evaluation. The error bars represent the statistical 1-sigma confidence interval of the binned data, determined from the $\chi ^2$ statistical evaluation. The blue solid line represents Eq. (20) fitted to the data expressed in terms of $\sigma _\tau ^2$ with weights according to the statistical uncertainties. The data are expected to deviate from the model by more than the statistical expectation due to the different composition of individual bins in terms of contributing sources and daily conditions

Radio Source Structure Radio source structure effects specific to our broadband observations are discussed in Sects. 4.2 and 4.3. The observed sources were selected from ICRF3 (the 3rd realization of the International Celestial Reference Frame) (Charlot et al. 2020). The source structure effect in S/X band VLBI is evaluated by the continuous structure index (SI) (Ma et al. 2009), which describes the median of group delay magnitude from ground-based VLBI observations as

$$\begin{aligned} \mathrm{SI} = 1 +2 \cdot \log (\tau _{\mathrm{med}} / {\mathrm {ps}}). \end{aligned}$$

(22)

SI values of the sources are provided by IERS Technical Note No. 35 and the Bordeaux VLBI Image Database³ and range from 1.83 to 4.12 for the 35 sources in the sessions, with an average of 2.87. The corresponding $\tau _{\mathrm{med}}$ are in the range 0.3–36.3 ps. Since this delay deviation varies symmetrically around zero with baseline orientation (Charlot 1990), we approximate the source structure delay observed at a random time as a Gaussian distribution with zero mean and a standard deviation according to a quadratic mean ${\bar{\tau }}_{\mathrm{med}}$. Weighting according to the number of scans for each source yields ${\bar{\tau }}_{\mathrm{med}}$ = 16.7 ps. While this accounts for the source structure contribution to the group delay instability as observed in the S/X band, broadband VLBI has an additional sensitivity to frequency-dependent source position variation. From surveying imaging data collected in the Bordeaux VLBI Image Database, we estimate that such effects may appear as a 0.1 mas (and no larger than 0.2 mas) source position variation over the 8780 km baseline and over the frequency range of 6–13.3 GHz in our observations. This corresponds to an instability of 14.2 ps, with a 28.5 ps worst case estimate. Adding both contributions in quadrature, we find an expected group delay instability contribution due to source structure of 22.0 ps (33.0 ps worst-case). Although Sect. 4.2 presents an investigation into source structure directly based on our broadband measurements, this presently provides insufficient insight into the resulting group delays: Even for a structured source, symmetry can cause the group delay error to vanish, since the imaginary term [Eq. (11) of Charlot (1990)] diminishes. Total instability of broadband group delay All contributions discussed above are summarized in Table 6, which separately lists expected and pessimistic values. The total expected group delay instability, representing single scans taken during the same session, is then 27.5 ps, with a pessimistic value of 40.7 ps. The obtained values are in reasonable agreement with the WRMS of the delay residuals during the observation sessions reported in Table 4, which vary from 22 to 39 ps. We can also relate them to an approximate uncertainty in a clock comparison as 27.5 ps/30 h = $2.5 \times 10^{-16}$ for a 30 h session (with a pessimistic value of 3.8$\times 10^{ -16}$). The baseline repeatability discussed in Sect. 3.2 depends on averages over full sessions and the variation of systematic delays over longer intervals of time. Additional difficulties arise from relating delay and baseline uncertainties, a process that depends on many parameters of observation session and analysis, not limited to scan interval, source selection, and weighting of the data in analysis. Many contributions are similar in magnitude, however, and it is nonetheless of interest to compare c $\cdot$ 27.5 ps = 8.2 mm (12.2 mm for the pessimistic case) to the observed baseline repeatability of 15.0 mm presented in Sect. 3.2.

4.2 Correlation amplitude in uv-track

The influence of radio source structure on group delays has been investigated through closure delays by Xu et al. (2016, 2019). This provides delay components related to the radio source while eliminating the contribution from the local observation sites. In our measurement, we did not achieve sufficient SNR on the baseline between the 2.4 m diameter antennas to obtain closure delay data. We find that the correlation data output of the KASHIM34 (VH)–MARBLE1M (V) baseline after the polarization synthesis (see Sect. 2.5) nevertheless provides interesting insight into radio source structure. The thirteen VLBI sessions in the period from October 2018 to February 2019 (Table 4) were performed with the same array of frequency bands (6.0, 8.5, 10.4, and 13.3 GHz). Some high latitude sources have been observed for longer than half a day. The left side of Fig. 20 displays the color contours of correlation amplitudes mapped on the uv-track for each of the four frequency bands. The right-side column shows correlation amplitude variation for each source as a function of position angle (counted clockwise from North) of the projected baseline. This uv-track-dependent amplitude variation is attributed to the source structure.

Fig. 20

Correlation amplitude contour on uv-track (left) and correlation amplitude versus position angle of the projected baseline (right) of KASHIM34-MARBLE1M. The plots show merged data of 13 sessions between 14 October 2018 and 14 February 2019 for five selected sources. Individual tracks on the left correspond to the frequency bands, with lower frequency close to the center. In the right panels, red downwards triangles (‘ ’), brown squares ( ), blue circles (‘ ’), and black upwards triangles ( ) represent observation frequencies of 6.0 GHz, 8.5 GHz, 10.4 GHz, and 13.3 GHz, respectively

Source 0059+581 is classified as a defining source of ICRF3 with small source structure effects (e.g., Fey and Charlot 1997; Anderson and Xu 2018). Although it is reported by Fey and Charlot (1997) to be composed of multiple weak sub-components within 1 mas distance from the main component, resulting in structure indices of 0.81 for S-band and 2.42 for X-band, these components are not resolved in our uv-track data. We observe an amplitude modulation as expected for a simple elliptical shape. Our data show similar patterns of amplitude variations for all the frequency bands, suggesting similar source shapes in the four observing frequencies.

Source 1044+719 is another source with small structure effects. Fey and Charlot (2000) evaluated it to have structure indices of 0.59 and 0.7 for S and X-band. Anderson and Xu (2018) assigned it to a group of sources with minimal source structure effect based on X-band data of the recent CONT14 experiment. Our data show almost constant amplitude in each frequency band, suggesting small structure effects in not only single band, but also broadband group delay.

Sources 0212+735 and 3C418 are classified as exhibiting large structure effects by Anderson and Xu (2018). This is clearly confirmed in our observation data. A particular point to note here is that for source 0212+735, the data for the 6 GHz band show an amplitude change with position that differs significantly from the other bands. The coupling of frequency-dependent source structure to the estimated ionospheric correction may then result in additional errors as discussed below.

Table 7 Model parameters of the source 1928+738 for each band obtained by fitting to correlation amplitude data

Fig. 21

Two-dimensional ellipsoid Gaussian model (lines) fitted to correlation amplitude data (symbols) of the source 1928+738 for each band. The panel on the right-hand side shows the modeled ellipsoids, corresponding to the parameters listed in Table 7. Contour lines mark 50% of peak flux for each band, as indicated by the line pattern and color

4.3 An interpretation of correlation amplitude for 1928+738

The large amplitude variation for source 1928+738 in Fig. 20 implies that the source is compact in one direction and larger than half of the fringe interval in the other. We apply a Gaussian source model as introduced by Charlot (1990) to describe this visibility dependence on the projected observing baseline.

Hereafter, our discussion will focus on the amplitude. The interferometric visibility differs from the correlation coefficient by calibration with the system noise temperature and the antenna gain curve. For simplicity, no correction is applied for elevation-dependent gain change, which is on the order of few percent for the Kashima 34 m antenna and expected to be negligible for the small MBL1 antenna. We similarly neglect changes of the system noise temperature during source tracking. Hereafter, we regard the correlation amplitude as proportional to the visibility amplitude to within 10 percent of error.

To model the source distribution, the expected visibility amplitude for an ellipsoidal 2D-Gaussian source model is evaluated on the uv-track for each frequency band. The ellipsoid axis ratio determines the ratio between peak and bottom of the variation. The inclination of the ellipsoid corresponds to the phase of the amplitude changes with respect to the position angle. Source model parameters obtained by manual optimization are summarized in Table 7, and correlation amplitude data and model are displayed in Fig. 21. A contour map of the derived source models is also shown. Since the models are derived independently for each frequency, their relative positions are not obtained in this analysis. The orientation and size of the models are roughly consistent with the results by Fey and Charlot (2000) for X-band data, with a deviation observed mainly for the axis ratio. Fey and Charlot (2000) remark that source models could not be uniquely determined due to incomplete sampling in the uv-plane in their observation with the VLBA (Very Long Baseline Array). Our modeling as a single component may also contribute to the difference, but time evolution of the structure is another intriguing possibility.

Despite the limited sensitivity of the 2.4 m antenna, the spatial resolution is sufficient to resolve sub-milliarcsecond structure. The data in Fig. 21 show a larger deviation from the model at higher observation frequencies, suggesting that these frequencies better resolve the structural complexity of source 1928+738.

It is notable that while the symmetric components of our model do not cause any delay contribution (e.g., Charlot 1990), more complex shapes generally will. The observation of band-dependent source structure has additional implications for broadband VLBI.

Any displacement of the brightness distribution between the bands will affect the estimate of the ionospheric contribution from assuming $\phi = A \cdot {\mathrm {TEC}}/f$ for the correlation coefficient in the WBWS process. A perturbation added to the phase response by the inter-band source structure will then affect the extracted group delay observable as well. Thus, the source structure effect depends not only on the observing frequency and baseline vector as a function of time, but also depends on the technique of ionospheric delay estimation in data analysis. Future surveys will be needed to identify sources suitable for geodetic and metrological broadband VLBI.

5 Conclusion

A broadband VLBI system has been developed inspired by the VGOS concept. The K6/GALAS system provides high-speed data acquisition using the RF-direct sampling technique, leading to increased phase stability over time and across frequency bands. This enables precise absolute group delay measurements from broadband data. A technique based on cross-spectrum observation data of natural radio sources provides phase calibration in the absence of a suitable Pcal system.

Our novel broadband feed can be fitted to conventional Cassegrain high-gain antennas. In addition to the stationary 34 m antenna at Kashima, we performed experiments with two stations that can easily be moved to desirable locations. By combining measurements over the two baselines involving the highly sensitive antenna, a virtual delay observable between the transportable stations can then be obtained. This new node-hub style VLBI scheme also provides cancelation of common delay errors associated with large dish sizes as used in the hub station, a role that can be fulfilled by any standard VGOS antenna. We hope that the use of smaller node antennas in combination with VGOS telescopes will contribute to geodetic VLBI, particularly in the southern hemisphere. Choosing the location of the hub will provide opportunities to manipulate and possibly attenuate source structure effects.

Measurements over both short and intercontinental distances confirmed the performance of the full system. NHS VLBI observations with the 2.4 m antenna pair revealed baseline repeatability at the same level as IVS-R1 and R4 legacy S/X-VLBI observations. The system also enabled long-distance frequency comparisons where the transportable node stations are near the standards to be compared (Pizzocaro et al. 2021). Here, we evaluated the absolute delay uncertainty over a 8780 km baseline between the 2.4 m antenna pair to be less than 50 ps, except for a yet-to-be-calibrated fixed timing offset. This opens new avenues for a future intercontinental network serving time and frequency metrology (Haas 2020).

In the broadband VGOS era, high group delay precision comes with increased sensitivity to frequency-dependent source shape. Even the single baseline between a 2.4 m diameter antenna and the Kashima 34 m antenna revealed information on radio source structure. Further studies and astronomical surveys targeting this subject should be encouraged. The development of new techniques is required for calibration and for assigning an extended structure index for VGOS.

Appendices

Relation between the VLBI clock rate and the frequency difference between frequency standards at VLBI stations

We work in a geocentric reference system (in practice realized by the ITRF), using terrestrial time TT as a coordinate time, hereafter noted t. For a signal from a distant radio source arriving at station 1/2 at coordinate time $t_1/t_2$, the virtual delay observable $\tau ^{{{{\mathrm {obs}}} *}}$ [Eq. (14)] of NHS VLBI contains the clock offset $\tau ^{\mathrm {clk}}$ in the same way as $\tau ^{\mathrm {obs}}$ in regular VLBI (e.g., Cannon 1999; Schuh and Behrend 2012):

$$\begin{aligned} \tau ^{{{{\mathrm {obs}}} *}} = \tau ^{\mathrm {geo}} + \cdots + \tau ^{\mathrm {clk}} \end{aligned}$$

(23)

where the measured delay $\tau ^{\mathrm {obs}} = \tau _2 - \tau _1$ is the difference of the reading of clock 2 at time $t_2$ and the reading of clock 1 at time $t_1$, and where the VLBI model reproduces the coordinate time difference as $\tau ^{\mathrm {geo}} + \cdots = t_2 - t_1$ (Petit and Luzum 2010). Therefore, we have

$$\begin{aligned} \tau ^{\mathrm {clk}} = (\tau _2(t_2) - \tau _1(t_1)) - (t_2 - t_1). \end{aligned}$$

(24)

With $r_i = \frac{{\mathrm{d}} \tau _{\mathrm{i}}({\mathrm{t}})}{{\mathrm{d}} t} - 1$ as the coordinate rate of clock i, we have

$$\begin{aligned} \tau _i(t) = \tau _i(t_0) + (1 + r_i) \cdot (t - t_0) \end{aligned}$$

(25)

where $t_0$ is an arbitrary time reference in an interval over which the rates may be considered constant. Equation (24) then gives

$$\begin{aligned} \tau ^{\mathrm {clk}} = {\varDelta }^0 + (r_2 - r_1) \cdot (t_1 - t_0) + r_2 \cdot (t_2 - t_1) \end{aligned}$$

(26)

where ${\varDelta }^0$ is the desynchronization of the clocks at $t_0$. Differentiating Eq. (26), we get

$$\begin{aligned} \frac{{\mathrm{d}} \tau ^{\mathrm {clk}}}{\mathrm{d} t_1} = (r_2 - r_1) + r_2 \cdot \frac{{\mathrm{d}}(t_2 - t_1)}{{\mathrm{d}} t_1}. \end{aligned}$$

(27)

In a standard VLBI station, the frequency of the local hydrogen maser is adjusted to the rate of UTC (and thus to TT) within about $10^{-12}$ such that $r_i$ and the first term of Eq. (27) are of order $10^{-12}$. For two stations on Earth, the VLBI delay rate is typically of order $10^{-6}$, and thus, the last term of Eq. (27) is of order $10^{-18}$. The VLBI determined clock rate is also equivalent to $\frac{{\mathrm{d}} \tau ^{\mathrm {clk}}}{{\mathrm{d}} t_1}$ to within about $10^{-12}$ in relative value. Thus, in practice the VLBI determined clock rate is equivalent to the difference $(r_2 - r_1)$ within about $10^{-18}$.

If the goal of the VLBI experiment is to compare other frequency standards operating close to the hydrogen masers, their coordinate rate difference may be obtained from $(r_2 - r_1)$ by local frequency comparison. If the standards are primary frequency standards (Cs fountains) or secondary frequency standards (e.g., optical lattice clocks), which aim at generating an ideal proper time at their location, an adequate relativistic transformation then provides the comparison of the proper times generated by the two remote standards.

Phase calibration and wideband bandwidth synthesis

VLBI correlation processing begins with congregating the streams of observation data from the two stations to find the cross-spectrum $S_k$ representing the complex correlation coefficients in the frequency domain. During this process, geometrical delay and the Doppler frequency shift resulting from its change due to Earth’s rotation are compensated by using an accurate geometrical delay model. The cross-spectrum is then interpreted in terms of the phase residuals as

$$\begin{aligned} S_k&= S_{k,0} \exp \left[ -j \left\{ 2 \pi f (\tau _{g,k} - \tau _{g,\mathrm{apr}}) \right. \right. \nonumber \\&\quad \left. \left. -A \frac{{\mathrm {TEC}}}{f} + \phi _{\mathrm{ins},k} (f) + \psi \right\} \right] , \end{aligned}$$

(28)

where $\tau _{g,k}$ is the group delay component in the frequency band k, $\tau _{g,\mathrm{apr}}$ is an a priori delay model, ${\mathrm {TEC}}$ is the differential total electron content representing the ionospheric contribution, and $\psi$ absorbs other phase contributions including radio source structure. The physical constant A is the same as in Eq. (3) and $j= \sqrt{-1}$. For phase calibration with a radio source (RSPcal), we first define the reference cross-spectrum and its phase as

$$\begin{aligned} S_k^{\mathrm {ref}}= & {} S_{k,0}^{\mathrm {ref}} \exp \big ( -j \phi _k^{\mathrm {ref}}\big ) \end{aligned}$$

(29)

$$\begin{aligned} \phi _k^{\mathrm {ref}}= & {} 2 \pi f {\Delta } \tau _{g,k}^{\mathrm {ref}} +\phi _k^{\mathrm {ref}}(f)+\psi ^{\mathrm {ref}}. \end{aligned}$$

(30)

The parameters ${\Delta } \tau _{g,k}^{\mathrm {ref}}$ and $\phi _k^{\mathrm {ref}}(f)$ can be determined from a reference scan of a source with ideally negligible structure effects. They include misalignment of sample timing between signal channels, instrumental phase characteristic, and the ionospheric phase contributions of the reference scan. In similar form as Eq. (28),

$$\begin{aligned} {\Delta } \tau _{g,k}^{\mathrm {ref}}= & {} \tau _{g,k}^{\mathrm {ref}} - \tau ^{\mathrm {ref}}_{g,\mathrm{apr}} \end{aligned}$$

(31)

$$\begin{aligned} \phi _k^{\mathrm {ref}}(f)= & {} -A \frac{{\mathrm {TEC}}^{\mathrm {ref}}}{f} + \phi _{\mathrm{ins},k}(f). \end{aligned}$$

(32)

Phase calibration by RSPcal assumes that instrumental phase $\phi _{\mathrm{ins},k}(f)$ is invariant with time and differential phase contribution of radio source structure is negligible. Applying RSPcal to the cross-spectrum of a scan i yields

$$\begin{aligned} S'^i_k= & {} S_k^i \cdot \exp (j \phi _k^{\mathrm {ref}}) =S_{k,0} \cdot \exp \big (-j {\Delta } \phi _k^i \big ) \end{aligned}$$

(33)

$$\begin{aligned} {\Delta } \phi _k^{i}= & {} 2\pi ({\Delta } \tau _{g,k}^i - {\Delta } \tau _{g,k}^{\mathrm {ref}} ) -A \frac{{\mathrm {TEC}}^i - {\mathrm {TEC}}^{\mathrm {ref}}}{f} + {\Delta } \psi , \nonumber \\ \end{aligned}$$

(34)

$$\begin{aligned} {\Delta } \tau ^i_{ g,k}= & {} \tau ^i_{ g, k} - \tau ^i_{ g,\mathrm{apr}}. \end{aligned}$$

(35)

Equation (33) represents the calibration process and defines the phase content ${\Delta } \phi ^i_k$. Equation (34) is a model for its interpretation. Instrumental phase is eliminated by the calibration, but the differential quantities of group delay and ionospheric ${\mathrm {TEC}}$ are retained, although the latter appears as ${\Delta } {\mathrm {TEC}} = {\mathrm {TEC}}^i - {\mathrm {TEC}}^{\mathrm {ref}}$. ${\Delta } \psi$ is a fixed phase constant.

The algorithmic details of the wideband bandwidth synthesis have been presented in (Kondo and Takefuji 2016, 2019) and are beyond the scope of this article. In summary, the WBWS is a process to find a set of parameters (${\Delta } \tau , {\Delta } \dot{\tau }, {\Delta } {\mathrm {TEC}}$) that maximize the magnitude of coherent integration F of the cross-spectrum.

$$\begin{aligned}&F({\Delta } \tau , {\Delta } \dot{\tau }, {\Delta } {\mathrm {TEC}}) = \left| \sum _{k=1}^n \int _{0}^{T} \int _{-B/2}^{B/2} \right. \nonumber \\&\quad \left. S'^i_{k} \cdot \exp j \left[ 2 \pi (f_k+f) \left\{ {\Delta } \tau + {\Delta } \dot{\tau } (t-t_0)\right\} -A \cdot \frac{{\Delta } {\mathrm {TEC}} }{f_k + f} \right] \mathrm{d}f \mathrm{d}t \right| ,\nonumber \\ \end{aligned}$$

(36)

where $f_k$ is the center frequency of the band k and B is the bandwidth. The parameters of optimization then estimate the differential total electron content ${\Delta } {\mathrm {TEC}}$, the residual group delay ${\Delta } \tau$, and delay rate ${\Delta } \dot{\tau }$. To obtain the delay observable of WBWS, we add back the a priori delay. This can then be analyzed with standard VLBI delay analysis software, even though it represents a differential delay quantity with regard to the reference scan.