EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 696 (April 2025) A&A, 696 (2025) A169 Full HTML
Open Access
Issue
A&A
Volume 696, April 2025
Article Number A169
Number of page(s) 11
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/202452038
Published online 21 April 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model.

Open Access funding provided by Max Planck Society.

1 Introduction

Continued improvements in instruments and algorithms for very long baseline interferometry (VLBI) enable us to image and analyze radio emission from near to supermassive black holes (SMBH) and inside the formation zone of relativistic jet in active galactic nuclei (AGNs). The recent Event Horizon Telescope (EHT) observations captured the black hole “shadow” in M87* and SgrA* (Event Horizon Telescope Collaboration 2019a, 2022), and the Global mm-VLBI Array (GMVA) and Very Long Baseline Array (VLBA) observations provided high fidelity images of the extended jet emission in radio-loud AGNs (Walker et al. 2018; Kim et al. 2018; Okino et al. 2022). The nearby galaxy M87, with a SMBH mass of 6.5 × 109 M⊙ (Event Horizon Telescope Collaboration 2019b) and a redshift of z = 0.004283 (Cappellari et al. 2011), provides a unique opportunity to study the black hole accretion disk and jet launching mechanism simultaneously due to the large angular scales of 0.08 pc or 260 Rg (where Rg = GM/c2) per milliarcsecond. So far, the black hole shadow and relativistic jet of M87 at radio frequencies have been observed independently due to instrumental limitations and the characteristics of the source in different frequency regimes.

Observations of M87 with the GMVA, the phased Atacama Large Millimeter/Submillimeter Array (ALMA), and the Greenland Telescope (GLT) in 2018 enabled the central ring-like structure and the extended jet to be imaged simultaneously. Lu et al. (2023) employed a CLEAN (Högbom 1974; Clark 1980) self-calibration and imaging method in order to reconstruct the ring and jet structures in M87 with a large field of view. Additionally, the ring in M87 was reconstructed with a smaller field of view using a regularized maximum likelihood (RML)-based imaging software SMILI (Akiyama et al. 2017). GMVA observations at 86 GHz play a significant role in resolving the core and extended jet emission of M87. Prior to the EHT M87 observation in 2017 (Event Horizon Telescope Collaboration 2019a), the GMVA M87 observation in 2014 and 2015 (Kim et al. 2018) provided the core and edge-brightened jet of M87 image down to 13 Rg. The GMVA observations with ALMA in 2018 facilitated the detection of the ring with extended jet emission since ALMA provides the longest north-south baselines with improved sensitivity. For the first time, these finding connected the central ring-like feature, which presumably corresponds to the compact accretion flow, with the innermost jet. This work was accompanied by recent groundbreaking observations of the jet in M87 on a variety of scales, ranging from the (full polarimetric results on) horizon scales (Event Horizon Telescope Collaboration 2019a, 2021a, 2023, 2024), to the strongly edge-brightened innermost jet (Walker et al. 2018; Kim et al. 2018, 2023), as well as to the triple-peaked, helical structure and dynamics in the large scale jet (Asada et al. 2016; Hada 2017; Nikonov et al. 2023; Cui et al. 2023). It is imperative to connect the compact scale structures of M87* to the large scale structures to understand AGN jet formation. The observations presented by Lu et al. (2023) play a significant role in connecting horizon scales to jet scales.

However, the data reduction process of GMVA+ALMA is challenging due to the sparsity of the UV coverage, tropospheric phase corruption at millimeter wavelengths, a low signal-to-noise ratio (S/N), and inhomogeneous antenna statistics due to the different sensitivity of antennas. This difficulty jeopardizes the interpretation of some of the features in the reconstruction, namely the ring-like structure (two brighter spots) and the jet (inner ridge line). Imaging with the conventional CLEAN method is not able to reconstruct a robust ring-like structure due to suboptimal resolution. Furthermore, the simple assumption in CLEAN that the sky brightness distribution is a collection of point sources is not an optimal prior to describe the extended continuous jet emission. An independent assessment of the robustness of these features is needed to facilitate a much-awaited scientific interpretation.

The latest advancement in forward modeling imaging algorithms enables us to generate more robust results from sparse millimeter-VLBI data sets. As an example, Bayesian imaging is a probabilistic approach that reconstructs the posterior distribution using Bayes’ theorem. Hence, Bayesian imaging is able to use posterior samples to estimate the uncertainty of parameters, such as image features and instrumental gains. However, image reconstruction is computationally demanding compared to other imaging methods since a collection of possible images is reconstructed instead of a single image. RML methods reconstruct an image by minimizing an objective function, which consists of a data fidelity term and regularizers that favor certain image structures such as smoothness or sparsity. Contrary to CLEAN, forward modeling methods fit the model to the data directly in the visibility domain. We can fit closure quantities directly or calibration can be incorporated in the image reconstruction. Furthermore, we are able to encode knowledge about the source and measurement setup explicitly in the prior distribution and regularizers. As a result, both of the imaging approaches outperform the traditional inverse modeling CLEAN algorithm and we can reconstruct reproducible images with improved resolution in a less supervised fashion. A detailed comparison between CLEAN and forward modeling approaches can be found in Chael et al. (2016); Arras et al. (2021); Müller et al. (2024).

In this work, we reconstruct images using two imaging algorithms: we perform Bayesian self-calibration and imaging jointly using the Bayesian imaging software resolve (Junklewitz et al. 2016; Arras et al. 2022; Roth et al. 2023; Kim et al. 2024) and we reconstruct an image with closure amplitudes and closure phases only (Chael et al. 2018) using the RML-based DoG-HiT software (Müller & Lobanov 2022, 2023a,b). These two independent imaging methods are utilized in order to estimate the robustness of the M87 ring and extended jet emission from the GMVA+ALMA observation in 2018. They quantify the robustness of the recovered features from two alternative, supplementary perspectives, by self-calibration and imaging from a probabilistic point of view and by imaging with closure quantities only without a potential self-calibration bias.

This article is structured as follows. In Sect. 2, we explain the resolve and DoG-HiT image reconstruction methods. In Sect. 3, we show two image reconstruction results and then analyze the robustness of the M87 ring structure and jet emission. In Sect. 4, we summarize our results.

2 Method

2.1 Bayesian self-calibration and imaging by resolve

resolve1 is an open-source Bayesian imaging software for radio interferometric data (Junklewitz et al. 2016; Arras et al. 2022; Roth et al. 2023; Kim et al. 2024). In resolve, samples of potential images and antenna-based gain solutions that are consistent with the data are reconstructed by Bayes’ theorem in a variational inference sense (Blei et al. 2016; Knollmüller & Enßlin 2019; Frank et al. 2021). In this paper, we used the metric Gaussian variational inference method (Knollmüller & Enßlin 2019, MGVI) to estimate the posterior distribution of the sky brightness distribution and antenna-based gains. The MGVI method enables us to perform high-dimensional Bayesian inference with affordable computational resources. The probabilistic approach could be advantageous for the M87 GMVA+ALMA observations due to the sparse UV coverage, and large and heterogeneous data uncertainties. For details of the M87 GMVA+ALMA data, we refer to Sects. 1 and 2 of the supplementary information in Lu et al. (2023).

We reconstructed the resolve image in Fig. 1 with a spatial domain of 2048 × 1024 pixels and a field of view of 4 mas × 2 mas from a priori calibrated (without self-calibration) data. The Bayesian self-calibration is performed simultaneously with the imaging. The number of posterior samples was 100 and the reduced χ2 value of the final result was 1.1. The wall-clock time for the resolve reconstruction was 5.5 hours on a single node of the MPIfR cluster with 25 MPI (Message Passing Interface) tasks.

Since the GMVA+ALMA array is highly inhomogeneous, it is a reasonable assumption that each antenna gain has a different temporal correlation structure. In resolve, we utilized a Gaussian process prior with a nonparametric correlation kernel in the NIFTy software2 for the sky brightness distribution (image) and gain prior model. The spatial correlation between image pixels and temporal correlation between gains can therefore be inferred from the data without manual steering of the gain solution interval constraints (see Fig. A.2). The amplitude gain prior is assumed to be correlated and different correlation kernels are inferred per antenna. Gain amplitudes for right-hand circular polarization (RCP) and left-hand circular polarization (LCP) are assumed to have the same correlation structure to stabilize the self-calibration and image reconstruction. The phase gain is assumed to be uncorrelated since the phase coherence time in GMVA observations is comparable with the data averaging time (ten seconds) and it can be even shorter under poor weather conditions.

Furthermore, the posterior distribution of model parameters, such as each pixel in the image and antenna-based gain solutions, can be explored in resolve. In other words, the reliability of the reconstructed parameters and the image features, such as the ring structure and extended jet emission, can be quantified by estimated uncertainties from posterior samples. As an example, if one antenna is problematic, then it would result in a high uncertainty in its gain solution. As a result, in Bayesian imaging, the high uncertainty of these data points can be self-consistently taken into account in the image reconstruction. More details about the Bayesian self-calibration and imaging method for VLBI data and validation with synthetic data can be found in Kim et al. (2024).

thumbnail Fig. 1

Image reconstructions of GMVA+ALMA M87 at 86 GHz. Each image presents results obtained by a different algorithm. The top posterior mean image was reconstructed using the resolve Bayesian self-calibration and imaging method. The middle image was obtained using DoG-HiT closure amplitudes and closure phases only imaging. The bottom image was obtained using the CLEAN self-calibration and imaging in Lu et al. (2023). The CLEAN image is convolved with an elliptical beam, which is represented as an ellipse with sizes 79 × 37 μas, PA = −63∘ in the bottom-left corner.

2.2 Imaging with closure quantities by DoG-HiT

DoG-HiT is a regularized maximum likelihood (RML) imaging algorithm that reconstructs the image using multiscalar wavelet basis functions (Müller & Lobanov 2022). The basis functions are fitted to the UV coverage, offering a neat separation between covered and noncovered Fourier coefficients, that is, gaps in the UV coverage (for more details on the wavelets, we refer the reader to Müller & Lobanov 2023a). The image is recovered by a sparsity promoting forward-backward splitting framework that effectively calculates the multiresolution support. In other words, the image is represented by the set of all statistically significant wavelet scales. It has been demonstrated that the multiresolution support is beneficial prior information that enables reconstruction even for sparse and weakly constrained settings (Müller & Lobanov 2022, 2023b).

In this work, we aim to validate the results presented in Lu et al. (2023) using closure-only imaging pioneered by Chael et al. (2018). Closure-only imaging is a self-calibration independent technique. As a result, we are able to estimate the robustness of the recovered features against the gain calibration. We represent the DoG-HiT image (see Fig. 1) in a square field of view of 4096 μas by 512 × 512 pixels. The reconstruction was run on a CPU (11th generation Intel core i7) with 32 GB RAM for roughly thirty minutes. The reduced χ2 to the closure phases was 1.36 and to the closure amplitudes was 1.1. For the reconstruction with DoG-HiT, we first select a set of wavelet basis functions, called a dictionary, fitted to the UV coverage of the observation. Then we run DoG-HiT with all large scale wavelets that were significant to fit the extended, diffuse jet emission. Then we use this image as an initial guess and add all small-scale wavelets that are relevant to represent the central ring, and minimize the χ2-metric to the closure phases and closure amplitudes.

In this framework, we directly fit to the self-calibration independent closure amplitudes and closure phases. For a given number of antennas, more closure quantities can be constructed than visibilities. However, not all the measurements are independent since they can be represented as linear combinations of other closure triangles and quadrilaterals (Twiss et al. 1960; Blackburn et al. 2020; Thyagarajan et al. 2022). Since the number of statistically independent closure phases and amplitudes is therefore smaller than the number of independent visibilities, which effectively leads to a number of degeneracies such as lost information of the total flux density and the absolute source position, we have to account for the larger freedom in the models. There are two opposite strategies to achieve this. We could either try to explore the multimodality and degeneracies inherent to closure quantities, for example using recently proposed multiobjective optimization schemes (Müller et al. 2023; Mus et al. 2024a,b), or we could utilize a more constraining piece of prior information that resolves the degeneracies. For the purpose of the latter approach, multiresolution support proved successful and was implemented in DoG-HiT.

3 Results

The M87 GMVA+ALMA images at 86 GHz by resolve, DoG-HiT, and CLEAN (Lu et al. 2023) are shown in Fig. 1. resolve and DoG-HiT image fits and results are archived in zenodo3.

thumbnail Fig. 2

Ring of the M87 image reconstructions at 86 GHz from Fig. 1. The color shows intensity in Jy mas−2 according to the linear color bar located at the right of the figure. Each image presents results obtained by a different algorithm, whose names are indicated in the lower right corners. The left image is the resolve posterior mean image using the Bayesian self-calibration and imaging method. The left-center image is the DoG-HiT reconstruction using closures only imaging. The right-center image represents the CLEAN reconstruction with the over-resolved 37 μas circular beam in Lu et al. (2023). The right image shows the SMILI reconstruction in Lu et al. (2023). All images were processed by a Gaussian interpolation.

3.1 Estimation of ring-like features in M87

A zoom-in on the central compact emission region is presented in Fig. 2. In the visibility domain, the presence of the visibility null amplitude at around 2.3 Gλ and the phase jump around the null amplitude (see Fig. S2 in Lu et al. 2023) is analogous to EHT M87 observations in 2017 and 2018 (Event Horizon Telescope Collaboration 2019a, 2024), which is strong evidence of the M87 ring structure. Figure 5 (middle panel) shows the visibility domain representation of resolve and DoG-HiT images, namely the posterior mean model visibility amplitudes by resolve and the amplitudes of the Fourier-transformed DoG-HiT image. The visibility null amplitudes in resolve and DoG-HiT are located at around 2.3 Gλ, which is consistent with the results in Lu et al. (2023). The visibility null amplitudes are shorter than the M87 EHT observation in 2017 and 2018 (at around 3.4 Gλ).

The disagreement of visibility null amplitudes between 3 mm and 1 mm implies that the diameter of the observed ring at 3 mm is larger than that of the ring at 1 mm, since the baseline location of the visibility null amplitude scales inversely with the ring diameter (Lu et al. 2023). For comparison, in Sagittarius A*, the intrinsic (descattered) size of the compact horizon-scale source changes by about a factor of 2 between 1 mm and 3 mm emission (Issaoun et al. 2019; Event Horizon Telescope Collaboration 2022). The M87 ring diameter at 86 GHz was estimated in Lu et al. (2023) to be (64−8+4$64_{-8}^{+4}$ μas), based on the circular ring fitting with the optimal set of SMILI images and visibility domain model fitting. Furthermore, they found that the thick ring is preferred over a thin ring using image analysis and model fitting.

In this article, the posterior distribution of the M87 ring diameter, width, and ellipticity is estimated using the variational image domain analysis (VIDA) software (Tiede et al. 2022). To validate the robustness of the ring structure, resolve provides more reliable error estimates than the previous report since the uncertainty is estimated from posterior sample images, and not from the selected top-set of images. VIDA is an image feature extraction tool that treats each image as a probability distribution and compares the image to the geometrical model image (template) by utilizing Kullback-Leiber (KL) or Bhattacharyya (Bh) divergences as the objective function. Using the corresponding template, the ring features of each posterior sample image are estimated. In this analysis, we used the CosineRingwFloor template, and the azimuthal intensity distribution is described by a cosine expansion.

Figure 3 shows that the M87 ring diameter estimated from the resolve images is 60.9 ± 2.2 μas and is 61.0 μas from the DoG-HiT image. The estimated ring diameter using resolve and DoG-HiT is within the errors of the estimation (64−8+4μas$64_{-8}^{+4} \mu \mathrm{as}$) in Lu et al. (2023). The width is 16.0 ± 0.9 μas using resolve and 20.6 μas using DoG-HiT. The discrepancy in the ring width between resolve and DoG-HiT results from the sparse UV coverage beyond the first visibility null amplitude. The visibility domain model fitting (see Fig. S10 in Lu et al. 2023) shows that the UV coverage beyond the first null is crucial to determine the best fitting model. However, the thick ring and thin ring models show similar reduced χ2 due to the limited UV coverage beyond the first null. As a result, the ring width is not well constrained by the data. Furthermore, we note that the ring diameter and width show anti-correlation, which is due to the finite resolution of the telescope array. The effective radius of the ring decreases with a larger ring width (see Appendix G of Event Horizon Telescope Collaboration 2019c), which explains the dependence between the estimated diameter and width. This anti-correlation is also shown by SMILI reconstructions (see Fig. S14 in Lu et al. 2023). The ellipticity of the ring, τ, is defined as τ = 1 − b/a, where a is the semimajor axis lengths and b is the semiminor axis lengths of the elliptical ring. The τ is 0.06 ± 0.04 by resolve and 0.04 by DoG-HiT, which means there is no significant ellipticity of the M87 ring. A detailed description of the VIDA.jl software can be found in Tiede et al. (2022).

The observation by GMVA+ALMA M87 on 2018, April 14 was conducted a week before the EHT M87 observation (2018, April 21 and 25), which is shorter than the expected timescale for decorrelation of the emission pattern (Georgiev et al. 2022). The flux density of the compact region 200 μas × 200 μas) at 1 mm is 0.5 ± 0.1 Jy by DIFMAP, eht-imaging, and SMILI softwares (see Table 2 in Event Horizon Telescope Collaboration 2024). We note that the compact region flux density constraints from EHT 2018 data were challenging due to the lack of short baseline coverage (Event Horizon Telescope Collaboration 2024). The total flux density at 3 mm is 0.57 ± 0.03 Jy on mas scales, and the flux density of the compact region4 at 3 mm is 0.33 ± 0.02 Jy by resolve. DoG-HiT images estimate 0.43 Jy in the compact field of view. As a result, the spectral index α of the M87 compact region is slightly positive α ∼ 0.4. This implies a mixed optical depth in the core, under a caveat that the emitting region at 3 mm is larger than at 1 mm, following the ring diameter analysis. The observed ratio of flux densities is reasonably consistent with the predictions of the numerical models, which typically indicate an inhomogeneous optical depth in the compact region (Event Horizon Telescope Collaboration 2019d, 2021b; Palumbo et al. 2024).

thumbnail Fig. 3

Posterior distribution of the M87 ring diameter, thickness, and ellipticity estimated from 100 resolve posterior sample images. Contours show 1σ and 2σ cumulative regions. The probability density function is obtained using Gaussian kernel density estimation. Each black dot marks the estimated ring parameter from the resolve posterior sample images. The red marks correspond to the estimation of the ring diameter and ring thickness obtained by the DoG-HiT reconstruction (diameter: 61.0 μas, thickness: 20.6 μas, ellipticity: 0.04).

3.2 Extended jet emission of M87

Our reconstructions of the limb-brightened M87 jet structure broadly agree with the one described in Lu et al. (2023): we see an edge-brightened jet anchored to the vicinity of the ring-like feature. One peculiar feature in the image reported by Lu et al. (2023) is the presence of a bright spine along the jet axis. This central ridge line may be related to the triple-helix structure in the jet of M87 that has been observed at larger scales (Nikonov et al. 2023).

However, it is questionable whether this feature of the image represents a real structure on-sky, or appears as a consequence of imaging artifacts. Particularly, it has been recently demonstrated that CLEAN deconvolution errors are prone to produce inner ridge lines for edge-brightened jet configurations (Pashchenko et al. 2023). In fact, the central spine is less prominent in the DoG-HiT and resolve reconstructions.

DoG-HiT and resolve allow the robustness of the central spine to be quantified from two independent perspectives: by calibration-independent imaging with a minimal human bias in DoG-HiT, and by uncertainty estimation in resolve. For DoG-HiT, the small scale structure and the large scale (diffuse) structure are represented by different wavelets, which are ultimately expressed by the multiresolution support. This allows us to estimate the robustness of the central ridge by a jackknife test, that is, we cut the diffuse emission at the location of the central spine and recalculate the fit statistics to the (calibration-independent) closure quantities. To this end, we applied the following strategy. We flagged long baselines of more than 2 Gλ (to focus the analysis on the diffuse emission), then we fitted the closure phases and closure amplitudes with DoG-HiT, only varying coefficients in the multiresolution support, and calculated the updated fit statistics. Next, we masked out the diffuse, central spine from the multiresolution support, and refit the closure quantities. Finally, we compared the scoring with a central spine and without a central spine.

This strategy resembles a standard strategy in VLBI, which is often applied in the discussion of the existence of a counter-jet. However, we note some key advantages of the strategy applied by us, compared to CLEAN. First, DoG-HiT directly fits closure quantities, hence the conclusions that we can draw are less dependent on the phase and amplitude self-calibration. Second, we perform the jackknife test on a multiscalar domain, which allows us to divide the emission more clearly into small and large scale structures. Finally, DoG-HiT directly fits a model to the data that is composed of diffuse and compact emission and does not need a final convolution with the beam, as opposed to the point source model of CLEAN (e.g., see the discussion in Müller & Lobanov 2023a). Hence, we compare the fit quality of the approximated on-sky representation rather than an unphysical list of CLEAN components (which would question the interpretability of the χ2-statistics).

We obtain χcph 2=1.104,χcla 2=1.494$\chi_{\text {cph}}^{2}=1.104, \chi_{\text {cla}}^{2}=1.494$ when fitting the data with a central spine, and χcph2=1.061,χcla2=1.592$\chi_{\mathrm{cph}}^{2}=1.061, \chi_{\text{cla}}^{2}=1.592$ when fitting without a central spine. Neither mode is strongly favored. From this study, we cannot report the conclusive detection of a central spine in the image.

An alternative perspective on the robustness of image features is offered by the built-in uncertainty quantification in resolve. In Fig. 4, the transverse flux intensity profiles of the M87 jet emission are depicted. The mean and standard deviation of the transverse jet profile can be obtained from 100 posterior sample images by resolve. The intensity profile of the jet at 0.25, 0.5, 0.75 mas shows that the edge-brightened structure (two peaks) is prominent, however a significant central spine structure is not seen. In Fig. 4, the standard deviation of the pixel fluxes at the central spine is not particularly higher than the limb-brightened feature. In Fig. 1, the DoG-HiT image shows a central spine, which however is fainter than that of the CLEAN reconstruction. Therefore, the central spine in the images obtained by CLEAN in Lu et al. (2023) may be a consequence of CLEAN artifacts resulting from the CLEAN windows and sparsity-promoting CLEAN sky model. Further observations with additional short baseline antennas are required to conclude the detection of the central spine. The edge-brightened morphology that we recover is both consistent with the earlier observations at 86 GHz (Kim et al. 2018) and well-motivated theoretically (e.g., Yang et al. 2024). The counter jet is not detected consistently in the three images, it is therefore not discussed further in this work.

thumbnail Fig. 4

Galaxy M87 jet transverse profiles and intensity map obtained by resolve. Intensity in the map is represented by false color according to the color bar used in Fig. 1 in logarithmic scale. The image is rotated 18∘ clockwise. Intensity plots at the bottom of the figure show flux density profiles of the jet at 0.25, 0.5 and 0.75 mas from the phase center. Each vertical line corresponds to a location where profiles were extracted.

3.3 Substructure of the M87 ring

A feature that appears consistent across all four different imaging algorithms (SMILI, CLEAN, resolve, and DoG-HiT) are the two bright blobs in the ring. While images are recovered with various resolutions, the recovered ring emission always shows the double structure within the ring toward the top and the bottom, at consistent positions across different imaging methods. We note that the limb-brightened structures that are recovered by resolve close to the ring seem to connect to the ring exactly at the positions where the brighter blobs within the ring occur consistently for all four reconstructions. Hence, it may be natural to interpret the double pattern in the ring as a physical phenomenon. In this subsection we present some discussion on how real these features may be.

First, it is noticeable that these brighter regions in the ring form a double structure that is point-symmetric to the center of the ring. That supports the interpretation of these structures as an imaging artifact, especially due to the sparse coverage at long baselines. In fact, Lu et al. (2023) have tested SMILI and CLEAN reconstructions on synthetic data and found that similar structures are artificially introduced by the imaging procedure (compare Sect. 4 and particularly Fig. S8 in Lu et al. 2023).

Multiple well-understood artifacts may cause such a double structure. Here, we discuss the three most natural scenarios. First, it could be caused by specific choices of the regularization assumption inherent to the respective imaging algorithm. Second, the structure could be introduced artificially by residual gain effects. Finally, the structure may be described by a residual sidelobe structure, which would make it essentially a consequence of the sparsity of the UV coverage. In what follows, we will discuss each of these concerns individually.

Four imaging methods that were utilized here approach the image reconstruction from four vastly different perspectives: CLEAN recovers the structure in an inverse modeling framework essentially processing a sparsity-promoting regularization approach (sky brightness distribution is represented by a collection of delta components) (Lannes et al. 1997), SMILI approaches the image by a weighted sum of multiple handcrafted data and regularization terms in a RML framework (Akiyama et al. 2017); DoG-HiT processes multiscale functions in the context of compressive sensing (Müller & Lobanov 2022); and resolve estimates the posterior distribution of the image and gains from the prior model encoding source and instrument information and likelihood (i.e., the data; Kim et al. 2024). We note that the prior in Bayesian imaging can be interpreted as regularizers in RML methods, and vice versa (Kim et al. 2024). The multiscale functions in DoG-HiT and the Gaussian process prior in resolve are flexible and do not ask for the double structure in the ring as prior knowledge. The fact that all four independent methods that use a variety of prior information (regularization) lead to a similar structure challenges the interpretation of the double structure as an artifact from the assumptions and prior information applied by the imaging procedure.

The structure may be a consequence of unsolved gain residuals. In fact, it is a possible issue that the alternating self-calibration and cleaning procedure produce “phantom” structures point symmetric to the origin that has been reported in practice for a long time. We note, however, that the double feature also appears in the DoG-HiT reconstruction that are independent from gain corruptions (closure-only imaging), which makes a potential cause by the calibration of the phases less likely. This claim is supported by the resolve reconstruction which solves for the self-calibration with imaging simultaneously in a probabilistic setup.

Finally, we ask whether the double structure could be introduced by the sparsity of UV coverage. That is a possibility that can never be eradicated completely, simply because the observation fundamentally misses relevant visibilities through the gaps in the UV coverage. Consequently, there is missing information. In fact, the Fourier domain representation of a double source is a fringe pattern (compare Fig. 5), which is exactly the kind of artifact that not fully cleaned residuals may introduce. This issue has been identified through synthetic data tests by Lu et al. (2023). They generated synthetic data sets with artificial baselines from symmetric ring ground truth images to understand the importance of the long west-east baselines. CLEAN and SMILI reconstructions from the synthetic data with additional baselines were able to reconstruct the ground truth. From the synthetic data with original UV coverage, SMILI reconstructions could recover symmetric ring ground truth images reasonably well, but CLEAN reconstructions tended to generate asymmetric ring structures. The occurrence of the double structure may be related to the lack of long west-east baselines. Nevertheless, we argue that the recovered double structure is represented by the measured visibilities. In Fig. 5, we show the central ring image (top panels) and the amplitudes of the full Fourier transform of the reconstructions (middle panels) for DoG-HiT (left panels) and resolve (right panels). A ring feature is identified in the Fourier domain by a first zero that is clearly covered by observations (compare, e.g., to the model fitting discussions in Lu et al. 2023). Moreover, the visibility domain representations of DoG-HiT and resolve show the ellipticity of the null visibility amplitude points in Fig. 5 (middle panel). We note that the location of the null visibility points in the UV domain is elliptical, with an elongation in the direction of the jet. This may result in bright spots in the image aligned perpendicularly to the jet, related to the stretch of the image domain ring. The data from multiple antennas (EF, KP, OV, PV, and YS) at the longest west-east baselines imply that the ellipticity of the null visibility points originates from the data.

To analyze the UV-coverage sparsity that corresponds to the M87 ring substructure, we subtract a uniform ring feature from the images (defined by 60% of the respective emission peak) to extract the double feature on top of the ring. The resulting double patterns are depicted in the bottom panels of Fig. 5, and their respective amplitudes in the Fourier domain in the last row. The observed UV points are overplotted with red crosses. Observations span the main fringe and the first side-lobe uniformly, the fringes are not produced exclusively in the gaps in the UV coverage. These findings, as well as the striking similarity between multiple imaging approaches, constitute a convincing argument for the physical nature of bright spots along the ring. However, a definitive answer cannot be given with the quality of the existing data set, particularly in the presence of the synthetic data tests performed by Lu et al. (2023), and follow-up observations are needed.

If the emission that forms the ring image is dominated by the accretion disk, elongation in the direction perpendicular to the ring may be a simple geometric effect for the inclined observer. While the 230 GHz ring image is strongly asymmetric with respect to the jet axis (Event Horizon Telescope Collaboration 2024), our 86 GHz reconstructions exhibit a high degree of symmetry. In numerical general relativistic magnetohydrodynamic (GRMHD) simulations of accretion within the EHT library (Dhruv et al. 2025), consistent behavior appears for some retrograde accretion models (negative black hole spins). In that case, 230 GHz image asymmetry is driven primarily by the spin effects (Event Horizon Telescope Collaboration 2019d), which are dominant very near the event horizon. In the 86 GHz image, which is formed by a more extended emission, these effects are balanced by the Doppler boost that enhances brightness on the opposite side of the black hole, which results in a ring structure elongated perpendicularly to the jet axis, with bright spots on both sides of the central brightness depression.

thumbnail Fig. 5

Reconstructions of the central feature by DoG-HiT (left panels) and resolve (right panels). Top panels: reconstructions of the compact emission. Top-middle panels: the central ring feature. Middle panels: The recovered amplitudes of the central feature, with the UV coverage overplotted (red points). Middle-bottom panels: the central ring feature, cut at 60% of the respective peak brightness, which showcases the double blob pattern on top of the ring. Bottom panels: the amplitudes and the UV coverage of the double pattern alone.

4 Conclusions

Lu et al. (2023) present observations of the core and jet in M87 observed with the GMVA+ALMA at 86 GHz. The image contains a ring-like feature that looks similar to the one reported by the EHT (Event Horizon Telescope Collaboration 2019a, 2024), but with an approximately 1.5 times larger ring diameter resulting from the null visibility points and phase jump at 2.3 Gλ. For the first time, this ring feature is connected to an innermost jet in the same image, possibly providing constraints on the launching mechanism of the jet. Furthermore, the recovered image contains several fainter features that may be of great importance for a scientific interpretation, especially in the context of recent works on the large scale jet structure in M87 (Kim et al. 2018; Cui et al. 2023; Nikonov et al. 2023; Kim et al. 2023).

To validate the results reported in Lu et al. (2023), we apply two more imaging algorithms specially designed to study the robustness of the recovered features: by Bayesian self-calibration and imaging with resolve, and by closure-only imaging with DoG-HiT. The distinctive features of resolve and Dog-HiT, namely the probabilistic approach and the multiscale wavelet-based deconvolution algorithm, allow us to quantify the robustness of the recovered M87 ring and extended jet emission. We obtained the posterior distribution of the ring diameter, width, and ellipticity by analysis of the resolve posterior sample images. We confirm the M87 ring-like structure at 86 GHz with a diameter of 60.9 ± 2.2 μas, a thickness of 16.0 ± 0.9 μas, and an ellipticity of 0.06 ± 0.04 by resolve, and a diameter of 61.0 μas, a thickness of 20.6 μas, and an ellipticity of 0.04 by DoG-HiT. The estimated ring diameter is consistent with the estimate (64−8+4μas$64_{-8}^{+4} \mu \mathrm{as}$) in Lu et al. (2023).

Furthermore, image reconstructions using resolve and DoG-HiT show that the ring is embedded in a strongly edge-brightened large scale jet structure, which agrees with the findings reported in Lu et al. (2023). Among the upper and lower arm of the edge-brightened jet, the CLEAN reconstruction presented in Lu et al. (2023) features a third, central spine that may be interpreted together with a triple-helix structure at larger scales (Nikonov et al. 2023). However, the central spine structure in the resolve and DoG-HiT reconstructions is less prominent compared to the CLEAN reconstruction. Our analysis shows that this central spine is neither necessary to fit the data, nor supported by the uncertainty quantification by resolve in the image domain. The validation by two independent imaging methods implies that the central spine is fainter than previously reported.

Finally, all utilized imaging algorithms coincide on the same substructure in the ring consisting of two bright spots to the north and the south, co-located with inner anchor points on the edge-brightened jet to the horizon scale. Stemming from a coincidence of this phenomenon across a variety of imaging algorithms, and its representation in the Fourier domain, we argue that is more likely that this substructure of the ring results from the data, and not an imaging or self-calibration artifact, although an artifact that originates from sparse UV coverage is not ruled out. The potential physical origin of these structures is unclear. If real, they might be transient emission structures during the observational period, however, their alignment perpendicular to the jet suggests otherwise. They might also be permanent structures potentially linked to the disk-jet transition, as their location within the disk seems to coincide with the position the jet edges point to. Finally, they may be a result of an interplay of Doppler boost and black hole spin effects, as seen in numerical simulations of retrograde accretion. Investigating the size and asymmetry of the 86 GHz emission ring in the EHT M87 simulation library could be a fruitful path toward understanding the detailed physics responsible for the morphology of our image.

Acknowledgements

We thank Jae-Young Kim for helpful suggestions and feedback on drafts of the manuscript, Jack Livingston and Thomas Krichbaum for comments and informative discussions. J. K. and A. N. received financial support for this research from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. This work was supported by the M2FINDERS project funded by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant Agreement No. 101018682). R.-S.L. is supported by the National Science Fund for Distinguished Young Scholars of China (Grant No. 12325302) and the Shanghai Pilot Program for Basic Research, CAS, Shanghai Branch (JCYJ-SHFY-2021-013). M.W. is supported by a Ramón y Cajal grant RYC2023-042988-I from the Spanish Ministry of Science and Innovation. This research has made use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM, Onsala, Metsahovi, Yebes, the Korean VLBI Network, the Greenland Telescope, the Green Bank Observatory and the Very Long Baseline Array (VLBA). The VLBA and the GBT are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Greenland Telescope (GLT) is operated by the Academia Sinica Institute of Astronomy and Astrophysics (ASIAA) and the Smithsonian Astrophysical Observatory (SAO). The data were correlated at the correlator of the MPIfR in Bonn, Germany.

Appendix A Uncertainty estimation by resolve

Figure A.1 shows the pixel-wise relative uncertainty of the resolve image in Fig. 1. The relative uncertainty is defined as the sky brightness posterior standard deviation normalized by the posterior mean from 100 posterior sample images. Lower relative uncertainty values of the M87 ring emission and limbbrightened are estimated, that means they are well-constrained by the data. However, the relative uncertainty values in the central spine are higher compared to the limb-brightened jet emission. Figure A.2 depicts the resolve posterior amplitude gains. ALMA (AA) amplitude gain solutions show smooth behavior with lower posterior standard deviations compared to other arrays due to the high sensitivity of ALMA array. The gain amplitudes with high uncertainties result from the absence of the data (BR: <4h, EF: >6h, GB: <4h, KP: <3h, LA: <3h, PT: <3h, and YS: RCP gains). For instance, the RCP gain amplitudes for YS antenna have uniformly distributed high uncertainties since only single polarization mode (LCP) was observed for YS antenna.

thumbnail Fig. A.1

resolve image pixel-wise relative uncertainty, which is the sky brightness posterior standard deviation normalized by the posterior mean from the resolve reconstruction in the top panel of Fig. 1.
thumbnail Fig. A.2

resolve posterior amplitude gains. The gain as a function of time is illustrated as a thin line with a semi-transparent standard deviation. The left and right columns of the figure show gains from the right (RCP) and left (LCP) circular polarizations correspondingly. Each row represents an individual antenna, whose abbreviated name is indicated in the bottom left corner of each LCP plot.

Appendix B Contour plots

thumbnail Fig. B.1

Same as Fig. 2, but overlaid with successive contours that increase by a factor of 2$\sqrt{2}$, starting from 20% of the peak.
thumbnail Fig. B.2

Same as Fig. 1, but overlaid with successive contours that increase by a factor of 2, starting from 0.2% of the peak. The starting level was chosen based on the structure of the CLEAN image.

References

  1. Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017, ApJ, 838, 1 [NASA ADS] [CrossRef] [Google Scholar]
  2. Arras, P., Bester, H. L., Perley, R. A., et al. 2021, A&A, 646, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Arras, P., Frank, P., Haim, P., et al. 2022, Nat. Astron., 6, 259 [NASA ADS] [CrossRef] [Google Scholar]
  4. Asada, K., Nakamura, M., & Pu, H.-Y. 2016, ApJ, 833, 56 [NASA ADS] [CrossRef] [Google Scholar]
  5. Blackburn, L., Pesce, D. W., Johnson, M. D., et al. 2020, ApJ, 894, 31 [Google Scholar]
  6. Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. 2016, arXiv e-prints, [arXiv:1601.00670] [Google Scholar]
  7. Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 413, 813 [Google Scholar]
  8. Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11 [Google Scholar]
  9. Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23 [Google Scholar]
  10. Clark, B. G. 1980, A&A, 89, 377 [NASA ADS] [Google Scholar]
  11. Cui, Y., Hada, K., Kawashima, T., et al. 2023, Nature, 621, 711 [NASA ADS] [CrossRef] [Google Scholar]
  12. Dhruv, V., Prather, B., Wong, G., & Gammie, C. F. 2025, ApJS, 277, 16 [Google Scholar]
  13. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019a, ApJ, 875, L1 [Google Scholar]
  14. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875 L6 [Google Scholar]
  15. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019c, ApJ, 875, L4 [Google Scholar]
  16. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019d, ApJ, 875, L5 [Google Scholar]
  17. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021a, ApJ, 910, L12 [Google Scholar]
  18. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021b, ApJ, 910, L13 [Google Scholar]
  19. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022, ApJ, 930, L12 [NASA ADS] [CrossRef] [Google Scholar]
  20. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2023, ApJ, 957, L20 [NASA ADS] [CrossRef] [Google Scholar]
  21. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2024, A&A, 681, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Frank, P., Leike, R., & Enßlin, T. A. 2021, Entropy, 23, 853 [NASA ADS] [CrossRef] [Google Scholar]
  23. Georgiev, B., Pesce, D. W., Broderick, A. E., et al. 2022, ApJ, 930, L20 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hada, K. 2017, Galaxies, 5, 2 [Google Scholar]
  25. Högbom, J. A. 1974, A&AS, 15, 417 [Google Scholar]
  26. Issaoun, S., Johnson, M. D., Blackburn, L., et al. 2019, ApJ, 871, 30 [Google Scholar]
  27. Junklewitz, H., Bell, M. R., Selig, M., & Enßlin, T. A. 2016, A&A, 586, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018, A&A, 616, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Kim, J.-Y., Savolainen, T., Voitsik, P., et al. 2023, ApJ, 952, 34 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kim, J.-S., Nikonov, A. S., Roth, J., et al. 2024, A&A, 690, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Knollmüller, J., & Enßlin, T. A. 2019, arXiv e-prints [arXiv:1901.11033] [Google Scholar]
  32. Lannes, A., Anterrieu, E., & Marechal, P. 1997, A&AS, 123, 183 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Lu, R.-S., Asada, K., Krichbaum, T. P., et al. 2023, Nature, 616, 686 [CrossRef] [Google Scholar]
  34. Müller, H., & Lobanov, A. P. 2022, A&A, 666, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Müller, H., & Lobanov, A. P. 2023a, A&A, 673, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Müller, H., & Lobanov, A. P. 2023b, A&A, 672, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Müller, H., Mus, A., & Lobanov, A. 2023, A&A, 675, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Müller, H., Massa, P., Mus, A., Kim, J.-S., & Perracchione, E. 2024, A&A, 684 A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Mus, A., Müller, H., & Lobanov, A. 2024a, A&A, 688, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Mus, A., Müller, H., Martí-Vidal, I., & Lobanov, A. 2024b, A&A, 684, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Nikonov, A. S., Kovalev, Y. Y., Kravchenko, E. V., Pashchenko, I. N., & Lobanov, A. P. 2023, MNRAS, 526, 5949 [NASA ADS] [CrossRef] [Google Scholar]
  42. Okino, H., Akiyama, K., Asada, K., et al. 2022, ApJ, 940, 65 [NASA ADS] [CrossRef] [Google Scholar]
  43. Palumbo, D. C. M., Bauböck, M., & Gammie, C. F. 2024, ApJ, 970, 151 [Google Scholar]
  44. Pashchenko, I. N., Kravchenko, E. V., Nokhrina, E. E., & Nikonov, A. S. 2023, MNRAS, 523, 1247 [NASA ADS] [CrossRef] [Google Scholar]
  45. Roth, J., Arras, P., Reinecke, M., et al. 2023, A&A, 678, A177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Thyagarajan, N., Nityananda, R., & Samuel, J. 2022, Phys. Rev. D, 105, 043019 [NASA ADS] [CrossRef] [Google Scholar]
  47. Tiede, P., Broderick, A. E., & Palumbo, D. C. M. 2022, ApJ, 925, 122 [CrossRef] [Google Scholar]
  48. Twiss, R. Q., Carter, A. W. L., & Little, A. G. 1960, The Observatory, 80, 153 [NASA ADS] [Google Scholar]
  49. Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128 [Google Scholar]
  50. Yang, H., Yuan, F., Li, H., et al. 2024, Sci. Adv., 10, eadn3544 [NASA ADS] [CrossRef] [Google Scholar]
4

Defined as the 200 μas (or about 50 Rg) field of view centered on the ring.

All Figures

thumbnail Fig. 1

Image reconstructions of GMVA+ALMA M87 at 86 GHz. Each image presents results obtained by a different algorithm. The top posterior mean image was reconstructed using the resolve Bayesian self-calibration and imaging method. The middle image was obtained using DoG-HiT closure amplitudes and closure phases only imaging. The bottom image was obtained using the CLEAN self-calibration and imaging in Lu et al. (2023). The CLEAN image is convolved with an elliptical beam, which is represented as an ellipse with sizes 79 × 37 μas, PA = −63∘ in the bottom-left corner.
In the text
thumbnail Fig. 2

Ring of the M87 image reconstructions at 86 GHz from Fig. 1. The color shows intensity in Jy mas−2 according to the linear color bar located at the right of the figure. Each image presents results obtained by a different algorithm, whose names are indicated in the lower right corners. The left image is the resolve posterior mean image using the Bayesian self-calibration and imaging method. The left-center image is the DoG-HiT reconstruction using closures only imaging. The right-center image represents the CLEAN reconstruction with the over-resolved 37 μas circular beam in Lu et al. (2023). The right image shows the SMILI reconstruction in Lu et al. (2023). All images were processed by a Gaussian interpolation.
In the text
thumbnail Fig. 3

Posterior distribution of the M87 ring diameter, thickness, and ellipticity estimated from 100 resolve posterior sample images. Contours show 1σ and 2σ cumulative regions. The probability density function is obtained using Gaussian kernel density estimation. Each black dot marks the estimated ring parameter from the resolve posterior sample images. The red marks correspond to the estimation of the ring diameter and ring thickness obtained by the DoG-HiT reconstruction (diameter: 61.0 μas, thickness: 20.6 μas, ellipticity: 0.04).
In the text
thumbnail Fig. 4

Galaxy M87 jet transverse profiles and intensity map obtained by resolve. Intensity in the map is represented by false color according to the color bar used in Fig. 1 in logarithmic scale. The image is rotated 18∘ clockwise. Intensity plots at the bottom of the figure show flux density profiles of the jet at 0.25, 0.5 and 0.75 mas from the phase center. Each vertical line corresponds to a location where profiles were extracted.
In the text
thumbnail Fig. 5

Reconstructions of the central feature by DoG-HiT (left panels) and resolve (right panels). Top panels: reconstructions of the compact emission. Top-middle panels: the central ring feature. Middle panels: The recovered amplitudes of the central feature, with the UV coverage overplotted (red points). Middle-bottom panels: the central ring feature, cut at 60% of the respective peak brightness, which showcases the double blob pattern on top of the ring. Bottom panels: the amplitudes and the UV coverage of the double pattern alone.
In the text
thumbnail Fig. A.1

resolve image pixel-wise relative uncertainty, which is the sky brightness posterior standard deviation normalized by the posterior mean from the resolve reconstruction in the top panel of Fig. 1.
In the text
thumbnail Fig. A.2

resolve posterior amplitude gains. The gain as a function of time is illustrated as a thin line with a semi-transparent standard deviation. The left and right columns of the figure show gains from the right (RCP) and left (LCP) circular polarizations correspondingly. Each row represents an individual antenna, whose abbreviated name is indicated in the bottom left corner of each LCP plot.
In the text
thumbnail Fig. B.1

Same as Fig. 2, but overlaid with successive contours that increase by a factor of 2$\sqrt{2}$, starting from 20% of the peak.
In the text
thumbnail Fig. B.2

Same as Fig. 1, but overlaid with successive contours that increase by a factor of 2, starting from 0.2% of the peak. The starting level was chosen based on the structure of the CLEAN image.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.