© The Authors 2025

1. Introduction

For a long time following its initial discovery, the giant elliptical galaxy M87 remained merely an entry in astronomical catalogues (Messier 1781). More than a century later, observations at the Lick Observatory led to the discovery of a ‘curious straight ray’ superimposed on the diffuse emission of the galaxy (Curtis 1918), which decades later was identified as the relativistic jet emanating from the region close to the central supermassive black hole (SMBH). With the advent of radio astronomy and the growing scientific interest in active galactic nuclei (AGNs), M87 became a prime target for observations across the electromagnetic spectrum during the 20th century (see e.g. EHT MWL Science Working Group 2021 [hereafter M87 MWL2017], EHT MWL Science Working Group 2024 [hereafter M87 MWL2018], and Hada et al. 2024 for a review).

At the core of M87 lies the SMBH M87*, and its radio properties have been studied for decades across various frequencies (e.g. Reid et al. 1989; Junor et al. 1999; Doeleman et al. 2012; Hada et al. 2016; Mertens et al. 2016; Walker et al. 2018; Kim et al. 2018; Lu et al. 2023). In 2019, the Event Horizon Telescope Collaboration (EHTC) produced the first total intensity image of M87*’s shadow, using data collected during its initial observing campaign in 2017 (EHTC 2019a,b,c,d,e,f). This was followed by imaging of the linearly polarized emission (EHTC 2021a,b) and an analysis of the circular polarization near the event horizon (EHTC 2023).

The total intensity image of M87* revealed a ring of (42 ± 3) μas diameter that is brighter in the south (EHTC 2019d,f). Using results from theoretical simulations of M87*’s accretion (EHTC 2019e), it was determined that the ring size of M87* corresponds to a central black hole with a mass of (6.5 ± 0.7)×10⁹ M_⊙. These horizon-scale mass measurements are consistent with the mass inferred on much larger scales from stellar velocity dispersion measurements (Gebhardt & Thomas 2009; Gebhardt et al. 2011; EHTC 2019f; Liepold et al. 2023; Simon et al. 2024).

Follow-up observations by the Event Horizon Telescope (EHT) in 2018 (EHTC 2024b) verified the ring size, measuring a diameter of ${42}_{-3}^{+2}\mathrm{\mu }\mathrm{as}$, and confirmed the original interpretation of the ring being due to lensed emission around a SMBH. However, while the ring diameter was stable, the azimuthal structure of the ring evolved significantly. Namely, the angle of the peak brightness shifted by 30°  anti-clockwise in 2018. This rotation is consistent with expectations from numerical simulations of M87* (EHTC 2019e), which show temporal variation in the angle of peak brightness because of intrinsic variability in the accretion flow (EHTC 2025). Analysis of observations from 2009 to 2013 with prototype EHT arrays also indicated that the structure of M87* was consistent with a stable ring with the peak brightness position angle varying from year to year (Wielgus et al. 2020).

Further evidence of the emission seen from M87* being due to a hot magnetized accretion flow was provided by the linear polarization maps produced by the EHT in 2021 (EHTC 2021a,b). The inner ring was found to be linearly polarized. Most of the linear polarization was concentrated in the south-western portion of the ring in 2017, with a polarization fraction reaching ∼15% (EHTC 2021a). The observed polarization fraction is consistent with simulations in which the Faraday rotation internal to the emission region causes the de-polarization of synchrotron radiation (EHTC 2021b).

To probe the magnetic field structure of the ring, the EHT reconstructed its electric vector position angle (EVPA) pattern, observing a largely azimuthal but slightly twisted structure. This pattern is consistent with semi-analytical models that have a strong poloidal magnetic field component (Narayan et al. 2021; EHTC 2021b) and, ignoring Faraday effects, is predicted for magnetically dominated accretion flows (Chael et al. 2023). However, considerable uncertainty remains about internal and external Faraday rotation in M87*, which has been studied at much larger scales using Atacama Large Millimeter/submillimeter Array (ALMA) measurements (Goddi et al. 2021).

General relativistic magneto-hydrodynamic (GRMHD) simulations predict that the polarized emission should be dynamic around M87*, on timescales as short as weeks. The first hints of variable polarimetric properties were detected by EHTC (2021a). While the observed EVPA was largely stable, mild fluctuations in the fractional polarization were detected. Furthermore, the peak of the polarized emission shifted by ∼25° anti-clockwise between the images taken on April 5 and 11, 2017. Comparing these polarimetric properties averaged over the four observation days in 2017, EHTC (2021b) were able to constrain various GRMHD models. From this analysis, it was found that weakly magnetized accretion models performed worse than magnetically dominated ones. Assuming M87* is similar to the magnetically arrested disk (MAD) models used in EHTC (2021b), it is predicted that the polarization fraction should remain approximately stable for prograde MAD simulations and increase for retrograde systems. Furthermore, ignoring the influence of external Faraday effects, Chael et al. (2023) note that the distribution of specific azimuthal polarization modes from Palumbo et al. (2020) may depend on M87*’s black hole spin for the preferred models from EHTC (2021b). Therefore, studying the stability of M87*’s EVPA pattern across multiple years will provide valuable insights into the nature of its accretion flow and the central black hole’s properties.

This work presents the first 230 GHz multi-epoch study of M87*’s polarimetric variability. We analyse the polarimetric properties of M87* in 2018 and the total intensity and polarimetric properties of M87* in 2021. We compare these new results with the properties of M87* in 2017 to better understand its polarimetric variability. In Sect. 2 we provide a brief overview of the 2021 EHT observation campaign and its data properties. Section 3 provides the basics of polarimetric very long baseline interferometry (VLBI) imaging and details the different calibration and imaging pipelines used in this work. Section 4 describes the interferometric data properties in 2017, 2018, and 2021. In Sect. 5 we validate the polarimetric imaging pipelines used in this work, demonstrating their reliability for the different array configurations in 2017, 2018, and 2021. The first polarimetric images of M87* in 2018 and 2021, along with their multi-epoch properties, are presented in Sect. 6. Additionally, taking Northern Extended Millimetre Array (NOEMA) and Kitt Peak (KP) data from 2021 into account, we discuss the detection of non-trivial closure phases on scales between 200 μas and 1 mas, providing the first measurements of M87*’s extended structure at 230 GHz. Finally, the interpretation of the results is presented in Sect. 7, and our conclusions are given in Sect. 8.

2. Observations

The new M87* data described and analysed in full polarization in this work were collected as part of the April 9–19 EHT 2021 observing campaign on April 13 and 18. Furthermore, we analysed the M87* data from April 21, 2018, in full polarization. Previously, these 2018 data were analysed only in total intensity (Stokes ℐ) in EHTC (2024b). The 2018 M87* observations were carried out by ALMA, the Atacama Pathfinder Experiment (APEX), the Greenland Telescope (GLT), the IRAM 30 m Telescope at Pico Veleta (PV), the James Clerk Maxwell Telescope (JCMT), the Alfonso Serrano Large Millimeter Telescope (LMT), the Submillimeter Array (SMA), and the Submillimeter Telescope (SMT).

The 12 m Kitt Peak Telescope and NOEMA¹ joined EHT observations for the first time in 2021. The LMT did not participate in 2021 but resumed regular EHT observations in 2022. The 2021 array is displayed in Fig. 1. The April 18 observations are the focus of this analysis, as they include NOEMA, which did not participate in the April 13 observations due to bad weather. The April 13 data are used for consistency checks, and results are presented in Appendix A.4. 64 VLBI scans were recorded on M87* between 19:20 UT, April 17 and 11:25 UT, April 18 in 2021. Similarly, we have 64 M87* scans between 19:40 UT, April 12 and 11:40 UT, April 13. In each scan, we integrated for five minutes on the source.

Fig. 1. EHT in its 2021 configuration. Compared to the original 2017 array, GLT was added in 2018, and KP and NOEMA joined the EHT for the 2021 campaign (indicated in blue). Baselines from SPT and LMT are greyed out since SPT cannot observe M87* (only its calibrator 3C 279), and LMT did not observe in 2021.

In 2018, the EHT data recording rate was upgraded from 32 gigabits per second (Gbps) to 64 Gbps, except for the GLT, which used a 32 Gbps rate in its inaugural observations (Chen et al. 2023). As of 2021, the GLT also uses four Mark-6 units, and thus the new EHT data described in this work marks the first 64 Gbps observations with the GLT. Furthermore, due to photogrammetry and panel adjustments, the GLT 230 GHz aperture efficiency increased from 22% to 66% after the EHT 2018 observations (Koay et al. 2020, 2023a; Chen et al. 2023). In this work we utilized the upper sideband, and the two new lower sideband frequencies will be analysed in future work. 2021 is also the first EHT observations where the JCMT observed in dual polarization thanks to the new 2 sideband Namakanui receiver (Mizuno et al. 2020; Mizuno & Han 2021). The station sensitivities and metadata used for the flux density calibration of the 2021 data are described in Romero-Cañizales et al. (2025).

The baseband data from both receiver sidebands recorded by each telescope are correlated over four frequency bands centred on 213.1, 215.1, 227.1, and 229.1 GHz; each with a bandwidth of 1875 MHz, which we refer to as bands 1–4, respectively. For the 2017 observations, the 227.1 GHz and 229.1 GHz bands were referred to as low and high bands, respectively. To accommodate different frequency recording setups at different stations while maintaining a fixed visibility frequency grid with 32 sub-bands, each with 116 channels per band, the DiFX (Deller et al. 2007, 2011) software was used for correlation with the outputbands mode². This work analyses the two upper sideband frequency bands: the low band (band 3) around 227.1 GHz and the high band (band 4) around 229.1 GHz. Following a linear to circular PolConvert (Martí-Vidal et al. 2016) process, we formed visibilities in a full-polarization circular feed basis: RR^*, RL^*, LR^*, and LL^*. Combinations of these four correlation products can be used to form the four Stokes parameters, as explained in the next section. The raw post-correlation visibilities after PolConvert are made publicly available as L1 releases on the EHT website³.

3. Methods

This study made use of two calibration pipelines for cross verification: rPICARD (Janssen et al. 2018, 2019) and EHT-HOPS (Blackburn et al. 2019). Additionally, we employed seven imaging algorithms to ensure the robustness of the results: two CLEAN-based imaging algorithms, GPCAL (Park et al. 2021) and LPCAL (Shepherd 1997; EHTC 2019d); three regularized maximum likelihood (RML) algorithms, DoG-HIT (Müller & Lobanov 2022), ehtim (eht-imaging; Chael et al. 2016, 2018; EHTC 2019d), and MOEA/D (Müller et al. 2023; Mus et al. 2024a); and two Bayesian methods, Comrade (Tiede 2022) and THEMIS (Broderick et al. 2020a).

3.1. Overview

The pairwise correlations between the electric field measurements from two idealized dual-polarized circular feeds (R_jR_k^*, R_jL_k^*, L_jR_k^*, and L_jL_k^*) are related to the Stokes visibilities ($\stackrel{\sim }{\mathcal{I}}$, $\stackrel{\sim }{\mathcal{Q}}$, $\stackrel{\sim }{\mathcal{U}}$, and $\stackrel{\sim }{\mathcal{V}}$) through linear algebraic transformations. For a baseline between two stations j and k, the true-source coherency matrix can be expressed as

$$\begin{array}{c}\hfill {\rho }_{\mathit{jk}}=\left(\begin{array}{cc}{\stackrel{\sim }{\mathcal{I}}}_{\mathit{jk}}+{\stackrel{\sim }{\mathcal{V}}}_{\mathit{jk}}& {\stackrel{\sim }{\mathcal{Q}}}_{\mathit{jk}}+i{\stackrel{\sim }{\mathcal{U}}}_{\mathit{jk}}\\ {\stackrel{\sim }{\mathcal{Q}}}_{\mathit{jk}}-i{\stackrel{\sim }{\mathcal{U}}}_{\mathit{jk}}& {\stackrel{\sim }{\mathcal{I}}}_{\mathit{jk}}-{\stackrel{\sim }{\mathcal{V}}}_{\mathit{jk}}\end{array}\right)=\left(\begin{array}{cc}{R}_j{R}_k^{\ast }& R_j{L}_k^{\ast }\\ L_j{R}_k^{\ast }& L_j{L}_k^{\ast }\end{array}\right).\end{array}$$(1)

Inverting this equation yields the Stokes visibilities:

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\stackrel{\sim }{\mathcal{I}}}_{\mathit{jk}}\\ {\stackrel{\sim }{\mathcal{Q}}}_{\mathit{jk}}\\ {\stackrel{\sim }{\mathcal{U}}}_{\mathit{jk}}\\ {\stackrel{\sim }{\mathcal{V}}}_{\mathit{jk}}\end{array}\right)=\frac{1}{2}\left(\begin{array}{c}{R}_j{R}_k^{\ast }+L_j{L}_k^{\ast }\\ R_j{L}_k^{\ast }+L_j{R}_k^{\ast }\\ -i(R_j{L}_k^{\ast }-L_j{R}_k^{\ast })\\ R_j{R}_k^{\ast }-L_j{L}_k^{\ast }\end{array}\right).\end{array}$$(2)

Unfortunately, this simple relation no longer holds exactly for realistic interferometers with atmospheric and instrumental effects. We used the radio interferometric measurement equation (RIME) formalism to incorporate instrumental and atmospheric effects. The RIME formalism provides a mathematical framework that relates observed visibilities to the true sky brightness distribution while accounting for instrumental and propagation effects (Hamaker et al. 1996; Smirnov 2011). The basic form of RIME for a quasi-monochromatic signal is expressed as

$$\begin{array}{c}\hfill {\rho }_{\mathit{jk}}^{\prime }={\mathit{J}}_j{\rho }_{\mathit{jk}}{\mathit{J}}_k^{†},\end{array}$$(3)

where ρ′_jk represents the measured visibility or coherency matrix, which is a 2 × 2 matrix that encapsulates four correlations between the voltage signals received from two stations, j and k, using dual-polarized feeds, ρ_jk is the true-source visibility or coherency matrix that describes the inherent brightness of the source, J is the Jones matrix that characterizes the linear transformations that an incoming signal undergoes due to propagation and instrumental effects (Jones 1941), and the † symbol represents the conjugate transpose. Different instrumental and propagation effects are represented by distinct Jones matrices, which are multiplied in a sequence (called a Jones chain; Smirnov 2011) that reflects the physical order of these effects along the signal path.

For EHT polarimetric data, the Jones matrix formalism incorporates effects particularly critical for polarization calibration (Thompson et al. 2017). After instrumental calibration and post-processing described above, we parameterized our Jones matrices by

$$\begin{array}{c}\hfill {\mathit{J}}_j={({\mathbf{\Phi }}^{-1}\mathit{G}\mathit{D}\mathbf{\Phi })}_j.\end{array}$$(4)

Here, G is the time-dependent residual instrumental gains matrix, Φ_j is the instrumental feed rotation matrix, ϕ(t) is the feed-rotation angle and D is the constant instrumental polarization leakage matrix:

$$\begin{array}{cc}\hfill {\mathit{G}}_j& =g_{j,R}(t)·\mathrm{diag}[1,g_{j,L}(t)/g_{j,R}(t)],\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill {\mathbf{\Phi }}_j& =\mathrm{diag}[e^{-i{\varphi }_i(t)},e^{i{\varphi }_i(t)}],\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\mathit{D}}_j& =\left(\begin{array}{cc}1& d_{j,R}\\ d_{j,L}& 1\end{array}\right).\hfill \end{array}$$(7)

In the following, we refer to g_j, R(t) as the station-based, time-variable gains, and g_j, L(t)/g_j, R(t) as the station-based, time-variable gain ratios. The feed rotation angle at site j depends on the source elevation f_j, el, the parallactic angle f_j, par, and an offset ϕ_j, offset,

$$\begin{array}{c}\hfill {\varphi }_j(t)=f_{j,\mathrm{el}}{\theta }_{j,\mathrm{el}}(t)+f_{j,\mathrm{par}}{\psi }_{j,\mathrm{par}}(t)+{\varphi }_{j,\mathrm{off}}.\end{array}$$(8)

We analysed the conjugate closure trace products, 𝒞 (Broderick et al. 2020a) for each observation to assess the presence of polarized emission in the data. Conjugate closure trace products are independent of any time-dependent stationized instrumental effects that can be represented as a Jones matrix, including residual instrumental gain errors (Eq. (5)) and instrumental polarization leakage errors (Eq. (7)) and are defined on quadrangles of four stations, j, k, l, and m, as

$$\begin{array}{c}\hfill {\mathcal{C}}_{\mathit{jklm}}={\mathcal{T}}_{\mathit{jklm}}{\mathcal{T}}_{\mathit{jmlk}},\end{array}$$(9)

where ${\mathcal{T}}_{\mathit{jklm}}=\frac{1}{2}\mathrm{tr}\left({\rho }_{\mathit{jk}}{\rho }_{\mathit{lk}}^{-1}{\rho }_{\mathit{lm}}{\rho }_{\mathit{jm}}^{-1}\right)$ is the closure trace. If no polarized emission is present, arg(𝒞) will be zero.

Image reconstruction refers to inferring from the measured coherency matrices/visibilities, the set of Stokes parameter maps ℐ(x), 𝒬(x), 𝒰(x), 𝒱(x), which fully characterize the polarized state of electromagnetic radiation at a given spatial coordinate x = (x, y). ℐ(x) gives total intensity, 𝒬(x) measures the difference between horizontal and vertical linear polarization, 𝒰(x) quantifies the difference between light polarized at 45° and −45°, and 𝒱(x) represents the level of circular polarization. The Stokes parameters are related to the Stokes visibilities through the Fourier transform (van Cittert 1934; Zernike 1938),

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathcal{I}}}_{\mathit{jk}}={\displaystyle \int \mathcal{I}}(\mathit{x})e^{2\pi i{\mathit{u}}_{\mathit{jk}}·\mathit{x}}\mathrm{d}x\mathrm{d}y,\end{array}$$(10)

where u_jk is the projected baseline between stations j and k. In addition, we also defined the total intensity closure phases (CPh), ψ_jkl, and log-closure amplitudes (LCA), 𝒜_jklm, which are insensitive to overall gain corruptions (see e.g. Rogers et al. 1974; Blackburn et al. 2020):

$$\begin{array}{cc}\hfill {\psi }_{\mathit{jkl}}& =arg\left[{\stackrel{\sim }{\mathcal{I}}}_{\mathit{jk}}^{\prime }{\stackrel{\sim }{\mathcal{I}}}_{\mathit{kl}}^{\prime }{\stackrel{\sim }{\mathcal{I}}}_{\mathit{lj}}^{\prime }\right]\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\mathcal{A}}_{\mathit{jklm}}& =log\left[\frac{|{\stackrel{\sim }{\mathcal{I}}}_{\mathit{jk}}^{\prime }||{\stackrel{\sim }{\mathcal{I}}}_{\mathit{lm}}^{\prime }|}{|{\stackrel{\sim }{\mathcal{I}}}_{\mathit{jl}}^{\prime }||{\stackrel{\sim }{\mathcal{I}}}_{\mathit{km}}^{\prime }|}\right],\hfill \end{array}$$(12)

where $\stackrel{\sim }{\mathcal{I}}{\prime }_{\mathit{jk}}=\frac{1}{2}\mathrm{tr}(\rho {\prime }_{\mathit{jk}})$ the approximate Stokes I visibility for baseline j,  k.

For polarization, the complex linear polarization is defined as

$$\begin{array}{c}\hfill \mathcal{P}\equiv \mathcal{Q}+i\mathcal{U}=\mathcal{I}|m|e^{2i\chi },\end{array}$$(13)

where m = (𝒬 + i𝒰)/ℐ is the linear polarization fraction, and χ = 0.5arg(𝒫) is the EVPA measured east of north on the sky. The circular polarization fraction is given by v = 𝒱/ℐ. In the visibility domain, we can define similar interferometric quantities (Johnson et al. 2015),

$$\begin{array}{cc}\hfill \breve{m}& =\frac{\stackrel{\sim }{\mathcal{Q}}+i\stackrel{\sim }{\mathcal{U}}}{\stackrel{\sim }{\mathcal{I}}}=\frac{2R{L}^{\ast }}{R{R}^{\ast }+L{L}^{\ast }},\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill \breve{\chi }& =\frac{1}{2}arg\breve{m}\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill \breve{v}& =\frac{\stackrel{\sim }{\mathcal{V}}}{\stackrel{\sim }{\mathcal{I}}}=\frac{R{R}^{\ast }-L{L}^{\ast }}{R{R}^{\ast }+L{L}^{\ast }}·\hfill \end{array}$$(16)

Note that $\breve{m},\breve{v}$ are not the Fourier transforms of m,  v.

The unresolved (image-integrated) linear and circular polarization fractions, as well as their resolved (image-averaged) counterparts, with ∑_i indicating a sum over all image pixels, are given by

$$\begin{array}{cc}\hfill {|m|}_{\mathrm{net}}& =\frac{1}{{\sum }_i{\mathcal{I}}_i}{[{\left({\displaystyle \sum _i}{\mathcal{Q}}_i\right)}^2+{\left({\displaystyle \sum _i}{\mathcal{U}}_i\right)}^2]}^{1/2},\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill v_{\mathrm{net}}& =\frac{1}{{\sum }_i{\mathcal{I}}_i}{\displaystyle \sum _i}{\mathcal{V}}_i,\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill ⟨|m|⟩& =\frac{1}{{\sum }_i{\mathcal{I}}_i}{\displaystyle \sum _i}{({\mathcal{Q}}_i^2+{\mathcal{U}}_i^2)}^{1/2},\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill ⟨|v|⟩& =\frac{1}{{\sum }_i{\mathcal{I}}_i}{\displaystyle \sum _i}\left(\left|\frac{{\mathcal{V}}_i}{{\mathcal{I}}_i}\right|{\mathcal{I}}_i\right).\hfill \end{array}$$(20)

Note that ⟨|m|⟩ and ⟨|v|⟩ are sensitive to image resolution, i.e. the restoring beam size, while |m|_net and v_net remain unaffected by convolution.

Another useful parameter for quantifying polarization structures and comparing polarimetric images is the complex coefficient of the azimuthal mode decomposition of 𝒫 given by (Palumbo et al. 2020)

$$\begin{array}{c}\hfill {\beta }_m=\frac{1}{{\mathit{I}}_{\mathrm{ann}}}{\displaystyle {\int }_{{\rho }_{\min }}^{{\rho }_{\max }}{\int }_0^{2\pi }}\mathcal{P}(\rho ,\phi )e^{-im\phi }\rho d\phi d\rho ,\end{array}$$(21)

where (ρ, φ) represent the polar coordinates on the image plane, I_ann is the Stokes ℐ flux density contained within the annulus between the minimum and maximum radii ρ_{min,  max}. In this paper, ρ_min is set to zero, while ρ_max is set to 45 μas to focus on the compact core emission. Note that we defined the EVPA helicity of the ring’s polarization pattern as the sign of ∠β₂ following conventions from Palumbo et al. (2020) and Chael et al. (2023). In addition to |m|_net, |v|_net, ⟨|m|⟩, and ⟨|v|⟩, the amplitude and phase of the second coefficient, |β₂| and ∠β₂ are useful parameters to score the GRMHD accretion models against the EHT data. This parameter is useful for distinguishing between accretion states (Palumbo et al. 2020) and will be used in a companion paper focusing on the theoretical interpretation of our results (Chael et al., in prep.).

3.2. Calibration and post-processing

Two pipelines are used for the post-correlation signal stabilization and data reduction (Janssen et al. 2022): rPICARD described in Janssen et al. (2018, 2019) and EHT-HOPS described in Blackburn et al. (2019). The updated data processing steps needed to address unique data issues encountered in the 2021 observations are described in Appendix B.

3.3. Polarized imaging

In this study seven different polarized imaging algorithms from three different frameworks were used. Here, we present a concise overview of the methods (see also Table 1 for a summary) and provide a more complete description of the methods in Appendix C.

For the two CLEAN-based pipelines, the first involved automated imaging with Difmap and polarimetric calibration procedures using GPCAL (see Appendix C.1.1 for more details) and collectively will be denoted as GPCAL. The second method involved another CLEAN procedure with Difmap but using different hyperparameter settings and instrumental polarization calibration using LPCAL (see Appendix C.1.2 for more details) and collectively will be denoted as LPCAL. Both methods achieved the final total intensity images through iterative CLEAN and self-calibration. The primary difference between these methods lies in their polarized imaging approach: GPCAL employed an automated parameter survey to produce a set of images, as seen in previous EHT imaging studies of M87 (EHTC 2019d, 2024b), and selected the representative image based on closure chi-squares. Similarly, LPCAL involved imaging involving a small survey over different hyperparameters where the representative image was chosen by minimizing the closure chi-square. However, the specific hyperparameter survey and the overall imaging procedure differed from GPCAL. This allowed us to test the sensitivity of the CLEAN reconstructions to differentassumptions.

For leakage calibration, GPCAL derived initial leakage solutions using the ‘similarity approximation’, which assumes that the linear polarization structure is proportional to the total intensity structure within each sub-model. The solutions were then refined through iterative linear polarization imaging, leakage solution estimation, and correction. LPCAL used the standard CALIB and LPCAL AIPS tasks to estimate the leakages and apply the corrections to the data, and then used Difmap to make the final ℐ, 𝒬, and 𝒰 maps. While both GPCAL and LPCAL used CLEAN to obtain final polarized images, their differing assumptions in deriving polarimetric leakages introduce a measure of uncertainty into the CLEAN polarimetric reconstructions. Finally, for both the GPCAL and LPCAL pipelines, the fiducial images were created by taking the final set of cleancomponents and blurring them with a 20 μas Gaussian beam similar to the nominal resolution of the EHT array (see Fig. 2).

Fig. 2. (u, v) coverage of M87* during the 2017 (left), 2018 (middle), and 2021 (right) campaigns for the band 3 (227.1 GHz) observations. The ticks show each year’s interferometric EVPA and the colour of the observed interferometric fractional linear polarization, after de-biasing for thermal noise and applying leakage corrections. The leakage terms used are the fiducial values from EHTC (2021a) in 2017 and the THEMIS leakage solutions in 2018 and 2021. The clusters of high polarization fractions in 2018 come from the GLT, which was underperforming in 2018 due to an incomplete commissioning at that time, as described in Koay et al. (2023a).

The RML framework minimizes a weighted sum of multiple objectives, data fidelity functionals (loss functionals, χ²), and regularization terms, R:

$$\begin{array}{c}\hfill \widehat{I}\in {\mathrm{argmin}}_I{\displaystyle \sum _d}{\alpha }_d{\chi }_d^2(I)+{\displaystyle \sum _x}{\beta }_x{R}_x(I).\end{array}$$(22)

The common minimization of data terms and regularization terms ensures a solution that matches observed data and is favourable with respect to the hand-crafted regularization terms. This framework has been realized in two different methods: ehtim (Chael et al. 2016, 2018) and DoG-HIT (Müller & Lobanov 2022, 2023a). ehtim has been used in all previous EHT studies on M87* (EHTC 2019d, 2021a, 2023, 2024b) and Sgr A* (EHTC 2022, 2024c), as well as for imaging various AGN sources (Kim et al. 2020; Janssen et al. 2021; Issaoun et al. 2022; Jorstad et al. 2023). DoG-HIT has also been applied to Sgr A* (EHTC 2024c) and lower frequency M87* data (Kim et al. 2025). This paper is the first time MOEA/D was applied to one of the primary EHT targets. The three RML approaches differ in the regularizer terms used, the optimality concept applied, and the minimization procedure, but use similar calibration procedures. For each RML pipeline, the representative image is the optimal image according to its loss function. For more information, see Appendix C.2.

To select the regularizer hyperparameters for ehtim, we used a combination that was consistent with the top set results from EHTC (2019d, 2024b) for total intensity and EHTC (2021a) for linear polarization. In Appendix C.2 we verify that this combination performs well on synthetic data. However, unlike previous EHT publications (EHTC 2019d, 2024b), no attempt was made to create the ‘top set’ of hyperparameters within the ehtim pipeline. A reduced parameter search has been performed as spot check instead, as described in Appendix C.2. Given that there are seven different imaging algorithms, including two Bayesian frameworks, there was sufficient diversity in our imaging algorithms to explore model uncertainty. Furthermore, we found that the uncertainty we reported from the 2017 image reconstructions (see Sect. 6) was consistent with those previously inferred using the ‘top set’ approach in EHTC (2019d).

Finally, two Bayesian polarized imaging pipelines are used: THEMIS and Comrade. Both methods jointly solve for all four Stokes parameters and instrumental terms, such as gains and leakage corrections. Both methods used a rasterized image in all four Stokes parameters, but their priors differed to test the robustness of M87*’s image. For both methods, every station in the array assumed gains that vary independently for each scan and frequency band. Leakages were assumed to be constant for each observing day, but could differ across frequency bands. Finally, the gain ratios were handled differently in Comrade and THEMIS. For each site in the array, Comrade fit an independent gain ratio for each scan and frequency band, while THEMIS enforced that the gain ratios were unity. For both Bayesian methods, the representative image was given by the mean image computed from their respective posterior samples.

3.4. Feature extraction

To evaluate the images, we focused on each polarized imaging algorithm’s ability to accurately measure the key image structural parameters and integrated polarization quantities. Quantifying the reliability of these quantities was critical since they are used to score the theoretical simulations in Chael et al. (in prep.) as well as the overall image quality across different sets of synthetic data.

For the Stokes ℐ images, we extracted the following parameters: the zero-width or δ ring equivalent diameter, $\widehat{d}$, (defined below), brightness asymmetry, A, brightness position angle, η, and compact ring flux density, F_com. We used the template matching algorithm from Tiede et al. (2022) to extract these parameters. The radial profile of the template was a Gaussian distribution whose diameter, d₀, is defined as the diameter of the peak brightness. To harmonize the measured ring size across different pipelines, we converted d₀ to the δ ring equivalentdiameter:

$$\begin{array}{c}\hfill \widehat{d}=d_0+\frac{1}{4ln(2)}\frac{w^2}{d_0},\end{array}$$(23)

which is derived in EHTC (2019d, Appendix G). This equation removes the approximate effect of blurring a δ ring with a Gaussian beam of size w, which we defined as the full width at half maximum (FWHM) of the Gaussian ring.

The template azimuthal brightness profile was assumed to be (Tiede et al. 2022),

$$\begin{array}{c}\hfill S(\varphi )=1-2{\displaystyle \sum _{k=1}^4}{A}_{k}cos[k(\varphi -{\eta }_k)].\end{array}$$(24)

We restricted A_k ∈ [0.0, 0.5] to limit the amount of negative brightness in the image. Following EHTC (2019d), the ring brightness asymmetry, A, is defined as A₁, and the ring brightness position angle, η, is defined as η₁.

To determine the optimal template, we optimized the cross-correlation coefficient

$$\begin{array}{c}\hfill {\rho }_{\mathcal{I}}(\mathit{I},\mathit{T})=\frac{⟨\mathit{I}-⟨I⟩,\mathit{T}-⟨T⟩⟩}{||\mathit{I}-⟨I⟩||||\mathit{T}-⟨T⟩)||},\end{array}$$(25)

where T is the template image and I is the image reconstruction. Given the optimal template, the compact ring flux density, F_com, is defined as the total flux density within a 90 μas diameter disk, whose centre matches the fitted ring centre.

Finally, since each polarized imaging algorithm has a different image centre, field of view, and intrinsic resolution, we first normalized the image reconstructions across each method based on their total intensity images. To create a uniform image resolution, we selected a reference image. For the synthetic data tests described below, we used the ground-truth image blurred with a 20 μas FWHM Gaussian beam as the reference. This choice matched the conventions in EHTC (2019d, 2021a). For the M87* reconstructions, we used LPCAL as our reference image. Note that the two CLEAN pipelines gave consistent results. However, we chose LPCAL because it is an established pipeline in VLBI. Given the reference image, the reconstructions for each method are blurred to maximize their total intensity cross-correlation, Eq. (25). Note that this blurring may differ from the intrinsic resolution of the linear polarization maps. Given these harmonized images, we used the template matching procedure described above to estimate the total intensity parameters, $\{\widehat{d},w,A,\eta \}$, and compute the centre of the ring. We then calculated |m_net|, ⟨|m|⟩, β_1, 2 and F_com by integrating radially about the fitted ring centre to a radius of 45 μas and azimuthally over 2π radians. Finally, to estimate the global polarization fidelity of the polarization reconstructions, we also computed the linear polarization cross-correlation between a linear polarization map 𝒫 and a reference map 𝒫₀ following EHTC (2021a),

$$\begin{array}{c}\hfill {\rho }_{\mathcal{P}}(\mathcal{P},{\mathcal{P}}_0)=\frac{\mathrm{Re}[⟨\mathcal{P},{\mathcal{P}}_0⟩]}{||\mathcal{P}||||{\mathcal{P}}_0||},\end{array}$$(26)

where ⟨ ⋅ , ⋅ ⟩ is the complex inner product, and ∥ ⋅ ∥ is the complex norm. If 𝒫 and 𝒫₀ are co-linear, then ρ_𝒫 is unity, and if they are not, then ρ_𝒫 < 1 by the Cauchy-Schwarz inequality.

4. Data properties

Figure 2 shows the (u, v) coverage for EHT observations in 2017, 2018, and 2021. The ticks denote the interferometric EVPA $\breve{\chi }$ and the linear polarization fraction $|\breve{m}|$ as defined in Eq. (13) after correcting for polarization leakage and amplitude noise bias. The 2018 data show some points with high $|\breve{m}|$, which are all related to the low signal-to-noise ratio data at baselines to the GLT. These points should be interpreted with care, as the GLT had a lower aperture efficiency in 2018 than in 2021. In particular, when considering only common regions in the (u, v) space across the years, 2017 typically displays a higher $|\breve{m}|$.

The M87* total intensity visibility amplitudes and phases as a function of (u, v) distance are shown in Fig. 3. The April 18, 2021, data have the best (u, v) coverage of any EHT observation to date. Additionally, we achieved a high calibration quality with very coherent phases and few outliers in amplitudes that resulted from residual gain errors. The data show the characteristic secondary peak beyond a deep amplitude minimum of a ring-like structure at ∼3.4 Gλ. Finally, the addition of the NOEMA station implies that the secondary visibility null located at ≈8.3 Gλ is probed by two baselines, PV-Hawaii and NOEMA-Hawaii, forming the first non-trivial closure triangle at such resolution scales.

Fig. 3. Band 3 (227.1 GHz) M87* total intensity amplitude and phase data measured in 2017, 2018, and 2021. The calibration applied is the pipeline-based signal stabilization and a priori flux density calibration without image-based self-calibration. The 2017 and 2018 data were produced by EHT-HOPS, and the 2021 data by rPICARD. The grey bands around 3.4 Gλ indicate the approximate location of the first visibility minima.

Compared to 2017 and 2018, the absence of LMT in the 2021 array implies that information from intermediate spatial scales from the ∼1.5 Gλ to ∼2 Gλ (∼100 μas) LMT-SMT baseline are poorly constrained. At the same time, the addition of NOEMA provides a new ∼700 Mλ (∼350 μas) baseline to PV, and KP-SMT provides a new short ∼100 Mλ (∼2 mas) baseline. As in previous years, the ∼3.5 Gλ GLT-PV baseline probes the source asymmetry; the addition of NOEMA allows us to constrain the asymmetry with a new high-fidelity triangle (GLT-PV-NOEMA).

Previous EHT observations of M87* have shown significant polarization structures. Since polarization signals are particularly prone to calibration errors, we inspected the conjugate closure trace phases arg(𝒞) defined in Eq. (9). The conjugate closure trace quadrangle ALMA-APEX-SMA-SMT is present in all three years, which we show in the top panel of Fig. 4. The bottom left and middle panels show the conjugate closure trace for the quadrangle ALMA-APEX-LMT-SMT for 2017 and 2018, while the bottom right panel shows the quadrangle ALMA-PV-NOEMA-GLT for 2021. Significant deviations from zero are visible in 2017 and 2021, but not in 2018. This result indicated the presence of significant and changing polarization structures throughout the years and provides calibration-independent evidence that M87* was de-polarized in 2018.

Fig. 4. Comparisons of the conjugate closure trace phases across the three observing campaigns for two quadrangles. Deviations from arg(𝒞) = 0 indicate the presence of significant polarization structure, independent of leakage and gain calibration. Closed and open points show low- and high-band data, respectively, offset in time for clarity. The model predictions from the low-band reconstructions described in Sect. 3 are shown for reference (for the Bayesian methods, THEMIS and Comrade, the 95-percentile range is shown along with five draws from the posterior). Note the variation in the vertical axis ranges between panels. Because LMT did not observe during the 2021 campaign, we show the ALMA-PV-NOEMA-GLT quadrangle in 2021, for which clear deviations from zero are apparent.

5. Image validation

We conducted two separate synthetic polarized imaging and calibration tests to validate our imaging methods. These tests ensured that the features we extracted from M87* were not significantly biased. Furthermore, the two tests were designed to investigate different potential systematics in the data, source model, and imaging methodologies.

The first set was developed to test the robustness of the different polarized imaging methods to accretion turbulence and small-scale structure below the intrinsic EHT resolution. A significant component was ensuring that the image reconstructions were of high quality and that the parameter estimates for the total intensity and polarized quantities were accurately determined. The ground truth images and the data generation process were blinded for all imaging teams to prevent human biases from affecting the image reconstructions. However, an issue was identified during evaluation of the MOEA/D results after the data had been un-blinded. This issue significantly impacted the MOEA/D results and was not due to any assumptions in the imaging pipeline but to a bug in one of its polarized regularizers. After fixing this bug, the MOEA/D results were rerun, which is what is shown below. As a result, the MOEA/D results are notconsidered blinded. Note that this problem did not impact the six other imaging methods.

The second set of tests focused on checking whether the polarized image reconstructions could be biased by over-resolved polarized emission from M87*’s extended jet. For this test, we created two versions of synthetic data that shared the same core component but differed in their polarization properties within the extended jet component. Note that these tests were not blinded during the data generation process. In the main text of this paper, we present results from the blinded synthetic data tests. Results from the extended emission test are presented in Appendix D.3.

5.1. Data generation

We generated synthetic data for each source model with the software ehtim using the (u, v) coverage from the April 11, 2017, April 21, 2018, and April 18, 2021, observations. These synthetic data were generated for band 3 for all three years. The complex visibilities were generated by sampling the Fourier transform of each ground truth image with the observation’s (u, v) coverage. The sampled complex visibilities were then corrupted with (i) time-stable polarimetric leakage terms, (ii) station-based, time-variable gains in the visibility phase and amplitude, (iii) station-based, time-variable complex gain ratios of the right-left circular polarization gains, and (iv) baseline-based thermal noise estimated from the corresponding EHT observations.

This generation process is identical to that of EHTC (2023) for the 2017 synthetic data. For 2018 and 2021, we changed the generation process to accommodate the changes in data properties throughout the three years of observations. These changes included adding or removing stations present during each year’s observations, modifying station mount types, and replicating the significant JCMT polarimetric leakage observed during 2018 (Appendix B).

For the blinded test source models, three snapshots from a polarized MAD GRMHD simulation with i = 17° , a = −0.5, R_low = 10, and R_high = 40 were used as the ground-truth images. This GRMHD simulation was chosen because it passed the EHT 2017 and 2018 multi-year theoretical constraints from EHTC (2025). For each year of the synthetic data, we used a different random snapshot. For more details, see Appendix D.

5.2. Synthetic data results

The blinded parameter estimation results are shown in Fig. 5. We found that for every year, the total intensity cross-correlation with the truth is > 0.975 for all methods, demonstrating their high quality. Each method recovered the true δ ring diameter to within 2 μas, the brightness asymmetry to within ∼0.05, and the ring position angle estimate, η, to within 5 − 10° . For the compact flux density, we found a larger spread across methods. This is not unexpected due to the lack of intermediate baselines in the EHT array, making compact flux constraints sensitive to the image priors and gain solutions as well as priors, as discussed in more detail below.

Fig. 5. Comparison of extracted parameters for the blinded synthetic data test. All methods are blurred to match the ground truth GRMHD, blurred to 20 μas before the parameters are estimated. The markers show the results for each method. For the two Bayesian methods, Comrade and THEMIS, the markers denote the median, and the error bars the 95% credible interval. Since the synthetic data only consider a single frequency, the non-Bayesian methods only produce a single image. The dashed black line shows the ground truth, estimated from the true on-sky image blurred with a 20 μas Gaussian beam.

For the linear polarization reconstructions, the polarized cross-correlation improved dramatically from 2017 to 2021. In 2017, 0.7 < ρ_𝒫 < 0.9, while in 2018 and 2021, ρ_𝒫 > 0.9 for all methods. This demonstrated the improved polarized imaging capabilities of the 2018 and 2021 EHT arrays. The image reconstructions and more details can be found in Appendix D.2 and Fig. D.1.

Analysing the derived image-averaged polarimetric quantities in Fig. 5, we found that the true |m_net| is contained within the spread across methods each year. Furthermore, Comrade’s posterior contained the truth within its 95% contours every year, while the estimates of LPCAL and ehtim are within 0.25% of the truth each year. Similarly to |m_net|, the image-integrated EVPA χ is also recovered each year, and all methods recovered the truth to within 10°.

For ⟨|m|⟩, we found a slightly more complicated result. Unlike |m_net| where the reconstructions tended to be distributed around the truth, we found that ⟨|m|⟩ tends to be biased low for all methods except THEMIS. This bias reflects the difficulty of measuring ⟨|m|⟩ due to its sensitivity to the linearly polarized resolution, and field of view of the reconstructions (see Appendix G of EHTC 2023, for a related discussion). Although the images were blurred to match their resolution, this was based on total intensity, which may differ from the resolution of the polarized maps. The magnitude of this bias tended to decrease from 2017 to 2021, and in 2021, the results of four of the seven methods were close (< 0.5%) to the true value. Therefore, after combining the results across all methods, our estimates of ⟨|m|⟩ recovered the true value.

Analysing the azimuthal structure of the ring reconstruction, the phases of the first two β modes are recovered each year. Furthermore, the improved coverage in 2018 and 2021, compared to 2017, increased the precision and accuracy of the measurements of ∠β₂. That is, we found a 10°  dispersion around the truth for all methods. The amplitudes of the first two β modes are also recovered, although the spread across methods is more pronounced. Specifically, for |β₁|, we found that the estimates from DoG-HIT and MOEA/D exhibited significant deviations from the true values. Furthermore, we found that, similar to ⟨|m|⟩, some methods tended to be biased towards lower values, specifically ehtim and MOEA/D for |β_1, 2| and DoG-HIT for |β₂| across all years. Similarly, Comrade was biased low for |β₁| in 2018 and |β₂| in 2017, although the truth is contained within the 99% credible interval in both cases. Since ⟨|m|⟩, and |β_i| are roughly proportional, the bias in the β mode amplitudes is likely of a similar origin.

In summary, even if individual methods were sometimes biased, the combined estimates across all methods contained the truth. Therefore, to estimate the parameters of M87*, we combined the estimates from all methods. The estimates are combined by concatenating the results across all methods, weighted by the inverse number of samples produced. Specifically, for the non-Bayesian methods, the inverse weights are equal to two since each method produced an estimate for the band 3 and band 4 data. For the Bayesian methods, the inverse weight is given by the number of posterior samples from each method and frequency band. Given this set of samples, we then computed the 95% percentile range to estimate the combined uncertainty.

In Appendix D and Fig. D.2, we also compare the leakage recovered from each method to the true value. In general, we found that the Bayesian methods recovered the leakage for every station to within 1% for the blinded synthetic data. For the non-Bayesian methods, we observed more scatter, especially for GLT in 2018. This can be understood in light of the relatively small parallactic angle coverage of GLT for M87* (∼15° ) compared to other EHT sites (50–100°) over a whole observation track. As a result, the Bayesian methods reported a relatively large leakage uncertainty for GLT compared to the other sites, but the true value was still recovered. However, the non-Bayesian method leakage estimates are more prone to local minima since they report a single value rather than characterizing the parameter space, making them more susceptible to biases. Regardless of these discrepancies, the different leakage estimates for GLT did not impact our results.

Finally, for the non-blinded extended emission tests, we found similar results. Namely, each method’s image reconstructions were not significantly impacted by the presence of extended emission. For more detailed information, see Appendix D.3.

6. Results

Figure 6 presents the total intensity and linear polarization maps of M87* in 2017, 2018, and 2021. Unless stated otherwise, the 2017 results were obtained from the EHT-HOPS band 3, April 11 data from EHTC (2021a). Band 3 was chosen for consistency with EHTC (2021a). The 2018 results were obtained from the HOPS reduction of bands 3 and 4 on April 21, 2018, in EHTC (2024b). The 2021 results from bands 3 and 4 on April 18 were produced using the rPICARD reduction. rPICARD was chosen in 2021 due to its ability to handle NOEMA’s phase jumps (see Appendix B and von Fellenberg et al. 2025).

Fig. 6. Fiducial images for 2017 (band 3), 2018 (bands 3 and 4), and 2021 (bands 3 and 4). The images were produced by averaging the reconstructions over the methods described in Sect. 3. Each method’s image reconstruction has been blurred to match the resolution of the representative LPCAL image. Top row: Polarization ‘field lines’ overlaid on the total intensity image. Second row: Total intensity image in grey scale with the contours showing the 22.5%, 45%, 67.5%, and 90% peak brightness levels, overlaid with polarization ticks. The polarization ticks indicate the EVPA, the tick length is proportional to the linear polarization intensity, and their colour indicates the linear polarization fraction. Polarization ticks are only shown in regions where the total intensity is > 10% of the maximum brightness and the linear polarization brightness is > 10% of the peak linear polarized brightness. Bottom row: Total linear polarized brightness, |𝒫|.

For each year, the images shown in Fig. 6 are produced by averaging all the imaging methods, blurred to a common resolution. For individual results, we refer the reader to Appendix F. Generally, the individual methods blurred to a common resolution are consistent with the averaged results, especially for the 2021 data, where the improved array strongly constrains the ring emission.

The parameter and feature extraction results in Table 2 show the EHT constraints averaged over all methods for each year. As mentioned above, the ranges are computed by averaging the results from each method, weighted equally, and calculating the 95% percentile range.

For 2017, the values in this work are consistent with the original values reported in EHTC (2021a). However, the uncertainty of our estimates is smaller than the original values EHTC (2021a). This reduction is expected since EHTC (2021a) analysed all four days of observations, while this paper only considered April 11, 2017. Therefore, we used the reported polarization values from EHTC (2021a) and ring orientation from EHTC (2019d,f) in Table 2 for a conservative estimate of M87*’s appearance in 2017. Note that the quantitative measurements of M87*’s polarized emission using the methods in this paper are described below in Sect. 6.2.

Finally, note that in addition to estimating the linear polarization structure, every method also estimated the leakage for each station. Unlike EHTC (2021a), we did not attempt multisource fitting to constrain the leakages for the EHT array. Instead, each imaging method estimated the polarization leakage for each site as a byproduct of polarimetric imaging. In general, we found consistent D-terms between imaging methods in 2018 and 2021, except for GLT in 2018. The measured GLT leakage displayed considerable scatter among the non-Bayesian methods. However, this scatter was also observed in the synthetic data tests and is explained by GLT’s poor parallactic angle coverage and their leakage reconstruction approach as we describe inSect. 5.

6.1. Stokes I results

Figure 6 demonstrates that M87* consistently shows a central depression – the black hole shadow – across all observations. The ring around this central depression is consistently asymmetric in brightness, always appearing brighter in the south. This aligns with expectations from simulations of accretion physics, given the location of the large-scale jet observed at lower frequencies. The size of the ring remains stable over the years and is consistent with previous results (EHTC 2019d, 2024b).

Specifically, Fig. 7 shows that the diameter of the ring is consistent between methods, although 2017 has the most scatter due to its sparser coverage. After averaging over methods, we found a diameter estimate of [42.0 μas,  46.4 μas] for 2017, [40.7 μas,  44.4 μas] for 2018, and [43.1 μas,  44.4 μas] for 2021. Therefore, as predicted by our theoretical understanding of accretion physics, the diameter of the ring is stable across all three observing epochs, confirming the findings in EHTC (2024b). Note that in EHTC (2024b) two estimates of M87*’s diameter were provided. Our imaging results are consistent with their method averaged estimate of $\sim 43.3_{-3.1}^{+1.5}\mathrm{\mu }\mathrm{as}$. This similarity demonstrates the robustness of the EHT ring diameter estimates, given that our results utilized different imaging and analysis pipelines.

Fig. 7. Extracted parameters of M87* across the three observations and averaged over band 3 and band 4. Each panel shows the results from April 11, 2017, April 21, 2018, and April 18, 2021. For the non-Bayesian methods, the error bars show the spread between the high- and low-band estimates and are not a measure of the statistical uncertainty of the image reconstruction. The error bars for the Bayesian methods show the 95% credible interval around the median and measure the statistical uncertainty of the reconstructions due to thermal noise, instrumental effects, incomplete coverage, and frequency dependence. The solid black line is the median, and the grey band is the 95-percentile range from all image reconstructions, with each method weighted inversely to the number of images they produced. Note that each image reconstruction has been blurred to match the resolution of the LPCAL reconstruction before the parameters were estimated.

Unlike the size of the ring, we found that its azimuthal structure changed from 2017 to 2018/2021. As first noted in EHTC (2024b) in 2018, the brightness peak changed from 2017 to 2018, moving anti-clockwise by ∼45° . From Fig. 6 this shift appears to be due to the bright region in the eastern part of the ring disappearing. In 2021, we found that M87*’s brightness profile is remarkably similar to the 2018 image. Quantitatively, calculating the cross-correlation between the fiducial 2017 and 2018 images gives 0.948, while the cross-correlation between 2018 and 2021 gives 0.992. This value is similar to the average cross-correlation between all pairs of imaging methods in 2021.

Analysing the total intensity ring properties in Fig. 7 and Table 2, we found that the ring brightness asymmetry, A, is consistent in 2018 and 2021 and marginally discrepant in 2017. The difference in azimuthal brightness distribution is seen most clearly in the evolution of the position angle η from 2017 to 2018/2021. The measured position angle was [166° , 175° ] in 2017, while in 2018 and 2021, η was between 200°  and 218°  over the two years. The stability of η from 2018 to 2021 is notable, as it aligns with the expected brightness maximum predicted by theoretical simulations, given the location of the low-frequency jet (EHTC 2019e, 2024b, 2025).

Finally, we measured significant disagreements in the compact flux density in the 2021 images. ehtim measured a compact flux density of around 0.5 Jy, similar to the range found in 2017 and 2018, while the other methods found a compact flux between 0.7 Jy and 0.95 Jy. A similar disagreement in compact flux was reported in EHTC (2024b), where a detailed analysis revealed that the compact flux measured by EHT was very prior-dependent, producing values ranging from 0.3 to 1.1 Jy. A significant factor in this uncertainty, beyond residual calibration issues, is the lack of intermediate baselines in the 2021 EHT array. For the 2017 and 2018 EHT array, the LMT-SMT baselines directly probed the emission on ∼100 μas scales. In 2021, the EHT gained a new intermediate baseline, NOEMA-PV (sensitive to ∼300 μas scales) and a short baseline KP-SMT (sensitive to ∼2-3 mas scales). Unfortunately, the LMT did not participate in the 2021 campaign, limiting our ability to directly constrain the flux on scales ≲100 μas⁴. Therefore, constraining the absolute emission on ≲100 μas scales strongly depends on the image and gain priors. Taking this uncertainty into account, we report a conservative range for the compact flux of M87* in 2021 of F_com = [0.5 Jy,  0.9 Jy]. This constraint is similar to that reported in EHTC (2024b) for M87* in 2018. Note, the 2021 total flux constraint could be improved in the future by utilizing multifrequency information, for example, studies of the core M87* at 86 GHz similar to the analysis in EHTC (2024b).

6.2. Polarization results

Unlike total intensity, where significant evolution in the ring parameters was only observed for the PA, M87*’s linear polarization emission appeared distinct each year. Figure 6 shows how the linear polarization emission of M87* differs yearly in the EVPA pattern and the total linear polarization brightness. In 2017, a peak linear polarization fraction of ∼15% was found near the brightest region of the ring. Similarly, the absolute linear polarized brightness peaked in the south-western region of the ring. In contrast, in 2018 and 2021, the total intensity brightness maximum is de-polarized with a measured linear polarization fraction of ≲5%. Moreover, in 2018, the ring is almost entirely de-polarized, except for a single region in the western part of the image that has a linear polarization fraction of ∼5–10%. While the ring is more polarized in 2021, it is still measured to be less than the previous 2017 estimates, and the peak polarization fraction of the ring never exceeds 10% after blurring all methods to a common resolution of 20 μas.

Examining the image-integrated non-structural parameters in Fig. 7 and Table 2, a similar pattern is found. In 2017, using the methods described in this paper, we measured the image resolved fractional polarization to be ⟨|m|⟩ ∈ [5.7%,7.8%] on April 11, 2017, consistent with the 5.7–10.7% found in EHTC (2021a) after averaging over the four 2017 observations. This is significantly higher than the values we found for M87* 2018 (⟨|m|⟩ ∈ [2.0%,3.6%]) and 2021 (⟨|m|⟩ ∈ [3.0%,4.7%]). This result quantitatively demonstrates the relatively low polarization state of M87* in 2018 and 2021. The unresolved fractional linear polarization in 2017, |m_net|, while typically higher (1.0–3.7%), is consistent with the values found in 2018 (0.2–2.1%) and 2021 (0.3–1.6%). For the image-integrated EVPA, χ, we measured significant variability every year, where χ ∈ [ − 33° , −3° ] in 2017, [9° , 44° ] in 2018, and [ − 24° , 12° ] in 2021, reflecting the changing EVPA pattern in the image reconstructions.

Analysing the properties of the polarized ring, we found differences between the behaviour of β₁ and β₂ from 2017 to 2021. Interestingly, ∠β₁ was stable from 2017 to 2021, with ∠β₁ ranging from 71° to 171° for all years and methods, except for DoG-HIT in 2017. Similarly, |β₁| is consistent every year, although most methods found a larger β₁ in 2017 than in the other years. The exception to this is DoG-HIT. However, DoG-HIT struggled to recover |β₁| in the blinded synthetic data test (see Fig. 5); it was discrepant from the truth by 2–3% in an absolute sense in 2017 and 2018.

Unlike β₁, β₂ evolves substantially from 2017 to 2021 in both amplitude and phase. Using the methods in this paper, we found |β₂|∈[4.5%,5.5%] on April 11, 2017, consistent with the values in Table 2 taken from EHTC (2021a), which are averaged over the four observations in 2017. Similarly to ⟨|m|⟩, |β₂| decreases significantly in 2018 to 0.5–2.0%, but recovers slightly in 2021 to 1.8–3.9%. The observed de-polarization is also evident from the polarization calibration-insensitive closure traces, as noted above.

Unexpectedly, we found that ∠β₂ evolved substantially from 2017 to 2021. Using the methods in this paper, we found that ∠β₂ ∈ [ − 141° , −128° ] on April 11, 2017, which is consistent with the values reported in EHTC (2021a), with ∠β₂ ≈ [ − 163° , −127° ] averaged over all four days of observations in 2017. Although very little β₂ is measured in 2018, we found a similar value to that of 2017, although with greater uncertainty, [ − 156° , −99° ]. However, this consistency was broken in 2021. That is, we found that ∠β₂ rotated by about 60°  and flips sign, signaling a change in the polarization helicity. The most apparent visual changes in the EVPA patterns between 2017 and 2021 are in the dimmer northern half of the ring. While the EVPA patterns look more similar in the brighter south-west region, there is still an apparent shift in EVPA in the SW; since the overall image β₂ is intensity weighted Eq. (21), changes in the polarization structure in the brighter south-western region in fact account for most of the overall ≈60°  shift in ∠β₂ across the image. These results motivate further studies into the polarization structure in different parts of the emission ring and their relation to changes in the underlying magnetic field morphology or in the Faraday rotation along different lines of sight.

Note that from the synthetic data tests, we do not believe that the significant changes in ∠β₂ or the other polarimetric properties from 2017 to 2021 are due to changes in the coverage of M87*. The synthetic data results in Fig. 5 demonstrate that all algorithms recover the EVPA pattern to a high degree of certainty. Specifically, ∠β₂ is consistently one of the more robust quantities we measured in the synthetic data tests. Its values in 2021 are consistent on both April 13 and 18 (Fig. A.5) and for both reduction pipelines rPICARD and EHT-HOPS (Fig. A.3). Additionally, while the image-integrated polarization quantities ⟨|m|⟩, |β₁|, and |β₂| are biased low for some methods in the synthetic data tests, this bias decreases for the Bayesian methods in later years. Therefore, the fact that 2017 had a higher overall polarization is likely not a result of beam de-polarization or other instrumental effects. We discuss the interpretation of these changes in Sect. 7.

Finally, we mention that, while most polarized imaging algorithms include circular polarization maps, we do not report them in this paper. In Appendix B and Fig. A.2, we inspect the right and left circular polarization closure phases differences in 2017–2021 and found weak signals in the 2018 and 2021 observations. However, we found significant discrepancies in the circular polarization maps across methods due to residual instrumental systematics, such as right-left gain ratios. Like in past works (EHTC 2023), we could not robustly recover horizon-scale circular polarization structure for this reason.

6.3. Constraints on the extended non-ring emission

M87* has a parsec-scale jet that is consistently detected at longer radio wavelengths over multiple magnitudes of spatial scales (Kim et al. 2018; Lu et al. 2023; Cui et al. 2023; Walker et al. 2018; Lister et al. 2018; Nikonov et al. 2023), including the inner few hundred microarcseconds with 86 GHz GMVA observations (Kim et al. 2018; Lu et al. 2023; Kim et al. 2025). At 230 GHz, ALMA observations revealed the jet is pointed 288°  east of north on scales of ∼100 mas (Goddi et al. 2021). EHT observations in all three years found missing flux between the intra-site and the next shortest baselines (see e.g. the discussion in EHTC 2024b for details). Before 2021, the EHT lacked baselines ≳200 μas, limiting sensitivity to structures at jet-launching scales (e.g. Roelofs et al. 2020, 2023; Broderick et al. 2022; EHTC 2024b, 2019d). In 2021, the inclusion of KP and NOEMA added baselines beyond previous capabilities: NOEMA–PV with a baseline length of ∼1100 km probes scales of approximately ∼340 μas, and KP–SMT (∼100 km) probes spatial scales ∼2700 μas. To analyse the characteristics of the emission on these scales, we inspected three sets of closure phases: ALMA-APEX closures, ALMA-SMT-KP closures, and ALMA-PV-NOEMA, which we show in Fig. 8.

Fig. 8. Closure phases (band 4) on the three closure triangles that probe extended non-ring emission trivial closures – to AA-AX-PV (top), ALMA-PV-NOEMA (middle), and ALMA-SMT-KP (bottom) – with fits from the final image reconstructions for the types of imaging algorithms – Bayesian (Comrade), CLEAN (GPCAL), and the RML method (ehtim). Other methods yield qualitatively similar fits where available. Strategies to fit the extended non-ring emission to the 2021 closure phases are presented in Georgiev et al. (2025) and Saurabh et al. (2025).

We found that on ALMA-APEX triangles, the closure phases displayed a systematic negative offset of a few degrees. This offset is present in all years, but it is most notable in 2021. Focusing on the 2021 observation, the ALMA-SMT-KP triangle displayed a consistent offset of −6°  to −7° . Furthermore, the ALMA-PV-NOEMA triangle displayed a ∼5°  bump near 16 GMST. Interestingly, none of the 2021 reconstructions from the inner 100 μas are consistent with these closure phases when compared to the data under the assumption of small ≲1% systematic errors for these triangles.

Unfortunately, even in 2021, we found that directly imaging the structure of M87* on these scales with EHT data is still highly uncertain and model-dependent due to the sparse intermediate baseline coverage. A companion paper, Saurabh et al. (2025), explores the effects of including a Gaussian component with core ≲100 μas ring image on the ALMA-SMT-KP and ALMA-PV-NOEMA closure phases. They reported a preference for a faint (≤100 mJy) component located ∼300 μas to the west of the ring. However, the details of the locations are likely dependent on the morphology chosen for the extended non-ring emission.

While the addition of the Gaussian component improves the closure phase fits on the ALMA-SMT-KP and ALMA-PV-NOEMA triangles, the offset in the ALMA-APEX triangles still exists. If this offset is not due to instrumental effects, it implies the existence of non-trivial emission on scales ≳300 μas. This conclusion of extended emission is also supported by the measured compact flux of the ring and additional Gaussian component, which only accounts for 40–70% of the measured flux on the ALMA-APEX and JCMT-SMA baselines for 2017, 2018, and 2021. To constrain the emission on larger scales, Georgiev et al. (2025) analysed the offset in the ALMA-APEX closure phases. They found that the offset can be explained by extended emission whose centroid is located about 1 mas north-west of the ring. Encouragingly, the inferred centroid direction is consistent with the extended jet observed by ALMA (Goddi et al. 2021). This larger-scale emission does not impact the smaller-scale ALMA-KP-SMT or ALMA-PV-NOEMA triangles. Note that residual issues in the EHT data could be caused by systematic non-closing errors (see e.g. EHTC 2019c, for a list) that may impact the recovered properties of this extended ∼1 mas emission. We refer to Georgiev et al. (2025) for a thorough accounting of these systematics.

In conclusion, with the 2021 EHT array and the supporting analyses in Saurabh et al. (2025). and Georgiev et al. (2025) we found that emission on at least three spatial scales is necessary to fit M87* data: a ring-like structure for scales ≳1 Gλ, larger diffuse emission located about a few hundred microarcseconds from the core to match the KP-SMT and PV-NOEMA baselines, and an even larger component on ∼1 mas scales to explain the offsets on ALMA-APEX triangles. We anticipate that the exact nature of this emission will be imageable in newer EHT observations, such as the 2022 array, which regained the LMT.

7. Discussion

7.1. Total intensity

The general consistency between the images of M87* reconstructed from the 2017, 2018, and 2021 observations is as expected from theoretical models and GRMHD simulations. All datasets produce a ring-like structure of ∼43 μas diameter, despite differences in the EHT array configuration and the resulting (u, v) coverage. Nonetheless, well-constrained differences exist between the 2017 images and the subsequent ones; namely, an offset in the total intensity position angle and a higher source polarization in 2017. It is interesting to consider whether M87* was in a special state in 2017, characterized by the aforementioned differences.

In unresolved data, the 2018 observations show the most differences with other years, with a 3-day flare in very-high energy gamma-ray emission. Notably, there is no clear correlation between the EHT image or millimetre-flux variability and the very high-energy emission from M87*, which showed a ∼3 day flare in 2018 (M87 MWL2018). The similarity between the 2018 and 2021 horizon-scale images suggests that this flare originated downstream in the jet, as supported by simple, single-zone modelling. However, we detected no accompanying changes in the flux of the nearest knot (HST-1) or outer jets. On the other hand, more complex scenarios are also consistent with our observations, such as magnetic reconnection in a current sheet near the black hole. Reconnection could produce the high-energy emission without altering the millimetre properties (see Hakobyan et al. 2023; M87 MWL2018 for more discussion of these and other possibilities).

7.1.1. Ring diameter

With this work, we have further established the persistence of a ring-like emission feature on horizon scales in M87*. The stability of this feature further supports its interpretation as a black hole shadow from an optically thin accretion disk surrounding a Kerr black hole. After averaging the diameters from 2017, 2018, and 2021, we measured an average diameter of 43.9 ± 0.6 μas. Our work corroborates the findings from Wielgus et al. (2020) that measured consistent ring diameters between the 2017 EHT measurements and the ‘proto-EHT’ observations from 2009 to 2013. Additionally, given the stability of the ring, in a future work we will use the ring diameter measurements from the 2017, 2018, and 2021 observations to further improve the statistical error on the EHT’s mass-to-distance ratio (M/D) estimate of M87*.

We did not compare the EHT spectral index measurements to ALMA measurements for two reasons. First, we only considered the upper two EHT sidebands in this publication, and second, ALMA’s measurements of M87*’s spectral index have significant contributions from regions ≫200 μas. In future work, including new EHT observations with better (u, v) coverage, we will investigate the spectral evolution of M87* by computing spectral index maps between the 212–216 and 226–230 GHz receiver sidebands.

Finally, comparing our results to the 86 GHz ring diameter measurements in Lu et al. (2023) and Kim et al. (2025), the stability of the 230 GHz ring supports the view that the larger diameter measured at 86 GHz is due to synchrotron opacity effects. Moreover, future EHT analyses of the 260 GHz data recorded in 2024 and the 345 GHz data recorded in 2021, 2024, and 2025 will improve our understanding of the spectral nature of M87*’s ring emission.

7.1.2. Ring position angle

The position angle shift of the ring’s brightness peak between 2017 and 2018 is discussed in detail in EHTC (2025). As predicted in EHTC (2019e), η is expected to have a mean value around ∼203° –209° . This is about a 90°  offset from the black hole spin axis, which, in turn, we assume to be closely aligned with the 288°  ± 10°  large-scale jet direction (Walker et al. 2018; Cui et al. 2023). As the 2018 and 2021 measurements agree with this expectation and show a very similar source structure overall, it is imperative to understand the evolution and variability of the ring η from further EHT observing campaigns in 2022 and onwards. This variability, most likely related to the turbulent character of the flow, may be an informative observable for comparing theoretical models of M87* with observations (Wielgus et al. 2020; EHTC 2025). Comparing images from 2017, 2018, and 2021, it seems likely that M87*’s ring-like structure is not as variable as inferred in Wielgus et al. (2020), which may have been biased by the simple geometric model considered.

7.2. Polarization

The horizon-scale polarization measured by the EHT highly constrains numerical models of M87*, (EHTC 2021b, 2023) and Sgr A* (EHTC 2024a). Past EHT papers have extensively compared EHT images to model images generated from GRMHD simulations, ray-traced with general relativistic radiative transfer codes. While the total intensity image from the EHT’s 2017 observations only weakly constrained GRMHD models (EHTC 2019e), the addition of linear polarization constraints – particularly ⟨|m|⟩ and ∠β₂ – largely ruled out weak-magnetic field simulations in favour of MAD models as the preferred description of M87*’s accretion flow (EHTC 2021b). These results are consistent with the upper limits on the circular polarization fraction reported in EHTC (2023).

The qualitative differences in M87*’s polarized image in 2017, 2018, and 2021 naturally raise two questions: are MAD models still preferred in all three years of EHT observations? Can existing GRMHD simulations explain the observed changes in the polarized image from year to year?

We defer a comprehensive scoring analysis of GRMHD images to future work and focus on the potential implications of two major differences in the polarized observations from 2018 and 2021 compared to 2017: the ≈50% lower beam-scale polarization fraction ⟨|m|⟩ in both 2018 and 2021, and the ≈ − 60°  shift in ∠β₂ from 2017 to 2021, resulting in a change in the sign of ∠β₂ from negative to positive.

7.2.1. Changes in ⟨|m|⟩

Although both 2018 and 2021 EHT images of M87* are significantly more de-polarized than the 2017 image (⟨|m|⟩ ≈ 8% in 2017 vs ≈3 − 4% in 2018 and 2021), images from both years are still likely to be more consistent with MAD GRMHD models than weakly magnetized ‘standard and normal evolution’ (SANE) models. Although the observed changes in fractional polarization are significant, most SANE simulations are even more de-polarized than the 2018 and 2021 results would indicate. In particular, most SANE simulations have both ⟨|m|⟩ < 2.5% and |β₂|< 2.5% (EHTC 2021b, Figs. 7, 9). SANE disks are significantly more de-polarized than MAD disks, partly due to increased plasma turbulence but largely due to significantly higher Faraday depths through the disk from the higher density plasma needed (compared to MAD models) to produce the observed total flux density (Mościbrodzka et al. 2017). SANE models are also more likely to overproduce circular polarization and violate the EHT’s upper limits on the resolved circular polarization fraction ⟨|v|⟩ (Ricarte et al. 2021; EHTC 2023).

While SANE models are naturally de-polarized, most MAD models naturally produce more linear polarization than is observed in M87* (⟨|m|⟩ ≳ 10% after convolution with a 20 μas Gaussian). De-polarizing MAD models requires relatively cold electrons, with the EHT model scoring tentatively preferring an ion-to-electron temperature ratio R = T_i/T_e ≳ 40 to explain the low fractional linear polarization observed in 2017. Producing the decreased fractional linear polarization seen in 2018 and 2021 with MAD GRMHD models will require even colder electrons (larger R). While ideal GRMHD simulations like those used in EHTC (2019e, 2021b, 2023) require R to be assigned manually in post-processing, radiative two-temperature simulations can predict R by evolving separate ion and electron entropies self-consistently after adopting a sub-grid model for how the different species are heated. Recent surveys of two-temperature simulations found relatively low values of R ≈ 1 − 10 in the EHT emission region of GRMHD simulations (Mościbrodzka 2025; Chael 2025), with correspondingly large ⟨|m|⟩ ≈ 20 − 40%. The fact that EHT observations suggest that cold thermal electrons in the inner 5r_g are necessary to produce the observed low fractional polarization poses an intriguing tension with existing two-temperature models. This tension demonstrates that EHT observations can directly constrain plasmaheating mechanisms near the black hole.

7.2.2. Changes in ∠β₂

In the absence of Faraday rotation, the sign of ∠β₂ for face-on systems like M87* encodes the direction of electromagnetic energy flux; if magnetic fields are presumed to co-rotate with the plasma clockwise on the sky EHTC (2019e), the observed sign of ∠β₂ < 0°  in 2017 and 2018 is consistent with outward electromagnetic energy flux in analytic Blandford-Znajek monopole (Blandford & Znajek 1977) and GRMHD simulations (Chael et al. 2023). The observed ∠β₂ ∈ [161° , 166° ] in 2021 is not immediately consistent with either the measured value of ∠β₂ in 2017 and 2018 or this theoretical expectation.

The ≈ − 60°  shift in ∠β₂ from 2018 to 2021 could be caused by several different changes in the near-horizon emission region, including: (1) a change in the underlying magnetic field structure; (2) a change in the degree of Faraday rotation along the line of sight; and (3) evolving contributions from different emission regions along the line of sight (e.g. in the disk or jet); or some combination of all three.

As reported in Table 3, the measured Faraday rotation measure from unresolved ALMA observations, RM_ALMA, is largely consistent in magnitude across the three years of EHT observations reported here, though the rotation measure is known to vary in sign and magnitude on short timescales (e.g. changing sign from positive to negative between April 5 and 11, 2017). However, the unresolved ALMA observations are significantly affected by large-scale polarized emission from M87’s jet that does not contribute to EHT observations of the ∼40 μas scale core structure. Goddi et al. (2021) used a two-component Faraday screen model, constrained by ALMA and EHT observations, to attempt to disentangle the RM in the core and the extended jet; they found that the magnitude of the variable core RM in this model can significantly exceed the value measuredby ALMA.

A preliminary application of the two-component model to the 2017, 2018, and 2021 observations suggests an increase in the core Faraday rotation in 2021 sufficient to provide consistent de-rotated ∠β₂ values across all three years. However, comparing the de-rotated ∠β₂ values with simulation images requires careful accounting of both the potentially large uncertainties in the core rotation measure from the two-component screen model (Goddi et al. 2021) and the intrinsic internal rotation measure produced in GRMHD simulations (Ricarte et al. 2020). Therefore, we leave the detailed quantitative analysis to a future paper.

If changes in M87*’s Faraday screen are insufficient to explain the observed changes in ∠β₂, the likely alternative explanation is that the EHT is capturing changes in the emission region magnetic field or shifts in the location of the emission region(s) between 2017, 2018, and 2021. MAD GRMHD simulations have intrinsically variable ∠β₂ values (Palumbo et al. 2020), but the distributions of ∠β₂ for MAD simulations with a fixed black hole typically do not change sign when applying standard assumptions (e.g. electron heating prescriptions) in post-processing (EHTC 2021b; Chael et al. 2023). Further EHT monitoring of the polarization structure will indicate whether the observed variability in ∠β₂ is consistent with GRMHD variability, plasma propagation effects, or a large-scale change in the magnetic field topology between 2017 and 2021.

In follow-up works, we will further constrain the core Faraday rotation measure with the Goddi et al. (2021) two-component model and directly compare de-rotated ∠β₂ in the 2018 and 2021 EHT observations to GRMHD simulations, including uncertainties in the modelled core rotation measure. We will present a comprehensive comparison of the polarized EHT images in 2017, 2018, and 2021 to several GRMHD simulation libraries and discuss whether all three years of EHT observations still support the conclusion that M87* is in a MAD accretion state (EHTC 2021b, 2023). Moreover, we will discuss whether or not the observed variability in ∠β₂ after Faraday de-rotation is consistent with the MAD picture that ∠β₂ tracks field lines with outward Poynting flux and is sensitive to black hole spin (Palumbo et al. 2020; Chael et al. 2023).

8. Summary and conclusions

We have presented the first multi-epoch, polarimetric imaging of M87* on event-horizon scales using EHT observations from 2017, 2018, and 2021. This work includes the first results from the 2021 EHT observing campaign, which added the 12 m Kitt Peak Telescope and NOEMA, substantially improving the baseline coverage.

From the measured visibility data, we observed the same overall radial visibility amplitude plot across 2017, 2018, and 2021, with the 2021 data being of very high quality. For all three years, we found a visibility amplitude null at ∼3.4 Gλ and a second null at ∼8.3 Gλ, characteristic of a ring-like structure in the image. Analysing the polarized or cross-hand visibilities, we found that in 2018 and 2021, the data are de-polarized relative to 2017. Furthermore, we found that in 2018, there were no significant non-zero conjugate closure trace products, further evidence that M87* was largely de-polarized. Similar to 2017, we found small offsets in the differences between the right and left circular polarization closure phases in 2021. However, structural maps of the measured circular polarization are still dominated by instrumentalsystematics.

We employed seven distinct methods for polarimetric imaging and instrumental calibration, each validated via several synthetic data tests. When applied to EHT observations of M87*, each method produces images dominated by a ∼42 − 46 μas ring for each year, with a brightness maximum in the south. Despite this broad consistency, the images show variation from year to year in terms of the azimuthal brightness distribution of the ring. The 2017 images exhibit the largest differences from the two other years. Notably, the ring appears identical in 2018 and 2021, even though M87* had a recorded gamma-ray flare just prior to the 2018 EHT observation. Future EHT observations of M87*, including completed campaigns in 2022 and 2024 (see e.g. EHTC 2024), will further elucidate the typical state of its accretion flow and its turbulence statistics.

Unlike the total intensity, we found that the polarized structure of M87* changes dramatically every year. The ring observed in 2018 and 2021 is strongly de-polarized relative to 2017, with only a single region having a fractional linear polarization of ∼5–10%. Furthermore, the polarization pattern shows significant variations, with ∠β₂ rotating by ∼60°  and changing the overall EVPA helicity in 2021. These changes highlight the importance of the EHT multi-frequency observations at 345 GHz in 2021, 2024, and 2025 and at 260 GHz in 2024 to separate the internal polarization structure from external Faraday effects. The implications of these results on our understanding of the accretion flow around M87* will be analysed in a future publication and demonstrate the EHT’s capabilities to constrain plasma physics around SMBHs.

To clarify the location and process driving the 2018 very high-energy γ-ray flaring, a full theoretical study that simultaneously takes into account the changes in the EHT image flux and polarization, as well as the ring and inner jet PA between 2017–2021, is needed. However, such a study will be limited by the annual cadence of our observations, and the image quality in 2018 was overall worse than in 2017. The potential to link dynamics to particle acceleration in real time provides strong motivation for a high-cadence (every few days), longer-duration (over months)intensive monitoring campaign of M87* (EHTC 2024), together with multi-wavelength coverage, including the high-energy facilities. If another γ-ray flare is detected, such a campaign would drastically improve our chances of definitively identifying the physics driving at least one mechanism for very high-energy flaring in AGNs.

Finally, the 2021 data provide the first structural hints of M87*’s extended non-ring emission up to 1 mas thanks to the addition of the Kitt Peak and NOEMA telescopes. Using closure phases on triangles that include short EHT baselines, we demonstrate that the compact ring emission does not fully fit the EHT data. Specifically, all our imaging methods fail to fit the ALMA-PV-NOEMA and ALMA-SMT-KP closure phases when confined to ∼100 × 100 μas², showing residuals of the order of a few degrees. However, given the limited EHT coverage on these scales, we could not robustly image the extended emission. The measurements that probe the extended non-ring emission are best described by simple Gaussian components. A more involved modelling or imaging of the extended jet structure would overfit the data, particularly for pre-2021 EHT observations, when the important closure triangles were not present. Future EHT observations with improved intermediate baseline coverage, for example the 2022 EHT array, may enable the direct imaging of these extended structures at 230 GHz. These images would provide a direct view of the connection between the jet launch region and the black hole shadow (see e.g. EHTC 2024).

References

Appendix A: Data consistency and further checks

A.1. Coverage and systematic errors

Figure A.1 compares the M87* (u, v) coverage between the two 2021 observing days, April 13 and 18, analysed in this work. We focused on the April 18 data due to the substantially improved (u, v) coverage provided by NOEMA.

Fig. A.1. M87* (u, v) coverage of VLBI scans shown for the band 3 (227.1 GHz) data on April 13 and 18, 2021. The co-located ALMA, APEX, JCMT, and SMA stations are described by the Chile and Hawaii points, respectively. The labelling is done by antenna, using conjugate baselines to display the other station.

Table A.1 shows an overview of our systematic errors from observations over the years. 2021 has the best (u, v) coverage, allowing for more accurate calibration and overall the lowest systematics, as reflected by the data quality. Our sensitivity to non-zero trivial closure quantities is explored in Georgiev et al. (2025) and briefly discussed in Sect. 7.

A.2. Evidence of circular polarization

Figure A.2 shows the RR-LL closure phase differences, which are a robust estimate of the existence of circular polarization. Overall, smaller uncertainties in the 2021 data lead to small but visible non-zero Stokes 𝒱 signals on some closure triangles. Our polarized imaging tools pick up these biases, but residual gain ratio uncertainties are too significant to reconstruct a reliable circular polarization map from the data.

Fig. A.2. RR* − LL* closure phase differences that describe source circular polarization structure shown for two triangles in 2017, 2018, and 2021. The data are averaged over 12-minute-long segments.

A.3. 2021 HOPS comparison

In the main text, we only considered the rPICARD reduction for the 2021 data. The primary reason was that HOPS cannot currently correct for the NOEMA phase jumps noted in Appendix B.1. To remove the effects of the phase jumps, the EHT-HOPS pipeline instead selects the largest part of the data within each sub-band that exhibits a stable phase and flags the rest at the cost of sensitivity on NOEMA baselines.

To test whether the NOEMA differences and other differences between the EHT-HOPS and rPICARD pipelines could impact the results (see EHTC 2019c, for an explanation), we reimaged the M87* April 18, 2021, HOPS reduction using the Comrade pipeline and compared it with the results presented in the main text. Figure A.3 shows the image reconstructions from the rPICARD and EHT-HOPS pipelines. Likely due to residual issues in the HOPS NOEMA calibration, we found a large amplitude offset in the NOEMA gain. Due to this offset, we measured a smaller compact flux for the HOPS data. To fix this, we renormalized the total flux to match the total flux from the fiducial 2021 image. After this adjustment, we found very similar image reconstructions between the rPICARD and EHT-HOPS pipelines. This result demonstrates the robustness of our results to different calibration pipelines.

Fig. A.3. Image reconstructions averaged over bands 3 and 4 using the Comrade pipeline and the April 18, 2021, reduction from the standard rPICARD pipeline used in this paper (left) and the EHT-HOPS pipeline (right). Both images have been normalized to have the same compact flux due to residual issues with the EHT-HOPS calibration of NOEMA.

A.4. Cross-check of the 2021 polarized imaging results from April 13

In the main text, we focus on the results from April 18, 2021. This day was selected because it was the only track in 2021 where the entire 2021 EHT array observed M87* for a significant period. The next best day is April 13, 2021, in terms of coverage; however, NOEMA did not participate due to bad weather. As a result, April 13, 2021, lacks any baselines on the scales of ∼100 − 500 μas, weakening the constraining power of the EHT on intermediate scales. Given the issues with NOEMA in the 2021 observation described in Appendix B, April 13 provides an additional check that NOEMA is not dramatically altering the conclusions of April 18. Additionally, since the M87* observations are 5 days apart, this allowed us to test for any significant evolution in the source, which is expected to be minimal given the dynamical timescale of M87* (EHTC 2019d,e; Georgiev et al. 2022).

Figure A.4 shows the fiducial image reconstruction produced by averaging over the 7 methods in Sect. 3 on April 13, 2021 (top row) and April 18, 2021 (bottom row). From Fig. A.4, we see that the two images display a very similar EVPA pattern and that the overall fractional polarization appears qualitatively similar. Both images indicate that the ring is de-polarized in the brightest region, and the polarized emission is concentrated in the east-west portion of the ring. The similarity in the EVPA pattern is noteworthy given that the ALMA rotation measure differed by a factor of 2 (see Table 3) between April 13 and April 18.

Fig. A.4. Fiducial images from April 13, 2021, and April 18, 2021. The April 13 image is constructed in the same manner as Fig. 6.

Looking at the images on April 13 and 18 quantitatively, Fig. A.5 shows the estimated parameters for both days. After averaging over all methods, the grey bands indicated that the total intensity image parameters were consistent across the two days. We see a similar result for the linear polarization parameters, except for |β₁|, where there was a > 1σ difference in the method average, although a |β₁|≈1% is within the 2σ error bars. The measurement of ∠β₁ also showed some discrepancies between April 13 and 18, although again it is not significant. Finally, we also found that the leakage estimates on April 13 agree with those on April 18 across all methods and sites.

Fig. A.5. Comparison of imaging parameter estimation results from April 13 and 18, 2021. The conventions follow Fig. 7.

Appendix B: Data issues and mitigation strategies

B.1. NOEMA instrumental phase calibration

NOEMA experienced phase jumps at the edges of each of its 64 MHz sub-bands. Each jump amounted to a multiple of 90° . The phase jumps are mostly constant in time, but change in value when the NOEMA PolyFiX backend is reset during observations (see Piétu et al. 2025, for a description of NOEMA’s phasing system). During the April 18 observations, the backends were reset twice, leading to multiple changes in the instrumental phases. In addition, the EHT data are correlated using the 58 MHz sub-band outputbands mode, which results in phase jumps within the final EHT sub-bands. The CASA FRINGEFIT task used by rPICARD cannot fit for phase jumps within sub-bands (spectral windows in the CASA Measurement Set). Thus, we calibrated these phase jumps in a dedicated calibration step.

For the rPICARD-calibrated data, we extended our pipeline with a dedicated NOEMA instrumental phase calibration step. We re-grid the data virtually by defining visibility chunks in the native 64 MHz gridding used by the NOEMA backend. These are then fit separately for each scan and NOEMA frequency window. To increase robustness, we required that the determined phase jumps are multiples of 90° . We also require that the phase jumps are consistent for more than 5 scans, to ensure that only real changes due to the infrequent resets of the PolyFiX backend are accepted.

The current HOPS processing environment does not support non-linear phase bandpass correction within a 58 MHz sub-band at this moment⁵. To mitigate de-coherence from the uncorrected phase jump, the EHT-HOPS pipeline flags the smaller portion of each sub-band before or after the jump, which allows the remaining bandpass to be calibrated using standard corrections. This process removes ∼25% of the total integrated bandwidth for NOEMA baselines, reducing overall sensitivity.

In late 2023, the origin of the phase jumps was identified in the signal processing firmware of the PolyFiX correlator. A PolyFiX firmware update was implemented to remove these phase offsets in the NOEMA passband. The new firmware version was installed at the NOEMA observatory in January 2024 and has been successfully tested on-sky. Consequently, the dedicated phase calibration of the NOEMA visibilities will not be required for EHT data from 2024 onwards. Further details on the rPICARD calibration steps and the PolyFiX firmware update are described in von Fellenberg et al. (2025).

B.2. JCMT leakage

In the 2018 data, we infer that JCMT had an unusually large polarization leakage of ∼ 30 %. We traced the likely origin to a 17°  misalignment of the quarter-wave plate (QWP). The QWP is aligned regularly using data from a baseline between the JCMT and a single SMA antenna. The optimal angle is found by nulling the signal, i.e. minimizing the cross-correlation signal for a calibrator target (3C 84 was used in 2018). From this position, where JCMT and SMA record in opposite polarizations, the QWP is rotated back by 90°  to align the JCMT polarization with the SMA. In 2018, this procedure was complicated by noisy data from poor weather conditions, making QWP alignment difficult. To remove the impact of the large JCMT leakage in 2018, all imaging methods flagged JCMT baselines.

B.3. 2021 specific calibration steps

Three stations required dedicated calibration steps in 2021. PV experienced a change in instrumental R-L delay at 00:37 UT on April 18. To account for this, we allowed for dedicated, per-VLBI-scan instrumental delay calibration of PV. KP showed a small drift in the L-R phase over time, which was corrected through self-calibration during imaging for the Comrade pipeline.

On April 13, APEX erroneously executed observations with the activated ‘wobbler’ for the first ≈2 h, which resulted in the antenna pointing off-source for roughly 50% of the time, leading to a correspondingly halved visibility amplitude, without affecting the phase stability. The resultant time-dependent amplitude gain error led to poor network calibration. We thus decided to calibrate the resulting APEX gain error before post-processing by fitting a Gaussian component to the data in Difmap. We then amplitude self-calibrate the APEX station, keeping all other stations fixed. The resulting corrected data are exported and used for further post-processing (i.e. network calibration).

B.4. GLT system temperatures on April 13

The GLT did not record system temperatures on April 13. Using the methods outlined in Section 5.2.4 of Koay et al. (2023b), we instead used modelled T_sys values, which resulted in large GLT amplitude gain errors that need to be corrected through self-calibration.

Appendix C: Details of the polarized imaging methods

C.1. CLEAN methods

C.1.1. Difmap/GPCAL

We conducted imaging with CLEAN using Difmap (Shepherd 1997, 2011) based on the imaging pipeline utilized for the 2018 EHT images of M87* (EHTC 2024b). In summary, the pipeline first searches for the best geometric model among a circular Gaussian component with 15 μas FWHM (representing an unresolved symmetric model), a uniform disk with sizes ranging from 56 to 84 μas in steps of 4 μas, and a uniform ring with sizes ranging from 36 to 68 μas (also in steps of 4 μas, without width). The best model was selected based on the closure phase normalized χ², and phase-only self-calibration was conducted using the model. The pipeline surveys five parameters: the total assumed compact flux density, cleaning stopping condition, relative weight correction factor for ALMA in self-calibration, diameter of the CLEAN window, and the power-law scaling of the (u, v) density weighting function. The same parameter ranges used for the 2018 M87* imaging were used. The final representative image was selected based on the lowest overall χ_CPh² and χ_LCA². We did not perform the top-set image selection procedure used in the 2017 and 2018 EHT imaging of M87*.

We utilized GPCAL (Park et al. 2021, 2023a,b), a novel instrumental polarization calibration pipeline based on AIPS and Difmap, to correct the polarimetric leakages in the visibilities. The pipeline first derives the leakages of the Chile and Hawaii intra-site stations, ALMA, APEX, SMA, and JCMT, using the visibilities of the intra-site baselines (ALMA-APEX and SMA-JCMT) only. In the derivation, the source structure was assumed to be point-like. Then, the pipeline derives the leakages of the other stations by utilizing all the baselines except for the intra-site baselines and by fixing the leakages of the intra-site stations during the fitting procedure.

Linear polarization images were obtained with CLEAN using Difmap and the leakage-corrected visibilities, as done in the 2017 polarimetric imaging of M87* (EHTC 2021a). CLEAN was used for Stokes 𝒬 and 𝒰 images using the CLEAN windows utilized for the Stokes ℐ imaging, until the peak in the residual dirty image within the windows reached the RMS noise level of the dirty image.

C.1.2. Difmap/LPCAL

For the determination of the instrumental polarization of the antennas (so-called D terms), we also utilized the traditional AIPS-based task LPCAL (Leppanen et al. 1995), which has been used extensively in VLBI polarimetry during the past two decades. The total intensity image (Stokes I) was performed with Difmap in the same way as described in Appendix C.1.1 but with different hyperparameters, for example the number of iterations and data flags. The edited and fully self-calibrated data were then imported to AIPS, where the AIPS task CALIB was used to remove the gain ratios between the parallel-hand polarizations (RR^* and LL^*). We then applied the AIPS task CCEDT to split the δ components of the Stokes I image into up to 10 compact sub-regions of constant fractional polarization, which are used as input models for the AIPS task LPCAL. The antenna D-terms derived from LPCAL were then applied to the data using another round of amplitude and phase self-calibration with the option DOPOL in the CALIB AIPS task. Finally, the data are exported from AIPS back to a local data file, from which the final linear polarization images are obtained in the same manner as described in the last paragraph of the previous section.

C.2. Regularized maximum likelihood methods

Regularized maximum likelihood techniques balance data fidelity and regularization assumptions on the image, such as smoothness, sparsity, or entropy. This is achieved by minimizing a weighted sum of terms measuring data fit quality, and terms measuring its feasibility:

$$\begin{array}{c}\hfill \widehat{I}\in {\mathrm{argmin}}_I{\displaystyle \sum _d}{\alpha }_d{\chi }_d^2(I)+{\displaystyle \sum _x}{\beta }_x{R}_x(I).\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}(22\text{revisited})\end{array}$$

Throughout this manuscript, three algorithms are applied that utilize this concept. These are ehtim (Chael et al. 2016, 2018), DoG-HIT (Müller & Lobanov 2022, 2023a), and MOEA/D (Müller et al. 2023; Mus et al. 2024a). While all these algorithms apply similar heuristics, they differ significantly in the regularization assumptions applied and in the minimization problem in Eq. (22) is solved. A comprehensive overview is presented in Table C.1.

In particular, ehtim has been used intensively and successfully in previous EHT data analyses of black hole shadows in total intensity (EHTC 2019d, 2022, 2024b) and polarization (EHTC 2021a, 2024c). The robustness of and confidence in the ehtim imaging stemmed partly from thorough and detailed parameter surveys and studies on synthetic data. We did not perform the full evaluation and parameter search in this manuscript. However, we performed spot checks with hyperparameter configurations from previous years and accompanied the RML imaging procedure with two alternative approaches: DoG-HIT and MOEA/D. DoG-HIT is a lightweight compressive sensing approach that, due to its simple design, the small number of hyperparameters, and the data-driven selection of the regularization assumption (ie, the exact form of the sparsifying basis), is a robust and fast alternative well suited for the EHT. MOEA/D replaces parameter surveys on synthetic data with a multi-objective exploration of the optimization landscape based on the data itself.

C.2.1. DoG-HIT

DoG-HIT is a compressive sensing algorithm that models the image by wavelets (Müller & Lobanov 2022). In contrast to traditional compressive sensing approaches, the wavelets have been specifically designed to fit to the (u, v) coverage (Müller & Lobanov 2023b). This is done to offer an ideal separation between covered and non-covered Fourier coefficients, highlighting the features introduced by the former, and suppressing the latter, motivated by the exceptional sparse coverage of the EHT. All in all, DoG-HIT solves the optimization problem (Müller & Lobanov 2022):

$$\begin{array}{cc}\hfill & \widehat{\omega }\in {\mathrm{argmin}}_{\omega }\{{\chi }_{\mathrm{CPh}}^2(\mathrm{\Psi }\omega )+{\chi }_{\mathrm{LCA}}^2(\mathrm{\Psi }\omega )+{||\omega ||}_{l^0}+R_{\mathrm{flux}}(\mathrm{\Psi }\omega )\},\hfill \\ \hfill & \widehat{\mathcal{I}}=\mathrm{\Psi }\widehat{\omega }.\hfill \end{array}$$(C.1)

Here, Ψ denotes the wavelet dictionary, and ω the wavelet coefficients.

Due to its sparsifying approach, DoG-HIT computes the multi-resolution support (the spatial scales and positions of all wavelet coefficients that are statistically significant to represent the features of the image) as a byproduct of the total intensity imaging. Müller & Lobanov (2023a) noted that the multi-resolution support provides a powerful prior information for the reconstruction of channelized datasets (e.g. spectral channels, polarization channels, or time-dynamic reconstructions). Linear polarization and circular polarization were recovered by constrained minimization, i.e. we minimized χ_pvis² and χ_CPh², respectively, and only allowed coefficients in the multi-resolution support to vary. Since the form of the basis functions is rigidly set by the instrumental configuration and there is no manual selection, the number of free hyperparameters and the human bias in imaging are reduced.

DoG-HIT handles total intensity, linear polarization, and circular polarization separately. First, we recover an image in total intensity by closure-only imaging. In a second step, we calibrate the gains and leakages, and recover linear polarization with the above-mentioned strategy of only changing the parameters in the multi-resolution support. Finally, we recover circular polarization using a gradient descent approach.

C.2.2. ehtim

ehtim directly implements the minimization problem Eq. (22), with objective functionals summarized in Table C.1. Reconstruction with ehtim depends strongly on the regularization hyperparameters α_d, β_x, which balance the data fidelity and regularization assumptions. These parameters have been selected by parameter surveys in previous studies (e.g. EHTC 2019d), i.e. various parameter combinations have been tested against a variety of geometric models. ehtim can process visibilities in total intensity, linear, and circular polarization simultaneously, by adding all the objectives together. However, the combined parameter space would be very high-dimensional. Instead, we followed the approach taken in EHTC (2021a) and EHTC (2024c). We recovered total intensity first, then fixed the total intensity structure, and recovered linear and circular polarization afterwards.

Reflecting this splitting strategy, the parameter search has also been split into total intensity and polarization. A parameter survey has been recomputed in total intensity on the four original geometric models (double, disk, ring, and crescent) used to validate the first release of an image of the black hole shadow (EHTC 2019d). Additionally, the data and regularization terms for linear polarization and circular polarization have been surveyed for a number of crescent images with varying asymmetry, orientation of the EVPA pattern, and polarization fraction as a spot-check. The best parameter combination (in terms of cross-correlation to the ground truth) matches the one reported in EHTC (2021a) and has therefore been adapted for the three years for the sake of consistency. Our final weights are shown in Table C.1. Ultimately, the validation for this configuration results from the validation described in Sect. 5.

The full polarized imaging and calibration pipeline consists first of imaging the source in total intensity. The reconstruction is iteratively blurred, and the reconstruction is redone to avoid getting trapped in local minima. The observation is scan-averaged, and imaging optimization alternates with self-calibration steps iteratively. During the initial self-calibration steps, we assume that the gains for the right- and left-handed polarization feeds are identical. As a second step, we fix the total intensity image and then recover the linear polarization image. The linear polarization map is recovered in a loop, alternating between image reconstruction and leakage estimation. Finally, we self-calibrate the data again, allowing for different gains in the two polarization feeds, and solve for the circular polarization image.

C.2.3. MOEA/D

MOEA/D is a widely used optimization framework for multi-objective optimization originally proposed by Zhang & Li (2007) and Li & Zhang (2009). Recently, it has been transferred to VLBI imaging (Müller et al. 2023), and expanded to include polarization (Mus et al. 2024a). RML methods proceed by minimizing a weighted sum that balances data fidelity terms and regularizers against each other. MOEA/D, aims to find the ‘best-compromise’ solutions between the different regularizers. The notion of Pareto optimality identifies this: For a multi-objective vector of functionals [χ₁², χ₂², ..., R₁, R₂, ...], a solution is called ‘Pareto optimal’ if the further optimization along one of the functionals automatically has to worsen the scoring in another.

The set of all Pareto-optimal image structures is known as the Pareto front. MOEA/D approximates this front. Compared to ehtim, the solution to the optimization problem defined by MOEA/D is not a single image, but a space of images (the Pareto front), representing the best-compromise solutions among all objectives. Minimization is achieved through genetic evolution, where every genome in every generation represents a single image (Müller et al. 2023). The final population approximates the Pareto front.

It has been demonstrated both for total intensity (Müller et al. 2023) and polarization (Müller 2024) that the Pareto front of the calibration-independent structure separates into multiple disjoint clusters of solutions, interpreted as local minima of a potentially multimodal optimization landscape. In this work, we report the image closest to the utopian ideal point following the strategy outlined in Müller et al. (2023); Mus et al. (2024b), and leave the detailed analysis of the multi-modality of the problem to a subsequent work.

MOEA/D has been implemented in the MrBeam (Müller & Lobanov 2022; Müller et al. 2023; Mus et al. 2024b) software package, which calls ehtim-subroutines for the calibration and evaluation of the data terms. The full imaging and calibration procedure is therefore similar to the one adapted for ehtim, although not equivalent due to the difference in the surveying strategy (i.e. no survey on synthetic data are needed), the different concept of optimality, and the optimization procedure. The data fidelity and penalty terms are summarized in Table C.1.

Due to this proximity, MOEA/D’s imaging pipeline resembles the ehtim strategy. We first recovered the total intensity image by fitting only closure quantities. To help MOEA/D converge, we derived a starting point for the population using an un-regularized run of ehtim. Afterwards, we selected a representative cluster of solutions (proximity to the utopian) and performed self-calibration and leakage calibration. Then, we ran MOEA/D again to recover the linear and circular polarization images. In contrast to ehtim, we solved for linear polarization and circular polarization at the same time.

C.3. Bayesian methods

This section describes the two Bayesian polarized imaging algorithms, THEMIS and Comrade, used for image reconstruction, residual calibration, and leakage corrections. Both Bayesian codes aim to simultaneously solve for all four Stokes parameters in the image and calibration and leakage terms through forward modelling of the measurement process.

Both THEMIS and Comrade utilize a Cartesian grid of fluxes, where for each pixel, the Stokes vector is parameterized using the Poincaré representation,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\mathcal{I}}_{\mathit{ij}}\\ {\mathcal{Q}}_{\mathit{ij}}\\ {\mathcal{U}}_{\mathit{ij}}\\ {\mathcal{V}}_{\mathit{ij}}\end{array}\right)={\mathcal{I}}_{\mathit{ij}}\left(\begin{array}{c}1\\ p_{\mathit{ij}}{\widehat{\mathit{p}}}_{ij,x}\\ p_{\mathit{ij}}{\widehat{\mathit{p}}}_{ij,y}\\ p_{\mathit{ij}}{\widehat{\mathit{p}}}_{ij,z},\end{array}\right),\end{array}$$(C.2)

where p_ij is the total polarization fraction and $\widehat{\mathbf{p}}$ is a unit vector in ℝ³, which are all included as parameters in the forward model. Each pixel is modelled as a point source and then convolved with a kernel to create a continuous image function. The ideal visibilities are calculated using Eq. 1 based on the Stokes Images produced from each method. Site-based instrumental corruptions are applied using a RIME formalism, utilizing the Jones matrices for feed rotation, leakage, and residual gains from Eq. 4. The model visibilities are then compared to the data using a complex Gaussian likelihood, which is calculated as the product of all individual likelihoods for all parallel and cross-hand products. We now describe the individual details for each method.

C.3.1. Comrade

Comrade (Tiede 2022) is a Bayesian polarized imaging and calibration software suite written in the Julia programming language (Bezanson et al. 2017). Here, we briefly describe the polarized imaging process and the priors used; a complete explanation can be found in Tiede (2022) and Tiede et al. (in prep.).

To prepare the data, Comrade first coherently averaged the data over observation scans using the Pyehtim.jl scan_average function and then added 2% systematic noise to model residual calibration errors, such as loss of coherence due to time and frequency averaging. Additionally, due to the over-resolved flux on short baselines, Comrade flagged the ALMA-APEX, JCMT-SMA, and KP-SMT baselines, and fit the total flux of the compact emission.

Comrade’s image model is aligned with the equatorial coordinate axes and uses a 64 × 64 raster with a field of view of 200 μas. For the total intensity image raster ℐ, we use a first-order Gaussian Markov random field (GMRF) prior on the log-ratio transformed pixel fluxes r_ij, which are related to the total intensity fluxes ℐ_ij by

$$\begin{array}{c}\hfill {\mathcal{I}}_{\mathit{ij}}=F_0\frac{e^{r_{\mathit{ij}}}}{{\displaystyle {\sum }_{\mathit{ij}}}{e}^{r_{\mathit{ij}}}},\end{array}$$(C.3)

where F₀ is the total flux of the raster component. For the GMRF, the variance and correlation length of the random field are included as hyperparameters. For more information, see Tiede et al. (in prep.). The GMRF acts as multiplicative fluctuations in some mean image structures, which we assume to be given by a 60 μas Gaussian blob. However, during testing, we found that the size of the Gaussian did not appreciably change the resulting images.

The total fractional polarization amplitude, p_ij, is given by

$$\begin{array}{c}\hfill p_{\mathit{ij}}={[1+exp(\overline{\ell }+{\sigma }_p{\ell }_{\mathit{ij}})]}^{-1},\end{array}$$(C.4)

where $\overline{\ell }$ is the average logistic total fractional polarization, σ_p is the variance of the logistic fractional polarization, and ℓ_ij are the scaled fluctuations of the total fractional polarization, whose prior is another independent GMRF. This corresponds to modelling the total fractional polarization in logit space, with some constant mean fractional polarization and some variance. This parameterization ensures that 0 ≤ p_ij ≤ 1 for all values. For the total polarization direction on the Poincaré sphere, for each unit vector ${\widehat{p}}_{\mathit{ij}}$, we used the uniform distribution on the sphere by modelling each direction with an independent unit normal distribution and then dividing by the total length to ensure that the vector is normalized.

Gain and gain ratios are allowed to vary every scan, whereas the D-terms are fixed to constant values over each track. For each station, we use a Gaussian prior on the log-gain amplitudes with a mean of 0 and a standard deviation of 0.2, except for GLT, which uses a standard deviation of 1.0 to model the large-gain-amplitude fluctuations. For the gain phases and gain phase ratios, first lock the gain phase and gain phase ratio for ALMA to zero to model the a priori EVPA calibration included in the PolConvert procedure. If ALMA is absent in the scan, we alphabetically select a different reference site and set only the right-hand gain phase to zero while keeping the gain ratio phase as a model parameter. For the gain phases, we use a zero-mean von Mises prior with concentration parameter π⁻² that essentially creates a uniform prior in the interval [ − π, π]. For the gain ratio phases, we use a zero-mean von Mises prior with a concentration parameter 0.5² to model the a priori calibration that the rPICARD pipeline does to remove gain phase ratio.

During imaging, Comrade simultaneously models all four Stokes parameters and the instrumental terms. Since Comrade is a Bayesian imaging algorithm, it estimates the uncertainty for all model parameters, including the image and instrumental terms. Its output is a set of samples approximately drawn from the posterior using Markov chain Monte Carlo (MCMC). Specifically, Comrade uses the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2011) in the Julia sampling package AdvancedHMC.jl, utilizing the automatic differentiation package Enzyme.jl (Moses & Churavy 2020; Moses et al. 2021) to compute gradients. To initialize the chain, we ran ten rounds of optimization to find the approximate posterior maximum. From this location, we then ran NUTS for 10,000 adaptation steps to tune the sampler, followed by 10,000 sampling steps. From these 10,000 samples, we randomly selected 500 images for further analysis.

C.3.2. THEMIS

THEMIS provides a Bayesian framework for modelling and parameter estimation from EHT data (Broderick et al. 2020a). Imaging with THEMIS employs a model consisting of a rectilinear set of control points spanned via a bicubic spline to represent the intensity map, the polarization fraction map, EVPA map, and the Stokes V fraction (Broderick et al. 2020b; EHTC 2021a). The orientation of the raster and the field of view are free parameters, which allow the raster to expand, contract, and rotate to maximize the model’s freedom to fit the data. The resolution of the raster is typically determined by maximizing the Bayesian evidence on the raster dimension; this is small due to the limited number of EHT resolution elements across M87*. To facilitate cross-epoch comparisons, here we adopted a 5x5 raster based on the polarized imaging study of the 2017 observations (EHTC 2021a). Finally, a large-scale asymmetric Gaussian was included to capture potential discrepancies between the intra-site and VLBI baselines.

THEMIS reconstructions are fit directly to each band’s scan-averaged complex visibilities (RR^*, LL^*, RL^*, LR^*) independently. Site-specific leakages and complex gains are fit simultaneously with image generation as described in EHTC (2021a). THEMIS assumes that the right and left-hand gains are equal; a 3% systematic error is added to partially mitigate non-unity gain ratios.

The result of the THEMIS fits is an approximate posterior composed of a set of images used for Bayesian uncertainty quantification. THEMIS provides a number of posterior sampling methods, for which the most common output is a MCMC chain that supports subsequent Bayesian interpretation. To ensure efficient sampling of the posterior, we use the differential even-odd parallel tempering scheme, with each tempering level explored via the Hamiltonian Monte Carlo NUTS algorithm implemented by the Stan package (Syed et al. 2019; Carpenter et al. 2017). This sampler has been shown to effectively capture multimodal posteriors (see e.g. EHTC 2021a). Chain convergence is assessed by visual inspection of parameter traces and quantitative chain statistics, including integrated autocorrelation time, split-$\widehat{R}$, and parameter rank distributions (Vehtari et al. 2021), and typically requires ∼10⁵ MCMC steps. The number of tempering levels is chosen to ensure efficient communication between the highest- and lowest-temperature levels, typically set at 80 due to the complicated nature of the model.

Appendix D: Synthetic data validation

D.1. Models

Here, we describe the models used to generate the synthetic data. For the blinded GRMHD synthetic data, a KHARMA MAD simulation (EHTC 2025) with inclination i = 17° , spin a = −0.5, R_low = 10, and R_high = 40 was used. For each epoch, 2017, 2018, and 2021, a random snapshot was chosen with parameters that did not necessarily match the properties of the measured data for each epoch. The selected GRMHD snapshots are shown in the first column of Fig. D.1.

To test the impact of over-resolved polarized flux on ring reconstructions, we constructed a pair of geometric models consisting of an m-ring and an extended jet component. The jet consists of three polarized, elliptical, Gaussian components — this is the same tri-Gaussian jet model used in the synthetic data generation of (EHTC 2023), only polarized. The amount of polarized flux in the jet was constructed to roughly match the upper and lower levels of polarization observed on the short baseline ALMA-APEX, JCMT-SMA, and KP-SMT.

D.2. Blinded GRMHD results: Individual methods

The blinded synthetic data image reconstructions for each method are shown in Fig. D.1 with image reconstruction blurred to match the resolution of the truth convolved with a 20 μas circular Gaussian beam. For 2017, we see the most diversity across the methods as expected due to the overall poorer coverage of the 2017 data. This is similar to the diversity between methods found in EHTC (2021a). In 2018 and 2021, we found that once blurred to the same resolution, the image reconstructions show more consistency between methods, especially in the brighter regions of the ring.

Fig. D.1. Polarized imaging results from the blinded GRMHD synthetic data test. The leftmost column shows the ground truth GRMHD image blurred by a 20 μas Gaussian beam; the other columns show the imaging results from each method blurred to an equivalent resolution. From top to bottom, the rows correspond to the 2017, 2018, and 2021 results.

Figure D.2 shows the recovered leakage compared to the true values used when creating the synthetic data. The rows show the recovered leakage for each epoch of synthetic data. Generally, the two Bayesian methods produce the best leakage estimates each year and provide reliable estimates of the true leakage to within a few percent for all years.

Fig. D.2. D-term recovery results of blinded GRMHD synthetic data from each polarized imaging method (columns) and year (rows). The x-axis of each plot is the true leakage, and the y-axis is the measured leakage. Perfect recovery means that all points should live on the y = x line. The total L2 difference between the ground truth and recovered leakage averaged over all stations is shown in the bottom right of each plot. For the Bayesian methods, we also show a χ2 value computed as the square difference between the true value and the recovered leakage divided by the variance of the posterior estimate. For the Bayesian methods, the points denote the posterior median leakage, and error bars show the 95% credible interval about the median.

For the non-Bayesian methods, there is more scatter in the estimated leakage, but the error in the leakage estimates improves each year due to the increased number of stations in the array, making it easier to solve for leakage. In 2018 and 2021, we found that GLT’s leakage estimates tended to have the largest discrepancies from the truth. This error is expected given the relatively poor parallactic angle coverage ( < 15° ) of the site. Encouragingly, the Bayesian methods can still constrain GLT’s leakage, with the posterior for Comrade and THEMIS containing the truth within their 95% credible interval for each station.

D.3. Impact of polarized extended emission on the core reconstructions

The impact of polarized extended emission for each imaging pipeline is shown in Fig. D.3. The overall reconstruction quality of the low-polarized jet synthetic data is shown in the left column of Fig. D.3. For total intensity (top row), we found that all methods produced reliable estimates of the total intensity structure across all epochs with ρ_ℐ > 0.985. For linear polarization, we found more mixed results. Both Bayesian methods THEMIS and Comrade performed well every year, finding a ρ_𝒫 ≳ 0.95 each year. In addition, ehtim performed well for the low-polarization jet reconstructions. DoG-HIT performed well in 2017 and 2021, but its total intensity image got noticeably worse in 2018; however, the reconstruction was still of high quality overall. In 2018, GPCAL and MOEA/D struggled to recover the true linear polarized structure but recovered in 2021, suggesting it may have been a convergence issue in the algorithms.

Fig. D.3. Cross-correlation between the low polarized and high polarized extended jet geometric tests. The ring model for the two tests is identical, meaning that the cross-correlation in the core region would be equal to unity for perfect reconstructions. Left column: Total intensity (top row) and linear polarization cross-correlation (bottom row) between the low-polarized jet reconstructions. Right column: Same quantities but with the cross-correlation computed between the low and high polarized jet reconstructions.

To compare the reconstructions of the same compact features in the low- and high-polarized jet synthetic data, we then cross-correlated the reconstructions from each pipeline, which is shown in the right column of Fig. D.3. Note that we cross-correlated the mean images for the Bayesian methods to remove the impact of different thermal noise realizations on the image reconstructions. The total intensity cross correlation between the low and high polarized jet reconstructions was substantially closer to unity ( > 0.998) for all methods than the cross correlation between the low polarized reconstructions and the truth. A similar result was found for the linear polarization cross-correlation between the two synthetic datasets, albeit with a few exceptions. Namely, for 5/7 methods, we found a linear polarization cross-correlation of > 0.98, and this consistency is always greater than the consistency between the reconstruction and ground truth. The only exceptions to this are DoG-HIT. For DoG-HIT, the low-polarized jet reconstruction showed a polarized cross-correlation of almost unity (0.96). In contrast, the high-polarized jet reconstruction performed worse (0.88), suggesting that the method is weakly sensitive to extended emission. The Comrade, THEMIS, and ehtim are not affected by different levels of extended emission when reconstructing the compact ring structure.

Appendix E: Leakage estimates in 2018 and 2021

Individual leakage solutions are shown in Fig. E.1 for 2018 and 2021. For the 2017 leakage estimates, we refer the reader to EHTC (2021a) for a multisource analysis. Note that in general, the 2017 leakage estimates we found in this paper are consistent with those found in EHTC (2021a). We found that the leakage estimates produced by each pipeline are consistent to within 5%, with a few notable exceptions. For the 2018 reconstructions, GLT leakage estimates show significant discrepancies between methods. GLT’s leakage, being difficult to estimate, was also observed in the synthetic data tests, where non-Bayesian methods could not measure GLT’s leakage to better than 10%. The discrepancy in the estimated leakage of GLT is not surprising due to its limited parallactic angle coverage. Similarly, in 2018, SMA had limited parallactic angle coverage for M87*, with only ∼5°  evolution of the parallactic angle. As a result, the leakage estimate for SMA is highly uncertain from just the M87* coverage. In 2021, we found that all leakage estimates, besides GLT, are consistent to within ≲3%. GLT still has the largest spread between methods due to its poor parallactic angle coverage, although the spread improved by a factor of ∼2 compared to 2018

Fig. E.1. Leakage estimates across all methods for 2018 and 2021. In 2018, we flagged JCMT (Appendix B.2). For the Bayesian methods, the error bars show the 95% credible intervals about the median leakage denoted by the marker.

Finally, the uncertainty of GLT’s leakage in 2018 and 2021, and SMA’s in 2018, does not appear to dramatically alter the image reconstructions across pipelines (see Fig. F.1) or the polarized parameter estimates (see Fig. 7). Although we expect leakage estimates for GLT and SMA to improve from multisource fitting in the future, it is unlikely that they will significantly alter the image reconstructions or change our quantitative results.

Appendix F: Individual method reconstructions of M87*

Figure F.1 shows the results of each imaging method for the three epochs presented, blurred to a common resolution. In addition to the total flux differences, we found excellent agreement in the structure of M87* every year. For 2017, all methods found that the south-west region of the ring is the most polarized region and has a similar total intensity azimuthal profile. However, for the same year, there are some differences in the EVPA pattern, especially in the eastern and northern parts of the images. This uncertainty was also found in EHTC (2021a), and is a result of the sparseness of the 2017 coverage.

Fig. F.1. Reconstructions of M87* across 2017, 2018, and 2021 from all methods considered. All methods have been blurred to match the resolution of the two CLEAN reconstructions GPCAL and LPCAL.

Remarkably, the uncertainty in the linear polarization and total intensity images in 2018, especially in 2021, has been significantly reduced. All methods found that the position angle of the ring brightness has shifted in 2018 and 2021 compared to 2017. Furthermore, all methods found that M87* in 2018 is significantly de-polarized, where only the western part of the image has significant polarized emission across all methods. For the 2021 reconstructions, all methods found a very similar EVPA pattern, further demonstrating the robustness of the ∠β₂ rotation and the changes in the EVPA helicity of the ring. Furthermore, unlike 2017, the brightest part of the total intensity was de-polarized relative to other parts in 2021.

Finally, Fig. F.2 shows the uncertainty maps for the Comrade and THEMIS pipelines. We found that the central ring has a small fractional uncertainty over all three years, demonstrating its robustness. Furthermore, the regions of the ring that are polarized in 2017 (SW), 2018 (W), and 2021 (E and W) have a small fractional uncertainty in total linear polarization and overall uncertainty in EVPA.

Fig. F.2. Uncertainty maps for the 2017 (top), 2018 (middle), and 2021 (bottom) M87* reconstructions using the Comrade (left) and THEMIS (right) pipelines. Left columns: Fractional total intensity uncertainty. Middle columns: Fractional linear polarization uncertainty. Right columns: Absolute EVPA uncertainty, which is defined as the circular standard deviation in radians of each pixel. Any pixels that are outside the field of view of the raster are assigned an uncertainty of > 1.

In summary, all key results presented in Sect. 6 are seen in individual image reconstructions and do not depend on the particular choices of the imaging algorithms.

Appendix G: Acknowledgments

The Event Horizon Telescope Collaboration thanks the following organizations and programs: the Academia Sinica; the Academy of Finland (projects 274477, 284495, 312496, 315721); the Agencia Nacional de Investigación y Desarrollo (ANID), Chile via NCN19_058 (TITANs), Fondecyt 1221421 and BASAL FB210003; the Alexander von Humboldt Stiftung (including the Feodor Lynen Fellowship); an Alfred P. Sloan Research Fellowship; Allegro, the European ALMA Regional Centre node in the Netherlands, the NL astronomy research network NOVA and the astronomy institutes of the University of Amsterdam, Leiden University, and Radboud University; the ALMA North America Development Fund; the Astrophysics and High Energy Physics programme by MCIN (with funding from European Union NextGenerationEU, PRTR-C17I1); the Black Hole Initiative, which is funded by grants from the John Templeton Foundation (60477, 61497, 62286) and the Gordon and Betty Moore Foundation (Grant GBMF-8273) - although the opinions expressed in this work are those of the author and do not necessarily reflect the views of these Foundations; the Brinson Foundation; the Canada Research Chairs (CRC) program; Chandra DD7-18089X and TM6-17006X; the China Scholarship Council; the China Postdoctoral Science Foundation fellowships (2020M671266, 2022M712084); ANID through Fondecyt Postdoctorado (project 3250762); Conicyt through Fondecyt Postdoctorado (project 3220195); Consejo Nacional de Humanidades, Ciencia y Tecnología (CONAHCYT, Mexico, projects U0004-246083, U0004-259839, F0003-272050, M0037-279006, F0003-281692, 104497, 275201, 263356, CBF2023-2024-1102, 257435); the Colfuturo Scholarship; the Consejo Superior de Investigaciones Científicas (grant 2019AEP112); the Delaney Family via the Delaney Family John A. Wheeler Chair at Perimeter Institute; Dirección General de Asuntos del Personal Académico-Universidad Nacional Autónoma de México (DGAPA-UNAM, projects IN112820 and IN108324); the Dutch Research Council (NWO) for the VICI award (grant 639.043.513), the grant OCENW.KLEIN.113, and the Dutch Black Hole Consortium (with project No. NWA 1292.19.202) of the research programme the National Science Agenda; the Dutch National Supercomputers, Cartesius and Snellius (NWO grant 2021.013); the EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, Center for Astronomical Mega-Science, Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute; the European Research Council (ERC) Synergy Grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (grant 610058) and Synergy Grant “BlackHolistic: Colour Movies of Black Holes: Understanding Black Hole Astrophysics from the Event Horizon to Galactic Scales” (grant 10107164); the European Union Horizon 2020 research and innovation programme under grant agreements RadioNet (No. 730562), M2FINDERS (No. 101018682) and FunFiCO (No. 777740); the European Research Council for advanced grant “JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales” (grant No. 884631); the European Horizon Europe staff exchange (SE) programme HORIZON-MSCA-2021-SE-01 grant NewFunFiCO (No. 10108625); the Horizon ERC Grants 2021 programme under grant agreement No. 101040021; the FAPESP (Fundação de Amparo á Pesquisa do Estado de São Paulo) under grant 2021/01183-8; the Fondes de Recherche Nature et Technologies (FRQNT); the Fondo CAS-ANID folio CAS220010; the Generalitat Valenciana (grants APOSTD/2018/177 and ASFAE/2022/018) and GenT Program (project CIDEGENT/2018/021); the Gordon and Betty Moore Foundation (GBMF-3561, GBMF-5278, GBMF-10423); the Institute for Advanced Study; the ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU; the Istituto Nazionale di Fisica Nucleare (INFN) sezione di Napoli, iniziative specifiche TEONGRAV; the International Max Planck Research School for Astronomy and Astrophysics at the Universities of Bonn and Cologne; the Italian Ministry of University and Research (MUR)– Project CUP F53D23001260001, funded by the European Union – NextGenerationEU; Deutsche Forschungsgemeinschaft (DFG) research grant “Jet physics on horizon scales and beyond” (grant No. 443220636) and DFG research grant 443220636; Joint Columbia/Flatiron Postdoctoral Fellowship (research at the Flatiron Institute is supported by the Simons Foundation); the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT; grant JPMXP1020200109); the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS Research Fellowship (JP17J08829); the Joint Institute for Computational Fundamental Science, Japan; the Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS, grants QYZDJ-SSW-SLH057, QYZDJSSW-SYS008, ZDBS-LY-SLH011); the Leverhulme Trust Early Career Research Fellowship; the Max-Planck-Gesellschaft (MPG); the Max Planck Partner Group of the MPG and the CAS; the MEXT/JSPS KAKENHI (grants 18KK0090, JP21H01137, JP18H03721, JP18K13594, 18K03709, JP19K14761, 18H01245, 25120007, 19H01943, 21H01137, 21H04488, 22H00157, 23K03453); the MICINN Research Projects PID2019-108995GB-C22, PID2022-140888NB-C22; the MIT International Science and Technology Initiatives (MISTI) Funds; the Ministry of Science and Technology (MOST) of Taiwan (103-2119-M-001-010-MY2, 105-2112-M-001-025-MY3, 105-2119-M-001-042, 106-2112-M-001-011, 106-2119-M-001-013, 106-2119-M-001-027, 106-2923-M-001-005, 107-2119-M-001-017, 107-2119-M-001-020, 107-2119-M-001-041, 107-2119-M-110-005, 107-2923-M-001-009, 108-2112-M-001-048, 108-2112-M-001-051, 108-2923-M-001-002, 109-2112-M-001-025, 109-2124-M-001-005, 109-2923-M-001-001, 110-2112-M-001-033, 110-2124-M-001-007 and 110-2923-M-001-001); the National Science and Technology Council (NSTC) of Taiwan (111-2124-M-001-005, 112-2124-M-001-014, 112-2112-M-003-010-MY3, and 113-2124-M-001-008); the Ministry of Education (MoE) of Taiwan Yushan Young Scholar Program; the Physics Division, National Center for Theoretical Sciences of Taiwan; the National Aeronautics and Space Administration (NASA, Fermi Guest Investigator grant 80NSSC23K1508, NASA Astrophysics Theory Program grant 80NSSC20K0527, NASA NuSTAR award 80NSSC20K0645); NASA Hubble Fellowship Program Einstein Fellowship; NASA Hubble Fellowship grants HST-HF2-51431.001-A, HST-HF2-51482.001-A, HST-HF2-51539.001-A, HST-HF2-51552.001A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555; the National Institute of Natural Sciences (NINS) of Japan; the National Key Research and Development Program of China (grant 2016YFA0400704, 2017YFA0402703, 2016YFA0400702); the National Science and Technology Council (NSTC, grants NSTC 111-2112-M-001 -041, NSTC 111-2124-M-001-005, NSTC 112-2124-M-001-014); the US National Science Foundation (NSF, grants AST-0096454, AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433, OIA-1126433, AST-1140030, DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509, OPP-1248097, AST-1310896, AST-1440254, AST-1555365, AST-1614868, AST-1615796, AST-1715061, AST-1716327, AST-1726637, OISE-1743747, AST-1743747, AST-1816420, AST-1935980, AST-1952099, AST-2034306, AST-2205908, AST-2307887, AST-2407810); NSF Astronomy and Astrophysics Postdoctoral Fellowship (AST-1903847); the Natural Science Foundation of China (grants 11650110427, 10625314, 11721303, 11725312, 11873028, 11933007, 11991052, 11991053, 12192220, 12192223, 12273022, 12325302, 12303021); the Natural Sciences and Engineering Research Council of Canada (NSERC); the National Research Foundation of Korea (the Global PhD Fellowship Grant: grants NRF-2015H1A2A1033752; the Korea Research Fellowship Program: NRF-2015H1D3A1066561; Brain Pool Program: RS-2024-00407499; Basic Research Support Grant 2019R1F1A1059721, 2021R1A6A3A01086420, 2022R1C1C1005255, 2022R1F1A1075115); the POSCO Science Fellowship of the POSCO TJ Park Foundation; Netherlands Research School for Astronomy (NOVA) Virtual Institute of Accretion (VIA) postdoctoral fellowships; NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation; Onsala Space Observatory (OSO) national infrastructure, for the provisioning of its facilities/observational support (OSO receives funding through the Swedish Research Council under grant 2017-00648); the Perimeter Institute for Theoretical Physics (research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science); the Portuguese Foundation for Science and Technology (FCT) grants (Individual CEEC program - 5th edition, https://doi.org/10.54499/UIDB/04106/2020, https://doi.org/10.54499/UIDP/04106/2020, PTDC/FIS-AST/3041/2020, CERN/FIS-PAR/0024/2021, 2022.04560.PTDC); the Princeton Gravity Initiative; the Spanish Ministerio de Ciencia, Innovación y Universidades (grants PID2022-140888NB-C21, PID2022-140888NB-C22, PID2023-147883NB-C21, RYC2023-042988-I); the Severo Ochoa grant CEX2021-001131-S funded by MICIU/AEI/10.13039/501100011033; The European Union’s Horizon Europe research and innovation program under grant agreement No. 101093934 (RADIOBLOCKS); The European Union “NextGenerationEU", the Recovery, Transformation and Resilience Plan, the CUII of the Andalusian Regional Government and the Spanish CSIC through grant AST22_00001_Subproject_10; “la Caixa” Foundation (ID 100010434) through fellowship codes LCF/BQ/DI22/11940027 and LCF/BQ/DI22/11940030; the University of Pretoria for financial aid in the provision of the new Cluster Server nodes and SuperMicro (USA) for a SEEDING GRANT approved toward these nodes in 2020; the Shanghai Municipality orientation program of basic research for international scientists (grant no. 22JC1410600); the Shanghai Pilot Program for Basic Research, Chinese Academy of Science, Shanghai Branch (JCYJ-SHFY-2021-013); the Simons Foundation (grant 00001470); the Spanish Ministry for Science and Innovation grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033; the Spinoza Prize SPI 78-409; the South African Research Chairs Initiative, through the South African Radio Astronomy Observatory (SARAO, grant ID 77948), which is a facility of the National Research Foundation (NRF), an agency of the Department of Science and Innovation (DSI) of South Africa; the Swedish Research Council (VR); the Taplin Fellowship; the Toray Science Foundation; the UK Science and Technology Facilities Council (grant no. ST/X508329/1); the US Department of Energy (USDOE) through the Los Alamos National Laboratory (operated by Triad National Security, LLC, for the National Nuclear Security Administration of the USDOE, contract 89233218CNA000001); and the YCAA Prize Postdoctoral Fellowship. This work was also supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government(MSIT) (RS-2024-00449206). We acknowledge support from the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) of Brazil through PROEX grant number 88887.845378/2023-00. We acknowledge financial support from Millenium Nucleus NCN23_002 (TITANs) and Comité Mixto ESO-Chile. We thank the staff at the participating observatories, correlation centers, and institutions for their enthusiastic support. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00841.V and ADS/JAO.ALMA#2019.1.01797.V. ALMA is a partnership of the European Southern Observatory (ESO; Europe, representing its member states), NSF, and National Institutes of Natural Sciences of Japan, together with National Research Council (Canada), Ministry of Science and Technology (MOST; Taiwan), Academia Sinica Institute of Astronomy and Astrophysics (ASIAA; Taiwan), and Korea Astronomy and Space Science Institute (KASI; Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc. (AUI)/NRAO, and the National Astronomical Observatory of Japan (NAOJ). The NRAO is a facility of the NSF operated under cooperative agreement by AUI. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC05-00OR22725; the ASTROVIVES FEDER infrastructure, with project code IDIFEDER-2021-086; the computing cluster of Shanghai VLBI correlator supported by the Special Fund for Astronomy from the Ministry of Finance in China; We also thank the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported by FAPESP (Fundacao de Amparo a Pesquisa do Estado de Sao Paulo) under grant 2021/01183-8. APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (Germany), ESO, and the Onsala Space Observatory (Sweden). The SMA is a joint project between the SAO and ASIAA and is funded by the Smithsonian Institution and the Academia Sinica. The JCMT is operated by the East Asian Observatory on behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of China, Chinese Academy of Sciences, and the National Key Research and Development Program (No. 2017YFA0402700) of China and Natural Science Foundation of China grant 11873028. Additional funding support for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada. The LMT is a project operated by the Instituto Nacional de Astrófisica, Óptica, y Electrónica (Mexico) and the University of Massachusetts at Amherst (USA). The IRAM 30 m telescope on Pico Veleta, Spain and the NOEMA interferometer on Plateau de Bure, France are operated by IRAM and supported by CNRS (Centre National de la Recherche Scientifique, France), MPG (Max-Planck-Gesellschaft, Germany), and IGN (Instituto Geográfico Nacional, Spain). The SMT is operated by the Arizona Radio Observatory, a part of the Steward Observatory of the University of Arizona, with financial support of operations from the State of Arizona and financial support for instrumentation development from the NSF. Support for SPT participation in the EHT is provided by the National Science Foundation through award OPP-1852617 to the University of Chicago. Partial support is also provided by the Kavli Institute of Cosmological Physics at the University of Chicago. The SPT hydrogen maser was provided on loan from the GLT, courtesy of ASIAA. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442. XSEDE Stampede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N. XSEDE JetStream resource at PTI and TACC was allocated through AST170028. This research is part of the Frontera computing project at the Texas Advanced Computing Center through the Frontera Large-Scale Community Partnerships allocation AST20023. Frontera is made possible by National Science Foundation award OAC-1818253. This research was done using services provided by the OSG Consortium (Pordes et al. 2007; Sfiligoi et al. 2009), which is supported by the National Science Foundation award Nos. 2030508 and 1836650. Additional work used ABACUS2.0, which is part of the eScience center at Southern Denmark University, and the Kultrun Astronomy Hybrid Cluster (projects Conicyt Programa de Astronomia Fondo Quimal QUIMAL170001, Conicyt PIA ACT172033, Fondecyt Iniciacion 11170268, Quimal 220002). Simulations were also performed on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, on the HazelHen cluster at the HLRS in Stuttgart, and on the Pi2.0 and Siyuan Mark-I at Shanghai Jiao Tong University. The computer resources of the Finnish IT Center for Science (CSC) and the Finnish Computing Competence Infrastructure (FCCI) project are acknowledged. This research was enabled in part by support provided by Compute Ontario (http://computeontario.ca), Calcul Quebec (http://www.calculquebec.ca), and the Digital Research Alliance of Canada (https://alliancecan.ca/en). The EHTC has received generous donations of FPGA chips from Xilinx Inc., under the Xilinx University Program. The EHTC has benefited from technology shared under open-source license by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER). The EHT project is grateful to T4Science and Microsemi for their assistance with hydrogen masers. This research has made use of NASA’s Astrophysics Data System. We gratefully acknowledge the support provided by the extended staff of the ALMA, from the inception of the ALMA Phasing Project through the observational campaigns of 2017 and 2018. We would like to thank A. Deller and W. Brisken for EHT-specific support with the use of DiFX. We thank Martin Shepherd for the addition of extra features in the Difmap software that were used for the CLEAN imaging results presented in this paper. We acknowledge the significance that Maunakea, where the SMA and JCMT EHT stations are located, has for the indigenous Hawaiian people.

All Tables

All Figures

	Fig. 1. EHT in its 2021 configuration. Compared to the original 2017 array, GLT was added in 2018, and KP and NOEMA joined the EHT for the 2021 campaign (indicated in blue). Baselines from SPT and LMT are greyed out since SPT cannot observe M87* (only its calibrator 3C 279), and LMT did not observe in 2021.
In the text

Fig. 2. (u, v) coverage of M87* during the 2017 (left), 2018 (middle), and 2021 (right) campaigns for the band 3 (227.1 GHz) observations. The ticks show each year’s interferometric EVPA and the colour of the observed interferometric fractional linear polarization, after de-biasing for thermal noise and applying leakage corrections. The leakage terms used are the fiducial values from EHTC (2021a) in 2017 and the THEMIS leakage solutions in 2018 and 2021. The clusters of high polarization fractions in 2018 come from the GLT, which was underperforming in 2018 due to an incomplete commissioning at that time, as described in Koay et al. (2023a).

In the text

Fig. 3. Band 3 (227.1 GHz) M87* total intensity amplitude and phase data measured in 2017, 2018, and 2021. The calibration applied is the pipeline-based signal stabilization and a priori flux density calibration without image-based self-calibration. The 2017 and 2018 data were produced by EHT-HOPS, and the 2021 data by rPICARD. The grey bands around 3.4 Gλ indicate the approximate location of the first visibility minima.

In the text

Fig. 4. Comparisons of the conjugate closure trace phases across the three observing campaigns for two quadrangles. Deviations from arg(𝒞) = 0 indicate the presence of significant polarization structure, independent of leakage and gain calibration. Closed and open points show low- and high-band data, respectively, offset in time for clarity. The model predictions from the low-band reconstructions described in Sect. 3 are shown for reference (for the Bayesian methods, THEMIS and Comrade, the 95-percentile range is shown along with five draws from the posterior). Note the variation in the vertical axis ranges between panels. Because LMT did not observe during the 2021 campaign, we show the ALMA-PV-NOEMA-GLT quadrangle in 2021, for which clear deviations from zero are apparent.

In the text

Fig. 5. Comparison of extracted parameters for the blinded synthetic data test. All methods are blurred to match the ground truth GRMHD, blurred to 20 μas before the parameters are estimated. The markers show the results for each method. For the two Bayesian methods, Comrade and THEMIS, the markers denote the median, and the error bars the 95% credible interval. Since the synthetic data only consider a single frequency, the non-Bayesian methods only produce a single image. The dashed black line shows the ground truth, estimated from the true on-sky image blurred with a 20 μas Gaussian beam.

In the text

Fig. 6. Fiducial images for 2017 (band 3), 2018 (bands 3 and 4), and 2021 (bands 3 and 4). The images were produced by averaging the reconstructions over the methods described in Sect. 3. Each method’s image reconstruction has been blurred to match the resolution of the representative LPCAL image. Top row: Polarization ‘field lines’ overlaid on the total intensity image. Second row: Total intensity image in grey scale with the contours showing the 22.5%, 45%, 67.5%, and 90% peak brightness levels, overlaid with polarization ticks. The polarization ticks indicate the EVPA, the tick length is proportional to the linear polarization intensity, and their colour indicates the linear polarization fraction. Polarization ticks are only shown in regions where the total intensity is > 10% of the maximum brightness and the linear polarization brightness is > 10% of the peak linear polarized brightness. Bottom row: Total linear polarized brightness, |𝒫|.

In the text

Fig. 7. Extracted parameters of M87* across the three observations and averaged over band 3 and band 4. Each panel shows the results from April 11, 2017, April 21, 2018, and April 18, 2021. For the non-Bayesian methods, the error bars show the spread between the high- and low-band estimates and are not a measure of the statistical uncertainty of the image reconstruction. The error bars for the Bayesian methods show the 95% credible interval around the median and measure the statistical uncertainty of the reconstructions due to thermal noise, instrumental effects, incomplete coverage, and frequency dependence. The solid black line is the median, and the grey band is the 95-percentile range from all image reconstructions, with each method weighted inversely to the number of images they produced. Note that each image reconstruction has been blurred to match the resolution of the LPCAL reconstruction before the parameters were estimated.

In the text

Fig. 8. Closure phases (band 4) on the three closure triangles that probe extended non-ring emission trivial closures – to AA-AX-PV (top), ALMA-PV-NOEMA (middle), and ALMA-SMT-KP (bottom) – with fits from the final image reconstructions for the types of imaging algorithms – Bayesian (Comrade), CLEAN (GPCAL), and the RML method (ehtim). Other methods yield qualitatively similar fits where available. Strategies to fit the extended non-ring emission to the 2021 closure phases are presented in Georgiev et al. (2025) and Saurabh et al. (2025).

In the text

	Fig. A.1. M87* (u, v) coverage of VLBI scans shown for the band 3 (227.1 GHz) data on April 13 and 18, 2021. The co-located ALMA, APEX, JCMT, and SMA stations are described by the Chile and Hawaii points, respectively. The labelling is done by antenna, using conjugate baselines to display the other station.
In the text

	Fig. A.2. RR* − LL* closure phase differences that describe source circular polarization structure shown for two triangles in 2017, 2018, and 2021. The data are averaged over 12-minute-long segments.
In the text

	Fig. A.3. Image reconstructions averaged over bands 3 and 4 using the Comrade pipeline and the April 18, 2021, reduction from the standard rPICARD pipeline used in this paper (left) and the EHT-HOPS pipeline (right). Both images have been normalized to have the same compact flux due to residual issues with the EHT-HOPS calibration of NOEMA.
In the text

	Fig. A.4. Fiducial images from April 13, 2021, and April 18, 2021. The April 13 image is constructed in the same manner as Fig. 6.
In the text

	Fig. A.5. Comparison of imaging parameter estimation results from April 13 and 18, 2021. The conventions follow Fig. 7.
In the text

	Fig. D.1. Polarized imaging results from the blinded GRMHD synthetic data test. The leftmost column shows the ground truth GRMHD image blurred by a 20 μas Gaussian beam; the other columns show the imaging results from each method blurred to an equivalent resolution. From top to bottom, the rows correspond to the 2017, 2018, and 2021 results.
In the text

Fig. D.2. D-term recovery results of blinded GRMHD synthetic data from each polarized imaging method (columns) and year (rows). The x-axis of each plot is the true leakage, and the y-axis is the measured leakage. Perfect recovery means that all points should live on the y = x line. The total L2 difference between the ground truth and recovered leakage averaged over all stations is shown in the bottom right of each plot. For the Bayesian methods, we also show a χ2 value computed as the square difference between the true value and the recovered leakage divided by the variance of the posterior estimate. For the Bayesian methods, the points denote the posterior median leakage, and error bars show the 95% credible interval about the median.

In the text

Fig. D.3. Cross-correlation between the low polarized and high polarized extended jet geometric tests. The ring model for the two tests is identical, meaning that the cross-correlation in the core region would be equal to unity for perfect reconstructions. Left column: Total intensity (top row) and linear polarization cross-correlation (bottom row) between the low-polarized jet reconstructions. Right column: Same quantities but with the cross-correlation computed between the low and high polarized jet reconstructions.

In the text

	Fig. E.1. Leakage estimates across all methods for 2018 and 2021. In 2018, we flagged JCMT (Appendix B.2). For the Bayesian methods, the error bars show the 95% credible intervals about the median leakage denoted by the marker.
In the text

	Fig. F.1. Reconstructions of M87* across 2017, 2018, and 2021 from all methods considered. All methods have been blurred to match the resolution of the two CLEAN reconstructions GPCAL and LPCAL.
In the text

Fig. F.2. Uncertainty maps for the 2017 (top), 2018 (middle), and 2021 (bottom) M87* reconstructions using the Comrade (left) and THEMIS (right) pipelines. Left columns: Fractional total intensity uncertainty. Middle columns: Fractional linear polarization uncertainty. Right columns: Absolute EVPA uncertainty, which is defined as the circular standard deviation in radians of each pixel. Any pixels that are outside the field of view of the raster are assigned an uncertainty of > 1.

In the text