Pilot KaVA monitoring on the M87 jet: confirming the inner jet structure and superluminal motions at sub-pc scales

TL;DR

This study used the KaVA array to monitor the M87 jet at 22 GHz, revealing detailed inner jet structure, superluminal motions, and acceleration, confirming these features are fundamental to the jet’s behavior.

Contribution

First high-cadence KaVA monitoring of M87 jet at 22 GHz, demonstrating its capability to observe superluminal motions and acceleration at sub-parsec scales.

Findings

Detected superluminal motions at 0.1-2 pc scales.

Observed gradual acceleration of jet components.

Confirmed similar kinematic features across different instruments.

Abstract

We report the initial results of our high-cadence monitoring program on the radio jet in the active galaxy M87, obtained by the KVN and VERA Array (KaVA) at 22 GHz. This is a pilot study that preceded a larger KaVA-M87 monitoring program, which is currently ongoing. The pilot monitoring was mostly performed every two to three weeks from December 2013 to June 2014, at a recording rate of 1 Gbps, obtaining the data for a total of 10 epochs. We successfully obtained a sequence of good quality radio maps that revealed the rich structure of this jet from <~1 mas to 20 mas, corresponding to physical scales (projected) of ~0.1-2 pc (or ~140-2800 Schwarzschild radii). We detected superluminal motions at these scales, together with a trend of gradual acceleration. The first evidence for such fast motions and acceleration near the jet base were obtained from recent VLBA studies at 43 GHz, and the fact that very similar kinematics are seen at a different frequency and time with a different instrument suggests these properties are fundamental characteristics of this jet. This pilot program demonstrates that KaVA is a powerful VLBI array for studying the detailed structural evolution of the M87 jet and also other relativistic jets.

Full text

1 Introduction

Very long baseline interferometry (VLBI) is a powerful technique for probing compact nonthermal radio sources and phenomena in the Universe. In particular, VLBI plays a major role in studying relativistic jets in active galactic nuclei (AGN), since the ultra-high-resolution capability allows direct imaging of the sites of e.g., collimation, acceleration and production of high-energy $\gamma$-ray emission in the jets, which generally take place at pc-to-subpc scales or even less. Over the past decades, tremendous observational efforts have been made with various VLBI arrays around the world (e.g., VLBA, EVN and LBA) to elucidate the relativistic jets. Nevertheless, a number of fundamental questions on jet physics still remain to be answered, requiring further imaging and monitoring of the jet-formation scales with high-resolution VLBI.

The KVN and VERA Array (KaVA) has recently emerged as the first regularly operated international VLBI network in East Asia, consisting of three 21-meter dishes in Korea (KVN; e.g., [Lee et al. (2014)]) and four 20-meter dishes in Japan (VERA; [Kobayashi et al. (2003)]). Over the last years, we made efforts to combine the two arrays in collaboration among the Korea Astronomy and Space Science Institute (KASI), National Astronomical Observatory of Japan (NAOJ) and associated universities. The aim of this project is to form a single higher performance array by complementing each other. Indeed, the increases in the number of baselines (from 3/6 for KVN/VERA to 21 for KaVA) and the baseline coverage (from 300-400 km/1000-2300 km for KVN/VERA to 300-2300 km for KaVA) significantly improve the overall array sensitivity and the imaging performance, as first demonstrated by Matsumoto et al. (2014) and Niinuma et al. (2014).

Following a commissioning phase lasting several years, partial (risk-shared) regular operation of KaVA began in 2014. The common frequency bands of 22 GHz and 43 GHz are available, and a data-recording rate of 1 Gbps is currently supported. All data observed by KaVA are correlated at the Korea-Japan Joint VLBI Correlator installed at KASI (the so-called Daejeon hardware correlator; Oh et al. (2012); Lee et al. (2015a)). A feature of KaVA is that the array operates for a quasi-full year, except for the annual (6–8)-week maintenance period that starts from mid-June. This offers an ideal opportunity to monitor structural evolution of relativistic jets via detailed multi-epoch observations.

While some early KaVA monitoring results on bright blazars have already been published (Oh et al., 2015; An et al., 2016), here we focus on the radio galaxy M87. M87 is one of the nearest AGN with a powerful relativistic jet. Due to its proximity (16.7 Mpc; Blakeslee et al. (2009)) and the huge mass of the central black hole ((3–6)$\times 10^{9}M_{\odot}$; Gebhardt et al. (2011); Walsh et al. (2013)), this source offers an exceptional opportunity to resolve the jet-formation scales with VLBI (e.g., Junor et al. (1999); Hada et al. (2011)). Recent high-sensitivity VLBI observations at 86 GHz revealed a limb-brightened, wide-opening-angle jet launching at a scale of $\lesssim$10 $R_{\rm s}$ (Hada et al., 2016; Kim et al., 2016). Moreover, extensive analyses of multi-frequency VLBI images revealed a parabolic collimation over a distance from z$$\sim100 $R_{\rm s}$ to z$$\sim$$10^{5} $R_{\rm s}$ (Asada & Nakamura, 2012; Hada et al., 2013). These results suggest that magneto-hydrodynamic (MHD) processes play a key role in collimating the M87 jet (e.g., Nakamura & Asada (2013)).

Compared to the progress in morphology analysis, the structural evolution of this jet is still controversial. Early multi-epoch VLBI studies of M87 (at sparse intervals of months-years) often reported slow or quasi-stationary motions at the pc-to-subpc scales (Reid et al., 1989; Dodson et al., 2006; Ly et al., 2007; Kovalev et al., 2007). However, such low speeds would be insufficient to explain the known large jet-to-counter-jet brightness ratio at the same scales (typically $>$10). To date, the most detailed program was performed by Walker et al. in 2007–2008 with VLBA at 43 GHz (Walker et al., 2008, 2016), where M87 was monitored every three weeks from 2007 Jan-Aug and every $\sim$five days from 2008 Jan-Apr. In contrast to the previous knowledge, they first noted the existence of fast ($\gtrsim$(1–2) $c$) motions near the jet base. Using the same data, the recent sophisticated image analysis by Mertens et al. (2016) have further suggested the presence of flow acceleration at these scales. Nevertheless, it is still unclear whether such fast motions were just a temporal event or not, and also whether similar motions can be seen at other frequencies or not.

Therefore, there is a growing consensus on the necessity of a high-cadence (every a few weeks or less), continuous monitoring program to properly understand the kinematics of the M87 jet. Unfortunately, the VLBA’s massive program is triggered only occasionally along with unpredictable TeV events, preventing a monitoring of the source on a regular basis. Since M87 is one of the few jet sources that gravitational scales of $\ll 10^{3-4}\,R_{\rm s}$ are directly resolvable, having a regular monitoring program should be of great value for testing the magnetically driven jet paradigm by comparing with theories, as well as rapid follow-up of high-energy events. This is our primary motivation for launching a new monitoring program of this jet with KaVA.

In this paper we report on the initial results from our first KaVA–M87 monitoring program. This is a pilot study conducted at the early stage of our KaVA operation, and thus the present paper primarily focuses on the validation of KaVA’s imaging and monitoring performance on this jet. More detailed analyses (including newer data) and comparisons with some specific theoretical models are currently prepared in separate papers (Park et al. in prep; Ro et al. in prep.). In the next section we describe our observations and data reduction. In Sections 3 and 4, the results on imaging and kinematics are presented. In the final section we will summarize the paper and describe our future prospects on this project. Throughout the paper we adopt $D=16.7\,{\rm Mpc}$ and $M_{\rm BH}=6.0\times 10^{9}\,M_{\odot}$, corresponding to $1\,{\rm mas}=0.08\,{\rm pc}=140\,R_{\rm s}$.

2 Observations and Data Reduction

From December 2013 to June 2014 M87 was densely monitored with KaVA at 22 GHz. The period was in a transition phase between internal commissioning and regular science operation of KaVA, so the data used here are the combined results of these two stages. The monitoring was made at sampling intervals of two to three weeks. The observations were largely successful, but two particular sessions performed on 2016 February 3 and 25 were lost due to recording issues. In addition, the data taken on 2014 May 15 was poor in quality due to the loss of VERA-Ishigaki and severe weather conditions at VERA-Mizusawa, which prevented us from obtaining a reliable jet image. As summarized in Table 1, we eventually obtained our dataset for a total 10 epochs during the seven months. For the first three epochs (Dec/05, Dec/25, Jan/15), the on-source time of M87 was relatively short because the source was intermittently scanned in a multi-source survey program, while for the rest of the epochs M87 was observed with a quasi-continuous mode. Weather conditions were generally stable at the KVN and VERA-Mizusawa/Iriki stations, while the weather at the two southern stations VERA-Ishigaki/Ogasawara was a little more variable (e.g., rain showers). Left-hand circular polarization was received and sampled with a 2-bit quantization. At the VERA stations, a single-beam mode was used in KaVA observations. All the data were recorded at 1 Gbps (256 MHz bandwidth, 16 MHz $\times$ 16 subbands) and correlated by the Daejeon correlator.

The initial data calibration was performed using the National Radio Astronomy Observatory (NRAO) Astronomical Image Processing System (AIPS) based on the standard VLBI data reduction procedures. A priori amplitude calibration was applied using the measured system noise temperature and the elevation-gain curve of each antenna. We calibrated the amplitude part of bandpass characteristics at each station using the auto-correlation data. Following the amplitude calibration, fringe-fitting was performed to calibrate the visibility phases, and finally the data were averaged over each sub-band. Imaging and self-calibration were performed in the Difmap software (Shepherd et al., 1994) in the usual manner.

For KVN data correlated with the Daejeon correlator, there is a known “amplitude-loss” problem due to the accumulation of multiple losses through signal processing (Lee et al., 2015b). For such KVN data Lee et al. (2015b) specified an amplitude correction factor of 1.35 (for the data correlated before 2015 March) or 1.30 (for the data correlated after 2015 March). As for KaVA data, the corresponding values are currently being examined (Oh et al. in prep), but we may reasonably assume the same correction factors, since the backend system of VERA is similar to that of KVN. Therefore in the present analysis of our KaVA data, we rescaled the amplitude of the calibrated visibilities by a factor of 1.35. In Section 3.2 we show that this factor is indeed reasonable.

3 Results and Discussion

3.1 Images

The first KaVA image of M87 at 22 GHz was presented in Niinuma et al. (2014), where the observation was made at 128 Mbps with a relatively short ($\sim$60 min) integration time. While the early KaVA image nicely demonstrated the improvement of image quality over the only-VERA case, the image dynamic range was still limited to $\sim$1000 and the jet emission was barely detected down to $\lesssim$10 mas from the core.

To describe the typical imaging quality of our present 1 Gbps KaVA dataset, we show two representative images of M87 in Fig. 1, which are the same data taken on 2014/Mar/2 but different uv-weighting and restoring beams are applied.

The naturally weighted image (the top panel of Fig. 1) illustrates the overall structure of the M87 jet. Thanks to the improvement of our data-recording rate, KaVA detected the extended jet emission down to $\sim$20 mas very clearly (at SNR much higher than 5 $\sigma$). The resulting dynamic range of our KaVA 1 Gbps images was 1600-4000, depending on the number of stations, weather condition and on-source time.

The uniformly weighted image (the bottom panel of Fig. 1) improves the angular resolution by $\sim$20% ($\sim$1 mas on average), resolving the rich structure of the inner jet regions. The well-known limb-brightened structure was clearly resolved up to $\sim$2 mas near the core. With this image we measured the peak-to-peak width of the limb-brightened jet as a function of distance. The result is shown in Fig. 2, together with a contemporaneous VLBA 24 GHz result presented in Hada et al. (2016). The KaVA result is in good agreement with the VLBA one in the overlapping regime, and the slope is also consistent with the previously known parabolic collimation profile.

Note that M87 is known to have the counter jet (typically a-few-mas long) at the eastern side of the core (Ly et al., 2007; Kovalev et al., 2007). For the present KaVA 22 GHz data, we were also required to put some CLEAN components at the counter side through our deconvolution process, but the extent was limited only to $\lesssim$1 mas. However. the contemporaneous higher-resolution VLBA images at 86/43/24 GHz (Hada et al., 2016) found an exceptionally small extent of the counter jet around this period (visible only $\sim$0.5 mas from the core). Thus the bare detection of the counter jet with KaVA is consistent with the VLBA results.

3.2 Light curves

Investigating a detailed radio light curve of the core or the innermost jet of M87 gives us important insight into high-energy activities near the black hole (Acciari et al., 2009; Hada et al., 2012, 2014; Akiyama et al., 2015; Lee et al., 2016). To check the feasibility of studying radio variability with KaVA, in Fig. 3 we show light curves for the M87 jet obtained through this program. Here we present two light curves describing different spatial scales: one is for the peak flux of the core when convolved with a 1-mas-circular beam, while the other is for the total CLEANed flux including the extended jet. In addition, we also plot available contemporaneous VLBA 24 GHz light curves (the data presented in Hada et al. (2016)) for the corresponding regions. For each data 10% uncertainty is assumed in amplitude.

For both the compact and the extended scales, the KaVA and VLBA flux densities are in good agreement with each other. For the core, the smooth flux variation is well sampled with KaVA especially for the period of continuous triweekly monitoring (2014 Mar–Jun). On the other hand, the measurement of integrated flux is generally more sensitive to the overall image sensitivity and uv-coverage. Nevertheless, thanks to the KVN baselines as short as 17 M$\lambda$ (almost equivalent to the VLBA’s shortest baselines 15 M$\lambda$), KaVA consistently recovered the M87 extended emission (well within 10% uncertainty) that was imaged by the VLBA.

Therefore, these results demonstrate that the KaVA array is also very suitable for studying light curves of relativistic jets on various spatial scales.

3.3 Structural evolution of the jet

The structural evolution of the M87 jet is one of the major motivations of our monitoring program. On the left hand side of Fig. 4, we show multi-epoch CLEAN images of M87 between January and June 2014 (all the images are rotated by 20∘ clockwise). While more dedicated studies on the kinematics and structural evolution including two-dimensional velocity field are currently prepared in other papers (Park et al. in prep.; Ro et al. in prep.), here we focus only on the radial structural evolution of the jet based on a rather simple analysis, which can be quickly compared with various previous studies.

Since the M87 jet has a limb-brightened and smooth intensity distribution, rigorous component identification in the jet is non trivial (e.g., Mertens et al. (2016)). Nevertheless, here we applied the conventional visibility-based model-fitting technique (with modelfit in Difmap), and reconstructed the M87 jet with a set of circular Gaussian components. These images are shown on the right hand side of Fig. 4. While this simple modeling often failed to recover the weak extended emission at $\sim$20 mas, the inner jet within 10 mas was reasonably characterized. This is supported by the fact that the number of fitted Gaussians and the derived parameters (location, size, flux) of each component in this region were remarkably similar among all the densely sampled images, indicating that the jet does have some “characteristic patterns” that can be commonly/consistently identified among different epochs. Of course the actual jet structure is much more complex than modeled here, but such a detailed modeling is beyond the scope of the present paper. A more complete modeling and characterization of this jet will be presented in a forthcoming paper (Park et al. in prep.)

For each Gaussian component, the time evolution of the separation from the core is plotted in Fig. 5. Note that we grouped two Gaussian components at $\approx$3–5 mas from the core and treated them as a single component (component D), estimating the flux-density-weighted component position. This grouping was made because (i) there is only a single Gaussian component in the last epoch at a similar location and (ii) their sizes and distribution vary a lot in different epochs but the grouped positions show a gradual outward motion and are smoothly connected to the single Gaussian component in the last epoch (see the left panel of Fig. 4). Thus the component D represents a spatially averaged motion around the 3–5 mas region. The densely sampled dataset confirmed that all the components were moving outward. We derived an apparent velocity of each component by fitting a linear function to each dataset, and obtained a wide range of values from $\sim$0.16 $c$ to $\sim$1.60 $c$. These values are then plotted in Fig. 6, as a function of angular separation from the black hole, by using the result of the previous core-shift measurement by Hada et al. (2011)111The result of Hada et al. (2011) indicates that the convergence point of the core-shift is located at 70 $\mu$as upstream of the 22 GHz radio core. Here we assume that the convergence point is identical to the black hole location and the relative separation between the black hole and the 22 GHz core is stationary with time..

As seen in Fig. 6, the apparent velocity profile obtained by KaVA suggests a trend of acceleration within 10 mas from the black hole, showing a transition from sub-to-super-luminal motion around $\sim$5 mas. For comparison, we also plot previous proper motion results at the same scales from the literature. The KaVA result is largely in agreement with those reported by Ly et al. (2007), Walker et al. (2008), Mertens et al. (2016), Hada et al. (2016) and Walker et al. (2016), where they typically found slow motions at $\lesssim$1 mas and fast motions beyond $\sim$1 mas. As an additional test to check the validity of our measurement, we also applied the Wavelet-based Image Segmentation and Evaluation (WISE) technique (Mertens & Lobanov, 2015; Mertens et al., 2016) to our KaVA data, and obtained a consistent velocity profile with that of our model fitting (see Appendix for more detail). This is additional support that our simple method is valid and the acceleration and superluminal motions obtained by KaVA are indeed real.

Note that, besides the moving components, our multi-epoch KaVA images suggest the presence of a quasi-stationary component at $\sim$20 mas from the jet base. Interestingly, the M87 jet has a bright standing component at this distance since 1980’s (Reid et al., 1989; Dodson et al., 2006), implying a fundamentally different origin from the temporal moving features. We are investigating this component in more detail, and a complete analysis will be presented in a forthcoming paper (Park et al. in prep.).

In Fig. 7, we show an extended plot of M87 jet kinematics including large-scale jet features studied in the literature (Biretta et al. (1999); Cheung et al. (2007); Giroletti et al. (2012); Meyer et al. (2013); Asada et al. (2014)). Adopting a viewing angle of 15 deg (e.g., Perlman et al. (2011)) and a BH mass of $6\times 10^{9}M_{\odot}$, the vertical and horizontal axes are described by Lorentz factor and deprojected distance from BH in $R_{\rm s}$ unit, respectively. The proper motion profile measured by KaVA ($v_{\rm app}\propto z^{0.86\pm 0.10}$) is translated into $\Gamma\propto z^{0.22\pm 0.03}$. Interestingly, this relation can consistently reproduce the maximum $\Gamma$ at HST-1 within $\approx 1\sigma$, implying that a single power-law acceleration profile may hold over three orders of magnitude in distance. Note that recent jet-shape studies of M87 at the same scales have independently discovered a collimation profile with a single power-law $z\propto r^{1.7}$ (Asada & Nakamura (2012); Nakamura & Asada (2013); Hada et al. (2013)). These results are consistent with the essence of the magnetically driven jet scenario, where the collimation and the acceleration are co-spatial (Komissarov et al., 2009).

Of course the jet velocity profile obtained here is still rather coarse, and the exact profile may be more complex than suggested here, in fact as proposed by other detailed studies (Asada et al., 2014; Mertens et al., 2016). Nevertheless, the important thing here is that this pilot program demonstrates that KaVA has excellent potential to study the velocity field and structural evolution of this jet in further detail. More extensive and advanced studies on this jet will be presented in the framework of the KaVA Large Program (as described below).

4 Summary and Future Prospects

In this paper we have reported the basic performance of KaVA imaging and monitoring capability for M87 based on a pilot multi-epoch program at 22 GHz. We have shown that KaVA 1 Gbps observations of M87 at 22 GHz can well image the extended jet structure down to $\sim$20 mas from the core. An image dynamic range of $\sim$4000 was achieved, but further improvement of the KaVA imaging performance is expected depending on weather condition and integration time. The limb-brightened, parabolic jet structure at pc-to-subpc scales can also be imaged in agreement with VLBA images. This indicates that we can deeply study the morphology of this jet from a few 100 $R_{\rm s}$ to $\sim$6000 $R_{\rm s}$ with KaVA (for a viewing angle of 15*∘*).

The densely sampled multi-epoch KaVA images have traced the detailed structural evolution of this jet. In agreement with the recent finding based on the VLBA 43 GHz data by Mertens et al. (2016), our KaVA 22 GHz data have also detected the superluminal outward motions together with linear acceleration within 10 mas from the jet base. The independent confirmation at different time, at different frequency, and with different instruments suggests that the rapid jet acceleration near the jet base is not episodic but a fundamental property of the M87 jet. The fast jet speeds near the jet base may consistently explain why the counter jet of M87 is so weak.

Finally, below we briefly describe the future prospects on our M87 program.

The KaVA Large Program: While the kinematics result presented here is still rather coarse, this pilot study demonstrated that the KaVA has excellent potential to reveal the velocity profile of this jet in more detail. Following this, since 2016 we have upgraded our KaVA monitoring program of M87 in the framework of the KaVA Large Program (Kino et al., 2015). With the new program, every year we monitor M87 mostly biweekly over $>$6 months at both 22 and 43 GHz quasi-simultaneously. Depending on the progress of our ongoing KaVA system upgrade, we will conduct observations at a higher recording rate ($\gtrsim$2 Gbps). This strategy will greatly improve the quality of our proper motion analysis with KaVA and set a much stronger constraint on the velocity field of the M87 jet. More detailed studies about the jet structural evolution (both radially and transversally) are now actively ongoing by including the newer data. Also, the dual-frequency multi-epoch data will enable us to derive a set of accurate spectral index maps and their detailed evolution with time, which will provide important insights into the underlying particle energetics, density and magnetic fields. Furthermore, monitoring the light curve of the radio core on various time scales (from days to years) may tell us a hint on the accretion rate onto the black hole (Park et al., 2017), a key parameter to understand the jet production.

The East Asian VLBI Network: Concurrently, the international VLBI collaboration in East Asia is rapidly growing. Besides the successful regular operation of the KVN and VERA array, the KaVA has started joint experiments with telescopes in China as well as some more stations in Korea and JVN. This will ultimately form the East Asian VLBI Network (EAVN), a huge VLBI array with up to $\sim$20 stations distributed throughout these countries (Wajima et al., 2016). The commissioning of EAVN is currently actively ongoing. In particular, in 2017, a part of our KaVA Large Program on M87 is expanded to invite some more stations: this includes the Tianma (Shanghai) 65 m (22/43 GHz), the Urumqi 25 m (22 GHz), the Nobeyama 45 m (43 GHz), the Takahagi 32 m (22 GHz), the Hitachi 32 m (22 GHz), the Kashima 34 m (22 GHz) and the Sejong 21 m (22/43 GHz). This will tremendously enhance the overall performance of our VLBI array such as angular resolution, sensitivity, uv-coverage and image dynamic range. With such a powerful capability of EAVN, we will be able to probe the acceleration and collimation properties of M87 over the whole distance from the jet base to HST-1 in unprecedented detail.

{ack}

We thank the anonymous referee for his/her careful review and suggestions that improved the manuscript. We acknowledge all staff members and students at KVN and VERA who supported the operation of the array and the correlation of the data. KVN is a facility operated at by the Korea Astronomy and Space Science Institute. VERA is a facility operated at National Astronomical Observatory of Japan in collaboration with associated universities in Japan. K.H. is supported by the Research Fellowship from the Japan Society of the Promotion of Science (JSPS). This work was partially supported by KAKENHI (26800109 and 15H00784). TJ and GYZ are supported by Korea Research Fellowship Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (NRF-2015H1D3A1066561). I.C. is financially supported by the National Research Foundation of Korea (NRF) via Global PhD Fellowship Grant (NRF-2015H1A2A1033752). Part of this work was achieved using the grant of Visiting Scholar Program supported by the Research Coordination Committee, National Astronomical Observatory of Japan (NAOJ).

Appendix: Application of the WISE technique to KaVA data

Mertens et al. (2016) showed that the detailed velocity field of the M87 jet can be investigated using the Wavelet-based Image Segmentation and Evaluation (WISE) method. This method provides an un-supervised scheme to detect and identify various characteristic structural patterns in different epochs and measure the jet velocities. We applied the WISE method to our data set and obtained velocity vectors at different regions of the jet on four spatial (decomposition) scales, 0.2, 0.3, 0.4, and 0.6 mas (Fig. 8, see e.g., Mertens & Lobanov (2015) for details)222We followed the command line tutorials for 3C 120 provided on the Web site; http://flomertens.github.io/wise/tutorials_cmd/Walkthrough3C120.html. Most of the velocity fields showed radial, outward motions except quasi-stationary ones at $\approx 0$ and $\approx 10$ mas from the core on all the four scales (green and red arrows, respectively), and another one at $\approx 3$ mas from the core on a scale of 0.6 mas (cyan arrow) showing back-and-forth motions. This behavior may reflect the complex nature in the 3–5 mas regions, as we needed a component grouping in this area in our modelfitting approach (see Section 3.3). In addition, we detected at most four to six velocity vectors only at specific locations in the jet. These numbers are quite smaller than those presented in Mertens et al. (2016). These are likely due to the limited angular resolution of KaVA, not to the limited sampling of our monitoring since the average interval between adjacent epochs is $\approx$ 19 days, which is comparable to that of Mertens et al. (2016).

In Fig. 9 we present the apparent velocities of the radial, outward motions derived from the WISE technique. Note that not all these values are independent of each other. The values obtained from the WISE method show a consistent acceleration trend with our modelfit analysis, although they seem to be slightly below our best-fit line. This small difference might be due to an under-estimation of the velocities in the WISE method or an over-estimation in our modelfit analysis. This could also be investigated in our future studies using the Large Program data in more detail (e.g., Park et al. in prep.).