First M87 Event Horizon Telescope Results. IX. Detection of Near-horizon Circular Polarization

Event Horizon Telescope (EHT) observations have revealed a bright ring of emission around the supermassive black hole at the center of the M87 galaxy. EHT images in linear polarization have further identified a coherent spiral pattern around the black hole, produced from ordered magnetic fields threading the emitting plasma. Here we present the first analysis of circular polarization using EHT data, acquired in 2017, which can potentially provide additional insights into the magnetic fields and plasma composition near the black hole. Interferometric closure quantities provide convincing evidence for the presence of circularly polarized emission on event-horizon scales. We produce images of the circular polarization using both traditional and newly developed methods. All methods find a moderate level of resolved circular polarization across the image (〈∣v∣〉 < 3.7%), consistent with the low image-integrated circular polarization fraction measured by the Atacama Large Millimeter/submillimeter Array (∣v int∣ < 1%). Despite this broad agreement, the methods show substantial variation in the morphology of the circularly polarized emission, indicating that our conclusions are strongly dependent on the imaging assumptions because of the limited baseline coverage, uncertain telescope gain calibration, and weakly polarized signal. We include this upper limit in an updated comparison to general relativistic magnetohydrodynamic simulation models. This analysis reinforces the previously reported preference for magnetically arrested accretion flow models. We find that most simulations naturally produce a low level of circular polarization consistent with our upper limit and that Faraday conversion is likely the dominant production mechanism for circular polarization at 230 GHz in M87*.


INTRODUCTION
The Event Horizon Telescope (EHT) has produced the first images of the event-horizon-scale millimeter emission around the supermassive black hole in the core of the massive elliptical galaxy M87 at the center of the Virgo Cluster.Using very-long-baseline interferometry (VLBI) at 230 GHz, these initial EHT observations from 2017 recovered a ring-like structure with a diameter similar to the predicted "black hole shadow" of a M BH ≈ 6.5 × 10 9 M ⊙ black hole at the distance of M87 * (D ≈ 16.8 Mpc).The resolved total intensity images of the ring were consistent with models of synchrotron emission from ultra-hot magnetized plasma near the event horizon (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f, hereafter Papers I-VI).
The EHT observes in full polarization, recording simultaneous data from orthogonally polarized feeds at each antenna.Images of the near-horizon linearly polarized radiation were published and analyzed in Event Horizon Telescope Collaboration et al. (2021a, hereafter Paper VII).These linear polarimetric images provide essential new information about the magnetic field struc-ture near the event horizon of M87's supermassive black hole, indicating that the near-horizon magnetic fields are ordered and dynamically important (Event Horizon Telescope Collaboration et al. 2021b, hereafter Paper VIII).
In this Letter we report on the search for resolved circularly polarized radiation (CP) on event-horizon scales in M87 * from EHT observations in 2017.The circularly polarized signal from synchrotron radiation near the black hole should contain unique information about the magnetic field and the nature of the radiating particles that cannot be inferred from total intensity or linear polarization alone.These include the possibility of directly measuring the strength of the magnetic field and determining whether the observed radiation is mainly from an electron-positron or an electron-ion plasma (e.g., Wardle et al. 1998).However, the circularly polarized signal is expected to be weaker than the linear polarization (Jones & O'Dell 1977;Jones 1988), requiring high sensitivity and accurate calibration of each antenna to be detected.
Previous radio and millimeter-wavelength observations of CP in M87 * are quite limited, whether at low or high angular resolution.Homan & Lister (2006), using the Very Long Baseline Array (VLBA) at 15 GHz, measured fractional circular polarization in the core of −0.49 ± 0.10 %.Interestingly, they found no linear polarization (LP < 0.07 %), reversing the Stokes parameter hierarchy described above.M87 is the only source in their list to show such behavior.On the other hand, using the Very Large Array (VLA) at 8.4 GHz, Bower et al. (2002) detected weak fractional linear polarization (LP = 1.74 ± 0.06 %) but no circular polarization (|CP| < 0.1 %).Single dish measurements of M87 with the Effelsberg 100m radio telescope at the same range of frequencies showed a similar trend with LP ∼ 1.5% and |CP| < 0.2% (Myserlis et al. 2018), while at 86 GHz the POLAMI monitoring program observed a fairly stable LP ∼ 5% and CP∼ −1.5% over a period of 12 years with the IRAM 30m telescope (Thum et al. 2018).Simultaneous observations during the 2017 EHT campaign with the Atacama Large Millimeter/submillimeter Array (ALMA) did not result in a significant detection of the unresolved fractional CP at 221 GHz (≈ −0.3 ± 0.6 %; Goddi et al. 2021).
Beyond M87 * , accurate measurements of circular polarization are generically difficult to obtain in VLBI.Homan & Lister (2006) detected CP at the level of 3σ or better in 17 sources out of their sample of 133 Active Galactic Nuclei (AGN) using a "gain transfer" technique (Homan & Wardle 1999), in which all sources are used for polarization calibration.These detections at 15 GHz all had fractional polarizations between 0.25 % and 0.70 %.Gabuzda et al. (2008) also detected circular polarization in eight AGN jets with 15 GHz VLBI measurements, and associated the observed CP signs with CP production by Faraday conversion in helical jet magnetic fields.At higher frequencies up to 43 GHz, Vitrishchak et al. (2008) found circular polarization fractions up to ∼ 1% in a sample of AGN cores resolved with VLBI.
This paper presents the details of the 2017 EHT observations and data calibration for circular polarization, procedures and results for circular polarimetric imaging, and their theoretical interpretation in constraining parameters in a library of simulation models.In Section 2, we summarize the EHT 2017 observations, describe evidence for the detection of circular polarization in M87 * , and describe our a priori calibration procedure.In Section 3, we describe our methods of circular polarimetric image reconstruction and gain calibration from EHT data.In Section 4, we present image reconstructions from all methods across all EHT 2017 datasets.We derive an upper limit on the average, resolved circular polarization fraction in M87 * at 230 GHz.
In Section 5, we examine circular polarization in a library of General Relativistic Magnetohydrodynamic (GRMHD) simulations of M87 * .We add upper limits to the circular polarization on event-horizon scales from our EHT observations to the list of constraints applied to theoretical models in Paper VIII.We discuss the effect of these new constraints on our preferred models for M87 * 's accretion flow, and we investigate the physical origin of circular polarization in passing GRMHD models.We summarize the work in Section 6.

EHT OBSERVATIONS OF M87 *
2.1.Conventions A radio interferometer, such as the EHT, samples the Fourier plane (visibility domain) of the brightness distribution (image) on the sky (e.g., Thompson et al. 2017).Each image domain Stokes parameter (I, Q, U, V) has a corresponding visibility, which we denote with a tilde (e.g., Ṽ).At most stations, the EHT data are recorded in two orthogonal circular polarizations, right and left (R, L).Interferometric visibilities are computed through the complex correlation between each pair of sites (j, k) and polarization (i.e., R j R * k , R j L * k , L j R * k , and L j L * k ).These measurements can be transformed to the Stokes representation through linear algebraic relations.In particular, the visibility domain circular polarization Ṽ on the j − k baseline is given by In practice, polarimetric measurements are corrupted through a multitude of systematic effects, which can be conveniently represented in a matrix formalism (e.g., Hamaker et al. 1996;Smirnov 2011): In this expression, ρ jk is a 2×2 matrix of the true visibilities on the j − k baseline, ρ ′ jk gives the measured visibilities, and J is a complex time-and site-dependent 2×2 matrix that describes the aggregate systematic effects.The latter can be further decomposed into a product of three terms, J = GDΦ, where G is a diagonal matrix that describes the time-dependent "gains" of the two feeds, D is a matrix with diagonal entries that are unity and off-diagonal entries that describe "leakage" between the feeds, and the matrix Φ describes the overall rotation of the feeds.For studies of circular polarization with circularly polarized feeds, the gain matrix is the most important source of contamination, predominantly through the gain ratio G R/L , In particular, while the gains G L and G R have rapid variations, especially in their phase, the gain ratio G R/L can be stable over timescales of multiple days.For addi-

Observations and data reduction
The EHT observed M87 * on 5, 6, 10, and 11 April 2017 with a VLBI array of 7 telescopes located at 5 geographical sites: 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 IRAM 30 m telescope (PV) on Pico Veleta in Spain; the Submillimeter Telescope (SMT) on Mt.Graham in Arizona (USA); and the Submillimeter Array (SMA) and the James Clerk Maxwell Telescope (JCMT) on Maunakea in Hawai ' i.The South Pole Telescope (SPT) also participated in the EHT observations but cannot see M87 * .
Each telescope recorded two frequency bands, each 2 GHz wide, centered at 229.1 GHz (high band, HI) and at 227.1 GHz (low band, LO).Most sites recorded right and left circular polarization simultaneously, except for the JCMT which recorded a single circular polarization each night, and ALMA, which recorded orthogonal linear polarizations that were subsequently converted to a circular basis (Martí-Vidal et al. 2016;Matthews et al. 2018;Goddi et al. 2019).Paper II provides a detailed description of the EHT array and observations.
The EHT data sets were first correlated and then calibrated using two independent pipelines, EHT-HOPS (Blackburn et al. 2019) and CASA rPICARD (Janssen et al. 2019), for the stabilization of the source signal in time and frequency (e.g., Janssen et al. 2022).The data presented in this paper correspond to the EHT-HOPS pipeline, following a verification of the inter-pipeline consistency (Paper III), although in Appendix A we also give a brief summary of the gain calibration in the CASA rPICARD pipeline, and we compare imaging results from EHT-HOPS and rPICARD pipeline data in Section 4. The amplitude flux density calibration was performed with custom-built EHT-HOPS post-processing scripts, based on the metadata provided by the participating telescopes.A more extensive description of other aspects of the EHT M87 * data set calibration was presented in Paper III.The data sets analyzed in this paper are identical as the ones used in Paper VII and Paper VIII, following the same polarimetric calibration procedures.The full-Stokes calibrated VLBI data from the 2017 EHT observations of M87 are publicly available through the EHT data portal1 under the code 2023-D01-01(Event Horizon Telescope Collaboration 2023).

R/L gain calibration
We calibrated the gain ratio G R/L for each site using multi-source and multi-day polynomial fits to the differences between RR * and LL * visibilities.Since ALMA observes in a linear-polarization basis, it provides the absolute electric vector position angle (EVPA) informa-tion and it is used as the reference station; its G R/L is fixed to unity (Martí-Vidal et al. 2016;Goddi et al. 2019). 2  For the other sites, the R/L gains were modeled separately for each band using visibilities on baselines to ALMA, using custom-built Python scripts.The use of multiple sources (10 sources were observed with ALMA during the EHT 2017 campaign) in the fitting procedure helps to distinguish between instrumental effects, which are largely source-independent, and contributions from circular polarization, which are source-dependent (Steel et al. 2019).
For most sites, G R/L could be successfully modeled using a constant value across the full observing campaign, separately for each frequency band.The sites that required a time-dependent G R/L were APEX, which showed a linear phase drift between RR * and LL * visibilities, and the SMA, which showed irregular phase variation that we modeled using third order polynomials specific to each day and frequency band.We calibrated the amplitudes |G R/L | assuming a constant value for each site.
In Appendix A, we provide more details on this strategy for relative gain calibration, as well as an example plot detailing how gains were estimated for the LMT (Figure 14).We also compare the above strategy for relative gain calibration used for the main results in this paper with the alternative strategy employed by the CASA-based rPICARD pipeline for EHT data reduction, which calibrates relative gains assuming the intrinsic V = 0.

Evidence for circular polarization
Under a perfect calibration, we could directly interpret the signal to noise (S/N) ratio of Ṽ, calculated from observables with Equation 1, to identify robust detections of CP.In Figure 1 we show the S/N of Ṽ as a function of a projected baseline length.Unfortunately, as circular polarization is encoded in the difference between RR * and LL * visibilities, residual errors in the calibration of G R/L will create spurious Ṽ signatures.As a consequence, the S/N will be inflated, and hence values shown in Figure 1 can only be treated as upper limits.This is further illustrated in Appendix A, where we discuss a more aggressive calibration strategy for G R/L , which reduces the S/N of Ṽ significantly.
There are, however, quantities that are robust against the effects of complex antenna gains.An example of such quantities are the closure phases (e.g., Thompson et al. 2017;Blackburn et al. 2020), defined for a triangle 2 A clockwise 45 degree rotation was also applied to all observations across the whole array (i.e., 90 degrees added to the G R/L phases of all stations), in order to account for the orientation of the ALMA Band 6 cartridges with respect to the antenna mounts (Paper VII). of antennas (A, B, and C) as where ϕ ij is the visibility phase measured by baseline i − j.One way to minimize the effects from possible antenna mis-calibrations and infer the presence of resolved source-intrinsic circular polarization in our observations is to compute the difference of closure phases between the RR * and LL * visibilities; by construction, the difference of these closures will not be affected by G R/L terms at any antenna. 3n Figure 2 (top panels), we show the closure phases, ψ, of the antenna triangle formed by ALMA, SMT, and PV, computed for the whole campaign (left column, empty markers for low band; right column, filled markers for high band).Closure phases for the RR * and LL * visibilities are shown in red and blue, respectively.It can be seen that there is an offset in the closure phases between the two polarization channels, which generates a non-zero closure difference (indicated in green in the bottom panels).This difference in the RR and LL closure phases is present in all epochs and in both bands.
For this antenna triangle, the average offset in closure phases between RR * and LL * (combining all epochs and assuming a constant residual value) is 6.7±1.3 deg in the low band and 7.9 ± 1.6 deg in the high band, indicated with green bands in the bottom panels.The offset is consistent despite each band being calibrated independently.Moreover, the measured offset is difficult to explain with the conservative systematic non-closing error budget discussed in Section 8.4 of Paper III, and hence it implies a tentative detection of a weak Stokes Ṽ at the level of S/N ∼ 5.
While this measurement implies the presence of a fractional CP reaching ∼ 3% somewhere in the visibility domain, the measurement can not be directly translated into a quantitative image domain constraint.The ALMA-SMT-PV triangle shown is the one that produces closure-phase differences with the most clear deviation from zero.In Appendix B, we show the results for all other triangles including ALMA.None of these other triangles shows an unambiguous detection of a nonzero closure phase difference like that seen on the ALMA-SMT-PV triangle, suggesting that SMT-PV is the baseline most sensitive to the CP signatures in M87 * .
Even though the RR * and LL * closure phases are robust against antenna gains, they may be affected by instrumental polarimetric leakage (antenna D-terms, the D matrix term in Equation 3).However, the effect of D-term uncertainties in the parallel-hand visibilities is much smaller than in the cross-hands (e.g., Smirnov 2011), which implies that Stokes V is much less affected by instrumental polarization than Stokes Q and U. To verify a negligible impact of the polarimetric leakage on the closure phase signal, we have compared the closurephase values between the data with and without D-term calibration.For all triangles related to ALMA, the effect of the D-terms on the RR * − LL * closure differences is always less than the standard deviation of the thermal noise on the closure phase difference.For the triangle shown in Figure 2, the maximum effect of the D-terms is only 0.42 σ.
Hence, we can conclude that the closure phase differences indicating the presence of circular polarization on EHT baselines, presented in Figure 2, are robust against both antenna gains and polarimetric leakage.In Appendix C, we discuss evidence for polarization in "closure trace" products (Broderick & Pesce 2020), quantities that are insensitive to all station-based systematic factors, including D-terms.Note-Each method produces images of M87 * in all four Stokes parameters.The circular polarization signal is weak and strongly depends on calibration assumptions, so we summarize the primary differences in calibration assumptions among the methods.For additional details (e.g., priors adopted for fitted values), see the more detailed descriptions in Appendix D.

POLARIMETRIC IMAGING METHODS
As discussed in Section 2, the circular polarization signal in M87 * is both weak and sensitive to calibration errors.Consequently, the inferred circular polarization structure can be sensitive to assumptions made about the residual calibration errors as well as choices made in producing images from the measured visibilities.To assess the potential impact of these effects, we produce fully polarized images of M87 * using five methods.We summarize these methods and their assumptions in Table 1.In Appendix D, we present more detail on each method's assumptions and procedures used for circular polarization image reconstructions.
The five methods we use can be divided into three general categories.The first category uses adaptations of the standard CLEAN imaging algorithm together with iterative "self-calibration" to solve for calibration errors; we use the software DIFMAP (Subsection D.1) and polsolve (Subsection D.2) as two implementations of this approach.The second category uses a Bayesian forward modeling approach, jointly solving for both a polarized image raster and for residual calibration errors and providing estimates for the posterior distributions of each; we use the software Themis and DMC (Subsection D.3).The third category also uses a Bayesian forward modeling approach, but it uses a simple geometric model for the sky image and only fits VLBI "closure" quantities to constrain the circular polarization structure; we use the software eht-imaging (Subsection D.4).
In addition to the differences inherent in each approach as a result of the underlying method (CLEAN vs raster fitting vs geometric modeling, Bayesian exploration vs fitting a single image), our methods face additional choices on how they deal with polarimetric leakage (D-terms) and residual gain ratios (G R/L ).Constraints on the D-terms were derived and discussed extensively in Paper VII; these constraints have been confirmed with analysis of EHT observations of AGN in Issaoun et al. (2022) and Jorstad et al. (2023).As a result, some methods (DIFMAP, and m-ring modeling) chose to directly apply the Paper VII D-term results to the data and not treat polarimetric leakage further; in contrast, DMC and Themis fully explored uncertainties in the D-terms under flat priors as part of their Bayesian posterior exploration.polsolve performed D-term self-calibration as a part of its imaging procedure using priors motivated by the Paper VII results.
Most significantly for circular polarization, each method had freedom to chose how to approach residual uncertainties in the gain ratios G R/L , in both amplitude and phase.Three imaging methodspolsolve, DIFMAP and DMC-solved for separate G R and G L terms as part of their self-calibration (polsolve and DIFMAP) or Bayesian forward modeling (DMC) approach.In contrast, Themis did not solve for separate right-and left-circular gains, but only solved for an overall gain term G = G R = G L , absorbing any uncertainty in the residual gain ratios into the recovered image structure.The geometric m-ring modeling method did not solve for any gain terms, as it directly fit gain-insensitive closure quantities in right-and left-circular polarizations independently.These choices in the treatment of the gain ratios can have a large impact on the recovered circular polarization structure and its uncertainty; for instance, while the methods are otherwise similar in their Bayesian approach and treatment of the D-terms, the weak priors on the gain ratios in DMC result in larger error bars in the recovered Stokes V structure as compared to Themis, which locks all gain ratios to unity.
Before applying our imaging methods to EHT observations of M87 * , we first tested each method on synthetic data taken from GRMHD simulation images of M87 * on EHT2017 baselines.We present these tests in Appendix E. Our results on synthetic data suggest that our imaging methods are generally not capable of unambiguously determining the horizon-scale structure of circular polarization in M87 * , given the low fractional polarization in the source (and GRMHD models), as well as the poor u − v coverage of the EHT in 2017.Indeed, when turning to real data in Section 4, we also find that our methods are unable to arrive at a single consistent image of the circular polarization in M87 * .However, we are able to use the images presented in Section 4 to derive an upper limit for the circular polarization fraction on horizon-scales.

Imaging results on individual days and bands
Each imaging method introduced in Section 3 was used to produce Stokes I, Q, U, and V images from the eight individual 2017 EHT data sets of M87 * , corresponding to the four observation days in both low and high band.Imaging methods were free to use data that were pre-calibrated for the zero-baseline D-terms derived in Paper VII or not, and to solve for residual G R/L gain errors or assume they are fixed to unity in selfcalibration.The choices made by each imaging method are summarized in Table 1.Of the images presented here, the DMC and Themis posteriors for the M87 * data are identical to those already presented for linear polarization in Paper VII.The polsolve results were generated with a similar procedure as the linear polarization results in Paper VII, but with an additional series of imaging and self-calibration for recovering CP and the G R/L offsets.The DIFMAP and m-ring modeling results are new to this work.
Figure 3 shows the April 11 low band reconstructions in both total intensity and linear polarization in the top row (in the style of Paper VII) and in an "ellipse-plot" representation of the total linear plus circular polarization in the bottom row.The "ellipse-plots" illustrate the degree of linear polarization relative to circular polarization by the eccentricity of small ellipses plotted across the image.As in Paper VII, all imaging methods recover consistent images of total intensity and linear polarization.The ≈ 40 µas diameter ring structure is recovered in all methods, as is the ≈ 15 % peak linear polarization fraction and predominantly azimuthal EVPA pattern in the south-west part of the image.Minor differences between images in the fractional linear polarization appear at the edges of the ring and in the m-ring pattern, which is limited by a small number of degrees of freedom in the m = 3 mode fit in linear polarization.The ellipse plots in the bottom row show that in all cases the circularly polarized brightness recovered is a small fraction (≲ 20 %) of the total polarized brightness, which is indicated by the large axis ratio/eccentricity of each ellipse.The colors of each ellipse show the sign of circular polarization recovered at each point.For April 11 low band there is a consistent negative sign of V at the total intensity maximum in the south-west, but across the rest of the image there is little consistency between methods in the sign of V.
In Figure 4 and Figure 5 we focus on the results of our EHT 2017 M87 * circular polarization imaging and modeling results by showing circular polarization images within contours indicating the I brightness.In the color map chosen, red corresponds to a positive sign of CP while blue corresponds to negative CP.As in the synthetic data tests, different methods show consistent structure in total intensity, but significantly different structure in circular polarization.The imaging results are most consistent for low-band data, where on most days most methods recover negative circular polarization at and near the total intensity maximum.In general, however, the Stokes V structures across the image are not consistent from method-to-method.Furthermore, the V images are not consistent between bands when imaged with the same method.An exception is mring modeling, which consistently indicates an approximately North-South asymmetry, with more negative circular polarization in the South.These results must be interpreted carefully, however, as strong assumptions on the source structure are imposed in the choice of a simple m = 1 model for fitting the circular polarization.The m-ring results are discussed in more detail in R23.
To test whether the observed inconsistency in the recovered Stokes V images is mitigated by a different calibration strategy, we produced images with a subset of methods using data on April 5 and April 11 reduced with the CASA-based rPICARD pipeline (Janssen et al. 2019).We made use of a new rPICARD calibration mode, where instrumental phase and delay offsets are solved initially to align the the RR * , RL * , LR * , and LL * phases of the high and low bands.Subsequently, all fringe-fitting and phase calibration solutions are obtained from the combined data of the four correlation products and the two frequency bands.Additionally, unlike the data reduced with the EHT-HOPS pipeline used for Figure 4 and Figure 5, where a multi-day and multi-source fit is employed, the CASA data utilized here is calibrated with G R/L assuming V = 0.The station-based amplitude gains are solved every few minutes and do not affect closures.We note that the R-L gain calibration strategies employed for the HOPS/CASA data will likely underfit/overfit instrumental gain offsets.The images produced form the CASA data still show inconsistency across imaging methods and observing days, and do not change our main conclusions based on the EHT-HOPS images, which we adopt as fiducial for the rest of the paper.We show full results of this test in Appendix F.
The general inconsistency among methods and across bands is in sharp contrast to the total intensity images (bottom two rows).Images of circular polarization on these consecutive days are expected to be nearly identical, as is seen in total intensity and linear polarization.The top/bottom row in each pair shows results from imaging the high/low-band data.In each panel, total intensity is indicated in the colored linear-scale contours, and the Stokes V brightness is indicated in the diverging colormap, with red/blue indicating a positive/negative sign.The colorbar ranges are fixed in both plots (and in Figure 5).For posterior exploration methods (DMC, Themis, m-ring fitting), the posterior-average image is shown.All images are blurred with the same 20 µas FWHM Gaussian, shown with the black inset circle in the upper-left panels.
presented in Paper IV and linear polarization images from Paper VII.This inconsistency is a result of the severe difficulties in recovering resolved circular polarization from sparse 2017 EHT observations with low S/N and G R/L calibration uncertainties.

Combining days and bands
Given the lack of consistency in the Stokes V reconstructions presented in Figure 4 and Figure 5 for individual EHT data sets, it is natural to wonder if by averaging data across frequency bands and in time we may increase S/N enough to more confidently recover circular polarization structure.A subset of our imaging methods tested this hypothesis.Figure 6 presents imaging results from DMC, Themis, and DIFMAP on datasets combining high and low band EHT observations of M87 * on two pairs of days: April 5 and 6 2017 (top row) and April 10 and 11 (bottom row).The source structure in M87 * evolved slightly over the week of observations in 2017, but was stable on each pair of days combined here (Paper IV,Paper VII).As in Figure 4 and Figure 5, we find no consistency in the reconstructions from the band-and day-averaged data.While the S/N is improved by a factor of two by averaging four data sets in each combined reconstruction, the low intrinsic circular polarization and residual G R/L uncertainty are still too large and result in inconsistent image reconstructions from the combined data.We focus on results derived from the individual day/band images in the rest of this work.

Upper limit on the resolved circular polarization fraction
We quantify the average circular polarization fraction in resolved EHT images by the image-averaged fractional circular polarization magnitude |V/I|, weighted by the Stokes I intensity: where the integral is over the whole area of the image A. The definition of ⟨|v|⟩ in Equation 5 is in close analogy with the average resolved linear polarization fraction ⟨|m|⟩ used in Paper VII and Paper VIII.We contrast this average polarization fraction on EHT scales with the unresolved circular polarization fraction By definition, ⟨|v|⟩ depends on the image resolution, while v net is invariant to convolution with a blurring kernel.The definition of ⟨|v|⟩ naturally down-weights contributions from regions where the total intensity image I is dim and noisy.However, ⟨|v|⟩ is by definition constrained to be positive and thus cannot indicate the predominant sign of circular polarization in an image.For weakly polarized/noisy images, ⟨|v|⟩ will also be biased to non-zero values (see Appendix G for a discussion).
We chose to use ⟨|v|⟩ instead of alternative metrics like V/I at the peak brightness location because our image reconstructions and GRMHD synthetic images often do not show a single peak in circular polarized intensity.Furthermore, image-integrated metrics are preferable for computing sensible posterior means and error bars from the results of Bayesian methods like DMC, Themis, and m-ring model fitting.
Figure 7 shows the image-integrated and resolved circular polarization fractions for all reconstruction methods, days, and bands of EHT 2017 M87 * data.For the posterior exploration methods (DMC, Themis, and mring modeling), we plot the posterior median of each quantity and error bars corresponding to 2.27th and 97.7th percentile of the posterior distribution (corresponding to 2σ error bars if the posterior were Gaussian).For DIFMAP and polsolve, we plot the single value corresponding to the image results in Figure 4 and Figure 5; we derive approximate 2σ error bars for these methods based on the measured off-source residuals in V and I using standard Gaussian error propagation.DMC has particularly large error bars on fractional circular polarization because of its permissive priors on the relative gain ratios G R/L .The values of the integrated circular polarization fraction v net in Figure 7 recovered by each imaging method are typically within the upper limits reported by Goddi et al. (2021) from ALMA observations, though some methods produce anomalously larger integrated polarization fractions on certain data sets (e.g.polsolve on both April 11 data sets, DIFMAP and Themis on April 6 low band data).Our methods do not recover a consistent sign of the integrated circular polarization fraction v net .
The values of the resolved circular polarization fraction at 20 µas resolution ⟨|v|⟩ for most reconstructed images in the right panel of Figure 7 are typically less than 4 %.The synthetic data results in Figure 20 indicate that for 2017 EHT coverage and sensitivity most methods will always produce ⟨|v|⟩ ≳ 1 %, even when the actual value is lower, as a consequence of uncertainties in the reconstruction in the low S/N regime, uncertainty in the residual gain ratios G R/L , and an upward bias on the quantity ⟨|v|⟩ (Appendix G).As a result, and because of the lack of agreement among methods in recovering a consistent source structure in V in Figure 4 and Figure 5, we use these results only to obtain a conservative upper limit for the resolved circular polarization fraction.
We estimate a combined upper limit on ⟨|v|⟩ in M87 * from each method's eight measurements across the four days and two observing bands.We average the measurements of ⟨|v|⟩ within each method across these eight data sets using standard inverse-variance weighting, then compute the 99th percentile of the resulting distribution (assumed to be Gaussian).These 99 % upper limits for each method are reported in Table 2.
We adopt the DMC results as the most conservative upper limit on ⟨|v|⟩.The upper limit computed from polsolve is slightly higher than that from DMC (3.8 % vs 3.7 %).Nonetheless, we adopt the DMC value as our fiducial upper limit in this work.DMC performs full Bayesian posterior exploration over the image intensity values and station gains.In contrast, polsolve computes only a single image with error bars computed using idealized Gaussian error propagation.None of our interpretation in Section 5 is affected by choosing 3.8 % vs 3.7 % as our fiducial upper limit.For these reasons, we adopt the DMC value ⟨|v|⟩ < 3.7 % going forward as our conservative upper limit on the average resolved horizon-scale circular polarization fraction in M87 * at 230 GHz.

THEORETICAL INTERPRETATION
On event horizon scales, circularly polarized images of hot accretion flows can encode valuable information about the plasma, including its magnetic field geometry and composition.While our imaging methods are unable to unambiguously determine the horizon-scale  2021) (used to constrain GRMHD models in Paper VIII) are indicated in the gray shaded region.In the right panel, the upper limit on ⟨|v|⟩ derived in this work (⟨|v|⟩ < 3.7 %, Table 2) is indicated by the dashed horizontal line.Note-We derive each upper limit using inverse-variance weighted averaging of the results in the right panel of Figure 7.We adopt the conservative upper limit, ⟨|v|⟩ < 3.7 % from DMC, as our fiducial value.Note that in computing the upper limit from polsolve, we exclude the anomalously highly polarized April 10 low band result, as it suffers from an anomalously large overall offset in the recovered G R/L .
structure of circular polarization, in Section 4 we establish an upper limit on the magnitude of the fractional circular polarization magnitude on scales of the EHT beam: ⟨|v|⟩ < 3.7%.This upper limit, combined with the existing limit on the unresolved, sourceintegrated circular polarization fraction from ALMA (|v net | < 0.8%), can be used to constrain models of the emitting plasma and accretion flow around M87 * , even without additional information on the structure of the circularly polarized emission.
Recent work has revealed that intrinsically circularly polarized synchrotron emission, Faraday conversion and rotation, and twisted field geometries are all important in generating the CP image that we observe in millimeter wavelengths (Tsunetoe et al. 2020;Mościbrodzka et al. 2021;Ricarte et al. 2021;Tsunetoe et al. 2021Tsunetoe et al. , 2022)).For the plasma parameters of interest in M87 * , millimeter synchrotron emission is intrinsically circularly polarized only at the ∼1 % level.However, Faraday conversion also plays a role in exchanging linear and circular polarization states, and in fact dominates the Stokes V production in many models.Stokes V generated by Faraday conversion is understood to be sensitive to the magnetic field geometry, which has been utilized to infer helical field structure in jets (e.g., Wardle & Homan 2003;Gabuzda et al. 2008).On event horizon scales, the connection between Stokes V and the magnetic field geometry is complicated by the geometries probed by geodesics that take non-trivial paths through the space-time, leading in particular to Stokes V sign flips of successive sub-images in some models (Mościbrodzka et al. 2021;Ricarte et al. 2021).
As with the previous papers in this series, our theoretical interpretation centers on ray-traced GRMHD simulations, which self-consistently generate the plasma that performs emission, absorption, and Faraday effects.Since images of total intensity and linear polarization were studied in detail in Paper V and Paper VIII respectively, here we focus mainly on the astrophysics governing the generation of circular polarization.In all cases, τρV > τρQ, indicating that with a constant field orientation in the emission region, circular polarization will dominate over linear polarization in these models.(Right) The average fractional circular polarization between 5 rg and 10 rg in one-zone models with a rotating magnetic field direction along the line of sight, as a function of the angular rotation frequency ω.We show three different models: a model with low Faraday conversion depth (blue), a model with median conversion depth (orange), and a model with high conversion depth (green).Dashed lines show corresponding results for one-zone models with no intrinsic emission of circular polarization, jV = 0.

One-zone models
In Paper V and Paper VIII, we used one-zone isothermal sphere models to derive order-of-magnitude estimates on the parameters of the synchrotron emitting plasma in M87*.In particular, applying constraints on the linear polarization, Paper VIII found that one-zone models for M87 * imply the dimensionless electron temperature Θ e = k B T e /m e c 2 lies in a mildly relativistic regime, 2 < Θ e < 20, the magnetic field is in the range 1 ≲ |B| ≲ 20 G, and the number density lies in the range 104 ≲ n e ≲ 10 7 cm −3 .While one-zone models neglect the critical effects of fluid velocity, gravitational redshift, and a non-uniform emitting region, they favor plasma parameters that are in general agreement with those found in favored GRMHD simulations (Paper VIII).
Here, we explore some implications for circular polarization in the one-zone model space selected in Paper VIII.We use the same model of a uniform R = 5r g isothermal sphere 4 with a magnetic field at a pitch angle θ B = π/3 to the line-of-sight.We determine polarized emissivities j I , j Q , j V , absorption coefficients α I , α Q , α U and Faraday rotation/conversion coefficients ρ V , ρ Q using the fitting functions for relativistic thermal electron distributions in Dexter (2016).
First, we consider the importance of Faraday conversion in producing the observed circular polarization in these models.The left panel of Figure 8 shows the distribution of the Faraday conversion optical depth zone models that pass the Paper VIII constraints.In all cases τ ρ Q > 1, indicating that Faraday conversion dominates intrinsic emission in producing circular polarization.
In the center panel of Figure 8, we consider the ratio of the Faraday rotation to conversion optical depth, τ ρ V /τ ρ Q .In all passing models, this ratio is greater than unity, as we expect for a mildly relativistic plasma.As a consequence, these one-zone models typically produce more circular polarization than linear polarization.This is because in one-zone models where Faraday effects are significant (τ ρ Q > 1, τ ρ V > 1) but absorption is insignificant (τ I < 1), the ratio of circular to linear polarized intensity in the limit of large Faraday depth approaches the ratio of the rotativities: Dexter 2016, for exact one-zone solutions in this limit).As a result, our one-zone models with uniform magnetic field orientations over-produce circular polarization relative to linear polarization, and we must invoke effects not considered in the one-zone models to explain the observed |V| < |P|.
We consider the effects of one such complicationa spatially twisted magnetic field -in the right panel of Figure 8. Faraday conversion only converts linear polarization in the U Stokes parameter into V, so a spatial rotation of the projected B-field can work with or against Faraday rotation to enhance or suppress the produced circular polarization.We add a constant rotation of the magnetic field direction at a spatial frequency ω to our one-zone models, and we solve for the final circular polarization fraction V/I as a function of ω in three models: one with a low conversion 20 µas For moderate rates of field rotation in the plane of the sky, only the model with the highest Faraday depth produces a constant circular polarization fraction at the values of ω considered.The other models produce V/I that varies rapidly and changes sign as the rotation rate ω passes through the critical frequency ω crit = ρ V /2.Thus, field twist through an inhomogenous emission region can have a significant impact on both the magnitude and sign of the observed circular polarization (e.g.Ricarte et al. 2021;Mościbrodzka et al. 2021).We also show results for the same models artificially setting the intrinsic circular polarization emission to zero (j V = 0, dashed lines).In all three cases, this change has a minor effect on the produced circular polarization fraction, again confirming that Faraday conversion dominates the production of circular polarization for plasma parameters appropriate for M87 * .
The fact that one-zone models overproduce V relative to linear polarization and the fact that a changing field geometry in the emission region has significant impact on V/I motivates consideration of more complex models (see, e.g., Goddi et al. 2021, for a two-zone model of M87 * 's Faraday rotation).We proceed next to consider circular polarization in GRMHD simulation images, which self-consistently include the effects of an inhomogenous emission region, field twist, and special and general relativistic redshift and parallel transport.

GRMHD image libraries
The 3D GRMHD simulations used in this work were first considered in Paper V and Paper VIII.We use the code ipole (Noble et al. 2007;Mościbrodzka & Gammie 2018) to perform general relativistic radiative transfer (GRRT) following the methodology outlined in Paper VIII and Wong et al. (2022).We briefly summarize our library generation here, and refer to Appendix H and these previous works for more details.
All images assume a fixed observing frequency of 230 GHz, BH mass of M BH = 6.2 × 10 9 M ⊙ , and distance of 16.9 Mpc. 5 The fluid density in the simulations is scaled to reproduce an average flux density of F ν = 0.5 Jy for each model, following the observed compact flux density in April 2017 (Paper IV).The observing inclination is tilted such that the approaching jet (parallel to the BH spin axis) is inclined at 17 • with respect to our line of sight.
Our simulation library probes 5 free parameters: (i) the magnetic field state, either a strong field "magnetically arrested disk" (MAD) (Bisnovatyi-Kogan & Ruzmaikin 1974; Igumenshchev et al. 2003;Narayan et al. 2003) or weak field "standard and normal evolution" (SANE) (Narayan et al. 2012;Sądowski et al. 2013), (ii) the BH spin a * ∈ {−0.94, −0.5, 0, 0.5, 0.94}, where 5 Note that the values of the M87 * BH mass and distance used in the GRMHD library (M BH = 6.2 × 10 9 M ⊙ ,D = 16.9Mpc) following Paper V; Paper VIII) are slightly different than the EHT's measured value from the ring diameter (M BH = 6.5 × 10 9 M ⊙ , D = 16.8Mpc, Paper VI) which we adopt in other Sections of this paper.Because we do not rely on the image size as a constraint on our models, we do not expect the 5% difference in mass and 0.5% difference in distance between the values adopted in our GRMHD simulations and Paper VI measurements to affect our interpretation.
a negative sign denotes retrograde accretion, (iii-iv) R high ∈ {1, 10, 20, 40, 80, 160} and R low ∈ {1, 10} which modulate the ion-to-electron temperature ratio in different regions, and (v) the magnetic field polarity, which is either aligned or anti-aligned with respect to the disk angular momentum on large scales.Each model is imaged at a cadence of 5 GM/c 3 for a total duration ranging from 2500 to 5000 GM/c 3 depending on the model.In total, we compute 184796 image snapshots from ten GRMHD simulations.For each snapshot we compute a set of polarimetric observables to compare with the data and score our models.Figure 9 displays a random selection of snapshots visualized in circular polarization as a function of spin, magnetic field state, and R high .All images have been rotated such that the oncoming jet is projected at a position angle of 288 • east of north.In the first three rows, we visualize the images in Stokes V at high resolution in symmetric logarithmic scale with three decades in dynamic range.We find a wide variety of morphologies in near-horizon Stokes V images, including sign flips in almost every snapshot.Some images, such as the SANE a * = 0 models, exhibit "noise-like" regions with rapid sign flips among adjacent pixels.This is equivalent to randomly oriented ticks in linear polarization: circular polarization in these regions is scrambled due to either large Faraday rotation6 or Faraday conversion depths.In the bottom three rows, we visualize these snapshots in linear scale with Stokes I contours, blurred with a Gaussian with a 20 µas FWHM.At EHT resolution, GRMHD simulations again predict a wide variety of morphologies including ubiquitous sign flips.

The importance of magnetic field polarity
This is the first paper in this series to consider the magnetic field polarity as a free parameter.We flip the magnetic field polarity in post-processing using the same GRMHD snapshots, since the equations of ideal GRMHD are invariant to a sign flip in the magnetic field vector.The equations of polarized GRRT are not invariant to sign of the magnetic field vector, however.Reversing the polarity of the magnetic field reverses the sign of ρ V (Faraday rotation) and j V (circularly polarized emission), but it does not affect ρ Q (Faraday conversion).Consequently, flipping the magnetic field direction does not necessarily simply reverse the sign of circular polarization across the image.
Throughout, we describe the magnetic field as "aligned" if its polarity is parallel to the angular momentum of the disk on large scales, or "reversed" if this polarity is anti-parallel.Note that the magnetic field structure can be complicated and turbulent in the nearhorizon emission region, especially in retrograde disks, so a magnetic field aligned with the angular momentum on large scales is not necessarily trivially aligned on event horizon scales.In Figure 10, we visualize a snapshot of the MAD a * = −0.94R low = 10 R high = 160 model with both magnetic field polarities.The top row is presented in linear scale blurred with a 20 µas FWHM Gaussian kernel, while the bottom row is presented in symmetric logarithmic scale with 2 decades of dynamic range in intensity.With a polarity aligned with the angular momentum of the disk as in previous work, the unresolved circular polarization fraction v net = 2.7 % vastly exceeds the upper limit of |v net | < 0.8 % from ALMA observations (Goddi et al. 2021).However, flipping the magnetic field polarity reverses the sign of circular polarization in a significant portion of the image, reducing v net to 0.7% and allowing the model to pass.In fact, this snapshot simultaneously passes all polarimetric constraints considered in this work with the reversed magnetic field configuration.
Figure 10 illustrates that it is not easy to predict which regions of a given image change upon a reversal of the magnetic field direction.We expect that regions dominated by intrinsic synchrotron emission should trivially flip sign, while regions dominated by Faraday conversion may remain unchanged unless Faraday rotation is significant along those geodesics.We further explore the effect of flipping the magnetic field polarity on linear polarization in Appendix I.While there are noticeable differences in the distribution of the ∠β 2 parameter (Palumbo et al. 2020) across all GRMHD models, the effect is less dramatic for linear polarization metrics than it is for circular.In Paper VIII, five observational metrics were used to score GRMHD models of M87 * 's accretion flow -(1) the unresolved linear polarization fraction |m| net ; (2) the unresolved circular polarization fraction |v| net ; (3) the image-average linear polarization fraction ⟨|m|⟩; (4) the amplitude |β 2 | and (5) the phase ∠β 2 of the second azimuthal Fourier coefficient of the linear polarization pattern (see Palumbo et al. 2020, Paper VIII).In this paper, we add one additional constraint to this set from our observational results in Section 4: an upper limit on the circular polarization fraction on EHT scales: ⟨|v|⟩ < 3.7 %.We summarize the observational constraints used to score GRMHD models in Table 3. 7In Paper VIII, we found that many GRMHD models which could satisfy both resolved and unresolved linear polarization constraints could also self-consistently satisfy the upper limit on the unresolved circular polarization, |v net | < 0.8 %.In the left panel of Figure 11, we plot histograms of v net for our models, now including flipped magnetic field configurations, which were not included in the original analysis of Paper VIII (their Figure 8).Both MAD and SANE models are capable of passing the upper limit on v net , but SANE models are more likely to fail.In the right panel of Figure 11, we plot histograms of the spatially resolved circular polarization ⟨|v|⟩, and we over-plot our allowed region from the imaging results in Section 4: ⟨|v|⟩ < 3.7 %.SANE models produce the largest values of ⟨|v|⟩, which extend to nearly 10 % in some cases.We find that 87 % of the images that fail our new upper limit on ⟨|v|⟩ are SANE.

GRMHD simulation scoring
We use two different methods for applying observational constraints to our GRMHD images, as in Paper VIII: 1. Simultaneous scoring: If a single model snapshot simultaneously passes every polarimetric constraint, then it passes.Otherwise, it is rejected.This method is relatively strict, as for a model to pass at least one image must satisfy all constraints.See Section 5.2 of Paper VIII for more detail on the simultaneous scoring procedure.
2. Joint scoring: We compute χ 2 j,k statistics for each of the six metrics j for each GRMHD snapshot k (Paper VIII, equation 17).We then compute a likelihood L j for each metric by calculating the fraction of snapshots where χ 2 j,k > χ 2 j,data .The final model likelihood is the product of the individual likelihoods from each metric:L = j L j .This method is relatively lenient, as no single image is required to satisfy all constraints simultaneously.See Section 5.3 of Paper VIII for more detail on the joint scoring procedure.
In Paper VIII, these methods produced slightly different results, but both methods favored MAD models.
Our new model scoring results from both methods are shown in Figure 12.This Figure is an update to Figure 12 of Paper VIII, updated for both new observational results (the addition of the upper limit on ⟨|v|⟩) and new theoretical models (the set of reversed-polarity GRMHD snapshots).As in Paper VIII, we find that MAD models are strongly favored over their SANE counterparts.
In the simultaneous scoring results (left panel of Figure 12), we find that a slight majority of the passing snapshots have aligned magnetic polarities with respect to the disk angular momentum vector.Furthermore, we find that 83 % of passing snapshots have magnetic field polarities aligned with the spin of M87 * , pointed away from us.This is due not to Stokes V constraints, but rather from ∠β 2 , for which the distributions shift slightly upon flipping the magnetic field polarity (see Appendix I).While potentially interesting, especially for constraining the origins of the magnetic field (e.g., Contopoulos & Kazanas 1998;Contopoulos et al. 2022), upon examining Figure 26 in Appendix I, it appears that this preference for an aligned field arises from the fact that ∠β 2 happens to be Faraday rotated out of the observed range more often in reversed field models than in aligned field models for the few spins that we sampled.It is possible that this effect may disappear if spin is sampled more densely, as ∠β 2 depends strongly on spin in GRMHD models (Palumbo et al. 2020;Emami et al. 2023a;Qiu et al. 2023).
Only a small improvement on the upper limit on ⟨|v|⟩ would have been necessary to rule out some currently passing snapshots.Among the passing snapshots, we find 0.30 % ≤ ⟨|v|⟩ ≤ 2.8 %, with a median value of 0.48 %.If the upper limit had instead been ⟨|v|⟩ < 1 %, then 204 snapshots would have passed, down from 288.If the upper limit had been ⟨|v|⟩ < 0.5 %, then only 35 snapshots would have passed, all of which would have been MAD.
In Figure 13, we visualize in circular polarization five of the snapshots that pass simultaneous constraints.Each of them are MAD models, of increasing spin, described in the figure caption.Some, but not all, of the Stokes V morphologies are dominated by a dipole at In Appendix J, we discuss the scoring results for both scoring strategies in more detail.We find that the new constraint on ⟨|v|⟩ has no effect on the results of the si-multaneous scoring method and only a slight effect on the results from the joint scoring method.The largest driver of quantitative differences in the scoring results from Paper VIII is the inclusion of reversed magnetic field models, while updated radiative transfer coefficients and slight differences in the snapshots included in ray-tracing also play a minor role.

Origins of Stokes V in passing models
To better understand the mechanisms by which Stokes V is generated in passing GRMHD models, we conducted a series of tests using our radiative transfer code ipole.By artificially switching off certain radiative transfer coefficients, we can isolate their effects on the results on the resulting images.
We present these tests in detail in Appendix K.In summary, we find that Faraday conversion is the dominant production mechanism for Stokes V in passing GRMHD models.Faraday rotation also plays a critical role, as the linear polarization that is later converted into circular polarization is scrambled on scales smaller than the EHT beam, leading to beam depolarization both in the linear and circular polarization images.Intrinsic or direct emission of circular polarization is subdominant to Faraday conversion, though it is not negligible in all models.

Extensions to GRMHD models
We caution that as with the previous papers in this series, constraints on astrophysical models from EHT results come with systematic uncertainties due to limits on our parameter space.As with previous theoretical interpretation papers in this series, we have limited our scoring to ideal GRMHD simulations with either perfectly aligned or anti-aligned disks, populated with ionelectron plasma with purely thermal distribution functions.
While our R-β prescription for setting the electron-toion temperature ratio broadly describes general trends as a function of plasma β, simulations with electron heating exhibit substantial zone-to-zone scatter and additional features that our prescription cannot reproduce (Mizuno et al. 2021;Dihingia et al. 2023;Jiang et al. 2023).On the other hand, simulations including radiative cooling produce denser and cooler accretion disk mid-planes, which may help produce Faraday rotation (Yoon et al. 2020).
Tilted disks are known to imprint signatures on the total intensity image (Chatterjee et al. 2020), and additional signatures may occur in polarization.Circularly polarized emission is also sensitive to the content of the plasma and its distribution function, which we briefly explore in Appendix L. While beyond the scope of our analysis, we expect that the polarized data published in this series will continue to constrain models probing additional aspects of BH and plasma astrophysics in future work.

SUMMARY AND CONCLUSION
We have presented an analysis of the circular polarization in M87 * at 1.3 mm wavelength from 2017 EHT observations, the first analysis of circular polarization in any source using EHT data.By examining the difference between EHT closure phases in right and left circular polarization measurements, we find firm evidence for a weak, astrophysical (non-instrumental) circular polarization signal in M87 * on event-horizon scales.
To further analyze these data, we developed five different approaches to image analysis, testing each method on a suite of synthetic data from GRMHD simulations.When applied to the observations of M87 * , these methods all find a moderate degree of image circular polarization ( < ∼ 4 %), which is consistent with expectations from the low degree of circular polarization ( < ∼ 1 %) seen in unresolved observations of M87 * with other facilities such as ALMA.The details of our reconstructed images vary considerably among the five methods, indicating that the circular polarization structure is sensitive to choices in the imaging and calibration.Overall, we find that the structure cannot be reliably inferred without additional data or stronger assumptions.From our set of image reconstructions, we establish an upper limit on the (unsigned) resolved circular polarization fraction at the EHT's resolution of 20µas: ⟨|v|⟩ < 3.7%.
Our results provide new constraints for models of the central supermassive black hole in M87 and its environment.We apply these constraints to a large library of images produced from ray-traced GRMHD simulations, which span a broad range of values for the black hole spin, weak (SANE) and strong (MAD) magnetic fields, updating our constraints relative to Paper VIII.We consider both prograde and retrograde accretion flows, and both aligned and anti-aligned poloidal magnetic fields relative to the angular momentum of the accretion disk.
We again find that strongly magnetized MAD models remain favored over weakly magnetized SANEs.As in Paper VIII and Qiu et al. (2023), our scoring prefers spin values roughly around -0.5 and 0.0, driven largely by our constraints on ∠β 2 .That is, the toroidal morphology of the linear polarization ticks favor models with substantial poloidal fields.Our scoring also favors models with larger ion-to-electron temperature ratio R high .That is, the relatively low polarization fraction favors models that contain relatively dense and cool electrons to perform Faraday depolarization.For the model snapshots which simultaneously pass observational constraints, we find that Faraday conversion is typically more important than the circular polarization inherent to synchrotron emission in generating the circular polarization that we observe.This is consistent with the conclusions of parsec-scale studies of other active galactic nuclei at lower frequencies (e.g.Jones 1988;Wardle et al. 1998;Bower et al. 1999;Gabuzda et al. 2008).Faraday rotation also serves an indirect role in limiting the circular polarization fraction in many models, scrambling much of the linear polarization that could otherwise be converted into circular.
Theoretical models of M87 * exhibit significant variability.In Appendix M, we explore this variability in circular polarization metrics in passing GRMHD models.Future M87 * observations may present more favorable conditions for Stokes V imaging (Figure 31).Thus, continued polarimetric monitoring of M87 * on horizon scales will allow us to place better constraints on physical models of the 230 GHz emitting region.Sagittarius A*, the black hole in the galactic center recently imaged by the EHT (Event Horizon Telescope Collaboration et al. 2022a), exhibits a much larger v net ≈ −1 % that may make prospects of Stokes V imaging more promising for this source (Bower et al. 2018;Goddi et al. 2021;Wielgus et al. 2022a), although its much more rapid time variability presents separate challenges (Event Horizon Telescope Collaboration et al. 2022b,c;Wielgus et al. 2022b).
Our results complete the analysis of EHT observations of M87 * in 2017, revealing the 230 GHz emission encircling the apparent shadow of the black hole in all four Stokes parameters.More recent EHT observations, with additional telescopes, significantly improved sensitivity, and higher observing frequencies, will provide crucial improvements to reveal unambiguous structure in the circular polarization and to characterize its variability.APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (Germany), ESO, and the Onsala Space Observatory (Sweden).The SMA is a joint project between the SAO and ASIAA and is funded by the Smithsonian Institution and the Academia Sinica.The JCMT is operated by the East Asian Observatory on behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of China, Chinese Academy of Sciences, and the National Key Research and Development Program (No. 2017YFA0402700) of China and Natural Science Foundation of China grant 11873028.Additional funding support for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada.The LMT is a project operated by the Instituto Nacional de Astrófisica, Óptica, y Electrónica (Mexico) and the University of Massachusetts at Amherst (USA).The IRAM 30-m telescope on Pico Veleta, Spain is operated by IRAM and supported by CNRS (Centre National de la Recherche Scientifique, France), MPG (Max-Planck-Gesellschaft, Germany), and IGN (Instituto Geográfico Nacional, Spain).The SMT is operated by the Arizona Radio Observatory, a part of the Steward Observatory of the University of Arizona, with financial support of operations from the State of Arizona and financial support for instrumentation development from the NSF.Support for SPT participation in the EHT is provided by the National Science Foundation through award OPP-1852617 to the University of Chicago.Partial support is also provided by the Kavli Institute of Cosmological Physics at the University of Chicago.The SPT hydrogen maser was provided on loan from the GLT, courtesy of ASIAA.This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442.XSEDE Stampede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N.XSEDE JetStream resource at PTI and TACC was allocated through AST170028.This research is part of the Frontera computing project at the Texas Advanced Computing Center through the Frontera Large-Scale Community Partnerships allocation AST20023.Frontera is made possible by National Science Foundation award OAC-1818253.This research was done using services provided by the OSG Consortium (Pordes et al. 2007;Sfiligoi et al. 2009) supported by the National Science Foundation award Nos.2030508 and 1836650.Additional work used ABACUS2.0, which is part of the eScience center at Southern Denmark University, and the Kultrun Astronomy Hybrid Cluster (projects Conicyt Programa de Astronomia Fondo Quimal QUIMAL170001, Conicyt PIA ACT172033, Fondecyt Iniciacion 11170268, Quimal 220002).Simulations were also performed on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, on the HazelHen cluster at the HLRS in Stuttgart, and on the Pi2.0 and Siyuan Mark-I at Shanghai Jiao Tong University.The computer resources of the Finnish IT Center for Science (CSC) and the Finnish Computing Competence Infrastructure (FCCI) project are acknowledged.This research was enabled in part by support provided by Compute Ontario (http://computeontario.ca), Calcul Quebec (http://www.calculquebec.ca), and Compute Canada (http://www.computecanada.ca).
The EHTC has received generous donations of FPGA chips from Xilinx Inc., under the Xilinx University Program.The EHTC has benefited from technology shared under open-source license by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER).The EHT project is grateful to T4Science and Microsemi for their assistance with hydrogen masers.This research has made use of NASA's Astrophysics Data System.We gratefully acknowledge the support provided by the extended staff of the ALMA, from the inception of the ALMA Phasing Project through the observational campaigns of 2017 and 2018.We would like to thank A. Deller and W. Brisken for EHT-specific support with the use of DiFX.We thank Martin Shepherd for the addition of extra features in the Difmap software that were used for the CLEAN imaging results presented in this paper.We acknowledge the significance that Maunakea, where the SMA and JCMT EHT stations are located, has for the indigenous Hawaiian people.

A. STRATEGIES FOR POLARIMETRIC GAIN CALIBRATION
The fiducial multi-source procedure for polarimetric gains calibration, described in Subsection 2.3, was employed for the whole series of EHT papers on the M87 * observations performed in 2017, as well as the multiyear M87 * variability study presented in Wielgus et al. (2020), and the 3C 279 quasar study (Kim et al. 2020).
Figure 14 shows an example of how we estimated the phase of G R/L for the LMT using this fiducial multisource approach.The figure shows visibility phase differences between RR * and LL * components recorded on the ALMA-LMT baseline, following the a priori correction for the field angle rotation.A near constant in time phase residual of about -157 • is clearly present, consistent between observing days and for multiple observed sources.This residual is interpreted as the negative phase of the LMT instrumental gain ratio G R/L and subsequently corrected in the data.Similarly, G R/L gains corresponding to other sites were corrected using the baseline from LMT to ALMA on multiple days and sources.
An alternative approach is based on calibrating complex G R/L to RR * = LL * , that is, to Ṽ = 0.The latter strategy, which is simpler but more aggressive was used in the CASA-based rPICARD pipeline (Janssen et al. 2019), where the phase difference between RR * and LL * visibilities is minimized through a polarimetric gains calibration.The Ṽ = 0 self-calibration constitutes a very robust approach that remains valid for the recovery of the remaining Stokes parameters as long as Ĩ ≫ Ṽ.In addition to its use in the rPICARD pipeline, we have made a choice to employ this calibration variant for the HOPS data in the EHT papers on Sagittarius A* (Event Horizon Telescope Collaboration et al. 2022a), Centaurus A (Janssen et al. 2021), J1924-2914 (Issaoun et al. 2022), and NRAO530 (Jorstad et al. 2023).
Note that self-calibration (multiplicative gaincalibration) to Ṽ = 0 is not equivalent to removing the baseline signatures of Ṽ = 0.In particular, this alternative calibration does not affect closure signatures, such as those presented in Figure 2.This is different than the EHT M87 * data set released with Paper III in 2019, 8 where all signatures of Ṽ were deliberately wiped out by replacing RR * and LL * on all baselines with their average, only retaining the Stokes Ĩ total intensity information intact.
Despite certain advantages, self-calibration to Ṽ = 0 may render circular polarization more difficult to recover For the actual estimate of a constant G R/L phase, 5 days and 10 sources were used.The origin of the small residual phases may be instrumental; however, they do not exhibit clear source-independent features.Alternatively, residuals may be caused by the presence of a small, source-specific circular polarization.and for that reason the multi-source polarimetric gains calibration procedure was selected as better suited for the conservative upper limits presented in this paper.In Figure 15 we present the version of the Figure 1, but for the Ṽ = 0 self-calibration data set.Significant decrease of Ṽ S/N is clearly visible, indicating, that high S/N detections seen in Figure 1 could be related to systematic errors of the polarimetric gain calibration, rather than to robust signatures of circular polarization.S/N of the Stokes Ĩ, on the other hand, is mostly insensitive to the G R/L calibration scheme.

B. CLOSURE PHASE DIFFERENCES ON ALMA TRIANGLES
Figure 2 in the main text shows the closure-phase differences for the triangle ALMA-SMT-PV.For this triangle, an average deviation from zero at the level of ∼ 5 σ is detected, which we interpret as evidence of the pres-ence of a spatially-resolved (and asymmetric) Stokes V brightness distribution in M 87*.
There is, though, a total of ten antenna triangles where ALMA is included, which are related to the closure phases with highest S/N.In Figure 16, we show the closure-phase differences of scan-averaged visibilities for all these triangles.There are hints of departures from zero in some cases and for some time ranges (e.g., ALMA-SMT-LMT toward the end of the observations, as well as ALMA-APEX-LMT and ALMA-LMT-PV), but there is no clear detection above 3 − 4 σ.The detection of Stokes V from all the ALMA-related antenna triangles is, therefore, only tentative if we use the closure phase differences.
C. CLOSURE TRACES Closure traces are a set of complex closure quantities defined on quadrangles that are insensitive to stationbased corruptions (see Broderick & Pesce 2020).Of particular relevance to the current discussion of Stokes V is the fact that this includes the right-and left-hand complex station gains and leakage terms.From these, a combination can be constructed on a single quadrangle, the conjugate closure trace product that differs from unity only in the presence of non-trivial polarization structure (see Broderick & Pesce 2020).For this reason, these conjugate closure trace products provide a robust direct test for source polarization.
Figure 17 compares the phases of the conjugate closure trace product generated on the APEX-ALMA-LMT-SMT and APEX-SMT-LMT-ALMA quadrangles from the image reconstructions in Figure 3 directly to the 2017 April 11 low-band data.9All reconstructions are broadly consistent with the conjugate closure trace products, similar to the results found in Appendix B of Paper VII.
Because the closure traces are invariant to rotations on the Poincaré sphere, the conjugate closure trace products cannot isolate the contribution from the Stokes V maps.However, we also show by dotted lines the conjugate closure trace products for each reconstruction in Figure 3 after setting the Stokes V to zero, indicating the gross magnitude of the impact of circular polarization.For all reconstruction methods, the difference attributable to the Stokes V map is small.For Themis and DMC, the difference is small in comparison to the range of conjugate closure trace product phases spanned by the image posterior, implying that reconstruction of the Stoke V maps is strongly dependent on the assumed calibration priors, consistent with what is found in Section 4.

D. METHOD SUMMARIES
In this Appendix, we describe each image reconstruction method introduced in Section 3 and summarized in Table 1.

D.1. DIFMAP
The CLEAN algorithm (Högbom 1974) is a traditional VLBI imaging method which performs imaging via inverse modeling in all four Stokes parameters, I, Q, U, and, V for visibility data with dual polarization.The method is implemented in the DIFMAP (Shepherd 2011) and Astronomical Imaging Process System (AIPS Greisen 2003) software packages.Imaging in DIFMAP is an iterative process using a set of CLEAN windows until the residual χ 2 between data and a model reach a minimum value.For total intensity imaging, the process involves self-calibration of phase and amplitude of visibilities with the model at different steps of iterations, which results in calibrated visibility data and a file containing a model of I delta functions characterized by the amplitude and position on the image.For linear and circular polarization imaging, no additional self-calibration procedure is performed, and DIFMAP generates models of Q,U, and V delta functions.
For EHT data, to obtain a Stokes V image we first perform an amplitude and phase calibration of R and L visibility data with the Stokes I image via the AIPS task CALIB, assuming that (RR * + LL * )/2 = Ĩ.Then we apply the D-term correction to the visibility data and repeat the amplitude and phase calibration with the image assuming again (RR * + LL * )/2 = Ĩ.For this study we use the M87 * D-terms, which are specified in Tables 3 & 5 of Appendix A in Paper VII, applying them through an antenna table via AIPS tasks TBOUT and TBIN.At this stage, the calibrated data files are written out to DIFMAP for visual inspection and further editing, if needed.
After reading the edited visibility data back to AIPS, we proceed with the final R/L gain calibration, which is necessary to remove any residual instrumental circular polarization.To do that, we apply an additional amplitude and phase calibration of R and L visibility data separately with the I image using the CALIB task, assuming that RR * = Ĩ and LL * = Ĩ.This is necessary because DIFMAP does not calibrate R and L polarization hands separately.The resulting solution (SN) table ideally includes only a small residual correction for R and L that represents their gain offset for each antenna in the array.This SN table is processed outside AIPS to edit the R and L amplitude and phases with the assumption that the R/L gains of the ALMA station equals unity and hence any R/L ratio seen by ALMA is ideally intrinsic to the source.Under this assumption, the R and L visibility data for all antennas are modified accordingly so that the average R/L ratios better match the R/L ratio seen by ALMA.We note that this approach corrects only for a R/L offset for each antenna and is insensitive to possible time variations of the R/L gains.The modified SN table is read back to AIPS with the task TBIN and used to produce the final data files (uvfits), using tasks CLCAL and SPLIT, which should be calibrated for both circular and linear polarization.Such calibrated uv-data are reloaded in DIFMAP and imaged in Stokes V using the set of windows employed for the Stokes I image.This is the procedure that we have used to produce Stokes V images for the 2017 M87 * data and synthetic data analyzed in this work.We have used the noise level (RMS) multiplied by a factor of 3 of Stokes I and V images as the uncertainties of the total and circular polarized intensities, respectively, to derive uncertainties of parameters presented in Appendix E and Section 4.

D.2. polsolve
The polsolve algorithm is also based on CLEAN, but for polarimetric calibration uses the full Measurement Equation (e.g., Smirnov 2011), including (second order) non-linear corrections to the instrumental polarization (see Martí-Vidal et al. 2021, for full details).While the original polsolve algorithm was not designed to include a Stokes V source model in the fitting of the instrumental polarization, the flexibility of CASA makes it possible to account for circular polarization by manually filling the model column of the measurement set and asking polsolve to use it in the fitting.This strategy corresponds to the "polarimetric self-calibration" approach described in detail in Martí-Vidal et al. (2021).
For the recovery of the Stokes V image with polsolve, we have followed this iterative polarimetric self-calibration approach, starting from the data calibrated as described in Paper VII. First, we have run the CASA-based CLEAN on the four Stokes parameters, by using a common mask based on the Stokes I image.The parameters used in this CLEAN are the same as those summarized for polsolve in Paper VII.Each run of (full-polarization) CASA-based CLEAN was followed by a (polarization-independent) self-calibration (i.e., gain solutions were computed based on Stokes Ĩ only, using RR * and LL * ).Then, a D-term fitting with polsolve was performed (this time, the full visibility matrix was used).This procedure was iterated until convergence (i.e., until the changes in the antenna gains and D-terms between consecutive iterations were negligible).After convergence, a final CLEAN image was generated.This approach was the same for all the data presented in this publication (i.e., all the synthetic data tests and all the EHT observations).

D.3. DMC and Themis
Themis (Broderick et al. 2020a) and DMC (Pesce 2021) produced Stokes V maps concurrently with the linear polarization maps presented in Paper VII, and details about the analyses are provided there.Here, we briefly summarize both methods and emphasize the differences that are most relevant for imaging circular polarization.
Themis and DMC formulate the imaging problem as one of Bayesian posterior exploration over the combined space of both image structures and station-based calibration quantities (i.e., complex gains and leakage terms).The output of both codes is a set of MCMC samples from the joint posterior distribution over both the full-Stokes image structure and the calibration quantities.
For the purposes of Stokes V reconstruction, the most salient difference between Themis and DMC is in their treatment of the relative right-and left-handed station gain quantities.Whereas Themis holds that the rightand left-handed station gains are identical, DMC explores potentially large differences between them.With the exception of a reference station -chosen here to be ALMA -that has both right-and left-hand gain phases fixed to be zero-valued at all times, DMC independently models the right-and left-hand complex gain quantities for each scan and each station.

D.4. Geometric modeling
Geometric modeling differs from imaging methods because the reconstructed source structure is restricted to the space defined by a family of simple geometric models.Geometric modeling has been used to infer total intensity source structure parameters for M87 * and Sgr A * , generally showing excellent consistency with the results of imaging (Event Horizon Telescope Collabora-tion et al. 2019f, 2022c).The restricted parameter space in geometric modeling is generally easier to constrain than that of many-parameter image reconstructions and their associated regularizers.Hence, geometric modeling provides a useful complement to imaging, especially to analyze observations with sparse uv-sampling.However, geometric modeling results will not be accurate if the true source structure is not represented in the model space, so care must be taken when modeling sources where the source structure is unknown and unconstrained from other methods.
In this work, we fit a full-Stokes m-ring model (R23) to EHT M87 * data from 2017.The total intensity m-ring model (Johnson et al. 2020;Event Horizon Telescope Collaboration et al. 2022c) consists of a ring with diameter d, with azimuthal brightness variations that are decomposed in Fourier modes.The image is given by Here, ρ and φ are polar sky coordinates centered on the ring, β k are the complex Fourier coefficients, F is the total flux density (we set β 0 ≡ 1), and the total intensity m-ring order m I is the maximum nonzero Fourier coefficient.By setting a higher m I , the total number of model parameters increases and an increasingly complex azimuthal structure can be modeled.Because the total intensity is real, The m-ring is given a finite width by convolving I(ρ, φ) with a circular Gaussian with FWHM α.A useful property of this model is that both the image and its corresponding visibility function are straightforward to calculate analytically (see, e.g., R23).
The m-ring model can be straightforwardly generalized to include all Stokes parameters.In linear polarization, the azimuthal variations are set by {β P,k }.Because the linear polarization image is complex, β P,k and β * P,k are not related by conjugation symmetry.In circular polarization, the azimuthal variations are set by {β V,k }.Because the circular polarization image is real, β V,k = −β * V,k .For our m-ring fits to synthetic and real EHT data, we use the geometric modeling module in eht-imaging (Chael et al. 2016(Chael et al. , 2018)).In this work, we fit an (m I = 3, m P = 3, m V = 1) m-ring to the parallelhand closure amplitudes and closure phases.Posterior exploration is done with the dynesty sampler (Speagle 2020).We first fit the total intensity and linear polarization structure, after which we fix the linear polarization parameters and fit the total intensity and circular polarization structure simultaneously.Because the closure products do not constrain the net circular polarization β V,0 ≡ v net , we fix this parameter to the ALMA-only estimate of −0.3 % (Goddi et al. 2021) for the EHT M87 * data, and to the ground-truth values for the synthetic datasets.Our fitting procedure and assumptions are motivated in more detail in R23, where additional tests are presented, including varying the circular polarization m-ring order m V , varying the assumed net circular polarization fraction (m c ≡ β V,0 ), and varying the fitted data products (e.g., using closure quantities versus using parallel-hand complex visibility ratios).
E. SYNTHETIC DATA TESTS In this Section, we test the ability of our image reconstruction and geometric modeling methods to recover resolved circular polarization from an idealized, highcircular-polarization "convention test" as well as representative images from GRMHD simulations of M87 * .

E.1. Synthetic data generation
As in Paper VII, we generated synthetic observations on the M87 * 2017 April 11 low-band (u, v) coverage using the eht-imaging software package.After sampling the image Fourier transforms on EHT baselines, the synthetic data were corrupted with time-variable, stationbased systematic phase and gain offsets, time-stable polarimetric leakage terms, and thermal noise.
Unlike in the Paper VII synthetic data, in this paper we do not assume the complex gains in the right and left circular polarizations are equal: G R ̸ = G L .Instead, we introduce time-variable offsets in both the amplitude and phase of the ratio G R/L .The magnitude of these offsets and their characteristic timescales were motivated by the a-priori limits described in Subsection 2.3.For most sites, we sample the R/L amplitude gain offset from a normal distribution centered at unity with a standard deviation of 20 %, and the R-L phase offset from a zero-centered normal distribution with a standard deviation of 10 deg.The exceptions are ALMA (where G R = G L ) and the SMA, where a priori limits on the phase offset motivate a use of a 40 deg standard deviation.The phase and amplitude gain offsets are correlated in time (as are the overall amplitude and phase gain terms), and we set the correlation timescales for the offsets to be larger than the correlation timescales for the overall amplitude (2 hr vs 1 hr) and phase gains (24 hr vs 0.5 hr).

E.2. Convention test results
Before imaging the GRMHD models presented in Subsection E.3, we first tested our methods on an optimistic "convention test"; a strongly circularly polarized (v ≈ 10 %) double Gaussian source with 1.2 Jy total flux density.In this test, instead of a GRMHD snapshot, a synthetic source image was constructed from two Gaussians with FWHM≈ 20µas separated by a distance d ≈ 50µas.The total flux density of the double source was 1.2 Jy, and the fractional circular polarization of each Gaussian component was high, V pk /I ≈ 10%.We show the results of this initial test in Figure 18.
All participating methods (the m−ring modeling group did not participate in this test) were able to re- Images are plotted in the same manner as in Figure 19.
cover the two-component structure as well as the approximate degree of circular polarization in each component, even with systematic gain offsets G R/L present in the synthetic data at similar levels to the GRMHD data sets discussed in the following section.This test demonstrates that our analysis methods would have been able to recover consistent, accurate resolved Stokes V images from 2017 EHT observations of M87 * if the intrinsic circular polarization brightness was ≈ 5 − 10 times higher than observed in 2017.

E.3. GRMHD results
To test our imaging methods on more realistic source models for M87 * , we selected three GRMHD images from Paper VIII that pass all the observational constraints from Paper VII.We sorted these images by the value of ⟨|v|⟩ after blurring to 20 µas resolution.We then chose three snapshot images from across this distribution to represent weak, average, and near-maximal ⟨|v|⟩ from the distribution of passing Paper VIII images.All of these images have low image integrated circular polarization |v net | < 0.3 %; images that have stronger polarization at 20 µas scales have significant regions of both positive and negative circular polarization that cancel in unresolved observations.
In Table 4, we summarize these three models, including the magnetic field configuration and spin of their underlying GRMHD simulation, the net 230 GHz circular polarization fraction v net , and the average resolved fraction ⟨|v|⟩.
Figure 19 shows the GRMHD ground truth models and reconstructions from our three synthetic datasets (Appendix E), for each of the five imaging or model- Table 4. Summary of synthetic data models.
ing methods described in Section 3.Because all three GRMHD synthetic data models were chosen from the set of passing Paper VIII models, all feature relatively similar total linear polarization fractions and EVPA structure, similar to those observed in M87 * by the EHT in Paper VII.For all models, the total intensity and linear polarization structure (top rows) are generally recovered well by all methods.Because the synthetic data included both large overall amplitude gains G on certain stations, the total flux density is not completely constrained by these datasets (as is the case for the real EHT 2017 M87 * data).The Stokes V reconstructions (bottom rows) show significant variation between the methods, and most methods have difficulty reconstructing the groundtruth features in circular polarization.
All methods struggle the most with model 01, the weakly polarized case (⟨|v|⟩ ≈ 0.5 %).In this model, no method correctly recovers the spatial distribution of the sign of circular polarization, and due to the presence of residual relative gains G R/L in the data some methods significantly overproduce the total Stokes V brightness in the image.In contrast, for the "high-polarization" model 03, (⟨|v|⟩ ≈ 4 %), all methods are able to recover the resolved distribution of V somewhat more accurately.Themis, polsolve, and m-ring modeling all recover the correct dipolar Stokes V structure of the image; DMC and DIFMAP produce mostly one sign of V corresponding to the more strongly polarized negative  component in the image.The results of model Model 02 are more mixed, with some methods (DMC, Themis, and DIFMAP) accurately recovering the sign and approximate magnitude of V at the total intensity peak, while m-ring modeling and polsolve see an overall dipole in the same orientation as the ground truth.No image reconstruction accurately captures the weak circular polarization in the northern half of the image.
In all three synthetic models, m-ring modeling does the best job of consistently recovering the level of V in the image and the orientation of the dipolar structure across the ring.The limits on structural complexity baked into the m V = 1 model likely helps the method lock on to the overall orientation of the source in these synthetic data tests, which all feature an underlying dipolar structure in V.While many GRMHD images show this dipolar structure at EHT resolution, some show more complicated patterns (see Section 5).More details and other model choices applied to m-ring modeling of these synthetic data sources are presented in R23).
Figure 20 shows the recovery of the image-integrated and resolved circular polarization fractions for all datasets and methods.The low image-integrated circular polarization fraction is recovered decently well by all methods on the three sources.An exception is DIFMAP, which recovers a relatively large negative integrated circular polarization fraction.The m-ring modeling method is constrained by construction to produce the ground-truth v net .The recovered resolved circular polarization ⟨|v|⟩ shows larger differences between methods and larger deviations from the ground truth.For the high circular polarization case 03, both DMC and Themis have the ground truth value of ⟨|v|⟩ within their 2σ error bars, while the polsolve and DIFMAP single-image results are off by ≈ 1 %.In the moderate and low polarization models 02 and 01, all methods except m-ring modeling produce results biased significantly upward from the ground truth.This could indicate that residual uncertainty in the G R/L gains is entering into over-producing the resolved circular polarization in these methods with weak intrinsic V and low V S/N in the data, even after self-calibration during imaging.Furthermore, ⟨|v|⟩ is biased upwards by uncertainty in the underlying image reconstructions (see Appendix G).Interestingly, m-ring modeling performs the best on datasets 01 and 02, accurately recovering ⟨|v|⟩ in both cases, while overproducing ⟨|v|⟩ in the "easier" case 03.This difference in recovery accuracy could be driven by the limited freedom of the m-ring model in the source structures it can recover; the ground-truth structure for model 03 is less dipolar and ring-like than for models 01 and 02, with both total intensity and polarized flux concentrated in the South.
Figure 21 shows the recovered amplitude gain ratios |G R/L | for selected stations on the three synthetic datasets.Note that Themis does not fit for these gain ratios, and m-ring modeling does not fit for any gains as it directly fits to closure quantities.The trends of |G R/L | as a function of time are generally recovered by all methods, although systematic offsets up to ∼ 10 % do occur.DMC is the only method that returns uncertainty estimates on these gain ratios; the DMC results are nearly all within 2σ of the ground truth.The inability of all methods to accurately recover these gain ratios even in the model with the most circular polarization (03) indicates that residual uncertainty in the solution for G R/L may significantly affect our V imaging results for M87 * .and 03 (bottom).From left to right -gain ratios are reported for ALMA, the LMT, the SMT, and the SMA sites.The actual applied gain ratios in the synthetic data set (motivated by the limits discussed in Subsection 2.3) are displayed in black.DMC gains are reported with 2σ error bars.
F. M87* RESULTS FROM CASA-CALIBRATED DATA Figure 22 presents Stokes V imaging results using data reduced with the CASA rPICARD pipeline (Janssen et al. 2019).The CASA-reduced data calibrated the amplitude gain ratios G R/L assuming V = 0, while the fiducial HOPS-reduced data used to produce the results in Figure 4 and Figure 5 used a multi-day and multi-source fit to calibrate G R/L .By producing images with these data, we test if our main conclusions regarding the level of the resolved V and (in)consistency of the resolved structure in different reconstructions depend strongly on the choices and assumptions made in data reduction, particularly in the way G R/L are calibrated before imaging.
Three of the imaging pipelines used in this paper (DIFMAP, Themis, and m-ring model fitting) produced Stokes V images using the CASA data for both high and low band on April 5 and April 11.The results presented in Figure 22 show that the resulting Stokes V images from the CASA-calibrated are still inconsistent across reconstruction methods and observing days.Interestingly, the reconstructions using CASA data show somewhat more consistency across the high and low band images than in the reconstructions using the fiducial EHT-HOPS data (Figure 4, Figure 5).This may be a result of the fact that the rPICARD pipeline combines information across the bands for the signal stabilization (Janssen et al. 2022) of the data, while the bands are treated entirely separately in the EHT-HOPS reduction.R23 presents an analysis of the CASA data m-ring reconstructions compared to the EHT-HOPS data reconstructions in more detail and without any G R/L amplitude gain calibration.

G. BIASES IN RESOLVED CIRCULAR POLARIZATION FRACTION
For an image composed of N res pixels or resolution elements, each of which has angular area ∆A, then Equation 5 can be expressed as a sum The mean of this distribution is µ i = σ i 2/π, which is nonzero.The mean of ⟨|v|⟩ will then be given by If each pixel or resolution element has a similar noise variance σ 2 i , then we can see that E [⟨|v|⟩] will be proportional to N res .Precisely, the value of the resolved circular polarization fraction depends on the number of pixels or resolution elements contained in the image, which means that it will be sensitive to things such as the size of a blurring kernel or the field of view of the image.

H. ADDITIONAL GRMHD LIBRARY GENERATION DETAILS
The GRMHD models in this paper were computed using the code HARM (Gammie et al. 2003;Noble et al. 2006) and were first presented in Paper V. Simulations were initialized with an idealized Fishbone-Moncrief torus of gas in hydrostatic equilibrium (Fishbone & Moncrief 1976), within which the magnetorotational instability (Balbus & Hawley 1991, 1998) naturally develops, driving accretion.Our GRMHD simulations are mainly characterized by two parameters.The first is the black hole spin a * ranging from −1 < a * < 1 with positive/negative values representing alignment/anti-alignment between the engines angular momenta.We only consider models where the BH angular momentum and the disk angular momentum on large scales is either perfectly aligned or anti-aligned, where anti-alignment is denoted with a negative sign.The second is the absolute magnetic flux in its dimensionless form ϕ ≡ Φ BH ( Ṁ r 2 g c) −1/2 , where Φ BH is the magnetic flux' magnitude traversing one hemisphere of the event horizon (see, e.g.Porth et al. 2019) and the mass accretion rate Ṁ set over the event horizon.The relative strength of the magnetic flux near the horizon is determined by the dimensionless quantity ϕ and affects the dynamics of the accretion flow solutions.This led to a division of the solutions into two categories: the Magnetically Arrested Disk (MAD) state which highly affects the accretion dynamics ϕ ≥ 50 (e.g.Tchekhovskoy et al. 2011;McKinney et al. 2012;Narayan et al. 2022), and the Standard and Normal Evolution (SANE) state ϕ ≈ 1 − 5 (in Gaussian units).SANE simulations were produced on a 288 × 128 × 128 grid, a fluid adiabatic index γ = 4/3, and a boundary of r out = 50 r g .MAD simulations adopt a 384×192×192 resolution, γ = 13/9, and a domain of r out = 10 3 r g .
To produce images via general relativistic radiative transfer (GRRT), we use the code ipole (Mościbrodzka & Gammie 2018).First, null geodesics are calculated backwards from the camera through the source.Then, radiative transfer is calculated forwards accounting for polarized synchrotron emission and self-absorption, as well as Faraday rotation and conversion.We use radiative transfer coefficients appropriate for a thermal relativistic plasma (Leung et al. 2011), with the updated fitting functions from Marszewski et al. (2021).The fast-light approximation is made during ray-tracing.
In ideal GRMHD, the fluid density, magnetic field strength, and internal energy can be rescaled during post-processing, and we do so to achieve an average flux density of F ν ≈ 0.5 Jy for our images at ν = 230 GHz (e.g., Event Horizon Telescope Collaboration et al. 2019e).In addition to the fluid scaling and BH parameters, there are three other free parameters we explore at this stage.Two of them, R low and R high are associated with the ion-to-electron temperature ratio, where the electron temperature is defined as In order for the description in Equation (H6) to be valid, the ions are assumed to be non-relativistic while electrons be relativistic.Adopting a prescription in Mościbrodzka et al. ( 2016), the ion-to-electron temperature ratio R is given by where the components describe the local plasma β = p gas /p mag with p gas = (γ − 1)u gas , p mag = B 2 /8π, and further u gas being the total internal energy density.To regulate the ratio in Equation (H7), our models explore R low = 1, 10 and R high = 1, 10,20,40,80,160 as in Paper VIII.Such a parameterization to treat the model as thermal flow is motivated by the fact that mean free path is much larger than the Debye-Length and Larmor radius of the jet's plasma.This assumption together with Equation (H7) mimic collisionless plasma properties and allows to associate the electron heating with its magnetic properties (see Paper VIII).

I. FIELD REVERSAL AND LINEAR POLARIZATION
As we discuss in 5.2, reversing the polarity of the magnetic field has significant effects on circular polarization, both on resolved and unresolved scales.Our previous studies including linear polarization of this source only included models wherein the poloidal magnetic field and the angular momentum of the disk are aligned.Here, we briefly explore the effect of flipping the magnetic field on linear polarization metrics.Overall, we find the effect is minor and does not have a strong impact on linear polarization constraints.
In Figures 23-26, we plot histograms of observable quantities considered in Paper VIII for the aligned case (as in previous work) on the left, and for the anti-aligned case on the right.As expected, most of the linear polarization metrics are insensitive to this flip.The quantity which shows the most significant, albeit minor, effect is ∠β 2 .As explored in more detail in Emami et al. (2023a), on average, Faraday rotation imprints a small shift in the EVPA pattern corresponding to the line-of-sight magnetic field direction.For example, since we are viewing M87 * nearly face-on, a poloidal field pointed towards the observer will on average Faraday rotate EVPA ticks in the counter-clockwise direction, and β 2 is also rotated counter-clockwise.Similarly, a poloidal field pointed away from the observer will on average rotate ticks β 2 clockwise.More complicated behavior can even arise if differences in the Faraday depth among different emitting regions cause summation or cancellation depending on the degree of Faraday rotation.
Ultimately, although flipping the polarity of the magnetic field has an effect on the distributions of ∠β 2 , it has little effect on our linear polarization constraints.SANE models do not appear to change significantly in this metric.The MAD models appear to shift slightly in ∠β 2 when the field polarity is reversed.Prograde models have a slight increase, whereas retrograde models have a slight decrease in ∠β 2 .This does not affect the simultaneous snapshot scoring model results, but causes an increased preference for retrograde MAD models in the joint scoring method.With multi-frequency data, a sign flip in the magnetic field polarity should naturally lead to a sign flip in the rotation measure (e.g.Contopoulos et al. 2022), but we do not consider rotation measure in this work.polarimetric constraints already naturally pass this constraint.This is perhaps not surprising, since the upper limit on v net reported by Goddi et al. (2021) and already considered in Paper VIII is more stringent than that on ⟨|v|⟩ obtained in this work.Only 261/184796 snapshots (0.14 %) fail the ⟨|v|⟩ constraint but pass the v net constraint.Producing a low v net but a high ⟨|v|⟩ requires regions of high resolved circular polarization fraction of both positive and negative sign to nearly cancel across an image.Ultimately when combined with the linear polarization constraints, ⟨|v|⟩ adds no discerning power.Consequently, model predictions for M87 * remain consistent with Paper VIII.
The minor differences in the simultaneous scoring histogram (left panel of Figure 12) compared to Paper VIII can be attributed not to the additional constraint on ⟨|v|⟩, but rather to the additional parameter probed in this paper, the magnetic field polarity, as well as slight differences in time sampling of the images and updates to radiative transfer coefficients.
To better understand how Stokes V alone affects our model scoring, Figure 27 plots the fraction of snapshots of a given model that pass Stokes V constraints without considering the others.We find that Stokes V constraints are more likely to rule out SANE models than MADs, since SANE models are more likely to overproduce v net .The ⟨|v|⟩ constraint alone does not rule out many models.Passing fractions appear similar for both field configurations, and thus they have been combined in the figure.

J.2. Joint Scoring
As shown in Fig. 12, the joint scoring method favors MAD R low = 10 models that are not highly spinning.The inclusion of reversed-field images greatly increases the relative likelihood for MAD spin −0.5 models, particularly for ∠β 2,LP .The joint scoring results appear qualitatively different from that given in Paper VIII, which is likely due to differences in the image library sets rather than just the ⟨|v|⟩ constraint.
The inclusion of the ⟨|v|⟩ constraint removes small likelihoods for MAD spin 0,−0.5, R low = 1 models but otherwise does not qualitatively change the likelihoods.Note that the midpoint for the ⟨|v|⟩ allowed range is 1.85 %, or half of the conservative upper limit.The results from the joint scoring method are more sensitive to this choice than is the case in simultaneous scoring.Using an aggressive upper limit on ⟨|v|⟩ of 2 % from Themis (Table 2), we find that qualitatively the joint scoring results remain largely the same, but some R low = 1 models now are included as those models have lower ⟨|v|⟩.A lower limit on ⟨|v|⟩ at 1 % tends to increase the likelihood of the MAD spin +0.5, R low = 10 models.
Comparisons between the image library sets between this paper and Paper VIII can yield insight into the differences in the joint scoring method.For example, the relatively lower likelihoods for the MAD spin +0.5 models in Fig. 12 can in part be attributed to only scoring images at 17 • or 163 • inclination, while Paper VIII included models at inclinations ±5 • from 17 • or 163 • .Models from Paper VIII scored using only 17 • or 163 • inclination cause the relative likelihood for MAD spin +0.5 models to decrease.The time sampling in the library sets can also influence the model scoring.In this paper, we sample uniformly for each model at 5M, while in Paper VIII the number of snapshots per model is kept fixed to 200, with the sampling cadence varying.Using the sampling in Paper VIII causes a slight increase to the likelihoods for the prograde MAD models.Moreover, the library set in Paper VIII excludes emission from within 5% of the event horizon, while the current image library consists of emission from further inside.While this only affects models that are very optically thin (MAD models), it can change the EVPA structure and thus the β 2,LP modes.The β 2,LP modes appear to be the most discrepant between the two library sets, and most linear and circular polarization constraints otherwise remain consistent.
Here, we present in more detail the test summarized in Subsection 5.5 to determine the physical origin of Stokes V emission in physical models.We investigate the full set of 288 snapshots that simultaneously pass all polarimetric constraints.In ipole, our GRRT code, it is straightforward to probe the effects of each physical effect by setting the appropriate radiative transfer coefficient to 0. Here, the most relevant coefficients are intrinsic circularly polarized emission (j V ), Faraday conversion (ρ Q ), and Faraday rotation (ρ V ).For details of these coefficients, we refer the reader to Mościbrodzka & Gammie (2018).
In this experiment, we re-image all the passing snapshots but each time, turn off one of these coefficients.Selecting one representative snapshot per passing model, we compare both v net and ⟨|v|⟩ against the full radiative transfer solution and plot the results in Figure 28.Since significant cancellation occurs in the image plane for Stokes V, large differences to the image do not necessarily correspond to large changes in v net .Therefore, ⟨|v|⟩ is somewhat more informative, and robust to cancellations, although cancellations within the Gaussian kernel of 20 µas still occur.
From Figure 28 we note 3 main results:     Each column corresponds to a comparison between the full radiative transfer (RT) simulation and the same calculation with one coefficient turned off: jV =0 (intrinsic emission) for column 1, ρQ=0 (Faraday conversion) for column 2, ρV =0 (Faraday rotation) for column 3.Each tick on the y axis corresponds to a model for which at least one snapshot passes all the polarimetric constraints.A single snapshot from each of these models has been plotted at each tick.
important not only for reducing the linear polarization fraction, but also for indirectly reducing the circular polarization fraction.In these models, Faraday rotation partly scrambles the linear polarization that is converted into circular.
In summary, Faraday effects are critical for Stokes V generation in our models, while intrinsic emission plays a sub-dominant, though non-negligible, role.
L. ALTERNATIVE PLASMA MODELS Throughout this work, we have modeled a pure ionelectron plasma with thermal electron distributions.Here, we explore our sensitivity to the details of the plasma by considering the existence of electron-positron pairs as well as non-thermal electron distribution functions.In brief, our Stokes V images are indeed sensitive to these details, motivating future studies in these areas.O'Dell 1977;Jones 1988;Wardle et al. 1998).For ionic plasma models, all of the radiative transfer coefficients are believed to be important for producing the Stokes V image.On EHT scales, the few recent studies on this topic find that increasing the pair fraction can alter both resolved Stokes V images at 230 GHz as well as the evolution as a function of frequency (Anantua et al. 2020;Emami et al. 2021Emami et al. , 2023b)).
In Figure 29, we test a single snapshot of the MAD a = +0.5 R high = 80 R low = 10 aligned field model, performing polarized ray-tracing with an increasing positronto-electron number density ratio, denoted as f in each panel.To obtain a given value of f , we inject electronpositron pairs to each cell, representing a simplistic scenario in which a large fraction of pairs are produced at number densities directly proportional to the number densities of pre-existing electrons.Since this would drastically increase the number density of emitting leptons in each successive panel, we also find a new value of the density normalization for each value of f , such that the total flux is always 0.5 Jy as in the original snapshot.
For a given value of f , the radiative transfer coefficients are modified via (e.g., MacDonald & Marscher 2018;Emami et al. 2021) capturing the effects described above.Figure 29 illustrates that the circularly polarized morphology of this snapshot clearly and strongly evolves with f , consistent with previous works.However, given the inherent diversity of Stokes V structures among the library, there is no known signature in Stokes V that clearly indicates the presence of pairs.Interestingly, we find that ⟨|v|⟩ increases monotonically in the perfect-resolution images, which can be attributed to drastically falling Faraday rotation depth (to be discussed in Anantua et al. in prep.).However, due to cancellations within the EHT beam, ⟨|v|⟩ does not necessarily increase monotonically at EHT resolution.In summary, we are sensitive to the content of the emitting plasma, but we are unaware of a clear signature distinguishing a pair plasma from an ionic plasma on event horizon scales where both intrinsic emission and Faraday conversion are important.While not explored here, images originating from GRMHD are also surprisingly sensitive to the atomic composition of a given ionic plasma (Wong & Gammie 2022).Identifying clear signatures to distinguish different models of plasma content on event horizon scales remains a topic of ongoing research.

L.2. Non-thermal electrons
Throughout this work, we consider only models with thermal electron distribution functions (eDF).However, it is believed that non-thermal particle acceleration can occur due to magnetic reconnection, MHD turbulence and collective plasma modes.At present, a consensus model for the presence of non-thermal electrons has not be found, but here we explore two physically motivated implementations.For a single snapshot, we compare their thermal eDF images to images produced assuming "kappa" models, which feature a thermal core and a highenergy non-thermal tail (Vasyliunas 1968).Such distributions have been observed in the solar wind (Decker & Krimigis 2003;Pierrard & Lazar 2010) and occur naturally in particle-in-cell (PIC) simulations of particle acceleration in magnetized plasmas (Kunz et al. 2016).We try both a constant κ = 5 model, where κ − 1 = 4 is the power law index of the high-energy tail, and variable kappa model where κ = κ(σ, β) as fit by the PIC simulations of Ball et al. (2016).Due to limitations of the fitting functions used, any time the prescription would assign κ > 7, the code instead uses a radiative transfer coefficient appropriate for a thermal eDF.
In Figure 30, we spot-check our sensitivity to our assumption of thermal electrons by ray-tracing a GRMHD snapshot with three different assumptions for the electron distribution function.A new plasma density scale is found for each image to match the total flux of the image with thermal electrons, 0.7 Jy.This snapshot corresponds to the MAD a = +0.5 R high = 80 R low = 10 aligned field model.In the top row, we plot images blurred to EHT resolution, while in the bottom row we plot perfect-resolution images of Stokes V in symmetric logarithmic scale.Compared to the leftmost image assuming thermal electrons, we find very little changes in the variable kappa model shown in the central panel.This is broadly consistent with our findings in Event Horizon Telescope Collaboration et al. (2022d) for Sgr A * in total intensity.Thus, at least for this snapshot, exchanging a thermal distribution for a physically motivated non-thermal electron distribution has very little effect on Stokes V.However, we report dramatic differences when switching to a constant kappa model with a value of κ = 5.As found in many other studies, non-thermal electrons make the image noticeably larger and more diffuse (e.g., Özel et al. 2000;Mao et al. 2017;Fromm et al. 2022;Ricarte et al. 2023).Intriguingly, although the plasma density is decreased only by a factor of 4 compared to the thermal model, the Stokes V signal almost entirely vanishes.This may be due to the fact that Faraday conversion is caused by the coldest relativistic electrons, which occur at a smaller fraction in κ models by definition.The standard deviations of models with different magnetic field polarities have been averaged together.Different magnetic field states are plotted in different rows, different spins are plotted in different columns, R high is plotted on the x-axis of each panel, and the color and shape of the markers encodes R low .In particular, models with R low = 10 are more variable than those with R low = 1, and a similar but weaker trend occurs with R high .

M. STOKES V VARIABILITY IN PASSING MODELS
In Figure 31, we consider each model's distribution of v net over time and plot its standard deviation, which we denote as σ vnet .We average together the standard deviations for each magnetic field polarity, since their distributions can be disjoint.We find that models with larger R low are more variable in v net , with a similar but weaker trend in R high .Although not shown, we find that the variability of ⟨|v|⟩ is qualitatively similar.Some models have σ vnet in excess of our upper limit of 0.008, suggesting that future observations during more favorable conditions may result in higher S/N detections of Stokes V on EHT baselines and hence more robust image reconstructions of the source.

Figure 1 .
Figure 1.Stokes Ṽ visibility S/N in scan-averaged EHT 2017 data, as a function of the projected baseline length.Data from both frequency bands, and all observing days are shown.No systematic uncertainties, like the imperfect calibration of the gain ratios G R/L , were accounted for.Hence, the plotted S/N represents upper limits on the Ṽ detections.Gray points in the background indicate the S/N of Stokes Ĩ detections on the corresponding baselines.

Figure 2 .
Figure 2. Closure phases observed on the triangle ALMA-SMT-PV triangle during M87 * observations on 2017 Apr 5-11 in low band (left column) and in high band (right column).Top: closure phases of scan-averaged visibilities for all epochs, RR * in red; LL * in blue.Bottom: difference of closure phases between RR * and LL * .The zero level of the closure difference (i.e., no Ṽ detected) is marked with a black dashed line.A light green band shows the RR * − LL * difference, inferred by-band, under the constant difference assumption.

Figure 3 .
Figure 3. Reconstructions of 2017 EHT M87 * data from April 11, low band.The top row shows total intensity images from all reconstruction methods in grayscale and fractional linear polarization in colored ticks as in Paper VII.The second row shows the same grayscale total intensity image overlaid with colored ellipses indicating the total polarization fraction |mtot| = √ Q 2 + U 2 + V 2 /I.The size of each ellipse indicates the total polarized brightness; the orientation of each ellipse indicates the linear EVPA, and axis ratio indicates the relative fraction of circular polarization.The color of each ellipse indicates the sign of circular polarization.

Figure 4 .
Figure 4. Circular polarization imaging results from 2017 EHT observations of M87 * on April 5 (top two rows) and April 6(bottom two rows).Images of circular polarization on these consecutive days are expected to be nearly identical, as is seen in total intensity and linear polarization.The top/bottom row in each pair shows results from imaging the high/low-band data.In each panel, total intensity is indicated in the colored linear-scale contours, and the Stokes V brightness is indicated in the diverging colormap, with red/blue indicating a positive/negative sign.The colorbar ranges are fixed in both plots (and in Figure5).For posterior exploration methods (DMC, Themis, m-ring fitting), the posterior-average image is shown.All images are blurred with the same 20 µas FWHM Gaussian, shown with the black inset circle in the upper-left panels.

Figure 5 .
Figure 5.The same as Figure 4, but for 2017 EHT observations of M87 * on April 10 (top two rows) and April 11 (bottom two rows).

Figure 6 .
Figure 6.M87 * imaging results combining days and bands.The top row shows results from three methods on a data set combining April 5 and 6 observations, both low and high band.The bottom row shows corresponding results from combining April 10 and 11 observations, low and high band.For Themis and DMC we show the posterior-average image.Images are plotted in the same manner as in Figure 4.

Figure 7 .
Figure7.Image-integrated statistics from M87 * images.The left panel shows the net circular polarization fraction vnet computed from each method for the eight EHT M87 * datasets and the right panel shows the average resolved circular polarization fraction ⟨|v|⟩ computed after blurring each image with a 20 µas kernel.The results from the posterior exploration methods are presented with the median value and 2σ error bars (Note that m-ring modeling strictly enforces |vnet| = 0).The DIFMAP and polsolve results are derived from the single fitted image and standard error propagation from measurements of the off-source image RMS in V and I.In the left panel, the limits on vnet from ALMA observations reported inGoddi et al. (2021)  (used to constrain GRMHD models in Paper VIII) are indicated in the gray shaded region.In the right panel, the upper limit on ⟨|v|⟩ derived in this work (⟨|v|⟩ < 3.7 %, Table2) is indicated by the dashed horizontal line.

Figure 8 .
Figure8.Circular polarization properties of passing Paper VIII one-zone models.(Left) Distribution of the Faraday conversion optical depth τρQ in passing models.In all passing models τρQ > 1, indicating that most circular polarization is likely produced by Faraday conversion.(Center) Distribution of the ratio of the Faraday rotation to Faraday conversion optical depths.In all cases, τρV > τρQ, indicating that with a constant field orientation in the emission region, circular polarization will dominate over linear polarization in these models.(Right) The average fractional circular polarization between 5 rg and 10 rg in one-zone models with a rotating magnetic field direction along the line of sight, as a function of the angular rotation frequency ω.We show three different models: a model with low Faraday conversion depth (blue), a model with median conversion depth (orange), and a model with high conversion depth (green).Dashed lines show corresponding results for one-zone models with no intrinsic emission of circular polarization, jV = 0.

Figure 9 .
Figure 9.A random selection of representative snapshots from our GRMHD image library.The color scales for each snapshot are normalized individually.The first three rows are presented at native resolution in symmetric logarithmic scale with three decades in dynamic range shown to better visualize faint features.The bottom three rows plot Stokes I in contours and Stokes V in color after blurring with a 20 µas FWHM Gaussian, both in linear scale.Models exhibit a wide variety of morphologies and almost always show sign reversals, both at perfect resolution and EHT resolution.depth(τ ρ Q ≈ 1 : |B| = 11 G, Θ e = 7, n e = 10 4 cm −3 ); one with a median conversion depth (τ ρ Q ≈ 5 : |B| = 5 G, Θ e = 8, n e = 10 5 cm −3 ); and one with a large conversion depth (τ ρ Q ≈ 200 : |B| = 9 G, Θ e = 3, n e = 4 × 10 6 cm −3 ).

Figure 10 .
Figure 10.Example GRMHD snapshot (MAD a * = −0.94R low = 10 R high = 160) plotted with both magnetic field configurations, aligned (left) and reversed (right).The top panels show the images blurred to EHT resolution, and the bottom panels show the images at their native resolution.As shown in the left panel, this snapshot fails simultaneous polarimetric constraints with the aligned field configuration, overproducing vnet.However, as shown in the right panel, flipping the magnetic field polarity produces some oppositely-signed regions that reduce |vnet|.Flipping the field has no effect on the total intensity image.

Figure 11 .
Figure 11.Number density distribution of net circular polarization vnet fraction (left panel) and image-averaged fractional circular polarization ⟨|v|⟩ (right panel) with an aligned and anti-aligned magnetic field respecting all images in the M87 * library.Allowed inferred ranges for ALMA-only data (vnet) and measured values of reconstructed polarimetric images of M87 * reconstructions (⟨|v|⟩) are limited by the dashed lines.

Figure 12 .Figure 13 .
Figure 12.Scoring results for all the models using the simultaneous scoring method (left panel) and the joint distribution method (right panel).Passing fraction or relative likelihoods are summed over R high and ⃗ B field alignment.MADs remain favored across both methods.

8Figure 14 .
Figure14.Example of a multi-source, multi-day estimation of the G R/L phase for the LMT, using RR * and LL * high band visibilities measured on the ALMA-LMT baseline.For the actual estimate of a constant G R/L phase, 5 days and 10 sources were used.The origin of the small residual phases may be instrumental; however, they do not exhibit clear source-independent features.Alternatively, residuals may be caused by the presence of a small, source-specific

Figure 15 .Figure 16 .
Figure 15.S/N of the Ṽ observations as a function of the projected baseline length, for the data set self-calibrated to Ṽ = 0. Gray points in the background correspond to the S/N of Stokes Ĩ detections.

Figure 17 .
Figure 17.Phases of the conjugate closure trace product constructed on the APEX-ALMA-LMT-SMT and APEX-SMT-LMT-ALMA quadrangles from the 2017 April 11 low-band observations.Colored lines show the same quantity computed from the mean images in Figure 3, with the 2σ regions for Themis and DMC indicated by bands.For comparison, the conjugate closure trace products when the Stokes V maps are ignored are shown by the corresponding dotted lines for each method.

Figure 18 .
Figure 18.Imaging results from the convention test data set, consisting of a highly polarized (V pk /I ≈ 10%) double source.Images are plotted in the same manner as in Figure 19.

Figure 19 .
Figure 19.Imaging results from the three synthetic data sets in Table 4.The first row shows the total intensity in grayscale and the fractional linear polarization in the colored ticks.The left column is the ground truth simulation image.From left to right, the next columns show the posterior average images from DMC and Themis, the final CLEAN images from polsolve and DIFMAP, and the posterior average of the m−ring model fits.All images are blurred with same 20µas FWHM circular Gaussian beam.The second row shows the total intensity image with 8 contour levels on a linear scale and the Stokes V structure in a a diverging colormap.Rows one and two show the results for model 01 (low resolved circular polarization), rows three and four show the results for model 02 (moderate resolved circular polarization), and rows five and six show the results for model 03 (high resolved circular polarization).Note that the colorbar for V has different maximum values in rows two, four, and six as the GRMHD simulation images become more polarized.

Figure 20 .
Figure 20.Image-integrated statistics from the results of the synthetic data tests presented in Figure19.The left panel shows the net circular polarization fraction vnet computed from each method for the three synthetic datasets and the right panel shows the average resolved circular polarization fraction ⟨|v|⟩ computed after blurring each image with a 20 µas kernel.The results from the posterior exploration methods are presented with the median value and 2 σ error bars (note that m-ring modeling strictly enforces that vnet is equal to the ground truth).The CLEAN based methods produce only one image and so only one value of each metric is reported, but 2 σ error bars for DIFMAP and polsolve are shown from standard error propagation using the off-source residuals in V and I.The ground truth value of each statistic is indicated with the black star and horizontal dashed line.

Figure 21 .
Figure 21.Derived gain ratios |GR|/|GL| from different imaging methods applied to synthetic dataset 01 (top), 02 (middle),and 03 (bottom).From left to right -gain ratios are reported for ALMA, the LMT, the SMT, and the SMA sites.The actual applied gain ratios in the synthetic data set (motivated by the limits discussed in Subsection 2.3) are displayed in black.DMC gains are reported with 2σ error bars.

Figure 22 .
Figure 22.M87 * Stokes V imaging results using data reduced with the CASA-based rPICARD pipeline (Janssen et al. 2019) and pre-calibrated with the assumption V = 0.The top two rows show results from three methods (DIFMAP, Themis, and m-ring geometric modeling) on April 5 observations, both low and high band.The bottom row shows corresponding results from April 11 observations.For Themis and m-ring modeling we show the posterior-average image.Images are plotted in the same manner as in Figure 4. where F is the total integrated Stokes I flux density in the image, and i is a subscript that indicates a particular pixel or resolution element.If I is everywhere positive in the image, then |V i /I i | I i = |V i | and we can simplify the expression for ⟨|v|⟩ to be ⟨|v|⟩ = ∆A F Nres

Figure 23 .Figure 24 .
Figure 23.Distribution of image-integrated net linear polarization fraction with an aligned (left panel) and anti-aligned (right panel) magnetic field respecting all images in the M87 * library.Allowed inferred ranges of EHT image reconstructions are limited by the dashed lines.

Figure 25 .Figure 26 .
Figure 25.Distribution of β2 amplitude for an aligned (left panel) and anti-aligned (right panel) magnetic field representing values taken from all images in the EHT M87 * library blurred with a 20 µas beam.Allowed inferred ranges of EHT image reconstructions are limited by the dashed lines.

Figure 27 .
Figure 27.Passing fraction of the vnet and ⟨|v|⟩ constraints for all the models in the simulation library.Aligned and anti-aligned models are scored combined since the individual passing fractions are similar.

Figure 28 .
Figure28.Exploration of the effects of various CP production mechanisms on vnet (top panels) and ⟨|v|⟩ (bottom panels).Each column corresponds to a comparison between the full radiative transfer (RT) simulation and the same calculation with one coefficient turned off: jV =0 (intrinsic emission) for column 1, ρQ=0 (Faraday conversion) for column 2, ρV =0 (Faraday rotation) for column 3.Each tick on the y axis corresponds to a model for which at least one snapshot passes all the polarimetric constraints.A single snapshot from each of these models has been plotted at each tick.

Figure 31 .
Figure31.Standard deviations of the distributions of vnet in GRMHD models over time.The standard deviations of models with different magnetic field polarities have been averaged together.Different magnetic field states are plotted in different rows, different spins are plotted in different columns, R high is plotted on the x-axis of each panel, and the color and shape of the markers encodes R low .In particular, models with R low = 10 are more variable than those with R low = 1, and a similar but weaker trend occurs with R high .

Table 1 .
Summary of the Imaging and Modeling Methods Used polsolve Priors from Paper VII Self-calibration (gains & D-terms).ALMA G R/L ≡ 1 DIFMAP Paper VII values fixed Self-calibration, assuming V = 0. ALMA G R/L ≡ 1

Table 3 .
Observational Constraints Applied to our GRMHD Image Library Note-Most of these constraints are inherited from Paper VII and were previously used to constrain models in Paper VIII.This work adds the new upper limit on ⟨|v|⟩.
NSF, and National Institutes of Natural Sciences of Japan, together with National Research Council (Canada), Ministry of Science and Technology (MOST; Taiwan), Academia Sinica Institute of Astronomy and Astrophysics (ASIAA; Taiwan), and Korea Astronomy and Space Science Institute (KASI; Republic of Korea), in cooperation with the Republic of Chile.The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc. (AUI)/NRAO, and the National Astronomical Observatory of Japan (NAOJ).The NRAO is a facility of the NSF operated under cooperative agreement by AUI.This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC05-00OR22725; the ASTROVIVES FEDER infrastructure, with project code IDIFEDER-2021-086; the computing cluster of Shanghai VLBI correlator supported by the Special Fund for Astronomy from the Ministry of Finance in China; We also thank the Center for Computational Astrophysics, National Astronomical Observatory of Japan.This work was supported by FAPESP (Fundacao de Amparo a Pesquisa do Estado de Sao Paulo) under grant 2021/01183-8.