The Jet and Resolved Features of the Central Supermassive Black Hole of M87 Observed with the Event Horizon Telescope (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 reliable for the precise calibration of VLBI data (Pearson & Readhead 1984; Readhead & Wilkinson 1978; Schwab 1980). VLBI systems are not so stable in phase and amplitude as compared to those of connected radio interferometers. In addition, millimeter- and submillimeter-wave observations are strongly 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, ¹¹ Greisen 2003).

3.2.1. Solutions of Self-calibration Using a Point-source Model

As a first step in this process, a single point source (located at the origin) was used as the first image model to obtain a solution for the visibility phase calibration from the self-calibration. The parameters used for the task CALIB are listed in Table 1.

Table 1. Parameters of CALIB for the First Self-calibration

Parameters
SOLTYPE	“L1”
SOLMODE	“P” (phase only)
SMODEL	1,0 (1 Jy single point)
REFANT	1 (ALMA)
SOLINT (solution interval)	0.15 (minute)
APARM(1)	1
APARM(7) (S/N cutoff)	3

Download table as:  ASCIITypeset image

As mentioned in Section 2, the coherence time of EHT public data is very short. The solution interval (SOLINT) was set to 0.15 minutes. We set the S/N cutoff = 3 for safety. This S/N cutoff value is larger than what many researchers use in the end. Solutions that did not meet the criteria (S/N cutoff) were flagged and abandoned.

Figure 2 shows the phase solution for the first step. Because phase is a relative quantity, the phase of the ALMA station (AA) is used here as a reference. For all stations, the nonzero and time-varying phase values were calculated by self-calibration. The four stations, APEX (AP), The Submillimeter Telescope at the Arizona Radio Observatory (SMT, hereafter AZ), Large Millimeter Telescope in México (LM), and Pico Veleta (PV), always show the same respective trends over the four days of observations, suggesting that there are errors in station positions (a sinusoidal curve with a period of one sidereal day is observed at all stations when there is an error in the position of the observed object).

Zoom In Zoom Out Reset image size

Another feature is the phase difference that occurs between IF1 and IF2, which is almost fixed for all stations except James Clerk Maxwell telescope (JCMT, hereafter JC) and AA (phase reference station), respectively.

If the EHT public data were sufficiently calibrated, the above two phenomena should not appear. In conclusion, the “calibrated” data published by the EHTC is not yet sufficiently calibrated. In order to obtain reliable images, the EHT’s public data needs to be further calibrated.

3.2.2. The First CLEAN Map

Figure 3 shows the CLEAN (Clark 1980; Högbom 1974) image obtained from the data after applying the first phase calibration solution shown in Figure 2. The parameters of the imaging of IMAGR are shown in Table 2. (In all figures showing the imaging results, the x-axis indicates relative R.A. and the y-axis relative decl.) The purpose of this imaging is to find the second brightest component following the central brightest peak. As will be explained in Section 5.1, despite the fact that the EHTC has utilized a large number of stations, the u−v coverage of the EHT array is formed by only 7 stations, or actually 5 stations if we exclude the very short baselines. The synthesized beam (dirty beam) of the EHT is not so sharp as those of multielement interferometers such as ALMA and Very Large Array (VLA). It is not easy to find the complex brightness distribution of the observed sources from a tentative map composed of such a scattered dirty beam (PSF). (We show the dirty beam and dirty map in Figures 20 and 21.) Therefore, we performed CLEAN, specified by the parameters shown in Table 2. This method is effective when the structure of the observed source is not point symmetric. We set the loop gain (GAIN) to 1.0 and extracted all of the brightest peak from the dirty map in the first CLEAN subtraction. Next, this component was replaced by a sharp Gaussian restoring beam and combined with the brightness distribution of the remaining dirty map. The image in Figure 3 was created in this way.

Zoom In Zoom Out Reset image size

Table 2. Parameters of IMAGR for the First Trial Imaging

Parameters
DOCALIB	2
CELLSIZE	1.22011 × 10−6, 1.22011 × 10−6 (arcsec)
FLDSIZE	8192, 8192 (pix)
ROBUST	0
NITER	1
GAIN	1.0

Download table as:  ASCIITypeset image

This has the effect of removing the bright but scattered PSF shape of the brightest point that dominates the dirty map, and clarifying the presence of the second brightest component in the image. Note that if the data is not properly calibrated and the actual PSF corresponding to a point source differs from the theoretically calculated PSF shape, the brightness distribution caused by the brightest point may remain in the afterimage. However, such remaining brightness distribution also shows a point-symmetric structure with respect to the location of the brightest point (practically the same location as the center of the map), and does not contribute to the asymmetric structure of the image. Therefore, if there is an asymmetric structure in the image, it is not related to the brightest component, but is due to another bright point source. So, by searching for the asymmetric structure, we can find the second component of the observed source.

This image (Figure 3 shows its central 600 μas square) has a nearly point-symmetric structure with respect to the center of the map. The overall feature is a series of multiple ridges in the PA = 55° direction. This structure is due to the nonuniformity of the u−v coverage. In addition to the central P, there are several other bright features. The peak brightness of these features is shown in Table 3. Features a, b, c, d, e, f, g, h, i, j, and k have corresponding features located at their symmetry points (denoted as a*, b*, c*, d*, e*, f*, g*, h*, i*, j*, and k*). Curiously, the features located in the upper right from the center are always brighter than their counterparts located in the lower left, i.e., the brightness ratio is greater than one. This may indicate the existence of a large-scale asymmetric brightness distribution in the observed object, extending from the center to the upper right. This is roughly consistent with the M87 jet propagation direction PA = −72° (Walker et al. 2018). Examining the brightness ratio of each pair, we find that the pair c & c* (ratio = 1.146) is the largest, followed by the pair a & a* (ratio = 1.122). Regarding the absolute value of brightness, the brightness of a (2.396 × 10⁻² Jy Beam⁻¹) is larger than that of c (1.982 × 10⁻² Jy Beam⁻¹). In addition, there is a ridge extending from the central bright point (P) toward feature a. The direction of this ridge (PA = −45°) is completely different from the direction of multiple ridges seen in the entire image, and no other ridge shows the same direction. Based on these characteristics, we decided to continue the hybrid mapping by adopting two points for the next image model. In other words, we chose to use a model with a point at each of the two locations, the center and the location of feature a.

Table 3. Peak Brightness of the Features That Appeared in the First Step Image

Name	Brightness (Jy beam−1)	Name	Brightness (Jy beam−1)	Ratio	Order
P	8.736 × 10−2
a	2.396 × 10−2	a*	2.136 × 10−2	1.122	2
b	2.438 × 10−2	b*	2.187 × 10−2	1.115	3
c	1.982 × 10−2	c*	1.729 × 10−2	1.146	1
d	2.099 × 10−2	d*	1.889 × 10−2	1.111	4
e	2.780 × 10−2	e*	2.534 × 10−2	1.097	5
f	2.420 × 10−2	f*	2.257 × 10−2	1.072	7
g	2.368 × 10−2	g*	2.205 × 10−2	1.074	6
h	2.334 × 10−2	h*	2.282 × 10−2	1.023	9
i	2.364 × 10−2	i*	2.337 × 10−2	1.012	10
j	2.475 × 10−2	j*	2.404 × 10−2	1.030	8
k	2.328 × 10−2	k*	2.319 × 10−2	1.004	11

Note. The 11 features near the center are shown. P is the peak in the center that was replaced by the restoring beam. Features other than P are as shown in the brightness distribution of the residual dirty map.

Download table as:  ASCIITypeset image

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 The Event Horizon Telescope Collaboration et al. (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 M87, 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 for a wrong image model. It is safe to try the self-calibration of the amplitude solution after a 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 the 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 shows 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 μas are merged into one component. In the images of the last two days, adjacent CLEAN components within 5 μas are merged into one component.

Zoom In Zoom Out Reset image size

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 FOV) 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 ± 37°.93 (the residual of the EHTC ring is −0°.73 ± 45°.33), and the residual of our map for the last two days of data is 1°.79 ± 42°.11 (the residual of the EHTC ring is 4°.22 ± 45°.74). Here, we integrated the 5 minutes 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 Figure 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 the 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).

Zoom In Zoom Out Reset image size

Figure 5 shows the phase solution. Compared to the initial phase solution shown in Figure 2, there is no significant change in the overall structure. 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 the Submillimeter Array in Hawaii (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, the errors in amplitude are due to noise in the atmosphere and in the receiving system. In addition, the 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.

Zoom In Zoom Out Reset image size

For the 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 station clocks and atmospheric variations, which are more difficult to measure correctly than other error factors.

Amplitude solutions for AA, AP, PV, and SM stations are within 50% fluctuation (∼1.0 ± 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 JC station shows a large amplitude value at T ∼ 4.25 hr on the last observation day (101 days). On the other hand, for the LMT station, the amplitude is large for several times as follows:

(a) T ∼ 1 hr and T ∼ 2.6 hr on the first observing day (095 days),

(b) T ∼ 1 hr on the second observing day (096 days),

(c) T ∼ 2.25 hr and T ∼ 6.25 hr on the third observing day (100 days).

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 The Event Horizon Telescope Collaboration et al. 2019d).

These large amplitude solutions might suggest that the resultant image is seriously 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 for the resultant image quality, the results obtained by both the EHTC and our work should be rejected. The EHTC considered the fact that some stations require large-amplitude corrections during data analysis. EHTC then analyzed the data from 3C 279, which was observed with M87, 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 M87 (The Event Horizon Telescope Collaboration et al. 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. Additionally, the fact that the EHTC obtained a nonring 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 M87 data, does not mean that the ring image of the EHTC is the correct image of M87.

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 μas square). As noted in Section 3.4, since the data calibration is not yet complete, our final images show the sidelobe structures around actual features. This is a common phenomenon in synthesis imaging with radio interferometers with only a small number of stations. The images in Figure 7 show that “the unresolved VLBI core” in M87 has finally resolved into substructures.

Zoom In Zoom Out Reset image size

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 μas apart from the core C. The flux densities from the obtained CLEAN components are F_C ∼ 60 mJy for the core (C), F_K ∼ 40 mJy for the knot (K), and F_W ∼ 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 μas in diameter, and calculate the brightness temperatures (lower limit). The average brightness of feature C is T_b = 1.1 × 10¹⁰ K. Feature K has a brightness of T_b = 7.1 × 10⁹ K. Feature W is T_b = 4.7 × 10⁹ K. Thus, we have detected the central features with brightness temperatures higher than the EHTC ring (up to T_b ∼ 6 × 10⁹ 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 μas in diameter, the brightness temperature of core C reaches 10¹¹ K. If this is the case, the brightness temperature is an order of magnitude higher than the previous measurement cases (Akiyama et al. 2015; Hada et al. 2016; Kim et al. 2018). 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 AGNs (Horiuchi et al. 2004; Homan et al. 2006), 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 M87 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 radiatively inefficient accretion flow (RIAF) disks (Kato et al. 2008; Nakamura et al. 1997). Our observational results are consistent with those of previous studies, supporting the existence of the RIAF disk in the M87 core (Di Matteo et al. 2003).

There is a clear difference between the two images observed over the five days. According to The Event Horizon Telescope Collaboration et al. (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 (The Event Horizon Telescope Collaboration et al. 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 the 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 Δα = −0.4 μas, Δδ = −4.5 μas in 5 days, and the proper motion is 0.33 mas yr⁻¹ (v ∼ 0.1c). 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 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° as a whole. The position of feature C has hardly changed, except for 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° 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.  
  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 M87. Assuming that the position of SMBH is approximately equal to the position of feature C (the distance between the core of M87 and SMBH is 41 μas or 6 R_S; Hada et al. 2011), feature W would be located (or at least projected) at 12 R_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 toward us (black hole echo; Saida 2017; Virbhadra 2009; Virbhadra & Ellis 2000). 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.  
  Another central black hole. Feature C is the primary SMBH of M87, and feature W is a secondary SMBH orbiting the primary SMBH. If there is a binary SMBH in M87, it can be permanently observed with the EHT array. Based on these observations, we calculated two possible orbits.

  • (a)  
    The proper motion of feature W (μ = 0.34 mas/yr, v ∼ 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 μas from feature C. In this case, the orbital period T is ∼1.5 yr, If the real radius is R = 1.4 × 10³ au, the mass of central object M_c is only 1.2 × 10⁶ M_⊙. Since the estimated mass is too small as compared to those of the previous M87 studies, this assumption must be rejected.
  • (b)  
    We assume that the observed proper motions and structure change of feature W are only due to the 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 a SMBH with a mass of 6 × 10⁹ M_⊙ and that the orbital radius of feature W is the distance between features C and W. The distance between them is 83 μas (1.4 × 10³ au or 11.9 R_S). It is consistent with the 86 GHz core size of ∼80 μas at 86 GHz observed in 2014 (Hada et al. 2016), suggesting that the two features C and W are not transient. Also the sinusoidal oscillations of the position angle of the jet were observed with a period of roughly 8 to 10 yr (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 × 10³ 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 (∼17°), the orbital period of the binary can be ∼10 yr.
• 3.  
  Unstable initial knot. Feature W is another knot moving toward a different direction. The jet of M87 is known to have a wide opening angle at scales well below 1 mas (Junor et al. 1999). Furthermore, Walker et al. (2018) found evidence from 43 GHz observations that the initial opening angle θ_app is ∼70°. We found the angle ∠KCW is 70°, and further that the line of the average jet axis (PA = −72°) divides this angle almost evenly into 34°, and 36°. Furthermore, the lines CK and CW extend in the directions of PA = −38° and PA = 252°, respectively. These directions are very similar to those of the ridges observed at 43 GHz from where Walker et al. (2018) measured the initial opening angle. We guess that not only feature K but also feature W are initial knots that have just emerged from the core; and still, the shape is very unstable and shows large and rapid changes.

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 μ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 the so-called jet launching region, which is a few milliarcseconds away from the core (Section 4.2.2).

4.2.1. The Overall Structure

Figure 8 shows the overall structure of the M87 we obtained. In order to make the emission obvious, we used a restoring beam of 200 μas circular Gaussian, 10 times larger than the default beam size. As already mentioned, it is found that the emission at 230 GHz comes not only from the central point source but also from other regions.

Zoom In Zoom Out Reset image size

The EHTC apparently assumed or concluded that there is no bright source outside the narrow region (128 μas in diameter) where the ring was found. However, we found that the emission was not from such a narrow range, but from a wide range of more than a few milliarcseconds. This is consistent with the results of VLBI observations at 43 GHz and 86 GHz.

Our final image shows a similar structure to the average image in the 43 GHz band (the inset of Figure 8). There are two main similarities.

First, as in the 43 GHz image, our image shows that the jet has an extended structure leading to the core. Then, up to a few milliarcseconds away from the core along the jet axis, both edges are bright, as in the previous observations.

Second, the brightness distribution of the jet in the 230 GHz is consistent with the trend of those obtained from lower-frequency observations. The core is vastly brighter than the jet structure. Within the radius of 0.25 mas (250 μas) from the center, 63% (the first two days’ image) and 75% (the last two days’ image) of all obtained flux densities are concentrated. However, features C, K, and W (several tens millijanskys at most, see Table 5) do not occupy them, rather the flux densities are distributed over a wider area. In contrast, the EHTC rings have a total of about 500 mJy that is contained entirely within a diameter of only a few tens of microarcseconds.

Table 5. Properties of Main Features: Positional Offsets from the Map Phase Center in μas, Flux Densities in mJy, and Minimum Brightness Temperatures in Kelvin Calculated with the Assumption That the Emission Area Is 15 μas in Diameter

	The First Two Days’		The Last Two Days’
Peak Position	R.A. (μas)	δ (μas)	R.A. (μas)	δ (μas)
Core	−1.8 ± 0.6	−1.5 ± 0.6	−2.3 ± 0.5	−1.5 ± 0.5
Knot	−22.2 ± 0.5	27.0 ± 0.5	−23.1 ± 0.3	22.5 ± 0.3
West	−78.1 ± 0.3	−33.0 ± 0.3	−77.2 ± 0.3	−37.5 ± 0.3
			Δ R.A. (μas)	Δδ (μas)
Core			−0.5	0.0
Knot			−0.9	−4.5
West			0.9	−4.5
Position of intensity increase			R.A. (μas)	δ (μas)
Core area a			−22.4 ± 0.4	15.0 ± 0.4
Knot area a			−29.8 ± 0.7	40.5 ± 0.7
b			−40.8 ± 0.3	33.0 ± 0.3
c			−48.0 ± 0.4	43.5 ± 0.4
West area a			−74.7 ± 0.7	−43.5 ± 0.7
b			−79.9 ± 0.5	−15.0 ± 0.5
c			−99.8 ± 0.7	−25.5 ± 0.7
Integrated intensity	(mJy)		(mJy)
Core	55.6 ± 5.2		66.1 ± 4.7
Knot	33.5 ± 2.7		44.9 ± 2.3
West	22.5 ± 1.2		30.2 ± 1.4
Brightness temperature	(K)		(K)
Core	>1.0 × 1010		>1.2 × 1010
Knot	>6.0 × 10 9		>8.1 × 10 9
West	>4.0 × 10 9		>5.4 × 10 9
Grouped-CLEAN components	(mJy)		(mJy)
FGCC > 0.1 mJy	707.4 (n = 1151)		1032.6 (n = 1657)
All	767.8 (n = 2824)		1154.6 (n = 7844)

Note. Positions of intensity increase in features are also shown. At the bottom, the sum-up intensities and the numbers of the grouped-CLEAN components in the whole images are shown.

Download table as:  ASCIITypeset image

The results of this observation at 230 GHz show that the brightness in the jet region is orders of magnitude lower than that in the core region. In addition, the decay of the flux density along the jet is more rapid at 230 GHz than in the lower-frequency observations. Compared to the peak luminosity of the core, the relative intensities are 6.6 × 10⁻² at 0.25 mas from the core, 9.9 × 10⁻³ at 0.5 mas, and 2.3 × 10⁻³ at 1 mas (Figure 9). While, in the observation at 43 GHz, the decreases of intensity are 2.8 × 10⁻¹ at 0.25 mas from the core peak, 8.5 × 10⁻² at 0.5 mas, and 2.5 × 10⁻² at 1 mas (from the upper panel of Figure 6 in Walker et al. 2018). At 1 mas position, the intensity of the jet structure is 2.5% of the core peak at 43 GHz; however the intensity at 230 GHz is only 0.2% of the core peak, namely the intensity of the jet structure is greatly attenuated at 230 GHz. However, with respect to the structure of the brightness and intensity distribution of both edges, the trend is in good agreement with the previous results of the M87 jet.

Zoom In Zoom Out Reset image size

The total flux density measured by the EHTC (The Event Horizon Telescope Collaboration et al. 2019c) was 1.12 and 1.18 Jy, for the first two days and the last two days, respectively. In contrast, the total flux density of the CLEAN component in our analysis is 767.8 and 1154.6 mJy, respectively. That is, there are missing flux densities of 353 and 25.4 mJy, respectively. ¹²

The difference between the total flux density of our image and the single-dish flux density is most likely due to the presence of extended emission somewhere, which the present EHT array cannot detect. As shown in Appendix A, the EHT array cannot detect a smoothed emission feature (like Gaussian brightness distributions) with size larger than 30 μas.

4.2.2. The Jet Launching Region

In this section, we present the structure within a few milliarcseconds from the core. We have found emission belonging to the jet component that was not detected by the EHTC.

In Figure 10, the brightness distribution in this region is represented in three ways. The logarithmic pseudo-color (arbitrary unit) is used to represent the large differences in brightness distribution. The upper panels (a) and (b) are composed by restoring beams of a circular Gaussian with half-power beam width (HPBW) of 20 μas, which corresponds to the size of the spatial resolution of the EHT array for M87 observations. The middle panels (c) and (d) are composed by restoring beams of circular Gaussian with HPBW of 200 μas. In order to facilitate a comparison with previous results, the beam size is close to the spatial resolution of previous lower-frequency observations (43, 86 GHz). Panels (e) and (f) in the bottom row show the image by the large restoring beam overlaid with the 43 GHz averaged image by the VLBA (Walker et al. 2018). Note that the 43 GHz image is time-averaged over 17 yr, so the knot-like features have been averaged out. It can be seen that the brightness distribution at 230 GHz is consistent with that at 43 GHz. The left panels (a), (c), and (e) show images of the data from the first two days. Panels (b), (d), and (f) on the right show images of the data from the last two days. Obviously, they are different from each other. However, the differences seen in the regions of a few milliarcseconds cannot be attributed to the intrinsic variability of the source that occurred during the five days. Rather, it seems to be mainly dependent on the observational conditions.

Zoom In Zoom Out Reset image size

The emission areas of our 230 GHz results are consistent with that of the 43 GHz average image. They also show that both edges of the jet are brightened, which have been observed in 43 and 86 GHz. Based on our data analysis, it seems that the detection of emissions in the range of several milliarcseconds from the core of the M87 has been successful to some extent.

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.

4.3.1. From Sensitivity

The Event Horizon Telescope Collaboration et al. (2019c) shows that the typical sensitivity of a baseline connected to ALMA is ∼ 1 mJy. We estimate that this sensitivity is for an integration of about 5 minutes. For an on-source time of 2 days, the attainable sensitivity reaches close to 0.1 mJy (ALMA-LMT baseline, S/N = 4, assuming a point source). It is difficult to estimate the practical sensitivity of the synthesized image of an interferometer composed of antennas with different performances, such as the EHT array. However, it is unlikely that the image sensitivity will not be worse than the baseline sensitivity noted above.

Here, we consider 0.1 mJy to be a reliable detection limit for our final images. Figure 11 shows the distribution of the grouped-CLEAN components with flux densities larger than 0.1 mJy. Almost all of the components are concentrated within a few milliarcseconds of the core. (The remaining components are located in a false bright spot created outside the range of this figure, about 20 mas west of the center.) The image from the core to a few milliarcseconds along the average jet axis seems to be reliable in terms of detection limits.

Zoom In Zoom Out Reset image size

A large number of grouped-CLEAN components with flux densities larger than 0.1 mJy are found in our final images. The number of grouped-CLEAN components with flux densities larger than 0.1 mJy is 1151 from the images of the first two days (with a sum of flux density of 707.4 mJy), and 1657 from the images of the last two days (with a sum of flux density of 1032.6 mJy; Table 5).

4.3.2. The Robustness of Our Final Image

In this section, we discuss another property of our final images: the robustness of the image structure.

If the data is not yet completely calibrated, BOX technique is effective. As well as the EHTC, we also used the BOX technique to limit the imaging area (FOV). This technique has the potential to produce good images even if the calibration is incomplete. On the other hand, it may create structures that do not actually exist. In fact, the bright spot on the right-hand side of our final image (Figure 8) is such an example. Therefore, care must be taken when using the BOX technique, because a false structure will be created in the BOX area, and the real structure outside of the BOX area will be removed from the image.

We examine how the image is affected when we change the BOX parameters. In other words, we investigate the stability of the image. We compare the images obtained by changing the size of the BOX. The panels in Figure 12 show the four cases. The upper left panel (a) is the image without BOX (the FOV is 24.576 mas square). The top right panel (b) shows the image with the same 8 BOXes that we used to obtain our final images. The lower right panel (c) is the image with a small BOX (circle with diameter D = 5 mas is used). The lower right panel (d) is the image with a very narrow BOX (circle with diameter D = 128 μas that corresponds to the FOV the EHTC used). These four CLEAN images were produced using data of the entire four days.

Zoom In Zoom Out Reset image size

In all panels, the emission can be seen at the positions of features C (core) and K (knots). On the other hand, feature W disappears in the case where the BOX setting is omitted (no BOX case). In the case of the EHTC FOV, no emission is seen at the position of feature W because the position of feature W is outside the BOX setting.

Without the BOX setting, the S/N of the image is degraded. From the comparison between panel (a) and the other panels, we can see that the BOX setting compensates for the lack of calibration and improves the image quality. Thus, the presence or absence of the BOX setting seems to have an effect on the image quality. Another noteworthy point is that, in the case of the very narrow BOX setting (panel (d)), several different bright spots newly appear in the BOX. Moreover, some of them are located at the boundaries of the BOX. In such a case, other actual brightness distributions could exist outside the BOX setting.

Since feature W disappears in the CLEAN image without the BOX setting, feature W is considered to be less reliable than features C and K. As mentioned at the end of Section 4.1, there are other reasons to consider that feature W has a real existence.

4.3.3. Self-consistency of Our Imaging as Compared to Those of the EHTC Images

At the end of this section on image reliability, we present the degree of matching between the visibility and the image model. Here, we compare the results with those of the EHTC ring.

• 1.  
  Relations between the visibility amplitude and u−v distance (projected baseline length). The amplitude of visibility obtained by the inverse Fourier transform of the image model is compared with those of the observed visibility data. Figures 13, and 14 correspond to Figure 12 in The Event Horizon Telescope Collaboration et al. (2019d). This kind of comparison of visibilities is often performed to check the reliability of an image. However, here, the observed visibility data are calibrated by self-calibration solutions using the image model. Therefore, it is important to note that the amplitudes of the observed visibility data and those from the image model are no longer independent of each other. What can be safely determined from this comparison is the internal consistency of the imaging and calibration process. Figure 13 shows the data from the first two days, and Figure 14 shows the data from the last two days. The top row of each figure shows the variation of the visibility amplitude with respect to the projected baseline length. The red dots are those of the image model. The middle and bottom panels show the normalized residual amplitudes between the image model and the calibrated observation data. The plotted points are the calibrated raw data that have not been time-integrated, and we can see that the scatter is much larger than that in their Figure 12 (The Event Horizon Telescope Collaboration et al. 2019d), where the time-integrated points have been plotted. We can see that the average and standard deviation of the normalized residuals of our final image are much smaller than that of any of the EHTC ring images in Figure 15. As an example, we show the normalized residual values for t = 180 s integration.For the data of the first two days, our image shows NR_ours = 0.030 ± 0.539, while the EHTC images show NR_EHTC = 0.148 ± 0.933.For the data of the last two days, our image shows NR_ours = 0.127 ± 1.259, while the EHTC images show NR_EHTC = 0.589 ± 2.370. Here, in the case of EHTC, we used the simple averages of those values for the four EHTC images. One thing that interests us is the large discrepancy in amplitude of the EHTC ring image cases at the longest baseline lengths over 8 × 10⁹ λ. It is three times larger than those of our final image cases. Since the EHTC ring images are very compact, if the images are really correct, the amplitude residuals at the longer baseline should become small at least. Another thing that interests us is the amplitudes at the very shorter baselines that are nearly zero λ. They contain the components of the extended structure that are resolved by high spatial resolution by EHT, so it is not surprising if they do not match. Our images reproduce the amplitudes of the very short baselines well. The differences are more significant in the cases of the EHTC rings. In our cases, the normalized residuals are 4 at most, but in the cases of the EHTC rings, they are distributed widely in the range of 0–15, which is not surprising since the EHTC rings are compact and have no extended components. However, in Figure 12 in The Event Horizon Telescope Collaboration et al. (2019d), the maximum is 4, as if the result shows good self-consistency. The EHTC Figure also shows the same results for the normalized residuals at the longest baselines. This is not consistent with our own analysis. Perhaps a different integration time of the data may cause this apparent discrepancy. (There is no explanation for the integration time of the data points in Figure 12 in The Event Horizon Telescope Collaboration et al. (2019d). Since the scatter of data points is affected by thermal noise, its value changes depending on the integration time of the data. Therefore, we examined the amounts of normalized residuals by changing the integration time. Figure 15 shows the average and standard deviation of the obtained normalized residuals. It can be seen that, at any integration time, our final image always shows smaller values than that of the EHTC ring image. The averages of the 10 s integrations and integrations over 180 s differ by a factor of 3, which explains the discrepancy above.The diagram of the visibility amplitude and u−v distance shows that our final images, both the first and the last days, show a better consistency of imaging and calibration processing compared to the cases of the EHTC ring images. The diagrams indicate that our images, not those of EHTC, are supported by the data.
• 2.  
  Closure phase variations. Following Figure 13 of The Event Horizon Telescope Collaboration et al. (2019d), we show the closure phases of the observed data and those derived from image models for some triangles in Figure 16. We added the closure phases of ALMA-LMT-SMA and LMT-PV-SMA to the three triangles (ALMA-LMT-PV, ALMA-SMT-LMT, SMT-LMT-SMA) shown by the EHTC. Closure phase is a quantity that is free from systematic phase errors and reflects the observed source structure. All panels in Figure 16 show large phase variations, which correspond not to time variation in the structure of the observed source, but to time variation of the shape of the triangle composed by the three stations as seen from the observed source. The green dots are the closure phase corresponding to the EHTC ring image, and the red dots correspond to our image. The dots of our image (green dots) appear to be better aligned with the observed data than those of the EHTC ring image. Our image is more complex than the EHTC ring, resulting in short-term small closure phase variations.The three panels from the top right toward the bottom correspond to the panels shown in Figure 13 of The Event Horizon Telescope Collaboration et al. (2019d). In the case of the two triangles ALMA-LMT-PV and SMT-LMT-SMA, our results are consistent with those of EHTC. However, our results for the closure phase in ALMA-SMT-LMT triangle differ from those of EHTC. In our case, the closure phase shows an increase from +25° to +85°, while that of the EHTC shows a decrease from −25° to −80°. Both the first and last values and the amount of change are opposite in positive and negative. All triangles were examined, but none were identical to the closure phase variation shown by the EHTC for the ALMA-SMT-LMT triangle. Our closure phase values are from the AIPS task, CLPLT. All triangles were examined, and none showed the variation similar to that the EHTC showed for the triangle. Also, there are no significant closure phase discrepancies between the real and model data. There seems nothing wrong with the CLPLT calculations. In our analysis, there is no clear difference in closure phase matching between our images and the EHTC rings. A notable difference is in the case of the LMT-PV-SMA triangle, where the closure phase in the EHTC rings is beyond ±3σ error bars, whereas in our final images, it manages to fall within it.The values of the closure phases also change depending on the integration time of the data; however, even when the integration time is changed, the residuals of either of them do not become overwhelmingly small. Figure 17 shows the statistics of the closure phase residuals for all triangles. As far as the closure phase residuals are concerned, between our images and the EHTC rings, there is no significant difference. Our image of the core-knot structure shows the same magnitude of closure phase residuals as those of the EHTC ring image. As an example, we show the standard deviations of the closure phase residuals at t = 180 s integration. For the data of the first two days, our image shows σ_ours = 40°.5, while the EHTC image shows σ_EHTC = 38°.5. For the data of the last two days, our image shows σ_ours = 43°.2, while the EHTC image shows σ_EHTC = 43°.7. As for closure phase residual, there is no significant difference between ours and the EHTC rings. If we claim that the EHTC ring image is correct due to the closure phase residual, then our images are also correct. ¹³ If these residuals are due to thermal noise, they should decrease inversely proportional to the root square of the integration time T, but as Figure 16 shows, they do not decrease in that way. It means that both images of ours and those of the EHTC still have differences from the true image.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We here investigate the u−v sampling of the EHT array for M87 observations. The u−v coverage itself is shown in the EHTC Paper I (The Event Horizon Telescope Collaboration et al. 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 M87 during the four observational days. The number of sampled data is plotted against the fringe spacings in μas unit. We can see that there are no samples for ranges d = 26–28 μas, d = 44 −46 μas, and d = 95–100 μas. ¹⁴ The number of sampled data for the size of the EHT ring (d = 42 ± 3 μas) is quite small, and for spacing between d = 44–46 μas, there are no sampled data.

Zoom In Zoom Out Reset image size

The few numbers of sampling around 42 μ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 examine the dirty beam calculated from the u−v coverage. Dirty beam is the term used in the field of radio interferometry, and its meaning is the same as that of PSF 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 δ-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 M87. The FWHM of the main beam is approximately 20 μ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 μas. ¹⁵ 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 ± 3 μas). Thus, it is important to clarify whether or not this PSF structure has affected the ring image or not. In addition, in Appendix B we show 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 ∼40 μas spacings.

Zoom In Zoom Out Reset image size

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 clear that, in the case of multielement 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 of 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 that look like the ghosts of this central ring-like structure.

Zoom In Zoom Out Reset image size

This is not a blurred intrinsic image of M87, 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 μas radii centered at (+2 μas, +22 μas) from the map center, which is the BOX setting as used in the EHTC-DIFMAP pipeline (Section 5.5).

Zoom In Zoom Out Reset image size

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 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 μ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 ∼40 μ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 central symmetry, one might think 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 CLEAN map 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 noncenter symmetry it, however, 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 M87, 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 simulated 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 M87 observations. As for the 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 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 μas (model B). Model B is a pair of points of 0.5 Jy located at (0 μas, 0 μas) and (0 μas, +45 μ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 μas radii centered at (+2 μas, +22 μ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.

Zoom In Zoom Out Reset image size

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 two points appearing show small shifts from the positions in model B. The center point shifted toward the east, while the upper point shifted toward 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 toward 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 the u−v coverage is sufficient, 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 ∼40 μas. We have shown that the shape of the PSF (essentially the effect of the sampling bias that makes it easy to create ∼40 μas structures) has the power to create ∼40 μas ring structures even from data with only one central point structure (centrosymmetry). However, the most powerful factor is not the 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 ∼40 μ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. ¹⁶ Among the three methods the EHTC used for the imaging, DIFMAP is the closest to the usual procedure ¹⁷ in the above web page. The 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 analysis 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 the 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 one-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 that 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. This technique may give us a good image, even when the calibration is incomplete. As already mentioned, we used the BOX technique in our analysis. The EHTC-DIFMAP analysis did the same.

Now we investigate the EHTC-DIFMAP procedure. The EHTC-DIFMAP imaging analysis 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 ¹⁸ in the first imaging trial. The CLEAN subtraction area was restricted to within the circle of 60 μas diameter specified by the BOX.

The BOX is roughly a circle of diameter 60 μ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 μas in the y-direction (δ direction). This offset of 22 μ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 ∼40 μas ring.

Zoom In Zoom Out Reset image size

The PSF structure within the BOX explains not only the diameter of the EHTC ring, which is ∼40 μ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 ∼40 μas ring is the result of such enhancement of ∼40 μ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 extended 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° and PA = 225° directions from the center. These are probably sidelobes, but it is very interesting that one emission peak in the PA = 225° 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.

Zoom In Zoom Out Reset image size

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 direction (left panels), (2) a two-point model aligned in north–south direction (center panels), and (3) a 40 μas diameter ring, but the ring width is 2 μas, which is different from that of the EHTC ring (right panels).

Zoom In Zoom Out Reset image size

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 μ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 FOV 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 we have performed many simulations and what are shown in Figure 25 are selected samples. The summary of our simulations is as follows.

• 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.  
  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.  
  Core-knot model. Two-points model consisting of a point of 1 Jy at the center (0 μas, 0 μas) and a point of 0.25 Jy at (−25 μas, +20 μas). The distance between them is about 32 μas, and they are located along a line in the direction PA = −38°.7. 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.  
  Ring images with diameters around 40 μas. The ring width is 2 μas, and the ring diameters are 30, 35, 40, and 45 μ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 μ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.  
  Two-point model (shift in the north–south direction). Place a point of 1 Jy at the center (0, 0 μ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 μas. Using IMAGR in AIPS, two points cannot be separated and resolved when the spacing is very narrow. At a separation of 10 μas, the image is elongated in the north–south direction. If the separation is larger than 25 μas, the two points are reproduced separately. Using the EHTC-DIFMAP pipeline, when the separation is larger than 45 μas, the two points are reproduced correctly, but otherwise the north point is not reproduced in the correct position.
• 6.  
  Two-point model (east–west shift). A point of 1 Jy is placed at the center (0 μas, 0 μas), and a point 0.5 Jy is placed at δ = 42 μas. The R.A. position of the second point is shifted every 5 μas from −20 to 20 μ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 μas.

Thus, in most cases, the EHTC-DIFMAP pipeline did not reproduce the intrinsic image of the simulation data. Figure 10 of The Event Horizon Telescope Collaboration et al. (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 reproduces the input model images reasonably well, 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 early 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 analysis case). This is an insufficient description of data calibrations because it is rare that only a-priori calbrations 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 the 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 definitely correct, 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 the sensitivity of the ring structure obtained by the EHTC to that from 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 μas, (b) 80 μas, (c) 100 μas, (d) 120 μas, (e) 240 μas, and (f) 450 μas. Their centers are located at the same position as that of the EHTC-DIFMAP pipeline case. Other parameters of the CLEANs are the same as those used for the right image in Figure 26.

Zoom In Zoom Out Reset image size

The results are shown in Figures 29 and 30. Figure 29 shows views of the entire imaging area of 2 mas square, and Figure 30 shows enlarged views of the central 256 μas square. The BOXes used are indicated by blue circles. The ring structure in panel (a) (D = 60 μ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 μas restoring beam is convolved to form the CLEAN map. In the cases of (b) D = 80 μas and (c) D = 100 μ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 μ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 μas and (f) D = 450 μas, there is no ring. Instead, one elongated bright spot appears. We also tried larger BOX cases, D = 900 μas, and D = 1300 μas, and the resulted image is similar to the case of D = 450 μas. They show a bright spot at the center of the map.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If the FOV is smaller than that for 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 μas. It is an odd coincidence that the EHTC imaging analysis set the FOV to be 128 μas (The Event Horizon Telescope Collaboration et al. 2019d), just around the boundary FOV. (It is 60 μas in the case of the EHTC-DIFMAP pipeline as already mentioned.) ¹⁹

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 The Event Horizon Telescope Collaboration et al. (2019d), in which we found some of the reasons why they created an artifact. The EHTC conducted various surveys, but their methods were not objective and were biased from the very beginning of their analysis.

5.8.1. Attention to Sampling Bias

5.8.2. The EHTC Simulations

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 M87, i.e., the missing ∼40 μas spatial Fourier components (Section 5.1). The sampling bias tends to create an artificial structure of ∼40 μas in imaging. We confirmed that the artificial ∼40 μ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 doing so. Also, 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 ∼40 μas of the sampling bias origin. In addition, the u−v sampling bias effects are independent of the need for data calibration. Therefore, even with the 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 ∼40 μas in size from the beginning, and thus performed the image synthesis with a very narrow FOV (128 μas in size). A large FOV setting will not enhance the artificial ring image of ∼40 μas size. On the other hand, the very small FOV set by the EHTC enhances the sampling bias effect to create the ∼40 μ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 the initial analysis of the EHTC.