The jet and resolved features of the central supermassive black hole
of M 87 observed with EHT

In the analysis of VLBI data, the hybrid mapping method is widely used to obtain a calibration solution for the data and to reconstruct the brightness distribution. Hybrid mapping, which consists of repeatedly assuming one image model, performing self-calibration, obtaining a trial solution for calibration, and improving the image model for the next self-calibration, is the only method that is essential for precise calibration of VLBI data (Pearson & Readhead, 1984; Readhead & Wilkinson, 1978; Schwab, 1980). VLBI systems are not so stable in phase and amplitude as connected radio interferometers. In addition, millimeter- and submillimeter-wave observations are more affected by atmospheric variations. Therefore, the hybrid mapping method is becoming more and more important in the calibration of high frequency VLBI data such as the EHT observations. We performed a standard hybrid mapping process using the tasks CALIB and IMAGR in AIPS (the NRAO Astronomical Image Processing System ⁸⁸8http://www.aips.nrao.edu/index.shtm, Greisen (2003)).

After obtaining the image models for the two points, more than 100 iterations, including trials and errors, were performed in the hybrid mapping process. Most of the CLEAN images were run with the parameters listed in Table 4.

As described in Event Horizon Telescope collaboration (2019d), care must be taken in the choice of FOV, as incorrect restrictions will result in incorrect image structures. Considering the well-known structure of M 87, we restricted the imaging region by eight BOXes where emission could be detected. For the self-calibration, the selected CLEAN components were used as the next imaging model and the parameters in Table 1 were used. By repeating the phase-only self-calibration in this way, we were able to find better images and calibration solutions. This is because the method of simultaneously solving the amplitude solution with hybrid mapping has the risk of going in the wrong direction. Self-calibration of phase-only can be done with a certain degree of safety. Even if the wrong model is used and a completely wrong solution is obtained, the closed phase is automatically preserved and convergence to a completely wrong image can usually be avoided. On the other hand, the amplitude solution from self-calibration can be infinitely large or small by using the wrong image model. It is safe to try the self-calibration of the amplitude solution after a confident correct image model is obtained from the phase self-calibration.

In the latter half of the imaging process, we selected several candidates for the final image by comparing the difference in closure phase between the image and the real data. Furthermore, we created several image models using the CLEAN components of the candidate images to find the best image with the minimum closure residuals. Since we found that the source structure has time variation, we divided the data into the data of the first two days and the data of the last two days. As a result, we obtained the best images as shown in Figure 4. The data from the first two days was a composite of seven raw CLEAN maps, and the data from the last two days was a composite of nine raw CLEAN maps. Both consist of CLEAN components only. In the images of the first two days, adjacent CLEAN components within $2~{}\mu\rm as$ are merged into one component. In the images of the last two days, adjacent CLEAN components within $5~{}\mu\rm as$ are merged into one component.

Since the eight BOXes cover a large area, the resulting CLEAN component contains three types of emission: real emission, associated diffraction (sidelobes), and false acquisitions. For example, the bright emission on the right edge is not real. Such unreal bright spots often appear when the VLBI data is analyzed and imaged. When such unreal brightness appears, there may actually be strong emissions outside the BOXes. We produced a large image (30 arcsec field of view) using very short baseline (SM-JC and AA-AP) visibility data, but could not detect any new strong features.

To get a more complete image, we need to select only the real ones from these CLEAN components and perform CLEAN imaging again with each narrow BOX to cover the selected ones. However, since our data did not have enough u-v coverage and quality to select the correct ones among the CLEAN components, we gave up the task of extracting the CLEAN components this time. Nevertheless, the quality of the final image does not seem to be too bad. The closure residuals of the resulting image show a small variance comparable to the EHTC ring image. The residual of our map for the first two days of data is $-4.90~{}\pm~{}37.93^{\circ}$ (the residual of the EHTC ring is $-0.73~{}\pm~{}45.33^{\circ}$) and the residual of our map for the last two days of data is $1.79~{}\pm~{}42.11^{\circ}$ (the residual of the EHTC ring is $4.22~{}\pm~{}45.74^{\circ}$). Here, we integrated the 5-minute data points and calculated the closure phase of all triangles. For more information on closure phase residuals, see Section 4.3.3. The EHTC ring images for comparison were generated using the EHTC-DIFMAP pipeline. It should be emphasized that the final image clearly contains an unrealistic CLEAN component, but still shows the same level of closure phase residuals as the image of the EHTC ring.

We also used this final image model to attempt self-calibration of the amplitude and phase for the solution, performed CLEAN, and obtained new images. However, the residuals of the closure phase in the new image were not improved. Therefore, we terminated the hybrid mapping without amplitude calibration. The images shown in Figure 4 are our best final images. The two upper panels in Figures 7 and the panels in Figure 10 are partial extracts from the final images. Figure 8 shows the full image of the last two days of data.

Here we show solutions of the self-calibration for both amplitude and phase using one of the final images, although we did not apply it to the data calibration. The reason we dare to show the unused solutions here is that we believe this study provides insight into the quality of EHT public data and the reliability of our images. We performed our self-calibration using CALIB in the AIPS task in ”A&P” mode with the image (CLEAN components). Other parameters of CALIB are the same as in the first self-calibration (Table 1). Figures 5 and 6 show the self-calibration solutions (the total amount of solutions to be applied for data calibration).

Figure 5 shows the phase solution. Compared to the initial phase solution shown in Figure 2, there is no significant change in the overall perspective. This can be attributed to the fact that the brightness distribution of the observed source is concentrated in the center and does not deviate significantly from the self-calibration solution assuming a single point source. However, there are offsets in the phase changes of JC, LM, PV, and SM stations. In addition, the phases of the two stations in Hawaii, JC and SM, show a larger phase dispersion than the solution of self-calibration assuming a point source on the 100 and 101 days’ data. Although small in comparison, the phases of JC and SM on 95 and 96 days’ data also show phase scatters on the same hours.

The amplitude solution is shown in Figure 6. In general, errors in amplitude are due to noise in the atmosphere and in the receiving system. In addition, changes in aperture efficiency depending on the elevation angle of antenna often cause systematic errors in amplitude. These effects can be measured by an auxiliary method.

For large aperture antennas, gain loss due to offset tracking of the target source from the narrow main beam angle may occur, which is difficult to calibrate.

Furthermore, coherency (phase stability) loss is observed due to the variations in each station clock and atmospheric variations, which is more difficult to measure correctly than other error factors.

Amplitude solutions for AA, AP, PV and SM stations are within $50\%$ fluctuation ($\sim 1.0\pm 0.5$). Such values are often found in amplitude solutions of most of the self-calibrations of VLBI data. On the other hand, JC and LM stations occasionally show large amplitude solutions reaching 10 and 30, respectively. The JCMT station shows a large amplitude value at $T~{}\sim 4.25~{}$hour on the last observation day (101 day). On the other hand, for the LMT station, the amplitude is large for several times as follows.
 (a) $T\sim 1~{}hour$ and $T\sim 2.6~{}hour$ on the first observing day (095 day)
 (b) $T\sim 1~{}hour$ on the second observing day (096 day)
 (c) $T\sim 2.25~{}hour$ and $T\sim 6.25~{}hour$ on the 3rd observing day (100 day).
The EHTC did not show the overall calibrations to be applied, but noticed the sudden large-amplitude errors at the LMT station (Figure 21 in Event Horizon Telescope collaboration (2019d)).

These large amplitude solutions may have implied that the resultant image is significantly wrong. For comparisons, we examined the solutions of self-calibration in the case of the EHTC ring image. Consequently, we found that the self-calibration solutions by the EHTC ring image also demonstrate large amplitudes occasionally, similar to those of our image (Section  5.6). Therefore, if such large amplitudes found in self-calibration solutions are negative signs against the resultant image quality, the results obtained by both the EHTC and our work should be rejected. The EHTC is concerned and considering the fact that some stations require large-amplitude corrections during data analysis. EHTC then analyzed the data from 3C 279, which was observed with M 87, and obtained consistent imaging results from all imaging methods. At the same time, the EHTC found that the amplitude correction was consistent with that obtained with M 87 (Event Horizon Telescope collaboration (2019d)). The amplitude corrections they found are also consistent with those we showed above. In other words, it is natural to consider that such a large amount of amplitude variation actually occurred. To add, the fact that the EHTC obtained a non-ring structure from the 3C 279 data and that the amount of error corrections the EHTC obtained at that time were consistent with those obtained from the M 87 data does not mean that the ring image of the EHTC is the correct image of M 87. The large amplitude solutions from the self-calibrations indicate that the ”calibrated” data released by the EHTC are not of high quality with respect to the amplitude.

In the central core region, we could not find the ring structure reported by the EHTC, but found a core-knot structure. Figure 7 shows the images of the central region ($300~{}\mu\rm as$ square). As noted in Section 3.4, since the data calibration is not yet complete, our final images show the sidelobe structure around actual features. This is a common phenomenon in synthesis imaging with radio interferometers with only a few stations. The images in Figure 7 show that “the unresolved VLBI core” in M  87 has finally resolved into substructures.

The high spatial resolution of the EHT array clearly shows the presence of two bright peaks, i.e., the core and knot structure. The core is indicated by ”C” and the knot by ”K” in the upper left panel of Figure 7. In addition, we found a feature, ”W”, located west about $83~{}\mu\rm as$  apart from the core C. The flux densities from the obtained CLEAN components are $F_{C}\sim 60~{}mJy$ for the core (C), $F_{K}\sim 40~{}mJy$ for the knot (K), and $F_{W}\sim 25~{}mJy$ for the west feature (W). In this observation, the solid angle of features was not so clear. Here, we assume that the solid angle of the emission is $15~{}\mu\rm as$ in diameter, and calculate the brightness temperatures (lower limit). The average brightness of feature C is $T_{\rm b}=1.1~{}\times 10^{10}$ K. Feature K has a brightness of $T_{\rm b}=7.1~{}\times 10^{9}$ K.  Feature W is $T_{\rm b}=4.7~{}\times 10^{9}$ K. Thus, we have detected central features with brightness temperatures higher than the EHTC ring (up to $T_{\rm b}\sim 6\times 10^{9}$ K). The solid angle assumed here is the maximum size of a single, smoothed object that the EHT array can detect in the 230 GHz (Appendix A). Therefore, the actual brightness temperature is likely to be much higher. If the solid angle of the emission is $5~{}\mu\rm as$ in diameter, the brightness temperature of core C reaches $10^{11}~{}K$. If this is the case, the brightness temperature is an order of magnitude higher than the previous measurement cases (Kim et al., 2018; Hada et al., 2016; Akiyama et al., 2015). This is mainly because the size of the emitting region has been identified as smaller due to the higher spatial resolution.

High brightness temperatures were often detected from some active galactic nuclei (Horiuchi et al., 2004; Homan et al., 2011), and can be explained by the Doppler boosting effect of relativistic motion of jet approaching toward us. Previous observations found no high velocity movement in the M 87 central core. Therefore, the brightness temperatures are not due to such Doppler boosting effects. If they actually reflect the physical temperatures, they can be explained easily by the simple RIAF (Radiatively Inefficient Accretion Flow) disks (Kato, Fukue, & Mineshige , 2008; Nakamura et al., 1997). Our observational results are consistent with those of previous studies, supporting the existence of the RIAF disk in the M 87 core (Di Matteo et al., 2003).

There is a clear difference between the two images observed over the five days. According to Event Horizon Telescope collaboration (2019c), they found a change in the closure phase between data sets from the first two days and the last two days. In other words, there was a clear time variation. However, the EHTC could not clearly identify from the structure of that ring where that change occurred (Event Horizon Telescope collaboration , 2019d). We identified the change in the closure phase as due to a change in the core knot structure (features C and K). In particular, the change in the position of feature K was also seen in the trial images during the hybrid mapping process.

Assuming that features are single components, we fitted a Gaussian brightness distribution to each feature and measured the central position and displacement over five days. Relative to the position of feature C, the change in position of feature K is $\Delta\alpha=-0.4~{}\mu\rm as$, $\Delta\delta=-4.5~{}\mu\rm as$ in 5 days, and the proper motion is $0.33~{}\rm mas/yr~{}(v\sim 0.1~{}c)$.  Feature K appears to be approaching feature C as if showing an inflow motion. However, if we look at the differences in the brightness distributions shown in the lower left of the Figure 7, we can see that the changes in the brightness distribution of feature K occur in three places, all at the north end of feature K. In the latter measurement of the position of feature K, feature K appears to be moving south because the brightness distribution of feature C affects the measurement; K is moving north on the line of $PA=-38^{\circ}$ as a whole. The position of feature C has hardly changed, but the location of ”a” where the intensity increased is at the northwest side of feature C. In other words, the structure of features C and K and their time variation can be interpreted as an outflow emanating in the direction of $PA=-38^{\circ}$ from the origin. There has never been a measurement of the motion of a knot so close to the central core.

In comparison, it is difficult to interpret what feature W is. Three hypotheses are presented below.

1. 1.

   Gravitational lensing image. Feature W is morphologically similar to feature K, and the pattern of brightness variation is also similar. This can be attributed to the formation of the gravitational lens image of feature K due to the strong gravitational field of the SMBH in M 87. Assuming that the position of SMBH is approximately equal to the position of feature C (the distance between the core of M 87 and SMBH is $41~{}\mu\rm as$ or 6 $R_{\rm S}$; Hada et al. (2011)), feature W would be located (or at least projected) at 12 $R_{\rm S}$ from the SMBH. There is a possibility that the radio waves emitting from feature K to the far side, orbiting in the strong gravity field of the SMBH, and being changed the propagation direction, come towards us (Black Hole Echo; Saida (2017)). If feature W is such a lensing image caused by the strong gravity field of the SMBH, it should be the image of the backside of feature K, so it is most likely a mirror image of feature K. However, the shape of feature W does not look like such an image. Needless to say, there are many possibilities for a gravitational lensing due to a strong gravity field, so detailed calculations are required to deny it completely; however, the possibility that feature W is a gravitational lensing image is not very high.

2. 2.

   Another central black hole. Feature C is the primary SMBH of M 87, and feature W is a secondary SMBH orbiting the primary SMBH. If there is a binary SMBH in M 87, it can be permanently observed with the EHT array. Based on these observations, we calculated two possible orbits.

   1. (a)

      The proper motion of feature W ($\mu=0.34~{}\rm mas/yr,~{}\rm v\sim 0.09~{}c$) is assumed to be due to a circular orbit motion, and its orbital radius is assumed to be the separation $R=83~{}\mu\rm as$ from feature C. In this case, the orbital period T is $\sim 1.5$ years. If the real radius is $R=1.4\times~{}10^{3}~{}au$, the mass of central object $M_{c}$ is only $1.2\times~{}10^{6}~{}M_{\odot}$.  Since the estimated mass is too small as compared to those of the previous M 87 studies, this assumption must be rejected.

   2. (b)

      We assume that the observed proper motions and structure change of feature W are only due to changes in surrounding matters, and that the measured proper motions of feature W have nothing to do with its orbital motion. In other words, we assume here that no change in the position of the center of gravity of feature W is observed. Also, it is assumed that feature C has an SMBH with a mass of $6\times 10^{9}~{}M_{\odot}$ and that the orbital radius of feature W is the distance between features C and W. The distance between them is $83~{}\mu\rm as$ ($1.4\times 10^{3}$ au or $11.9~{}R_{\rm S}$).  It is consistent with the 86 GHz core size of $\sim 80~{}\mu\rm as$ at 86 GHz observed in 2014 (Hada et al., 2016), suggesting that the two features C and W are not transient. Also sinusoidal oscillations of the position angle of the jet were observed with a period of roughly 8 to 10 years (Walker et al., 2018). If the two features C, and W compose a binary of black holes, its orbital motion can cause such jet oscillations. Certainly, the apparent separation of approximately $1.4\times 10^{3}$  au is too short to explain the observed period of the jet oscillation. However, if the real distance is longer by a factor of about 3.42, which is the correction factor of the viewing angle of the jet axis from us ($\sim 17~{}^{\circ}$), the orbital period of the binary can be $\sim 10$ years.

3. 3.

   Unstable initial knot. Feature W is another knot moving toward a different direction. The jet of M 87 is known to have a wide opening angle at scales well below $1\rm~{}mas$ (Junor, Biretta & Livio, 1999). Furthermore, Walker et al. (2018) found evidence from 43 GHz observations that the initial opening angle $\theta_{app}$ is $\sim 70~{}^{\circ}$. We found the angle $\angle KCW$ is $70~{}^{\circ}$, and further that the line of the average jet axis ($PA=-72^{\circ}$) divides its angle almost evenly into $34^{\circ},~{}36^{\circ}$. Furthermore, the lines CK and CW extend in the directions of $PA=-38^{\circ}$ and $PA=252^{\circ}$, respectively. These directions are very similar to the ridges observed at 43 GHz from where Walker et al. (2018) measured the initial opening angle. We guess that not only feature K, feature W is also an initial knot that has just emerged from the core, and still, the shape is very unstable and shows changes significantly.

Adopting the most conservative hypothesis, feature W, like feature K, can be understood to be a knot that represents the initial jet structure.

As we will discuss in Section 4.3.2, the core-knot structure (features C and K) is robust in the sense that it can be obtained with different imaging parameters. On the other hand, the  $40~{}\mu\rm as$ ring of the EHTC is sensitive to BOX parameters and can be easily destroyed, even if it can be created as shown in Section 5.7. Due to the robustness of the core-knot structure, we consider it to be a real structure. On the other hand, the feature W is sensitive to the BOX size, so its detection is not as reliable as it could be. However, the structure corresponding to feature W had already appeared in the first imaging results (Section 3.2.2). That is, feature c in Table 3 is in a similar position to feature W and also shows the largest asymmetry. Also, if we run the EHTC-DIFMAP pipeline without its BOX setting, an emission feature appears at the position close to feature W (lower right panel of Figure 24). These suggest that the feature W is also a real structure.

Here, we show the overall brightness distribution (Section 4.2.1) and that of so-called the jet launching region, which is a few mas away from the core (Section 4.2.2).

As mentioned in Section  3, our calibration method was limited to self-calibration because the public the EHTC data do not contain raw data.

We also had to give up on the amplitude self-calibration because the closure phase residuals were not reduced as compared to the case when only phase calibration was performed. Therefore, the calibration is not yet fully complete. As clues to the reliability, we describe the properties of the final images from three aspects: detection limit (Section 4.3.1), robustness (Section 4.3.2), and the self-consistency of our imaging (Section 4.3.3), where our images show better self-consistency than those of the EHTC.

We here investigate the u-v sampling of the EHT array for M 87 observations. The u-v coverage itself is shown in the EHTC paper I (Event Horizon Telescope collaboration , 2019a). It shows a good u-v coverage as an interferometer of practically 5 stations and 10 baselines.  We looked at the number of data samples from another point of view. Figure 18 shows the distribution of spatial Fourier components (fringe spacings) sampled by EHT for M 87 during the four observational days. The number of sampled data is plotted against the fringe spacings in $\mu\rm as$ unit. We can see that there are no sample for ranges $d=26-28~{}\mu\rm as$, $d=44-46~{}\mu\rm as$, and $d=95-100~{}\mu\rm as$ ¹¹¹¹11 Fringe spacings larger than $100~{}\mu\rm as$ are omitted in the histogram. We are interested in the samplings contributing to very high spatial resolutions less than $100~{}\mu\rm as$. . The number of sampled data for the size of the EHT ring ($d=42\pm~{}3~{}\mu\rm as$) is quite small, and for spacing between $d=44-46~{}\mu\rm as$, there are no sampled data.

The few numbers of sampling around $42~{}\mu\rm as$ fringe spacings should affect the imaging performance. In the following subsections, we investigate how this lack of spatial Fourier components can affect the imaging result.

The EHTC show no description of the dirty beam structure. We, then, examine the dirty beam calculated from the u-v coverage. Dirty beam is the term used in the field of radio interferometer and its meaning is the same as that of PSF (Point Spread Function) in optics. In other words, it is the diffraction image of a single point when the data calibration is perfect. If we can obtain all spatial Fourier components, the dirty beam becomes a two-dimensional $\delta$-function.

In practice, we can obtain only a limited sample of spatial Fourier components, and thus the dirty beam has a complex shape. Figure 19 shows the dirty beam of the EHT array on the first day observation of M 87. The FWHM of the main beam is approximately 20$~{}\mu\rm as$. We can see a point-symmetric pattern around the main beam. This pattern includes peaks lower than the main beam. Such peaks are often called “sidelobes” in radio astronomy. We can see that the distances between the main peak and the first peaks of sidelobes are about 45$~{}\mu\rm as$. ¹²¹²12 Further, the first sidelobe levels reach more than 70 % of the height of the main beam. The sidelobe level is still extremely high even when changing u-v weights. Natural and uniform weighting cases, both show that the first sidelobe levels reach more than 60 % of the height of the main beam (see Appendix B). The separation between them corresponds to the radial range for which the spatial Fourier components are missing, and is very close to the diameter of the EHTC ring ($42~{}\pm~{}3~{}\mu\rm as$). Thus, it is important to clarify whether or not this PSF structure has affected the ring image or not. In addition, Appendix B shows other PSF shapes by different weightings of u-v points. Even in these PSF shapes, the separations between the main beam and adjacent sidelobes do not change much. They also show substructures with $\sim 40~{}\mu\rm as$ spacings.

Here, we show what happens if we are not very careful in using the dirty map when applying the hybrid mapping process. A dirty map is an image created by simply performing an inverse Fourier transform on the obtained spatial Fourier components. Therefore, it is influenced by the data sampling bias, that is, the structure of PSF is convoluted into the dirty map. Furthermore, if the data calibration is insufficient, it will be reflected in the dirty map, and the obtained structure in such a dirty map can be far from the actual image. Therefore, it is dangerous to perform self-calibration using the dirty map as the image model.

It is certain that in the case of multi-element radio interferometers like ALMA, or VLA, the dirty beams have sharp main beams and very low sidelobes. Only in such cases, it is not so dangerous to estimate the true image by the dirty map. While, in the case of VLBI observations, the dirty beam is comparatively dirty, usually, we do not estimate the true image from them and do not use them as the model image for self-calibration.

We, however, try that here and show what happens. The left-side panels of Figure 20 show the dirty maps from the data of the first day observations by EHT. To obtain these maps, we applied the phase-only calibration by the self-calibration technique using a point source as the image model. We can see a ring-like structure at the center, and many images which look like the ghosts of this central ring-like structure.

This is not a blurred intrinsic image of M  87, but a strong reflection of the substructure in the dirty beam (PSF) of the EHT array. Figure 21 shows the dirty map and the dirty beam. Not surprisingly, they agree rather well, and that means the central ring-like structure we showed in Figure 21 is just the reflection of the shape of the dirty beam and not a physical reality. If we were to believe there should be a ring, it would be very natural to select the partial image from what looks like a ring in the dirty map. If we do so, what will happen? To answer this question, we tried the calibration of phase and amplitude using this ring-like structure at the center as the image model. We calculated the amount of calibration of the amplitude and phase by the self-calibration method. After the calibration is applied, we made an image using the CLEAN algorithm assuming that the source is single and compact. Here, the CLEAN subtraction area is limited by a circle with $30~{}\mu\rm as$ radii centered at ($+2~{}\mu\rm as$, $+22~{}\mu\rm as$) from the map center, which is the BOX setting as used in the EHTC-DIFMAP pipeline (Section 5.5).

The obtained ring image is shown in the right panel of Figure 20. The size of the ring is close to that of the EHTC ring. What we also would like to emphasize here is that in order to create a ring structure, we need a narrow BOX setup. Without a narrow BOX, a ring structure cannot be created. If we change the position or size of the BOX, the ring deforms. When we use BOX offset setting of $22~{}\mu\rm as$ from the center, we get a ring very similar to the EHTC ring. In other words, the shape strongly depends on the location of the imaging region. As we will show in Section  5.4, the EHTC BOX setting limits the FOV so narrowly that it can produce the EHTC $\sim 40~{}\mu\rm as$ ring shape even from simulated data that does not contain the ring shape.

Close inspection reveals that the ring-like structure of the dirty map is at an offset position from the center. Since the structure of the PSF always shows center symmetry, we can suppose that the offset ring-like structure is not due to the PSF, but really exists at the offset position. However, the shape of Figure 3 shows that this is not the case. Figure 3 shows the structure of the dirty beam shape after deconvolution from only the brightest point at the center of the dirty map. There is no ring-like structure here. We can see that the ring-like structure of the dirty map (Figure 20) is created by convolution of the dirty beam to the brightest point in the center. The ring-like structure is offset because the true source image inherent in the data does indeed exhibit non-center symmetry, but it is not an offset ring, but a core-knot structure as shown in Figure 7.

In this section, we demonstrate that, with the u-v coverage of EHT for M 87, one could “observe” a ring even if the physical sources have different structures.  As the physical source structure, we consider a single point with flux density of 1 Jy at the center. We made the simulation visibility data by applying the Fourier transformation to the single point image (We used UVMOD, a task in AIPS.). The u-v coverage is exactly the same as that of EHT for M 87 observations. As for noise, we consider two cases, one with relatively large noise (case 1) and the other with no noise (case 2). For case 1, we added thermal noise. The noise level is proportional to the weight of each data noted in original FITS public data. The signal-to-noise ratio (S/N) is 0.01 on average. We then perform experiments of calibration and imaging. For calibration, we obtain solutions from self-calibration with incorrect image models that are different from the true structures of the source, then apply the solutions to the data. We inspect what kinds of images appear from CLEAN imaging. For the self-calibration, we assume models of two incorrect images. One is the EHTC ring image (model A) and the other is a pair of two points separated by $45~{}\mu\rm as$ (model B). Model B is a pair of points of 0.5 Jy located at ($0~{}\mu\rm as$, $0~{}\mu\rm as$) and ($0~{}\mu\rm as$, $+45~{}\mu\rm as$) respectively. The locations of these two points roughly correspond to the positions of the main peak and the first sidelobe peak in dirty beam. We get calibration solutions of amplitude and phase from the self-calibrations (We used CALIB, a task in AIPS.).

By using IMAGR, a task in AIPS, we performed the CLEAN imagings. We limit the CLEAN subtraction area by a narrow BOX setting; a circle with $30~{}\mu\rm as$ radii centered at ($+2~{}\mu\rm as$, $+22~{}\mu\rm as$) from the map center, which is a mimic BOX setting as used in the EHTC-DIFMAP pipeline (Section 5.5). The left panel in Figure 22 shows the ring image obtained from the large noise data (case 1). The obtained image from CLEAN is identical to the model image A used in self-calibration. When the S/N is very low, the CLEAN image after self-calibration can be identical to the model image used. Therefore, the result of case 1 does not tell much.

The results of case 2 are surprising. The middle panel in Figure 22 shows the ring image obtained from no noise data. The obtained image from CLEAN is very close to the model image A used in self-calibration. Even when there is no noise, if we start with the assumption that there is a ring in the self-calibration image model, with the u-v coverage of EHT, we obtain a ring similar to that obtained by the EHTC.

The right panel in Figure 22 shows the result of using model image B for self-calibration. Though we used the data of a single point with no noise, we obtained a two-point image, which is very close to model B used in self-calibration. In addition, the appeared two points show small shifts from the positions in model B. The center point shifted towards the east, while the upper point shifted towards the west. These shifts seem to be the influence of the structure of the dirty beam. The shift of the center point is along the ridge elongated from the main beam in the dirty beam, while the shift of the upper point is towards the peak position of the first sidelobe.

The results of these experiments show that the reproduced images are very close to the image model used in the self-calibration. Such a result does not occur when there is sufficient u-v coverage, as in the case of VLA and ALMA. Unfortunately, the u-v coverage of EHT is not sufficient, and it has a special bias that makes it prone to producing shapes (including rings) with a size of $\sim 40~{}\mu\rm as$. We have shown that the shape of the PSF (essentially the effect of the sampling bias that makes it easy to create $\sim 40~{}\mu\rm as$  structures) has the power to create $\sim 40~{}\mu\rm as$ ring structures even from data with only one central point structure (centrosymmetry). However, the most powerful factor is not data correction by self-calibration with wrong model images. Here and in Section 5.3, we showed that the very narrow FOV setting of the BOX has great power to create a $\sim 40~{}\mu\rm as$ ring structure. In Section  5.5, we will discuss the setting of BOX in the EHT-DIFMAP pipeline.

In this section, we investigate the data processing pipeline that the EHTC used for imaging the EHTC ring with DIFMAP ¹³¹³13https://github.com/eventhorizontelescope/2019-D01-02. Among the three methods the EHTC used for the imaging, DIFMAP is the closest to the usual procedure ¹⁴¹⁴14difmap/EHT_Difmap in the above web page. Other two are rather new, and difficult to study here. Since all three methods gave essentially the same image, we believe it is sufficient to analyze one of the three methods to find out why the ring structure was formed. The EHTC-DIFMAP imaging team used the hybrid mapping technique to calibrate the data and obtained the image. The hybrid mapping method is the standard method of VLBI calibration and imaging.

Let us first summarize the standard procedure of VLBI data calibration and imaging. The radio interferometer samples the spatial Fourier components of the brightness distribution of the observed source (often called visibility). Theoretically, the brightness distribution of the observed source can be obtained by collecting the samples and performing the inverse Fourier transform. However, there are actually two problems as follows. (a) The sampled visibility data contain errors, so they are not the correct spatial Fourier components. Therefore, it is necessary to calibrate them. Self-calibration is a method of obtaining calibration solutions by using an assumed image. If the assumed image is not correct, the correct calibration solution cannot be obtained. (b) The number of observed samples is limited. The inverse Fourier transform alone does not produce the correct image. The PSF does not become a point. Therefore, the PSF shape must be deconvolved. For this purpose, image processing is performed. As a method to obtain correct calibration solutions and to attain clear imaging results, hybrid mapping, which alternates between the above two tasks, is usually used in VLBI imaging. We use self-calibration for calculating the calibration solution using a tentative image model, and CLEAN for performing deconvolution of the scattered PSF shape.

Hybrid mapping is the following algorithm. (a) Assume an image model for self-calibration. For the first one, a 1-point model is usually used. (b) Calculate tentative calibration solutions from self-calibration using the image model. (c) Calibrate visibility data using the tentative calibration solutions. (d) Obtain the next tentative image using CLEAN from the calibrated visibility data. CLEAN image is composed of point sources (CLEAN components). (e) Make the next image model for self-calibration by picking up reliable CLEAN components from the image. (f) Go back to (b) and repeat the iteration from (b) to (e) until a satisfactory image is obtained.

The visibility (spatial Fourier component) is a complex quantity and so it has amplitude and phase. It is safe first to repeat self-calibrations only for phase solutions until a nearly satisfactory image is obtained and to perform self-calibration for amplitude and phase by using that image. This is because self-calibration for amplitude tends to give a solution which is too close to the model image used in self-calibration (Cornwell & Fomalont, 1999).

However, the repetition of amplitude self-calibration is often performed under the careful check of the practical quality of the data, and then, in many cases, good calibration solutions can be obtained.

In addition, the BOX technique is an auxiliary means used for imaging with CLEAN. (It is essentially unrelated to the hybrid mapping method.) We can limit the area of CLEAN subtraction intentionally by the BOX setting. The technique may give us a good image, despite an incomplete calibration. As already mentioned, we used the BOX technique in our analysis. Also the EHTC-DIFMAP team did.

Now we investigate the EHTC-DIFMAP procedure. The EHTC-DIFMAP imaging team performed the self-calibration only for phase solutions with a single point source model, following the general procedure of the hybrid mapping method. However, after the first phase self-calibration, they used the BOX technique with a very narrow BOX¹⁵¹⁵15BOXCircMask_r30_x-0.002_y0.022.win in the web in the first imaging trial. The CLEAN subtraction area was restricted to within the circle of $60~{}\mu\rm as$ diameter specified by the BOX. This means that from the beginning of their data analysis, the EHTC-DIFMAP team assumed that the image is a very compact single one, and excluded other possibilities.

The BOX is roughly a circle of diameter $60~{}\mu\rm as$ composed of 30 small rectangles. By comparing the BOX shape and the dirty beam we found that the BOX covers the main beam and the first sidelobe, but leaves out the second and other sidelobes. The BOX is not located at the phase center but offset by $+22~{}\mu\rm as$ in the y-direction ($\delta$ direction). This offset of $22~{}\mu\rm as$ coincides with the radius of the EHTC ring. Figure 23 shows the EHTC BOX position on the dirty beam (PSF) aligned by their phase centers. The area of the EHTC BOX covers (1) the main beam peak, (2) the ridge extending to the left (east) from the main beam peak, (3) a part of the first sidelobe located north of the main beam, and (4) an area of negative strength near the center of the box area. The structure of the dirty beam within the EHTC BOX is close to the shape of the EHTC $\sim 40~{}\mu\rm as$ ring.

The PSF structure within the BOX explains not only the diameter of the EHTC ring, which is $\sim 40~{}\mu\rm as$, but also the asymmetric structure of the EHTC ring. The EHTC ring has an asymmetric structure with a bright south side. The brighter south side corresponds to the main beam (the brightest in the box), while the darker north side of the ring corresponds to the first sidelobe in the north (the less bright peak in the box).

If such a narrow BOX with the offset is used from the beginning to the end of the hybrid mapping process, the effect of the dirty beam can be enhanced as we showed in Section 5.3. It is quite clear that the EHTC $\sim 40~{}\mu\rm as$ ring is the result of such enhancement of $\sim 40~{}\mu\rm as$ substructures in PSF.

We actually ran the EHTC-DIFMAP pipeline to investigate its behavior and performance. The left panel of Figure 24 shows the resultant ring image from running the EHTC-DIFMAP pipeline with the default parameter settings. On the other hand, the right panel of Figure 24 shows the image obtained by removing the BOX setting (Figure 23) from the EHTC-DIFMAP pipeline. The ring in the BOX setting region has been destroyed and the brightness distribution has been moved to outside the BOX setting region. The new brightness distribution shows a core-knot structure in the center, similar to the results of our analysis. In addition, several emission peaks appear in the $PA=45^{\circ}$ and $PA=225^{\circ}$ directions from the center. These are probably sidelobes, but it is very interesting that one emission peak in the $PA=225^{\circ}$ direction corresponds to the feature W we found in our analysis. This result implies that the very narrow BOX setting of the EHTC-DIFMAP pipeline is the main factor that creates the EHTC ring image.

We also created simulated data and applied the EHTC-DIFMAP pipeline to it to investigate what kind of images are obtained. The purpose of this is to check if the EHTC-DIFMAP pipeline always calibrates the data correctly and reproduces the true picture of the data. We created simulated data with the task UVMOD in AIPS and acquired images with both the task IMAGR in AIPS and the EHTC-DIFMAP pipeline for comparison. Figure 25 shows examples of the simulation results. Three cases are shown. (1) a core-knot model aligned in about $PA=-38.7^{\circ}$ direction (left panels), (2) a two-point model aligned North-South direction (center panels), and (3) a $40~{}\mu\rm as$ diameter ring, but the ring width is $2~{}\mu\rm as$, different from that of the EHTC ring (right panels).

The top panels show images obtained from data simulated using the AIPS IMAGR task with the same box settings as the EHTC-DIFMAP pipeline. The IMAGR produced the expected image.

The bottom panels of Figure 25 show the results of the EHTC-DIFMAP pipeline. The original images of the simulated data could not be reproduced. In the case of the core-knot model (lower left panel), the knot component is lost. Also, in the case of the two-point model (bottom center), the weak north point has disappeared. Even more surprisingly, in the case of the $40~{}\mu\rm as$ diameter ring model (bottom right panel), the image quality is degraded and the ring structure is obscured.

Since these images are compact and fit within the BOX setting of the EHTC-DIFMAP pipeline, this phenomenon is not due to a narrow field of view caused by the BOX setting. This is probably because the performance of the hybrid mapping process of the EHTC-DIFMAP pipeline strongly depends on the image structure and noise structure of the input data.

Note that our simulation this time is not only that shown in Figure 25. The summary of the simulation is as follows.

1. 1.

   A single point source model of 1 Jy with no noise. The image was reproduced from both IMAGR in AIPS and the EHTC-DIFMAP pipeline.

2. 2.

   Complete noise model. No specific image was detected from the IMAGR in AIPS, but a single point was obtained from the EHTC-DIFMAP pipeline. This is probably due to the fact that the first self-calibration was performed using the image model of a single point source. This is a common phenomenon that occurs when self-calibration is applied to data sets with very low S/N.

3. 3.

   Core knot model. Two points model consisting of a point of 1 Jy at the center ($0~{}\mu\rm as$, $0~{}\mu\rm as$) and a point of 0.25 Jy at ($-25~{}\mu\rm as$, $+20~{}\mu\rm as$). The distance between them is about $32~{}\mu\rm as$ s and they are located along a line in the direction $PA=-38.7^{\circ}$. 12 simulation data with noise ranging from zero to 640 Jy/weight were generated. IMAGR in AIPS detected the core-knot structure as expected. In the EHTC-DIFMAP pipeline, the knot component disappeared in 11 of the 12 cases.

4. 4.

   Ring images with diameters around $40~{}\mu\rm as$. The ring width is $2~{}\mu\rm as$, and the ring diameters are 30, 35, 40, and 45 $\mu\rm as$. The flux density of the image is 1 Jy. No noise was added. Using IMAGR in AIPS, the hole in the center of the ring is obscured when the ring diameter is smaller than $30~{}\mu\rm as$. When the ring diameter is larger, the ring can be recognized. On the other hand, using the EHTC-DIFMAP pipeline, the structure of the ring became ambiguous in all images, like the example shown in Figure 25.

5. 5.

   Two-point model (shift in the north-south direction). Place a point of 1Jy at the center ($0~{}\mu\rm as$, $0~{}\mu\rm as$) and a point of 0.5 Jy north of it. There are 11 intervals: 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, and 50 $~{}\mu\rm as$. Using IMAGR in AIPS, two points cannot be separated and resolved when the spacing is very narrow. At a separation of 10$~{}\mu\rm as$, the image is elongated in the north-south direction. If the separation is larger than 25$~{}\mu\rm as$, the two points are reproduced separately. Using the EHTC-DIFMAP pipeline, when the separation is larger than 45$~{}\mu\rm as$, the two points are reproduced correctly, but otherwise the north point is not reproduced in the correct position.

6. 6.

   Two-point model (east-west shift). A point of 1Jy is placed at the center ($0~{}\mu\rm as$, $0~{}\mu\rm as$) and a point 0.5-Jy is placed at $\delta=42~{}\mu\rm as$. The R. A. position of the second point is shifted every $5~{}\mu\rm as$ from -20$~{}\mu\rm as$ to 20$~{}\mu\rm as$ to create 9 simulation data. The positions of these two points are reproduced from IMAGR in AIPS. the EHTC-DIFMAP pipeline reproduces these two points in four cases. In the other four cases, the position of the north point deviated from the correct position by more than 10 $\mu\rm as$.

Thus, in most cases, the EHTC-DIFMAP pipeline did not reproduce the intrinsic image of the simulation data. Figure 10 of Event Horizon Telescope collaboration (2019d) shows that their procedure successfully reproduces images from simulated data (ring, crescent, disk, and double source), but their results are inconsistent with ours (at least for the EHT DIFMAP case). In our simulations, noise is either not added at all, or if added, it is thermal noise. IMAGR in AIPS almost reproduces the input model images, while the EHTC-DIFMAP pipeline does not. The difference between the two methods is the hybrid mapping process. The EHTC-DIFMAP pipeline lacks the performance of the general calibration of the data described in Section  5.8.3.

The EHTC papers give the detailed description of the data calibrations at the pre-stage process, but not much about the calibrations performed during the imaging process (through the self-calibration in the hybrid mapping process of the EHTC-DIFMAP imaging team case). This is an insufficient description of data calibrations because it is rare that only the pre-calibrations are sufficient for obtaining VLBI fine images. We, therefore, estimate and show the amounts of calibrations performed by the EHTC during their imaging process. We attempted to reconstruct them in the following procedure. First, following the EHTC open procedure, we reproduced the EHTC ring image of the first observing day (the left panel in Figure 26). Then, using the CLEAN components of the ring image as the model image, we performed self-calibration by the task CALIB in AIPS and got solutions for both phase and amplitude. The parameters of CALIB are shown in Table 6. In Figures 27 and 28, we show the total amount of calibrations both of phase and amplitude in the case of the EHTC ring image model. Note that we performed CLEAN imaging to verify if the solutions from the self-calibrations are consistent with the original the EHTC ring image. First, we flagged out data points where the amplitude solutions were more than 4 by using the task SNCOR in AIPS. Here we followed the standard VLBI teaching, which states that ”discarding bad data will make a better image than keeping them”. Then, we applied the solutions to the selected data and obtained an image using the task IMAGR in AIPS for CLEAN. The right panel of Figure 26 shows the resultant image whose quality is similar to that of the original EHTC ring. Thus, “our” solutions of self-calibrations successfully reproduced the EHTC ring image.

The phase and amplitude solutions for the EHTC ring image are very similar to, or worse than our results in Figures 5 and 6. Therefore, if such large amplitudes and their rapid variations found in the self-calibration solutions are negative signs against the resultant image quality, both the results of our final images and the EHTC ring images should be rejected. This implies that the ”calibrated” data released by the EHTC are not of high quality (as is often the case with VLBI archival data).

In Section 5.6, we reproduced the ring structure obtained by the EHTC, by using the very narrow box used by the EHTC. If the ring image is the definitely correct image, it should only weakly depend on the BOX size. In fact, we studied the effect of changing the BOX in Section 4.3.2, and found that the main structure of our final image (features C & K) does not disappear. In this section, we investigate how sensitive the ring structure obtained by the EHTC to the BOX method they used.

We calibrated the data by applying the calibration solutions obtained from the self-calibration using the EHTC ring as the image model, which are shown in Figures 27 and 28. We performed CLEANs using BOXes of circles with diameters of (a) $60~{}\mu\rm as$, (b) $80~{}\mu\rm as$, (c) $100~{}\mu\rm as$, (d) $120~{}\mu\rm as$, (e) $240~{}\mu\rm as$, and (f) $450~{}\mu\rm as$. Their centers are located at the same position as that of the EHTC-DIFMAP case. Other parameters of the CLEANs are the same as those used for the right image in Figure 26.

The results are shown in Figures 29 and 30. Figure 29 shows views of the entire imaging area of $\rm 2~{}mas$ square, and Figure 30 shows enlarged views of the central $256~{}\mu\rm as$ square. The BOXes used are indicated by blue circles. The ring structure in panel (a) ($D=60~{}\mu\rm as$) is the same as that in the right panel of Figure 26. Here, we can see the ring structure similar to that obtained by the EHTC. The ring extends beyond the BOX area. This is because the obtained CLEAN components are near the boundary of the BOX, and a $20~{}\mu\rm as$ restoring beam is convolved to form the CLEAN map. In the cases of (b) $D=80~{}\mu\rm as$ and (c) $D=100~{}\mu\rm as$, the ring is still recognizable. However, we can see CLEAN components not on the ring but near the boundary of the BOX when the BOX becomes larger. In the case of (d) $D=120~{}\mu\rm as$, there are several CLEAN components almost on the boundary of the BOX, though the ring structure is still recognizable. In the cases of (e) $D=240~{}\mu\rm as$ and (f) $D=450~{}\mu\rm as$, there is no ring. Instead, one elongated bright spot appears. We also tried larger BOX cases, $D=900~{}\mu\rm as$, and $D=1300~{}\mu\rm as$, and the resulted image is similar to the case of $D=450~{}\mu\rm as$. They show a bright spot at the map center.

If the field of view (FOV) is smaller than the case (d), a ring image will appear, and if FOV is larger than that, the ring image will be destroyed. In other words, the boundary FOV is around $120~{}\mu\rm as$. It is a queer coincidence that the EHTC imaging teams set the FOV to be $128~{}\mu\rm as$ (Event Horizon Telescope collaboration , 2019d), just around the boundary FOV. (It is $60~{}\mu\rm as$ in the case of the EHTC-DIFMAP pipeline as already mentioned) ¹⁶¹⁶16Figure 6 in Event Horizon Telescope collaboration (2019d) shows similar test results to ours, but they only show the results with BOX sizes narrower than $100~{}\mu\rm as$ in diameter..

Thus, we can conclude that the EHTC ring appears only when a very narrow BOX is applied.

Here we explain why the EHTC created the image of the ring. The EHTC described the details of their imaging process in Event Horizon Telescope collaboration (2019d), in which we found some of the reasons why they created the artifact. The EHTC conducted various surveys, but their methods were not objective and biased towards their desires from the very beginning of their analysis. They also failed to perform the basic data checking that VLBI experts always do. These points are noted below.

We here summarize the reasons why the EHTC obtained the artifact image unintentionally. The fundamental reason is the sampling bias in EHT’s u-v coverage for M 87, i.e., the missing $\sim 40~{}\mu\rm as$ spatial Fourier components (Section 5.1). The sampling bias tends to create an artificial structure of $\sim 40~{}\mu\rm as$ in imaging. We confirmed that the artificial $\sim 40~{}\mu\rm as$ sized structures appear in the dirty beams (PSF) (Section 5.2, Appendix B) and in our CLEAN imaging results (Appendix C).

The CLEAN algorithm was originally designed to remove false structures such as sidelobes that appear in PSF, but it does not always succeed in removing them. As well, other imaging methods of the EHTC, if no special countermeasures are taken against u-v sampling bias, will produce artificial structures in the imaging results with sizes of $\sim 40~{}\mu\rm as$ of the sampling bias origin. Also, u-v sampling bias effects are independent of the need for data calibration. Therefore, even with imaging methods based on observables that are free from systematic errors such as closure phase, we have no choice but to be aware of the effects. The sampling bias effect does not produce serious errors if careful data analysis is performed. However, the EHTC expected a ring-like structure of $\sim 40~{}\mu\rm as$ in size from the beginning, and thus performed image synthesis with a very narrow FOV ($128~{}\mu\rm as$ in size). A large FOV setting will not enhance the artificial ring image of $\sim 40~{}\mu\rm as$ size. On the other hand, the very small FOV set by the EHTC enhances the sampling bias effect to create the $\sim 40~{}\mu\rm as$ ring structure (Section 5.7). There were large amplitude errors in the observed data possibly originating from atmospheric variations (Section 3.5). This could have played a role in confusing the initial analysis of the EHTC.