The galaxy M87 is one of the prime targets for high resolution radio imaging pursuing the ringlike shadow of its supermassive black hole the innermost regions of accretion flow and the formation of the relativistic jet. However it remains challenging to observe both jointly. Only recently global mmVLBI array GMVAALMA observations at 86 GHz in 2018 were able to reconstruct the M87 black hole shadow and the extended jet emission simultaneously. In order to analyze the ring and jet of M87 conventional CLEAN algorithms were mainly employed alongside the RML method SMILI in the previous work. To test the robustness of the reconstructed structures of M87 GMVAALMA observations at 86GHz we estimate the ring diameter width and the extended jet emission with the possible central spine by two different novel imaging algorithms resolve and DoGHiT. Overall reconstructions are consistent with the results reported in the previous paper. The ring structure of the M87 is resolved at higher resolution and the posterior distribution of M87 ring features is explored. The resolve images show that the ring diameter is 60.9  2.2 muas and width is 16.0  0.9 muas. The ring diameter is 61.0 muas and width is 20.6 muas by DoGHiT. The ring diameter is therefore in agreement with the estimation 6448 muas by SMILI and the geometrical modeling. Two bright spots in the ring are reconstructed by four independent imaging methods the substructure in the ring is therefore most likely originated from the data. A consistent limbbrightened jet structure is reconstructed by resolve and DoGHiT albeit with a less pronounced central spine. Modern datadriven imaging methods confirm the ring and jet structure in M87 complementing traditional VLBI methods with novel perspectives on the significance of recovered features. They confirm the result of the previous report.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Imaging the black hole shadow and extended jet of M87
{'Imaging the black hole shadow and extended jet of M87': "Jong-Seo Kim 1, ⋆ , Hendrik Müller 1 , 2 , Aleksei S. Nikonov 1 , Ru-Sen Lu 1 , 3 , 4 , Jakob Knollmüller 5 , Torsten A. Enßlin 6 , 7 , Maciek Wielgus 1 , and Andrei P. Lobanov 1 , 8  \n- 1 Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany\n- 2 Jansky Fellow of National Radio Astronomy Observatory, 1011 Lopezville Rd, Socorro, NM 87801, USA\n- 3 Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai, People's Republic of China\n- 4 Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, Nanjing, People's Republic of China\n- 5 Radboud University, Heyendaalseweg 135, 6525 AJ, Nijmegen, The Netherlands\n- 6 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany\n- 7 Ludwig-Maximilians-Universität, Geschwister-Scholl-Platz 1, 80539 Munich, Germany\n- 8 Institut für Experimentalphysik, Universität Hamburg, Luruper Chaussee 149, 22761, Hamburg, Germany  \nReceived August 29, 2024", 'ABSTRACT': 'Context. The galaxy M87 is one of the prime targets for high resolution radio imaging pursuing the ring-like "shadow" of its supermassive black hole, the innermost regions of accretion flow, and the formation of the relativistic jet. However, despite the physical process of ring-like structure and jet are closely interconnected, and this connection may unravel the physics behind the jet-launching mechanism, it remains challenging to observe both jointly. Only recently, global mm-VLBI array (GMVA) + ALMA observations at 86 GHz in 2018 were able to reconstruct the M87 black hole shadow and the extended jet emission simultaneously. In order to analyze the ring and jet of M87, conventional CLEAN algorithms were mainly employed alongside the regularized maximum likelihood method SMILI in the previous work.  \nAims. To test the robustness of the reconstructed structures of M87 GMVA + ALMA observations at 86GHz, we estimate the ring diameter, width, and the extended jet emission with the possible central spine by two di ff erent novel imaging algorithms: resolve and DoG-HiT .  \nMethods. We performed Bayesian self-calibration and imaging with uncertainty estimation by resolve . We reconstructed the image using closure amplitudes and phases only by DoG-HiT .  \nResults. Overall reconstructions are consistent with the results reported in the previous paper. The ring structure of the M87 is resolved at higher resolution and the posterior distribution of M87 ring features is explored. The resolve images show that the ring diameter is 60 . 9 ± 2 . 2 µ as and width is 16 . 0 ± 0 . 9 µ as. The ring diameter is 61 . 0 µ as and width is 20 . 6 µ as by DoG-HiT . The ring diameter is therefore in agreement with the estimation (64 + 4 -8 µ as) by SMILI and the geometrical modeling. Two bright spots in the ring are reconstructed by four independent imaging methods, the substructure in the ring is therefore most likely originated from the data. A consistent limb-brightened jet structure is reconstructed by resolve and DoG-HiT , albeit with a less pronounced central spine.  \nConclusions. Modern data-driven imaging methods confirm the ring and jet structure in M87, complementing traditional VLBI methods with novel perspectives on the significance of recovered features. They confirm the result of the previous report.  \nKey words. techniques: interferometric - techniques: image processing - techniques: high angular resolution - methods: statistical galaxies: active - galaxies: jets - galaxies: individual (M87)', '1. Introduction': 'Continued improvements of sensitivity and uv-coverage of very long baseline interferometry (VLBI) allow to image and analyze radio emission from the close vicinity of supermassive black holes (SMBH) and inside the formation zone of relativistic jet in active galactic nuclei (AGNs). The recent Event Horizon Telescope (EHT) observations have captured the black hole "shadow" in M87* and SgrA* (Event Horizon Telescope Collaboration et al. 2019a, 2022), and the Global mm-VLBI Array (GMVA) and Very Long Baseline Array (VLBA) observations have 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 × 10 9 M ⊙ (Event Horizon Telescope Collaboration et al. 2019b) and redshift z = 0 . 004283 (Cappellari et al. 2011) pro-  \nes the 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 R g (where R g = GM / c 2 ) per milliarcsecond. So far, the black hole shadow and relativistic jet of M87 at radio frequencies are observed independently due to instrumental limitations and characteristics of the source in di ff erent frequency regimes.  \nOnly lately, observations of M87 with the GMVA and the phased Atacama Large Millimeter / Submillimeter Array (ALMA) in 2018 was able to detect the black hole shadow and extended jet emission altogether. 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 called SMILI (Akiyama et al. 2017). The global Millimeter VLBI Array (GMVA) observations at 86 GHz play a significant role to  \nresolve the core and extended jet emission of M87. Prior to EHT M87 observation in 2017 (Event Horizon Telescope Collaboration et al. 2019a), GMVA M87 observation in 2014 and 2015 (Kim et al. 2018) provided the core and edge-brightened jet of M87 image down to 13 R g . In 2018, GMVA M87 observation was conducted jointly with the phased ALMA. The joint observation with ALMA facilitates the detection of the ring with extended jet emission since ALMA provides the longest northsouth baselines with improved sensitivity. These findings, for the first time, connected the central ring-like phenomenon which is interpreted as the black hole shadow with the innermost jet. This work is accompanied with 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 et al. 2019a, 2021, 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 connected horizon scales to jet scales.  \nHowever, the data reduction process of GMVA + ALMA is challenging due to the Fourier coverage sparsity of uv-coverage, tropospheric phase corruption at mm wavelength, low signal-tonoise ratio (SNR), and inhomogeneous antenna statistics due to the di ff erent sensitivity of antennas and characteristics of polarization receivers between ALMA (linear) and the rest of the antennas (circular). This di ffi culty jeopardizes the interpretation of some of the features in the reconstruction, namely the substructure in the ring (two brighter spots) and the jet (inner ridge line). An independent assessment of the robustness of these features is needed, to facilitate the much awaited scientific interpretation of these features.  \nThe latest advancement of forward modeling imaging algorithms enables us to generate more robust results from sparse mm-VLBI data sets. As an example, Bayesian imaging is a probabilistic approach reconstructing the posterior distribution using the Bayes\' theorem. Hence, Bayesian imaging is able to estimate uncertainty of parameters, such as image features and instrumental gains, however the image reconstruction is comparatively computationally demanding. RML methods reconstruct an image by minimizing the data fidelity term with regularizers. The forward modeling approaches fit the model to the data directly in visibility domain and the model in Bayesian imaging and RML methods tend to be flexible and can be modified easily. 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 traditional inverse modeling CLEAN algorithm and we can reconstruct reproducible images with improved resolution in less supervised fashion. A detailed comparison between CLEAN and forward modeling approaches can be found in Arras et al. (2021); Müller et al. (2024). In this work, we reconstruct images by two imaging algorithms: we perform Bayesian self-calibration and imaging jointly by the Bayesian imaging software resolve (Junklewitz et al. 2016; Arras et al. 2022; Roth et al. 2023; Kim et al. 2024) and reconstruct an image with closure amplitudes and closure phases only (Chael et al. 2018) using RML-based DoG-HiT software (Müller & Lobanov 2022, 2023a,b). Those 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 robust-  \nss of the recovered features from two alternative, supplementary perspectives, by self-calibration and imaging in a probabilistic point of view and by imaging only with closure quantities without potential self-calibration bias.  \nThis manuscript is structured as follows. In section 2, we explain resolve and DoG-HiT image reconstruction methods. In section 3, we show two image reconstruction results and then analyze the robustness of the M87 ring structure and jet emission. In section 4, we summarize our results.', '2.1. Bayesian self-calibration and imaging by resolve': "resolve 1 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 which are consistent with the data are reconstructed by Bayes' theorem in variational inference sense (Blei et al. 2016; Knollmüller & Enßlin 2019; Frank et al. 2021). In this paper, we used 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 a ff ordable computational resources. The probabilistic approach can be advantageous for the M87 GMVA + ALMA observation due to the sparse uv-coverage, large and heterogeneous data uncertainties. For details of the M87 GMVA + ALMA data, we refer to the section 1 and 2 of the supplementary information in Lu et al. (2023).  \nWe reconstructed the resolve image in Figure 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 are 100 and the χ 2 n value of the final result is 1 . 1. The wall-clock time for the resolve reconstruction is 5 . 5 hours on a single node of MPIfR cluster with 25 MPI (Message Passing Interface) tasks.  \nSince GMVA + ALMAarray is highly inhomogeneous, it is a reasonable assumption that each antenna gain has di ff erent temporal correlation structure. In resolve , we utilized Gaussian process prior with non-parametric correlation kernel in NIFTy software 2 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 gain solution interval constraints (see Figure A.2). The amplitude gain prior is assumed to be correlated and di ff erent correlation kernels are inferred per antenna. Gain amplitude for right-hand circular polarization (RCP) and left-hand circular polarization (LCP) are assumed to have the same correlation structure to stabilize the selfcalibration and image reconstruction. The phase gain is assumed to be uncorrelated since the phase coherence time in GMVA observation is comparable with data averaging time (10 seconds) and it can be even shorter under poor weather conditions.  \nFurthermore, the posterior distribution of desired parameters, such as each pixel in the image and antenna-based gain solutions can be estimated 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  \nby estimated uncertainties from posterior samples. As an example, if one antennna is problematic, then it would result in 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 Bayesian self-calibration and imaging method for VLBI data and validation with synthetic data can be found in Kim et al. (2024).", '2.2. Imaging with closure quantities by DoG-HiT': 'DoG-HiT is a Regularized Maximum Likelihood (RML) imaging algorithm that models the image by multiscalar wavelet basis functions (Müller & Lobanov 2022). The basis functions are fitted to the uv-coverage, o ff ering a neat separation between covered and non-covered Fourier coe ffi cients, i.e. gaps in the uvcoverage (for more details on the wavelets, we refer to Müller & Lobanov 2023a). The image is recovered by a sparsity promoting forward-backward splitting framework which e ff ectively calculates the multiresolution support, i.e. the set of all statistically significant wavelet scales to represent the image. It has been demonstrated that the multiresolution support is a beneficial prior information that allows for the reconstruction even for sparse and weakly constrained settings (Müller & Lobanov 2022, 2023b).  \nIn this work, we aim to validate the results presented in Lu et al. (2023) by a self-calibration independent technique, i.e. by estimating the robustness of the estimate against the gaincalibration. Therefore, we perform closure-only imaging as pioneered by Chael et al. (2018). We represent the DoG-HiT image (see Figure 1) in a square field of view of 4096 µ as by 512 × 512 pixels. The reconstruction was run on a standard notebook for roughly thirty minutes. The χ 2 n to the closure phases was 1 . 36 and to the closure amplitudes 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, di ff use jet emission. Then we use this image as an initial guess and add all small-scale wavelets relevant to represent the central ring, and minimize the χ 2 -metric to the closure phases and closure amplitudes.  \nIn 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 / quadrilaterals (Twiss et al. 1960; Blackburn et al. 2020; Thyagarajan et al. 2022). Since the number of statistically independent closure phases / amplitudes is therefore smaller than the number of independent visibilities, e ff ectively leading to a number of degeneracies such as the 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 / degeneracies inherent to closure quantities, e.g. by recently proposed multiobjective optimization schemes (Müller et al. 2023; Mus et al. 2024, 2023), or we utilize a more constraining prior information that resolves the degeneracies. For the purpose of the latter approach, the multiresolution support proved successful and has been implemented in DoG-HiT .', '3. Results': 'The M87 GMVA + ALMA images at 86GHz by resolve , DoG-HiT , and CLEAN (Lu et al. 2023) are shown in Figure 1. resolve and DoG-HiT image fits files and results are archived in zenodo 3 .', '3.1. Estimation of ring-like features in M87': 'Azoom-in into 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) are analogous to EHT M87 observations in 2017, 2018 (Event Horizon Telescope Collaboration et al. 2019a, 2024), which is the 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 amplitude by resolve and amplitude of 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 amplitude null is shorter than M87 EHT observation in 2017, 2018 (at around 3.4 G λ ). It implies that the M87 ring at 3 mm is larger than the ring at 1 mm since the baseline location of the visibility null amplitude scales inversely with the ring diameter (Lu et al. 2023). The M87 ring diameter at 86 GHz (64 + 4 -8 µ as) is estimated in Lu et al. (2023) by geometric modeling and SMILI images. Furthermore, they found that the thick ring is preferred over a thin ring by imaging analysis and model fitting.  \nIn this manuscript, we estimate the M87 ring features, such as diameter and width, by using two independent imaging methods resolve and DoG-HiT . To validate the robustness of the ring structure, posterior distribution of the M87 ring diameter, width, and ellipticity is estimated by VIDA software (Tiede et al. 2022) from 100 resolve posterior sample images. VIDA is an image feature extraction tool treating each image as a probability distribution and comparing the image to the geometrical model image (template) by utilizing Kullback-Leiber (KL) or Bhattacharyya (Bh) divergences as the objective function. From the corresponding template, the ring features of each posterior samples image are estimated. In this analysis, a flexible ring model ( CosineRingwFloor ) is used.  \nFigure 3 shows that the M87 ring diameter by resolve is 60 . 9 ± 2 . 2 µ as and 61 . 0 µ as by DoG-HiT . The estimated ring diameter by resolve and DoG-HiT is within errors of the estimation (64 + 4 -8 µ as) in (Lu et al. 2023). The width is 16 . 0 ± 0 . 9 µ as by resolve and 20 . 6 µ as by DoG-HiT . The discrepancy of ring width from resolve and DoG-HiT results from the sparse uvcoverage beyond the visibility null amplitude. We note that the ring diameter and width show anti-correlation, which is due to the finite resolution of the telescope array. The e ff ective radius of the ring decreases with a larger ring width (see the Appendix Gof Event Horizon Telescope Collaboration et al. 2019c), which explains the dependence between the estimated diameter and width. This anti-correlation is also shown by SMILI reconstructions (see Figure S14 in Lu et al. 2023). The ellipticity of the ring, τ , is defined as τ = 1 -b / a , where a is the semi-major axis lengths and b is the semi-minor axis lengths of elliptical ring. The τ is 0 . 06 ± 0 . 04 µ as by resolve and 0 . 04 by DoG-HiT , which means there is no significant ellipticity of the M87 ring.  \n2  \nRelative RA [mas]  \nFig. 1. GMVA + ALMA M87 image reconstructions at 86 GHz. Each image presents results obtained by a di ff erent algorithm. The top posterior mean image was reconstructed using 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 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, P.A. = -63 · in the bottom left corner.  \n<!-- image -->  \nDetailed description of the VIDA.jl software can be found in Tiede et al. (2022).  \nGMVA + ALMA M87 observation (2018, April 14) was conducted a week before 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  \nJy by DIFMAP , eht-imaging , and SMILI softwares (see Table 2 in Event Horizon Telescope Collaboration et al. (2024)). We note that the compact region flux density constraints from EHT 2018 data was challenging due to the lack of short baseline coverage (Event Horizon Telescope Collaboration et al. 2024). The total flux density at 3 mm is 0 . 57 ± 0 . 03 Jy on mas scales, and the flux density of the compact region at 3 mm is 0 . 33 ± 0 . 02 Jy by  \n0  \n2  \nFig. 2. The ring of the M87 image reconstruction at 86 GHz from the Figure 1. The color shows intensity in Jy mas -2 according to the linear color bar located at the top of the figure. Each image presents results obtained by a di ff erent algorithm, whose names are indicated in the lower right corners. The top left image is the resolve posterior mean image using the Bayesian self-calibration and imaging method. The top right image is the DoG-HiT reconstruction using closures only imaging. The bottom left image represents the CLEAN image with the over-resolved 37 µ as circular beam in Lu et al. (2023). The bottom right shows the SMILI image in Lu et al. (2023). All images were processed by a Gaussian interpolation.  \n<!-- image -->  \n-  \n-  \n-  \n-  \nRelative RA [ µ as]  \nresolve . 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 predictions of the numerical models, that typically indicate inhomogeneous optical depth in the compact region.', '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 bright spine along the jet axis. This central ridge-line may be related to the triple-helix structure in  \nthe jet of M87 that has been observed at larger scales (Nikonov et al. 2023).  \nHowever, 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.  \nDoG-HiT and resolve allow for the quantification of the robustness of the central spine 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 di ff erent wavelets, ultimately expressed by the multiresolution support. This allows us to estimate the robustness of the central ridge by a jackknife test, i.e.  \nFig. 3. The posterior distribution of the M87 ring diameter, thickness, and ellipticity estimated from 100 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).  \n<!-- image -->  \nwe cut the di ff use 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 (to focus the analysis on the di ff use emission), then we fitted the closure phases and closure amplitudes with DoG-HiT only varying coe ffi cients in the multiresolution support, and calculate the updated fit statistics. Next, we mask out the di ff use, central spine from the multiresolution support, and refit the closure quantities. Finally, we compare the scoring with a central spine and without.  \nThis strategy resembles a standard strategy in VLBI, 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, allowing 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 physically reasonable, opposed to 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).  \nWe obtain χ 2 cph = 1 . 104, χ 2 cla = 1 . 494 when fitting the data with a central spine, and χ 2 cph = 1 . 061, χ 2 cla = 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.  \nAn alternative perspective on the robustness of image features is o ff ered by the built-in uncertainty quantification in resolve . In Figure 4, the transverse flux intensity profiles of the M87 jet emission are depicted. The mean and standard deviation  \nof 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. The standard deviation of the pixel fluxes at the central spine is not particularly higher than the limb-brightened feature. Therefore, the noticeable central spine of the M87 jet is not detected in resolve image. The central spine in the images obtained by CLEAN in Lu et al. (2023) may be a consequence of CLEAN artifacts resulting from 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.', '3.3. Substructure of the M87 ring': "A feature that appears consistent across all four di ff erent 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 towards the top and the bottom, at consistent positions across methods. Note that the limbbrightened structure 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.  \nFirst, 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 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 section 4 and particularly figure S8 in Lu et al. (2023).  \nMultiple 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 e ff ects. Finally, the structure may be described by a residual sidelobe structure, hence essentially a consequence of the sparsity of the uv-coverage. In what follows, we will discuss each of these concerns individually.  \nThe double structure appears for all four imaging methods. Four imaging methods that were utilized here, approach the image reconstruction from four vastly di ff erent perspectives: CLEAN recovers the structure in an inverse modelling framework essentially processing a sparsity promoting regularization approach (Lannes et al. 1997), SMILI approaches the image by a weighted sum of multiple handcrafted data and regularization terms as a RMLalgorithm (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 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). Multiscale functions in DoG-HiT and Gaussian pro-  \n5  \nFig. 4. M87 jet transverse profiles and intensity map obtained by resolve . Intensity in the map is represented by false-color according to the colorbar used in Figure 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.  \n<!-- image -->  \ncess prior in resolve are flexible and does not ask for the double structure in the ring as a prior knowledge. The fact that all four independent methods using a variety of regularization / prior information 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.  \nThe structure may be a consequence of unsolved gain residuals. In fact, it is a possible issue of the alternating self-calibration and cleaning procedure to 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 DoG-HiT reconstructions which are independent against gain corruptions (closure-only imaging), making a potential cause by the calibration of the phases less likely, supported by resolve which solves for the self-calibration with imaging simultaneously.  \nFinally, could the double structure be introduced by the sparsity of the uv-coverage? That is a possibility that never can be eradicated completely, simply since the observation fundamentally misses relevant visibilities in the gaps of the uv-coverage. Consequently, there are missing information. In fact, the Fourier domain representation of a double source is a fringe (compare Fig. 5), which is exactly the kind of artifacts not fully cleaned residuals may introduce. This issue has been identified by synthetic data tests by Lu et al. (2023). Moreover, the occurrence of the double structure may be related to the lack of long westeast baselines. Nevertheless, we argue that the recovered double structure is represented by the measured visibilities. We show in Fig. 5 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  \nfitting discussions in Lu et al. 2023). Moreover, the visibility domain representation of DoG-HiT and resolve show the ellipticity of the null visibility amplitude points in Figure 5 (middle panel). We note that the location of the null visibility points in uv-domain is elliptical, with a 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, YS) at the longest west-east baselines implies that the ellipticity of the null visibility points is originated from the data.  \nTo analyze the uv-coverage sparsity corresponding 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 Figure 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 of the uvcoverage. 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. A definitive answer however may not 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.  \nIf the emission forming the ring image is dominated by the accretion disk, elongation in the direction perpendicular to the ring may be a simple geometric e ff ect for the inclined observer. While the 230 GHz ring image is strongly asymmetric with respect to the jet axis (Event Horizon Telescope Collaboration et al. 2024), our 86 GHz reconstructions exhibit high degree of symmetry. In numerical general relativistic magnetohydrodynamic (GRMHD) simulations of accretion, consistent behavior appears for retrograde accretion (negative black hole spins). In that case 230 GHz image asymmetry is driven primarily by the spin e ff ects (Event Horizon Telescope Collaboration et al. 2019d), dominant very near the event horizon. In the 86 GHz image, formed by a more extended emission, these e ff ects are balanced by the Doppler boost enhancing brightness on the opposite side of the black hole, resulting in a ring structure elongated perpendicularly to the jet axis, with bright spots on both sides of the central brightness depression.", '4. Conclusions': "Lu et al. (2023) presents observations of the core and jet in M87 observed with the GMVA + ALMA at 86GHz. The image contains a ring-like feature that looks similar to the one reported by the EHT (Event Horizon Telescope Collaboration et al. 2019a, 2024), but with 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 a 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 the 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).  \nTo validate the results reported in Lu et al. (2023), we apply two more imaging algorithms specially designed to study the robustness of 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 ,  \nnamely the probabilistic approach and the multiscale waveletbased deconvolution algorithm, allow us to quantify the robustness of recovered M87 ring and extended jet emission. We obtained the posterior distribution of ring diameter, width, and ellipticity by analysis of the resolve posterior sample images. We confirm the M87 ring-like structure at 86 GHz with the diameter of 60 . 9 ± 2 . 2 µ as, the thickness of 16 . 0 ± 0 . 9 µ as, and the ellipticity is 0 . 06 ± 0 . 04 by resolve and the diameter of 61 . 0 µ as, the thickness of 20 . 6 µ as, and the ellipticity is 0 . 04 by DoG-HiT . The estimated ring diameter is consistent with the estimate (64 + 4 -8 µ as) in Lu et al. (2023).  \nFurthermore, image reconstructions by resolve and DoG-HiT show that the ring is embedded in a strongly edgebrightened large scale jet structure agreeing 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 by resolve and DoG-HiT reconstructions is less prominant compared to 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 image domain. The validation by two independent imaging methods implies that the central spine is fainter than previous report.  \nFinally, 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 of 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 artifact originated 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 also might 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 e ff ects, as seen in numerical simulations of retrograde accretion.  \nAcknowledgements. 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). 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 data were correlated at the correlator of the MPIfR in Bonn, Germany.", 'References': 'Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017, ApJ, 838, 1 Arras, P., Bester, H. L., Perley, R. A., et al. 2021, A&A, 646, A84 Arras, P., Frank, P., Haim, P., et al. 2022, Nature Astronomy, 6, 259 Asada, K., Nakamura, M., & Pu, H.-Y. 2016, ApJ, 833, 56 Blackburn, L., Pesce, D. W., Johnson, M. D., et al. 2020, ApJ, 894, 31  \nKim et. al.: Imaging the ring and extended jet of M87  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nJy  \n/  \nmas  \n2  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nJy  \n/  \nmas  \n<!-- image -->  \nFig. 5. Reconstructions of the central feature by DoG-HiT (left panels), and resolve (right panels). Top panels: reconstructions of the compact emission. Top-medium panels: The central ring feature. Medium panels: The recovered amplitudes of the central feature, with the uv-coverage overplotted (red points). Medium-bottom panels: The central ring feature, cutted at 60% of the respective peak brightness, showcasing the double blob pattern on top of the ring. Bottom panels: The amplitudes and the uv-coverage of the double pattern alone.  \n<!-- image -->  \n0.8  \nFig. A.1. resolve image pixel-wise relative uncertainty, which is the sky brightness posterior standard deviation normalized by the posterior mean by the resolve reconstruction from the top panel of Figure 1.  \n<!-- image -->', 'Appendix A: Uncertainty estimation by resolve': 'Figure A.1 shows the pixel-wise relative uncertainty of the resolve image in Figure 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 limb-brightened are estimated. However, the relative uncertainty values in the central spine are higher compared to the limbbrightened jet emission. Figure A.2 depicts the posterior amplitude gains by resolve . 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 amplitude for YS antenna has uniformly distributed high uncertainties since only single polarization mode (LCP) was observed for YS antenna.  \nFig. 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.  \n<!-- image -->'}