Super-Resolved Image of M87 Observed with East Asian VLBI Network

Abstract

Obtaining high-resolution images at centimeter-or-longer wavelengths is vital for understanding the physics of jets. We reconstructed images from the M87 22 GHz data observed with the East Asian VLBI Network (EAVN) by using the regularized maximum likelihood (RML) method, which is different from the conventional imaging method CLEAN. Consequently, a bright core and jet extending about 30 mas to the northwest were detected with a higher resolution than in the CLEAN image. The width of the jet was 0.5 mas at 0.3 mas from the core, consistent with the width measured in the 86 GHz image in the previous study. In addition, three ridges were able to be detected at around 8 mas from the core, even though the peak-to-peak separation was only 1.0 mas. This indicates that the RML image’s spatial resolution is at least 30% higher than that of the CLEAN image. This study is an important step for future multi-frequency and high-cadence observations of the EAVN to discuss the more detailed structure of the jet and its time variability.

1. Introduction

M87, located at 16.8 Mpc from Earth in the constellation Virgo, is a giant elliptical galaxy with a super-massive black hole of $6.5\times {10}^9{M}_{\odot }$ at its center [1]. This close proximity to a huge black hole allows for a linear resolution of 1 milliarcsecond (mas) = 0.08 pc = 130 Schwarzschild radii ($R_{\mathrm{s}}$). The relativistic jet erupts from the central core, emitting radiation ranging from radio to X-rays and gamma rays [2]. Even considering only the radio band, the jet’s appearance is quite different depending on the wavelength. For VLBI at centimeter-or-longer wavelengths, M87 jets extending to ∼900 mas at 150 mm [3] and ∼20 mas at 13 mm [4] have been observed; however, as the wavelength becomes shorter, the downstream of the jet becomes less visible, extending to ∼10 mas at 7 mm [5] and ∼3 mas at 3.5 mm [6,7]. The 1.3 mm data observed with the EHT imaged the ring-like structure in the very vicinity of the black hole [1]. The extended jet at the 1.3 mm waveband is so faint that it was not detected in the 2017 EHT observation. However, by improving short-baseline coverage, the dynamic range of the image can be increased, and the faint, extended jet can be recovered. How far downstream the jet can be detected depends on the sensitivity of the telescope, but there is no doubt that observations at centimeter-or-longer wavelengths are still more suitable for studying the downstream of the jet.

It is more difficult to obtain high-resolution images downstream of a jet than upstream because shorter wavelengths have a higher spatial resolution for the same telescope aperture. One method for obtaining high-resolution images of the downstream of a jet is space VLBI. By connecting ground and space telescopes, VLBI can be constructed to obtain a long baseline length. This effectively means that it plays the same role as a giant telescope that is larger than the size of the Earth. For example, the VLBI Space Observatory Programme (VSOP) satellite, led by the Institute of Space and Astronautical Science (ISAS) in JAXA in collaboration with the National Astronomical Observatory of Japan (NAOJ), was launched in 1997 and renamed HALCA, which carries a radio telescope with an 8-meter diameter [8]. High spatial resolution observations at 1.6 GHz and 5 GHz using VSOP successfully resolved the jet into three ridges at the mas scale [9]. Investigating how the complex internal structure of the jet evolves over time is critical to understanding the jet’s physics. However, as space VLBI observation time is limited, it is difficult to conduct studies that closely monitor the jet and reveal component motions in detail.

Therefore, to image the jet’s fine structure further from the central region, it is necessary to observe the jet at centimeter-or-longer wavelengths, where the jet is the brightest, and to reconstruct a high-resolution image. The regularized maximum likelihood (RML) method applied in this study is an imaging technique for interferometric data that was primarily developed for imaging the data obtained with the Event Horizon Telescope (EHT). While CLEAN, a conventional method, creates an image model after constructing a dirty map, resulting in images with resolutions limited by the beam size, RML directly constructs an image model that fits the observed data, resulting in higher-resolution images than the nominal diffraction limit determined by the uv coverage. As the ideal spatial resolution of VLBI is several times smaller than beam size [10], the use of RML may yield images with several-times-higher spatial resolutions than CLEAN images.

The East Asian VLBI Network (EAVN) is a VLBI system in East Asia, and currently consists of up to 16 telescopes with a longest baseline of 5078 km [11]. Three science working groups (SWGs), active galactic nuclei (AGNs), evolved stars, and star-forming regions are operated in the EAVN, of which M87 has continued to be monitored since 2017 at 13 and 7 mm wavelengths [12] as a target of the AGN SWG. This is a major monitoring project that has been carried out since the era of the KVN and VERA Array (KaVA; VLBI network in Japan and Korea; [13]), and the results of KaVA’s M87 observations are summarized in [14], which found where the jet accelerated and transitioned from subluminal to superluminal speeds. These high-cadence monitoring observations provide precise measurements of jet motion, which is an important step toward understanding the physics and evolution of jets.

2. Observations and Data Reductions

The data treated in this paper were observed with the EAVN on 18 March 2017, and have already been published in [12]. The paper [12] presents observations of four AGNs in two frequency bands of the EAVN, namely, 22 GHz and 43 GHz. They used the Astronomical Image Processing System (AIPS; [15]) provided by the National Radio Astronomy Observatory (NRAO) for an initial calibration of the complex visibility. After that, CLEAN imaging and self-calibration were performed from the calibrated data using DIFMAP software [16]. On the other hand, the data treated in this paper are only the M87 22 GHz observation, which is imaged with a different imaging method (see Section 3 for the details) using the initial calibrated data prior to imaging. In this section, a brief overview of the observational data used in this study is given.

The telescopes participating in this observation were the KVN and VERA Array (KaVA), Tianma 65 m Radio Telescope (TMRT), Nanshan 26 m Radio Telescope (NSRT), and Hitachi (HIT). VERA consists of four stations in Japan, namely, Mizusawa (MIZ), Iriki (IRK), Ogasawara (OGA), and Ishigakijima (ISG). KVN stands for Korean VLBI Network, which has three stations, namely, Yonsei (KYS), Ulsan (KUS), and Tamna (KTN), though KYS was unable to participate in this observation due to an issue at the site. See Table 1 in [12] for specifications of each participating telescope. Figure 1 shows that NSRT contributes to filling the outer part of the uv plane. The longest baseline is 5100 km between NSRT and Ogasawara, which corresponds to the smallest synthetic beam size of 0.6 mas in the northwest–southeast direction.

M87 was observed along with 3C 273, 1219 + 044, and M84. Within the 7-hour session, M87 was observed for 6 scans at 47 min per scan. The recording rate was 1 Gbps (2-bit sampling), where a total bandwidth of 256 MHz was divided into eight 32 MHz intermediate-frequency (IF) bands. Only left-hand circular polarization was received. All the data were correlated at the Daejeon hardware correlator installed at the Korea Astronomy and Space Science Institute (KASI). The initial calibration of visibility amplitude, bandpass, and phase was performed in the standard manner by using the AIPS software package.

3. Imaging

We reconstructed M87 images from the EAVN data set with the regularized maximum likelihood (RML) method implemented in the Sparse Modeling Imaging Library for Interferometry (SMILI; [17,18]), which is developed primarily for imaging EHT data [19]. The VLBI observation fills the uv-plane sampling depending on the telescope placement and the position of the source, as shown in Figure 1. To obtain an image from these data, the conventional method CLEAN draws a dirty map by inverse Fourier transform with zeros in the empty uv plane and reconstructs a set of point-source models from the dirty map. RML, on the other hand, uses regularization to ensure a plausible image consistent with the data without performing an inverse Fourier transform. In this imaging, we use the regularization terms of the weighted-L1 (wL1), total variation (TV), total squared variation (TSV), and maximum entropy method (MEM) (see Appendix A in [19] for the definition of each term), which describe assumptions about the source structure, such as that only limited regions in the field of view have brightness or that the source structure is smooth. Each term contains a hyperparameter, which adjusts the relative weighting of the regularization term to the data. If the weight of the regularization term is too strong, the image will be inconsistent with the data. On the other hand, if the weight of the data is too strong, the image will not reflect the features of the image expressed by the regularization term. This allows the production of dirty beam-free images with RML, which can recover finer structures than CLEAN.

We used visibility amplitudes, log closure amplitudes, and closure phases for image reconstruction. The visibility amplitude of TMRT shows some offsets compared with that of KaVA. Due to the uncertainty of a priori calibration, the visibility amplitudes of NSRT and HIT are very low compared with KaVA. These systematic errors in visibility amplitude do not affect the log closure amplitude because they are offset in the process of calculating the log closure amplitude. Therefore, the baselines, including TMRT, NSRT, and HIT, are excluded only from the visibility amplitude for the first imaging, and they are applied after the visibility amplitude is corrected by self-calibration.

wL1 regularization applies a pixel intensity penalty so that the dark areas in a prior image are also dark in the restored image, reducing the noise of the background region. The prior image is obtained by convolving the CLEAN map with 1.5 mas circular Gaussian. Since only loose constraints on the intensity distribution of the core and jet are required, we convolved the map with a large-enough circular Gaussian; therefore, our images are not affected by the detailed CLEAN model. We applied a total flux of 2.0 Jy based on the CLEAN results. The image properties are set to a pixel size of 80 $\mathsf{\mu }$as and a field of view of 41 mas × 20 mas, referring to Figure 9 of [12]. In addition, the image window, in which the intensity is calculated in imaging, was set as shown by the white circles in Figure 2. To remove noise within the image window and outside of the source structure, the brightness outside of the yellow-circles region in Figure 2 was corrected to zero at the end of each image. As the jet’s structure is restored in an area much smaller than these windows, it can be assumed that the window settings are appropriate and that the subsequent image reconstruction is successful.

In the first step, after excluding the HIT, NSRT, and TMRT baselines, we added a 200% error to the visibility amplitudes to reduce the weight of the visibility amplitudes and reconstruct images with a greater emphasis on the closure quantities information. We started the first imaging with the initial image of a circular Gaussian with the full width at a half maximum of 0.1 mas. The regularization parameters of wL1, TV, TSV, and MEM were set to 1, 10, 10, and 0.0001, respectively. We performed the first self-calibration by using the obtained image as a model and obtained a calibrated uv data set to proceed to the next step of the parameter survey. Self-calibration restored the depressed visibility amplitudes, and all baselines were used hereafter.

In the next step, iterative imaging was performed using the self-calibrated uv data. The iterative pipeline of imaging and self-calibration was created to investigate how the image changes or is consistent with the data depending on the parameter set. The following parameter combinations were prepared: additional errors of visibility amplitudes (err) = [0.1, 0.01, 0], regularization parameters of MEM (${\lambda }_{\mathrm{MEM}}$) = [0.01, 0.001, 0.0001], regularization parameters of wL1 (${\lambda }_{\mathrm{wL}1}$) = [10, 1, 0.1], and regularization parameters of TSV (${\lambda }_{\mathrm{TSV}}$) = [100, 10, 1]. After the iterative imaging of 81 combinations of the parameter sets, we selected the final images with good fits to the data, i.e., those with reduced ${\chi }^2$ around 1. The selection criteria of reduced ${\chi }^2$ for the closure phases and log closure amplitudes are less than 1.2 and 1.3, respectively.

4. Image Properties

In total, 48 final images, which fit the data better than the selection criteria, were selected among the 81 images. There are slight differences among the 48 selected images depending on the parameter set. For example, MEM had a particularly pronounced effect on the images, with the larger ${\lambda }_{\mathrm{MEM}}$ resulting in a blurrier image. Ex. 1, 2, and 3 in Figure 3 show images with ${\lambda }_{\mathrm{MEM}}$ set to 0.0001, 0.001, and 0.01, respectively, all with a reduced ${\chi }^2$ near to 1. The top image in Figure 3 is the average of 48 such images that are slightly different but meet the reduced ${\chi }^2$ selection criteria. All images have a bright central core from which a fainter jet extends in a northwesterly direction. Furthermore, the jet dims at about 15 mas from the core and brightens again at about 20 mas, which is consistent with the CLEAN image [12].

A counter-jet-like structure is also visible, though it is very faint and short compared with the approaching jet. We tried to calculate the noise level of the image at a sufficient distance from the source structure using DIFMAP, assuming that the pixel values of the averaged RML image are the CLEAN model. Using natural weighting, we estimated rms = 0.3 mJy/beam with the beam size (FWHM) of $1.39\times 0.605$ mas in 11.6 degrees. The counter-jet-like structure was about six times brighter than the noise rms and is considered to be significantly detected. The brightest part of the counter feature in the average RML image, which is located 1.7 mas southeast of the core, has about 23% of the brightness of the main jet symmetrically located relative to the core. This means that the brightness ratio (BR) of the main jet to the counter jet is about 4.3. However, there is no conspicuous blob to the northwest of the core that fully corresponds to the blob to the southeast. In previous studies, the BR of the M87 jet at 1 mas from the core was about 5–25 [6], and at 0.5–3.1 mas, 10–15 [20]. In [21], the BR integrated over a region was 14.4, but the peak-to-peak ratio was 3.4. Although this study is roughly consistent with the results of previous studies, more precise BR measurements are needed to estimate the jet velocity.

To evaluate the jet shape, all 48 images were processed as follows. We rotated all images 19 degrees clockwise so that the jet extended horizontally rightward from the bright core. The vertical profile of this image was obtained to analyze the intensity distribution in the direction perpendicular to the jet axis. The profiles were located 0.3 mas and 0.5 mas from the core, and then in 0.5 mas increments up to 8.0 mas, they were averaged by ±0.2 mas horizontally, respectively. The jet width was defined using the full width at half-maximum intensity (FWHM). If there are multiple peaks in the profile, the FWHM of the peaks at both ends and the peak-to-peak distance between them were used to calculate the width of the entire jet. Jet widths were measured in all 48 images, and the average values are plotted in Figure 4 with red diamond marks. The error bars are their standard deviations. The widths of 0.3 mas and 8 mas from the core were 0.5 mas and 2.7 mas, respectively. The jet width’s dependence on the core distance can be fitted with a power law with an exponent of $0.54\pm 0.09$ (95% confidence level). These widths are consistent with the jet width of the 43 and 86 GHz images obtained with the CLEAN method in [6], and they are shown with green and blue circle marks in Figure 4.

To see if we can resolve the ridge structure in the images, we investigated the profile perpendicular to the jet axis at a distance of 8 mas along the jet from the core (Figure 5), where we see three peaks in each image. The best-fit parameters were obtained by fitting the profile with three Gaussian components. The separation between the central and southern peaks was 1.1–1.4 mas, and the separation between the northern and central peaks was even narrower at 0.9–1.1 mas. On the other hand, no triple ridges were detected in the CLEAN image obtained using the same observation data [12]. As the beam size used for the convolution of the CLEAN image was 1.35 mas perpendicular to the jet, by assuming it to be the spatial resolution of the CLEAN image in this direction, the RML image with a peak separation of 1.0 mas shows at least 30% higher spatial resolution than the CLEAN image. The ridges with the brightest peak intensities were, in order, the northern, the southern, and the central ridges. The central ridge seems to be systematically thicker than the northern and southern ridges; however, for some images, the difference in width is so small that it is unclear whether there is a significant difference. Previous observations at 1.6 GHz with VSOP [9] and the ultra-deep image of 2 Gbps VLBA with the phased VLA at 15 GHz [22] detected a triple-ridge structure more than 5 mas downstream from the core. The triple ridge at 8 mas from the core detected in this study is also consistent in position and width with these, so it is likely that the same phenomenon was captured.

5. Summary

We reconstructed M87 images from the EAVN 22 GHz data set with the RML method implemented in the SMILI tool using sparse modeling techniques. Visibility amplitudes, log closure amplitudes, and closure phases are utilized for the image reconstruction. By setting various parameters in an iterative pipeline of imaging and self-calibration, 81 images were obtained. Finally, 48 images with a reduced ${\chi }^2$ of about 1 were selected as the final images. Although there are slight differences among the selected images, depending on the parameter set, all images have a bright central core from which a fainter jet extends in a northwesterly direction. A counter-jet-like structure is also visible with the detection at 6$\sigma$, although it is very faint and short compared with the approaching jet. We successfully measured the jet width of the region of a few mas at the root of the jet, which is consistent with the results of previous studies measured at 43 GHz and 86 GHz. We observed three peaks in the profile perpendicular to the jet axis at a distance of 8 mas from the core in each image, which is consistent with the previous studies from a space VLBI by VSOP and an ultra-deep image by the 2 Gbps VLBA observation with the phased VLA. The image of the M87 jet obtained by RML allows us to examine the differences in the brightness and thickness of the three ridges in detail, as well as their distance from each other, features which could be lost by convolving with nominal beam sizes, as in the CLEAN method.

In the future, we must image the M87 jet at a centimeter-or-longer wavelength with higher resolution and sensitivity to investigate the entire jet from the root to the downstream with more precise profiles. This will help us to understand the acceleration and collimation mechanisms of the jet. If we observe the same region by using multiple wavelengths, we should be able to determine the physical properties, such as the optical depth and magnetic field. Furthermore, if monitoring observations are made at a resolution high enough to distinguish the three ridges, it is possible to study the time variation of the peculiar structure. Expanding the VLBI array will be one of the keys to achieving these goals. An attempt to create a global array centered on East Asia has already begun with EATING VLBI, a joint effort between the EAVN and the Italian VLBI Network [23]. Collaboration between East Asia and Australia has also started as a joint observation between the Tidbinbilla 70 m radio telescope and the EAVN, and it will include other Australian stations in the future. In addition, a radio telescope in Thailand is scheduled to join the EAVN in the future. These array extensions will play a pilot role for the next generation of VLBI, such as SKA-VLBI [24] and ngVLA [25].