First M87 Event Horizon Telescope Results. VII. Polarization of the Ring

The giant elliptical galaxy Messier 87 (M87, NGC 4486) is the central member of the Virgo cluster of galaxies and hosts a low-luminosity radio source (Virgo A, 3C 274, B1228+126). M87 is nearby and bright, and at its center is one of the best-studied active galactic nuclei (AGNs). M87 was the first galaxy in which an extragalactic jet (first described as a "narrow ray") extending from the nucleus was discovered (Curtis 1918). This kiloparsec-scale jet is visible, with remarkably similar morphology, at all wavelengths from radio to X-ray. The optical radiation from the jet on kpc scales was found to be linearly polarized by Baade (1956), which was confirmed by Hiltner (1959), suggesting that the emission mechanism is synchrotron radiation.

The central engine that powers the jet contains one of the most massive black holes known, measured from the central stellar velocity dispersion (Gebhardt et al. 2011; M =(6.6 ± 0.4) × 10⁹ M_⊙) and directly from the size of the observed emitting region surrounding the black hole shadow (Paper VI; M = (6.5 ± 0.7) × 10⁹ M_⊙). For this mass, the Schwarzschild radius is R_s = 2GM/c² = 1.8 × 10¹⁵ cm. At the distance of M87, ${16.8}_{-0.7}^{+0.8}\,$ Mpc (Blakeslee et al. 2009; Bird et al. 2010; Cantiello et al. 2018, Paper VI), the EHT resolution of about 20 micro-arcseconds (μas) translates into a linear scale of 0.0016 pc = 2.5 R_s.

The M87 jet has been imaged at subarcsecond resolution in both total intensity and linear polarization at optical wavelengths with the Hubble Space Telescope (Thomson et al. 1995; Capetti et al. 1997), and at radio wavelengths with the Very Large Array (e.g., Owen et al. 1989). Observing the launching region of the jet closer to the black hole and the region surrounding the black hole requires milliarcsecond (mas) resolution or better, and hence very-long-baseline interferometry (VLBI) techniques used at the highest frequencies (e.g., Boccardi et al. 2017 and references therein).

Milliarcsecond-scale VLBI observations show that the core itself is unpolarized even at millimeter wavelengths. Zavala & Taylor (2002), observing at 8, 12, and 15 GHz, set upper limits on the fractional polarization of the compact core of m < 0.1%. About 20 mas downstream from the core, patchy linear polarization starts to become visible in the jet at the level of 5%–10%, although no large-scale coherent pattern to the electric-vector position angles (EVPAs) χ is apparent. However, at each patch in the downstream jet, the EVPAs exhibit a linear change with λ², allowing the rotation measures (RMs) to be estimated. These RMs range from −4000 rad m⁻² to 9000 rad m⁻² (Zavala & Taylor 2002). The linear dependence of EVPA on λ² over several radians is important, as it shows that the Faraday-rotating plasma in the jet cannot be mixed in with the relativistic emitting particles (Burn 1966) but must be in a cooler (sub-relativistic) foreground screen.

On kiloparsec scales, Owen et al. (1990) found a complex distribution of RM. Over most of the source the RM is typically of order 1000 rad m⁻², but there are patches where values as high as 8000 rad m⁻² are found.

More recently, Park et al. (2019) studied Faraday RMs in the jet using multifrequency Very Long Baseline Array (VLBA) data at ≲8 GHz. They found that the RM magnitude systematically decreases with increasing distance from 5,000 to 200,000 R_s. The observed large (≳45°) EVPA rotations at various locations of the jet suggest that the dominant Faraday screen in this distance range would be external to the jet, similar to the conclusion of Zavala & Taylor (2002). Homan & Lister (2006), also observing at 15 GHz with the VLBA (as part of the MOJAVE project) found a tight upper limit on the fractional linear polarization of the core of <0.07%. They also detect circular polarization of (−0.49 ± 0.10)%.

At 43 GHz, Walker et al. (2018) presented results from 17 years of VLBA observations of M87, with polarimetric images presented at two epochs. These show significant polarization (up to 4%) in the jet near the 43 GHz core, but at the position of the total intensity peaks the fractional polarizations are only 1.5% and 1.1%. They interpret these fractions as coming from a mix of emission from the unresolved, unpolarized core and a more polarized inner jet.

Hada et al. (2016) showed images at four epochs at 86 GHz made with the VLBA and the Green Bank Telescope. At this frequency, the resolution is about (0.4 × 0.1) mas, corresponding to (56 × 14) R_s. Again, the core is unpolarized with no linear polarization detected at the position of the total intensity emission's peak, while there is a small patch of significant (3.5%) polarization located 0.1 mas downstream. At 0.4 mas downstream from the peak, there is another patch of significant polarization (20%). These results indicate that there are regions of significantly ordered magnetic field very close to the central engine.

Very recently, new observations by Kravchenko et al. (2020) using the VLBA at 22 and 43 GHz show two components of linear polarization and a smooth rotation of EVPA around the 43 GHz core. Comparison with earlier observations show that the global polarization pattern in the jet is largely stable over an 11 year timescale. They suggest that the polarization pattern is associated with the magnetic structure in a confining magnetohydrodynamic wind, which is also the source of the observed Faraday rotation.

The EHT presently observes at ∼230 GHz and has previously reported polarimetric measurements only for Sagittarius A* (Sgr A*; Johnson et al. 2015). The only previous polarimetric measurements of M87 at this frequency were done by Kuo et al. (2014) using the Submillimeter Array (SMA) on Maunakea, Hawai'i, USA. The SMA is a compact array with a (1.2 × 0.8) arcsec beam, 10000 times larger than the EHT beam. Li et al. (2016) used the value from this work to calculate a limit on the accretion rate onto the M87 black hole. Most recently, Goddi et al. (2021) reported results on M87 around 230 GHz as part of the Atacama Large Millimeter/submillimeter Array (ALMA) interferometric connected-element array portion of the EHT observations in 2017. The ALMA-only 230 GHz observations (with a FWHM synthesized beam in the range 1''–2'', depending on the day) resolve the M87 inner region into a compact central core and a kpc-scale jet across approximately 25''. It has been found that the 230 GHz core at these scales has a total flux density of ∼1.3 Jy, a low linear polarization fraction ∣m∣ ∼ 2.7%, and even less circular polarization, ∣v∣ < 0.3%. Notably, ALMA-only observations show strong variability in the RM estimated based on four frequencies within ALMA Band 6 (four spectral windows centered at 213, 215, 227, and 229 GHz; Matthews et al. 2018). The RM difference is clear between the start of the EHT observing campaign on 2017 April 5 (RM ≈ 0.6 × 10⁵ rad m⁻²) and the end on 2017 April 11 (RM ≈ − 0.4 × 10⁵ rad m⁻²). Because these measurements were taken simultaneously with the EHT VLBI observations presented here, the ALMA-only linear polarization fraction measurements can be used as a point of reference, and we discuss possible implications of the strong RM evolution on the EHT polarimetric images of M87.

This Letter presents the details of the polarimetric data calibration, the procedures for polarimetric imaging, and the resulting images of the M87 core. In Section 2, we briefly overview the basics of polarimetric VLBI. In Section 3, we summarize the EHT 2017 observations, describe the initial data calibration procedure and validation tests, and describe the basic properties of the polarimetric data. In Section 4, we describe our methods, strategy, and test suite for our polarimetric calibration and imaging. In Section 5, we present and analyze the polarimetric images of the M87 ring and examine the calibration's impact on the polarimetric image. We discuss the results and summarize the work in Sections 6 and 7.

This Letter is supplemented with a number of appendices supporting our analysis and results. The appendices summarize: polarimetric data issues (Appendix A); novel VLBI closure data products (Appendix B); details of calibration and imaging methods (Appendix C); validation of polarimetric calibration for telescopes with an intra-site partner (Appendix D); fiducial leakage D-terms from M87 imaging (Appendix E); preliminary results of polarimetric imaging of M87 (Appendix F); polarimetric imaging scoring procedures (Appendix G); details of Monte Carlo D-term simulations (Appendix H); consistency of low- and high-band results for M87 (Appendix I); comparison to polarimetric properties of calibrator sources (Appendix J); and validations of assumptions made in polarimetric imaging of the main target and the calibrators (Appendix K).

Eight observatories at six geographical locations participated in the 2017 EHT observing campaign: ALMA and the Atacama Pathfinder Experiment (APEX) in the Atacama Desert in Chile; the Large Millimeter Telescope Alfonso Serrano (LMT) on the Volcán Sierra Negra in Mexico; the South Pole Telescope (SPT) at the geographic south pole; the IRAM 30 m telescope (PV) on Pico Veleta in Spain; the Submillimeter Telescope (SMT) on Mt. Graham in Arizona, USA; SMA and the James Clerk Maxwell Telescope (JCMT) on Maunakea in Hawai'i, USA. ¹³¹

The EHT observations were carried out on five nights between 2017 April 5 and 11. M87 was observed on April 5, 6, 10, and 11. Along with the main EHT targets M87 and Sgr A*, several other AGN sources were observed as science targets and calibrators.

Observations were conducted using two contiguous frequency bands of 2 GHz bandwidth each, centered at frequencies of 227.1 and 229.1 GHz, hereby referred to as low and high band, respectively. The observations were arranged in scans alternating different sources, with durations lasting between 3 and 7 minutes. Apart from the JCMT, which observed only a single polarization (right-circular polarization on 2017 April 5–7 and left-circular polarization on 2017 April 10–11), all stations observed in full polarization mode. ALMA is the only station to natively record data in a linear polarization basis. Visibilities measured on baselines to ALMA were converted from a mixed linear-circular basis to circular polarization after correlation using the PolConvert software (Martí-Vidal et al. 2016; Matthews et al. 2018; Goddi et al. 2019). A technical description of the EHT array is presented in Paper II and a summary of the 2017 observations and data reduction is presented in Paper III.

After the sky signal received at each telescope was mixed to baseband, digitized, and recorded directly to hard disk, the data from each station were sent to MIT Haystack Observatory and the Max-Planck-Institut für Radioastronomie (MPIfR) for correlation using the DiFX software correlators (Deller et al. 2011). The accumulation period adopted at correlation is 0.4 s, with a frequency resolution of 0.5 MHz. The clock model used during correlation to align the wavefronts arriving at different telescopes is imperfect, owing to an approximate a priori model for Earth's geometry as well as rapid stochastic variations in path length due to local atmospheric turbulence (Paper III). Before the data can be averaged coherently to build up signal-to-noise ratio (S/N), these effects must be accurately measured and corrected. This process, referred to as fringe fitting, was conducted using three independent software packages: the Haystack Observatory Processing System (HOPS; Whitney et al. 2004; Blackburn et al. 2019); the Common Astronomy Software Applications package (CASA; McMullin et al. 2007; Janssen et al. 2019a); and the NRAO Astronomical Image Processing System (AIPS; Greisen 2003, Paper III). Automated reduction pipelines were designed specifically to address the unique challenges related to the heterogeneity, wide bandwidth, and high observing frequency of EHT data. The field rotation angle is corrected with Equations (4)–(5), using coefficients given in Table 1. Flux density (amplitude) calibration is applied via a common post-processing framework for all pipelines (Blackburn et al. 2019; Paper III), taking into account estimated station sensitivities (Issaoun et al. 2017; Janssen et al. 2019b). Under the assumption of zero circular polarization of the primary (solar system) calibrator sources, elevation-independent station gains possess independent statistical uncertainties for the right-hand-circular polarization (RCP) and left-hand-circular polarization (LCP) signal paths, estimated to be ∼20% for the LMT and ∼10% for all other stations (Janssen et al. 2019b).

To remove the instrumental amplitude mismatch between the LL^* and RR^* visibility components (the R–L phases are correctly calibrated in all scans by using ALMA as the reference station), calibration of the complex polarimetric gain ratios (the ratios of the G_R and G_L terms in the G matrices) is performed. This is done by fitting global (multi-source, multi-days) piecewise polynomial gain ratios as functions of time. The aim of this approach is to preserve differences in LL^* and RR^* visibilities intrinsic to the source (Steel et al. 2019). After this step, preliminary polarimetric Stokes visibilities $\tilde{{ \mathcal I }},\tilde{{ \mathcal Q }},\tilde{{ \mathcal U }},\tilde{{ \mathcal V }}$ can be constructed. However, the gain calibration requires significant additional improvements. The final calibration of the station phase and amplitude gains takes place in a self-calibration step as part of imaging or modeling the Stokes ${ \mathcal I }$ brightness distribution, preserving the complex polarimetric gain ratios (e.g., Papers IV, VI). Fully calibrating the D-terms requires modeling the polarized emission.

The Stokes ${ \mathcal I }$ (total intensity) analysis of a subset of the 2017 observations (Science Release 1 (SR1)), including M87, was the subject of Papers I–VI. The quality of these Stokes ${ \mathcal I }$ data was assured by a series of tests covering self-consistency over bands and parallel hand polarizations, and consistency of trivial closure quantities (Wielgus et al. 2019). Constraints on the residual non-closing errors were found to be at a 2% level.

For additional information on the calibration, data reduction, and validation procedures for EHT, see Paper III. Information about accessing SR1 data and the software used for analysis can be found on the EHT website's data portal. ¹³² In this Letter, we utilize the HOPS pipeline full-polarization band-averaged (i.e., averaged over frequency within each band) and 10-second averaged data set from the same reduction path as SR1, but containing a larger sample of calibrator sources for polarimetric leakage studies. In addition, the ALMA linear-polarization observing mode allows us to measure and recover the absolute EVPA in the calibrated VLBI visibilities (Martí-Vidal et al. 2016; Goddi et al. 2019). Other minor subtleties in the handling of polarimetric data are presented in Appendix A.

In Figure 1 (top row), we show the (u, v) coverage and low-band interferometric polarization of our main target M87 as a function of the baseline (u, v) after the initial calibration stage but before D-term calibration. The colors code the scan-averaged amplitude of the complex fractional polarization $\breve{m}$ (i.e., the fractional polarization in visibility space; for analysis of $\breve{m}$ in another source, Sgr A*, see Johnson et al. 2015). M87 is weakly polarized on most baselines, $| \breve{m}| \lesssim 0.5$. Several data points on SMA–SMT baselines have very high fractional polarization $| \breve{m}(u,v)| \sim 2$ that occur at (u, v) spacings where the Stokes ${ \mathcal I }$ visibility amplitude enters a deep minimum. The fractional polarization $\breve{m}$ of the M87 core is broadly consistent across the four days of observations and between low- and high-frequency bands, therefore high-band results are omitted in the display.

Zoom In Zoom Out Reset image size

In Figure 1 (bottom row) we show the field rotation angles ϕ for each station observing M87 on the four observing days. The data are corrected for this angle during the initial calibration stage, but the precision of the leakage calibration depends on how well this angle is covered and on the difference in the field angles at the two stations forming a baseline. In the M87 data the field rotation for stations forming long baselines (LMT, SMT, and PV) is frequently larger than 100° except for April 10, for which the (u, v) tracks are shorter.

In addition to the M87 data, a number of calibrators are utilized in this Letter for leakage calibration studies. To estimate D-terms for each of the EHT stations we use several EHT targets observed near in time to M87. In VLBI, weakly polarized sources are more sensitive to polarimetric calibration errors so they are preferred calibrators. For full-array leakage calibration, we focus on two additional sources: J1924–2914 and NRAO 530 (calibrators for the second EHT primary target, Sgr A*), which are compact and relatively weakly polarized. The main calibrator for M87 in total intensity, 3C 279 (Kim et al. 2020; Paper IV), is bright and strongly polarized on longer baselines and is not used in this work. The properties and analysis of the calibrators are discussed in more detail in Appendices J and K.

The closure traces for M87 and the calibrators can be used both to probe the data for uncalibrated systematic effects (see Appendix B.2) and to ascertain the presence of polarized flux density in a calibration-insensitive manner (see Appendix B.3).

Unless otherwise stated, the following analysis is focused on the low-band half of the data sets.

Producing an image of the linearly polarized emission requires both solving for the sky distribution of Stokes parameters ${ \mathcal Q }$ and ${ \mathcal U }$ and for the instrumental polarization of the antennas in the EHT array. In this work, we use several distinct methods to accomplish these tasks. Our approaches can be classified into three main categories: imaging via sub-component fitting; imaging via regularized maximum likelihood; and imaging as posterior exploration. In this section we only briefly describe each method; fuller descriptions are presented in Appendix C.

The calibration of the instrumental polarization by sub-component fitting was performed using three different codes (LPCAL, GPCAL, and polsolve) that depend on two standard software packages for interferometric data analysis: AIPS ¹³³ and CASA. ¹³⁴ In all of these methods, the Stokes ${ \mathcal I }$ imaging step is performed using the CLEAN algorithm (Högbom 1974), and sub-components with constant complex fractional polarization are then constructed from collections of the total intensity CLEAN components and fit to the data. In AIPS, two algorithms for D-term calibration are available: LPCAL (extensively used in VLBI polarimetry for more than 20 years; Leppänen et al. 1995) and GPCAL ¹³⁵ (Park et al. 2021). In CASA, we use the polsolve algorithm (Martí-Vidal et al. 2021), which uses data from multiple calibrator sources simultaneously to fit polarimetric sub-components and allows for D-terms to be frequency dependent (see Appendix D). In all sub-component fitting and imaging methods, we assume that Stokes ${ \mathcal V }=0$. Further details on LPCAL, GPCAL, and polsolve can be found in Appendix C.1.

Image reconstruction via the Regularized Maximum Likelihood (RML) method was used in Paper IV along with CLEAN to produce the first total intensity images of the 230 GHz core in M87. RML algorithms find an image that maximizes an objective function composed of a likelihood term and regularizer terms that penalize or favor certain image features. In this work, we use the RML method implemented in the eht-imaging ¹³⁶ software library (Chael et al. 2016, 2018) to solve for images in both total intensity and linear polarization. Like the CLEAN-based methods, eht-imaging does not solve for Stokes ${ \mathcal V }$. Details on the specific imaging methods in eht-imaging used in the reconstructions presented in this work can be found in Appendix C.2.

Imaging as posterior exploration is carried out using two independent Markov chain Monte Carlo (MCMC) schemes: D-term Modeling Code (DMC) and Themis. Both codes simultaneously explore the posterior space of the full Stokes image (including Stokes ${ \mathcal V }$) alongside the complex gains and leakages at every station; station gains are permitted to vary independently on every scan, while leakage parameters are modeled as constant in time throughout an observation. We provide more detailed model specifications for both codes in Appendix C.3 and in separate publications (Pesce 2021; A. E. Broderick et al. 2021, in preparation).

Hereafter, we often refer to eht-imaging, polsolve, and LPCAL methods as imaging methods/pipelines and to DMC and Themis methods as posterior exploration methods/pipelines.

In the imaging methods we divide the polarimetric calibration procedure for EHT data into two steps. In the first step, we calibrate the stations with an intra-site partner (ALMA–APEX, SMA–JCMT) using the assumption that sources are unresolved on intra-site baselines, where the brightness distribution can be approximated with a simple point source model. In the imaging pipelines we apply the D-terms for ALMA, APEX, and SMA to the data before polarimetric imaging and D-term calibration of the remaining stations. Baselines to the JCMT (which are redundant with SMA baselines) are removed from the data sets, to reduce complications from handling single-polarization data. The ALMA, APEX, and SMA D-terms are fixed in imaging with eht-imaging and polsolve; because LPCAL is unable to fix D-terms of specific stations to zero, it derives a residual leakage for these stations, which remains small. ¹³⁷ In the second step, we perform simultaneous imaging of the source brightness distribution and D-term calibration of stations for which only long, source-resolving baselines are available. In contrast, the posterior exploration pipelines do not use the D-terms derived using the intra-site baseline approach and instead solve for all D-terms (and station gains) starting with the base data product described in Section 3.2.

The point source assumption adopted in the imaging method intra-site baseline D-term calibration step is an extension to the intra-site redundancies already exploited in the EHT network calibration (Paper III), allowing us to obtain a model-independent gain calibration for ALMA, APEX, SMA, and JCMT. For an unresolved, slowly evolving source we can assume the true parameters of the coherency matrix ρ_jk in Equation (6) to be constant throughout a day of observations, as very low spatial frequencies u are sampled, ρ_jk ≈ ρ_jk ( u = 0). Hence, only four intrinsic visibility components of ρ_jk per source and four complex D-terms (two for each station) need to be determined from all the data on an available baseline.

We fit the D-terms of ALMA, APEX, JCMT, and SMA for each day using the multi-source feature of polsolve, combining band-averaged observations of multiple sources (3C 279, M87, J1924–2914, NRAO 530, 3C 273, 1055+018, OJ287, and Cen A as shown in Appendix D) on each day in one single fit per band. The results of these fits per station, polarization, day, and band are presented in Figure 2 (left panel), where we also plot the mean and standard deviation of the D-terms across all days and both bands for each station and polarization. In Appendix D, we provide tables with D-term values and further discuss the time and frequency dependence of D-terms and JCMT single polarization handling. In Appendix D we also present several validation tests of our intra-site baseline D-term estimation method carried out to motivate the use of band-averaged data products, comparisons to independent polarimetric source properties measured from simultaneous interferometric-ALMA observations (Goddi et al. 2021) near-in-time interferometric-SMA leakage estimates, and comparisons to results from a model fitting approach.

Zoom In Zoom Out Reset image size

In addition to intra-site baseline D-term calibration in the imaging pipelines, we also account for residual station-based amplitude gain errors by calibrating the data to pre-determined fiducial Stokes ${ \mathcal I }$ images of chosen calibrator sources. Given the extreme resolving power of the EHT array, all available calibrators are resolved on long baselines. Therefore, we must select sources that are best imaged by the EHT array; these are compact non-variable sources with sufficient (u, v) coverage. Four targets in the EHT 2017 observations fit these criteria: M87, 3C 279, J1924–2914, and NRAO 530. Stokes ${ \mathcal I }$ images of M87 and 3C 279 have been published (Papers I–VI; Kim et al. 2020). Final Stokes ${ \mathcal I }$ images for the Sgr A* calibrators NRAO 530 and J1924–2914 will be presented in upcoming publications (S. Issaoun et al. 2021, in preparation; S. Jorstad et al. 2021, in preparation) but the best available preliminary images are used to self-calibrate our visibility data for D-term comparisons in this Letter (Appendix J).

For M87, although multiple imaging packages and pipelines were utilized in the Stokes ${ \mathcal I }$ imaging process, the resulting final "fiducial" images from each method are highly consistent at the EHT instrumental resolution (e.g., Paper IV, Figure 15). Therefore, we selected a set of Stokes ${ \mathcal I }$ images for self-calibration from the RML-based SMILI imaging software pipeline (Akiyama et al. 2017a, 2017b; Paper IV). The images we use for self-calibration are at SMILI's native imaging resolution (∼10 μas), which provide the best fits to the data and are not convolved with any restoring beam. We self-calibrate our visibility data to these images, thereby accounting for residual station gain variations in the data that make imaging challenging. Using these self-calibrated data sets allows the imaging methods to focus on accurate reconstructions of the polarimetric Stokes ${ \mathcal Q }$ and ${ \mathcal U }$ brightness distributions and D-term estimation.

Preliminary D-terms estimated by the three imaging methods before testing and optimizing imaging parameters on synthetic data are reported in Appendix F. The right panels of Figure 2 show the final D-terms for LMT, PV, and SMT derived from the imaging and posterior modeling methods after optimization on synthetic data (see Section 4.3). To quantify the agreement (or distance in the complex plane) between D-term estimates from different methods we calculate L₁ norms. The L₁ norms averaged over left and right (also real and imaginary) D-term components, over all stations and over the four observing days, are all less than 1% for each pair of imaging methods (see Figure 20 in Appendix E). The mean values of the D-terms from the posterior exploration methods correlate well with the D-terms estimated by the imaging methods. For each combination of imaging and posterior exploration method the station averaged L₁ norms range from 1.5% to 1.89%.

DMC is the only method that solves for independent left- and right-circular station gains. The ratio of these gains derived by DMC on April 11 is shown in Figure 3. As expected, the assumption made by all of the imaging pipelines and one of the posterior exploration pipelines (Themis), that right- and left-hand gains are equal for all stations at all times, holds. For verification purposes, we also estimate D-terms using data of several calibrator sources. We find that the D-terms derived by polarimetric imaging of these other sources are consistent with those of M87 (Appendix J). Finally, we note our estimated SMT D-terms are similar to those computed previously using early EHT observations of Sgr A* (Johnson et al. 2015).

Zoom In Zoom Out Reset image size

Each imaging and leakage calibration method has free parameters that must be set by the user before the optimization or posterior exploration takes place. Some of these parameters (e.g., field of view, number of pixels) are common to all methods, but many are unique to each method (e.g., the sub-component definitions in LPCAL or polsolve, or the regularizer weights in eht-imaging). In VLBI imaging, these parameters are often simply set by the user given their experience on similar data sets, or based on what appears to produce an image that is a good fit to the data and free of noticeable imaging artifacts. In this work, we follow Paper IV in choosing the method parameters that we use in our final image reconstructions more objectively by surveying a portion of the parameter space available to each method.

We perform surveys over the different free parameters available to each method and attempt to choose an optimal set of parameters based on their performance in recovering the source structure and input D-terms from several synthetic data models. Appendix G provides more detail on the individual parameter surveys performed by each method. The parameter set that performs best on the synthetic data for each method is considered our "fiducial" parameter set for imaging M87 with that method. ¹³⁸ The corresponding images reconstructed from various data sets using these parameters are the method's "fiducial images."

The synthetic data sets that we used for scoring the imaging parameter combinations consist of six synthetic EHT observations using the M87 2017 April 11 equivalent low-band (u, v) coverage. The source structure models used in the six sets vary from complex images generated using general relativistic magnetohydrodynamic (GRMHD) simulations of M87's core and jet base (Models 1 and 2 from Chael et al. 2019) to simple geometrical models (a filled disk, Model 3, and simple rings with differing EVPA patterns, Models 4–6). The synthetic source models have varying degrees of fractional polarization and diverse EVPA structures. The synthetic source models blurred to the EHT nominal resolution are displayed in the first column of Figure 4.

Zoom In Zoom Out Reset image size

All M87 synthetic data sets were generated using the synthetic data generation routines in eht-imaging. We followed the synthetic data generation procedure in Appendix C.2 of Paper IV, but with models featuring complex polarization structure. The synthetic visibilities sampled on EHT baselines are corrupted with thermal noise, phase and gain offsets, and polarimetric leakage terms. Mock D-terms for the SMT, LMT, and PV stations were chosen to be similar to those found by the initial exploration of the M87 EHT 2017 data reported in Appendix F. Random residual D-terms for ALMA, APEX, JCMT, and SMA (reflecting possible errors in the intra-site baseline calibration procedure) were drawn from normal distributions with 1% standard deviation. After generation, the phase and amplitude gains in the synthetic data were calibrated for use in imaging pipelines in the same way as the real M87 data; that is, they were self-calibrated to a Stokes ${ \mathcal I }$ image reconstructed via the SMILI fiducial script for M87 developed in Paper IV.

In Figure 4, we present our fiducial set of images (in a uniform scale) from synthetic data surveys carried within each method. In each panel we report a correlation coefficient $\left\langle I\cdot {I}_{0}\right\rangle$ between recovered Stokes ${ \mathcal I }$ and the ground-truth ${{ \mathcal I }}_{0}$ images,

$$\begin{eqnarray}&&\langle I\cdot {I}_{0}\rangle =\displaystyle \frac{\langle \,({ \mathcal I }-\overline{{ \mathcal I }})({{ \mathcal I }}_{0}-{\overline{{ \mathcal I }}}_{0})\,\rangle }{\sqrt{\langle \,{\left({ \mathcal I }-\overline{{ \mathcal I }}\right)}^{2}\,\rangle }\sqrt{\langle \,{\left({{ \mathcal I }}_{0}-{\overline{{ \mathcal I }}}_{0}\right)}^{2}\,\rangle }}.\end{eqnarray} \tag{ 8 }$$

This reflects the dot product of the two mean-subtracted images when treated as unit vectors. We also calculate a correlation coefficient for the reconstructed linear polarization image ${ \mathcal P }\equiv { \mathcal Q }+i{ \mathcal U }$,

$$\begin{eqnarray} &&\langle \vec{P}\cdot {\vec{P}}_{0} \rangle =\displaystyle \frac{{\rm{Re}}[ \langle \,{{ \mathcal P }}_{}^{}{{ \mathcal P }}_{0}^{\ast }\, \rangle ]}{\sqrt{ \langle \,{{ \mathcal P }}^{}{{ \mathcal P }}^{\ast }\, \rangle }\sqrt{ \langle \,{{ \mathcal P }}_{0}^{}{{ \mathcal P }}_{0}^{\ast }\, \rangle }}.\end{eqnarray} \tag{ 9 }$$

The real part is chosen to measure the degree of alignment of the polarization vectors $({ \mathcal Q },{ \mathcal U })$. In both cases, images are first shifted to give the maximum correlation coefficient for Stokes ${ \mathcal I }$. Because Stokes ${ \mathcal I }$ image reconstructions are tightly constrained by an a priori known total image flux density, the Stokes ${ \mathcal I }$ correlation coefficients are mean subtracted to increase the dynamic range of the comparison. This introduces a field-of-view dependence to the metric, as only spatial frequencies above (field of view)⁻¹ are considered; up to the beam resolution. There is no such dependence in the linear polarization coefficient, which is not mean subtracted.

The correlation is equally strong independently of the employed method. The polarization structure is more difficult to recover for models with high or complex extended polarization (Models 1 and 2) for which correlation of the recovered polarization vectors is strong to moderate. In Figure 5 we present a uniform comparison of the recovered D-terms and the ground-truth D-terms for all synthetic data sets and methods. For all methods the recovered D-terms show a strong correlation with the model D-terms. To quantify the agreement (or distance in the complex plane) between D-term estimates and the ground-truth values D_Truth in each approach, we calculate the L₁ ≡ ∣D_i − D_Truth∣ norm, where D_i is a D-term component derived within a method i. Overall, for the fiducial set of parameters the agreement between the ground truth and the recovered D-terms in synthetic data measured using the L₁ norm is ≤1.3% on average (when averaging is done over stations, D-term components, and models). The reported averaged L₁ norms give us a sense of the expected discrepancies in D-terms between employed methods for their fiducial set of parameters. However, we notice again that the discrepancies do depend on source structure. For example, in models with no polarization substructure (e.g., Model 3) all methods had difficulty in recovering D-terms for PV (visible as large error bars for the station), a station forming only very long baselines on a short (u, v) track. If we exclude PV from the L₁ metrics the expected L₁ norms for LMT and SMT alone for all methods are L₁ ∼ 0.6%–0.8% when averaged over models.

Zoom In Zoom Out Reset image size

In Figure 6, we present the fiducial M87 linear-polarimetric images produced by each method from the low-band data on all four observing days. The fiducial images from each method are broadly consistent with those from the preliminary imaging stage shown in Figure 21 of Appendix F.

Zoom In Zoom Out Reset image size

Unless otherwise explicitly indicated, we display low-band results in the main text. The high-band results are given in Appendix I. We decided to keep the analysis of the high- and low-band data separate for several reasons. First, the main limitations in the dynamic range and image fidelity in EHT reconstructions arise from the sparse sampling of spatial frequencies, not the data S/N. Increasing the S/N by performing band averaging does not improve the dynamic range of the reconstructed images. Second, treating each band separately minimizes any potential chromatic effects that might add extra limitations to the dynamic range, such as intra-field differential Faraday rotation. Finally, separating the bands in the analysis allowed us to use the high-band results as a consistency check on the calibration of the instrumental polarization and image reconstruction for the low-band data. We perform this comparison of the results obtained at both bands in Appendix I. We conclude that both the recovered D-terms and main image structures are broadly consistent between the low and high bands.

The different reconstruction methods have different intrinsic resolution scales; for instance, the CLEAN reconstruction methods model the data as an array of point sources, while the RML and MCMC methods have a resolution scale set by the pixel size. In Figure 6, we display the fiducial images from each method at the same resolution scale by convolving each with a circular Gaussian kernel with a different FWHM. The FWHM for each method is set by maximizing the normalized cross-correlation of the blurred Stokes ${ \mathcal I }$ image with the April 11 "consensus" image presented in Figure 15 of Paper IV. The blurring kernel FWHMs selected by this method are 19 μas for eht-imaging, DMC, and Themis, 20 μas for LPCAL, and 23 μas for polsolve.

The M87 emission ring is polarized only in its southwest region and the peak fractional polarization at ≈20 μas resolution is at the level of about 15%. The residual rms in linear polarization (as estimated from the CLEAN images) is between 1.10–1.30 mJy/beam in all epochs, which implies a polarization dynamic range of ∼10. The nearly azimuthal EVPA pattern is a robust feature evident in all our reconstructions across time, frequency, and imaging method. The images show slight differences in the polarization structure between the first two days, 2017 April 5/6 and the last two, 2017 April 10/11. Notably, the southern part of the ring appears less polarized on the later days. This evolution in the polarized brightness is consistent with the evolution in the Stokes ${ \mathcal I }$ image apparent in the underlying closure phase data (Paper III, Figure 14; Paper IV, Figure 23). However, as with the Stokes ${ \mathcal I }$ image, the structural changes in the polarization images with time over this short timescale (6 days ≈16 GM/c³) are relatively small, and it is difficult to disentangle which differences in the polarized images are robust and which are influenced by differences in the interferometric (u, v) coverage between April 5 and 11 (Paper IV, Section 8.3).

In Figure 7, we show the simple average of the five equivalently blurred fiducial images (one per method) for each of the four observed days. The averaging is done independently for each Stokes intensity distribution. These method-averaged images are consistent with the EHT closure traces, as shown in Figure 13 in Appendix B. We adopt the images in Figure 7 as a conservative representation of our final M87 polarimetric imaging results.

Zoom In Zoom Out Reset image size

While the overall pattern of the linearly polarized emission from M87 is consistent from method to method, the details of the emission pattern can depend sensitively on the remaining statistical uncertainties in our leakage calibration. In addition, the different assumptions and parameters used in each reconstruction method affect the recovered polarized intensity pattern, introducing an additional source of systematic uncertainty in our recovered images. In this section, we assess the consistency of the recovered polarized images across different D-term calibration solutions within and between methods.

We explore the consistency of our image reconstructions against the uncertainties in the calibrated D-terms by generating a sample of 1000 images for each method, each generated with a different D-term solution. For the imaging methods, we define complex normal distributions for each D-term based on the scatter in recovered D-terms in Figure 2 and reconstruct images after calibrating to each set of random D-terms without additional calibration. This procedure is explained in detail in Appendix H. For the posterior exploration methods we simply draw 1000 images from the posterior for each observing day.

In each method's set of 1000 image samples covering a range of D-term calibration solutions, we study the azimuthal distribution of the polarization brightness (p) and EVPA (χ) by performing intensity-weighted averages of these quantities over different angular sections along the ring. The width of the angular sections used in the averaging is set to Δφ = 10° and the averages are computed from a position angle φ = 0° to φ = 360°, in steps of 1°.

Comparing angular averages of these quantities with a small moving window Δφ avoids spurious features from the different pixel scales used in the different image reconstruction methods. The pixel coordinates of the image center are estimated (for each method) from the peak of the cross-correlation between the ${ \mathcal I }$ images and the representative images of M87 used in the self-calibration. To avoid the effects of phase wrapping in the averaging (which biases the results for values of χ around ±90°), the quantity 〈χ〉 is computed coherently within each angular section, i.e., the averages are defined as

$$\begin{eqnarray} &&\langle { \mathcal P } \rangle \equiv \displaystyle \frac{ \langle { \mathcal I }\sqrt{{{ \mathcal Q }}^{2}+{{ \mathcal U }}^{2}} \rangle }{ \langle { \mathcal I } \rangle },\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\langle \chi \rangle \equiv \displaystyle \frac{1}{2}\arctan \left(\displaystyle \frac{\langle { \mathcal Q }\times { \mathcal I }\rangle }{\langle { \mathcal U }\times { \mathcal I }\rangle }\right).\end{eqnarray} \tag{ 11 }$$

In Figure 8, we show histograms of these quantities for two days, 2017 April 5 and 11, as a function of the orientation of the angular section used in the averaging (i.e., the position angle around the ring). We consider these two days because they have the best (u, v) coverage and span the full observation window; these results will thus include any effects of intrinsic source evolution in the recovered parameters. From Figure 8, it is evident that the difference in 〈p〉 between methods is larger than the widths of the 〈p〉 histograms in each method. This means that effects related to the residual instrumental polarization, giving rise to the dispersion seen in the histograms, are smaller than artifacts related to the deconvolution algorithms. In other words, the 〈p〉 images are limited by the image fidelity due to the sparse (u, v) coverage rather than by the D-terms.

Zoom In Zoom Out Reset image size

Even though there are differences among methods in the p azimuthal distribution, some features are common to all our image reconstructions. The peak in the polarization brightness is located near the southwest on 2017 April 5 (at a position angle of 199° ± 11°, averaged among all methods) and close to the west on 2017 April 11 (position angle of 244° ± 10°). That is, the polarization peak appears to rotate counter-clockwise between the two observing days (see the dotted lines in Figure 8). On both days, the region of high polarization brightness is relatively wide, covering a large fraction of the southern portion of the image (position angles from around 100°–300°).

In the azimuthal distribution of 〈χ〉, all methods produce very similar values in the part of the image with the highest polarized brightness (the southwest region, between position angles of 180° and 270°). The EVPA varies almost linearly, from around 〈χ〉 = − 80° (in the south) up to around 〈χ〉 = 30° (in the east). The EVPAs on 2017 April 11 are slightly higher (i.e., rotated counter-clockwise) compared to those on 2017 April 5. This difference is clearly seen for eht-imaging, polsolve, and Themis, though the difference is smaller for DMC and LPCAL. We notice, though, that the differences in the EVPAs between days could also be affected by small shifts in the estimates of the image center on each day. Outside of the region with high polarization, the EVPA distributions for all methods start to depart from each other. There is a hint of a constant EVPA 〈χ〉 ∼ 0° in the northern region (i.e., position angles around 0°–50°) in polsolve and LPCAL on both days, but the other methods show larger uncertainties in this region.

The discrepancies in EVPA among all methods only appear in the regions with low brightness (i.e., around the northern part of the ring). Therefore, polarization quantities defined from intensity-weighted image averages, discussed in the next sections, will be dominated by the regions with higher brightness, for which all methods produce similar results. Image-averaged quantities are somewhat more robust to differences in the calibration and image reconstruction algorithms, though they are not immune to systematic errors.

In comparing polarimetric images of M87, we are most interested in identifying acceptable ranges of three image-averaged parameters that are used to distinguish between different accretion models in Paper VIII: the net linear polarization fraction of the image ∣m∣_net, the average polarization fraction in the resolved image at 20 μas resolution 〈∣m∣〉, and the m = 2 coefficient of the azimuthal mode decomposition of the polarized brightness β₂. These parameters are defined below.

First, the net linear polarization fraction of the image is

$$\begin{eqnarray}&&| m{| }_{\mathrm{net}}=\displaystyle \frac{\sqrt{{\left({\displaystyle \sum }_{i}{Q}_{i}\right)}^{2}+{\left({\displaystyle \sum }_{i}{U}_{i}\right)}^{2}}}{{\displaystyle \sum }_{i}{I}_{i}},\end{eqnarray} \tag{ 12 }$$

where the sum is over the pixels indexed by i. ALMA measured ∣m∣_net = 2.7% on 2017 April 11 (Goddi et al. 2021), but this measurement includes emission at large scales outside of the 120 μas field of view of the EHT images. We also consider the intensity-weighted average polarization fraction across the resolved EHT image:

$$\begin{eqnarray}&&\langle | m| \rangle =\displaystyle \frac{{\displaystyle \sum }_{i}\sqrt{{Q}_{i}^{2}+{U}_{i}^{2}}}{{\displaystyle \sum }_{i}{I}_{i}}.\end{eqnarray} \tag{ 13 }$$

The value of 〈∣m∣〉 is determined by the intensity of the polarized emission at each point in the image, and thus it is sensitive to the resolution of the image and the choice of restoring beam. Specifically, images restored with beams of larger FWHM will tend to be more locally depolarized and thus have lower 〈∣m∣〉 than images restored with beams of smaller FHWM. In contrast, the integrated polarization fraction ∣m∣_net is insensitive to convolution.

We quantify the polarization structure with a decomposition into azimuthal modes. In particular, Paper VIII considers the complex amplitude β₂ of the m = 2 mode defined in Palumbo et al. (2020), who found this mode to be the most important in distinguishing different modes of accretion from 230 GHz images produced by different GRMHD simulations. The β₂ azimuthal mode decomposition coefficient is defined as

$$\begin{eqnarray}&&{\beta }_{2}=\displaystyle \frac{1}{{I}_{\mathrm{ring}}}{\int }_{{\rho }_{\min }}^{{\rho }_{\max }}{\int }_{0}^{2\pi }P(\rho ,\varphi )\,{e}^{-2i\varphi }\ \rho d\varphi d\rho ,\end{eqnarray} \tag{ 14 }$$

where (ρ, φ) are polar coordinates in the image plane, and I_ring is the Stokes ${ \mathcal I }$ flux density in the ring between the minimum radius ${\rho }_{\min }$ and the maximum radius ${\rho }_{\max }$. Because our reconstruction methods recover no significant extended brightness off the main ring, we take ${\rho }_{\min }=0$ and extend ${\rho }_{\max }$ to encompass the full-image field of view, so I_ring is equal to the total Stokes ${ \mathcal I }$ flux density in the image.

Note that both the amplitude ∣β₂∣ and phase ∠β₂ depend on the choice of image center, image resolution, and restoring beam size. In the comparisons that follow, we convolve images from each method with circular Gaussian beams with the FWHMs specified in Section 5.1 chosen to bring all images to the same resolution scale. Furthermore, we center each reconstruction by finding the pixel offset that maximizes the cross-correlation between the blurred Stokes ${ \mathcal I }$ image and the 2017 April 11 consensus Stokes ${ \mathcal I }$ image from Paper IV. In general, we find these offsets to be small, and our results do not change significantly if we do not apply any centering procedure in calculating β₂ from our reconstructed images.

From the sets of 1000 images generated by each method to explore variations of the image structure with the D-term solution, we compute distributions of each key metric—∣m∣_net, 〈∣m∣〉, ∣β₂∣, and ∠β₂—that is used in Paper VIII for theoretical interpretation. These distributions are summarized in Figure 9, which displays the mean points and 1σ error bars from different D-term realizations for all four methods on both 2017 April 5 and 11. We present a more complete look at these distributions with histograms for each quantity from each method/day pair in Appendix H, Figures 25 and 26. Figure 9 shows results for low-band images only; we compare these results to results derived from the high-band data in Figure 29 in Appendix I. Because we derived and vetted our imaging procedures for the low-band data, we use only the low-band results in determining our final parameter measurements and use the high-band results in Appendix I as a consistency check.

Zoom In Zoom Out Reset image size

On each observation day, the distributions of ∣m∣_net appear consistent between most pairs of reconstruction methods, with some notable exceptions. Many of the distributions of ∣m∣_net peak around the ALMA measured value of 2.7%, but the LPCAL distributions on both days and the eht-imaging distributions on 2017 April 11 are peaked closer to 1%. The distributions of 〈∣m∣〉 are peaked between 6% and 11% for all five methods across both days. On both days, the 〈∣m∣〉 distributions for eht-imaging, DMC, and Themis are peaked at values 2%–3% higher than the corresponding LPCAL or polsolve distributions. This systematic shift may indicate residual issues with bringing the reconstruction methods to the same resolution scale; in particular, the same circular Gaussian kernel was used to blur Stokes ${ \mathcal I },{ \mathcal Q }$ and ${ \mathcal U }$ in each method, while the intrinsic resolution of the reconstruction in ${ \mathcal Q }$ and ${ \mathcal U }$ may be lower than in total intensity. In each method, there appears to be a decrease in 〈∣m∣〉 of ≈1%–2% between 2017 April 5 and 11. Note that, because it is constrained to be positive and defined as the ratio of two uncertain flux densities, the mean of a distribution of fractional polarization m can have a positive bias. In Figure 25, we see that the distributions of ∣m∣_net and 〈∣m∣〉 both can have long tails; this is most evident on 2017 April 10, when the image reconstructions are the most uncertain due to poor (u, v) coverage. On 2017 April 5 and 11, we do not see prominent tails in the distributions of ∣m∣_net and 〈∣m∣〉. Furthermore, we expect any bias in the mean of these quantities in the measurement from a single method to be overwhelmed by the systematic uncertainty between different reconstruction methods.

The mean of the ∣β₂∣ distribution is peaked between 0.04 and 0.07 for all methods on both days; however, the ∣β₂∣ distributions from eht-imaging, DMC, and Themis have larger mean values on both days than the corresponding distributions for polsolve and LPCAL. Again, because ∣β₂∣ is sensitive to the restoring beam size, this may be due to residual errors in bringing the polarized images to the same resolution scale. Similarly to the distributions of 〈∣m∣〉, there are indications of a shift downward in ∣β₂∣ by an absolute value of ≈0.01 in all four methods between 2017 April 5 and 11. The distributions of the phase ∠β₂ are consistent between most pairs of methods with no obvious systematic difference between the sub-component methods (LPCAL, polsolve) and those that use a continuous image representation (eht-imaging, DMC, and Themis). Furthermore, there is no apparent systematic difference in the ∠β₂ results between 2017 April 5 and 11.

To score different accretion models from GRMHD simulations against constraints from the EHT data, Paper VIII uses a range for each quantity that incorporates both the uncertainties in the parameters from the D-term calibration process (the error bars for each method in Figure 9) and the systematic uncertainty across methods (the scatter in the points across methods). The final ranges used in Paper VIII for each parameter were set by taking the minimum/maximum of the 10 mean values minus/plus the 1σ error bars across both days and all methods. These parameter ranges are denoted by colored bands in Figure 9 and are presented in Table 2.

Phased ALMA records the VLBI signals in a basis of linear polarization, which need a special treatment after the correlation (Martí-Vidal et al. 2016; Matthews et al. 2018). The post-correlation conversion of the ALMA data from a linear basis into a circular basis has implications for the kind of instrumental polarization left after fringe fitting. As discussed in Goddi et al. (2019), any offset in the estimate of the phase difference between the X and Y signals of the ALMA antenna used as the phasing reference (an offset likely related to the presence of a non-zero Stokes ${ \mathcal V }$ in the polarization calibrator) maps into a post-conversion polarization leakage that can be modeled as a symmetric, pure-imaginary D-term matrix (i.e., D_R = D_L = iΔ). The amplitude of the ALMA D-terms, Δ, can be approximated (to a first order) as the value of the phase offset between X and Y in radians (Goddi et al. 2019). Hence, we expect the D_R and D_L estimates for ALMA to be found along the imaginary axis and to be of similar amplitude.

Furthermore, the ALMA feeds in Band 6 (the frequency band used in the EHT observations) are rotated by 45° with respect to their projection on the focal plane. This introduces a phase offset between the RCP and LCP post-converted signals that has to be corrected after the fringe fitting. This offset can be applied as a global phase added (subtracted) to the RL^* (LR^*) correlation products in all baselines (because ALMA has been used as the reference antenna in the construction of the global fringe-fitting solutions). We have applied this 45° rotation to all the visibilities before performing the analysis described in this Letter. Hence, the absolute position angles of the electric vectors (EVPA) derived from our EHT observations are properly rotated into the sky frame. This property of the ALMA–VLBI observations (see Appendix D) gives us absolute EVPA values instantaneously.

The LMT shows an unexpectedly high leakage signal with a large delay of ∼1.5 ns, which affects the cross-polarization phase spectra of the baselines related to the LMT. As a consequence, all the baselines related to the LMT show spurious instrumental fringes in the RL^* and LR^* correlations, with amplitudes similar to (and even higher than, for the case of sources with low intrinsic polarization) that of the main fringe. These instrumental fringes are smallest in the parallel-hand correlations (RR^* and LL^*), but relatively high in the cross-polarization hands and are related to strong polarization leakage likely due to reflections in the optical setup of the LMT receiver used in 2017 (Paper III). For the EHT observations on year 2018 and beyond, the special-purpose interim receiver used at the LMT was replaced by a dual-polarization sideband-separating 1.3 mm receiver, with better stability and full 64 Gbps sampling as for the rest of the EHT (Paper II), so future polarimetry analyses of the EHT may be free of this instrumental effect from the LMT.

If we take the frequency average over all intermediate-frequency (IF) sub-bands (the results presented in Paper I–Paper VI are based on this averaging), the effect of this leaked fringe is smeared out, as the average is equivalent to taking the value of the visibility at the peak of the main fringe. This main peak is only affected by the sidelobe of the delayed leaked fringe, with a relative amplitude that we estimate to be of 10%–20% of the cross-polarization main fringe. Therefore, the effect of the leaked fringe is small in comparison to the contribution from the ordinary instrumental polarization, which can especially dominate the cross-polarization signal for observations of sources with low polarization like M87, and can be ignored.

The dual-polarization observations performed by the SMA use two independent receivers at each antenna to register the RCP and LCP signals. However, the visibility matrices of the baselines related to the SMA are built from the combination of the RCP and LCP streams as if they were registered with one single receiver. Therefore, some of the assumptions made in the RIME (see Equation (6)) for the polarimetry calibration (e.g., stable relative phases and amplitude between polarizations) may not apply for the SMA-related visibilities. However, the fringe fitting of the parallel-hand correlations related to the SMA, as well as the absolute amplitude calibration (both described in Paper III) did account for the drifts in cross-polarization phase and amplitude between the SMA receivers, which makes it possible to model the instrumental polarization using ordinary leakage matrices.

One extra correction that has to be applied to the D-terms of the SMA is a phase rotation between the RCP and LCP leakages, to account for the 45° rotation of the antenna feed with respect to the mount axes. The D-terms shown in Section 4.2 and in Appendix D are corrected for this rotation.

The JCMT was equipped with a single-polarization receiver for these observations, so that only one of the two polarizations can be used at each epoch. Therefore, only one of the two cross-polarization correlations can be computed in all the baselines related to the JCMT; depending on which product is computed, we can only solve for one of the two D-terms of the JCMT (i.e., D_L if RCP is recorded; D_R otherwise).

As explained in Martí-Vidal et al. (2016), a byproduct of the use of polconvert in VLBI is the calibration of the absolute cross-polarization delays and phases in the stations with polconverted data, which allow for the reconstruction of the absolute EVPAs of the observed sources. The only condition to have this absolute R/L delay and phase calibration is to use the polconverted station (i.e., ALMA, in the case of the EHT) as the reference antenna in the fringe fitting.

In Figure 11, we show the difference of multi-band delays between RL^* and LR^* after the GFF calibration described in Paper III. All source scans and baselines with an S/N higher than 10 are shown for times when ALMA was participating in the observations. According to Martí-Vidal et al. (2016), the delay difference between RL^* and LR^* should be around zero when ALMA is the reference antenna. We see, though, hints of a small global residual delay difference after the GFF calibration (the points are not symmetrically distributed around zero). The weighted average of all the delay differences shown in Figure 11 is Δτ = − 28 ± 1 ps. This is a very small delay in absolute value (the amplitude losses due to this delay in each correlation product is lower than 1%), but still detectable at the remarkably high S/N level of the EHT observations.

Zoom In Zoom Out Reset image size

Closure traces (Broderick & Pesce 2020) are calibration-insensitive quantities constructed on station quadrangles from the coherency matrices ρ_jk defined in Equation (2):

$$\begin{eqnarray}&&{{ \mathcal T }}_{{ijkl}}=\displaystyle \frac{1}{2}\mathrm{tr}\left({\rho }_{{ij}}{\rho }_{{kj}}^{-1}{\rho }_{{kl}}{\rho }_{{il}}^{-1}\right).\end{eqnarray} \tag{ B1 }$$

These data products are a superset of the more familiar closure quantities (closure phases and closure amplitudes), with the additional property that they are independent of instrumental polarization. The closure traces are also independent of any other station-based effects that can be described in a Jones matrix formalism, including the definition of the polarization basis (e.g., the representation of the polarized quantities in terms of linear or circular feeds). The closure traces thus provide a powerful tool with which to make calibration-independent statements regarding polarimetric data and intrinsic source structure.

By analogy with trivial closure phases (see Paper III), trivial closure traces may be constructed on "boomerang" quadrangles, i.e., quadrangles in which a station is effectively repeated in such a way as to make the quadrangle area vanish (Broderick & Pesce 2020). Given two co-located stations i and $i^{\prime}$, the closure trace ${{ \mathcal T }}_{{iji}^{\prime} k}$ reduces to unity.

Each quadrangle ijkl has an associated "conjugate" quadrangle ilkj, constructed by reordering the baselines within the coherency matrix product. ¹³⁹ Conjugate closure trace products can be expressed as

$$\begin{eqnarray}\,\begin{array}{rcl}{{ \mathcal C }}_{{ijkl}} & \equiv & {{ \mathcal T }}_{{ijkl}}{{ \mathcal T }}_{{ilkj}}=1+\ {\left({\breve{q}}_{{ij}}-{\breve{q}}_{{kj}}+{\breve{q}}_{{kl}}-{\breve{q}}_{{il}}\right)}^{2}\\ & & \,+\ {\left({\breve{u}}_{{ij}}-{\breve{u}}_{{kj}}+{\breve{u}}_{{kl}}-{\breve{u}}_{{il}}\right)}^{2}\\ & & \,+\ {\left({\breve{v}}_{{ij}}-{\breve{v}}_{{kj}}+{\breve{v}}_{{kl}}-{\breve{v}}_{{il}}\right)}^{2}\\ & & \,+\ { \mathcal O }\left({\breve{q}}^{3},{\breve{u}}^{3},{\breve{v}}^{3}\right),\end{array}\,\end{eqnarray} \tag{ B2 }$$

where ${\breve{q}}_{{ij}}\equiv {\tilde{{ \mathcal Q }}}_{{ij}}/{\tilde{I}}_{{ij}}$, ${\breve{u}}_{{ij}}\equiv {\tilde{{ \mathcal U }}}_{{ij}}/{\tilde{I}}_{{ij}}$, and ${\breve{v}}_{{ij}}\equiv {\tilde{{ \mathcal V }}}_{{ij}}/{\tilde{I}}_{{ij}}$. The ${{ \mathcal C }}_{{ijkl}}$ are identically unity in the absence of intrinsic source polarization and for point sources. Deviations from unity require non-constant interferometric polarization fractions on baselines in the quadrangle ijkl, and therefore closure trace products are a robust indicator of polarized source structures (Broderick & Pesce 2020).

For EHT observations, boomerang quadrangles are formed using the redundant baselines presented by ALMA and APEX. ¹⁴⁰ In the top panel of Figure 12, the phases of all of the trivial closure traces are shown for M87 (blue) and the calibrators (J1924–2914, NRAO 530, 3C 279; gray), constructed from scan-averaged visibility data. The values of these phases are clustered about zero, consistent with the expectation that the trivial closure traces are unity.

Zoom In Zoom Out Reset image size

The distribution of the normalized residuals provide a direct assessment of the systematic error budget of the polarimetric data independent of the gain and leakage calibration. These residuals are shown in the bottom panel of Figure 12 for both M87 and calibrators. We find that the data match the anticipated unit-variance Gaussian, which is consistent with an absence of unidentified systematic uncertainties in the polarimetric data.

Figure 13 shows the phase of the conjugate closure trace product, ${{ \mathcal C }}_{{ijkl}}$, for a quadrangle pair ALMA–APEX–LMT–SMT and ALMA–SMT–LMT–APEX. The presence of non-zero ${{ \mathcal C }}_{{ijkl}}$ is a calibration-insensitive indicator of significant polarized structures in the Stokes map. Because the uncertainties of the closure traces on conjugate quadrangles are correlated, the resulting uncertainty in the conjugate closure trace product is typically smaller than would be estimated from assuming independent uncertainties in the individual closure traces. The errors shown have been estimated using Monte Carlo sampling of the constituent visibilities.

Zoom In Zoom Out Reset image size

Closure trace phases are shown on a handful of non-trivial quadrangles in Figure 14 for M87. These phases are clearly non-zero and exhibit variations throughout the observing night, which is consistent with non-trivial source structure. The behavior of the closure trace evolution is similar across neighboring observation days (e.g., 2017 April 5/6, April 10/11) and consistent between quadrangles constructed using ALMA (filled markers) and APEX (open markers). The behavior of the closure trace evolution is dissimilar between the 2017 April 5/6 and April 10/11 observations, providing direct evidence for an evolving source structure in M87 independent of all station-based corrupting effects, including polarization leakage. Because ALMA, LMT, and SMT are nearly co-linear as seen from M87 for much of the observations, and because the closure traces are presumably tracing primarily the Stokes ${ \mathcal I }$ emission, the closure trace phases in Figure 14 are very similar to the Stokes ${ \mathcal I }$ closure phases shown in Figure 14 of Paper III. Note that the points plotted in Figure 14 have been averaged across both the high and low bands.

Zoom In Zoom Out Reset image size

Conjugate closure trace product phases are shown in Figure 15 for each observation day for the ALMA–PV–LMT–SMT and ALMA–SMT–LMT–PV quadrangle pair. There is the appearance of temporal evolution from April 5/6 to April 10/11, with an attendant implication for an evolution in the polarization map of M87 between those periods. However, the paucity of quadrangles exhibiting significant evolution renders this conclusion suggestive at best. Note that the points plotted in Figure 15 have been averaged across both the high and low bands.

Zoom In Zoom Out Reset image size

In the sub-component fitting method for polarimetric modeling, the instrumental polarization (i.e., the complex D-terms) and the source polarized brightness distribution are estimated simultaneously from the interferometric observables and a fixed estimate of the source brightness distribution ${ \mathcal I }({\boldsymbol{x}})$. The sub-component fitting calibration algorithms estimate the D-terms from Equations (4) and (7) by modeling the polarized source structure as a disjoint set of N "polarization sub-components," ${{ \mathcal I }}_{i}({\boldsymbol{x}})$, such that:

$$\begin{eqnarray}&&{ \mathcal I }({\boldsymbol{x}})=\displaystyle \sum _{i}^{N}{{ \mathcal I }}_{i}({\boldsymbol{x}}).\end{eqnarray} \tag{ C1 }$$

The fractional polarization of each sub-component is assumed to be constant, implying that:

$$\begin{eqnarray}\,\begin{array}{rcl}{ \mathcal Q }({\boldsymbol{x}}) & = & \displaystyle \sum _{i}^{N}{q}_{i}\,{{ \mathcal I }}_{i}({\boldsymbol{x}})\\ { \mathcal U }({\boldsymbol{x}}) & = & \displaystyle \sum _{i}^{N}{u}_{i}\,{{ \mathcal I }}_{i}({\boldsymbol{x}}),\end{array}\,\end{eqnarray} \tag{ C2 }$$

where q_i and u_i are real-valued constants. In the sub-component fitting method, we therefore assume that the polarized brightness is exactly proportional to ${{ \mathcal I }}_{i}$ for each source sub-component. This condition is known as the "similarity approximation" and may produce inaccurate estimates of the instrumental polarization for cases of strongly polarized and resolved calibrators (e.g., Cotton 1993) and/or if the sub-division of ${ \mathcal I }$ into sub-components is not performed properly. Discussions about the self-similarity assumption can be found in Appendix K.

The polsolve, LPCAL, and GPCAL algorithms determine which values of q_i , u_i , D_R , and D_L minimize the difference between the calibrated visibility matrix and the Fourier-transformed model brightness matrix (Equations (6) and (4)). The total number of parameters used in this fit is equal to two times the the number of source sub-components (i.e., 2N, which correspond to q_i and u_i in Equation (C2)) plus four times the number of antennas (i.e., 4N_a , accounting for the real and imaginary parts of the D_R and D_L of each antenna). The error function to be minimized is the sum of the χ² values computed for the cross-polarization matrix elements of the RIME, i.e.,

$$\begin{eqnarray}\begin{array}{c}\begin{array}{c}\,\begin{array}{ccc}{\chi }^{2} & = & \displaystyle \sum _{m=1}^{{N}_{v}}{w}_{m}{\left|{R{L}^{\ast }}_{{kl},m}^{c}-{\left(\tilde{{ \mathcal Q }}+i\tilde{{ \mathcal U }}\right)}_{m}\right|}^{2}\\ & & +\,\displaystyle \sum _{m=1}^{{N}_{v}}{w}_{m}{\left|{L{R}^{\ast }}_{{kl},m}^{c}-{\left(\tilde{{ \mathcal Q }}-i\tilde{{ \mathcal U }}\right)}_{m}\right|}^{2},\end{array}\,\end{array}\end{array}\end{eqnarray} \tag{ C3 }$$

where w_m is the weight of the mth visibility, the index c stands for calibrated visibilities (corrected both for station gains and for instrumental polarization using the current estimate of the D-terms) and N_v is the number of visibilities.

The calibrated visibilities, ${R{L}^{\ast }}_{{kl},m}^{c}$ and ${L{R}^{\ast }}_{{kl},m}^{c}$, depend on ${D}_{R}^{k}$, ${D}_{L}^{k}$, ${D}_{R}^{l}$ and ${D}_{L}^{l}$ (Equation (4)), whereas $\tilde{{ \mathcal Q }}$ and $\tilde{{ \mathcal U }}$ depend on q_i and u_i . The χ² minimization solves for the intrinsic source Stokes parameters and the instrumental polarization simultaneously. We note that the effects of instrumental polarization are constant in the frame of the antenna feed for Cassegrain-mounted feeds, whereas the intrinsic source polarization is defined in the sky frame; as a consequence, the changing feed angle of each antenna across the observations (i.e., the Earth rotation during the extent of the observations) allows the model fitting to decouple the antenna D-terms from the Stokes parameters of the source sub-components. For Nasmyth-mounted feeds, there is an additional rotation between cross-polarization introduced by the feed/optics and the telescope itself. In this case, we assume a minimal contribution from the antennas themselves, a reasonable assumption for these on-axis telescopes. Equation (C3) implies that there are several implicit assumptions in the polarimetric modeling of polsolve and LPCAL. First, RL*^c and LR*^c are computed by setting ${ \mathcal V }=0$ (i.e., any circular polarization in the calibrators is neglected, compared to Stokes ${ \mathcal I }$). Furthermore, the real and the imaginary parts of the residual visibilities (i.e., either RL*^c or LR*^c minus the Fourier transforms of the corresponding model brightness distributions) are assumed to be statistically independent.

If the linear polarization structures of calibrators are not similar to their total intensity structures, a breakdown of the similarity approximation can occur. This can be a source of uncertainties in D-term estimation. In Appendix K, we discuss this effect on our results for different polarization calibrators reported in this Letter.

The package eht-imaging (Chael et al. 2016, 2018) implements polarimetric image reconstruction via RML. eht-imaging solves for an image X by minimizing an objective function via gradient descent. The objective function J( X ) is a weighted sum of data-consistency log-likelihood terms and regularizer terms that favor or penalize certain image features. That is, to find an image (in either total intensity, polarization, or both) we minimize

$$\begin{eqnarray}&&J({\boldsymbol{X}})=\displaystyle \sum _{\mathrm{data}\,\mathrm{terms}\,i}{\alpha }_{i}{\chi }_{i}^{2}({\boldsymbol{X}})+\displaystyle \sum _{\mathrm{regularizers}\,j}{\beta }_{j}{S}_{j}({\boldsymbol{X}}).\end{eqnarray} \tag{ C4 }$$

Picking optimal values of the "hyperparameter" weights α_i and β_j in Equation (C4) is an essential task in RML imaging. Here we describe the data terms and regularizers that we use for polarimetric imaging, and in Appendix G we describe our method for determining the hyperparameters using parameter surveys.

For polarized image reconstructions, we follow the method laid out in Chael et al. (2016), with the addition of iterative self-calibration of any uncorrected station D-terms. We start with data that has had the overall time-dependent station amplitude and phase gains calibrated using the SMILI fiducial image from Paper IV, and the ALMA, APEX, SMA, and JCMT D-terms have been corrected using the zero-baseline solutions described in Section 4.2. The data are scan-averaged. We then reconstruct the Stokes ${ \mathcal I }$ using the same fiducial imaging script for eht-imaging developed in Paper IV. We fix the image field of view at 120 μas and solve for a grid of 64 × 64 pixels. In the Stokes ${ \mathcal I }$ imaging, the total flux density is constrained to be 0.6 Jy. We next (re)self-calibrate the station amplitude and phase gains (assuming G_R = G_L ) to our final Stokes I image. Having extensively explored the imaging parameter space for Stokes ${ \mathcal I }$ imaging in Paper IV, we do not vary these parameters in our polarimetric imaging surveys. After self-calibrating to our final total intensity image, we drop zero baselines for the polarimetric imaging stage.

In defining an objective function of the form in Equation (C4) for the polarized image reconstruction, we consider two log-likelihood χ² terms; one computed using the RL* polarimetric visibility ${ \mathcal P }=\tilde{{ \mathcal Q }}+i\tilde{{ \mathcal U }}$, and one using the visibility domain polarimetric ratio $\breve{m}={ \mathcal P }/\tilde{{ \mathcal I }}$. ${\chi }_{\breve{m}}^{2}$ is immune to any residual station gain error left over from Stokes ${ \mathcal I }$ imaging, while ${\chi }_{\tilde{p}}^{2}$ is not. We use two regularizers on the polarized intensity. First, the Holdaway-Wardle (Holdaway & Wardle 1990) regularizer S_HW (Equation (13) of Chael et al. 2016) acts like an entropy term that prefers image pixels take a value less than ${m}_{\max }=0.75$. This regularizer encourages image pixels to stay below the theoretical maximum polarization fraction for synchrotron radiation, but it is not a hard limit. Second, the total variation (TV) regularizer S_TV (Rudin et al. 1992) acts to minimize pixel-to-pixel image gradients in both the real and imaginary part of the complex polarization brightness distribution (Equation (15) of Chael et al. 2016).

Taken together, the objective function we minimize in polarimetric imaging is

$$\begin{eqnarray}&&{J}_{\mathrm{pol}}({ \mathcal Q },{ \mathcal U })={\alpha }_{p}{\chi }_{\tilde{p}}^{2}={\alpha }_{m}{\chi }_{\breve{m}}^{2}-{\beta }_{\mathrm{HW}}{S}_{\mathrm{HW}}-{\beta }_{\mathrm{TV}}{S}_{\mathrm{TV}}.\end{eqnarray} \tag{ C5 }$$

The relative weighting between the data constraints and the regularizer terms is set by the four hyperparameters α_P , α_m , β_HW, and β_TV.

We solve for the polarized intensity distribution that minimizes Equation (C5) parameterized by the fractional polarization m and EVPA ξ in each pixel. The Stokes ${ \mathcal I }$ image is fixed in the polarimetric imaging step and defines the region where polarized emission is allowed. To ensure our solution respects ${{ \mathcal Q }}^{2}+{{ \mathcal U }}^{2}\lt {{ \mathcal I }}^{2}$ everywhere, we transform the fractional polarization m in each pixel from the range m ∈ [0, 1] to κ ∈ (−∞ , ∞ ) and solve for κ (See Appendix D of Chael et al. 2016). In the eht-imaging script for EHT M87 observations, we solve for the pixel values of κ and ξ that minimize the objective function by gradient descent, and we then transform $({ \mathcal I },\kappa ,\xi )\to ({ \mathcal I },{ \mathcal Q },{ \mathcal U })$. We often restart the gradient descent process several times, using the output of the previous round of imaging blurred by a 20 μas Gaussian kernel as the new initial point.

In between rounds of polarimetric imaging with eht-imaging, we iteratively solve for the remaining D-terms by minimizing the χ² between the real (gain-calibrated) data and sampled data from the current image reconstruction corrupted with Jones matrices (Equation (4)). We do not use any linearized approximations of the effects of the Jones matrices when solving for the D-terms, but throughout we assume the model image has no circular polarization (${ \mathcal V }=0$). The eht-imaging pipeline thus alternates between rounds of polarimetric imaging and D-term calibration; often it takes many successive rounds of imaging and D-term calibration (n_iter ≈ 50–100) for the process to converge on a stable D-term solution.

In this section we describe two MCMC schemes developed for polarimetric imaging. Both MCMC codes model the polarized emission structure on a Cartesian grid of intensity points, with the Stokes vector parameterized using a spherical (Poincaré) representation,

$$\begin{eqnarray}&&\left(\begin{array}{c}{{ \mathcal I }}_{i}\\ {{ \mathcal Q }}_{i}\\ {{ \mathcal U }}_{i}\\ {{ \mathcal V }}_{i}\end{array}\right)={{ \mathcal I }}_{i}\left(\begin{array}{c}1\\ {{\ell }}_{i}\cos ({\xi }_{i})\sin ({\varsigma }_{i})\\ {{\ell }}_{i}\sin ({\xi }_{i})\sin ({\varsigma }_{i})\\ {{\ell }}_{i}\cos ({\varsigma }_{i})\end{array}\right),\end{eqnarray} \tag{ C6 }$$

where the index i runs over individual grid points, where we solve for the total intensiy ${{ \mathcal I }}_{i}$, the fractional polarization ℓ_i , and the two angles ξ_i , ς_i . Stokes visibilities are generated from the gridded emission structure via a direct Fourier transform (i.e., treating each grid point as a point source), and the visibilities are then multiplied with a smoothing kernel to impose image continuity. The parallel- and cross-hand visibilities on each baseline are then computed from the Stokes visibilities using Equation (7), and the gains and leakage terms are applied to the model visibilities using a Jones matrix formalism (see Equation (6)). The model and data visibilities are ultimately compared via complex Gaussian likelihood functions for each of the parallel- and cross-hand data products independently, with the total likelihood taken to be the product of the individual likelihoods for all parallel- and cross-hand data products.

C.3.1. DMC

C.3.2.  Themis

C.3.3.  Themis−DMC Leakage Posterior Comparison

In Figure 16 we show a comparison between the leakage posteriors for all stations, as determined by DMC and Themis fits to the 2017 April 11 observations of M87. Despite the various different assumptions and model specifications, we find excellent agreement in both the means and shapes of the posterior distributions recovered from both methods. The modest discrepancies between the posteriors shown in Figure 16 are associated with the different treatment of systematic uncertainty and right/left gain ratios between the two methods; when these model choices are homogenized, the DMC and Themis fits to both synthetic data sets and to the M87 data return indistinguishable posteriors. Notably, both model treatments of the Stokes map appear to be comparably capable of capturing the source structure.

Zoom In Zoom Out Reset image size

The polarimetric imaging procedure alternates between imaging via minimization of the objective function (Equation (C5)) and D-term calibration, as described in Section C.2. In the imaging stage, the critical parameters that influence the final reconstruction include the four hyperparameters α_P , α_m , β_HW , and β_TV that set the relative weighting in the objective function between the different data constraints and regularizing terms. In surveying different parameters in the eht-imaging survey, we fix α_m = 1 and vary the other three hyperparameter weights.

In addition to the hyperparameter weights, an additional free parameter in our objective function is the amount of additional systematic noise to add to the data as a budget for non-leakage sources of systematic error (see Papers III, IV). To account for these systematics, we add a term equal to ${f}_{\mathrm{sys}}\times | \tilde{{ \mathcal I }}|$ in quadrature to our baseline thermal noise estimates, where f_sys is an overall multiplicative factor. Finally, we also include as a parameter in our surveys the number of iterations n_iter of alternating between imaging and calibrating the station D-terms. The full list of parameters that we vary in the eht-imaging parameter survey is listed in Table 6, with the fiducial parameters used in reconstructing images of M87 denoted in bold.

To select a fiducial set of parameters that performs best on all six synthetic data tests, we assign each parameter set p two scores on each synthetic data set a; one scoring the fidelity of the polarized image reconstruction s_fid,a ( p ), and one scoring the accuracy of the D-term estimation s_dterm,a ( p ). First, we score the image fidelity by computing the normalized cross-correlation ρ_NX between the reconstructed and ground-truth polarimetric intensity distribution. That is, we use Equation (15) of Paper IV on the images of $\sqrt{{{ \mathcal Q }}^{2}+{{ \mathcal U }}^{2}}$. Then the fidelity score for the parameter set p on the synthetic data set a is simply

$$\begin{eqnarray}&&{s}_{\mathrm{fid},a}({\boldsymbol{p}})={\rho }_{\mathrm{NX}}.\end{eqnarray} \tag{ G1 }$$

We compute the D-term estimation accuracy metric by first calculating the average ℓ₂ distance d_D between the reconstructed D-terms and the ground truth for the data set. For M87, we average this distance over the three stations that we calibrate at this stage: SMT, LMT, and PV. Then we transform this distance to a score between 0 (bad) and 1 (good) by using a sigmoid function;

$$\begin{eqnarray}&&{s}_{{\rm{dterm}},a}({\boldsymbol{p}})=1-{\rm{Erf}}({d}_{{\rm{D}}}/\sqrt{2}{d}_{{\rm{tol}}}),\end{eqnarray} \tag{ G2 }$$

where d_tol is a threshold for the average distance between the ground truth and recovered D-terms beyond which we begin to heavily penalize the reconstruction. We set d_tol = 5%.

Finally, having computed the two scores s_fid,a and s_dterm,a for the parameter set p on the synthetic data set a, we compute a final score s( p ) for the parameter set by multiplying these individual scores together on all six synthetic data sets a:

$$\begin{eqnarray}&&s({\boldsymbol{p}})=\displaystyle \prod _{a}{s}_{\mathrm{fid},a}({\boldsymbol{p}}){s}_{\mathrm{dterm},a}({\boldsymbol{p}}).\end{eqnarray} \tag{ G3 }$$

We then have a final score s for each parameter set incorporating its performance in accurately reconstructing the polarized flux distribution and input D-terms on six synthetic data sets. We take the parameter set p with the highest score as our fiducial parameter set. The fiducial set is indicated in bold in Table 6.

The polsolve algorithm is characterized by several degrees of freedom: the pixel's angular size (smaller values increase the astrometric accuracy of the CLEAN components and, hence, the quality of the deconvolution), the field of view (the effect of this parameter is minimum if a CLEAN masking is applied), the visibility weighting (mainly defined with the Briggs "robustness" parameter, r; Briggs 1995), and the method for division of the Stokes ${ \mathcal I }$ model into sub-components of constant fractional polarization (see Section C.1).

The first step in the polsolve procedure is to generate a first version of the ${ \mathcal I }$ image (using the CASA task clean). Several iterations of phase and amplitude self-calibration (using tasks gaincal and applycal) may be applied to the data, in order to optimize the dynamic range of the ${ \mathcal I }$ model. The self-calibration gains are forced to be equal for the R and L polarizations at all antennas. ¹⁴⁴ Then, the ${ \mathcal I }$ model is split into several sub-components and formatted for its use in polsolve, using the CASA task CCextract. ¹⁴⁵ Finally, polsolve estimates the D-terms of LMT, SMT, and PV, together with the fractional polarizations and EVPAs of all the source sub-components. The estimated D-terms are applied to the data and a final run of clean is peformed, to generate the final version of the full-polarization images.

In the polsolve parameter survey, the ${ \mathcal I }$ division is done in two ways. In one approach, a centered square mask of 50 × 50 μas is created and divided into a regular grid of n × n cells (in this case, cells that do not contain CLEAN components are not used in the fit). Alternatively, a centered circular mask of 40–80 μas diameter is created and divided azimuthally into a regular set of n pieces. The full parameter survey with polsolve consists of an exploration of the following.

• 1.  
  Both types of model sub-divisions (i.e., either regular grid or azimuthal cuts), with n running from 1 to 12.
• 2.  
  Robustness parameter, from r = − 2 to r = 2 (i.e., from uniform to natural weighting) in steps of 0.2.
• 3.  
  Relative weight of the ALMA antenna (which affects the shape of the point-spread function (PSF) of the instrument, especially for values of r far from −2), running from 0.1 to 1.0 in steps of 0.1.
• 4.  
  In all images, a circular CLEAN mask of diameter varying from 40 to 80 μas (in steps of 5 μas) is used. The size of the CLEAN mask is the same as the mask to define the ${ \mathcal I }$ model sub-division.

The pixel size is fixed to 1 μas and the image size covers 256 × 256 μas. We note that for models with extended emission (i.e., Models 1 and 2, see Figure 4) an additional CLEAN mask is added at the southern part of the image.

For each combination of parameters in the survey, we compute the L₁ norm of the differences between true and estimated D-terms, as well as the correlation coefficients, ρ_s (for each Stokes parameter, s) between the CLEAN image reconstructions and true-source structures, properly convolved with the same beam. These quantities can be used to select the best combination of parameters for D-term calibration and image reconstruction (the fiducial imaging parameters for polsolve). Depending on the relative weight that is given to L₁ and the average image correlation, $\rho =({\rho }_{{ \mathcal I }}+{\rho }_{{ \mathcal Q }}+{\rho }_{{ \mathcal U }})/3$, we obtain slightly different fiducial parameters.

In Figure 22, we show an example plot from our polsolve parameter survey for Model 5 (see Figure 4; similar results are found for the rest of models in the survey). The chosen figure of merit is (1 − ρ)L₁, which we show as a function of the robustness parameter r and the number of slices n in the azimuthal sub-division. From Figure 22, we note that the dependence of the figure of merit with n becomes weak for values of n larger than 3–4 and robustness parameters between −1 and −2. This behavior also happens if the regular gridding is used to generate the ${{ \mathcal I }}_{i}$ sub-components. A qualitative explanation of this effect may be that large values of n translate into sub-components of sizes smaller than the synthesized resolution. Therefore, increasing the number of sub-components does not improve the fit, because the small separations between neighboring sub-components correspond to spatial frequencies that are not sampled by the interferometer. Conversely, the fitted fractional polarizations of close by sub-components become highly correlated in the fit, and the L₁ norm of the D-terms saturate around a minimum value.

Zoom In Zoom Out Reset image size

Based on the combined analysis of all the synthetic data sets, we determine the fiducial parameters for polsolve: a robustness parameter of −1 (though −2 produces similar results, especially for models 4 to 6, see Figure 4), a circular slicing with n = 8 sub-components (which produces results similar to 9–10 sub-components, and also similar to those from a regular gridding, n × n, with n = 3–5), relative ALMA weights of 1.0 (which produces similar results as for values between 0.5–1.0), and a circular CLEAN mask of 50 μas diameter.

The standard procedure for D-term calibration using LPCAL is as follows. 1.) Select a calibrator that has either a low fractional polarization or a simple polarization structure, and a wide range of parallactic angle coverage. This is M87 in our case. 2.) Produce a total intensity CLEAN map of the calibrator with e.g., Difmap. 3.) Split the CLEAN model into several sub-models with the task CCEDT in AIPS. LPCAL assumes that each sub-model has a constant fractional polarization and EVPA. Therefore, the more sub-models that we use to divide the Stokes I image, the more degrees of freedom we have for modeling the source linear polarization. 4.) Run LPCAL using the sub-models.

We follow this standard procedure for the D-term calibration. In this work, we consider an additional parameter, the ALMA weight-scaling factor. Down-weighting the ALMA data can be useful when there are significant systematic uncertainties in the ALMA visibilities. In this case, the solutions for other stations can be distorted as the fitting would be dominated by the ALMA baselines due to its high sensitivity and the corresponding smaller formal error bars. In addition to the ALMA weight-scaling factor, we consider the number of CLEAN sub-models as the main parameter that may significantly affect D-term estimation with LPCAL.

We first performed a manual parameter survey using the synthetic data. We reconstructed D-terms with LPCAL by using different numbers of sub-models and ALMA weight-scaling factors, and compared with the ground-truth values. We conclude that using a relatively large number of sub-models (≳10) gives better reconstructions, while the results do not change much when more than 10 sub-models are used. Also, we find that the results are not sensitive to the ALMA weight re-scaling. The strategy and parameters that we adopted could reproduce the ground-truth D-terms in the synthetic data within an accuracy of ∼1% (Figure 5).

Next, we analyzed the low-band M87 data. We used a common approach as for the synthetic data tests, but we allowed several users to independently calibrate and image the real data with different parameter settings to test the robustness of the method. First, each user reconstructed the Stokes ${ \mathcal I }$ image with CLEAN in Difmap. The CLEAN components were divided into several sub-models by the task CCEDT in AIPS. Each user then used LPCAL in AIPS to solve for D-terms for all antennas. This includes the possible residual D-terms of ALMA, APEX, and SMA (as LPCAL cannot set any station D-terms to zero). We let each user test different schemes for CLEAN and use different parameters that were not seen as sensitively impacting results in the synthetic data survey. However, all users split their CLEAN models into a number of sub-models (≳10), in accordance with our synthetic data parameter survey results. Different (u, v) weighting parameters, CLEAN windows, and CLEAN cutoffs were used. Users could choose to downweight ALMA baselines through ParselTongue (Kettenis et al. 2006), a Python interface to AIPS, or average the data in time for LPCAL. We selected fiducial D-terms for the LPCAL pipeline by taking the median of the real and imaginary parts of the surveyed D-terms of each station. This approach allows us to take into account the uncertainties in LPCAL that may be associated with the Stokes ${ \mathcal I }$ image reconstruction and the parameters used for the different tasks. To get the final M87 image for LPCAL, we applied these fiducial D-terms to the data and imaged in. For the high-band results (Appendix I), we obtained D-term solutions and images from a single, best-bet pipeline using representative parameters (15 CLEAN sub-models, ALMA weight = 1.0) from our survey on the low-band data.

Additionally, we investigated the effects of the parameter selection on the real data with GPCAL, using the same parameters as we identify for LPCAL, ¹⁴⁶ to examine the improvements of the fit statistics when changing parameters. Figure 23 shows the distribution of ${\chi }_{\mathrm{red}}^{2}$ for the two parameters from the survey on the M87 data for 2017 April 11. We explored the number of sub-models from five to 20 with an increment of one and the ALMA weight-scaling factor of (0.1, 0.2, 0.4, 0.7, 1.0, 2.0, 5.0, 10.0). Similar to our conclusions from the synthetic data analysis, we found that the fit statistics on the M87 data improve with a larger number of sub-models up to about 10 sub-models. The result is insensitive to changing the ALMA weight-scaling factors. This trend was seen for the M87 data for the other observation days as well. Therefore, we conclude that the parameters that we used for the real data analysis with LPCAL are reasonable.

Zoom In Zoom Out Reset image size

For the DMC image reconstructions, we surveyed two hyperparameters: the pixel separation and the image field of view. Because the DMC method fits for a systematic uncertainty term in addition to the image and calibration parameters, all fits are formally "good" from the perspective of, e.g., a χ² metric. We thus determine an acceptable fit for a particular data set to be the one that minimizes the number of model parameters (i.e., a combination of largest pixel separation and smallest image field of view) while recovering the expected level of systematic uncertainty (i.e., 0% for the synthetic data sets and 2% for the M87 data sets; see Paper III) within some threshold (taken to be the 3σ bounds determined by the posterior distribution).

Associated with Themis reconstructions are two hyperparameters corresponding to the two-dimensional number of control points. These have natural values set by the number of independent modes that may be reconstructed; for the EHT, 5 × 5. Because this resulted in formally acceptable fits, i.e., reduced-χ² near unity, and based on similar considerations from Stokes ${ \mathcal I }$ image reconstructions (see, e.g., Broderick et al. 2020b), no additional exploration was required for M87. For the Model 1 and 2 synthetic data reconstructions, the number of control points were incrementally increased until acceptable fits were obtained.

To sample the effects of uncertainties in the D-terms on the reconstructed images in the posterior exploration methods (DMC and Themis), we can simply draw our 1000 image sample randomly directly from the full posterior. For the three non-MCMC imaging pipelines (eht-imaging, LPCAL, polsolve), our method is to use a simple Monte Carlo approach, similar to the analysis of Martí-Vidal et al. (2012). For each method, we draw 1000 random sets of D-terms from normal distributions with means and (co)variances determined by considering the fiducial results from Table 5 across the four observing days. We assume for this test that the uncertainties in the D-terms are uncorrelated from station to station and between LCP and RCP. This represents a worst-case scenario test, as correlations between the D-terms would reduce the volume of the D-term parameter space surveyed by each method. The full sample of 1000 D-term sets sampled for SMT, LMT, and PV on 2017 April 11 for the eht-imaging, polsolve, and LPCAL pipelines are shown in Figure 24. In addition to the distributions shown in Figure 24, we also include D-terms on ALMA, APEX, SMA, and JCMT drawn from circular complex Gaussians with 1% standard deviation; these represent residual uncertainties left over from the zero-baseline D-term calibration of these stations.

Zoom In Zoom Out Reset image size

After drawing a given set of random D-terms, we apply this calibration solution to the data and reconstruct a polarized image using the same fiducial imaging procedure described for each method in Appendix G. Our imaging scripts in this stage do not involve any simultaneous leakage calibration, but only reconstruct the Stokes ${ \mathcal Q }$ and ${ \mathcal U }$ from the visibilities with the assumed calibration solution already applied. ¹⁴⁷ That is, in this procedure we draw a set of random D-terms from distributions reflecting our uncertainty in the recovered D-terms from the earlier stage of calibration and imaging, and then we reconstruct an image assuming that this D-term calibration is perfect with no need for further leakage calibration.

In Figures 25 and 26, we show histograms over each imaging method's sample of 1000 images with different D-term calibration solutions on all four observing days of the four image-integrated quantities used in Section 5.3. Figure 25 shows histograms of the image net polarization fraction ∣m∣_net (Equation (12), plotted in red) and the intensity-weighted average polarization fraction at 20 μas resolution 〈∣m∣〉 (Equation (13), plotted in green). Figure 26 shows the amplitude ∣β₂∣ (plotted in brown) and phase ∠β₂ (plotted in purple) of the β₂ coefficient of the azimuthal decomposition defined in Equation (14). Because the observations on 2017 April 5 and 11 have the highest-quality (u, v) coverage and bracket the observed time evolution of the source, we choose to define acceptable ranges for these parameters (the shaded bars in Figures 25, 26, taken from Table 2) using only these two days. In particular, the poor quality of the 2017 April 10 (u, v) coverage leads to broader distributions of the four key quantities with large systematic uncertainties between imaging methods (third columns of Figures 25, 26). The distributions on 2017 April 5 and 11 are summarized with mean and 1σ error bars in the main text Figure 9, and are discussed in Section 5.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, Table 7 presents ranges of the image-integrated Stokes parameters ${ \mathcal I },{ \mathcal Q },{ \mathcal U }$ derived from the surveys over different D-term calibration solutions from each day of observations. The ranges in Table 7 were calculated by taking the minimum mean − 1σ and maximum mean + 1σ point from the five individual method surveys on each day.

Table 7. Ranges of Image-integrated Stokes ${ \mathcal I },{ \mathcal Q },{ \mathcal U }$ Obtained from Each Method's D-term Calibration Survey on the Low-band Data

Stokes (Jy)	Min	Max	Min	Max
	2017 April 5		2017 April 6
${ \mathcal I }$	0.419	0.512	0.376	0.508
${ \mathcal Q }$	−0.0136	−0.0002	−0.0056	0.0075
${ \mathcal U }$	−0.0157	−0.0055	−0.0167	−0.0063
	2017 April 10		2017 April 11
${ \mathcal I }$	0.381	0.545	0.407	0.565
${ \mathcal Q }$	−0.0067	0.0158	0.0030	0.0180
${ \mathcal U }$	−0.0074	0.0115	−0.0113	0.0007

Download table as:  ASCIITypeset image