Stochastic Averaging of Radiative Transfer Coefficients for Relativistic Electrons

Synchrotron emissivities, absorptivities, and Faraday rotation and conversion coefficients are needed in modeling a variety of astrophysical sources, including Event Horizon Telescope (EHT) sources. We develop a method for estimating transfer coefficients that exploits their linear dependence on the electron distribution function, decomposing the distribution function into a sum of parts each of whose emissivity can be calculated easily. We refer to this procedure as stochastic averaging and apply it in two contexts. First, we use it to estimate the emissivity of an isotropic κ distribution function with a high-energy cutoff. The resulting coefficients can be evaluated efficiently enough to be used directly in ray-tracing calculations, and we provide an example calculation. Second, we use stochastic averaging to assess the effect of subgrid turbulence on the volume-averaged emissivity and along the way provide a prescription for a turbulent emissivity. We find that for parameters appropriate to EHT sources turbulence reduces the emissivity slightly. In the infrared, turbulence can dramatically increase the emissivity.

1. INTRODUCTION Semi-analytic and numerical models of black hole accretion flows, including general relativistic magnetohydrodynamics simulations, are now routinely used for physical interpretation of observational data (e.g., Event Horizon Telescope Collaboration et al. 2019Collaboration et al. , 2022)).In many cases the dominant radiative process at the frequencies of interest is synchrotron; scattering is typically negligible.Synchrotron photons are emitted, absorbed, and their polarization is modified by interaction with the plasma.In the radiative transfer equation these interactions are encoded in the transfer coefficients for emission, absorption, and Faraday rotation (rotation of the polarization) and conversion (conversion of circular to linear polarization and vice versa).The transfer coefficients depend on the electron distribution function f .Calculation of synchrotron transfer coefficients from f is notoriously difficult (Schwinger 1949).Methods for numerically evaluating coefficients for isotropic distribution functions now exist (Marszewski et al. 2021), or even some coefficients for anisotropic distribution function (Verscharen et al. 2018;Galishnikova et al. 2023).Although the direct numerical evaluation of coefficients is possible, direct evaluation is not practical in numerical radiative transfer calculations.
Instead simple, easy to evaluate analytic approximations are typically used for a small set of common model distribution functions.For example, useful expressions exist for the relativistic isotropic thermal (Maxwell-Jüttner) and non-thermal (power-law or κ) electron distribution functions (e.g., Shcherbakov 2008;Huang & Shcherbakov 2011;Leung et al. 2011;Dexter 2016;Pandya et al. 2016Pandya et al. , 2018;;Marszewski et al. 2021).Nevertheless, one might want to evaluate a model for a new or slightly modified distribution function -for example a power law distribution with an imposed exponential cutoff.In general this would require an expensive numerical re-evaluation and re-fitting of the coefficients.
Corresponding author: Monika Mościbrodzka m.moscibrodzka@astro.ru.nl Here we propose a method for evaluating synchrotron transfer coefficients semi-analytically.Our method is based on the observation that the coefficients depend linearly on the distribution function f , so that, for example, the emissivity j ν (f 1 + f 2 ) = j ν (f 1 ) + j ν (f 2 ).This observation has been made before in the context of accreting black holes by Mao et al. (2017).It implies that if we can write the distribution function as a weighted sum of distribution functions for which accurate but approximate analytic transfer coefficients are known, then we can also evaluate the transfer coefficients as a weighted sum.Here we apply this notion in two contexts.
First, we develop an accurate, efficient procedure for estimating transfer coefficients for the κ distribution function, which has been used in studies of black hole accretion and jets to model electron acceleration processes (e.g., Davelaar et al. 2018;Event Horizon Telescope Collaboration et al. 2019;Davelaar et al. 2023).The κ distribution smoothly connects a quasi-thermal core at low momentum to a power-law tail at high momentum.It is known to be a good fit to electron populations in the solar wind (Vasyliunas 1968); for a recent review see Pierrard et al. (2022).We also show how our procedure makes it easy to incorporate a cutoff in the distribution function.
The application to κ distributions is inspired by work on "stochastic averaging" by Schwadron et al. (2010), and we adopt this term to describe our method of taking weighted averages of emission coefficients.
Second, we use stochastic averaging to write down turbulent transfer coefficients, which may be useful in estimating the emergent radiation from numerical simulations.Large eddy simulations of turbulence are common in astrophysics.Each resolution element, or zone, represents a region with unresolved turbulent substructure.The effect of that substructure on the radiation field can be modeled by replacing the emission associated with each zone's temperature, density, velocity, and magnetic field strength by a weighted average of emissivities over a distribution of temperatures, densities, etc.We provide a general expression for turbulent emissivity that depends on the covariances of the emissivity parameters in a turbulent flow, and estimate the magnitude of the change in emissivities for models of Sgr A*.
This paper is organized as follows.Section 2 describes various distribution functions and develops the necessary notation.Section 3 applies stochastic averaging to the κ distribution function.Section 4 applies stochastic averaging to turbulent transfer coefficients.Section 5 summarizes and identifies directions for future research.

PRELIMINARIES
The polarized radiative transfer equation for the intensities in Stokes parameters I, Q, U, V along a coordinate s is . (1) where we have neglected scattering.Here j designates an emissivity, α an absorptivity, and ρ a rotativity, i.e. a Faraday rotation or conversion coefficient.Altogether there are 11 transfer coefficients.Following Leung et al. (2011) and Pandya et al. (2016) the polarized synchrotron emissivity coefficients for an isotropic distribution are where p is the momentum and f (p) ≡ (1/n e )dn e /d 3 p is the distribution function per particle rather than per unit volume; similarly f (γ) ≡ (1/n e )dn e /dγ, where γ is the electron Lorentz factor.Here S = I, Q, U, V , δ is the Dirac delta function and K S and y c functions are given in e.g.Pandya et al. (2016).Evidently the emissivities depend linearly on the distribution function.We will work explicit examples for j ν later on, but the stochastic averaging procedure can be applied to all transfer coefficients since all depend linearly on f , which must in general be derived from kinetic theory.One commonly used model for f is the thermal or Maxwell-Jüttner distribution where Θ e = k B T e /(m e c 2 ) is the dimensionless plasma temperature.Notice that if Θ e ≫ 1 then K 2 (1/Θ e ) ≈ 2Θ 2 e .Accurate approximate expressions for the thermal transfer coefficients are available in Marszewski et al. (2021).
A second commonly used model for f is the κ distribution.A relativistic κ distribution proposed by Xiao (2006) reads: where κ and w are parameters, and N is the normalization constant.For κw ≫ 1, N (κ, w) ≈ (κ − 2)(κ − 1)/2κ 2 w 3 .Evidently f κ smoothly connects a thermal core at low γ to a power-law tail at high γ, which is more convenient than stitching a power-law tail to a Maxwellian core (e.g., Özel et al. 2000;Zhao et al. 2023).In the limit κ → ∞, f κ → f th with Θ e ≡ k B T e /m e c 2 = w.Notice that N goes to zero as κ → 2 from above, i.e. the number of electrons in the unnormalized distribution diverges.Notice also that the mean energy per electron diverges for κ ≤ 3. Now consider a distribution function with a continuous parameter P, f (γ; P) with a known emissivity j ν (P).For example, f might be the thermal distribution with P = Θ e .Suppose that we want the emissivity associated with a new distribution that can be written as a stochastic average over f : where dPF (P) = 1, so that F is a probability distribution over P. Then evidently This defines the stochastic averaging procedure.The averaging procedure can be applied to the emissivities because they are linear in f (Equation [2]).The remaining transfer coefficients -absorptivities and rotativities -for a general distribution function can be obtained from a linear transformation of the susceptibility tensor, and the susceptibility tensor can be written as an integral over df (γ)/dγ (see e.g., Pandya et al. 2018).The absorptivities and rotativities are therefore also linear in f and can be stochastically averaged.

SYNCHROTRON EMISSION FOR κ DISTRIBUTIONS
The synchrotron transfer coefficients for the κ function cannot, to our knowledge, be calculated analytically.Numerical results can be fit with complicated analytic functions that are restricted, due to computational cost, to specific ranges or values of the κ parameter.Here we develop an efficient method for calculating κ transfer coefficients from a stochastic average of thermal distribution functions.
So far all we have done is reproduce f κ in approximate form.Now suppose that we introduce a cutoff in energy by setting F = 0 for λ < λ min .Then where now Γ denotes the upper incomplete Γ function and λ ′ min = λ min wκ.Then and this immediately yields an approximate expression for the distribution function (which can be expanded in the limit of small λ ′ min ) as well as a method for obtaining consistent, efficient-to-evaluate transfer coefficients.
The results can be compared to transfer coefficients evaluated numerically by the above authors assuming the κ electron distribution function.Figure 1 shows examples of thermal, original κ and our averaged distribution functions together with their synchrotron emissivities, absorptivities (Stokes I is shown as an example) and rotativities (Stokes Q is shown as an example).The figure shows stochastically averaged emissivities, absorptivities and rotativities in comparison to κ fit functions from Pandya et al. 2016 andMarszewski et al. 2021.There is a good agreement between these two.For emissivities/absorptivities the difference is at most 6%/7%.Given that we use here only fit functions and our averaged eDF is not exactly the κ function, this agreement is remarkable.There is also good agreement between the ρ Q coefficient obtained from stochastic averaging of thermal coefficient from Dexter 2016 and ρ Q for κ distribution function from Marszewski et al. 2021.
While we find a very good agreement between transfer coefficients for κ DF and for thermal DF averaged without the cutoff, calculating transfer coefficients with a cutoff (λ min =0.001) shows that applying a cutoff to the DF is not equivalent to multiplying the transfer coefficients by the same exponential (or other) frequency cutoff factor (e.g., Davelaar et al. 2019), which may dominate the errors if the cutoffs are known precisely.The method presented above offers a convenient solution to this limitation.

Emission from accelerated electrons near black hole event horizon
The main application of the presented scheme for computing transfer coefficients is modeling emission (images as well as broadband spectral energy distributions) from non-thermal electrons in relativistic accretion flows and jets from black holes.Given the general transfer coefficients, which are no longer limited to certain values or ranges of the κ parameter, we can introduce more realistic (smooth) model for electron acceleration.In this exercise, the underlying model for plasma density, velocity, and magnetic fields is provided by a 3D general relativistic magnetohydrodynamics (GRMHD) simulation of magnetically arrested disk accretion around non-rotating black hole evolved for 30,000 GM/c 3 using GRMHD code ebhlight (Ryan et al. 2015).The details of the model setup are comparable to simulations presented in Event Horizon Telescope Collaboration et al. (2022).
The simulations are ray-traced with polarized relativistic radiative transfer code ipole (Mościbrodzka & Gammie 2018), where the scheme for averaging emissivities, absorptivities and rotativities is implemented.When ray-tracing, the GRMHD simulations are scaled to Galactic Center supermassive black hole system, Sgr A*.It has been long thought that flaring behaviour of this accreting black hole is a manifestation of electron acceleration (e.g.Özel et al. 2000).Here we present example ipole images of Sgr A* at single frequency of 230 GHz which roughly corresponds to the peak of the synchrotron emission and the observing frequency of Event Horizon Telescope (EHT).The purpose of this calculation is to illustrate the flexibility and performance of a ray-tracing scheme equipped with stochastically averaged transfer coefficients.The radiative transfer calculation makes the following assumptions about the electron distribution functions, which are not followed by the GRMHD simulation.The basic model assumes thermal electrons where electron temperatures Θ e are given by the R − β formula from Mościbrodzka et al. 2016, where β ≡ P gas /P mag is provided by the underlying GRMHD simulation.The electron temperatures are derived from formula: where R high = 160 and R low = 1 and proton temperatures, T p , are followed in the GRMHD simulation.Each nonthermal DF is centered around w = Θ e where Θ e is given by the thermal model but we allow non-thermal electron distribution function only above some threshold temperature Θ e > 5. Apart from this floor for the electron acceleration we also set a ceiling for acceleration energies by setting λ min = 0.001, corresponding to a maximum Lorentz factor of electrons γ ≈ 10 4 .Finally, we assume that the κ parameter is variable and, similarly to Θ e and w parameters, is a smooth function of β: where κ low and κ high are the model free parameters.We set κ high = 30 and κ low = 3.5 for which there is a smooth transition between nearly-thermal electrons in the high-β regions (usually closer to the disk equatorial plane) and nonthermal electrons in the low-β regions (usually in the jet regions).Together these conditions guarantee that κw ≫ 1 everywhere where a non-thermal DF is allowed.Notice that equation 18 would be difficult to implement with existing fits to the κ distribution, which extend only to κ = 7.
Figure 2 shows images of the GRMHD simulations assuming: purely thermal electrons (top panels) and non-thermal electrons (bottom panels).Notice that the thermal and non-thermal models have been scaled with slightly different mass scaling unit M = 8 × 10 17 (non-thermal) and M = 10 × 10 17 (thermal) so that in both cases the flux density is ≈ 2 Jy, characteristic for Sgr A*.The differences between purely thermal and non-thermal images at the chosen frequency are only subtle at this observing frequency.1In particular non-thermal Stokes I images are slightly more extended in size compared to the thermal ones.This is consistent with findings of Mao et al. (2017).Interestingly, the polarimetric (Stokes Q,U,V) images are also slightly more extended.With more extended images the ratio of direct emission and lensed (into a photon ring) emission fluxes changes; this may have an impact on the variability of the light curves synthesized from the images.
The integration of the image with stochastically averaged transfer coefficients took only approximately 1.5 times longer compared to the purely thermal image.It is therefore feasible to carry out more extensive parameter surveys of the non-thermal models, including models with different assumptions than those presented here.

APPLICATION TO COARSE-GRAINING OF TURBULENT EMISSION
The finite resolution of GRMHD models inevitably cuts off turbulent structures close to and below the grid scale l .The unresolved structures contain a distribution of temperatures, densities, field strengths, frequency of emission, and field orientations (collectively, the parameter vector P) that may cause the radiative transfer coefficients averaged over a region the size of a cell to differ from the coefficients associated with the average state, since the transfer coefficients are nonlinear functions of P.Here we apply the notion of stochastic averaging to explore the implications of turbulent sub-structure for the total intensity emissivity alone; similar considerations apply to the other transfer coefficients.
Using stochastic averaging, we define the turbulent emissivity where F is the distribution of the n parameters.
If the variations in emissivity parameters within a zone are small, and we do not need to capture the instantaneous distribution within each zone but only the ensemble-averaged distribution, then we may not do too badly with the ansatz that the parameters P are distributed like a multivariate Gaussian, i.e.
where Σ ij = C ij σ i σ j is the covariance matrix, C ij is the dimensionless correlation matrix, and ∆P i ≡ P i − ⟨P i ⟩.
Expanding to second order in ∆P, where j ν and its derivatives are evaluated at ⟨P⟩, and repeated indices i, j are summed over.The first order term vanishes because ⟨∆P i ⟩ = 0 by definition.Notice however that variation of emissivity with n e enters the turbulent emissivity if density varies in a correlated way with other parameters, as might be the case if density and temperature vary but pressure is approximately constant.
To go further, we need (1) the standard deviation σ i for each emissivity parameter P i ; (2) the correlation matrix C ij ; and (3) the Hessian of the emissivity with respect to its parameters.
We start by estimating the standard deviation σ P for each parameter P i inside a zone of size l .This is related to the structure function of the turbulence at scale l , i.e. ⟨(P i − ⟨P i ⟩) 2 ⟩ 1/2 .For a passive scalar P i in turbulence obeying Kolmogorov statistics with outer scale L and variance at the outer scale σ P (L), For a thermal distribution the emissivity parameters are: electron density n e ; electron temperature Θ e ; magnetic field strength B; the angle between the line of sight and the magnetic field θ (the observer angle); and the frequency on emission ν, which is Doppler shifted to the mean frequency of emission in a zone.
Next consider the correlation matrix C ij ≡ Σ ij /(σ i σ j ).This is a property of turbulence and thus of the flow.It will vary with the global flow structure and also with lengthscale.
For definiteness we evaluate two examples using a set of GRMHD simulations of aligned accretion around a black hole with spin a = 0.75.The simulations have adiabatic index γ = 5/3 and resolution 384 × 192 × 192.We compute the covariance in a narrow annulus within 0.1π rad of the midplane, with 3.8 < rc 2 /(GM ) < 4.2.For simplicity we ignore the velocity (frequency) variations and the variations in field direction, and consider only the gas temperature Θ = P/(ρc 2 ) rather than using a model for the electron temperature, which -because electron temperature assignment models commonly depend on plasma β -can alter the correlations.Thus we consider correlations between ρ, Θ, and B, averaged over azimuth and time.
The diagonal elements of the correlation matrix are 1 and it is symmetric, so there are 3 nontrivial entries.Finally consider the Hessian of j ν,I with respect to its parameters.To make an estimate of the importance of turbulence in changing the emissivity we must choose an emissivity and therefore a distribution function.We use a thermal distribution.Defining

A ≡
n e e 3 B sin θ 27 √ 2m e c 2 ; X ≡ 9πm e cν Θ 2 e eB sin θ ; (23) and working in the limit X ≫ 1, The general expression for the Hessian is readily evaluated but not very interesting.Again in the interests of getting to an estimate we adopt values for the parameters appropriate to the EHT source Sgr A*: and consider only variations in the parameters n e , Θ e , and B, assuming that Θ e varies like Θ. Notice that for these parameters X ∼ 10 2 , justifying our use of Equation ( 24).Then at last we find for the MAD models 1 2 That is, turbulence reduces the volume-averaged emissivity by 6%.The largest contribution to the sum in (26) is from the (i, j) = (Θ e , Θ e ) term, and is negative because the emissivity is concave with respect to Θ e .The second largest term is from the (n e , B) terms, and arises only because these quantities are anticorrelated.To sum up: for the EHT-relevant parameters considered here the emissivity correction due to turbulence is small and negative.Small negative corrections due to turbulence will alter the Sgr A* model normalization: a higher accretion rate, and therefore stronger magnetic fields, are needed to achieve the same millimeter flux density.This may be a point of concern for EHT modelers in both Sgr A* and M87*, since it will shift the accretion rate and jet power up, and potentially change the polarization properties.
Finally, notice that there are circumstances in which the turbulent correction is large.Consider the correction in the near infrared where the emissivity is dropping exponentially with increasing frequency.Then the dominant term in Equation 26is the (Θ e , Θ e ) term, which is ∼ 14.The second largest is the (B, B) term, which is ∼ 4. Evidently if the Hessian is large then the emergent radiation can be very sensitive to turbulent corrections.

CONCLUDING REMARKS
We have developed a notion of stochastic averaging of radiative transfer coefficients.This procedure, which is just a weighted sum of transfer coefficients over a distribution of parameters, can provide a quick, practical method for evaluating transfer coefficients.
In particular we have shown that one can apply stochastic averaging to construct radiative transfer coefficients for the κ distribution function from those for a thermal distribution.The calculation presented here is inaccurate for mildly relativistic plasmas, though mildly relativistic regions may not contribute significant emission anyway.
Stochastic averaging of transfer coefficients is valid for a non-relativistic electron distribution function, as shown by Schwadron et al. (2010); the non-relativistic electrons distributed into κ function will not emit synchrotron photons but can still contribute to Faraday rotation.
We have also applied stochastic averaging to estimate corrections to the emitted radiation from subgrid scale turbulence.We find that for parameters appropriate to Sgr A* in the millimeter, the turbulent corrections are small and negative.In the near-infrared the turbulent corrections are large (an order of magnitude) and positive.The turbulent emissivity calculation is valid for relativistic, mildly-relativistic and sub-relativistic plasmas.
The applications presented in this work stochastically average over a thermal distribution, and are therefore trivially applicable only to isotropic distributions.It would be interesting to generalize to anisotropic distributions (e.g.Galishnikova et al. 2023).
As a next step, one might also incorporate the stochastic averaging scheme into Compton scattering kernels to predict high energy emission from scattering of synchrotron photons on non-thermal electrons or scattering on turbulent substructures which are unresolved in the GRMHD models.Scattering off turbulent Compton kernel can be thought of as a parametric variability model which can be then tested observationally e.g. by X-ray observations.
Finally, we emphasize that the covariance between emissivity parameters needed in calculating the turbulent emissivity is a comparatively unstudied characteristic of turbulence (as of this writing we are not aware of other calculations of the covariance for MHD turbulence).It would be interesting to understand how this covariance depends on lengthscale and on flow parameters.
The work presented in this paper was initiated at workshop Modeling Plasmas Around Black Holes hosted by Lorentz Center in Leiden in fall of 2023.MM acknowledges support from Dutch Research Council (NWO), grant no.OCENW.KLEIN.113and NWO Athena Award.Regarding calculations in section 3.3, authors gratefully acknowledge the HPC RIVR consortium (www.hpc-rivr.si)and EuroHPC JU (eurohpc-ju.europa.eu)for funding this research by providing computing resources of the HPC system Vega at the Institute of Information Science (www.izum.si).We are grateful to Ben Prather and Vedant Dhruv for providing some of the calculations used in Section 4. This work was supported by NSF AST 20-34306.This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.CFG was supported in part by the IBM Einstein Fellow Fund at the Institute for Advanced Study, and also in part by grant NSF PHY-2309135 and the Gordon and Betty Moore Foundation Grant No. 2919.02 to the Kavli Institute for Theoretical Physics (KITP).We thank the referee for helpful comments that improved the quality of this paper.

Figure 2 .
Figure 2. 230 GHz polarimetric images of the magnetically arrested accretion flow around non-spinning supermassive black holes produced by ray-tracing simulations.Purely thermal and thermal+non-thermal models are shown in the top and bottom panels, respectively.In the bottom panels, the emission regions are slightly more extended.