First Sagittarius A* Event Horizon Telescope Results. VIII. Physical Interpretation of the Polarized Ring

Synchrotron emission from the plasma near supermassive black holes (BHs) provides a crucial source of insight into the physical processes that drive accretion and outflow in galactic cores. It is intrinsically polarized, and both linear polarization and circular polarization provide information about the emitting plasma's density, temperature, composition, and magnetic field. In the rest frame of the emitting fluid, the linear polarization direction is orthogonal to the local magnetic fields, so images of linear polarization capture the projected magnetic field structure perpendicular to the line of sight. Any magnetized plasma along the line of sight imparts additional polarimetric effects via Faraday rotation, which rotates the plane of linear polarization with a λ² dependence, where λ is the observing wavelength, and Faraday conversion, which exchanges linear and circular polarization states. Finally, for emission near a BH, the polarization is subject to achromatic rotation from propagation in a curved spacetime.

Recently, the Event Horizon Telescope (EHT) Collaboration published images of the supermassive BH at the Galactic center, Sagittarius A* (Sgr A*; Event Horizon Telescope Collaboration et al. 2022a, 2022b, 2022c, 2022d, 2022e, 2022f, hereafter Papers I–VI). These images revealed a bright emission ring encircling a central brightness depression (the "apparent shadow"), consistent with the expected appearance of a Kerr BH with a mass M ≈ 4 × 10⁶ M_⊙ that is only accreting a trickle of material relative to that captured at the Bondi radius in a radiatively inefficient manner (e.g., Hilbert 1917; Bardeen 1973; Luminet 1979; Jaroszynski & Kurpiewski 1997; Falcke et al. 2000). Comparisons of the EHT measurements with numerical simulations provide estimates of the mass accretion rate $\dot{M}\sim {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\sim {10}^{-3}{\dot{M}}_{{\rm{B}}}$ and a luminosity that is L ≲ 10³⁶ erg s⁻¹ ∼ 10⁻⁹ L_Edd (see, e.g., Paper V, and references therein). Here ${\dot{M}}_{{\rm{B}}}$ is the Bondi mass accretion rate and L_Edd ≡ 4π GMc m_p /σ_T is the Eddington luminosity, with G, c, m_p , and σ_T being the gravitational constant, speed of light, proton mass, and Thomson cross section, respectively. Previously, measurements of linearly polarized emission near Sgr A* gave strong evidence for this low accretion state (e.g., Agol 2000; Quataert & Gruzinov 2000). In addition, the emission ring morphology, including the lack of a pronounced brightness asymmetry in EHT images, favors a viewing angle in Sgr A* that is at a low to moderate inclination (≲50°) relative to the angular momentum of the inner accretion flow (see, e.g., Figure 9 in Paper V).

Event Horizon Telescope Collaboration et al. (2024, hereafter Paper VII) reports the first polarized images of Sgr A*, using EHT observations at 230 GHz taken in 2017. These images show a prominent spiral polarization pattern in the emission ring that is temporally stable, strongly linearly polarized (≈25%), and dominated by azimuthally symmetric structure. Both the image-averaged polarization fraction (m_net ∼ 5%) and the resolved polarization fraction (〈∣m∣〉 ≈ 25%) are significantly higher in Sgr A* than in the EHT's observations of M87* (Event Horizon Telescope Collaboration et al. 2021a, hereafter M87* Paper VII). In M87*, this polarization pattern was explained by coherent and dynamically important magnetic fields, depolarized by Faraday effects (Event Horizon Telescope Collaboration et al. 2021b, hereafter M87* Paper VIII).

In this paper, we provide the theoretical modeling and interpretation to accompany Paper VII. In Section 2, we summarize the new polarimetric observational constraints on Sgr A*. In Section 3, we provide general arguments about what these constraints imply for Sgr A* through comparison with three simple models: one-zone physical models to evaluate the plasma properties, geometrical ring models to evaluate the degree of coherence in the polarized image, and semianalytic emission models to evaluate the interplay between spacetime and emission parameters in determining polarized image structure. In Section 4, we describe a large library of general relativistic magnetohydrodynamic (GRMHD) simulations for Sgr A*. In Section 5, we evaluate which of these GRMHD models are compatible with the observational constraints. In Section 6, we summarize our findings and describe the prospects for improved constraints from future observations of Sgr A*.

In Paper VII, static polarimetric images are constructed from the Sgr A* EHT data taken on 2017 April 6th and 7 between 226.1 and 230.1 GHz (see Section 2 of Paper VII for more details). For theoretical interpretation, we adopt eight observational constraints derived from images generated by the THEMIS and the m-ring reconstruction methods (note that "m" is the azimuthal/angular mode number here, not polarization fraction; see Johnson et al. 2020). Of the four methods included in Paper VII, these are the only methods that provide Bayesian posteriors, from which we compute 90% confidence intervals. These methods make drastically different assumptions and, in a sense, bracket the possible spatial and temporal variability. In brief, the m-ring method fits a ring model to each snapshot independently, but the allowed spatial variability is very limited by construction (m ≤ 2 for total intensity, m ≤ 3 for linear polarization, and m ≤ 2 for circular polarization). In contrast, Themis attempts to optimize a single static image most consistent with the full data over time, with a noise budget attributed to time variability. Despite the vast differences between these models, they recover key image quantities with similar accuracy in synthetic data tests and arrive at mostly consistent observables (Paper VII).

Throughout this work, the large and time-variable rotation measure (RM) of Sgr A* poses a significant systematic uncertainty. Defined as RM ≡ Δχ/Δλ², where χ is the electric vector position angle (EVPA), the RM of Sgr A* may originate from Faraday rotation internal to the emitting region, an external screen, changes in the plasma probed as a function of optical depth, or a combination of these effects. Examining the polarized light curves for the same 2 days as our EHT observations, Wielgus et al. (2024) arrive at $\langle \mathrm{RM}\rangle =-{4.65}_{-1.18}^{+1.25}\times {10}^{5}\ \mathrm{rad}\,{{\rm{m}}}^{-2}$. We reserve a lengthy discussion of the RM of Sgr A* in both observations and theory for Appendix C. In summary, the fraction of the RM that can be attributed to an external Faraday screen is currently unresolved. Thus, throughout this work we consider the recovered image statistics both with and without RM derotation. Derotating the image corresponds to an interpretation where the time-averaged RM is attributed to a relatively stable external Faraday screen, separate from our models, which can be corrected for. Refraining from doing so corresponds to an interpretation in which all of the RM is generated internally, within our models. Our GRMHD simulations can reproduce the intraday variability of the RM, but not its stability of sign (see Appendix C).

For each of these methods, eight observational constraints explored in this paper are computed, listed in Table 1. To generate these ranges, a large quantity of images consistent with the data were generated from each method's posterior distribution. We computed the relevant observables for each of these images and then inferred 90% confidence regions. The m-ring method does not provide independent values of v_net, which is fixed to the mean ALMA-inferred value for circular polarization analysis (see Paper VII). When combining the two methods for theoretical interpretation, we adopt the minimum and maximum of the union of both 90% confidence regions (see Figure 10 in Paper VII for a visualization).

Table 2. Summary of the Sgr A* GRMHD Simulation Library Used in This Work

Setup	GRMHD	GRRT	a*	Mode	Γad	tfinal	rout	Resolution
Torus	KHARMA	ipole	0, ±0.5, ±0.94	MAD/SANE	$\tfrac{4}{3}/\tfrac{4}{3}$	50,000	1000	288 × 128 × 128
Torus	BHAC	RAPTOR	0, ±0.5, ±0.94	MAD/SANE	$\tfrac{4}{3}/\tfrac{4}{3}$	30,000	3333	512 × 192 × 192
Torus	H-AMR	ipole	0, ±0.5, ±0.94	MAD/SANE	$\tfrac{13}{9}/\tfrac{5}{3}$	35,000	1000/200	348/240 × 192 × 192

Note. The last column is N₁ × N₂ × N₃, with coordinate x₁ monotonic in radius, x₂ monotonic in colatitude θ, and x₃ proportional to longitude ϕ. Times are given in units of t_g and radii in units of r_g . Different settings may be adopted for MAD models compared to SANE ones, as denoted by a /.

Download table as:  ASCIITypeset image

The quantities m_net and v_net correspond to the net linear and circular polarization that would be inferred from a spatially unresolved measurement for the time-averaged image. These are given by

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

$$\begin{eqnarray}&&{v}_{\mathrm{net}}=\displaystyle \frac{{\sum }_{i}{{ \mathcal V }}_{i}}{{\sum }_{i}{{ \mathcal I }}_{i}},\end{eqnarray} \tag{ 2 }$$

where ∑_i denotes a summation over each pixel i. For the time-resolved light curves, which are distinct from the values inferred from our static image reconstructions, Wielgus et al. (2022b, 2024) find 2.6% < m_net < 11% and −2.1% < v_net < − 0.7%, respectively, where we quote the central 90% of the values observed during the same 2 days of observation. Interestingly, we find that the m-ring method arrives at much lower values of m_net than Themis, which may be attributable to temporal cancellations of fluctuating EVPA patterns.

The remainder of our constraints are structural quantities, beginning with 〈∣m∣〉 and 〈∣v∣〉, the image-averaged linear and circular polarization fraction. These are given by

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

$$\begin{eqnarray}&&\langle | v| \rangle =\displaystyle \frac{{\sum }_{i}\left|{{ \mathcal V }}_{i}/{{ \mathcal I }}_{i}\right|{{ \mathcal I }}_{i}}{{\sum }_{i}{{ \mathcal I }}_{i}}.\end{eqnarray} \tag{ 4 }$$

Note that these quantities depend on the effective resolution of our images. Throughout this work we quote values from our simulations corresponding to 20 μas resolution to mimic EHT resolution. We treat the resolved circular polarization fraction 〈∣v∣〉 as an upper limit, and thus the combined range extends to 0 in Table 1. This is due to the fact that the circularly polarized images presented in Paper VII show structural differences that we attribute to noise (see also Event Horizon Telescope Collaboration et al. 2023, hereafter M87* Paper IX). Because of the absolute magnitude inherent to the definition of this quantity, it is biased high when the signal-to-noise ratio is too low.

Complex β_m modes correspond to Fourier decompositions of the linear polarization structure, where m refers to the number of times that an EVPA tick rotates with azimuth (Palumbo et al. 2020). These coefficients are defined by

$$\begin{eqnarray}&&{\beta }_{m}=\displaystyle \frac{1}{{I}_{\mathrm{tot}}}{\int }_{0}^{\infty }{\int }_{0}^{2\pi }P(\rho ,\varphi )\,{e}^{-{im}\varphi }\ \rho d\varphi d\rho ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{I}_{\mathrm{tot}}={\int }_{0}^{\infty }{\int }_{0}^{2\pi }I(\rho ,\varphi )\ \rho d\varphi d\rho \end{eqnarray} \tag{ 6 }$$

where ρ and φ correspond to polar coordinates in the image and P = Q + iU. The rotationally invariant mode, β₂, has natural connections to what we believe are azimuthally symmetric disk/jet structures, in particular the magnetic field geometry. Its amplitude encodes the strength of this mode, while its phase encodes the pitch angle and handedness of EVPA ticks. We observe ∠β₂ closer to ±180° than 0°, which corresponds to tick patterns that are more toroidal than radial.

When considering observational constraints without RM derotation, we simply adopt the range of ∠β₂ as observed on the sky. When considering observational constraints with RM derotation, we derotate ∠β₂ assuming that there is an external Faraday screen between us and the emitting region that we can characterize by the mean RM over time. Since ∠β₂ depends on twice the EVPA, we therefore add −2〈RM〉λ² to ∠β₂, where 〈RM〉 is the mean RM observed on April 6 and 7. Therefore, the range on ∠β₂ had been significantly shifted by the Faraday screen by $2\mathrm{RM}{\lambda }^{2}=-{92.0}_{-23.4}^{+24.7}$ deg. Applying this derotation both shifts and broadens the constraint.

Mean images from the posterior distributions generated by each method are plotted in Figure 1. Two sets of linearly polarized images are shown, corresponding to images without and with derotation, respectively. Note that derotation reverses the handedness of the polarization spiral, which has important implications for the flow structure. In the first two rows, total intensity is shown in gray scale, with contours drawn at 25%, 50%, and 75% of the peak brightness. These same contours are repeated in the bottom row. In the top and middle rows, the colored ticks encode linear polarization, where the length scales with the total linearly polarized intensity and the color scales with the fractional polarization. The dashed white contours plot the linearly polarized intensity rather than the total intensity.

Zoom In Zoom Out Reset image size

Finally, we also compute the simplest nonrotationally symmetric mode, β₁, as a probe of polarization asymmetry. Again, ∣β₁∣ encodes the strength of this mode, and we use ∣β₂∣/∣β₁∣ as a probe of rotational symmetry. Since there is no clear axis (such as the spin axis) to define ∠β₁ = 0°, we do not study ∠β₁. We also refrain from computing higher-order β_m modes, which are more likely to be sensitive to smaller-scale noise fluctuations.

As discussed in the previous section, the linearly polarized image of Sgr A* exhibits three salient features:

• 1.  
  It has a large resolved polarization fraction of 24%–28%, with a peak of ∼40%, much higher than M87*.
• 2.  
  The linear polarization structure is highly ordered.
• 3.  
  The ordered structure exhibits a high degree of rotational symmetry, which appears to spiral inward with counterclockwise handedness after derotating by the apparent RM, or clockwise without derotating.

Before exploring more physically complete GRMHD models, we demonstrate that each of these features can be understood in the context of simple analytic models.

3.1. One-zone Modeling

We use the basic assumptions described in Paper V that Sgr A* is an accreting BH with extremely small Eddington ratio and follow M87* Paper VIII to include polarimetry. This polarized one-zone model validates the more complicated numerical models shown later in this paper and offers a natural explanation for the high polarization fraction of Sgr A* relative to M87*.

We model the accretion flow around Sgr A* as a uniform sphere of plasma with radius r = 5 r_g , where r_g = GM/c², comparable to the observed size of Sgr A* at 230 GHz (Papers III and IV), with uniform magnetic field oriented at a fiducial 60° inclination relative to the line of sight. The outcomes of our one-zone model depend only weakly on the field orientation. Note that the plasma velocity and the gravitational redshift are neglected.

In Paper V, we assumed that the plasma is optically thin, the ion−electron temperature ratio is 3, the ions are subvirial by a factor of 3, and plasma β ≡ P_gas/P_mag = 1. Adopting the observational flux constraint F_ν = 2.4 Jy (Wielgus et al. 2022a), we obtained the self-consistent solution n_e ≃ 10⁶ cm⁻³ and B ≃ 29 G. Using this solution, we can estimate the strength of the Faraday rotation at 230 GHz with the optical depth to Faraday rotation ${\tau }_{{\rho }_{V}}$:

$$\begin{eqnarray}&&{\tau }_{{\rho }_{V}}\approx r\times {\rho }_{V}\simeq 0.98\left(\displaystyle \frac{r}{5\ {r}_{g}}\right),\end{eqnarray} \tag{ 7 }$$

where ρ_V is the Faraday rotation coefficient (e.g., Jones & Hardee 1979). In contrast, similar modeling arrived at ${\tau }_{{\rho }_{V}}\sim 5.2(r/5\ {r}_{g})$ for M87* (M87* Paper VIII). The value inferred for Sgr A* suggests that the internal Faraday rotation may not be negligible (see also Wielgus et al. 2024), but it also may not necessarily lead to substantial depolarization.

By including optical depth effects and using the Dexter (2016) polarized synchrotron emission and transfer coefficients, we relax some assumptions, such as ion−electron temperature ratios and virial factor, and plot the allowed parameter space as in M87* Paper VIII. Specifically,

• 1.  
  we relax the flux constraint to 2 Jy < F_ν < 3 Jy to include the effect of variability; and
• 2.  
  we require the same assumption that Sgr A* is optically thin, i.e., τ < 1.

The above requirements are marked by the gray and green regions in Figure 2, respectively. The magnetic field strengths are shown as red dotted contour lines, and the different panels assume different plasma β. In blue, we plot the contour corresponding to ${\tau }_{{\rho }_{V}}\gt 2\pi$, beyond which internal Faraday depolarization becomes increasingly important. Unlike for M87* (see Figure 2 of M87* Paper VIII), we find that the regions where the total flux and optically thin constraints are satisfied only occur in Faraday thin regions of parameter space. We note that this is compatible with multifrequency RM measurements that suggest ${\tau }_{{\rho }_{V}}\sim 1$ (Wielgus et al. 2024). Again, this is enough to noticeably rotate the EVPA pattern, but not enough to cause substantial depolarization.

Zoom In Zoom Out Reset image size

In summary, the total flux and optical depth constraints of Sgr A* naturally require small Faraday depths, which explains the large inferred values of 〈∣m∣〉.

3.2. Ordered Polarization: Ordered Fields

Because beam depolarization can only decrease the observed polarization fraction, measurements of the linear polarization at varying angular scales provide information about the degree of order in the underlying polarization. A priori, it could be possible that the the underlying magnetic field is significantly tangled on scales much smaller than the beam. However, the combination of unresolved (m_net ≈ 0.07) and EHT-resolved (〈∣m∣〉 ≈ 0.25) linear polarization measurements constrains the degree of order in the true, underlying polarization pattern on scales smaller than our beam size, disallowing significant spatially unresolved disorder.

As a simple toy model, we analyzed a thin, circular ring with polarization confined to two azimuthal Fourier modes, labeled with index ℓ. ¹⁵⁵ First, we include a constant (ℓ = 0) mode that defines m_net. We fix the amplitude of this mode to be 0.07 to match unresolved observations of Sgr A*. Next, we add a second mode with varying index ℓ > 0 and an amplitude of 0.7, similar to the peak fractional polarization expected for synchrotron emission. By varying ℓ, we can crudely assess the allowed degree of coherence in the polarization of Sgr A*.

Figure 3 shows the resolved fractional polarization 〈∣m∣〉 at an angular resolution of 20 μas as a function of the secondary mode index ℓ. Both a perfectly ordered polarization field (ℓ = 0) and a highly disordered polarization field (ℓ ≫ 1) will have m_net ≈ 〈∣m∣〉. For the former, there is no beam depolarization; for the latter, the beam depolarization eliminates all small-scale polarized power, even at the resolution of the EHT. Hence, the high value of 〈∣m∣〉 relative to m_net that we observe is a powerful diagnostic of coherent polarized structure.

Zoom In Zoom Out Reset image size

As expected, small values of ℓ produce resolved polarization fractions that are too high, while large values of ℓ produce resolved polarization fractions that are too low. Many effects that are not included in this toy model could further decrease the resolved fractional polarization—the amplitude of the small-scale polarization structure could be significantly less than the synchrotron maximum (e.g., from optical depth or Faraday depolarization), there could be a mix of more than two modes, and there could be radial polarization structure that causes beam depolarization. Hence, this example provides a conservative lower limit on the scale of coherent polarized structure. To be consistent with our measurements of Sgr A*, we require ℓ ≲ 4, corresponding to structure on angular scales of $\theta \approx \tfrac{\pi }{{\ell }}5{\theta }_{{\rm{g}}}\approx 20\,\mu \mathrm{as}$. Here θ_g = r_g /d, where d is the distance and 5θ_g is the approximate radius of the emission ring in Sgr A*. Hence, even without detailed modeling, we anticipate that the underlying polarization in Sgr A* is highly ordered, with significant power on azimuthal scales of θ ≈ 4M or more. That is, the large resolved polarization fraction implies relative order of the magnetic field pattern on scales below the beam size.

3.3. Decoding the Polarization Morphology

Semianalytic models enable computationally inexpensive investigation of the effects of model parameters on images. For example, semianalytic models of radiatively inefficient accretion flows have been used for decades to gain tractable yet physically motivated insights into accretion flows (Bromley et al. 2001; Broderick et al. 2009, 2011, 2014, 2016; Pu et al. 2016; Pu & Broderick 2018; Vincent et al. 2022). Here we explore a very simple model, KerrBAM (or Kerr Bayesian Accretion Modeling), a semianalytic model for equatorial, axisymmetric synchrotron emission around a Kerr BH (Palumbo et al. 2022). This modeling framework carries out ray-tracing in a Kerr spacetime to produce a model image assuming an equatorial ring of emission with a specified fluid velocity, magnetic field geometry, and radial emission profile. Here we use this simple model to illustrate the effects of inclination and spin on polarized image structure.

As our starting point, we average ¹⁵⁶ magnetic fields and velocity fields in three KHARMA GRMHD simulations (to be discussed in Section 4) in both time and azimuth. We specify a ring of emission centered at a radius of 6r_g and use the values of the fluid velocity and magnetic field extracted from the GRMHD midplane at this radius. ¹⁵⁷ To give the emission ring a realistically finite width, the emission is spread in a Gaussian spanning approximately 4r_g –8r_g , keeping the velocity and magnetic field vectors constant. With these values, KerrBAM is able to capture the effects of beaming, frame dragging, and lensing on the resultant image. Note that this model excludes the likely contribution of emission off the midplane (e.g., Falcke et al. 1993; Markoff et al. 2007).

For three different magnetically arrested disks (MADs) with spins of 0, +0.5, and +0.94, we plot several polarimetric quantities of interest (leftmost column) and their model images (subsequent columns) in Figure 4. Along with the polarimetric observables, we overlay our constraints in gray, where for ∠β₂ the range without RM derotation is shown as a hatched region. Since this model places emission exactly at the midplane by construction, images produced at inclinations too close to 90° are misleading and therefore not included. The KerrBAM prescription does not include Faraday effects, only crudely modeling optical depth (in this case applying a midplane-normal crossing optical depth τ_⊥ = 0.5 applied uniformly to ${ \mathcal I }$, ${ \mathcal Q }$, and ${ \mathcal U }$) and assuming a prespecified emission model confined to the midplane, so detailed agreement with the GRMHD models is neither expected nor achieved. Nevertheless, this model is useful for understanding several qualitative trends in our GRMHD library that are successfully reproduced.

Zoom In Zoom Out Reset image size

First, the net polarization is minimized at low inclination, since the symmetry of the accretion flow causes cancellation of polarization in the integrated image. The amplitude of the rotationally invariant mode β₂ is always high, due to the underlying azimuthal symmetry of the system. Meanwhile, the amplitude of ∣β₁∣ is stronger at higher inclinations, as it is sensitive to asymmetries in the polarized image. Finally, we highlight the spin dependence of ∠β₂, which this modeling demonstrates is driven by the evolution of the magnetic field and velocity structure in the GRMHD simulations due to frame dragging (see also Ricarte et al. 2022; Chael et al. 2023; Emami et al. 2023b). The a_* = 0 model has ∠β₂ ∼ −180°, corresponding to a very toroidal EVPA pattern and thus radial magnetic fields. Meanwhile, the higher spin models acquire −180 ≲ ∠β₂ ≲ 0 owing to their more spiral EVPA structures. Interestingly, ∠β₂ remains strikingly stable with inclination, although the overall image structure appears to evolve substantially by eye.

This exploration shows that some of the most salient qualitative features of the polarized image can be traced back to fundamental properties of the fluid and spacetime (magnetic field geometry and spin) without necessarily invoking more uncertain aspects of GRMHD models such as Faraday rotation, the electron-to-ion temperature ratio, and the electron distribution function. However, more physically complete calculations with GRMHD simulations that include these effects self-consistently are still necessary for quantitative comparison.

While semianalytic models provide qualitative insights and intuition about BH accretion flows, they do not enforce conservation laws or capture time-dependent phenomena such as turbulence and shocks that play a crucial role in determining the detailed system structure. Thus, we generate dynamical source models using numerical ideal GRMHD simulations. A fluid approximation would appear to conflict with the fact that the rate of Coulomb collisions is small, leading to mean free paths well exceeding the system size, implying that a collisionless kinetic treatment of the plasma may be necessary (Mahadevan & Quataert 1997). However, kinetic instabilities can produce small-scale inhomogeneities in the magnetic field that produce an effective collisionality through particle–wave interactions (Kunz et al. 2014; Riquelme et al. 2015; Sironi & Narayan 2015; Meyrand et al. 2019). We implicitly assume that radiative effects like cooling are not dynamically important for the fluid evolution. This assumption is well motivated given the low accretion rate of Sgr A*, $\dot{M}\lesssim {10}^{-6}{\dot{M}}_{\mathrm{Edd}}$, for which the radiative cooling timescale is long compared to the accretion timescale (Dibi et al. 2012; Ryan et al. 2017; Chael et al. 2018; Porth et al. 2019; but see also Yoon et al. 2020).

In Paper V, to compare with total intensity EHT and multiwavelength constraints, we generated a suite of GRMHD-derived images sampling a range of initial conditions and parameterizations of the electron temperature and distribution function. We simplify our exploration in this work, limiting ourselves to simulations with untilted torus-like initial conditions, relativistic thermal electron distribution functions (eDFs) lacking nonthermal contributions, and electron temperatures prescribed via the Mościbrodzka et al. (2016) R − β prescription (see Equation (8) below). The properties of our GRMHD simulations are summarized in Table 2. Radiative transfer is integrated within a radius of 100r_g , explicitly ignoring material in highly magnetized regions with σ ≡ b²/ρ > 1, within which mass density is artificially injected to keep the simulation stable. We briefly test the impact of our choices of outer integration radius, the σ cut, and the eDF in Appendices D–F, respectively. While departures from these assumptions are both interesting and physically justified, we defer a thorough investigation of these topics to future work.

Our GRMHD library samples a five-dimensional parameter space. The first parameter is the magnetic field state, either an MAD model (Bisnovatyi-Kogan & Ruzmaikin 1976; Igumenshchev et al. 2003; Narayan et al. 2003; Tchekhovskoy et al. 2011) or a standard and normal evolution (SANE) model (De Villiers et al. 2003; Gammie et al. 2003; Narayan et al. 2012; Sądowski et al. 2013). These describe models in which the magnetic flux threading the horizon for a given accretion rate has saturated and become dynamically important (MAD) or not (SANE). The second is the BH spin, which we denote as a_* ∈ [ − 1, 1], where a negative sign indicates a retrograde disk with respect to the spin vector. Third is the inclination, which uniformly samples i ∈ [0°, 180°], instead of only i ∈ [0°, 90°] as probed in Paper V, because Faraday rotation and emission of circular polarization break the symmetry when polarization is considered. Our fourth parameter is R_high, which sets the asymptotic value of the ion-to-electron temperature ratio as plasma β → ∞ (Mościbrodzka et al. 2016). Specifically,

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{i}}{{T}_{e}}={R}_{\mathrm{low}}\displaystyle \frac{1}{1+{\beta }^{2}}+{R}_{\mathrm{high}}\displaystyle \frac{{\beta }^{2}}{1+{\beta }^{2}},\end{eqnarray} \tag{ 8 }$$

where T_i and T_e are the ion and electron temperatures, respectively. While the potential importance of electron cooling for M87* motivated models with cooler electrons, R_low = 10, here we only consider R_low = 1 owing to the much smaller Eddington ratio of Sgr A*. Finally, our fifth parameter is the magnetic field polarity with respect to the angular momentum vector of the disk, either aligned or reversed, which affects the direction of Faraday rotation and the handedness of circularly polarized emission. This last degree of freedom only matters for polarized radiative transfer and was ignored in Paper V. We produce a library of images for each combination of these parameters, tabulated in Table 3.

We retain the use of multiple codes to assess numerical systematic differences. For scoring, we generate libraries spanning 15,000t_g (t_g ≡ r_g /c), equivalent to about 10 8 hr nights of observation for the parameter combinations listed in Table 3 using two code combinations: KHARMA ¹⁵⁸ (Prather et al. 2021) + IPOLE ¹⁵⁹ (Mościbrodzka & Gammie 2018) and BHAC ¹⁶⁰ (Porth et al. 2017; Olivares et al. 2019) + RAPTOR ¹⁶¹ (Bronzwaer et al. 2018, 2020), where the first and second codes in each pair correspond to GRMHD and GRRT, respectively. As a further consistency check, a third set is generated with H-AMR ¹⁶² (Liska et al. 2022) + IPOLE for a subset of parameter space (only i ≤ 90°, aligned fields, and 5000t_g ) that we do not use for scoring.

Each simulation is initialized with a torus of gas in constant angular momentum hydrodynamic equilibrium (Fishbone & Moncrief 1976). These tori are perturbed with a weak, poloidal magnetic field. The simulations vary in their initial radius of maximum pressure (from ∼15r_g to 40r_g ) and adiabatic index, Γ_ad. Codes differ in their choice of Γ_ad because Γ_ad = 4/3 applies to a fluid of relativistic electrons and Γ_ad = 5/3 applies to a fluid of nonrelativistic ions, but only one fluid is evolved in these models. Depending on the torus size and initial magnetic field configuration, the simulations develop into an MAD or SANE state (see, e.g., Wong et al. 2022).

In Figure 5, we plot a selection of time-averaged GRMHD snapshots from our library, blurred to EHT resolution using a Gaussian convolution kernel with an FWHM of 20 μas. In the left panel of each set we plot total intensity in gray scale and the resolved linear polarization as colored ticks. In the right panel of each set, we plot the circular polarization from blue to red with total intensity contours. Each panel is individually normalized such that the color maps span from 0 to the $\max ({ \mathcal I })$ on the left and $\pm \max (| { \mathcal V }| )$ on the right. Each of these models is an MAD a_* = 0.94, R_high = 40 aligned field simulation, computed with different codes as indicated above.

Zoom In Zoom Out Reset image size

The codes exhibit agreement in terms of total intensity and polarized morphology but differ somewhat in the degree of polarization. As the inclination grows, the total intensity image becomes more asymmetric owing to Doppler beaming (e.g., Falcke et al. 2000; Medeiros et al. 2022; Paper V). The same holds true for the polarization, which is further affected by a Faraday depolarization gradient (see Appendix A.3). The magnetic field geometry as sampled by deflected light rays is encoded in the image of circular polarization. In particular, edge-on images in circular polarization exhibit sign inversions along both a horizontal and vertical axis due to flips in the line-of-sight magnetic field direction, and this signal disappears as the viewing angle decreases (Ricarte et al. 2021; Tsunetoe et al. 2021).

We introduce a novel methodology to score each of our GRMHD models using the eight polarimetric constraints in Table 1. Our new scoring scheme acts on time-averaged GRMHD images and attempts to accommodate variations between codes. Note that we only include quantities inferred from our polarimetric images in these constraints, but we will discuss comparisons with total intensity and multifrequency constraints derived in Paper V.

• 1.  
  First, each model time series of images is split into 10 windows, each with 1500 M duration. Within each window, we produce a time-averaged image by averaging each of the Stokes parameters. Then, we blur the average image with a Gaussian kernel with an FWHM of 20 μas and compute each of the eight observables for scoring.
• 2.  
  For each combination of parameters, we combine the values of the observables predicted by the KHARMA and BHAC codes. Since there are 10 windows and two sets of codes, this results in 20 different samples. From these values, we compute the 90% quantiles ¹⁶³ of each observable to capture the time variability.
• 3.  
  A model passes an individual observational constraint if there is overlap between its 90% quantile region and that of the observations. A model passes a set of observational constraints if it passes all of the constraints in the set simultaneously.

The most important differences compared to the scoring system utilized in Paper V are that this new system operates on time-averaged images and combines the results from multiple codes into a single theoretical range. We tested performing scoring using only one simulation set at a time. Since KHARMA model electron temperatures are assigned systematically hotter than those of the BHAC models (see Appendix H), KHARMA passes models with larger R_high. There is more disagreement between the codes for SANE models than for MAD models. The constraints with the most disagreement between the two codes are ∠β₂, ∣β₂∣/∣β₁∣, and m_net, with the KHARMA simulations ruling out more SANE models than the BHAC simulations in each case.

Each of the observational constraints has known connections with the underlying physics. For brevity, we defer a pedagogical exploration of how each of our free parameters is imprinted onto the observables to Appendix A. We study how each individual constraint affects model selection in Appendix B. Here we summarize the highest-level scoring results, first excluding ∠β₂ and then including ∠β₂ either as observed or after performing RM derotation.

5.1. Constraints Independent of RM

In Figure 6, we plot a pass/fail table combining all polarimetric constraints, with the exception of ∠β₂. These plots combine both polarities of the magnetic field, showing a pass as long as either polarity passes. These tables are slightly but not systematically different as a function of magnetic field polarity.

Zoom In Zoom Out Reset image size

We find that the tight constraint on 〈∣m∣〉 (24%–28%) is the most powerful, driving most of the trends shown in this figure. It is much more constraining on parameter space than m_net, for which a much larger range (2.0%–7.3%) is allowed. The ∣β₂∣ constraint rules out a few additional typically edge-on models, but it does not provide too much more additional constraining power because 〈∣m∣〉 and ∣β₂∣ are correlated. Without ∠β₂, Figure 6 reveals no significant preference between i > 90° and i < 90° models.

While our total intensity constraints generally favored larger values of R_high (due largely to multiwavelength constraints; Paper V), our polarimetric constraints usually prefer more moderate values. This is because larger values of R_high usually lead to larger internal Faraday rotation depths (see Appendix A.4), which is the most important physical driver of depolarization in our models. However, an interesting trend with respect to spin allows one of the best-bet models of Paper V to continue to pass with R_high = 160. This is the MAD a_* = 0.94, R_high = 160, i = 30°/150° model. MAD models with larger spin have smaller Faraday rotation depths (see Appendix H), allowing them to pass the 〈∣m∣〉 constraint for larger values of R_high. We refer readers to Appendix B for a more detailed breakdown of each constraint considered individually.

5.2. Constraints Including ∠β₂ without RM Derotation

First, we discuss the ∠β₂ constraint if RM derotation is not performed. It is possible that the RM may be attributed entirely to Faraday rotation captured within our simulation domain. GRMHD models are capable of producing the correct magnitude of RM from Faraday rotation on event horizon scales, but they tend to produce RM sign flips that are not consistent with decades of Sgr A* observations that produce negative values of the RM (Ricarte et al. 2020; M87* Paper VII; Wielgus et al. 2024). However, it is possible that this problem is related to the excess variability in our models identified in Paper V. We further discuss the uncertainties surrounding our interpretation of the RM in Appendix C.

If one attributes the RM entirely to internal Faraday rotation, then our constraint on ∠β₂ spans the interval (125°, 160°). Adding this constraint to Figure 6 results in Figure 7. A selection for i < 90° arises because the handedness of the polarization spiral is opposite that of the magnetic field, which inherits the handedness of the inflowing and emitting gas (see Section 3.3 and Appendix A.3). This corresponds to counterclockwise motion, which disagrees with hot spot interpretations of polarized flares both in the near-IR (NIR; GRAVITY Collaboration et al. 2018, 2020a, 2020b) and in the submillimeter (Vos et al. 2022; Wielgus et al. 2022b). That is, consistency with clockwise motion would require −180° < ∠β₂ < 0° if we assume that ∠β₂ traces magnetic field lines with outgoing Poynting flux (Chael et al. 2023), which does not agree with the linearly polarized morphology as observed on the sky.

Zoom In Zoom Out Reset image size

Without RM derotation, no model can simultaneously pass all total intensity and polarimetric constraints. This is because the a_* = 0.94 best-bet model of Paper V produces an EVPA pattern that is too radial (see Appendix A.2). All models that pass our polarization constraints in Figure 7 fail multiple constraints on the total intensity. In particular, all eight models shown in Figure 7 produce too much flux in the infrared to match observations, and all but the SANE model at a_* = 0.94 overproduce the X-ray flux (Paper V). Both of these are serious failures, as both the IR and X-ray fluxes estimated by our models are lower limits owing to our lack of nonthermal electrons and small simulation domain relative to the X-ray-emitting area. Five of the models additionally fail to match the observed size and flux of the source at 86 GHz (Issaoun et al. 2019). All of these models also fail at least one total intensity structural constraint (m-ring and visibility amplitude morphology tests in Paper V). In conclusion, we cannot find a concordance model of Sgr A* without RM derotation.

5.3. Constraints Including ∠β₂ with RM Derotation

Alternatively, in this section we interpret the mean RM as an external Faraday screen, motivating derotation. As discussed in Section 2, ∠β₂ depends on twice the RM, for which a mean value of $\langle \mathrm{RM}\rangle =-{4.65}_{-1.18}^{+1.25}\times {10}^{5}\ \mathrm{rad}\ {{\rm{m}}}^{-2}$ has been obtained. This potentially results in a shift in ∠β₂ of $2\langle \mathrm{RM}\rangle {\lambda }^{2}=-{92.0}_{-23.4}^{+24.7}$ deg if this RM is interpreted as an external Faraday screen. In this picture, a relatively stable external screen explains the constant sign of RM that has been observed for decades (nevertheless with variation on the order of ∼10⁵ rad m⁻²). Then, an additional component on event horizon scales, which is already included self-consistently in our models, explains the subhour time variability.

If one attributes the mean RM of a given day entirely to an external screen, then our constraint on ∠β₂ spans (−168°, −85°). Adding this constraint to Figure 6 results in Figure 8. Performing this cut requires inclination angles >90°, corresponding to clockwise motion on the sky, which now agrees with the aforementioned models of polarized NIR and submillimeter flares.

Zoom In Zoom Out Reset image size

With RM derotation, one of the best-bet models from our total intensity analysis passes all applied total intensity and polarimetric constraints. This is the MAD a_* = 0.94, R_high = 160, i = 150° aligned model. The second best-bet model from Paper V had a_* = 0.5 and otherwise identical parameters. This second model passes all constraints except 〈∣m∣〉, which it underproduces by ∼3%. In order for the a_* = 0.94 best-bet model to pass, at least 97% of the measured RM must arise from an external screen. Notably, the best-bet model fails if the smaller RM measured at 86 GHz a few days prior, (−2.14 ± 0.51) × 10⁵ rad m⁻² (Wielgus et al. 2024), is instead interpreted as the external screen.

In Figure 9, we visualize the best-bet model (BHAC shown) that survives with RM derotation. In the left two columns, we plot its full polarimetric image in the style of Figure 5. No blurring is applied in the leftmost column, and a 20 μas FWHM Gaussian kernel is convolved with the image in the second column to approximate EHT resolution. This model features a bright photon ring, and in our image without blurring, we omit total intensity contours from the circular polarization map to reveal a photon ring sign inversion (discussed in Mościbrodzka et al. 2021; Ricarte et al. 2021).

Zoom In Zoom Out Reset image size

On the right, we produce a map of the density of the observed emission in the equivalent KHARMA simulation (using Kerr–Schild coordinates). The emission density map is normalized such that its peak value is unity, and it is visualized in logarithmic scale with 3 orders of magnitude in dynamic range. Our line of sight is indicated by the green arrow, and a white contour encloses 90% of the total emission. This reveals that while the emission is peaked at small radius near the disk midplane, a substantial fraction of the emission originates from a more diffuse jet funnel region. Computing an emission-weighted characteristic emission radius $\bar{x}\equiv \int x\epsilon {dV}/\int \epsilon {dV}$, where is the emission density and x is the radius in cylindrical coordinates, we find $\bar{x}=7.3$. We note that our choices to include only thermal electron distribution functions and cut out regions with σ > 1 in this work minimize the potential contribution of a jet to the total emission (e.g., Figure 12 of Fromm et al. 2022). A significant jet component may be necessary to reproduce the flat spectral index at these frequencies (Falcke et al. 1993; Falcke & Markoff 2000; Mościbrodzka & Falcke 2013).

At a radius of 7.3r_g , we compute a mass-weighted average magnetic field strength of ${26}_{-4}^{+3}\ {\rm{G}}$, where the range quoted here corresponds to the 16th to 84th percentile values obtained in the time series. This value agrees reasonably well with the one-zone model discussed in Section 3.1, although we note that this value evolves substantially with radius, reaching ${67}_{-9}^{+8}\ {\rm{G}}$ at a radius of 4r_g and ${560}_{-80}^{+80}\ {\rm{G}}$ at the horizon.

This model produces an outflow power of 4 × 10³⁸ erg s⁻¹ and has an accretion rate of 5 × 10⁻⁹ M_⊙ yr⁻¹. This model has a very large jet efficiency of approximately 150% powered by the Blandford & Znajek (1977) mechanism. Yet despite its efficiently, the jet's power is not high enough to expect global effects on the evolution of our Galaxy (e.g., Su et al. 2021).

The first polarized image of Sgr A* on event horizon scales exhibits a high resolved polarization fraction of 24%–28% and an ordered, rotationally symmetric EVPA pattern. Through semianalytic arguments and comparisons to GRMHD simulations, we come to the following conclusions:

• 1.  
  The large resolved polarization fraction implies that the magnetic field on event horizon scales cannot be very tangled on scales smaller than beam, nor can Faraday rotation add too much additional disorder to the EVPA structure. The disparity between the spatially resolved (24%–28%) and unresolved (2.0%–7.3%) linear polarization fractions can be attributed to cancellations due to the symmetric nature of the image.
• 2.  
  Driven mostly by the spatially resolved polarization fraction, our constraints strongly favor MAD models over their SANE counterparts, as in M87* Paper VIII.
• 3.  
  If we rely on internal Faraday rotation to produce the observed RM and do not perform derotation, then there is no model that passes all total intensity and polarimetric constraints.
• 4.  
  On the other hand, if we assume that the RM can be attributed to an external screen and derotate the EVPA pattern, then we find one model that passes all applied total intensity and polarimetric constraints: MAD a_* = 0.94, R_high = 160, i = 150° aligned.

While our ideal GRMHD simulations containing only thermal electron distributions have done remarkably well at reproducing many of the observed quantities of Sgr A*, they nevertheless have many known imperfections. Most of these models overestimate time variability, including the best-bet model (Paper V), and we caution that the values inferred from our best-bet model should not be interpreted as measurements. Known areas where these simulations can be improved include the following:

• 1.  
  Initial Conditions : All of our simulations are initialized with tori that are either perfectly aligned or antialigned with the BH angular momentum axis. Simulations feeding the BH via stellar winds have different variability characteristics (Murchikova et al. 2022) and can self-consistently predict an external Faraday screen (Ressler et al. 2019, 2023). Tilted disk models (e.g., Fragile et al. 2007; Liska et al. 2018; Chatterjee et al. 2020) may lead to different Faraday rotation characteristics owing to their geometry at large radii.
• 2.  
  Electron Thermodynamics : The Mościbrodzka et al. (2016) prescription that we adopt to set the electron temperature broadly captures the trends seen in kinetic simulations that explicitly model heating and cooling (e.g., Chael et al. 2018; Dexter et al. 2020; Mizuno et al. 2021; Dihingia et al. 2023) but does not reproduce them in much detail. More generally, a nonthermal contribution to the electron distribution function is believed to be necessary to reproduce the spectral energy distribution (Özel et al. 2000; Markoff et al. 2001; Davelaar et al. 2018) and is naturally predicted by particle-in-cell simulations (Kunz et al. 2016; Ball et al. 2018). Nonthermal electron distribution functions can have significant impacts on both total intensity and polarized properties (e.g., Markoff et al. 2001; Mao et al. 2017; Davelaar et al. 2018; Cruz-Osorio et al. 2022; Fromm et al. 2022; Paper V) and are a promising avenue to continue theoretical exploration.
• 3.  
  Plasma Composition : Wong & Gammie (2022) demonstrate that models fed by helium rather than hydrogen may have substantially different emission morphologies, tending toward higher temperatures and lower densities and thus higher polarization fractions. Meanwhile, the presence of electron–positron pairs can significantly alter Faraday effects, leading to potential signatures in both linear and circular polarization that have not been fully explored (Anantua et al. 2020; Emami et al. 2021, 2023a; M87* Paper IX).

Several ongoing developments within the EHT will be impactful for testing our present interpretation, especially explorations in time and frequency. An effort is ongoing to produce dynamical movies of Sgr A*, despite the challenges of very sparse snapshot (u,v) coverage (Tiede et al. 2020; Farah et al. 2022; Levis et al. 2023). Measurements of the apparent angular velocity or potentially the motion of hot spots will provide additional constraints on spin and inclination (Wielgus et al. 2022b; Conroy et al. 2023). The dynamic reconstruction and geometric modeling of these data by Knollmüller et al. (2023) are consistent with the inferred inclination and clockwise motion of our best-bet model. On longer timescales (of years), it will be important to obtain averages of quantities such as ∠β₂, which varies little in our models owing to its tight link with BH spin.

In the frequency domain, future EHT data sets will include 345 GHz data. The wavelength dependence of the scattering screen toward the Galactic center inhibits imaging of Sgr A* at lower frequencies below 86 GHz (Johnson et al. 2018; Issaoun et al. 2019, 2021). On its own, a 345 GHz polarized image would already strongly mitigate one of our largest systematic uncertainties, the RM; the total EVPA rotation would decrease by a factor of 2, as it is proportional to ν⁻². These images will also have intrinsically higher resolution by a factor of 50%. Simultaneous dual-band observations could enable the production of RM maps, which would be our best tool for characterizing the Faraday screen and disambiguating our approach to derotation. If the RM truly originates from an external Faraday screen and the emission origin does not significantly change, then at 345 GHz we should observe a spatially uniform EVPA rotation of ∼20° clockwise relative to our 230 GHz image (roughly halfway between the top two rows in Figure 1). Meanwhile, RM due to internal Faraday rotation may exhibit more spatial variation and potentially sign flips owing to turbulence in the inner accretion flow (Ricarte et al. 2020).

Given the vastness of parameter and modeling space available to theoretical interpretation, we expect the polarized image of Sgr A* to continue to constrain models for many studies to come. This growing EHT data set will continue to challenge theoretical models and provide insights into the nature of BHs, accretion, and plasma physics.

The Event Horizon Telescope Collaboration thanks the following organizations and programs: the Academia Sinica; the Academy of Finland (projects 274477, 284495, 312496, 315721); the Agencia Nacional de Investigación y Desarrollo (ANID), Chile via NCN19_058 (TITANs), Fondecyt 1221421 and BASAL FB210003; the Alexander von Humboldt Stiftung; 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 program by MCIN (with funding from European Union NextGenerationEU, PRTR-C17I1); the Black Hole Initiative, which is funded by grants from the John Templeton Foundation (60477, 61497, 62286) and the Gordon and Betty Moore Foundation (Grant GBMF-8273)—although the opinions expressed in this work are those of the author and do not necessarily reflect the views of these Foundations; the Brinson Foundation; "la Caixa" Foundation (ID 100010434) through fellowship codes LCF/BQ/DI22/11940027 and LCF/BQ/DI22/11940030; 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 program 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's Horizon 2020 research and innovation program under grant agreements RadioNet (No. 730562) and M2FINDERS (No. 101018682); the Horizon ERC Grants 2021 program under grant agreement No. 101040021; 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 Fondo CAS-ANID folio CAS220010; the Generalitat Valenciana (grants APOSTD/2018/177 and ASFAE/2022/018) and GenT Program (project CIDEGENT/2018/021); the Gordon and Betty Moore Foundation (GBMF-3561, GBMF-5278, GBMF-10423); the Institute for Advanced Study; the 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. 443220636); Joint Columbia/Flatiron Postdoctoral Fellowship (research at the Flatiron Institute is supported by the Simons Foundation); the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT; grant JPMXP1020200109); the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS Research Fellowship (JP17J08829); the Joint Institute for Computational Fundamental Science, Japan; the Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS, grants QYZDJ-SSW-SLH057, QYZDJ-SSW-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 MICINN Research Projects PID2019-108995GB-C22 and PID2022-140888NB-C22; the MIT International Science and Technology Initiatives (MISTI) Funds; the Ministry of Science and Technology (MOST) of Taiwan (103-2119-M-001-010-MY2, 105-2112-M-001-025-MY3, 105-2119-M-001-042, 106-2112-M-001-011, 106-2119-M-001-013, 106-2119-M-001-027, 106-2923-M-001-005, 107-2119-M-001-017, 107-2119-M-001-020, 107-2119-M-001-041, 107-2119-M-110-005, 107-2923-M-001-009, 108-2112-M-001-048, 108-2112-M-001-051, 108-2923-M-001-002, 109-2112-M-001-025, 109-2124-M-001-005, 109-2923-M-001-001, 110-2112-M-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 and 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 and Technology Council (NSTC, grants NSTC 111-2112-M-001-041, NSTC 111-2124-M-001-005, NSTC 112-2124-M-001-014); the US National Science Foundation (NSF, grants AST-0096454, AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433, OIA-1126433, AST-1140030, DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509, OPP-1248097, AST-1310896, AST-1440254, AST-1555365, AST-1614868, AST-1615796, AST-1715061, AST-1716327, AST- 1726637,OISE-1743747, AST-1743747, AST-1816420, AST-1952099, AST-1935980, AST-2034306, AST-2205908, 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, 12325302, and 12303021); 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, 2022R1F1A1075115); Netherlands Research School for Astronomy (NOVA) Virtual Institute of Accretion (VIA) postdoctoral fellowships; NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation; Onsala Space Observatory (OSO) national infrastructure, for the provisioning of its facilities/observational support (OSO receives funding through the Swedish Research Council under grant 2017-00648); the Perimeter Institute for Theoretical Physics (research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science); the 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; and 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.

Software: eht-imaging (Chael et al. 2016), Numpy (Harris et al. 2020), Scipy (Jones et al. 2001), Pandas (McKinney 2010), Astropy (The Astropy Collaboration et al 2013, 2018), Jupyter (Kluyver et al. 2016), Matplotlib (Hunter 2007), Themis (Broderick et al. 2020), ipole (Noble et al. 2007; Mościbrodzka & Gammie 2018), KHARMA (Prather et al. 2021), BHAC (Porth et al. 2017; Olivares et al. 2019), H-AMR (Liska et al. 2022), RAPTOR (Bronzwaer et al. 2018, 2020), KerrBAM (Palumbo et al. 2022).

Using our GRMHD models, we explore a five-dimensional parameter space, constrained by eight observable aspects of the polarized image that we believe are tied to the models in physically understood ways. Below, we highlight the most salient trends in our simulated image library to explain their physical origins. We focus on illustrative examples in this appendix, but we provide exhaustive distributions of observables calculated from our GRMHD models in Appendix H.

A.1. Magnetic Field State

By construction, SANE models have weaker magnetic fields near the horizon than their MAD counterparts at a given accretion rate. As a result, once the fluid is rescaled to reproduce the observed millimeter flux, SANE models usually have larger mass densities. This translates directly to a larger Faraday rotation depth, which is directly implicated for scrambling/depolarizing EHT model images (Mościbrodzka et al. 2017; Jiménez-Rosales & Dexter 2018; Ricarte et al. 2020; M87* Paper IX). Our SANE models are also colder (Paper V), which further increases the efficiency of Faraday rotation (Jones & O'Dell 1977; Quataert & Gruzinov 2000).

In Figure 10, we explore the differences between our MAD and SANE models with fixed a_* = 0.5, R_high = 40, i = 50° and aligned fields. In the top panels, we plot the time-averaged KHARMA images in total intensity and linear polarization, blurred to a resolution of 20 μas. In the bottom panels, we compare differences in resolved linear and circular polarization fraction, Faraday rotation depth, and β₂. In these and the following plots in this section, we display theoretical error ranges that are calculated based on differences between codes, time variability, and nearest neighbors in parameter space. The details of how these theoretical error bars are calculated are provided in Section 5.

Zoom In Zoom Out Reset image size

Here we see that the SANE model has much lower linear polarization fraction (〈∣m∣〉 and ∣β₂∣), which can be attributed to a much larger Faraday depth ($\langle {\tau }_{{\rho }_{V}}\rangle$). ¹⁶⁴ Much larger Faraday depths in SANE models than their equivalent MADs drive most of the differences between these two classes of models. SANE models can also produce larger circular polarization (〈∣v∣〉) owing to Faraday conversion (M87* Paper IX). Palumbo et al. (2020) showed that ∣β₂∣ is a strong discriminant between MAD and SANE models of M87*. As expected, ∣β₂∣ is significantly larger for the MAD model than for the SANE model. Interestingly, while SANE models of M87* usually exhibit ∠β₂ ∼ 0, corresponding to radial EVPA patterns, the EVPA pattern in this SANE model acquires some twist owing to a tilted forward jet that we view in projection (top left).

A.2. Spin

The BH spin is a particularly interesting quantity to constrain owing to implications for its cosmic assembly and feedback processes. A number of EHT-related studies have recently explored signatures of spin, and resolved linear polarization structure has been shown to be one of the most promising and accessible probes (Palumbo et al. 2020; Chael et al. 2023; Emami et al. 2023b; Qiu et al. 2023; Ricarte et al. 2023).

In Figure 11, we plot the phase and amplitude of β₂ as a function of spin for the subset of the MAD R_high = 10, i = 30°/150° reversed models. The outer accretion disk rotates counterclockwise on the sky for i = 30° and clockwise on the sky for i = 150°, which is reflected by the sign of ∠β₂ (or rather, the sign of its imaginary component). As discussed in Section 3.3, ∠β₂ evolves with spin owing to frame dragging, which results in changes in the magnetic field and velocity structure (Palumbo et al. 2020; Event Horizon Telescope Collaboration et al. 2021b; Ricarte et al. 2022; Chael et al. 2023; Emami et al. 2023b; Qiu et al. 2023). The most highly spinning prograde models acquire a strong azimuthal magnetic field component, resulting in more radial EVPA patterns (∠β₂ closer to 0°). Finally, ∣β₂∣ is stronger for symmetric and ordered progrades than for their messier retrograde counterparts (see also Qiu et al. 2023).

Zoom In Zoom Out Reset image size

A.3. Inclination

The inclination of Sgr A* is of particular interest because its polarized flaring activity can be interpreted with a polarized hot spot model that favors a relatively face-on viewing angle (GRAVITY Collaboration et al. 2020a, 2020b; Wielgus et al. 2022b). In addition, it is of interest whether or not the accretion disk or BH angular momentum axes align with any structure in its environment.

Inclination is imprinted on the polarized image in a variety of ways, and we plot most of our polarimetric observables as a function of inclination in Figure 12. Here MAD a_* = 0.94, R_high = 10 models are considered. These models produce rotationally symmetric images when viewed face-on, and thus cancellation leads to opposite behavior of m_net and 〈∣m∣〉, the latter of which decreases with inclination owing to Faraday depolarization. Intuitively, ∣β₂∣, the amplitude of the rotationally invariant mode, is strongest for face-on viewing angles and weakest for edge-on viewing angles. Meanwhile, the asymmetric β₁ mode has the largest amplitude for intermediate inclinations. The handedness of the linear polarization spiral is directly encoded in $\mathrm{sign}(\mathrm{Im}({\beta }_{2}))$, and thus we see that ∠β₂ > 0° for i < 90° and ∠β₂ < 0° for i > 90°. Finally, v_net is sensitive to whether the poloidal field is pointed toward us or away from us, but note that it is not perfectly antisymmetric about i = 90° owing to contributions from Faraday conversion (Ricarte et al. 2021).

Zoom In Zoom Out Reset image size

A.4.  R_high (Electron Temperature)

As described in Section 4, R_high sets the ratio of ion-to-electron temperature as plasma β → ∞ (Mościbrodzka et al. 2016). Increasing R_high while fixing all other parameters makes the electrons of a given model cooler and less efficient emitters. Thus, models with larger R_high tend to have larger values of ${ \mathcal M }$ when rescaled to achieve the same target flux. As a result, increasing R_high indirectly increases the Faraday rotation depth (Mościbrodzka et al. 2017; Jiménez-Rosales & Dexter 2018; Ricarte et al. 2020; M87* Paper VIII). Increasing R_high also shifts emission away from the midplane and concentrates it toward the jet funnel region (Paper V; Wong et al. 2022). This effect is much weaker for MADs than for SANEs, since MAD models intrinsically have smaller plasma β on horizon scales.

In Figure 13, we plot time-averaged BHAC MAD a_* = 0.5, i = 130° aligned field models as a function of R_high, as well as several of their linear polarization observables. Increasing Faraday depolarization explains the declines in 〈∣m∣〉 and ∣β₂∣ with R_high. The polarization grows more asymmetric as R_high increases, because at this inclination the Faraday thick midplane is at the top half of the image. This, combined with increased Faraday rotation that slightly turns ticks clockwise, ¹⁶⁵ leads to a shift in ∠β₂. In addition, ∣β₂∣/∣β₁∣ decreases as the polarization grows more asymmetric.

Zoom In Zoom Out Reset image size

A.5. Magnetic Field Polarity

In ideal GRMHD, the equations governing the evolution of a magnetized fluid are invariant to a sign flip of the magnetic field direction. However, the equations of GRRT are not, leading to potential polarimetric signatures of the poloidal field direction. When performing radiative transfer, j_V (intrinsic circular polarization of emitted radiation) and ρ_V (Faraday rotation) are each sensitive to the direction of the field with respect to the photon wavevector. The historically negative Stokes V of Sgr A* is suggestive of a magnetic field oriented away from us. However, M87* Paper IX discusses how flipping the magnetic field direction can have nontrivial effects on the circularly polarized image (beyond a simple sign flip) and noticeable effects on ∠β₂ due to Faraday effects (see also Ricarte et al. 2021; Emami et al. 2023b).

In Figure 14, we highlight the differences between aligned and reversed field models for the time-averaged KHARMA MAD a_* = 0.5, R_high = 160, i = 130° models. Each model is blurred with a 20 μas FWHM Gaussian beam shown in total intensity and linear polarization ticks on the left and circular polarization and total intensity contours on the right. We write ∠β₂ and v_net for each model in the lower left corner, revealing significant and unpredictable differences, motivating independent ray-tracing for each magnetic field polarity. In linear polarization, the difference comes from reversing the direction that Faraday rotation shifts the EVPA pattern. The magnitude of this effect is larger than that reported in M87* Paper IX because M87* models are oriented almost completely face-on, viewed through an evacuated funnel (Ricarte et al. 2020). Models of Sgr A* can accumulate larger Faraday rotation depths as radiation passes through more of the disk at larger inclinations. In circular polarization, this particular model is mostly characterized by an overall sign flip, but this is not uniform across the image, leading to a small difference in v_net. This is because the coefficient of Faraday conversion, which exchanges linear and circular polarization, is invariant to a sign flip in the magnetic field direction.

Zoom In Zoom Out Reset image size

In Section 5, we included a limited selection of plots reflecting which of our models passed each of the polarimetric observational constraints on Sgr A*. Here we break down the impact of each constraint individually.

In Figure 15, we plot the impact of our 〈∣m∣〉 constraint, which we find is the most important for model selection. Compared to the other constraints, 〈∣m∣〉 is measured relatively precisely and the two methods agree very well. The Faraday rotation depth explains the trends in this figure (see Appendix H). More Faraday depolarization tends to occur if R_high is larger, if the inclination is larger, or if the model is SANE. Of the models that fail the 〈∣m∣〉 constraint, most are too depolarized, but some low-R_high, high-spin, face-on models are ruled out for predicting values of 〈∣m∣〉 that are too large. We find that 〈∣m∣〉 is much more constraining than m_net (Figure 16), which is measured much less precisely. Recall that m_net is substantially lower (and less consistent with the light curve) in the m-ring model than Themis. We find that if the higher and tighter m_net constraint from Themis had been adopted on its own, then this would have ruled out many face-on models (explained in Section 3.3 and Appendix A.3), including the a_* = 0.94 best-bet model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our circular polarization constraints are not very impactful. Our upper limit on 〈∣v∣〉 rules out no models (Figure 17), as all GRMHD models produce 〈∣v∣〉 lower than the upper limit (similar to M87* Paper IX). Our constraint on v_net is also not very impactful (Figure 18), but while not visible with our plotting scheme, it does rule out many retrograde models that have aligned fields. These models produce preferentially positive v_net, while decades of Sgr A* observations produce v_net < 0.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our constraints on ∣β₁∣ (Figure 19), ∣β₂∣ (Figure 20), and ∣β₂∣/∣β₁∣ (Figure 21) are impactful, but they are correlated with each other and 〈∣m∣〉. Compared to 〈∣m∣〉, ∣β₂∣ additionally rules out some i = 90° models. The ratio ∣β₂∣/∣β₁∣ is not very constraining, as most models naturally produce ∣β₂∣ > ∣β₁∣, in agreement with the observations. While some methods in Paper VII produced ratios up to ∼5, which would have pushed our selection toward more face-on inclinations, the two methods retained in this paper produced more modest values. Interestingly, a few face-on models are ruled out for being too dominated by the rotationally symmetric mode.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we consider the effect of ∠β₂ both with and without RM derotation in Figures 22 and 23, respectively. In either case, models with preferentially radial EVPA patterns are most likely to fail, such as face-on prograde MAD models (see Appendix A.2). With derotation, this constraint produces a preference for clockwise motion on the sky (i > 90°). Without derotation, the opposite is true, and more models fail outright since the constraint is tighter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The RM of Sgr A* is a significant systematic uncertainty in our work, affecting our interpretation of ∠β₂. The RM is defined as

$$\begin{eqnarray}&&\mathrm{RM}\equiv \displaystyle \frac{{\rm{\Delta }}\chi }{{\rm{\Delta }}{\lambda }^{2}},\end{eqnarray} \tag{ C1 }$$

where χ is the EVPA and λ is the wavelength. If the EVPA of the polarized emission does not intrinsically change with wavelength (due to optical depth), and the polarized emission is situated entirely behind a Faraday screen that is uniform relative to the size of the emitting region, then the RM is related to a path integral along the line of sight via

$$\begin{eqnarray}&&\mathrm{RM}=8.1\times {10}^{5}\ \mathrm{rad}\ {{\rm{m}}}^{-2}{\int }_{\mathrm{source}}^{\mathrm{observer}}{f}_{\mathrm{rel}}({{\rm{\Theta }}}_{e})\displaystyle \frac{{n}_{e}}{1\ {\mathrm{cm}}^{-3}}\displaystyle \frac{{B}_{| | }}{{\rm{G}}}\displaystyle \frac{{ds}}{\mathrm{pc}},\end{eqnarray} \tag{ C2 }$$

where n_e is the electron number density, B_∣∣ is the local magnetic field parallel to the photon wavevector, and ${f}_{\mathrm{rel}}({{\rm{\Theta }}}_{e})\approx \mathrm{log}({{\rm{\Theta }}}_{e})/(2{{\rm{\Theta }}}_{e}^{2})$, a factor causing lower efficiency as electrons become too relativistic (Jones & O'Dell 1977). If the two assumptions above are correct, then the "intrinsic" EVPA pattern can be easily recovered by derotating the EVPA by RMλ².

Sgr A* has exhibited a constant sign of RM for decades (Bower et al. 2018), which supports the interpretation of a stable external Faraday screen. GRMHD simulations including RM from event horizon scales predict ubiquitous sign flips on subhour timescales that are not observed (Ricarte et al. 2020; Ressler et al. 2023; Wielgus et al. 2024). On the other hand, Sgr A* exhibits non- λ² evolution of the EVPA when comparing the 86 GHz and 230 GHz bands. At 86 GHz, the RM on nearly simultaneous days to our observations is only −2 × 10⁵ rad m⁻² compared to −5 × 10⁵ rad m⁻² at 230 GHz (Wielgus et al. 2024). In addition to subhour time variability, this suggests that at least some of the RM must also come from internal Faraday rotation on event horizon scales.

Carefully predicting the RM directly for all of our GRMHD simulations would increase the computational cost by factors of a few (more than 2) with the software utilized in this work. This is because ray-tracing must be performed at different frequencies at nonuniform spacings to resolve potential phase wrapping and non-λ² behavior of the EVPA. Nevertheless, we check the RM for a few snapshots of our models in Figure 24, where the RM is estimated by ray-tracing at 213, 215, 227, and 229 GHz (emulating observations) and then fitting for the slope RM = d χ/d λ². MAD models are plotted in the top row, and SANE models are plotted in the bottom row. Three inclinations are shown: 30° in blue, 50° in orange, and 90° in gray. All models are at R_high = 40 and in an aligned field configuration. Note that these simulations only include material within 100r_g , but ab initio simulations of the accretion of Sgr A* from stellar winds suggest that a steady Faraday screen could potentially be situated at even larger radii (Ressler et al. 2019, 2023).

Zoom In Zoom Out Reset image size

We find that most of our models naturally produce ∣RM∣ ∼ 10⁵ rad m⁻² at at least one point in time, in rough agreement with the observed value. The SANE models, as well as the MADs at 90°, tend toward larger values, similar to models of M87* (Ricarte et al. 2020). However, as in previous works, the RM flips sign in every model at least once. Interestingly, we find similar order-of-magnitude values of RM if Faraday rotation is explicitly switched off during ray-tracing (ρ_V = 0) in some of these models. This suggests that evolving emission origin as a function of frequency may contribute to the inferred RM and its variability.

Our findings in Figure 24 are broadly consistent with an interpretation wherein the rapid time variability of RM is caused by variability on event horizon scales, but the stability of sign is maintained by an external Faraday screen along the line of sight, motivating derotation of ∠β₂. On the other hand, it may also be possible that all of the RM originates from event horizon scales, and our GRMHD models overpredict the variability in RM in the same way that they overpredict variability in total intensity (Paper V). To resolve this, 345 GHz imaging of Sgr A* will be critical; 345 GHz is less affected by Faraday rotation by a factor of (345/230)² ≈ 2. In addition, RM maps produced via simultaneous multifrequency imaging will help determine the nature of the Faraday screen.

Although we are confident that most of the emission in our models originates close to the event horizon (r ≲ 10r_g ), Faraday rotation can originate at much larger radius in our models, more so as the inclination increases (Dexter et al. 2020; Ricarte et al. 2020). This is especially problematic because material at these radii may not have had enough time in the simulation to reach equilibrium. This concern is more important for studies of Sgr A* than for M87* because we view M87* at an inclination of only 17° through an evacuated funnel.

We test the impact of the outer radiative transfer integration radius in Figure 25, where we ray-trace a few KHARMA snapshots at a variety of radii ranging from 30r_g to 300r_g . We focus on ∠β₂, which should be directly affected by Faraday rotation on large scales. Inclinations of both 50° and 90° are considered, with R_high values of both 10 and 160. Fortunately, we find that ∠β₂ appears to have converged for most of these models before 100r_g , where we perform the ray-tracing in this paper. We find that the models that do exhibit substantial evolution with outer integration radius all produce 〈∣m∣〉 lower than observed. Note that SANE models at 90° inclinations with R_high = 160 are the most Faraday-thick models in our library. Models at i = 90 and/or high R_high appear to have the most evolution with respect to the integration radius. This is consistent with the expectation that higher inclinations and higher R_high values will increase the amount of Faraday rotation owing to more photons traveling through dense, cold regions in the GRMHD domain.

Zoom In Zoom Out Reset image size

While ∠β₂ appears to show evolution for some models, the other polarimetry metrics are well converged and show minimal change for all models across integration radius. However, although we have checked the GRRT step, recall that our GRMHD models are only converged within r ≲ 30r_g owing to computational limitations. Exploration with simulations that are valid to larger radii that may produce an external Faraday screen self-consistently (e.g., Ressler et al. 2023) would be an interesting avenue for future analysis.

The polar funnel in the GRMHD simulations is filled with horizon-penetrating field lines and thought to contain plasma with orders of magnitude lower density than the accretion disk. By the same token, the funnel magnetization σ := B²/ρ is believed to be much larger than the magnetization in the disk. Since there are very few emitting particles in the funnel, its contribution to the overall image is expected to be negligible. In practice, to keep the numerical GRMHD evolution stable, σ is not allowed to assume realistic values but is instead capped at moderate values σ ≲ 50–100 by artificially injecting mass (e.g., Porth et al. 2019). Hence, we cannot trust the inflated mass density in this region. Assuming that emission in the σ ≫ 1 funnel should in reality be negligible, we follow the common practice and set all radiation transport coefficients to zero when the magnetization exceeds a critical value σ_cut = 1. This choice is only safe when no σ ≥ 1 regions form naturally in the disk and when the mixing of disk and funnel plasma at the jet wall is inefficient. In this case the gradient in magnetization is steep which means that whether we adopt σ_cut = 1 or, e.g., σ_cut = 25 does not affect the results. In reality, however, finite-resolution effects in the GRMHD simulations, resolved interchange instabilities, and potentially strong disk magnetization can cause a dependence on the adopted threshold value.

Using the BHAC/RAPTOR data, we have carried out spot checks with two "best-bet" models whereby we increase the threshold to σ_cut = 25: model 1 is MAD a_* = 0, R_high = 40, i = 150° aligned, and model 2 is MAD a_* = 0.94, R_high = 160, i = 30° aligned. In either case, the constraints change only by a few percent, e.g., in model 2 the average β₂ phase changed from 63° to 66° and the average net polarization went down from 2.9% to 2.7%. In model 1, the change in average β₂ phase is somewhat larger (going from −97° to −109°), but still small compared to the overall spread of the distributions. This shows that the results on polarized submillimeter emission are quite robust against change in the adopted value of the σ_cut and emission at or within the highly magnetized funnel does not dominate in the model.

Throughout this work, we have considered only thermal eDFs when performing GRRT. Here we briefly explore the impact of nonthermal electrons in the polarimetric properties of one GRMHD model: MAD a_* = 0, R_high = 40, i = 150° aligned. Two nonthermal prescriptions are explored:

• 1.  
  Variable κ: In each cell, a κ distribution (Vasyliunas 1968; Xiao 2006) is applied, using a κ(σ, β) prescription originating from particle-in-cell simulations (Ball et al. 2018; Davelaar et al. 2019).
• 2.  
  κ = 5: A κ distribution with a constant value of κ = 5 is applied globally (Davelaar et al. 2018).

We ray-trace 300 snapshots for each of these cases and compare with the thermal model snapshots. The accretion rate is kept fixed, but we find that the average flux density is 2.3 Jy for all cases. In Figure 26, we plot a selection of polarimetric quantities for these models. Each marker is placed at the median, and the error bars extend to the 16th and 84th percentiles. Overall, we find only subtle differences between these different eDF models. We find that 〈∣m∣〉 declines in the nonthermal eDF models, coincident with increases in the Faraday rotation depth (2.2, 4.2, and 6.3 for thermal, Variable κ, and κ = 5 models, respectively). Interestingly, v_net switches sign in the κ = 5 model, while ∠β₂ varies only slightly, due to its link with the underlying field geometry. Overall, images with nonthermal eDFs will be useful to study in future work.

Zoom In Zoom Out Reset image size

With our GRMHD models, we coarsely sample a five-dimensional parameter space. Here we investigate the possibility that this sparse sampling misses potentially passing models by performing scoring using expanded theoretical error bars. We conceptualize each combination of a_*, R_high, and i as a volume in three-dimensional parameter space. For each neighbor in parameter space, if the 90% quantiles of the neighbor do not overlap, we linearly interpolate the lower and upper ranges of each observable to the midpoints of their nearest neighbors. This scheme helps mitigate sparse sampling but, as we discuss, may lead to false positives if observables evolve rapidly between adjacent models. In addition, this methodology fails to consider correlated evolution between observables.

In Figures 27 and 28, we show the results of our interpolative scoring scheme considering all polarimetric constraints without and with RM derotation, respectively. As expected, many more models pass in both cases. The preference for clockwise motion with derotation or counterclockwise motion without derotation is less dramatic with this scheme. Without derotation, both best-bet models still fail. With derotation, the second best-bet model from Paper V, MAD a_* = 0.5, R_high = 160, i = 30°/150°, also passes in this scheme. Without interpolation, this model had only failed by producing too little 〈∣m∣〉.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This interpolative scoring scheme does not produce as clear of a preference for MAD over SANE models. We find that this difference is driven by a shortcoming of this method: SANE models evolve very rapidly with R_high, especially between R_high = 1 and R_high = 10, leading to very large theoretical error bars. We explore one example in Figure 29, where a set of KHARMA SANE a_* = − 0.5, i = 150° aligned field models are ray-traced at intermediate values of R_high ∈ {3, 5, 8}. Each of our eight polarimetric observables is plotted, and we better resolve the rapid evolution in these parameters with R_high. A noteworthy interaction occurs in our interpolation scheme with 〈∣m∣〉 and ∠β₂, two of our most constraining observables. We see that at R_high = 1 the model overproduces 〈∣m∣〉 but fails to reproduce ∠β₂, which is too radial. Meanwhile, SANE models with R_high = 10 have too low 〈∣m∣〉 and a uniformly distributed ∠β₂. Interpolation allows models in this region to pass because our scoring system suggests that there might be a model with intermediate R_high that has both a correct ∠β₂ and a correct 〈∣m∣〉. However, with better resolution in R_high, we do not find an individual model that would pass. Overall, this exercise shows that our main conclusions are not likely driven by our sparse sampling of parameter space.

Zoom In Zoom Out Reset image size

To visualize trends of our eight observables in the five-dimensional parameter space that we explore, we provide "violin" plots of our observables from our models as a figure set, the complete version of which is available in the online journal. In each figure, we consider one observable and one magnetic field state (either MAD or SANE models). One figure, the distributions of m_net for MAD models, is shown in Figure 30. Different spins are shown in different columns, and different values of R_high are shown in different rows. Within each panel, we plot distributions as a function of inclination, where only five of the nine inclinations ray-traced in this work are included to improve readability. Aligned field models are shown on the left, and reversed field models are shown on the right. The distributions with opposite magnetic field polarity are usually very similar, with the notable exceptions of v_net and, more subtly, ∠β₂. To display the relative agreement or disagreement between codes, we plot BHAC models in red and KHARMA models in blue. H-AMR models, which are ray-traced for a subset of models only for comparison here and not for scoring, are displayed as dashed distributions when available. Finally, the observational constraints are shown in gray, where, as usual, the allowed range for ∠β₂ without RM derotation is shown as a hatched region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (18 panels)

Tarball of figure set images

Our last set of plots, distributions of the Faraday rotation depth $\langle {\tau }_{{\rho }_{V}}\rangle$, are not directly observable but drive many of our physical trends, as well as differences between codes. For a detailed discussion of the physical trends present in these figures, we refer readers to Appendix A.

Differences between our KHARMA and BHAC models inflate our theoretical error bars in Section 5. We find that at least part of these differences arise from physical approximations regarding the assignment of electron temperature during the GRRT. One fluid with a single adiabatic index is evolved in our GRMHD codes, but it represents both relativistic electrons (with an adiabatic index of 4/3) and nonrelativistic ions (with an adiabatic index of 5/3). During the GRRT step of our calculations, only the electron temperature is relevant for the synchrotron emission that we observe. When assigning electron temperatures, RAPTOR adopts (see, e.g., Davelaar et al. 2018)

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{e}=\displaystyle \frac{u}{\rho }\displaystyle \frac{{m}_{p}}{{m}_{e}}\displaystyle \frac{1}{3(R+1)},\end{eqnarray} \tag{ H1 }$$

where Θ_e is the electron temperature, u is the internal energy, and R = T_i /T_e given by Equation (8). Meanwhile, IPOLE accounts for the difference in adiabatic indices by adopting

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Theta }}}_{e} & = & \displaystyle \frac{u}{\rho }\displaystyle \frac{{m}_{p}}{{m}_{e}}\displaystyle \frac{({\gamma }_{p}-1)({\gamma }_{e}-1)}{({\gamma }_{e}-1)R+({\gamma }_{p}-1)}\\ & = & \displaystyle \frac{u}{\rho }\displaystyle \frac{{m}_{p}}{{m}_{e}}\displaystyle \frac{2}{3(2+R)},\end{array}\end{eqnarray} \tag{ H2 }$$

where γ_e = 4/3 and γ_p = 5/3. Equation (H2) is physically justified, but it sacrifices internal consistency with the GRMHD simulations, where a single fluid with γ = 4/3 is evolved (Wong et al. 2022). When we set γ_e = γ_p = γ = 4/3 in Equation (H2), we recover Equation (H1) used by RAPTOR. Electron temperatures assigned by RAPTOR are systematically colder, 3/4 as hot as the ipole prescription at R = 1, and 1/2 as hot as R → ∞ . This explains the systematically larger Faraday depths in our BHAC models relative to both KHARMA and H-AMR, which are both ray-traced with ipole.

Larger differences are seen between SANE models than MADs. A unique SANE model is not believed to exist, and differences are known to occur at the GRMHD fluid level (Porth et al. 2019).