Balanced Turbulence and the Helicity Barrier in Black Hole Accretion

Horizon-scale observations from the Event Horizon Telescope (EHT) have enabled precision study of supermassive black hole accretion. Contemporary accretion modeling often treats the inflowing plasma as a single, thermal fluid, but microphysical kinetic effects can lead to significant deviations from this idealized picture. We investigate how the helicity barrier influences EHT-accessible electromagnetic observables by employing a simple model for electron heating based on kinetic physics and the cascade of energy and helicity in unbalanced turbulence. Although the helicity barrier plays only a minor role in regions with high plasma β, like in standard and normal evolution (SANE) disks, it may substantially impact regions with more ordered magnetic fields, such as the jet and its surrounding wind in SANE flows as well as throughout the entire domain in magnetically arrested disk (MAD) flows. In SANE flows, emission shifts from the funnel wall toward the lower-magnetization disk region; in MAD flows the emission morphology remains largely unchanged. Including the helicity barrier leads to characteristically lower electron temperatures, and neglecting it can lead to underestimated accretion rates and inferred jet powers. The corresponding higher plasma densities result in increased depolarization and Faraday depths thereby decreasing the amplitude of the β 2 coefficient while leaving its angle unchanged. Both the increased jet power and lower β2 may help alleviate outstanding tensions between modeling and EHT observations. We also find that the estimated ring diameter may be underestimated when the helicity barrier is neglected. Our results underscore the significance of the helicity barrier in shaping black hole observables and inferred accretion system parameters.


INTRODUCTION
Low luminosity active galactic nuclei (LLAGN) are usually modeled as radiatively inefficient accretion flows (RIAFs) onto supermassive black holes (e.g., Ichimaru 1977;Rees et al. 1982;Narayan & Yi 1995;Reynolds et al. 1996).In contrast to radiatively efficient thin disks model, RIAFs comprise geometrically thick, optically thin disks of hot plasma that circle the hole at subkeplerian speeds; in RIAFs, the gravitational binding energy of the inflowing plasma is converted into heat that cannot be radiated away before the plasma accretes down to the horizon.The excess heat provides a thermal pressure that supports a puffy disk.
Two of the LLAGN with the largest known sizes on the sky lie at the centers of our galaxy (Sgr A*) and the nearby elliptical galaxy Messier 87 (M87*).These two presumed RIAF sources are large enough to be directly observed by the Event Horizon Telescope (EHT), and the EHT's very long baseline interferometric experiment has produced radio images of the horizon-scale emission in spectacular detail (Doeleman et al. 2012;Akiyama et al. 2015Akiyama et al. , 2017;;Event Horizon Telescope Collaboration et al. 2019a, 2021a, 2022a).The observations can be used to probe plasma physics in these extreme environments while also providing constraints on key physical parameters of the accretion system, like the black hole mass and spin, the system accretion rate, and the amount of magnetic flux trapped on the horizon (Event Horizon Telescope Collaboration et al. 2019b, 2021b, 2022b).Ongoing measurements and the next generation of these experiments will help inform further parameter constraints and enable precision tests of our understanding of the physics governing these systems.The first observational results have already revealed tension between models and data in quantities like the resolved linear polarization fraction, which magnetized models often overproduce, and the predicted jet powers, which lie almost categorically at the lower end of the observational bounds.
Parameter estimation infers that M87* and Sgr A* are likely Coulomb collisionless since the path length to a Coulomb interaction greatly exceeds the size of the system.The electrons and ions that make up the infalling plasma therefore do not have time to equilibrate and relax to a thermal Maxwell-Jüttner distribution (e.g., Shapiro et al. 1976;Mahadevan & Quataert 1997;Quataert 1998;Sądowski et al. 2017;Chael et al. 2018;Ryan et al. 2018).Nevertheless, it may be that the ion-and electron-distribution functions are independently thermal, since intraspecies interactions due to kinetic plasma instabilities can drive particle-wave interactions that enable relaxation (see Kunz et al. 2014 and discussion therein).
The mechanisms that govern the heating and cooling of the ions and electrons are the subject of detailed study.Since the radio emission observed from Sgr A* and M87* is produced by the synchrotron process, the distribution function of the electrons plays a crucial role in determining the observational features of the sources.Accurately modeling particle acceleration is thus essential, since turbulent heating, reconnection, and shocks all yield different heating profiles and it is likely that different combinations of all heating mechanisms operate in different parts of the accretion flow.Interpreting the nonthermal features of the observations will also require a detailed understanding of the processes that determine the local particle distribution functions.
RIAFs have long been modeled with semianalytic (e.g., Ichimaru 1977;Rees et al. 1982;Narayan & Yi 1995;Özel et al. 2000) and numerical (e.g., Hawley 2000;De Villiers & Hawley 2003;McKinney & Gammie 2004) methods.The latter models are usually produced through general relativistic magnetohydrodynamics (GRMHD) simulation and are often favored over semianalytic models because they naturally incorporate properties of the turbulent dynamics, produce variability, and effectuate the connection between the accretion disk, wind, and jet, all of which may play an important role in determining the observational appearance of the system.The output of the fluid simulations is typically processed through general relativistic ray tracing (GRRT) codes to generate simulated observables like images and spectra.
In the standard modeling procedure, only the total energy of the fluid is tracked and evolved in the simulation.Since the electron distribution function is required to compute the radiative transfer coefficients, the typical modeling approach is to assume that the electrons are thermal and assign their temperature in a post-processing step by partitioning the total internal energy of the fluid into the ions and electrons following a prescription that depends on the ratio of the gas-to-magnetic pressure and the magnetization.These thermodynamic prescriptions are usually motivated by (kinetic) plasma theory for turbulent cascades, magnetic reconnection, collisionless shocks, and so on. 1 This approach introduces significant uncertainty and may well explain the model/data tension: the presence of a population of cold electrons would require an increased mass accretion rate, higher jet power, and would result in more depolarization from Faraday scrambling.
The turbulent cascade model is often invoked to quantify the energy partition: energy is injected at large scales (e.g., from the magnetorotational instability or large-scale torques) and cascades to higher wavenumbers until it is dissipated as thermal energy into the ions or electrons at their associated Larmor scales.But when the turbulence is imbalanced, not all of the energy injected at large scales can be treated the same way, and if plasma β ≡ P gas /P mag is small, conservation of helicity can inhibit energy flow in the cascade.
In low-β plasmas, the sense of the helicity cascade above and below the ion Larmor scale changes: as wavenumber increases into the kinetic scale, the fluid cross-helicity transforms conservatively into magnetic helicity and the direction of the helicity cascade inverts (Biskamp 2003;Meyrand et al. 2021;Squire et al. 2022).Since (generalized) helicity is conserved, the helicityendowed component of the turbulence is then unable to cascade below the ion Larmor scale and an effective helicity barrier is produced, limiting the fraction of the energy at large scales that can reach and heat the electrons.Observational evidence for the helicity barrier's operation has recently been found in the context of the solar wind via correlation of ion-cyclotron waves and electron-scale turbulence (Bowen et al. 2023).
In this paper, we use a simple model to study the effect of the helicity barrier on electron heating in radiatively inefficient accretion flows and probe its effect on the electromagnetic observables accessible to the horizon-scale radio observations.Since the presence of long-lived imbalanced turbulence is required for the helicity barrier to operate, one might expect the helicity barrier to be most important in the directed winds above the surface of the disk, but it is challenging to generate a predictive model for the quantitative details.In this study, we consider a limited subset of models and perform a preliminary study of the effect of the helicity barrier.We thus aim only to identify and describe broad qualitative trends to judge the importance of the effect and whether it may help explain contemporary questions raised by the data.We leave a more detailed study to future work.
In Section 2 we provide a brief overview of the problem of particle heating and its connection to the helicity barrier.We describe the details of our numerical methods in Section 3. In Section 4, we present the results of our numerical exploration, and we discuss model assumptions and limitations in Section 5. We conclude in Section 6.

PARTICLE HEATING IN COLLISIONLESS TURBULENCE
Robust interpretation of black hole images requires a detailed understanding of how emission is produced by the accretion flow.Most of the emission is due to synchrotron radiation (e.g., Yuan & Narayan 2014) and thus it is very sensitive to the electron momentum distribution function.The electron distribution function is determined by the details of the particle heating mechanisms that transform the gravitational energy released during accretion into thermal kinetic energy.The channels responsible for this conversion include dissipation in shocks (e.g., Tidman & Krall 1971;Blandford & Eichler 1987;Ghavamian et al. 2007;Mondal & Basu 2020;Tran & Sironi 2020), at reconnection sites (e.g., Bisnovatyi-Kogan & Lovelace 1997;Rowan et al. 2017;Werner et al. 2018;Rowan et al. 2019;Sironi & Beloborodov 2020;Werner & Uzdensky 2021), and in turbulent cascades.The relative importance of these channels is not well understood.In this paper, we focus on the turbulent cascade as the main source of energy for electrons.

The turbulent cascade
The conventional picture of the turbulent cascade involves a specified outer injection scale at which energy is supplied by large-scale processes (e.g., the MRI or largescale torques) and specified smaller dissipation scales at which the energy transforms into unordered kinetic motion (e.g., plasma kinetic scales or viscous/resistive scales in collisional systems).Solutions that bridge between these scales must conserve energy flux and are often assumed to only include interactions that are local in scale (e.g., only the eddies of similar sizes can efficiently interact with each other).The large separation between injection and dissipation scales often leads to assumption of "zeroth law of turbulence," which states that the large-scale behavior of the cascade does not depend on the physics responsible for its dissipation.The assumption of scale-independence has been very useful in constructing models for collisional turbulence (Kolmogorov 1941;Goldreich & Sridhar 1995;Boldyrev 2006).
The RIAF systems most relevant for the EHT are much better described as Coulomb collisionless, but when the ratio of gas to magnetic pressure is large (which is a likely description of much of the accretion flow), perturbations in the magnetic field can drive sufficient deviations from local thermodynamic equilibrium to trigger kinetic micro-instabilities, which are non-local in nature and increase the effective collisionality beyond what the naïve Coulomb collision picture implies (Kunz et al. 2014;Sironi & Narayan 2015;Melville et al. 2016;Riquelme et al. 2016;Squire et al. 2017;Bott et al. 2021Bott et al. , 2023;;Arzamasskiy et al. 2022;Ley et al. 2022;Tran et al. 2022).This enhanced collisionality can lead to considerable dissipation close to the injection scale due to pressure-anisotropic viscous stress (Arzamasskiy et al. 2022;Squire et al. 2023).
To avoid committing to a particular mechanism for dissipation, we adopt the sigmoidal R low -R high model (Mościbrodzka et al. 2016;Event Horizon Telescope Collaboration et al. 2019b), in which the ion-to-electron heating ratio smoothly transitions between two asymptotic values in regions with low and high plasma β = P gas /P mag .Although this form is quite simplified, it is straight-forward to implement and the sigmoidal shape is qualitatively supported by some studies of energy dissipation in collisionless plasmas (e.g., Quataert 1998;Rowan et al. 2019;Arzamasskiy et al. 2022).

The helicity barrier
When plasma β ≪ 1, further complexities arise if the energies of waves propagating in opposite directions are unequal, i.e., when the turbulence is imbalanced.This condition is typical in the case of the solar-wind plasma, but it may also occur in black hole accretion when strong outflows produce an imbalance biased along the outflow direction.How does imbalance alter the picture of a turbulent cascade?In imbalanced turbulence, the fluid is endowed with non-zero helicity, which must be conserved across the cascade in addition to the standard conservation of energy flux.In the inertial range (on scales k ⊥ ρ i ≪ 1 with ρ i the ion Larmor scale), both energy and helicity, which takes the form of a cross-helicity, can cascade simultaneously towards smaller scales; however, in the kinetic range, the cross-helicity is conservatively transformed into a magnetic helicity, the dispersion of the waves changes (Alfvén waves are converted into kinetic Alfvén waves, which are dispersive), and there is no solution that conserves the fluxes of both the energy and the generalized helicity.
The lack of solution results in an effective helicity barrier (Meyrand et al. 2021;Squire et al. 2022), as the unbalanced portion of the cascading energy is trapped at scales k ⊥ ρ i ∼ 1.The energy accumulates at that scale until other cascade directions are enabled (Squire et al. 2022 found that helicity barrier allows energy to enter a cascade of ion-cyclotron-waves, which eventually dissipate though ion-cyclotron heating).The imbalanced portion of the cascade thus only energizes the ions, and the maximum energy the electrons can receive is the balanced portion of the energy flux, which is itself divided between ions and electrons.The level of imbalance is quantified by the normalized quantity σ c ∈ [−1, 1], whose absolute value increases to unity as the level of imbalance grows.The primary effect of the helicity barrier that we consider in this paper is thus the reduction of electron heating, To compute σ c , it is useful to work in the Elsässer formulation of magnetohydrodynamics.In the relativistic context and written in terms of the fluid four-velocity u µ and magnetic field four-vector b µ = −u ν (⋆F ) µν , with ⋆F µν the Hodge dual of the electromagnetic Faraday tensor, the Elsässer variables are (Chandran et al. 2018) where the enthalpy is and where ρ is the rest-mass density of the fluid, u is its internal energy, and P is its pressure.The standard interpretation of z µ ± is that they describe the evolution of (pseudo-)Alfvén waves propagating through the equilibrium magnetic field.
Describing the fluid as a mean background with fluctuations, the fluctuations are just the differences between z µ ± and their locally time-averaged values, and it is easy to show that the reduced relativistic Elsässer equations, written in terms of these difference variables, reduce to the standard equations of Newtonian reduced magnetohydrodynamics.
The Elsässer variables can be used to compute two ideal pseudoenergy invariants, δz µ ± 2 , where we have The sum of the two pseudoenergies is the total energy in the system, and the difference measures the preference to generate waves in one direction or another (thus when the difference is non-zero, the system generates imbalanced turbulence).The normalized difference is the normalized cross-helicity: (5) Notice that there is ambiguity in how to perform the average in Equations 3 & 4: When the system is variable, the fluid frame changes and the part of the electromagnetic field that is seen as the magnetic field by the fluid, b µ , changes with time.The quantities z µ ± should represent the mean background flow; in regions where a characteristic background can be identified, performing a direct average of the four-vector components is then acceptable.In contrast, in regions that are highly variable, a mean background may not be readily identifiable and the meanings of δz µ ± become less clear.In such regions, our averaging procedure yields smaller values for σ c , which one might heuristically expect since the rapidly varying magnetic field and flow geometry do not allow helicity to accumulate along a particular direction over a sustained period of time.We compute z µ ± as an average over the full duration of the simulations beginning after the transient from the initial conditions dies out.We have verified that decreasing the averaging window by a factor of two or four does not qualitatively change our results.
The value of σ c in each fluid snapshot cannot be used directly to compute the effect of the helicity barrier, since the latter arises because of accumulated crosshelicity.The physical picture is as follows: the crosshelicity injected at large-scales cascades down to smaller scales on the eddy turnover timescale, and eventually cross-helicity (of a particular sign) accumulates at the ion Larmor scale.The helicity barrier is effective only when cross-helicity has time to build up at the Larmor scale-injections of negative and positive cross-helicity at the large scales do cascade, but they ultimately cancel out.
Because we perform our analysis in post-processing, we cannot track the cascade and injection of crosshelicity over time since that would require tracking the flow of non-zero cross-helicity fluid parcels as they evolve with the fluid.In this analysis, we instead approximate the buildup of cross-helicity by averaging the signed value of σ c over approximately a dynamical time.This signed average is a good proxy for the total amount of accumulated cross-helicity in an axisymmetric flow with small radial velocities; we adopt this procedure even though our simulations are three-dimensional for the sake of computational efficiency.
Figure 1 shows the signed value of σ c in a snapshot compared against the average of σ c over one dynamical time at r = 3 GM/c 3 where much of the observed emission is produced.The location of the disk can be identified by the plasma density, and the maximal extent of the emission region is bounded by magnetization σ = 1 contours.Evidently, the effect of averaging is to smooth out small-scale fluctuations in σ c while the broader, large-scale features are left mostly unchanged.The sign of σ c in the disk region is determined by the instantaneous flow properties.In MADs especially, crosshelicity of a particular sign may be long lived, as transient vertically asymmetric features are launched from large radii and fall through the event horizon.We discuss this smoothing procedure and compare between different averaging windows in the discussion section (see especially Figure 10).

NUMERICAL METHODS
We use the PATOKA pipeline to produce simulated images of RIAFs assuming the Kerr geometry (Wong et al. 2022).The pipeline comprises a fluid simulation step, in which the general relativistic magnetohydrodynamics (GRMHD) code iharm3d (Gammie et al. 2003;Prather et al. 2021) produces the time evolution of the accretion flow in full 3D, and a ray-tracing step, in which the general relativistic radiative transfer code ipole (Mościbrodzka & Gammie 2018) is used to compute the emission, extinction, and rotation of polarized light throughout the fluid simulation domain and track it to an observer at large distance.

The fluid model
The fluid evolution is obtained by solving the GRMHD equations, which take the form of a hyperbolic system of conservation laws along with the constraint Here, the plasma rest mass density is ρ and its fourvelocity is u µ .The magnetic field is represented by the b µ four-vector.The spacetime geometry enters through the metric g µν , its determinant g, and the Christoffel symbol Γ α βγ .The symmetric rank-2 tensor T µν represents the stress-energy of the fluid, which has contributions from both the fluid and the electromagnetic field where here u is the internal energy of the fluid and the fluid pressure P is related to its internal energy by a constant adiabatic index γ with P = (γ − 1) u.

Radiative transfer
The time series fluid data are processed into simulated images with a radiative transfer post-processing step using the ipole code (Mościbrodzka & Gammie 2018).Each simulated image comprises a square grid of square pixels defined by a field-of-view (or width) in units of GM/c 2 , distance from observer to source d src , and orientations with respect to the black hole spin axis and midplane (inclination and position angle).Pixels report the Stokes parameters I ν , Q ν , U ν , V ν at their centers.
To construct an image, ipole first traces photon trajectories backward from the camera into the simulation domain by solving the geodesic equations where Γ is a Christoffel symbol, λ is an affine parameter, and k α is the photon wavevector.ipole then integrates forward along each geodesic trajectory to solve the polarized radiative transfer equation, which in flat space can be written where we have neglected scattering as its effect is negligible at the radio frequencies we are interested in.Here, emissivities j ν , absorptivities α ν , and rotativities ρ ν are frame-dependent quantities (Chandrasekhar 1960).To compute the transfer coefficients, we use the thermal fits described in Marszewski et al. (2021).Further detail about ipole can be found in (Mościbrodzka & Gammie 2018;Wong et al. 2022).
Because GRMHD simulations introduce numerical floors in regions with high magnetization σ = b 2 /ρ, the plasma density and temperature are unreliable in such regions.To avoid contaminating the simulated images with numerical artifacts from the floors, we set the plasma density to zero in regions with σ > 1. Applying this σ-cutoff is reasonable, as the true density in highly magnetized regions like the jet is very small and therefore very little emission is produced there.

Computing the electron temperature
Since the fluid simulations only track the total internal energy of the total fluid, there is freedom in assigning the electron distribution function.For M87*, radio frequency emission is produced by the synchrotron process (e.g., Yuan & Narayan 2014), and for the relevant plasma parameters, the 230 GHz emission observed by the EHT likely comes predominantly from the thermal core of the distribution function.We thus assume that the electron population can be modeled as a relativistic thermal Maxwell-Jüttner distribution, which is characterized by a single temperature T e .
The problem is thus to determine T e given the total internal energy u of the fluid and the local fluid properties, which requires partitioning the total internal energy u into an ion component u i and an electron one u e .Schematically, the internal energy can be written as Here, we have used the subscript h (or z) to denote the part of the internal energy that can be related to heating via the helical (or zero-helicity) part of the turbulent fluctuations.When β is small, the helicity barrier stops any of u h from cascading below the ion Larmor scale and heating the electrons, so u e,h = 0 subject to the condition that β < β critical .We set β critical = 1 to be consistent with the physical derivation of the barrier, but we have found that varying this cutoff value above unity has negligible impact on our results.
When β < β critical , we assume that u z = (1 − σ c ) u, i.e., that energy imbalance is equal to injection imbalance.This equivalence is likely not true in general: Schekochihin (2022) finds that energy imbalance is larger than injection imbalance, although how well the quantitative details hold in non-idealized scenarios is uncertain.Nevertheless, under our assumption, the ion and electron energies are simply For apples-to-apples comparison, we fix R ≡ u i,z /u e,z regardless of σ c , which is reasonable under the approximation that the balanced component of the turbulent cascade is unaware of the imbalanced component.The ratio of total internal energies is then Finally, to compute the electron temperature, we must find the relationship between the ion-electron temperature ratio R T ≡ T i /T e and the energy ratio R u , which we do by assuming an ideal gas equation of state.Let the internal energies be Taking 1/y and 1/z to be the number of electrons and nucleons (protons + neutrons) per (unionized) atom, respectively, then n e = yρ/m p , n i = zρ/m p , and we have that n i = zn e /y.The ratio of energies is therefore Assuming fully ionized hydrogen,2 which has y = z = 1, if the ions are nonrelativistic γ i = 5/3 and the electrons are relativistic γ e = 4/3, then We have not yet specified R, the ratio of ion-toelectron energies in the zero-helicity fluid component, which is in reality determined by the microphysics.To parameterize over this uncertainty, we let R take any form allowed by the R low -R high prescription (see especially Mościbrodzka et al. 2016) which is motivated by models for electron heating in a turbulent collisionless plasma that preferentially heats the ions when the gas pressure exceeds the magnetic pressure.Here β R ≡ β/β 0 , and β 0 , r low , and r high are parameters that control the temperature ratio in regions of low (high) β where the plasma is dominated by gas (magnetic) pressure.The value of β 0 determines where the transition between r low and r high occurs.We adopt typical values for r low and β 0 and set them each to unity.Since plasma β is large in disk regions, models with large r high mostly have cooler disks and, by contrast, hotter coronae and funnel walls, and thus often produce more emission from regions off the midplane.

RESULTS
We now use a set of GRMHD simulations to study the importance of the helicity barrier in simulated polarized observations of RIAFs.For simplicity, we focus on the M87* accretion system and so set the black hole mass to M = 6.5 × 10 9 M ⊙ for consistency with observational results (see Table 1 of Event Horizon Telescope Collaboration et al. 2019c).This mass choice provides a physical length scale to the simulations.We use Equation 5and the averaging procedure described above to compute σ c across the domain and calculate electron temperatures.We image each snapshot of the fluid simulation twice, once using electron temperatures computed incorporating the helicity barrier and once with the effects of the helicity barrier turned off.To compare against observations, the camera must be assigned both an inclination and a position angle, which we define relative to the jet launched by the system.For our target M87*, there is clear evidence of a large-scale jet (see Walker et al. 2018) with a measured inclination angle of 17 • relative to the line of sight.We therefore orient our camera at either 17 • or 163 • relative to the axis of the jet in the simulation (this is coincident with the black hole spin axis) according to which parity reproduces the observed brightness asymmetry seen in Event Horizon Telescope Collaboration et al. (2019b), which manifests as a greater brightness temperature in the bottom half of the image.The statistical axisymmetry of the accretion flow means that rotating the image to align the position angle of the jet with its observed value does not influence the other images statistics.We thus fix the position angle of the images so that the jet lies in the vertical direction, as determined by the default PA = 0 setting for the simulations.

Model parameter space
The space of possible accretion configurations is high dimensional, covering the black hole mass and angular momentum parameters, the accretion rate of the system, boundary conditions and gas composition, and the magnetic field configuration.It is not computationally feasible to explore the full parameter space, so we focus on the subset corresponding to the canonical models used in the initial Event Horizon Telescope analysis of M87*.We thus aim to identify general trends and gauge the overall importance of the helicity barrier rather than make quantitatively precise predictions.
The magnetization of an accretion system can be used to differentiate flows according to whether the magnetic pressure near the horizon is strong enough to counterbalance the inward ram pressure of the fluid.When the magnetic pressure is high enough, the infalling motion of the plasma is arrested and the flow enters the magnetically arrested disk (MAD; Bisnovatyi-Kogan & Ruzmaikin 1974;Igumenshchev et al. 2003;Narayan et al. 2003) state.The alternative scenario is canonically referred to as standard and normal evolution (SANE; Narayan et al. 2012;Sądowski et al. 2013).SANE flows are turbulent but steady; in MAD flows, large tubes of magnetic flux arrest the inward motion of the flow and accretion is chaotic and mediated by transient filaments of hot plasma that thread the region between the hole and plasma at large radius.We consider both the MAD and SANE accretion states.
We express the black hole angular momentum in terms of the dimensionless spin parameter a * ≡ Jc/GM 2 with |a * | ≤ 1, where J is the magnitude of the angular momentum.By convention, we set a * < 0 when the angular momentum of the accretion flow and the spin of the black hole are anti-aligned.There is no reason that the angular momenta of the hole and the flow must be precisely aligned or anti-aligned.Tilted systems have recently gained broad attention; for simplicity, however, we restrict our focus to systems with no tilt.We con-sider five black hole spins a * = −15/16, −1/2, 0, 1/2, and 15/16 (hereafter written as −0.94, −0.5, 0, 0.5, 0.94 to be consistent with EHTC publications).
Although computing the radiative transfer coefficients requires choosing mass-density and length scales, since the GRMHD equations and Equation 5 are invariant under these rescalings, it is possible to measure the degree of cross-helicity directly from the scale-free fluid snapshot variables before restricting to a particular observer inclination or black hole accretion system.In Figure 2 we show the simulation-averaged values both for |σ c | and plasma β.For |σ c |, we have computed the timeaverage of the absolute value of the signed quantity σ c that has been calculated per fluid snapshot as described in Section 2.2 (and shown in the right-most panel of Figure 1).We show the average of the absolute value to account for the fact that the infall timescale is often shorter than the timescale over which σ c changes sign, since the magnitude of σ c controls the helicity barrier.Thus, the non-zero imbalance in the midplane of the MADs is due to spontaneous symmetry breaking that does not average out before the fluid parcels carrying the cross-helicity fall through the event horizon.Different choices for averaging windows are considered in the discussion (see especially Figure 10).
In regions where β is large, the accretion flow takes the form of a turbulent disk, fluctuations are large, and σ c is smallest.This effect is most prominently seen in the SANE flows and flows with small a * .Since MAD flows have more consistent magnetic fields, σ c keeps the same sign over longer timescales; this is reflected in the characteristically larger values of σ c in the MAD flows.
In all cases, the helicity barrier operates most strongly in regions with the most ordered magnetic field.In SANE flows the most ordered fields live within the jet and its enveloping wind.These regions have low plasma β and are approximately bounded by the magnetization σ = 1 contour.In MAD flows, the field maintains order throughout the domain and helicity builds up nearly equally everywhere.The "funnel" regions in the low-spin cases exhibit particularly ordered fields, which arise as accreted magnetic field lines build up near the horizon and are less perturbed by, e.g., the strong torques that a spinning black hole would impart on them due to frame dragging.

Images and emission source
We have thus far explored σ c in global accretion models from an observer-agnostic perspective.To understand how the helicity barrier influences observables, it is necessary to adopt an emission model, i.e., we must both choose thermodynamic flow parameters and set the The average for |σc| is computed from the value of |σc| calculated for each snapshot according to the signed averaging procedure described in §2.2.The accretion flow is concentrated in disks near the midplane, where β is largest.The magnetization σ = 1 surface is plotted as a white line in each panel; σ increases towards the poles and decreases in the midplane.In our models, the vast majority of emission is produced in regions with σ < 1, so from an observational perspective the value of σc between the two white lines matters the most.
observer inclination.We will focus on observations of the M87* accretion flow targeted by the EHT.
We use the ipole code to produce polarimetric images at the 230 GHz operational frequency of the EHT.In Figure 3, we show example images produced from the same single fluid snapshot shown in Figure 1 evaluated using thermodynamic models that either do or do not incorporate the influence of σ c on the ion-electron energy partition as described in Section 3.3.Columns show the full polarized properties of the light, including total intensity, degree of local linear polarization Q 2 + U 2 /I, electric vector position angle 1 2 arctan U/Q (measured east-of-north or counterclockwise-from-vertical on the sky), and the degree of circular polarization V /I, respectively.The bottom two rows of Figure 3 show the same images at the top two rows but blurred with a 20 µas Gaussian to simulate the effective resolution of the Event Horizon Telescope.This blurring is particularly important when considering observations of, e.g., resolved linear polarization, which may be high when the resolution element is smaller than the spatial correlation length of the EVPA but which is decreased dramatically when blurring over regions with a rapidly varying EVPA.
Although the images with and without the barrier are produced from the same fluid model, they correspond to different accretion rates selected such that the average flux density is 0.65 Jy to be consistent with observations.Thus, although the morphology of the fluid in the underlying accretion flows is the same for the dif- 2 arctan U/Q (measured east-of-north or counterclockwise-from-vertical on the sky), and the degree of circular polarization V /I, respectively.The bottom two rows show the images after being blurred with a 20 µas Gaussian beam to simulate the effective resolution of the Event Horizon Telescope.Since including the helicity barrier produces cooler electrons, the number density of the plasma must be increased to compensate, leading to more depolarization in the polarimetric signature (especially visible in the lower parts of the EVPA images).It also leads to less lensed emission coming through the disk in the photon ring, which produces less depolarization.While the detailed structure of the circular polarization can change significantly, the overall level of circular polarization does not tend to change.
ferent images, the number density and magnetic field strength differ.Figure 4 shows the factor by which the accretion rate must be increased for the flux to match observations.For internal consistency with the scalefree GRMHD equations, the increased accretion rate requires any local energy density quantity be increased by the same factor.In our case, the plasma number density, the fluid internal energy, and the square of the magnetic field strength must all be increased by the value shown in Figure 4. MADs and especially large r high have the largest required increase, as emission in those systems tends to be in regions with the largest imbalance and the greatest importance of the helicity barrier.Incorporating the helicity barrier yields higher estimates for jet power, since the jet power scales directly with the accretion rate. 3iven an emission model, it is possible to evaluate how the helicity barrier alters the source morphology.In The helicity barrier results in cooler electrons, which require higher local number densities to produce the same observed target flux.The higher number densities correspond to larger accretion rates.Systems with large r high tend to exhibit the greatest increase, as the effect of the helicity barrier is highest in the jet regions, where much of the emission is produced in large-r high models.MAD models typically have larger overall increases in ṁ, since σc is typically larger across their entire domains (see Figure 2).The spin of each model is given by the closest labeled value on the x-axis.Figure 6.Comparison of emission-weighted σc for library models with (solid) or without (dashed) the effect of the helicity barrier.The area under each curve has been normalized to unity.MADs have larger values of σc across their entire domains, and the influence of the helicity barrier is evident in the change of the curve shape for all models.In SANE models with low r high , emission primarily comes from the disk region, which already has low values of σc.Thus, the helicity barrier has little impact on the SANE r high = 1 models.As r high is increased in the SANEs, more emission comes from the funnel wall region and the effects of the helicity barrier's operation become more evident.
funnel region that lies at the interface between the jet core and the disk typically has larger values of σ than the disk.As expected, including the effects of the helicity barrier limits emission from regions with large σ c . Figure 2 shows that emission in MAD models does not change drastically while in SANE models the emission tends to shift away from the funnel wall and toward the lowermagnetization disk region.The right panels of Figure 5 show this trend as well: emission shifts from regions of high σ c to small σ c while the characteristic magnetization σ in the emission regions shifts from large values to small values in the SANE flow.
Figure 6 shows how the emission source changes across all library models.Emission in regions with large σ c decreases as expected.MADs typically produce emission throughout their infalling regions regardless of the thermodynamics prescription; as more of their domain has large values for σ c , the effect of the helicity barrier is very evident as emission in regions with large σ c drops significantly, altering the shape of each curve.SANE models with significant funnel-wall emission (i.e., models with large values of r high ) are often the most strongly affected.SANE models with low r high , i.e., models where the majority of the emission comes from the disk, are almost completely unaffected.

Polarization
We now evaluate how including the effect of the helicity barrier can affect the linear polarimetric β 2 ob-servable, which has been used by the EHT to gauge the strength of the horizon-scale magnetic field and differentiate between different accretion models (Palumbo et al. 2020;Event Horizon Telescope Collaboration et al. 2021b).The complex β 2 coefficient measures the power in (amplitude) and orientation of (argument) the azimuthally symmetric mode of the linear polarization vector across the image.The final value of the β 2 coefficient is determined by both the structure of the magnetic field in the emitting regions of the flow and the degree of depolarization due to differences in, e.g., Faraday rotation as the light propagates through the flow.Since the spin of the black hole influences the structure of the magnetic field, there is a trend of ∠β 2 with spin, with higher values of |a * | producing more toroidal fields and pushing ∠β 2 towards zero (radial linear polarization pattern).It is worthwhile to understand how plasma physics uncertainties might complicate this relationship.
Figure 7 shows how both the amplitude and argument of β 2 change for the different models in our library.Broadly, the amplitude of β 2 decreases with the inclusion of the helicity barrier while the argument of β 2 is mostly unaffected.The amplitude typically decreases the most in MAD models.Since σ c in MAD flows is mostly consistent across the domain, the regions that contribute to the image are mostly unchanged so that the general image structure persists.The differences are instead due to the reduced emission-per-particle due to lower temperatures, which must be compensated by an increased number density.This renormalization results Figure 8. Approximate increase in Faraday depth across all images from all models in our library when the helicity barrier is turned on.We compute the Faraday depth for each image as a polarization-weighted average over all image pixels of the Faraday depth along the entirety of each geodesic.Since MAD flows see large σc across their full domains, the plasma density in emitting regions increases significantly and is reflected in the increased Faraday depth.The effects of the helicity barrier in SANEs is more confined to the jet funnel region, and so the increase in Faraday depth is largest when emission comes from those regions, i.e., especially in the models with large r high .The box extends between the first and third quartiles and the whiskers extend the standard 1.5 × the interquartile range.Points in the lower-right triangle are smaller or less variable with the helicity barrier.For both statistics, MAD models are mostly unaffected since their emission regions remain similar.SANE models with larger r high tend to shift emission towards the disk when the helicity barrier is included, which results in more extended emission and larger apparent ring diameters.The variability for these models increases since it is controlled by the turbulent dynamics of the infalling plasma closer to the event horizon instead of in the funnel wall (notice that the points trend towards the value of the red r high = 1 SANE points).
in both an increased accretion rate but, more importantly, also increased depolarization since the differences in the increased column density of plasma along neighboring lines of sight lead to more extreme differing levels of rotation over the course of the light's propagation.
The differences that produce more scrambled images are quantitatively related to the Faraday depth along the geodesics as emission travels from its source to the observer.Figure 8 shows a proxy for the factor by which the Faraday depth increases when the effects of the he-licity barrier are included.Our proxy is computed by first evaluating the total Faraday depth along the full geodesic for each pixel in the image and then computing the polarization P = Q 2 + U 2 -weighted average of these values over all image pixels.The increase in Faraday depth is significant in MAD flows for all models; in SANE flows, the Faraday depth increases are more evident in models with large r high , where the emission is more likely to arise in the large-σ c jet funnel regions.
This scrambling effect can be seen in the EVPA panels of Figure 3, especially in the lower parts of the images where the bottom panel (with its lower density) has more coherent EVPA compared to the top panel).In the blurred linear polarimetric maps in the same figure, it is clear that the linear polarization in the bottomright (southwest) part of the image decreases because the overall intensity is canceled out by the near-random phases of the neighboring pixels' EVPAs.In SANE models the polarization pattern is often already highly scrambled because the magnetic fields in the emission region are highly disordered.Since the images start out scrambled, increasing the number density of the flow does not have as noticeable an effect.
The increased optical depth through the disk also means that the image of the lensed photon ring will appear less depolarized (see the blue ring that appears in the "without barrier" resolved linear polarization image of Figure 3).When the disk is optically thin, each pixel contains contributions from the direct image as well as the lensed secondary (and so on) images.The lensed images exhibit a conjugate polarization signature; the contributions from the lensed images cancel in part and the summed final polarization signal is decreased (for more detail see Himwich et al. 2020;Palumbo & Wong 2022).The increased column density due to the effect of the helicity barrier on the temperatures means that the secondary image is less prominent and thus less cancellation happens along the relevant trajectories.

Ring diameter & variability
Finally, we check whether disregarding the helicity barrier can bias several other parameters inferred by the EHT.Here, we focus on the ring diameter (Event Horizon Telescope Collaboration et al. 2019c;Psaltis et al. 2020;Event Horizon Telescope Collaboration et al. 2022b,c), which has been used to test consistency of the observational data with the theory of general relativity, and on the variability in the compact-flux light curve, which has demonstrated notable disagreement between models and the observational data Event Horizon Telescope Collaboration et al. (2022b).
We measure a ring diameter for each image in our library with the ring extractor (rex) method described in §9 of Event Horizon Telescope Collaboration et al. (2019d, see also Chael et al. 2022).rex makes its measurement from the algorithmically identified "center point" of each image-the point that is most equidistant from the peak intensity along each of 360 equally-spaced rays cast from itself.The left panel of Figure 9 shows a sampling of ring diameters taken from our M87*-like models, which are at low inclination where a ring diameter measurement is easiest to perform.As can be seen in the figure, disregarding the helicity barrier in MAD models tends to increase the measured ring diameter slightly although the overall measurements stay roughly consistent.In SANE models the helicity barrier alters the electron temperatures such that the measured ring diameter is roughly consistent or slightly larger.The largest increase in measured ring diameter occurs for the models with large r high , where emission in the jet funnel is suppressed and the disk contributes much more significantly to the image.
The right panel of Figure 9 shows how the measured modulation index changes when the effect of the helicity barrier is included in the electron thermodynamics calculation.The modulation index is where σ ∆T and µ ∆T are the standard deviation and mean of the time series, respectively, measured over some interval ∆T .We use ∆T = 553 GM/c 3 (≈ 6.5 months for M87*) to be consistent with the timescale used in the EHT analysis of the Galactic Center, which found inconsistencies between data and observation.MAD models are mostly unaffected since the geometric extent of the emission region does not change significantly with or without the helicity barrier.In contrast, in SANE models, especially with larger values of r high , the emission tends to shift towards the disk when the helicity barrier is included, increasing the relative imprint of the turbulent dynamics near the horizon (which is more variable than in the funnel wall).

DISCUSSION
Our approach is subject to several limitations.First, our GRMHD simulations do not dynamically evolve electron temperatures and instead only track the total energy of the ion-plus-electron fluid, leaving the electron distribution function to be prescribed in postprocessing.This procedure relies on the assumption that the ratio of heating rates can be directly mapped to the ratio of temperatures, which may be a reasonable assumption if the majority of internal energy is locally generated but need not be the case.Second, our base heating model does not depend on the structure of the accretion flow and does not correspond to any specific dissipation mechanism.Additionally, our simulations do not take into account the effects of pressure anisotropy on dynamics of the flow and thermodynamics of ions and electrons.Additional limitations are due to our implementation of the helicity barrier physics.Although it is necessary to compute the strength of local fluctuations in the Elsässer variables, there is no clear way to calculate the mean flow ⟨z ± ⟩ due to the global geometric structure and the relativistic nature of the problem.In this work, we use temporal averages rather than spatial ones, and our results depend on the details of the averaging procedure.To estimate the uncertainty due to averaging, in Figure 10 we compare the measured value of β 2 for different averaging windows for a representative set of models and find that the choice is relatively robust.Our model assumes that inhomogeneities do not allow the locally generated helicity to be transported away.
Finally, even though the barrier is expected to form in low-β regions of turbulence, we have applied the barrierinduced heating reductions across the entire domain.We do not expect this distinction to be qualitatively important, as high-β regions have relatively small amounts of electron heating in any case.Nevertheless, the heating reduction due to helicity barrier applies only in the regions of the flow where the main dissipation mechanism is turbulence.In our simulations, we assume that the entire domain is dissipated through turbulence, which is most likely incorrect.The relative importance of different dissipation channels is not yet well understood.
It is worthwhile to consider whether the effects of the helicity could be reproduced with modifications to the canonical electron temperature prescription of Equation 26.To first order, the helicity barrier produces cooler electrons across the domain and thereby increases the temperature ratio T i /T e everywhere; this change could be emulated by increasing the r low parameter by a factor of a few to order ten.Cool electron populations like the one that would result from this change have been invoked to explain disagreements between observations and model predictions for jet power and polarimetric properties (Event Horizon Telescope Collaboration et al. 2019b, 2021b).
Is it possible to do better?The strength of the helicity barrier depends on the normalized cross helicity σ c .Comparing the panels of Figures 2 and 5 shows that there is no clear relationship between σ c and other fluid parameters, like σ or β.Any modification to Equation 26 would at least need to be a function of some other locally calculable quantity that is not readily identifiable, so it is not clear how to modify the prescription without introducing an extra complexity comparable to directly evaluating σ c .
Thus, while such a global approach would produce the same qualitative effects as incorporating the helicity barrier, the complicated structures seen in Figure 1 suggest that any global approach would be inaccurate in detail.How the inaccuracies due to this approximation would compare to other modeling uncertainties is a different question, and we caution that the sensitive dependence of the observables on the details of electron distribution function makes it challenging to evaluate any kind of Jacobian.Performing a rigorous comparison is thus well beyond the scope of this paper.

SUMMARY
We have studied the effects of imbalanced turbulence and the resultant helicity barrier in the context of radiative inefficient black hole accretion.We have computed the degree of cross-helicity buildup in a suite of numerical accretion simulations covering both magnetically arrested disk (MAD) and standard and normal evolution (SANE) flows and over a range of black hole spins.We have also used results from local simulations of non-relativistic low-β turbulence (Meyrand et al. 2021;Squire et al. 2022) to explore how including (or not) the helicity barrier in the imaging procedure can affect predictions for 230 GHz horizon-scale black hole images relevant for Event Horizon Telescope analyses (Event Horizon Telescope Collaboration et al. 2019b, 2021b, 2022b).
The local level of sustained imbalance determines the importance of the helicity barrier, which in turn limits electron heating.We have found that the imbalance tends to be smaller in regions of the flow with high plasma β (commonly found in the disks of SANE flows and flows with low black hole spin).In contrast, in regions with ordered magnetic fields, such as in the jet and its surrounding wind in SANE flows as well as throughout much more of the domain in MAD flows, imbalance persists, helicity builds, and electron heating is more restricted.Accounting for the helicity barrier thus causes emission to shift away from the funnel wall towards the lower-magnetization disk region in SANE flows, while the emission morphology is largely unaffected in MADs.
When comparing to observations, the total emission produced by a candidate accretion flow must match its observed value, and the cooler electrons require larger plasma number densities and magnetic field strengths.Thus, neglecting the helicity barrier can lead to underestimated accretion rates and inferred jet powers by more than a factor of two.The higher plasma densities also lead to increased Faraday depths and depolarization, resulting in decreased amplitudes of the polarimetric β 2 observable.Finally, we find that the inferred ring diameter and light curve variability modulation index are mostly unchanged for MAD flows but may increase for SANE flows, especially with large values of r high .The increased jet powers and decreased coherent polarizations due to inclusion of the helicity barrier may help explain some qualitative differences between observed EHT data and contemporary modeling efforts (Event Horizon Telescope Collaboration et al. 2019b, 2021b).
The authors thank Michi Bauböck, Andrew Chael, Matt Kunz, Elias Most, Eliot Quataert, Jonathan Squire, Jim Stone, and Muni Zhou for useful discussions and suggestions.The authors also thank the anonymous referee for useful comments and suggestions.G.N.W. was supported by the Taplin Fellowship.Support for L.A. was provided by the Institute for Advanced Study.

Figure 1 .
Figure 1.Effect of averaging signed σc over a dynamical time in an example snapshot from a MAD accretion flow with a * = 0.5.Left panel: normalized plasma density plotted in linear scale with contour showing the magnetization σ = 1 surface, within which all emission is produced.Right panels: signed σc across the domain with overplotted σ contour.The center panel shows the instantaneously computed σc and the right-most panel shows σc averaged over one rotation period at r = 3 GM/c 3 , where emission tends to peak.Although averaging smooths out smaller features, the overall magnitude of the final |σc| remains non-negligible.The grey circle at the left center of each panel outlines the black hole event horizon.
For each time series of images, we rescale the mass density of the accreting plasma until the average of the 230 GHz flux density light curve matches the observed (instantaneous) value of F 230 GHz = 0.65 Jy (see Appendix B.1 of Event Horizon Telescope Collaboration et al. 2019d).See Appendix D of Wong et al. 2022 for caveats and more detail about the flux-fitting procedure.

Figure 2 .
Figure 2. Average value of |σc| and plasma β across domain for ten different GRMHD simulations of RIAFs across MAD and SANE states and five different black hole spins a * = −0.94,−0.5, 0, 0.5, and 0.94.The average for |σc| is computed from the value of |σc| calculated for each snapshot according to the signed averaging procedure described in §2.2.The accretion flow is concentrated in disks near the midplane, where β is largest.The magnetization σ = 1 surface is plotted as a white line in each panel; σ increases towards the poles and decreases in the midplane.In our models, the vast majority of emission is produced in regions with σ < 1, so from an observational perspective the value of σc between the two white lines matters the most.

Figure 3 .
Figure 3. Example image from a MAD a * = 0.5 r high = 40 model snapshot with and without including the effect of the helicity barrier on the electron heating model.Columns show the total intensity, linear polarization Q 2 + U 2 /I, electric vector position angle 12 arctan U/Q (measured east-of-north or counterclockwise-from-vertical on the sky), and the degree of circular polarization V /I, respectively.The bottom two rows show the images after being blurred with a 20 µas Gaussian beam to simulate the effective resolution of the Event Horizon Telescope.Since including the helicity barrier produces cooler electrons, the number density of the plasma must be increased to compensate, leading to more depolarization in the polarimetric signature (especially visible in the lower parts of the EVPA images).It also leads to less lensed emission coming through the disk in the photon ring, which produces less depolarization.While the detailed structure of the circular polarization can change significantly, the overall level of circular polarization does not tend to change.

Figure 4 .
Figure4.Relative in accretion rate ṁ for all library models after including effect of the helicity barrier.Colors correspond to different values of r high .The helicity barrier results in cooler electrons, which require higher local number densities to produce the same observed target flux.The higher number densities correspond to larger accretion rates.Systems with large r high tend to exhibit the greatest increase, as the effect of the helicity barrier is highest in the jet regions, where much of the emission is produced in large-r high models.MAD models typically have larger overall increases in ṁ, since σc is typically larger across their entire domains (see Figure2).The spin of each model is given by the closest labeled value on the x-axis.

Figure 5 .
Figure5.Location and plasma properties of emitting regions for two sample GRMHD simulation snapshots.Only emission that contributes to the final image is considered (seeWong et al. 2022  for more detail).The left panels show the effect of the helicity barrier on the geometry of the emission region.The right panels show histograms of where the emission originates by σc and magnetization σ.In MADs the value of σc is comparable across the domain so that emission everywhere is decreased but significantly change in morphology; in SANEs, σc is greater along the funnel walls and thus the relative contribution from near the funnel vs. the midplane decreases.The white dashed lines show the magnetization σ = 1 contour, outside of which emission has been disallowed.The accretion rate for each of the models has been adjusted so that the flux density at 230 GHz for each image is 0.65 Jy to be consistent with observations of M87*.right panels of the figure show the characteristic mag-netization σ and cross-helicity σ c of the emission.The

Figure 7 .
Figure7.The effect of including the helicity barrier in the electron temperature model on the polarimetric observable β2 across the library models.The amplitude of the β2 coefficient tends to decrease more in MAD models due to increased scrambling across the image from higher Faraday depths.In all models, the distribution of ∠β2 is less affected, although sometimes tends towards the more toroidal configuration.The trend of ∠β2 vs. a * remains unchanged.

Figure 9 .
Figure9.Exploration of how including the effect of the helicity barrier on electron temperatures changes rex-fit diameters (left) and variability (modulation index; right) across library models.The color of each point indicates the model r high .Points would lie on the black lines if the helicity barrier had no effect; the grey dashed lines in the rex panel show ±10% deviation.Points in the lower-right triangle are smaller or less variable with the helicity barrier.For both statistics, MAD models are mostly unaffected since their emission regions remain similar.SANE models with larger r high tend to shift emission towards the disk when the helicity barrier is included, which results in more extended emission and larger apparent ring diameters.The variability for these models increases since it is controlled by the turbulent dynamics of the infalling plasma closer to the event horizon instead of in the funnel wall (notice that the points trend towards the value of the red r high = 1 SANE points).

Figure 10 .
Figure 10.Comparison of measured β2 values from representative accretion models for two different choices for the averaging windows of the signed σc quantity.The blue ranges corresponds to an average over one dynamical time at the characteristic emission radius r = 3 GM/c 3 and is the choice adopted in the rest of the paper.The green ranges show the effect of removing this average, which tends to increase the absolute value of σc.The red ranges show the values measured before incorporating the helicity barrier.The cyan bar corresponds the values inferred from EHT observations of M87* (Event Horizon Telescope Collaboration et al. 2021a).