Synchrotron Polarization Signatures of Surface Waves in Supermassive Black Hole Jets

To model the dynamics of the accretion flow around a Kerr black hole, we use the Black Hole Accretion Code (BHAC; Porth et al. 2017; Olivares et al. 2019), which solves the ideal GRMHD equations in curved spacetime. The equation of state is assumed to be an ideal gas law, described via the specific enthalpy

$$\begin{eqnarray}&&h(\rho ,{P}_{\mathrm{gas}})=1+\displaystyle \frac{{\gamma }_{\mathrm{ad}}}{{\gamma }_{\mathrm{ad}}-1}\displaystyle \frac{{P}_{\mathrm{gas}}}{\rho },\end{eqnarray} \tag{ 1 }$$

with gas pressure P_gas and mass density ρ in the fluid frame, and the adiabatic index is set to γ_ad = 13/9.

We initialize our simulation with a Fishbone & Moncrief (1976) torus in spherical Boyer–Lindquist coordinates (t, r, θ, ϕ), where θ is the angle with respect to the spin axis of the black hole, and ϕ is the azimuthal angle around the black hole spin axis. The initial conditions are then transformed to CKS coordinates (t, x, y, z). The initial disk has an inner radius of r_in = 20r_g, with r_g the gravitational radius given by r_g = GM/c², where G is Newton's constant, M is the black hole mass, and c is the speed of light. We set the pressure maximum of the disk at ${r}_{\max }=40\,{r}_{{\rm{g}}}$. The initial magnetic field is given by the vector potential ${A}_{\phi }=\rho {(r/{r}_{\mathrm{in}}\sin \theta )}^{3}{e}^{-r/400}-0.2$, where r is the radial coordinate, and θ is the polar angle. The vector potential follows isocontours of density ρ. These initial conditions are chosen such that the simulation will reach the saturated state of a MAD accretion flow.

We set the dimensionless black hole spin parameter to a = 15/16, which results in a black hole horizon size of r_h ≈ 1.34r_g. The domain size is (−4000r_g, 4000r_g) in all three spatial Cartesian directions. Our base resolution is 192³ cells. We employ nine additional levels of static mesh refinement, resulting in an effective uniform Cartesian resolution of 98,304³. The highest-resolution grid is centered at the horizon and has a resolution of Δx_i = 0.08r_g. The simulation is run for 10,000r_g/c. We introduce a maximum for the cold magnetization parameter $\sigma =\tfrac{{b}^{2}}{\rho }$, where b is the magnetic field strength, and we inject mass to maintain σ ≤ 100. Additionally, we ensure that ${\beta }_{{\rm{p}}}^{-1}\leqslant {10}^{3}$, where β_p = 8π P_gas/b² is the plasma beta parameter. We use floor profiles for density as well as pressure, given by ρ_floor = 10⁻⁴ r⁻² and ${P}_{\mathrm{floor}}={10}^{-6}{\rho }_{\mathrm{floor}}^{{\gamma }_{\mathrm{ad}}}$. Due to our Cartesian grid, we do not have an inner radius where we can employ outflowing boundary conditions; we therefore introduce an artificial treatment of the fluid variables inside the black hole event horizon, known as excision in numerical relativity, to limit the accumulation of energy and density, which otherwise could numerically diffuse out of the event horizon when accumulated to too-high values. In our case, we introduce a ceiling on density (${\rho }_{\max }=6$), as well as pressure (${P}_{\max }=2$) when r < r_crit, where our critical radius is set to be r_crit = 1r_g. Given the location of the critical radius, we have four cells between the event horizon and the critical radius.

To generate synthetic polarized images of the accretion flow, we use the general relativistic ray-tracing code RAPTOR. RAPTOR solves the polarized radiative transport equations along null geodesics. The geodesic equation is solved starting from a virtual camera outside the simulation domain (r_cam = 10⁴ r_g). We employ an adaptive camera as described in Davelaar & Haiman (2022). The adaptive camera allows a varying resolution over the image plane, resulting in a computational benefit, since most of the resolution can be focused on the event horizon scale, which shows small-scale emission varying on short timescales, while the larger scale can be fully resolved with a lower resolution. We use a base resolution of 625² pixels and double the resolution four times, within 60r_g, 40r_g, 20r_g, and 10r_g, bringing the effective resolution to 10,000² pixels. ¹³ We use an adaptive Runga–Kutta–Fehlberg method to integrate the geodesic equation and a fourth-order finite difference method to compute the metric derivatives needed for the Christoffel coefficients. The step size in RAPTOR is estimated based on the wavevector in the previous step; see Davelaar et al. (2019). For this work, we developed an additional step-size criterion based on a Courant–Friedrichs–Lewy condition, where we require a minimum of eight steps per cell to ensure convergence of all Stokes parameters. We compute synthetic images between 80 and 100 GHz with a frequency spacing of 3 GHz. We also compute time-averaged spectra from the image-integrated total intensity at 25 frequencies with a logarithmic spacing between 10¹⁰ and 10¹⁵ Hz.

The ideal GRMHD equations are dimensionless and do not describe the evolution and thermodynamics of electrons. To this end, we need to introduce a mass scaling and a prescription for the electron temperature. To scale our simulation to a specific black hole, we use a unit of length ${ \mathcal L }={r}_{{\rm{g}}}$, a unit of time ${ \mathcal T }={r}_{{\rm{g}}}/c$, and a unit of mass ${ \mathcal M }$. The unit of mass is related to the mass accretion rate via $\dot{M}={\dot{M}}_{\mathrm{sim}}{ \mathcal M }/{ \mathcal T }$, where ${\dot{M}}_{\mathrm{sim}}$ is the accretion rate in simulation units. Combinations of these units are then used to scale all relevant fluid quantities. The density is scaled by ${\rho }_{0}={ \mathcal M }/{{ \mathcal L }}^{3}$, the internal energy by u₀ = ρ₀ c², and the magnetic fields by ${B}_{0}\,=\,c\sqrt{4\pi {\rho }_{0}}$ (where B₀ is expressed in Gaussian units). Since the black hole mass is often constrained observationally, the only free parameter in our system is ${ \mathcal M }$, which can be used to set the energy budget of the simulation such that the total emission produced equals a user-set target. We follow Mościbrodzka et al. (2016) to parameterize the ratio between electron temperature T_e and proton temperature T_p via

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{ratio}} & = & \displaystyle \frac{1}{1+{\beta }_{{\rm{p}}}^{2}}+R\displaystyle \frac{{\beta }_{{\rm{p}}}^{2}}{1+{\beta }_{{\rm{p}}}^{2}},\\ {{\rm{\Theta }}}_{{\rm{e}}} & = & {k}_{{\rm{b}}}{T}_{{\rm{e}}}/{m}_{{\rm{e}}}{c}^{2}=\displaystyle \frac{U({\gamma }_{\mathrm{ad}}-1){m}_{{\rm{p}}}/{m}_{{\rm{e}}}}{\rho \left(1+{T}_{\mathrm{ratio}}\right)},\end{array}\end{eqnarray} \tag{ 2 }$$

where U is the internal energy, m_p is the proton mass, m_e is the electron mass, and Θ_e is the dimensionless electron temperature. The variable R is a free parameter that sets the temperature ratio in regions where β_p ≫ 1, allowing for different emission morphology depending on the choice of R, e.g., disk-dominated if R = 1 or jet-dominated if R ≫ 1. Here, we will limit ourselves to R = 20, e.g., a more jet-dominated model. Note that MADs are relatively insensitive to the exact value of R given that most of the emission originates from a region where β_p ≲ 1. Finally, we must choose the electron DF's shape and spatial variation. We consider two models, one where the DF is limited to a thermal relativistic Maxwell–Jüttner DF (MJ-DF), and one where we combine the MJ-DF with a κ-DF. The κ-DF deviates from an MJ-DF by having a thermal core and a power law at high Lorentz factors. The κ-DF (Xiao 2006) is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{e}}{d\gamma }={n}_{e}N\gamma \sqrt{{\gamma }^{2}-1}{\left(1+\displaystyle \frac{\gamma -1}{\kappa w}\right)}^{-(\kappa +1)},\end{eqnarray} \tag{ 3 }$$

where the free parameters are w, which sets the width of the distribution, and κ, which sets the slope of the power law, and N is a normalization constant such that ${\int }_{1}^{\infty }\tfrac{{{dn}}_{e}}{d\gamma }d\gamma ={n}_{{\rm{e}}}$. The κ parameter is related to the power-law index p of the nonthermal tail of the DF via κ = p + 1, such that for γ ≫ 1, $\tfrac{{{dn}}_{e}}{d\gamma }\propto {\gamma }^{1-\kappa }$. For the width w, we follow Davelaar et al. (2019) by enforcing that the energy in the DF equals the total available internal energy of the electrons,

$$\begin{eqnarray}&&w=\displaystyle \frac{(\kappa -3)}{\kappa }{{\rm{\Theta }}}_{e},\end{eqnarray} \tag{ 4 }$$

where Θ_e is computed with Equation (2). Note that this formula requires κ > 3 (p > 2).

We then compute the emission coefficients c_ν (emission, absorption, and rotation coefficients) using a prescription introduced in Event Horizon Telescope Collaboration et al. (2022) that combines thermal and κ coefficients via

$$\begin{eqnarray}&&{c}_{\nu }=(1-\eta ({\beta }_{{\rm{p}}},\sigma )){c}_{\mathrm{thermal}}+\eta ({\beta }_{{\rm{p}}},\sigma ){c}_{\kappa },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\eta ({\beta }_{{\rm{p}}},\sigma )=\left(1-{e}^{-{\beta }_{{\rm{p}}}^{-2}}\right)\left(1-{e}^{-{\left(\sigma /{\sigma }_{0}\right)}^{2}}\right),\end{eqnarray} \tag{ 6 }$$

where σ₀ sets the transition point for the magnetization, which we set to σ₀ = 0.5. The function η(β_p, σ) smoothly transitions from thermal to nonthermal components if β_p < 1 and σ/σ₀ > 1, representing sites with a large reservoir of magnetic energy that can dissipate to accelerate particles, e.g., in the jet's shear layer. The polarized radiation transfer coefficients are computed via a fit formula. For the thermal DF, we use the emission and absorption coefficients presented in Dexter (2016) and the rotativities from Shcherbakov (2008); for the κ-DF, we follow Pandya et al. (2016) and Marszewski et al. (2021).

As an archetype LLAGN, we use model parameters consistent with M87*, meaning a black hole mass of M = 6.5 × 10⁹ M_⊙ and a distance of d = 16.8 Mpc (Event Horizon Telescope Collaboration et al. 2019a). We set the angle between the BH spin axis and the observer to i = 160°, following Walker et al. (2018). We set ${ \mathcal M }$ such that the average flux at 86 GHz is F_86GHz = 1 Jy, as observed by EHT MWL Science Working Group et al. (2021). The spectrum obtained by EHT MWL Science Working Group et al. (2021) also shows a spectral slope in the optically thin part of the spectrum in the near-infrared (NIR) of F_ν ∝ ν⁻¹. Given that for optically thin synchrotron emission, the flux follows F ∝ ν^−(p−1)/2 = ν^{−(κ−2)/2}, to match the observed spectral slope in the NIR, we need to use κ = 4.

To exclude the emission from the spine (interior of the jet, which in GRMHD simulations is typically dominated by artificial density floors), we exclude all of the emission if σ > 5. We expect our results to be insensitive to this choice for larger σ values. However, lower σ values would result in a smaller emission region at the jet–wind interface. We also exclude the larger-scale disk, $\sqrt{{x}^{2}\,+\,{y}^{2}}\gt 150{r}_{{\rm{g}}}$, which has not settled into a steady state for the runtime of our simulation. However, given the viewing angle of i = 160°, this choice does not strongly affect our results.

To assess when the simulation reaches the MAD state, we compute at the black hole's event horizon the mass accretion rate, $\dot{m}$, and magnetic flux threading the horizon, Φ_B, both in dimensionless units. The mass accretion rate is defined as $\dot{m}={\int }_{r={r}_{{\rm{h}}}}\rho {u}^{r}\sqrt{g}d\theta d\phi$, where u^r is the radial component of the velocity, and $\sqrt{g}$ is the determinant of the metric. The magnetic flux is defined as ${{\rm{\Phi }}}_{B}=0.5{\int }_{r={r}_{{\rm{h}}}}| {B}^{r}| \sqrt{\gamma }d\theta d\phi$, where B^r is the radial component of the magnetic field, and $\sqrt{\gamma }$ is the determinant of the spatial part of the metric. ¹⁴ We also define the MAD parameter ${\phi }_{\mathrm{mad}}={{\rm{\Phi }}}_{{\rm{B}}}/\sqrt{\dot{m}}$, which was introduced by Tchekhovskoy et al. (2011) to quantify that a simulation reaches the MAD state when ϕ_mad ∼ 15. All three quantities, $\dot{m}$, Φ_B, and ϕ_mad, are shown in Figure 1. The accretion flow reaches the MAD state for the first time at ≈3000r_g/c, when the magnetic flux on the horizon saturates as ϕ_mad ∼ 15. Our MAD simulation shows globally similar properties to standard simulations on spherical grids; see Appendix A for a comparison.

Zoom In Zoom Out Reset image size

In Figure 2, we show two-dimensional maps of the logarithm of density sliced along the spin axis (top row) or the equatorial plane (bottom row). Initially, the jet is more laminar in the left column, as flux is still building up on the horizon, and no eruptions have occurred. In the middle column, after the system has reached the MAD limit, flux tubes can be seen in the x–y plane, indicated by lower densities in the disk. The exhaust from magnetic flux eruptions generates these flux tubes. During the eruptions, magnetic energy is dissipated via large equatorial current sheets generated when the disk becomes magnetically arrested, and the northern and southern jets get in direct contact (Ripperda et al. 2022). The exhaust of these sheets forms flux tubes containing vertical magnetic fields that spiral outward in the disk before dissipating due to Rayleigh–Taylor mixing (Zhdankin et al. 2023).

Zoom In Zoom Out Reset image size

In the middle panel, small-amplitude waves propagate along the shear layer, interfacing the higher-density disk wind and low-density jet in the top panel. Large-amplitude waves propagate outward in the right column because the system produces strong magnetic flux eruptions; see Figure 1 at t = 8800r_g/c. The panels correspond to the vertical lines in the middle panel of Figure 1. The variability introduced by the flux eruptions at the base of the jet acts like a forced oscillator, introducing waves that propagate and grow from near the event horizon to the shear layer between the jet and the disk. The waves propagate to large scales, growing in size while shearing magnetic field lines and generating field reversals, shown in Figure 3.

Zoom In Zoom Out Reset image size

To assess the potential effect of the waves on the polarized emission, we show slices along the x – z plane in Figure 4 at t = 8850r_g/c (right panels of Figure 2) of several quantities relevant to the radiation transport. We find that distinct patches of magnetization close to σ ∼ 1 can be seen along the jet–wind shear layer coinciding with the location of the waves in the top right panel of Figure 2. The higher magnetization is important, since particle acceleration is more efficient at higher magnetization values (Sironi & Spitkovsky 2014). We use a passively advected tracer to study the acceleration of wind-based material to relativistic speeds at the jet–wind surface. The tracer is set to 0 when the density or pressure is set to floors and set to 1 in the initial torus. We then evolve this quantity as a passive scalar, tracing the advection of disk-based material into regions that were at least once set to floors. We find that matter from the disk (unfloored material) around the location of the jet–wind shear layer waves is accelerated to a high bulk Lorentz factor, and emission generated in the waves is emitted from unfloored matter originating from the accretion disk, shown by a nonzero tracer and Γ ≥ 2 in the top right panel of Figure 4. Looking at the bottom left panel, the waves also correlate with regions with a large fraction of nonthermal electrons since η ∼ 1, which is an evident result of the high magnetization (and low plasma-β) and our choice of electron distribution model (Equation (6)). Finally, the electrons are relativistically hot (i.e., Θ_e > 1), as shown in the bottom right panel. We associate the relativistic temperatures and the large fraction of nonthermal electrons with the heating of the plasma by the waves due to the dissipation of magnetic energy, as predicted in Sironi et al. (2021). In summary, the shear layer that is at the surface of the jet interfacing with the wind has a magnetization of order unity, has a high relativistic electron temperature, is moving at relativistic bulk Lorentz factors (Γ ∼ 3), and is likely dominated by nonthermal electrons.

Zoom In Zoom Out Reset image size

In Figure 5, we show the spectral DFs of our thermal-DF (thermal-jet) and κ-DF (κ-jet) models. The top panel shows the total intensity (Stokes ${ \mathcal I }$) as a function of frequency. The thermal-jet model recovers the low-frequency part of the spectrum accurately, e.g., ν < 10¹² GHz; however, at higher frequencies, it drops off too fast, which is consistent with Davelaar et al. (2019). In the case of the κ-jet model, the high-frequency emission is enhanced and obtains a spectral slope of α ≈ −1, consistent with the observations. In contrast, the thermal-jet model underestimates the NIR flux. The κ-jet model predicts that nonthermal electrons emit energetic photons predominantly in the jet boundary and are, therefore, a probe of dissipation of magnetic energy due to wave dynamics. To match the observed flux at 86 GHz, F_86GHz = 1Jy, we set the mass scaling to ${ \mathcal M }=1.5\,\times \,{10}^{25}$ g for both the thermal and κ-jet models. These units of mass correspond to a mass accretion rate of ≈10⁻⁴ M_⊙ yr⁻¹ and a jet power of ≈10^42–43 erg s⁻¹, both similar to values obtained in previous works (Collaboration et al. 2019b; Chael et al. 2019; Cruz-Osorio et al. 2022) and consistent with the jet powers inferred from observations of M87 (Prieto et al. 2016).

Zoom In Zoom Out Reset image size

The two bottom panels show LP and CP. For LP, the thermal jet achieves similar fluxes as the κ jet at a lower frequency (ν ≲ 10¹³ Hz), while at a higher frequency, a clear power law is visible in the κ-jet case. A similar power law is visible for CP, but at a lower frequency, the κ-jet model is comparable to the thermal case. Subsequent analysis is done at 86 GHz. This frequency was chosen because it probes emission structures at the base of the jet (Hada et al. 2016; Walker et al. 2018).

In the first row of Figure 6, we show total intensity maps for our κ-jet model at 86 GHz. The images correspond to t = 8500r_g/c and 9300r_g/c. At the core of the image, a darkening is visible, corresponding to the "black hole shadow" (Luminet 1979; Falcke et al. 2000) recently observed at 86 GHz by Lu et al. (2023). At the later stage (right panel), when the surface waves are prominently visible in the GRMHD simulations, helical wavelike substructures can be distinguished within the jet at larger scales.

Zoom In Zoom Out Reset image size

This subsection summarizes the LP results computed at 86 GHz. In unresolved LP fractions, we find that enhanced emission regions generate loops in a Stokes ${ \mathcal Q }-{ \mathcal U }$ diagram during the magnetic flux eruptions. The resolved LP fraction is then inversely proportional to the magnetic flux on the horizon, resulting in an enhancement of the fraction as the flux goes down. Furthermore, we find that the waves seen in the GRMHD simulation imprint themselves in LP maps as they lower the LP fraction. The effect of nonthermal κ-DF on the LP is minor, although we see a slight decrease compared to the thermal model for the LP fraction in the jet.

3.4.1. Core

In Figure 7 (top panel), we show a time series of the unresolved LP fraction m_net for both the thermal-jet (red lines) and κ-jet (blue lines) models, defined as

$$\begin{eqnarray}&&{m}_{\mathrm{net}}=\displaystyle \frac{\sqrt{{\left({\sum }_{\mathrm{pixels}}{S}_{Q}\right)}^{2}+{\left({\sum }_{\mathrm{pixels}}{S}_{U}\right)}^{2}}}{{\sum }_{\mathrm{pixels}}{S}_{I}}.\end{eqnarray} \tag{ 7 }$$

The image-integrated net LP fraction does not depend on telescope beam size, since it is an incoherent addition of the Stokes parameters. We obtain an average value of m_net = 0.026, consistent with the low values found by observations (Hada et al. 2016; Walker et al. 2018). The thermal and κ-jet models show almost identical values.

Zoom In Zoom Out Reset image size

In Figure 7 (bottom panel), we show the resolved LP fraction, m, defined as

$$\begin{eqnarray}&&m=\displaystyle \frac{{\sum }_{\mathrm{pixels}}\sqrt{{S}_{Q}^{2}+{S}_{U}^{2}}}{{\sum }_{\mathrm{pixels}}{S}_{I}}.\end{eqnarray} \tag{ 8 }$$

Here we follow the definition of resolved LP fraction from Equation (8) of Event Horizon Telescope Collaboration et al. (2021), which is an image-averaged LP fraction, taking into account some telescope beam size. As a first step, we assumed that the telescope fully resolves the image, taking the beam size to be much smaller than the intrinsic features we are interested in. Due to the coherent addition, the resulting resolved LP fraction m is substantially higher than the unresolved LP fraction m_net. The reason for this is that we preserve the sign of ${ \mathcal Q }$ and ${ \mathcal U }$ in the summation in the case of m_net, but the sign is dropped for m. In reality, ${ \mathcal Q }$ and ${ \mathcal U }$ will be convolved with the telescope's beam, resulting in incoherent addition. This effect can be seen in the bottom panel of Figure 7, where we compute the convolved LP fraction m_conv for varying telescope beam sizes of 20, 40, 80, and 160 μas. This computation is done by blurring the original images with a Gaussian filter where the full-width at half maximum represents the beam size. As the telescope beam size increases, the underlying substructure in both ${ \mathcal Q }$ and ${ \mathcal U }$ is being averaged out, resulting in an overall drop in LP fraction asymptotically approaching m_net as the beam becomes comparable to the core size.

In the fully resolved LP fraction m, as well as the 20–80 μas m_conv cases in Figure 7, at t = 9000r_g/c, an increase in LP fraction is visible. This increase coincides with a flux eruption at the horizon, visible from the overplotted Φ_B curve (solid black line; identical to the middle panel of Figure 1). The m shows a fractional increase of 20%, while Φ_B shows a fractional decrease of 20%, indicating an inverse relationship between the two quantities. A similar correlation can also be seen at the smaller eruption at t = 8400r_g/c. The κ-jet model reaches identical LP fractions as the thermal jet.

Close to the flux eruption starting around t ∼ 8700r_g/c, we show the unresolved Stokes parameters in a ${ \mathcal Q }/{ \mathcal I }-{ \mathcal U }/{ \mathcal I }$ diagram. After the onset of flux eruption, we see a clear clockwise loop with an LP excess of m_net ∼ 0.06 (Figure 8). The loop we find is similar to the loops found by Marrone et al. (2008) and Wielgus et al. (2022) at 230 GHz during X-ray flares for Sagittarius A*. Najafi-Ziyazi et al. (2023) finds evidence that these loops could be generated by enhanced polarized emission in the accretion disk by orbiting flux bundles ejected into the disk after a flux eruption. The enhanced emission leads to a local polarization excess; as the emission increases, the polarized emission increases. This enhanced emission region, often called a hot spot, orbits through the local magnetic field. Since the spot only lights up a small region of the accretion disk, it acts as a probe for the underlying magnetic field geometry. The magnetic field geometry close to the jet base is mostly poloidal; therefore, Stokes ${ \mathcal Q }$ and ${ \mathcal U }$ generate four quadrants, like a spoke wheel pattern, in the image plane, which alternate in sign; see, e.g., Narayan et al. (2021) and The GRAVITY Collaboration et al. (2023). Since the magnetic field orientation periodically varies, Stokes ${ \mathcal Q }$ and ${ \mathcal U }$ will show sinusoidal behavior in the total integrated ${ \mathcal Q }$ and ${ \mathcal U }$, since an excess of either positive or negative ${ \mathcal Q }$ and ${ \mathcal U }$ is seen; for more details, see Vos et al. (2022) and Najafi-Ziyazi et al. (2023).

Zoom In Zoom Out Reset image size

3.4.2. Jet

To study the large-scale jet only, we exclude the image plane's inner 30r_g; in other words, we exclude the near-horizon emission. The resulting m_net and m are shown in Figure 9. The unresolved LP fraction m_net ranges from 5% to 20%. The resolved fraction, m, reaches values of around 50%. The κ-jet model shows slightly lower values than the thermal-jet model in the case of m. In the second row of Figure 6, we show a map of m. In these maps, alternating regions of high and low LP are visible, coinciding with enhanced emission in the total intensity panels in the first row. This indicates a potential correlation between the LP fraction and the presence of the waves seen in the GRMHD simulation. This potential correlation will be investigated further in Section 3.8. The LP substructure seen in our simulation in Cartesian coordinates is absent in a low-resolution MAD simulation in spherical coordinates at an effective resolution of 192 × 96 × 96 cells in r, θ, and ϕ, respectively; see Appendix B.

Zoom In Zoom Out Reset image size

The jet stands out as a high-intensity emission region, while the accretion disk does not contribute to the emission at larger scales (r > 30r_g). Comparing the total intensity map with the LP fraction map, the foreground disk surrounding the jet shows high values of m (close to unity). To understand this behavior, we evaluate the asymptotic limit of our fitting formula for the emission coefficients ${J}_{{ \mathcal S }}$ (with ${ \mathcal S }$ indicating one of the Stokes parameters; see Equation (31) in Pandya et al. 2016). Given that the disk has weak magnetic fields and low temperatures, we have ν/ν_c ≪ 1, with ν_c = eB/(2π m_e c²), and Θ_e ≪ 1, and we find that ${J}_{{ \mathcal I }}/| {J}_{{ \mathcal Q }}| \to 1.0$. This makes physical sense because, due to the low temperature, the thermal DF is narrow, which means that there is a quasi-monoenergetic population of electrons that is causing the emission, so the polarization fraction should go to unity. This, however, does not alter the computation of our image-integrated m, m_conv, and m_net, since this outer region has low intensity, both polarized and Stokes ${ \mathcal I }$, so they do not contribute to the numerator and denominator of Equations (7) and (8). Additionally, we also exclude regions of very low intensity from our map, where we set the Stokes parameters of a pixel to 0 if ${S}_{{ \mathcal I }}/\max ({S}_{{ \mathcal I }})\lt {10}^{-6}$.

Lastly, we compute polarization angle maps as a function of beamwidth. These can be seen in Figure 10, where we overplotted ticks of the polarization vector on top of the total intensity map from the top left panel of Figure 6. The length of the ticks is set by the LP fraction, while the angle is set by the electric vector position angle (EVPA), as defined by $\chi =\tfrac{1}{2}{\tan }^{-1}({S}_{{ \mathcal Q }}/{S}_{{ \mathcal U }})$; here we also exclude LP fractions when ${S}_{{ \mathcal I }}/\max ({S}_{{ \mathcal I }})\lt {10}^{-6}$, so the tick lengths are set to 0. For the case of zero beam size, the polarization vector clearly follows the ridges of enhanced intensity. As the beam size increases, the correlation becomes weaker. However, reorientation of the polarization vector is also visible here; e.g., at x ∼ 0r_g and y ∼ 50r_g, the polarization pattern switches from diagonal, to horizontal, and back to diagonal. We only show beam sizes up to 60 μas, which is the expected resolution of the next-generation EHT (ng-EHT) at 86 GHz (Issaoun et al. 2023), assuming identical baselines to the current EHT array (using θ ∼ 20 μas (λ/1mm) with λ = 3 mm).

Zoom In Zoom Out Reset image size

In this subsection, we study the CP fractions computed at 86 GHz. We find that the resolved CP fraction decreases during a flux eruption as the inner accretion disk is ejected, resulting in a more dilute plasma to perform Faraday conversion. Overall, we find low unresolved and resolved CP fractions for our thermal and κ-jet models. The CP maps show sign reversals, indicative of an alternating magnetic field orientation that coincides with the features seen in the LP maps.

3.5.1. Core

In Figure 11, we show the unresolved and resolved CP fractions, v_net and v, defined as

$$\begin{eqnarray}&&{v}_{\mathrm{net}}=\displaystyle \frac{\sqrt{{\left({\sum }_{\mathrm{pixels}}{S}_{V}\right)}^{2}}}{{\sum }_{\mathrm{pixels}}{S}_{I}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&v=\displaystyle \frac{{\sum }_{\mathrm{pixels}}\sqrt{{S}_{V}^{2}}}{{\sum }_{\mathrm{pixels}}{S}_{I}}.\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

Both v_net and v are small in value, as expected from synchrotron radiation (Rybicki & Lightman 1979). In Figure 11 (bottom panel), during the strongest flux eruption at t = 9000r_g/c, both models show a slight decrease in CP fraction. The accretion disk enhances the amount of CP emission as Faraday conversion converts LP to CP. Due to the ejection of the inner part of the accretion disk during a flux eruption, there is more dilute plasma in this region, resulting in a drop in the CP fraction.

3.5.2. Jet

In Figure 12, we compute v_net and v but now also exclude the inner 30r_g to exclude the near-horizon emission and focus on the larger-scale jet only. This exclusion results in smaller fractions than the entire image-integrated values, since Faraday conversion happens in high-density regions. The third row of Figure 6 shows maps of v where the CP fraction is smaller at larger radii. In the map, we preserve the sign of Stokes ${ \mathcal V }$, which is set by the direction of the magnetic field along the line of sight. The image shows reversals in the sign of v and additional ridges of low CP fraction. Since Stokes ${ \mathcal V }$ carries information on the direction of the magnetic field along the line of sight, this indicates a potential orientation switch in the underlying magnetic field geometry. This reversal could be caused by the waves, as shown in Figure 3; we will further investigate the sign reversal in Section 3.8.

Zoom In Zoom Out Reset image size

Figure 6 (fourth row) shows maps of RMs computed between 80 and 100 GHz. The RM is defined as

$$\begin{eqnarray}&&\mathrm{RM}=\displaystyle \frac{{\chi }_{1}-{\chi }_{2}}{{\lambda }_{1}^{2}-{\lambda }_{2}^{2}},\end{eqnarray} \tag{ 11 }$$

where χ_ν is the EVPA computed at a specific frequency ν defined as $\chi =\tfrac{1}{2}{\tan }^{-1}({S}_{{ \mathcal Q }}/{S}_{{ \mathcal U }})$, and λ is the wavelength. Large RMs are visible along the jet, and in the core, we find values as large as 10⁴–10⁵ rad m^–2. However, the core dominates the total RM, which is expected, since the Faraday depth is larger due to higher density and lower temperatures. The relatively large value for the RM in the jet is somewhat surprising, given that the jet does not exhibit large Faraday depths. Given the small Faraday depth, the change in EVPA is caused not by Faraday rotation but rather by the transverse gradients in n_e, Θ_e, and B in the emitting shear layer. These gradients result in emission at different frequencies peaking at different depths in the shear layer, which have different plasma properties; e.g., different magnetic field orientations will result in different orientations of the EVPA, giving rise to nonzero RM values.

To test if the waves cause linear depolarization, we identified representative light rays showing high or low polarization fractions. The selected geodesics are indicated with the red, blue, and green dots in the top panel of Figure 13. We then compute the LP fraction as a function of the CKS $x^{\prime}$ coordinate, meaning smaller values of $x^{\prime}$ are closer to the spin axis. The result of this is shown in the bottom panel of Figure 13. We show the geodesics only as they approach the jet–wind surface and only show segments when the local magnetization σ < 5, meaning no radiation transfer is applied when the geodesic is inside the jet. The red ray is terminated early, meaning that for larger $x^{\prime}$, its net polarization fraction is higher, indicating that the shear layer at that point is thinner. The blue and green rays have a larger travel path, meaning that the polarization starts to drop. In the case of the green ray, the situation is even more interesting. The line of this ray is interrupted twice, indicating that it crossed into two regions of high magnetization but then left. We interpret this as the ray crossing through a rolling wave, similar to the wave seen at y = −100r_g in Figure 2. We test this by also computing the magnetization along the ray; this is also shown in the bottom panel of Figure 13 by a black line. This line crosses our magnetization threshold (σ > 5) twice. In general, the waves alter the thickness of the jet–wind shear layer, and their presence results in varying path lengths inside the shear layer for different rays, which affects the LP fraction.

Zoom In Zoom Out Reset image size

To connect the waves in GRMHD with the structures seen in the synthetic 86 GHz images, we compute emissivity-weighted averages of the magnetization σ, the pitch angle between the wavevector and the magnetic field orientation $\cos ({\theta }_{{\rm{B}}})$, the electron number density n_e, and the electron temperature Θ_e via

$$\begin{eqnarray}&&\langle q\rangle =\displaystyle \frac{\int {j}_{\nu }{qd}{\lambda }_{\mathrm{aff}}}{\int {j}_{\nu }d{\lambda }_{\mathrm{aff}}},\end{eqnarray} \tag{ 12 }$$

where q represents the weighted quantity; the result of this computation is shown in Figure 14. The top left panel, σ, shows that the images' higher-intensity features also have a larger magnetization. This agrees with the waves having a larger magnetization in Figure 4. Additionally, the same patterns are visible in the higher electron temperatures (top right panel of Figure 14) and the overdensities (bottom left panel of Figure 14), also in agreement with the properties of the waves, as shown in Section 3.1.

Zoom In Zoom Out Reset image size

For Stokes ${ \mathcal V }$, we finally compare the emissivity-weighted average of $\cos ({\theta }_{{\rm{B}}})$ (bottom right panel of Figure 14). The overall map of $\cos ({\theta }_{{\rm{B}}})$, where θ_B is the angle between the wavevector and the magnetic field vector, shows the same sign as Stokes ${ \mathcal V }$, implying that reversals in Stokes ${ \mathcal V }$ are caused by the shearing of the magnetic fields leading to reversed field orientations.