$v_{\mathrm{net}}\approx-1$ % that may make prospects of Stokes $\mathcal{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.

##  7 Acknowledgements

The Event Horizon Telescope Collaboration thanks the following organizations
and programs: the Academia Sinica; the Academy of Finland (projects 274477,
284495, 312496, 315721); the Agencia Nacional de Investigación y Desarrollo
(ANID), Chile via NCN$19\\_058$ (TITANs), Fondecyt 1221421 and BASAL FB210003;
the Alexander von Humboldt Stiftung; an Alfred P. Sloan Research Fellowship;
Allegro, the European ALMA Regional Centre node in the Netherlands, the NL
astronomy research network NOVA and the astronomy institutes of the University
of Amsterdam, Leiden University, and Radboud University; the ALMA North
America Development Fund; the Astrophysics and High Energy Physics programme
by MCIN (with funding from European Union NextGenerationEU, PRTR-C17I1); the
Black Hole Initiative, which is funded by grants from the John Templeton
Foundation and the Gordon and Betty Moore Foundation (although the opinions
expressed in this work are those of the author(s) and do not necessarily
reflect the views of these Foundations); the Brinson Foundation; “la Caixa”
Foundation (ID 100010434) through fellowship codes LCF/BQ/DI22/11940027 and
LCF/BQ/DI22/11940030; the Fondo CAS-ANID folio CAS220010; Chandra DD7-18089X
and TM6-17006X; the China Scholarship Council; the China Postdoctoral Science
Foundation fellowships (2020M671266, 2022M712084); Consejo Nacional de
Humanidades, Ciencia y Tecnología (CONAHCYT, Mexico, projects U0004-246083,
U0004-259839, F0003-272050, M0037-279006, F0003-281692, 104497, 275201,
263356); the Colfuturo Scholarship; the Consejería de Economía, Conocimiento,
Empresas y Universidad of the Junta de Andalucía (grant P18-FR-1769), the
Consejo Superior de Investigaciones Científicas (grant 2019AEP112); the
Delaney Family via the Delaney Family John A. Wheeler Chair at Perimeter
Institute; Dirección General de Asuntos del Personal Académico-Universidad
Nacional Autónoma de México (DGAPA-UNAM, projects IN112820 and IN108324); the
Dutch Organization for Scientific Research (NWO) for the VICI award (grant
639.043.513), the grant OCENW.KLEIN.113, and the Dutch Black Hole Consortium
(with project No. NWA 1292.19.202) of the research programme the National
Science Agenda; the Dutch National Supercomputers, Cartesius and Snellius (NWO
grant 2021.013); the EACOA Fellowship awarded by the East Asia Core
Observatories Association, which consists of the Academia Sinica Institute of
Astronomy and Astrophysics, the National Astronomical Observatory of Japan,
Center for Astronomical Mega-Science, Chinese Academy of Sciences, and the
Korea Astronomy and Space Science Institute; the European Research Council
(ERC) Synergy Grant “BlackHoleCam: Imaging the Event Horizon of Black Holes"
(grant 610058); the European Union Horizon 2020 research and innovation
programme under grant agreements RadioNet (No. 730562) and M2FINDERS (No.
101018682); the Horizon ERC Grants 2021 programme under grant agreement No.
101040021; the Generalitat Valenciana (grants APOSTD/2018/177 and
ASFAE/2022/018) and GenT Program (project CIDEGENT/2018/021); MICINN Research
Project PID2019-108995GB-C22; the European Research Council for advanced grant
‘JETSET: Launching, propagation and emission of relativistic jets from binary
mergers and across mass scales’ (grant No. 884631); the FAPESP (Fundação de
Amparo á Pesquisa do Estado de São Paulo) under grant 2021/01183-8; the
Institute for Advanced Study; the Istituto Nazionale di Fisica Nucleare (INFN)
sezione di Napoli, iniziative specifiche TEONGRAV; the International Max
Planck Research School for Astronomy and Astrophysics at the Universities of
Bonn and Cologne; DFG research grant “Jet physics on horizon scales and
beyond” (grant No. FR 4069/2-1); Joint Columbia/Flatiron Postdoctoral
Fellowship (research at the Flatiron Institute is supported by the Simons
Foundation); the Japan Ministry of Education, Culture, Sports, Science and
Technology (MEXT; grant JPMXP1020200109); the Japan Society for the Promotion
of Science (JSPS) Grant-in-Aid for JSPS Research Fellowship (JP17J08829); the
Joint Institute for Computational Fundamental Science, Japan; the Key Research
Program of Frontier Sciences, Chinese Academy of Sciences (CAS, grants QYZDJ-
SSW-SLH057, QYZDJSSW-SYS008, ZDBS-LY-SLH011); the Leverhulme Trust Early
Career Research Fellowship; the Max-Planck-Gesellschaft (MPG); the Max Planck
Partner Group of the MPG and the CAS; the MEXT/JSPS KAKENHI (grants 18KK0090,
JP21H01137, JP18H03721, JP18K13594, 18K03709, JP19K14761, 18H01245, 25120007,
23K03453); the Malaysian Fundamental Research Grant Scheme (FRGS)
FRGS/1/2019/STG02/UM/02/6; the MIT International Science and Technology
Initiatives (MISTI) Funds; the Ministry of Science and Technology (MOST) of
Taiwan (103-2119-M-001-010-MY2, 105-2112-M-001-025-MY3, 105-2119-M-001-042,
106-2112-M-001-011, 106-2119-M-001-013, 106-2119-M-001-027,
106-2923-M-001-005, 107-2119-M-001-017, 107-2119-M-001-020,
107-2119-M-001-041, 107-2119-M-110-005, 107-2923-M-001-009,
108-2112-M-001-048, 108-2112-M-001-051, 108-2923-M-001-002,
109-2112-M-001-025, 109-2124-M-001-005, 109-2923-M-001-001,
110-2112-M-003-007-MY2, 110-2112-M-001-033, 110-2124-M-001-007, and
110-2923-M-001-001); the Ministry of Education (MoE) of Taiwan Yushan Young
Scholar Program; the Physics Division, National Center for Theoretical
Sciences of Taiwan; the National Aeronautics and Space Administration (NASA,
Fermi Guest Investigator grant 80NSSC20K1567, NASA Astrophysics Theory Program
grant 80NSSC20K0527, NASA NuSTAR award 80NSSC20K0645); NASA Hubble Fellowship
grants HST-HF2-51431.001-A, HST-HF2-51482.001-A awarded by the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Inc., for NASA, under contract NAS5-26555; the National
Institute of Natural Sciences (NINS) of Japan; the National Key Research and
Development Program of China (grant 2016YFA0400704, 2017YFA0402703,
2016YFA0400702); the National Science Foundation (NSF, grants AST-0096454,
AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433,
AST-1140030, DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509,
OPP-1248097, AST-1310896, AST-1440254, AST-1555365, AST-1614868, AST-1615796,
AST-1715061, AST-1716327, OISE-1743747, AST-1816420, AST-1935980, AST-2034306,
AST-2307887); NSF Astronomy and Astrophysics Postdoctoral Fellowship
(AST-1903847); the Natural Science Foundation of China (grants 11650110427,
10625314, 11721303, 11725312, 11873028, 11933007, 11991052, 11991053,
12192220, 12192223, 12273022); the Natural Sciences and Engineering Research
Council of Canada (NSERC, including a Discovery Grant and the NSERC Alexander
Graham Bell Canada Graduate Scholarships-Doctoral Program); the National Youth
Thousand Talents Program of China; the National Research Foundation of Korea
(the Global PhD Fellowship Grant: grants NRF-2015H1A2A1033752, the Korea
Research Fellowship Program: NRF-2015H1D3A1066561, Brain Pool Program:
2019H1D3A1A01102564, Basic Research Support Grant 2019R1F1A1059721,
2021R1A6A3A01086420, 2022R1C1C1005255); Netherlands Research School for
Astronomy (NOVA) Virtual Institute of Accretion (VIA) postdoctoral
fellowships; Onsala Space Observatory (OSO) national infrastructure, for the
provisioning of its facilities/observational support (OSO receives funding
through the Swedish Research Council under grant 2017-00648); the Perimeter
Institute for Theoretical Physics (research at Perimeter Institute is
supported by the Government of Canada through the Department of Innovation,
Science and Economic Development and by the Province of Ontario through the
Ministry of Research, Innovation and Science); the Princeton Gravity
Initiative; the Spanish Ministerio de Ciencia e Innovación (grants
PGC2018-098915-B-C21, AYA2016-80889-P, PID2019-108995GB-C21,
PID2020-117404GB-C21); the University of Pretoria for financial aid in the
provision of the new Cluster Server nodes and SuperMicro (USA) for a SEEDING
GRANT approved toward these nodes in 2020; the Shanghai Municipality
orientation program of basic research for international scientists (grant no.
22JC1410600); the Shanghai Pilot Program for Basic Research, Chinese Academy
of Science, Shanghai Branch (JCYJ-SHFY-2021-013); the State Agency for
Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa”
award for the Instituto de Astrofísica de Andalucía (SEV-2017- 0709); the
Spanish Ministry for Science and Innovation grant CEX2021-001131-S funded by
MCIN/AEI/10.13039/501100011033; the Spinoza Prize SPI 78-409; the South
African Research Chairs Initiative, through the South African Radio Astronomy
Observatory (SARAO, grant ID 77948), which is a facility of the National
Research Foundation (NRF), an agency of the Department of Science and
Innovation (DSI) of South Africa; the Toray Science Foundation; the Swedish
Research Council (VR); the UK Science and Technology Facilities Council (grant
no. ST/X508329/1); the US Department of Energy (USDOE) through the Los Alamos
National Laboratory (operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the USDOE, contract
89233218CNA000001); and the YCAA Prize Postdoctoral Fellowship.

We thank the staff at the participating observatories, correlation centers,
and institutions for their enthusiastic support. This paper makes use of the
following ALMA data: ADS/JAO.ALMA#2016.1.01154.V. ALMA is a partnership of the
European Southern Observatory (ESO; Europe, representing its member states),
NSF, and National Institutes of Natural Sciences of Japan, together with
National Research Council (Canada), Ministry of Science and Technology (MOST;
Taiwan), Academia Sinica Institute of Astronomy and Astrophysics (ASIAA;
Taiwan), and Korea Astronomy and Space Science Institute (KASI; Republic of
Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory
is operated by ESO, Associated Universities, Inc. (AUI)/NRAO, and the National
Astronomical Observatory of Japan (NAOJ). The NRAO is a facility of the NSF
operated under cooperative agreement by AUI. This research used resources of
the Oak Ridge Leadership Computing Facility at the Oak Ridge National
Laboratory, which is supported by the Office of Science of the U.S. Department
of Energy under contract No. DE-AC05-00OR22725; the ASTROVIVES FEDER
infrastructure, with project code IDIFEDER-2021-086; the computing cluster of
Shanghai VLBI correlator supported by the Special Fund for Astronomy from the
Ministry of Finance in China; We also thank the Center for Computational
Astrophysics, National Astronomical Observatory of Japan. This work was
supported by FAPESP (Fundacao de Amparo a Pesquisa do Estado de Sao Paulo)
under grant 2021/01183-8.

APEX is a collaboration between the Max-Planck-Institut für Radioastronomie
(Germany), ESO, and the Onsala Space Observatory (Sweden). The SMA is a joint
project between the SAO and ASIAA and is funded by the Smithsonian Institution
and the Academia Sinica. The JCMT is operated by the East Asian Observatory on
behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of
China, Chinese Academy of Sciences, and the National Key Research and
Development Program (No. 2017YFA0402700) of China and Natural Science
Foundation of China grant 11873028. Additional funding support for the JCMT is
provided by the Science and Technologies Facility Council (UK) and
participating universities in the UK and Canada. The LMT is a project operated
by the Instituto Nacional de Astrófisica, Óptica, y Electrónica (Mexico) and
the University of Massachusetts at Amherst (USA). The IRAM 30-m telescope on
Pico Veleta, Spain 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.

##  Appendix 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 multi-year 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 multi-source 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 $\tilde{\mathcal{V}}=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 $\tilde{\mathcal{V}}=0$ self-calibration constitutes a very
robust approach that remains valid for the recovery of the remaining Stokes
parameters as long as $\tilde{\mathcal{I}}\gg\tilde{\mathcal{V}}$. 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).

Figure 14: 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

circular polarization.

Note that self-calibration (multiplicative gain-calibration) to
$\tilde{\mathcal{V}}=0$ is not equivalent to removing the baseline signatures
of $\tilde{\mathcal{V}}=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,888https://datacommons.cyverse.org/browse/iplant/home/shared/commons_repo/curated/EHTC_FirstM87Results_Apr2019
where all signatures of $\tilde{\mathcal{V}}$ were deliberately wiped out by
replacing $RR^{*}$ and $LL^{*}$ on all baselines with their average, only
retaining the Stokes $\tilde{\mathcal{I}}$ total intensity information intact.

Figure 15: S/N of the $\tilde{\mathcal{V}}$ observations as a function of the
projected baseline length, for the data set self-calibrated to
$\tilde{\mathcal{V}}=0$. Gray points in the background correspond to the S/N
of Stokes $\tilde{\mathcal{I}}$ detections.  Figure 16: Closure phase
differences between scan-averaged $RR^{*}$ and $LL^{*}$ visibilities, for all
the antenna triangles related to ALMA. All epochs and both bands (empty
symbols for high band; filled symbols for low band) are shown. The zero line
(i.e., no detection of spatially resolved Stokes $\mathcal{V}$) is shown as a
dotted line.

Despite certain advantages, self-calibration to $\tilde{\mathcal{V}}=0$ may
render circular polarization more difficult to recover 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 $\tilde{\mathcal{V}}=0$
self-calibration data set. Significant decrease of $\tilde{\mathcal{V}}$ 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
$\tilde{\mathcal{I}}$, on the other hand, is mostly insensitive to the
$G_{R/L}$ calibration scheme.

##  Appendix 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
$\sim 5\,\sigma$ is detected, which we interpret as evidence of the presence
of a spatially-resolved (and asymmetric) Stokes $\mathcal{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$
$\sigma$. The detection of Stokes $\mathcal{V}$ from all the ALMA-related
antenna triangles is, therefore, only tentative if we use the closure phase
differences.

##  Appendix C Closure Traces

Closure traces are a set of complex closure quantities defined on quadrangles
that are insensitive to station-based corruptions (see Broderick & Pesce,
2020). Of particular relevance to the current discussion of Stokes
$\mathcal{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.999The
insensitivity of the closure traces to the gain ratios $|G_{R/L}|$ results in
similar width $2\sigma$ bands in Figure 17 for Themis and DMC despite their
different treatments of the station gains, and thus very different error bars
in Figure 7. All 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 $\mathcal{V}$ maps. However, we also show by dotted lines the conjugate
closure trace products for each reconstruction in Figure 3 after setting the
Stokes $\mathcal{V}$ to zero, indicating the gross magnitude of the impact of
circular polarization. For all reconstruction methods, the difference
attributable to the Stokes $\mathcal{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 $\mathcal{V}$ maps is strongly dependent on the assumed calibration
priors, consistent with what is found in Section 4.

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\sigma$ regions for Themis and DMC
indicated by bands. For comparison, the conjugate closure trace products when
the Stokes $\mathcal{V}$ maps are ignored are shown by the corresponding
dotted lines for each method.

##  Appendix 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,
$\mathcal{I,Q,U,\text{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 $\chi^{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 $\mathcal{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 $\mathcal{Q}$,$\mathcal{U}$, and
$\mathcal{V}$ delta functions.

For EHT data, to obtain a Stokes $\mathcal{V}$ image we first perform an
amplitude and phase calibration of R and L visibility data with the Stokes
$\mathcal{I}$ image via the AIPS task CALIB, assuming that
$(RR^{*}+LL^{*})/2=\tilde{\mathcal{I}}$. 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=\tilde{\mathcal{I}}$. 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 $\mathcal{I}$ image using the CALIB task, assuming that
$RR^{*}=\tilde{\mathcal{I}}$ and $LL^{*}=\tilde{\mathcal{I}}$. 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 $\mathcal{V}$ using the set of windows employed
for the Stokes $\mathcal{I}$ image. This is the procedure that we have used to
produce Stokes $\mathcal{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 $\mathcal{I}$ and $\mathcal{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 $\mathcal{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 $\mathcal{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 $\mathcal{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 $\tilde{\mathcal{I}}$
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
$\mathcal{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 $\mathcal{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 right-
and 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 Collaboration
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

$\displaystyle\mathcal{I}(\rho,\varphi)$ $\displaystyle=\frac{F}{\pi
d}\delta\left(\rho-\frac{d}{2}\right)\sum_{k=-m_{\mathcal{I}}}^{m_{\mathcal{I}}}\beta_{k}e^{ik\varphi}.$
(D1)

Here, $\rho$ and $\varphi$ are polar sky coordinates centered on the ring,
$\beta_{k}$ are the complex Fourier coefficients, $F$ is the total flux
density (we set $\beta_{0}\equiv 1$), and the total intensity $m$-ring order
$m_{\mathcal{I}}$ is the maximum nonzero Fourier coefficient. By setting a
higher $m_{\mathcal{I}}$, the total number of model parameters increases and
an increasingly complex azimuthal structure can be modeled. Because the total
intensity is real, $\beta_{k}=-\beta_{k}^{*}$. The $m$-ring is given a finite
width by convolving $\mathcal{I}(\rho,\varphi)$ with a circular Gaussian with
FWHM $\alpha$. 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
$\left\\{\beta_{\mathcal{P},k}\right\\}$. Because the linear polarization
image is complex, $\beta_{\mathcal{P},k}$ and $\beta_{\mathcal{P},k}^{*}$ are
not related by conjugation symmetry. In circular polarization, the azimuthal
variations are set by $\left\\{\beta_{\mathcal{V},k}\right\\}$. Because the
circular polarization image is real,
$\beta_{\mathcal{V},k}=-\beta_{\mathcal{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, 2018). In this work, we
fit an ($m_{\mathcal{I}}=3,m_{\mathcal{P}}=3,m_{\mathcal{V}}=1$) $m$-ring to
the parallel-hand closure amplitudes and closure phases. Posterior exploration
is done with the ${\tt 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 $\beta_{\mathcal{V},0}\equiv v_{\rm
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_{\mathcal{V}}$, varying the
assumed net circular polarization fraction ($m_{\rm
c}\equiv\beta_{\mathcal{V},0}$), and varying the fitted data products (e.g.,
using closure quantities versus using parallel-hand complex visibility
ratios).

##  Appendix E Synthetic Data Tests

Figure 18: Imaging results from the convention test data set, consisting of a
highly polarized ($\mathcal{V}_{\rm pk}/\mathcal{I}\approx 10$%) double
source. Images are plotted in the same manner as in 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\mu$as FWHM circular Gaussian beam. The second row shows
the total intensity image with 8 contour levels on a linear scale and the
Stokes $\mathcal{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 $\mathcal{V}$ has
different maximum values in rows two, four, and six as the GRMHD simulation
images become more polarized.

In this Section, we test the ability of our image reconstruction and geometric
modeling methods to recover resolved circular polarization from an idealized,
high-circular-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, station-based 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}\neq 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\approx 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$\approx 20\mu$as separated by a
distance $d\approx 50\mu$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_{\rm pk}/I\approx 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 recover 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 $\mathcal{V}$ images from 2017 EHT
observations of M87∗ if the intrinsic circular polarization brightness was
$\approx 5-10$ times higher than observed in 2017.

###  E.3 GRMHD results

Figure 20: Image-integrated statistics from the results of the synthetic data
tests presented in Figure 19. The left panel shows the net circular
polarization fraction $v_{\rm net}$ computed from each method for the three
synthetic datasets and the right panel shows the average resolved circular
polarization fraction $\langle|v|\rangle$ computed after blurring each image
with a 20 $\mu$as kernel. The results from the posterior exploration methods
are presented with the median value and 2 $\sigma$ error bars (note that
$m$-ring modeling strictly enforces that $v_{\rm net}$ 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 $\sigma$ error bars for DIFMAP and polsolve
are shown from standard error propagation using the off-source residuals in
$\mathcal{V}$ and $\mathcal{I}$. The ground truth value of each statistic is
indicated with the black star and horizontal dashed line.  Figure 21: Derived
gain ratios $|G_{R}|/|G_{L}|$ 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\sigma$ error bars.

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
$\langle|v|\rangle$ after blurring to $20\,\mu$as resolution. We then chose
three snapshot images from across this distribution to represent weak,
average, and near-maximal $\langle|v|\rangle$ from the distribution of passing
Paper VIII images. All of these images have low image integrated circular
polarization $|v_{\rm net}|<0.3\,\%$; images that have stronger polarization
at $20\,\mu$as scales have significant regions of both positive and negative
circular polarization that cancel in unresolved observations.

Our three models are: the “weak polarization” case (model 01) with an average
circular polarization fraction at $20\,\mu$as resolution of
$\langle|v|\rangle$ $\approx 0.5\,\%$; the “moderate polarization” case (model
02) with $\langle|v|\rangle$ $\approx 2\,\%$; and the “strong polarization”
case (model 03) with $\langle|v|\rangle$ $\approx 4\,\%$. 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_{\rm net}$, and the average resolved fraction
$\langle|v|\rangle$.

Model |  $B-$Field | $a*$ |  $v_{\rm net}$ (%) |  $\langle|v|\rangle$ (%)
---|---|---|---|---
01 | MAD | -0.5 | -0.05 | 0.5
02 | MAD | 0.5 | +0.25 | 1.75
03 | SANE | 0 | +0.15 | 4.1
Table 4: Summary of synthetic data models.

Figure 19 shows the GRMHD ground truth models and reconstructions from our
three synthetic datasets (Appendix E), for each of the five imaging or
modeling 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 $\mathcal{V}$ reconstructions
(bottom rows) show significant variation between the methods, and most methods
have difficulty reconstructing the ground-truth features in circular
polarization.

All methods struggle the most with model 01, the weakly polarized case
$(\langle|v|\rangle\approx 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 $\mathcal{V}$ brightness in
the image. In contrast, for the “high-polarization” model 03,
$(\langle|v|\rangle\approx 4\,\%)$, all methods are able to recover the
resolved distribution of $\mathcal{V}$ somewhat more accurately. Themis,
polsolve, and $m$-ring modeling all recover the correct dipolar Stokes
$\mathcal{V}$ structure of the image; DMC and DIFMAP produce mostly one sign
of $\mathcal{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 $\mathcal{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 $\mathcal{V}$ in the image and the
orientation of the dipolar structure across the ring. The limits on structural
complexity baked into the $m_{\mathcal{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 $\mathcal{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_{\rm
net}$. The recovered resolved circular polarization $\langle|v|\rangle$ 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 $\langle|v|\rangle$ within their $2\sigma$ error
bars, while the polsolve and DIFMAP single-image results are off by $\approx
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 $\mathcal{V}$ and low $\mathcal{V}$ S/N in
the data, even after self-calibration during imaging. Furthermore,
$\langle|v|\rangle$ 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 $\langle|v|\rangle$ in
both cases, while overproducing $\langle|v|\rangle$ 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 $\sim~{}10\,\%$ do occur. DMC is the only method that returns uncertainty
estimates on these gain ratios; the DMC results are nearly all within
$2\sigma$ 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 $\mathcal{V}$ imaging results for M87∗.

Figure 22:  M87∗ Stokes $\mathcal{V}$ imaging results using data reduced with
the CASA-based rPICARD pipeline (Janssen et al., 2019) and pre-calibrated with
the assumption $\mathcal{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.

##  Appendix F M87* Results From CASA-calibrated Data

Figure 22 presents Stokes $\mathcal{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 $\mathcal{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 $\mathcal{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 $\mathcal{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 $\mathcal{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.

##  Appendix G Biases in Resolved Circular Polarization Fraction

For an image composed of $N_{\text{res}}$ pixels or resolution elements, each
of which has angular area $\Delta A$, then Equation 5 can be expressed as a
sum

$\langle|v|\rangle=\frac{\Delta
A}{F}\sum_{i=1}^{N_{\text{res}}}\left|\frac{\mathcal{V}_{i}}{\mathcal{I}_{i}}\right|\mathcal{I}_{i},$
(G1)

where $F$ is the total integrated Stokes $\mathcal{I}$ flux density in the
image, and $i$ is a subscript that indicates a particular pixel or resolution
element. If $\mathcal{I}$ is everywhere positive in the image, then
$\left|\mathcal{V}_{i}/\mathcal{I}_{i}\right|\mathcal{I}_{i}=|\mathcal{V}_{i}|$
and we can simplify the expression for $\langle|v|\rangle$ to be

$\langle|v|\rangle=\frac{\Delta
A}{F}\sum_{i=1}^{N_{\text{res}}}\left|\mathcal{V}_{i}\right|.$ (G2)

If $\mathcal{V}_{i}$ is purely noise-like and distributed according to a
normal distribution with mean zero and variance $\sigma_{i}^{2}$, then the
probability distribution function of $|\mathcal{V}_{i}|$ is given by

$|\mathcal{V}_{i}|\sim\begin{cases}\sqrt{\frac{2}{\pi\sigma_{i}^{2}}}\exp\left(-\frac{\mathcal{V}_{i}^{2}}{2\sigma_{i}^{2}}\right),&|\mathcal{V}_{i}|\geq
0\\\ 0,&|\mathcal{V}_{i}|<0\end{cases}.$ (G3)

The mean of this distribution is $\mu_{i}=\sigma_{i}\sqrt{2/\pi}$, which is
nonzero. The mean of $\langle|v|\rangle$ will then be given by

$E\left[\langle|v|\rangle\right]=\frac{\Delta
A}{F}\sqrt{\frac{2}{\pi}}\sum_{i=1}^{N_{\text{res}}}\sigma_{i}.$ (G4)

If each pixel or resolution element has a similar noise variance
$\sigma_{i}^{2}$, then we can see that $E\left[\langle|v|\rangle\right]$ will
be proportional to $N_{\text{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.

##  Appendix 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
magneto-rotational 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 $\phi\equiv\Phi_{\text{BH}}(\dot{M}r_{g}^{2}c)^{-1/2}$, where
$\Phi_{\text{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 $\dot{M}$ set over the event horizon. The relative strength of the
magnetic flux near the horizon is determined by the dimensionless quantity
$\phi$ 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 $\phi\geq 50$ (e.g.
Tchekhovskoy et al., 2011; McKinney et al., 2012; Narayan et al., 2022), and
the Standard and Normal Evolution (SANE) state $\phi\approx 1-5$ (in Gaussian
units). SANE simulations were produced on a $288\times 128\times 128$ grid, a
fluid adiabatic index $\gamma=4/3$, and a boundary of
$r_{\text{out}}=50\,r_{g}$. MAD simulations adopt a $384\times 192\times 192$
resolution, $\gamma=13/9$, and a domain of $r_{\text{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_{\nu}\approx 0.5\,$Jy for our images at
$\nu=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_{\mathrm{low}}$ and
$R_{\mathrm{high}}$ are associated with the ion-to-electron temperature ratio,
where the electron temperature is defined as

$\displaystyle T_{e}=\frac{2m_{e}u_{\text{gas}}}{3\rho
k_{B}\left(2+R\right)}.$ (H1)

In order for the description in Equation (H1) 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

$\displaystyle
R=\frac{T_{i}}{T_{e}}=R_{\text{high}}\frac{\beta^{2}}{1+\beta^{2}}+R_{\text{low}}\frac{1}{1+\beta^{2}}$
(H2)

where the components describe the local plasma
$\beta=p_{\text{gas}}/p_{\text{mag}}$ with
$p_{\text{gas}}=(\gamma-1)u_{\text{gas}}$, $p_{\text{mag}}=B^{2}/8\pi$, and
further $u_{\text{gas}}$ being the total internal energy density. To regulate
the ratio in Equation (H2), our models explore $R_{\text{low}}=1,10$ and
$R_{\text{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 (H2) mimic
collisionless plasma properties and allows to associate the electron heating
with its magnetic properties (see Paper VIII, ).

##  Appendix 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.

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 24: Distribution
of image-averaged fractional linear polarization $\langle|m|\rangle$ with an
aligned (left panel) and anti-aligned (right panel) magnetic field respecting
all images in the M87∗ library blurred with a 20$\,\mu$as FWHM Gaussian beam.
Allowed inferred ranges of EHT image reconstructions are limited by the dashed
lines. Figure 25: Distribution of $\beta_{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\,\mu$as beam.
Allowed inferred ranges of EHT image reconstructions are limited by the dashed
lines. Figure 26: Distribution of $\beta_{2}$ phase 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\,\mu$as beam.
Allowed inferred ranges of EHT image reconstructions are limited by the dashed
lines.

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 $\angle\beta_{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 $\beta_{2}$ is also rotated counter-
clockwise. Similarly, a poloidal field pointed away from the observer will on
average rotate ticks $\beta_{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 $\angle\beta_{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 $\angle\beta_{2}$ when
the field polarity is reversed. Prograde models have a slight increase,
whereas retrograde models have a slight decrease in $\angle\beta_{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.

##  Appendix J More Details on GRMHD Scoring Results

###  J.1 Simultaneous Snapshot Scoring

For simultaneous scoring, the upper limit (3.7 %) for $\langle|v|\rangle$
reported in this paper does not provide any additional information to
constrain the GRMHD image parameter space. That is, all models that pass the
other polarimetric constraints already naturally pass this constraint. This is
perhaps not surprising, since the upper limit on $v_{\mathrm{net}}$ reported
by Goddi et al. (2021) and already considered in Paper VIII is more stringent
than that on $\langle|v|\rangle$ obtained in this work. Only 261/184796
snapshots ($0.14\,\%$) fail the $\langle|v|\rangle$ constraint but pass the
$v_{\mathrm{net}}$ constraint. Producing a low $v_{\rm net}$ but a high
$\langle|v|\rangle$ 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,
$\langle|v|\rangle$ 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 $\langle|v|\rangle$, 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.

Figure 27: Passing fraction of the $v_{\mathrm{net}}$ and $\langle|v|\rangle$
constraints for all the models in the simulation library. Aligned and anti-
aligned models are scored combined since the individual passing fractions are
similar.

To better understand how Stokes $\mathcal{V}$ alone affects our model scoring,
Figure 27 plots the fraction of snapshots of a given model that pass Stokes
$\mathcal{V}$ constraints without considering the others. We find that Stokes
$\mathcal{V}$ constraints are more likely to rule out SANE models than MADs,
since SANE models are more likely to overproduce $v_{\mathrm{net}}$. The
$\langle|v|\rangle$ 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_{\mathrm{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 $\angle\beta_{2\mathrm{,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
$\langle|v|\rangle$ constraint.

The inclusion of the $\langle|v|\rangle$ constraint removes small likelihoods
for MAD spin $0$,$-0.5$, $R_{\mathrm{low}}=1$ models but otherwise does not
qualitatively change the likelihoods. Note that the midpoint for the
$\langle|v|\rangle$ 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 $\langle|v|\rangle$ of 2 % from Themis (Table 2), we find that
qualitatively the joint scoring results remain largely the same, but some
$R_{\mathrm{low}}=1$ models now are included as those models have lower
$\langle|v|\rangle$. A lower limit on $\langle|v|\rangle$ at 1 % tends to
increase the likelihood of the MAD spin $+0.5$, $R_{\mathrm{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^{\circ}$ or
$163^{\circ}$ inclination, while Paper VIII included models at inclinations
$\pm 5^{\circ}$ from $17^{\circ}$ or $163^{\circ}$. Models from Paper VIII
scored using only $17^{\circ}$ or $163^{\circ}$ 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
$\beta_{2\mathrm{,LP}}$ modes. The $\beta_{2\mathrm{,LP}}$ modes appear to be
the most discrepant between the two library sets, and most linear and circular
polarization constraints otherwise remain consistent.

##  Appendix K Stokes $\mathcal{V}$ Origin in Simulations.

Figure 28:  Exploration of the effects of various CP production mechanisms on
$v_{\mathrm{net}}$ (top panels) and $\langle|v|\rangle$ (bottom panels). Each
column corresponds to a comparison between the full radiative transfer (RT)
simulation and the same calculation with one coefficient turned off:
$j_{\mathcal{V}}$=0 (intrinsic emission) for column 1, $\rho_{\mathcal{Q}}$=0
(Faraday conversion) for column 2, $\rho_{\mathcal{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.

Here, we present in more detail the test summarized in Subsection 5.5 to
determine the physical origin of Stokes $\mathcal{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_{\mathcal{V}}$), Faraday
conversion ($\rho_{\mathcal{Q}}$), and Faraday rotation
($\rho_{\mathcal{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_{\mathrm{net}}$ and $\langle|v|\rangle$
against the full radiative transfer solution and plot the results in Figure
28. Since significant cancellation occurs in the image plane for Stokes
$\mathcal{V}$, large differences to the image do not necessarily correspond to
large changes in $v_{\mathrm{net}}$. Therefore, $\langle|v|\rangle$ is
somewhat more informative, and robust to cancellations, although cancellations
within the Gaussian kernel of $20\,\mu\mathrm{as}$ still occur.

From Figure 28 we note 3 main results:

  1. 1.

Turning off intrinsic Stokes $\mathcal{V}$ emission ($j_{\mathcal{V}}$) does
not appear to greatly influence either $v_{\mathrm{net}}$ or
$\langle|v|\rangle$ (first column of Figure 28). While this is not true for
every snapshot in every model, for most of the passing models, we find that
intrinsic emission is not the dominant origin of Stokes $\mathcal{V}$.

  2. 2.

On the other hand, turning off Faraday conversion ($\rho_{\mathcal{Q}}=0$)
greatly suppresses $\langle|v|\rangle$ (second row, second column in Figure
28), suggesting that Faraday conversion is the dominant production mechanism
for Stokes $\mathcal{V}$ in passing models.

  3. 3.

Interestingly, turning off Faraday rotation increases $\langle|v|\rangle$ for
some of the passing models. This implies that scrambling via Faraday rotation
is 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 $\mathcal{V}$ generation
in our models, while intrinsic emission plays a sub-dominant, though non-
negligible, role.

##  Appendix L Alternative Plasma Models

Throughout this work, we have modeled a pure ion-electron 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 $\mathcal{V}$
images are indeed sensitive to these details, motivating future studies in
these areas.

###  L.1 Ionic vs. pair plasmas

Figure 29: Test of plasma content, where we ray-trace a single snapshot of the
MAD $a=+0.5$ $R_{\mathrm{high}}=80$ $R_{\mathrm{low}}=10$ aligned field model
with an increasing positron-to-electron ratio, denoted as $f$ in the top left
of each panel. The Stokes $\mathcal{V}$ structure clearly evolves as $f$
increases, but we do not observe a clear discriminant of plasma content.
Figure 30: A single GRMHD snapshot (MAD $a=+0.5$ $R_{\mathrm{high}}=80$
$R_{\mathrm{low}}=10$ aligned field) ray traced with three different electron
distribution functions. The physically motivated variable kappa model produces
both a Stokes $\mathcal{I}$ and Stokes $\mathcal{V}$ image very similar to
thermal model. However, a model with fixed $\kappa=5$ produces a much more
diffuse Stokes $\mathcal{I}$ and extremely little Stokes $\mathcal{V}$.

Unlike an ionic plasma (assumed throughout this paper), a pure pair plasma
produces no circular polarization via synchrotron emission and has no Faraday
rotation, but can still perform Faraday conversion (Jones & 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 $\mathcal{V}$ image. On EHT scales, the few recent studies on this
topic find that increasing the pair fraction can alter both resolved Stokes
$\mathcal{V}$ images at 230 GHz as well as the evolution as a function of
frequency (Anantua et al., 2020; Emami et al., 2021, 2023b).

In Figure 29, we test a single snapshot of the MAD $a=+0.5$
$R_{\mathrm{high}}=80$ $R_{\mathrm{low}}=10$ aligned field model, performing
polarized ray-tracing with an increasing positron-to-electron number density
ratio, denoted as $f$ in each panel. To obtain a given value of $f$, we inject
electron-positron 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)

$\displaystyle j_{I,Q,U}$ $\displaystyle\to(1+f)j_{I,Q,U},$ $\displaystyle
j_{V}$ $\displaystyle\to(1-f)j_{V},$ $\displaystyle\alpha_{I,Q,U}$
$\displaystyle\to(1+f)\alpha_{I,Q,U},$ $\displaystyle\alpha_{V}$
$\displaystyle\to(1-f)\alpha_{V},$ $\displaystyle\rho_{Q,U}$
$\displaystyle\to(1+f)\rho_{Q,U},$ $\displaystyle\rho_{V}$
$\displaystyle\to(1-f)\rho_{V},$ (L1)

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 $\mathcal{V}$ structures among the library, there is no
known signature in Stokes $\mathcal{V}$ that clearly indicates the presence of
pairs. Interestingly, we find that $\langle|v|\rangle$ 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, $\langle|v|\rangle$ 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 high-energy 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 $\kappa=5$ model, where $\kappa-1=4$ is the power law index of the
high-energy tail, and variable kappa model where $\kappa=\kappa(\sigma,\beta)$
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 $\kappa>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_{\mathrm{high}}=80$
$R_{\mathrm{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 $\mathcal{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 $\mathcal{V}$. However, we
report dramatic differences when switching to a constant kappa model with a
value of $\kappa=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 $\mathcal{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 $\kappa$ models
by definition.

Figure 31: Standard deviations of the distributions of $v_{\mathrm{net}}$ 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_{\mathrm{high}}$ is plotted on the $x$-axis of each
panel, and the color and shape of the markers encodes $R_{\mathrm{low}}$. In
particular, models with $R_{\mathrm{low}}=10$ are more variable than those
with $R_{\mathrm{low}}=1$, and a similar but weaker trend occurs with
$R_{\mathrm{high}}$.

##  Appendix M Stokes $\mathcal{V}$ Variability in Passing Models

In Figure 31, we consider each model’s distribution of $v_{\mathrm{net}}$ over
time and plot its standard deviation, which we denote as
$\sigma_{v_{\mathrm{net}}}$. We average together the standard deviations for
each magnetic field polarity, since their distributions can be disjoint. We
find that models with larger $R_{\mathrm{low}}$ are more variable in
$v_{\mathrm{net}}$, with a similar but weaker trend in $R_{\mathrm{high}}$.
Although not shown, we find that the variability of $\langle|v|\rangle$ is
qualitatively similar. Some models have $\sigma_{v_{\mathrm{net}}}$ 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
$\mathcal{V}$ on EHT baselines and hence more robust image reconstructions of
the source.

## References

  * Anantua et al. (2020) Anantua, R., Emami, R., Loeb, A., & Chael, A. 2020, ApJ, 896, 30, doi: 10.3847/1538-4357/ab9103
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Ball et al. (2016) Ball, D., Özel, F., Psaltis, D., & Chan, C.-k. 2016, ApJ, 826, 77, doi: 10.3847/0004-637X/826/1/77
  * Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1974, Ap&SS, 28, 45, doi: 10.1007/BF00642237
  * Blackburn et al. (2020) Blackburn, L., Pesce, D. W., Johnson, M. D., et al. 2020, ApJ, 894, 31, doi: 10.3847/1538-4357/ab8469
  * Blackburn et al. (2019) Blackburn, L., Chan, C.-k., Crew, G. B., et al. 2019, ApJ, 882, 23, doi: 10.3847/1538-4357/ab328d
  * Bower et al. (1999) Bower, G. C., Falcke, H., & Backer, D. C. 1999, ApJ, 523, L29, doi: 10.1086/312246
  * Bower et al. (2002) Bower, G. C., Falcke, H., & Mellon, R. R. 2002, ApJ, 578, L103, doi: 10.1086/344607
  * Bower et al. (2018) Bower, G. C., Broderick, A., Dexter, J., et al. 2018, ApJ, 868, 101, doi: 10.3847/1538-4357/aae983
  * Broderick & Pesce (2020) Broderick, A. E., & Pesce, D. W. 2020, ApJ, 904, 126, doi: 10.3847/1538-4357/abbd9d
  * Broderick et al. (2020a) Broderick, A. E., Pesce, D. W., Tiede, P., Pu, H.-Y., & Gold, R. 2020a, ApJ, 898, 9, doi: 10.3847/1538-4357/ab9c1f
  * Broderick et al. (2020b) Broderick, A. E., Gold, R., Karami, M., et al. 2020b, ApJ, 897, 139, doi: 10.3847/1538-4357/ab91a4
  * Chael et al. (2018) Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23, doi: 10.3847/1538-4357/aab6a8
  * Chael et al. (2016) Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11, doi: 10.3847/0004-637X/829/1/11
  * Chatterjee et al. (2020) Chatterjee, K., Younsi, Z., Liska, M., et al. 2020, MNRAS, 499, 362, doi: 10.1093/mnras/staa2718
  * Contopoulos & Kazanas (1998) Contopoulos, I., & Kazanas, D. 1998, ApJ, 508, 859, doi: 10.1086/306426
  * Contopoulos et al. (2022) Contopoulos, I., Myserlis, I., Kazanas, D., & Nathanail, A. 2022, Galaxies, 10, 80, doi: 10.3390/galaxies10040080
  * Decker & Krimigis (2003) Decker, R. B., & Krimigis, S. M. 2003, Advances in Space Research, 32, 597, doi: 10.1016/S0273-1177(03)00356-9
  * Dexter (2016) Dexter, J. 2016, MNRAS, 462, 115, doi: 10.1093/mnras/stw1526
  * Dihingia et al. (2023) Dihingia, I. K., Mizuno, Y., Fromm, C. M., & Rezzolla, L. 2023, MNRAS, 518, 405, doi: 10.1093/mnras/stac3165
  * Emami et al. (2021) Emami, R., Anantua, R., Chael, A. A., & Loeb, A. 2021, ApJ, 923, 272, doi: 10.3847/1538-4357/ac2950
  * Emami et al. (2023a) Emami, R., Ricarte, A., Wong, G. N., et al. 2023a, ApJ, 950, 38, doi: 10.3847/1538-4357/acc8cd
  * Emami et al. (2023b) Emami, R., Anantua, R., Ricarte, A., et al. 2023b, Galaxies, 11, 11, doi: 10.3390/galaxies11010011
  * Event Horizon Telescope Collaboration (2023) Event Horizon Telescope Collaboration. 2023, EHT M87 Polarized Data, CyVerse Data Commons, doi: 10.25739/q46m-m857
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration, et al. 2019a, ApJ, 875, L1, doi: 10.3847/2041-8213/ab0ec7
  * Event Horizon Telescope Collaboration et al. (2019b) —. 2019b, ApJ, 875, L2, doi: 10.3847/2041-8213/ab0c96
  * Event Horizon Telescope Collaboration et al. (2019c) —. 2019c, ApJ, 875, L3, doi: 10.3847/2041-8213/ab0c57
  * Event Horizon Telescope Collaboration et al. (2019d) —. 2019d, ApJ, 875, L4, doi: 10.3847/2041-8213/ab0e85
  * Event Horizon Telescope Collaboration et al. (2019e) —. 2019e, ApJ, 875, L5, doi: 10.3847/2041-8213/ab0f43
  * Event Horizon Telescope Collaboration et al. (2019f) —. 2019f, ApJ, 875, L6, doi: 10.3847/2041-8213/ab1141
  * Event Horizon Telescope Collaboration et al. (2021a) —. 2021a, ApJ, 910, L12, doi: 10.3847/2041-8213/abe71d
  * Event Horizon Telescope Collaboration et al. (2021b) —. 2021b, ApJ, 910, L13, doi: 10.3847/2041-8213/abe4de
  * Event Horizon Telescope Collaboration et al. (2022a) —. 2022a, ApJ, 930, L12, doi: 10.3847/2041-8213/ac6674
  * Event Horizon Telescope Collaboration et al. (2022b) —. 2022b, ApJ, 930, L14, doi: 10.3847/2041-8213/ac6429
  * Event Horizon Telescope Collaboration et al. (2022c) —. 2022c, ApJ, 930, L15, doi: 10.3847/2041-8213/ac6736
  * Event Horizon Telescope Collaboration et al. (2022d) —. 2022d, ApJ, 930, L16, doi: 10.3847/2041-8213/ac6672
  * Fishbone & Moncrief (1976) Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962, doi: 10.1086/154565
  * Fromm et al. (2022) Fromm, C. M., Cruz-Osorio, A., Mizuno, Y., et al. 2022, A&A, 660, A107, doi: 10.1051/0004-6361/202142295
  * Gabuzda et al. (2008) Gabuzda, D. C., Vitrishchak, V. M., Mahmud, M., & O’Sullivan, S. P. 2008, MNRAS, 384, 1003, doi: 10.1111/j.1365-2966.2007.12773.x
  * Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444, doi: 10.1086/374594
  * Goddi et al. (2019) Goddi, C., Martí-Vidal, I., Messias, H., et al. 2019, PASP, 131, 075003, doi: 10.1088/1538-3873/ab136a
  * Goddi et al. (2021) —. 2021, ApJ, 910, L14, doi: 10.3847/2041-8213/abee6a
  * Greisen (2003) Greisen , E. W. 2003, in Astrophysics and Space Science Library, Vol. 285, Information Handling in Astronomy - Historical Vistas, ed. A. Heck, 109, doi: 10.1007/0-306-48080-8_7
  * Hamaker et al. (1996) Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Homan & Lister (2006) Homan, D. C., & Lister, M. L. 2006, AJ, 131, 1262, doi: 10.1086/500256
  * Homan & Wardle (1999) Homan, D. C., & Wardle, J. F. C. 1999, AJ, 118, 1942, doi: 10.1086/301108
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Igumenshchev et al. (2003) Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042, doi: 10.1086/375769
  * Issaoun et al. (2022) Issaoun, S., Wielgus, M., Jorstad, S., et al. 2022, ApJ, 934, 145, doi: 10.3847/1538-4357/ac7a40
  * Janssen et al. (2022) Janssen, M., Radcliffe, J. F., & Wagner, J. 2022, Universe, 8, 527, doi: 10.3390/universe8100527
  * Janssen et al. (2019) Janssen, M., Goddi, C., van Bemmel, I. M., et al. 2019, A&A, 626, A75, doi: 10.1051/0004-6361/201935181
  * Janssen et al. (2021) Janssen, M., Falcke, H., Kadler, M., et al. 2021, Nature Astronomy, 5, 1017, doi: 10.1038/s41550-021-01417-w
  * Jiang et al. (2023) Jiang, H.-X., Mizuno, Y., Fromm, C. M., & Nathanail, A. 2023, MNRAS, doi: 10.1093/mnras/stad1106
  * Johnson et al. (2020) Johnson, M. D., Lupsasca, A., Strominger, A., et al. 2020, Science Advances, 6, eaaz1310, doi: 10.1126/sciadv.aaz1310
  * Jones (1988) Jones, T. W. 1988, ApJ, 332, 678, doi: 10.1086/166685
  * Jones & O’Dell (1977) Jones, T. W., & O’Dell, S. L. 1977, ApJ, 214, 522, doi: 10.1086/155278
  * Jorstad et al. (2023) Jorstad, S., Wielgus, M., Lico, R., et al. 2023, ApJ, 943, 170, doi: 10.3847/1538-4357/acaea8
  * Kettenis et al. (2006) Kettenis, M., van Langevelde, H. J., Reynolds, C., & Cotton, B. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 351, Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, 497
  * Kim et al. (2020) Kim, J.-Y., Krichbaum, T. P., Broderick, A. E., et al. 2020, A&A, 640, A69, doi: 10.1051/0004-6361/202037493
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt (IOS Press), 87 – 90
  * Kunz et al. (2016) Kunz, M. W., Stone, J. M., & Quataert, E. 2016, Phys. Rev. Lett., 117, 235101, doi: 10.1103/PhysRevLett.117.235101
  * Leung et al. (2011) Leung, P. K., Gammie, C. F., & Noble, S. C. 2011, ApJ, 737, 21, doi: 10.1088/0004-637X/737/1/21
  * MacDonald & Marscher (2018) MacDonald, N. R., & Marscher, A. P. 2018, ApJ, 862, 58, doi: 10.3847/1538-4357/aacc62
  * Mao et al. (2017) Mao, S. A., Dexter, J., & Quataert, E. 2017, MNRAS, 466, 4307, doi: 10.1093/mnras/stw3324
  * Marszewski et al. (2021) Marszewski, A., Prather, B. S., Joshi, A. V., Pandya, A., & Gammie, C. F. 2021, ApJ, 921, 17, doi: 10.3847/1538-4357/ac1b28
  * Martí-Vidal & Marcaide (2008) Martí-Vidal, I., & Marcaide, J. M. 2008, A&A, 480, 289, doi: 10.1051/0004-6361:20078690
  * Martí-Vidal et al. (2021) Martí-Vidal, I., Mus, A., Janssen, M., de Vicente, P., & González, J. 2021, A&A, 646, A52, doi: 10.1051/0004-6361/202039527
  * Martí-Vidal et al. (2016) Martí-Vidal, I., Roy, A., Conway, J., & Zensus, A. J. 2016, A&A, 587, A143, doi: 10.1051/0004-6361/201526063
  * Matthews et al. (2018) Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018, PASP, 130, 015002, doi: 10.1088/1538-3873/aa9c3d
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083, doi: 10.1111/j.1365-2966.2012.21074.x
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Mizuno et al. (2021) Mizuno, Y., Fromm, C. M., Younsi, Z., et al. 2021, MNRAS, 506, 741, doi: 10.1093/mnras/stab1753
  * Mościbrodzka et al. (2016) Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38, doi: 10.1051/0004-6361/201526630
  * Mościbrodzka & Gammie (2018) Mościbrodzka, M., & Gammie, C. F. 2018, MNRAS, 475, 43, doi: 10.1093/mnras/stx3162
  * Mościbrodzka et al. (2021) Mościbrodzka, M., Janiuk, A., & De Laurentis, M. 2021, MNRAS, 508, 4282, doi: 10.1093/mnras/stab2790
  * Myserlis et al. (2018) Myserlis, I., Angelakis, E., Kraus, A., et al. 2018, A&A, 609, A68, doi: 10.1051/0004-6361/201630301
  * Narayan et al. (2022) Narayan, R., Chael, A., Chatterjee, K., Ricarte, A., & Curd, B. 2022, MNRAS, 511, 3795, doi: 10.1093/mnras/stac285
  * Narayan et al. (2003) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  * Narayan et al. (2012) Narayan, R., Sądowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241, doi: 10.1111/j.1365-2966.2012.22002.x
  * Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626, doi: 10.1086/500349
  * Noble et al. (2007) Noble, S. C., Leung, P. K., Gammie, C. F., & Book, L. G. 2007, Classical and Quantum Gravity, 24, S259, doi: 10.1088/0264-9381/24/12/S17
  * Özel et al. (2000) Özel, F., Psaltis, D., & Narayan, R. 2000, ApJ, 541, 234, doi: 10.1086/309396
  * Palumbo et al. (2020) Palumbo, D. C. M., Wong, G. N., & Prather, B. S. 2020, ApJ, 894, 156, doi: 10.3847/1538-4357/ab86ac
  * Pesce (2021) Pesce, D. W. 2021, AJ, 161, 178, doi: 10.3847/1538-3881/abe3f8
  * Pierrard & Lazar (2010) Pierrard, V., & Lazar, M. 2010, Sol. Phys., 267, 153, doi: 10.1007/s11207-010-9640-2
  * Pordes et al. (2007) Pordes, R., Petravick, D., Kramer, B., et al. 2007, in 78, Vol. 78, J. Phys. Conf. Ser., 012057, doi: 10.1088/1742-6596/78/1/012057
  * Porth et al. (2019) Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26, doi: 10.3847/1538-4365/ab29fd
  * Qiu et al. (2023) Qiu, R., Ricarte, A., Narayan, R., et al. 2023, MNRAS, doi: 10.1093/mnras/stad466
  * Ricarte et al. (2023) Ricarte, A., Gammie, C., Narayan, R., & Prather, B. S. 2023, MNRAS, 519, 4203, doi: 10.1093/mnras/stac3796
  * Ricarte et al. (2021) Ricarte, A., Qiu, R., & Narayan, R. 2021, MNRAS, 505, 523, doi: 10.1093/mnras/stab1289
  * Roelofs et al. (2023) Roelofs, F., Johnson, M. D., Chael, A., et al. 2023, ApJ, 957, L21, doi: 10.3847/2041-8213/acff6f
  * Sfiligoi et al. (2009) Sfiligoi, I., Bradley, D. C., Holzman, B., et al. 2009, in 2, Vol. 2, 2009 WRI World Congress on Computer Science and Information Engineering, 428–432, doi: 10.1109/CSIE.2009.950
  * Shepherd (2011) Shepherd, M. 2011, Difmap: Synthesis Imaging of Visibility Data, Astrophysics Source Code Library.  http://ascl.net/1103.001
  * Sądowski et al. (2013) Sądowski, A., Narayan, R., Penna, R., & Zhu, Y. 2013, MNRAS, 436, 3856, doi: 10.1093/mnras/stt1881
  * Smirnov (2011) Smirnov, O. M. 2011, A&A, 527, A106, doi: 10.1051/0004-6361/201016082
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa27810.48550/arXiv.1904.02180
  * Steel et al. (2019) Steel, S., Wielgus, M., Blackburn, L., Issaoun, S., & Johnson, M. 2019, EHT Memo Series, 2019-CE-03
  * Tange (2011) Tange, O. 2011, ;login: The USENIX Magazine, 36, 42.  http://www.gnu.org/s/parallel
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  * The Astropy Collaboration et al. (2013) The Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Thompson et al. (2017) Thompson, A. R., Moran, J. M., & Swenson, Jr., G. W. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd Edition (Springer International Publishing), doi: 10.1007/978-3-319-44431-4
  * Thum et al. (2018) Thum, C., Agudo, I., Molina, S. N., et al. 2018, MNRAS, 473, 2506, doi: 10.1093/mnras/stx2436
  * Tsunetoe et al. (2022) Tsunetoe, Y., Mineshige, S., Kawashima, T., et al. 2022, ApJ, 931, 25, doi: 10.3847/1538-4357/ac66dd
  * Tsunetoe et al. (2020) Tsunetoe, Y., Mineshige, S., Ohsuga, K., Kawashima, T., & Akiyama, K. 2020, PASJ, 72, 32, doi: 10.1093/pasj/psaa008
  * Tsunetoe et al. (2021) —. 2021, PASJ, 73, 912, doi: 10.1093/pasj/psab054
  * Vasyliunas (1968) Vasyliunas, V. M. 1968, J. Geophys. Res., 73, 2839, doi: 10.1029/JA073i009p02839
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vitrishchak et al. (2008) Vitrishchak, V. M., Gabuzda, D. C., Algaba, J. C., et al. 2008, MNRAS, 391, 124, doi: 10.1111/j.1365-2966.2008.13919.x
  * Wardle & Homan (2003) Wardle, J. F. C., & Homan, D. C. 2003, Ap&SS, 288, 143, doi: 10.1023/B:ASTR.0000005001.80514.0c
  * Wardle et al. (1998) Wardle, J. F. C., Homan, D. C., Ojha, R., & Roberts, D. H. 1998, Nature, 395, 457, doi: 10.1038/26675
  * Wielgus et al. (2020) Wielgus, M., Akiyama, K., Blackburn, L., et al. 2020, ApJ, 901, 67, doi: 10.3847/1538-4357/abac0d
  * Wielgus et al. (2022a) Wielgus, M., Moscibrodzka, M., Vos, J., et al. 2022a, A&A, 665, L6, doi: 10.1051/0004-6361/202244493
  * Wielgus et al. (2022b) Wielgus, M., Marchili, N., Martí-Vidal, I., et al. 2022b, ApJ, 930, L19, doi: 10.3847/2041-8213/ac6428
  * Wong & Gammie (2022) Wong, G. N., & Gammie, C. F. 2022, ApJ, 937, 60, doi: 10.3847/1538-4357/ac854d
  * Wong et al. (2022) Wong, G. N., Prather, B. S., Dhruv, V., et al. 2022, ApJS, 259, 64, doi: 10.3847/1538-4365/ac582e
  * Yoon et al. (2020) Yoon, D., Chatterjee, K., Markoff, S. B., et al. 2020, MNRAS, 499, 3178, doi: 10.1093/mnras/staa3031