**Dominic W. Pesce 1,2,\*, Daniel C. M. Palumbo 1,2, Angelo Ricarte 1,2, Avery E. Broderick 3,4,5, Michael D. Johnson 1,2, Neil M. Nagar 6, Priyamvada Natarajan 2,7,8 and José L. Gómez 9**

**Abstract:** We present estimates for the number of supermassive black holes (SMBHs) for which the next-generation Event Horizon Telescope (ngEHT) can identify the black hole "shadow", along with estimates for how many black hole masses and spins the ngEHT can expect to constrain using measurements of horizon-resolved emission structure. Building on prior theoretical studies of SMBH accretion flows and analyses carried out by the Event Horizon Telescope (EHT) collaboration, we construct a simple geometric model for the polarized emission structure around a black hole, and we associate parameters of this model with the three physical quantities of interest. We generate a large number of realistic synthetic ngEHT datasets across different assumed source sizes and flux densities, and we estimate the precision with which our defined proxies for physical parameters could be measured from these datasets.

**Keywords:** SMBHs; VLBI; ngEHT

#### **1. Introduction**

The Event Horizon Telescope (EHT) observations of the supermassive black holes (SMBHs) in M87 [1–8] and Sgr A\* [9–14] are the first in a new era of horizon-scale studies of black holes. The primary observational signature on horizon scales is the black hole "shadow", a ring-like emission structure surrounding a darker central region [15,16]. Simulations of accretion flows around SMBHs generically produce images that exhibit such shadows [5,13], which typically have a size comparable to that of the theoretical curve bounding the locus of impact parameters for photon trajectories that escape the black hole (i.e., the "apparent shape" of the black hole, from Bardeen [17]). A driving motivation for the EHT to pursue observations of M87\* and Sgr A\* was because these sources were anticipated to have the largest shadow sizes of all black holes on the sky [1].

The next-generation Event Horizon Telescope (ngEHT) will build on the capabilities of the EHT by improving (*u*, *v*)-coverage through the addition of more stations to the array, increasing baseline sensitivities by using wider observing bandwidths, and accessing finer angular resolution by observing at higher frequencies [18]. A natural question to ask is whether these improved capabilities will yield access to a larger pool of shadow-resolved

**Citation:** Pesce, D.W.; Palumbo, D.C.M.; Ricarte, A.; Broderick, A.E.; Johnson, M.D.; Nagar, N.M.; Natarajan, P.; Gómez, J.L. Expectations for Horizon-Scale Supermassive Black Hole Population Studies with the ngEHT. *Galaxies* **2022**, *10*, 109. https://doi.org/ 10.3390/galaxies10060109

Academic Editor: Luigina Feretti

Received: 6 November 2022 Accepted: 25 November 2022 Published: 2 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

SMBHs.  Access to a population of shadow-resolved SMBHs would thus provide an opportunity to make uniquely direct and self-consistent measurements of such spacetime properties, with attendant implications for studies of SMBH formation, growth, and co-evolution with host galaxies.

The suitability of any particular SMBH for shadow-resolving ngEHT observations depends primarily on three properties [19]:

The first of the above properties is set primarily by the mass of and distance to the black hole, while the latter two are more complex and depend also on the mass accretion rate and other physical conditions in the accretion flow. However, the detectability of horizon-scale structure from a SMBH does not guarantee the measurability of any particular quantity of interest; additional conditions must be met to ensure that, e.g., a black hole mass can be measured, or that the ring-like structure associated with the black hole shadow can be distinguished from other possible emission morphologies.

In this paper, we provide estimates for the number of SMBHs for which the ngEHT could plausibly make mass, spin, and shadow measurements. In Section 2, we define observational proxies for each of these quantities of interest that can be accessed from the horizon-scale emission structure. Section 3 describes our synthetic data generation procedure and our approach to estimating parameter measurement precision from ngEHT data. Our conditions for the measurability of each proxy are defined in Section 4, where we also report the number of objects expected to satisfy these conditions for each quantity of interest. We summarize and conclude in Section 5. Throughout this paper, we use the results from Pesce et al. [19] as our baseline for how many SMBHs satisfy the above three detection criteria as a function of *θ* and *Sν*. 1

#### **2. Measurable Proxies for Quantities of Interest**

For a given SMBH, the two primary quantities of scientific interest are its mass and spin, neither of which is directly observable by the ngEHT. Instead, analyses of ngEHT observations will need to identify and measure features of the emission structure that serve as proxies for the desired quantities, or else they will need to carry out some form of physical modeling to infer the SMBH mass and/or spin from the ngEHT data. For the proof-of-concept analyses presented in this paper, we pursue the former strategy.

#### *2.1. Proxy for SMBH Shadows*

One of the most generic predictions from simulated images of SMBHs is that the observed emission structure on event horizon scales should exhibit a ring-like morphology associated with the black hole shadow (e.g., [5,13]). Though it is possible for other processes to give rise to ring-like emission structures—e.g., the Einstein ring from a bright, compact emitter passing behind the black hole—in such cases the ring-like structure is expected to be transient. For the purposes of this paper, we thus consider the observation of a ring-like emission morphology to be a proxy for verifying the object's black hole nature. If we can determine from ngEHT observations that the emission structure from a particular object is ring-like—i.e., if we can discern the shadow—then we can identify that object as a black hole.

#### *2.2. Proxy for SMBH Masses*

The mass of a SMBH sets the physical scale for its associated spacetime metric, and all spacetime-sensitive length scales in the system should thus exhibit a proportionality with the gravitational radius,

$$\theta\_{\S} = \frac{GM}{c^2 D'} \tag{1}$$

with *M* the black hole mass and *D* its distance from Earth. The most observationally accessible length scale is the overall size of the emission region, which for a ring-like emission structure corresponds to the ring diameter, *d*. The EHT has demonstrated that black hole mass measurements for both M87\* and Sgr A\* can be made by calibrating the scaling relationship between *d* and *θg* using a large number of simulated images of the emission structure [6,12]. In this paper, we thus take *d* to be a proxy for *M*2; i.e., we assume that if *d* can be measured for a particular SMBH, then *M* can also be determined.

#### *2.3. Proxy for SMBH Spins*

The spin, *a*, of a SMBH has historically proven to be more difficult to measure than the mass; e.g., the EHT observations of M87\* and Sgr A\* have not yet yielded strong constraints on the spin of either SMBH [5,13]. There are a number of possible avenues for measuring *a* from horizon-scale images of SMBH systems (e.g., [23]), but the most observationally accessible of these approaches target the imprint of the SMBH spin on the horizon-scale magnetic field topology, which in turn can be accessed through observations of the linear polarization structure around the ring (e.g., [7,8]). Palumbo et al. [24] have developed a useful decomposition of the polarization structure in terms of a basis that captures the azimuthal behavior of the electric vector position angle (EVPA, i.e., the orientation of the linear polarization around the ring). This decomposition takes the form

$$\beta\_{\rm{fl}} = \frac{1}{S\_0} \int \int P(r, \phi) e^{-im\phi} r \mathrm{d}r \mathrm{d}\phi,\tag{2}$$

where (*r*, *φ*) are polar coordinates in the image, *P*(*r*, *φ*) = *Q*(*r*, *φ*) + *iU*(*r*, *φ*) is the complexvalued linear polarization field (with *Q* and *U* the standard Stokes intensities), and *S*0 is a flux normalization factor. When studying images of M87\* from GRMHD simulations, Palumbo et al. [24] found that the "twisted" morphology of the linear polarization pattern, quantified by the (complex-valued) *β*2 coefficient, is correlated with the spin of the black hole. It is now believed that this relation arises from a magnetic field geometry that evolves with the black hole spin: black holes with larger spins exhibit more frame dragging, and produce more strongly toroidal magnetic fields than lower-spin black holes [25]. Qiu et al. [26] further explored the connection between polarized image morphology and SMBH spin, finding that the asymmetry (*A*) of the Stokes *I* emission, the polarimetric *β*1 mode, and the modulus of the polarimetric *β*2 mode also encode spin information (though *β*2 continues to stand out as the most discriminating measurable parameter). In this paper, we thus take a joint measurement of *β*1, *β*2, and *A* to be our proxy for *a*.

#### **3. Synthetic Data Generation and Fitting Procedure**

To determine the region of the (*θ*, *Sν*) parameter space—and thus the number of SMBHs—for which the quantities of interest described in the previous section could be measured by the ngEHT, we carry out a series of model-fitting exercises using synthetic data. We use a model for the SMBH emission structure that captures the salient features relevant for measuring the physical quantities of interest. Per Section 2, these salient features include the diameter and thickness of the emitting ring, as well as the structure of the linear polarization pattern. As our parameterization of the SMBH emission structure, we thus use a polarized "m-ring" model [12,27] convolved with a circular Gaussian blurring kernel. This model is restricted to describing ring-like morphologies, but it can flexibly distribute both the total intensity and the linearly polarized flux about the ring using a

relatively small number of parameters. The emission structures produced by this model qualitatively match those expected from both simple analytic treatments (e.g., [28]) as well as numerical GRMHD simulations (e.g., [8,12]).

In our polarized source model, the Stokes I image structure is given by

$$I(r,\phi) = \left[\frac{S\_0}{\pi d}\delta\left(r - \frac{d}{2}\right)\sum\_{k=-m}^{m} a\_k e^{ik\phi}\right] \ast \left[\frac{4\ln(2)}{\pi l \mathcal{W}^2} \exp\left(-\frac{4\ln(2)r^2}{\mathcal{W}^2}\right)\right],\tag{3}$$

where ∗ denotes the convolution operation, *d* is the ring diameter, *W* is the FWHM ring width, and *δ* denotes the Dirac delta function. We enforce *α*0 = 1 so that *S*0 is the total flux density, and we also enforce *α*−*k* = *α*∗ *k* so that the image intensity is real-valued. We define *A* = |*α*1| to be the asymmetry parameter mentioned in Section 2 as potentially relevant for spin constraints.3 The linear polarization structure is similarly given by

$$P(r,\phi) = \left[\frac{1}{\pi d}\delta\left(r - \frac{d}{2}\right)\sum\_{k=-m}^{m}\beta\_k e^{ik\phi}\right] \ast \left[\frac{4\ln(2)}{\pi l \mathcal{W}^2} \exp\left(-\frac{4\ln(2)r^2}{\mathcal{W}^2}\right)\right],\tag{4}$$

where we now allow both *β*−*k* and *βk* to be free parameters because *P* is complex-valued in general.

We generate a number of synthetic SMBH images by gridding the (*d*, *S*0) parameter space, spanning [0.1, 100] *μ*as in *d* and [10−3, 1] Jy in *S*0, with both dimensions uniformly gridded on a logarithmic scale. We set *m* = 1 for the Stokes I emission, with both the real and imaginary parts of *α*1 uniformly sampled within [−0.5, 0.5]. For the polarized emission we set *m* = 2, with the real and imaginary parts of *β*0 and *β*−2 uniformly sampled within [−0.1, 0.1], the real and imaginary parts of *β*1 and *β*−1 uniformly sampled within [−0.05, 0.05], and the real and imaginary parts of *β*2 uniformly sampled within [−0.3, 0.3]. For all synthetic images, we enforce *W* = *d*/3. Though these choices are not unique, they cover a range of parameter values similar to that seen in the GRMHD simulations developed by the EHT collaboration [5,13]. An example polarized m-ring image generated using these specifications is shown in Figure 1.

**Figure 1.** Example polarized source model used for generating the synthetic data described in Section 3. The grayscale image shows the Stokes I emission, while the colored ticks mark the EVPA of the linear polarization structure. The length of each tick is proportional to the intensity of the linear polarization (i.e., |*P*|), while the color of each tick reflects the fractional polarization (i.e., |*P*|/*I*).

To generate synthetic ngEHT observations corresponding to the synthetic images, we use the ngehtsim4 package, which expands on the synthetic data generating functionality of the ehtim library [30,31]. We assume the observations are carried out at an observing frequency of 230 GHz and with 8 GHz of bandwidth using the "full" ngEHT Phase 1 array configuration from [32], which consists of the 2022 EHT array plus the OVRO 10.4 m dish, the Haystack 37 m dish, and three 6.1 m dishes located in Baja California (Mexico), Las Campanas Observatory (Chile), and the Canary Islands (Spain). We use historical weather data to determine appropriate system equivalent flux densities at each site following a procedure similar to that in Raymond et al. [33]. To emulate fringe-finding signal-to-noise ratio (SNR) thresholds, we flag any visibilities from baselines that contain a station not participating in at least one other baseline that achieves an SNR of 5 in a 10-s integration time. We add complex station gain corruptions at the level of 10% in amplitude and uniformly sampled within [0, 2*π*] in phase for all stations on every 300-s time interval, to emulate scans, and we assume that the data have been calibrated to remove polarimetric leakage effects.

We generate synthetic datasets across a grid in right ascension and declination, with spacings between grid points of 1 h in right ascension and 10 degrees in declination. To gather information on the performance of the array in different weather conditions and for different black hole structure realizations, we generate 100 instantiations of synthetic data at each grid location. We assume weather conditions typical for the month of April.

For each synthetic dataset, we estimate the precision with which the parameters of a polarized m-ring model fit to the data could be recovered. We compute these estimates using a Fisher matrix approach implemented within the ngEHTforecast5 package. This approach does not explicitly carry out fits of the model to the data; instead, it assumes that a "good" fit to the data has already been achieved, and it then provides an estimate of the uncertainty in each of the fitted parameters via a second-order expansion of the logarithmic probability density around the best-fit location. We compute parameter precision estimates assuming that the fits have been carried out using complex visibilities as the input data products, with broad priors on the station gain amplitudes and phases at every scan.

#### **4. Results: The Expected Number of Measurable SMBH Masses, Spins, and Shadows**

The results of the modeling exercises described in the previous section are summarized in Figure 2, which shows the sky density of sources expected to have measurable masses (top panel), spins (middle panel), and shadows (bottom panel). At each sky location, the plotted density corresponds to an average over 100 instantiations of weather conditions and source structure. Our criteria for determining whether a particular mass, spin, or shadow is deemed "measurable" are as follows:

Given the above measurability thresholds, we can see from Figure 3 that there is a hierarchy of measurement difficulty with increasing *Sν* and *θ*. The "easiest" quantity to measure is *d* (and thus the black hole mass), which can be recovered for ∼50 sources after integrating over the whole sky. The next most well-constrained quantities are those pertaining to the black hole spin, which we find can be recovered for ∼30 sources. The most difficult quantity to measure is *W* (and thus the black hole shadow), which can be recovered for ∼7 sources. The measurements are cumulative within this hierarchy: for all sources for which spin is measurable, mass is also measurable; for all sources for which the shadow is measurable, both spin and mass are also measurable.

**Figure 2.** Estimated sky density of SMBHs with measurable masses (top), spins (middle), and shadows (bottom), as a function of right ascension and declination. These estimates have been determined according to the criteria outlined in Section 4, and they assume an underlying distribution of observable SMBHs from Pesce et al. [19]. The stochastic variations seen from pixel to pixel are primarily the result of sampling noise. The location of M87\* is marked with a red star.

All three quantities of interest exhibit two regimes of non-measurability in Figure 3. For sources with flux densities below *Sν* 10 mJy, the source is too weak to be detected on most baselines, and there are thus simply insufficient data to enable significant constraints on the parameters of interest. For sources that are stronger than ∼10 mJy but smaller than several *μ*as, there can be many detected data points, but the source is insufficiently resolved to enable significant constraints on morphological parameters. In both regimes, all three quantities of interest exhibit a measurability tradeoff between *θ* and *Sν*. In the second regime, this tradeoff is such that it is possible to make a measurement for sources with smaller *θ* so long as they have sufficiently larger *Sν* (because increasing signal-to-noise ratio permits subtler features to be recovered), while in the first regime the tradeoff is reversed (because compact sources yield more detections—particularly on long baselines—than extended sources).

**Figure 3.** Approximate number density of SMBHs that are expected to satisfy different thresholds of measurability, assuming an observing frequency of 230 GHz. The background colorscale and contours mark the number density (per unit solid angle) of SMBHs that have flux densities greater than *Sν* and shadow diameters larger than *θ*, as a function of *Sν* and *θ* and assuming that sources are distributed isotropically on the sky [19]. The solid contours start with the thick contour indicating a count of 1 and then increase by factors of 10 towards the lower left, while the dashed contours each decrease by a factor of ten towards the upper right. The overplotted colored dashed contours indicate where various parameters of interest could be measurable for different combinations of (*θ*, *Sν*), assuming observations appropriate for the "full" ngEHT Phase 1 array observing at a declination of 10 degrees (i.e., averaged over right ascension). The red dashed contour marks the lower boundary of the region in which black hole mass can be measured, the green dashed contour marks the lower boundary of the region in which black hole spin can be measured, and the blue dashed contour marks the lower boundary of the region in which black hole shadow can be measured.

Figure 2 shows the sky distribution of objects with measurable masses, spins, and shadows, after averaging over weather and source structure instantiations. We find that the distribution is quite uniform, and that accessible objects can be located almost anywhere

in the sky; there is no strong dependence on right ascension. The only major gaps in accessibility are for sources having declinations within ∼30 degrees of the southern celestial pole, for which the (*u*, *v*)-coverage of the array is particularly poor. A modest increase in source density is seen around declinations of ∼30–40 degrees, where the (*u*, *v*)-coverage of the array is densest. For the shadow measurements, we also see a modest increase in source density around the northern celestial pole; northern polar observations provide the most complete long-baseline coverage, so this bump in density may indicate that the long baselines are the most constraining for the width parameter.

#### **5. Summary and Conclusions**

To date, the EHT has observed the horizon-scale emission structure around two SMBHs. The ngEHT aims to improve on the capabilities of the EHT by adding new dishes to the array, increasing the observing bandwidth, and expanding the frequency coverage, all of which will improve the sensitivity and fidelity of reconstructed images.

Motivated by the promise of the ngEHT for population studies of SMBHs, we have identified three scientific quantities of interest that the ngEHT can expect to measure for a number of SMBHs: the black hole mass, the black hole spin, and the black hole shadow. We construct a geometric ring model for the polarized emission structure around a SMBH, and we identify parameters of this model as observable proxies for the scientific quantities of interest. Specifically, we associate the diameter of the ring with measurements of the black hole mass, the thickness of the ring with measurements of the black hole shadow, and the linear polarization structure with measurements of the black hole spin.

Assuming a Phase 1 ngEHT array configuration observing in April conditions at a frequency of 230 GHz, we generate a large number of realistic synthetic observations spanning a range of source structure (i.e., flux density *Sν* and angular size *θ*) and site weather (i.e., opacity and atmospheric temperature) instantiations. For each synthetic dataset, we use a Fisher matrix formalism to estimate the precision with which each of the geometric ring model parameters of interest could be measured. We use the statistics of these measurement precision estimates (across all weather instantiations) to determine the corresponding number of SMBHs on the sky whose properties could be well-constrained as a function of *Sν* and *θ*. We carry out this procedure for synthetic observations covering a grid in right ascension and declination, finding that the sky density of measurable sources (in each parameter of interest) is approximately uniform for declinations above roughly −60◦.

Associating these measurable parameters with their corresponding physical quantities of interest, we present estimates for the number of SMBHs for which the Phase 1 ngEHT can expect to make measurements of these quantities.  The measurable SMBHs have characteristic observed flux densities of ∼30 mJy and angular sizes of ∼10 *μ*as; per Pesce et al. [19], we expect the bulk of these SMBHs to lie in the redshift range between *z* ≈ 0.01 and *z* ≈ 0.1. Our estimate for the number of measurable shadows is consistent with the predictions from Pesce et al. [19].

We note that our detection criteria for mass and spin are likely optimistic. A primary analysis limitation is that our model for the appearance of an SMBH does not include emission that extends much beyond the near-horizon region. Mass estimates for SMBHs of interest to the ngEHT may be complicated by additional image features—such as, e.g., AGN jets (e.g., [34]—that could limit the ability to accurately estimate the ring diameter when it is only marginally resolved or weakly detected. For spin, the situation is even more uncertain: the EHT has already produced tight estimates for the ring parameters *β*1, *β*2, and *A* for M87\*, but it has not yet claimed a corresponding measurement of the black hole spin [6–8]. A secure association between these ring parameters and spin will require a combination of continued observational and theoretical studies.

On the other hand, we have employed a simplified analysis that likely underestimates the number of accessible sources, given any particular set of detection criteria. For instance, the synthetic datasets used in this paper are currently limited to April weather conditions and an observing frequency of 230 GHz; a more comprehensive exploration of year-round weather conditions and the addition of a 345 GHz observing band would likely increase the number of accessible sources. Furthermore, the synthetic datasets generated for the analyses in this paper have assumed an EHT-like calibration procedure; more advanced calibration strategies that can bootstrap phase information across frequency bands (e.g., [35]) are also expected to increase the number of accessible sources. Addressing these shortcomings will be the focus of future work.

Observationally, the most critical next step is to identify a list of credible targets and start surveying them to determine flux densities and compactness for ngEHT followup. Ramakrishnan et al. [36] are compiling a comprehensive sample of all plausible ngEHT AGN targets, which is expected serve as a source catalog for pursuing SMBH population studies with the ngEHT.

**Author Contributions:** Conceptualization, J.L.G., P.N., D.W.P.; methodology, M.D.J., D.C.M.P., D.W.P., A.R.; software, A.E.B., M.D.J., D.W.P.; writing—original draft preparation, D.W.P.; writing—review and editing, A.E.B., J.L.G., M.D.J., N.M.N., P.N., D.C.M.P., D.W.P., A.R.; visualization, D.W.P.; supervision, J.L.G., P.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** Support for this work was provided by the NSF through grants AST-1440254, AST-1935980, and AST-2034306, and by the Gordon and Betty Moore Foundation through grant GBMF-10423. This work has been supported in part by the Black Hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University. NN acknowledges funding from TITANs NCN19-058 and Fondecyt 1221421.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Notes**

#### **References**

### *Article* **Multi-Wavelength and Multi-Messenger Studies Using the Next-Generation Event Horizon Telescope**

**Rocco Lico 1,2,\*, Svetlana G. Jorstad 3,4,\*, Alan P. Marscher 3, Jose L. Gómez 1, Ioannis Liodakis 5, Rohan Dahale 1, Antxon Alberdi 1, Roman Gold 6, Efthalia Traianou 1, Teresa Toscano 1 and Marianna Foschi 1**

**Abstract:** The next-generation Event Horizon Telescope (ngEHT) will provide us with the best opportunity to investigate supermassive black holes (SMBHs) at the highest possible resolution and sensitivity. With respect to the existing Event Horizon Telescope (EHT) array, the ngEHT will provide increased sensitivity and uv-coverage (with the addition of new stations), wider frequency coverage (from 86 GHz to 345 GHz and higher), finer resolution (<15 micro-arcseconds), and better monitoring capabilities. The ngEHT will offer a unique opportunity to deeply investigate the physics around SMBHs, such as the disk-jet connection, the mechanisms responsible for high-energy photon and neutrino events, and the role of magnetic fields in shaping relativistic jets, as well as the nature of binary SMBH systems. In this white paper we describe some ngEHT science cases in the context of multi-wavelength studies and synergies.

**Keywords:** very long baseline interferometry (VLBI); supermassive black holes; active galactic nuclei; multi-wavelength studies; relativistic jets

#### **1. Introduction**

Cosmic systems powered by accretion onto supermassive black holes (SMBHs) are the most intense sources of energy in the universe. The tremendous power they generate manifests in a variety of forms: electromagnetic waves from radio to *γ*-ray wavelengths, ultra-energetic particles, strong magnetic fields, and plasma jets propelled outwards at near-light speeds. There is now firm evidence that SMBHs lie at the hearts of nearly all galaxies. The first direct measurements of the sizes and masses of SMBHs have been obtained with the Event Horizon Telescope (EHT) at 230 GHz for the radio galaxy, M87, (MB*H* = 6.5 ± 0.7 × 109 M*Sun*, [1]) and the Milky Way black hole, SgrA\*, (MB*H* ≈ 4 × 106 M*Sun*, [2]). These measurements show consistency with the predictions of general relativity spanning over three orders of magnitude in central mass. Systems where matter falls into the black hole at a high rate create the phenomenon of active galactic nuclei (AGNs). The masses of black holes in AGNs range from less than 1 million to roughly 10 billion solar masses (e.g., [3,4]). The luminosity of the jets can be as high as ∼1015 times that of the sun.

AGNs with powerful jets are rare in the universe, so most lie at distances of hundreds of millions or billions of light-years from the Earth. Much of their light comes from regions of a similar size as our solar system. Their extremely compact nature, combined with the great distances involved, causes these regions to appear extremely small in the sky; therefore, observations designed to study AGNs require telescopes with extraordinary

**Citation:** Lico, R.; Jorstad, S.G.; Marscher, A.P.; Gómez, J.L.; Liodakis, I.; Dahale, R.; Alberdi, A.; Gold, R.; Traianou, E.; Toscano, T.; et al. Multi-Wavelength and Multi-Messenger Studies Using the Next-Generation Event Horizon Telescope. *Galaxies* **2023**, *11*, 17. https://doi.org/10.3390/ galaxies11010017

Academic Editor: Massimiliano De Pasquale

Received: 15 November 2022 Revised: 10 January 2023 Accepted: 12 January 2023 Published: 14 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

resolving powers. Furthermore, the extremely luminous emission, and thus the physical conditions that create it, are highly variable with respect to time. Due to these two factors, instruments that probe the nature of AGNs and their jets must provide ultra-fine resolution imaging, as well as the ability to monitor rapid changes in the images.

Such opportunities and capabilities will be offered in the near future by the nextgeneration Event Horizon Telescope (ngEHT). The main goal of the ngEHT is to upgrade and enhance the current capabilities of the EHT array by adding more than 10 new stations, improving the image dynamic range by ∼2 orders of magnitude, enabling simultaneous multi-frequency observations up to 345 GHz and possibly higher frequencies, including polarization capabilities at all frequencies, and quadrupling the observing bandwidth [5]. The ngEHT will have an angular resolution of ∼10–15 micro-arcseconds (*μas*) and reach linear scales that are a fraction of a light-year. Moreover, the development of new superresolution algorithms (e.g., [6–8]) will allow the smallest features resolved by the ngEHT to be further improved by a factor of up to ∼3. Observations at 86 GHz will also be included, which are expected to be important for understanding the connection between the photon ring and inner jet. The ngEHT project development and implementation will consist of two phases. The first phase is based on the instrument design and site selection, and it is expected to be concluded by the end of 2023. The second phase will consist of the actual construction and commissioning of the ngEHT array, which is expected to be fully operational by 2030, with a gradual deployment in the period 2026–2030.

#### **2. Science Cases**

The EHT has already demonstrated the power of multi-wavelength (MWL) synergies for characterizing the underlying physics and emission mechanisms in AGNs. This was by means of an intensive multi-wavelength observation campaign of M87 conducted in 2017 with 19 different observing facilities (Figure 1, EHT MWL Science Working Group et al. [9]). In this section, we describe some potential science cases for the ngEHT in the context of MWL studies.

#### *2.1. Disk-Jet Connection*

The ngEHT will allow us to determine a jet's profile starting at its origin in the vicinity of an SMBH, the location and its effect on the jet and length of the jet flow's acceleration and collimation zone (ACZ), and the location of the Bondi radius (inside of which the SMBH's gravitational influence is important). This will provide crucial insights toward understanding the connection and interplay between the jet, the accretion disk, and the black hole. Currently, the jet-disk connection has only been reported for two radio galaxies: 3C 120 [10–12] and 3C 111 [13,14]. Both these radio galaxies possess relativistic jets with apparent speeds of ∼5c [15]. They exhibit X-ray properties similar to Seyfert galaxies, which implies that their X-ray emissions are produced in the accretion disk (plus its "corona" of hot electrons) in the vicinity of the SMBH. Monitoring with the Very Long Baseline Array (VLBA) at 43 GHz and X-ray satellites (*RXTE, Swift, and Suzaku*), a connection was established between X-ray flux variations and the emergence of new superluminal emission features ("knots") in the jet: appearances of new knots are preceded by dips in the X-ray luminosity. This X-ray/radio connection conclusively demonstrates a link between changes in accretion disk structure and powerful ejection events. Moreover, there is a similarity with the behavior of X-ray binary systems ("microquasars"), allowing comparisons that are important for generalizing insights obtained regarding the origin of relativistic jets in both types of systems [16].

**Figure 1.** Composite image of M87 as seen across the entire electromagnetic spectrum during the 2017 EHT observation campaign, adapted from [9] under the terms of the CC BY 3.0 licence. For each image, the instrument, the observing wavelength, and the scale are shown next to it. The frequency– wavelength scale at the bottom indicates the location of each instrument in the electromagnetic spectrum. Image credit: the EHT Multi-wavelength Science Working Group; the EHT Collaboration; ALMA (ESO/NAOJ/NRAO); the EVN; the EAVN Collaboration; VLBA (NRAO); the GMVA; the Hubble Space Telescope; the Neil Gehrels Swift Observatory; the Chandra X-ray Observatory; the Nuclear Spectroscopic Telescope Array; the Fermi-LAT Collaboration; the MAGIC collaboration; the VERITAS collaboration; NASA and ESA.

The innermost jet profile and its evolution with respect to distance from the SMBH are key factors in understanding a jet's collimation and acceleration. Recently, a number of studies have been devoted to this topic. Kovalev et al. [17] performed a search for transitions in the jet shape in a sample of radio-loud AGNs (367 sources) at 15 and 1.4 GHz. The authors found that, for nearby AGNs (redshift *z* < 0.07), 10 out of 29 jets exhibited a transition from a parabolic to a conical shape. They concluded that the transition in geometry is perhaps a common feature of AGNs jets, but they can be observed only when sufficient linear resolution is achieved. The authors also indicated that a jet transition occurs when the particle kinetic energy becomes equal to the energy of the magnetic field, while the Bondi radius position is governed by ambient medium pressure and may not coincide with the transition break. Casadio et al. [18] have analyzed stacked VLBI images at 86 and 43 GHz of BL Lacertae (BL Lac, *z* = 0.069). They found that its jet expanded with a conical geometry, with a higher expansion rate between ∼5 and 10 pc (de-projected) from the BH. They proposed that the transition in the profile is connected with a change in the external pressure at the location of the Bondi radius (∼3.3 × 105 R*g*).

The capabilities of the high-resolution ngEHT will allow direct imaging of these regions, tracing the jet down to the expected accretion disk and the areas where the jet is ejected and collimated in several AGNs. Such images and the studies they allow in combination with MWL observations, will provide decisive information for constructing and testing theoretical models.

#### *2.2. Nature of the VLBI Core*

VLBI images of blazar jets, the most extreme objects in the family of AGNs, are dominated by a compact feature known as the "core". The nature of the core remains an open question, mostly because of optical depth effects and the limited angular resolution of VLBI arrays. At wavelengths longer than ∼7 mm, the core is optically thick and essentially unresolved, probably representing the area in the jet where the optical depth is ∼1. At 0.8–3 mm, where the resolution is higher and the core is optically thin, it should represent a physical structure in the jet's acceleration/collimation zone, such as a standing

shock [19], a kink in the jet flow, a site of magnetic reconnection (Figure 2), or some other phenomenon. The EHT Collaboration imaged the core region of the blazar 3C 279, finding it to be dynamic with a puzzling multi-component structure [20]. The ngEHT will provide the data needed to make MWL movies to follow the motions, and perhaps the formation, of the components. This will lead to an understanding of the core, the brightest feature seen in relativistic jets, and how it is connected to the SMBH and its accretion disk.

**Figure 2.** Sketch of the innermost jet region in blazars (image credit: A.Marscher) [21].

A particularly important diagnostic of the nature of the core is its polarization. For example, a standing conical shock is expected to display a radial linear polarization pattern for viewing angles very close to the jet axis [22], a pattern that can be distinguished by the ngEHT. A kink can develop from an instability in a helical field, which has the signature of a rotation measure gradient across the jet (e.g., [23]) that can be detected at the resolution of the ngEHT. Magnetic reconnection can occur where a jet's magnetic field is highly turbulent, with low polarization and chaotic position angles of polarization across the jet. They can also originate where the magnetic field lines form loops that are stretched parallel to the jet, in which case the linear polarization vectors should lie transversely to the jet.

#### *2.3. Spine-Sheath Structure*

High-frequency synchrotron peaked (HSP) blazars, such as Mrk 421 and Mrk 501, are associated with a long-standing problem: how to reconcile extreme luminosities of teraelectronvolt (TeV) gamma rays, which require high relativistic Doppler beaming factors, with low apparent speeds observed in microwave VLBI images (known as the "Doppler crisis"; e.g., [15,24]). One of the most favored models to resolve the Doppler crisis involves the presence of a transverse structure in the jet, with an extremely fast (Lorentz factor > 20) inner spine surrounded by a slower outer sheath (e.g., [25]). This produces distinctive observable signatures in VLBI images, such as limb brightening, especially in linearly polarized intensities, with electric-vector position angles (EVPAs) that are transverse to the jet. All of these signatures have already been observed in the jets of both Mrk 421 and Mrk 501 (e.g., [26–28]), as well as in other blazars (e.g., [23,29,30]). However, as it was previously mentioned, very-high angular resolutions and sensitivities are needed to resolve the transverse structure in the innermost jet regions; these will become routinely available with the ngEHT. In fact, for the case of Centaurus A, EHT observations at 228 GHz have revealed a highly collimated, asymmeytrically edge-brightened jet together with a fainter counterjet [31]. In addition, a spine-sheath structure was recently detected in the bright quasar 3C 273 with the fine resolution of space-VLBI *RadioAstron* observations at wavelengths of 18 and 6 cm [32]. At 18 cm, the jet exhibited limb-brightening that had not previously been seen in this source. At 6 cm, an emission from the central stream of plasma was imaged, with a spatial distribution that was complementary to the limb-brightened emission, which indicated the presence of a spine. Regular monitoring with high resolution

and sensitivity is needed to determine the jet stratification structure and its evolution. Such observations do not yet exist, but will be possible with the ngEHT, which will allow measurements of apparent speeds across the jets of several HSP blazars.

#### *2.4. Supermassive Binary Black Holes*

Binary SMBH systems represent a natural product of hierarchical galaxy formation and merger events. Currently, one of the best candidates for hosting a binary SMBH system is OJ 287, a bright blazar showing quasi-periodic variability patterns in its optical lightcurve that could be explained as the result of the secondary SMBH impacting the accretion disk of a primary one [33,34]. However, alternative scenarios to the binary SMBH system can also explain the variability pattern (e.g., [35]).

Moreover, according to current theoretical models, when the system reaches an orbital separation of 0.1–10 parsec the dynamical processes that drive the coalescence of the two SMBHs cease to be effective and inhibit the merging process. This is known as "the final parsec problem". The ngEHT, with a resolution finer than ∼15 *μas*, will represent a unique opportunity to detect and spatially resolve from tens to hundreds of sub-parsec binary SMBH systems (see Pesce et al. [36] for more details about the statistics and the expected detections). Spatially resolving such binary SMBH systems would provide the ultimate proof of their existence. With its monitoring capabilities, an ngEHT will allow us to determine and characterize the orbital parameters of such tight binary systems, whose orbital periods are expected to range from months to a few years, and fully resolve their orbital motion. This is currently being investigated in a set of numerical simulations by using the Black Hole Accretion Code (BHAC, [37]), see Figure 3.

**Figure 3.** Equatorial slice of a binary SMBH system in the corotating frame interpolated on the Black Hole Accretion Code (BHAC) adaptive mesh grid. The colorscale shows the spacetime lapse function and the arrows show the shift vector.

Ongoing and upcoming large surveys in the optical to X-ray bands (e.g., the European-Chinese X-ray mission Einstein Probe) will provide prime lists of new binary SMBH candidates, among which the best-suited ones can be followed up with the ngEHT to check the source association and measure the orbital parameters. Moreover, the ngEHT highresolution observations would also represent an important resource for future gravitational wave detectors, both ground-based ones similar to the pulsar timing arrays (PTAs) and space-based ones similar to the Laser Interferometer Space Antenna (LISA), which is sensitive to the low-frequency band of the gravitational wave spectrum expected from binary supermassive black holes.

#### *2.5. High-Energy and Neutrino Events*

Blazars are very efficient cosmic accelerators and have been proposed as a possible source of several past high-energy neutrino events, although the production mechanisms and the precise region(s) where they form are still unknown [38]. The identification of the astrophysical counterparts of high-energy neutrinos is not straightforward (see e.g., [39]), and, in this context, the ngEHT could be helpful in two ways. On one hand, the highresolution mm-wavelength ngEHT observations will allow us to identify and constrain the region(s) responsible for neutrino production by tracing the variations and delay between the radio emission from the innermost jet regions and the neutrino detection event. On the other hand, owing to its monitoring and observing agility capabilities, an ngEHT will help to more accurately localize neutrino emission sources by triggering prompt observations when a detection event occurs (see Y. Y. Kovalev et al., this issue of *Galaxies*).

The ngEHT will also provide essential details about the innermost regions of blazar jets, allowing us to disentangle the different regions where visible X-ray and *γ*-ray emissions originate (e.g., near the SMBH, in the jet's acceleration zone, in a standing shock, in a moving knot, in a kink, in regions of magnetic reconnection, etc.), and to better understand the mechanisms and physical processes responsible for the observed emission. The time-domain imaging capabilities of the ngEHT will allow events to be matched in the jet with outbursts of high-energy photons with high accuracy. Moreover, multi-epoch and multi-frequency ngEHT polarization observations will allow investigations into the possible connection of magnetic field variations to the high-energy emission events (e.g., [40]). This can be achieved by determining the time evolution of fundamental polarization parameters, such as the fractional polarization, rotation of the polarization plane, and Faraday rotation measure.

#### *2.6. Polarization and Magnetic Field Studies*

The ngEHT will allow the characterization of a polarized radio emission in the near vicinity of the SMBH and constrain the magnetic field geometry, the plasma temperature, and the underlying particle density. This will determine whether the magnetic field is helical, as predicted by most theories, turbulent, or compressed by shocks (e.g., [41–43]). The presence of an ionized medium with a magnetic field causes the direction of linear polarization to rotate as radio waves pass through it; this is also known as Faraday rotation. MWL observations with the ngEHT can measure this effect to explore the medium both inside and surrounding the jet near its vertex. This can determine whether the jet is confined by a wind, an accretion flow, or its own magnetic field. Moreover, with the ngEHT it will be possible to measure the velocity field at the base of a jet, allowing us to test the different magneto-hydrodynamic jet formation models on scales of hundreds to thousands of Schwarzschild radii, and to conclusively determine whether the jets are powered by the extraction of rotational energy from the spinning black hole (Blandford– Znajek process, [44]) or from the accretion flow (Blandford–Payne process, [45]).

Synergies with current (Imaging X-ray Polarimetry Explorer—IXPE [46]) and future (e.g., e-XTP [47], COSI [48]) high-energy polarimeters could also prove invaluable. The combination of broadband radio polarization from the ngEHT and high-energy polarization observations will allow us to distinguish between shock and magnetic reconnection models for particle energization in the jets (e.g., [19,49]), as well as accurately pinpoint the particle acceleration regions. This could have significant implications for models of jet propagation as well as the nature of the blazar core (see Section 2.2). It will also bring a unique view of the spatial evolution of the magnetic field configuration on diverse scales from the vicinity of the black hole to pc and kpc away.

#### **3. Technical Requirements**

In order to measure Faraday rotation and gradients in the magnetic field of the most compact regions of jets, and other *μ*as-scale emission regions, simultaneous multi-frequency ngEHT polarization observations are required (see A. Ricarte et al., this issue of *Galaxies*). This can be implemented in several ways: (1) by using antennas whose frequency can be switched quickly, as is the case, for example, for a VLBA; (2) by combining ngEHT observations with other VLBI arrays (e.g., VLBA or next-generation VLA, at 43 and 86 GHz, respectively), or with other future space observatories (e.g., Millimetron); (3) by splitting the ngEHT into sub-arrays operating at different frequencies.

To achieve the science goals presented in this paper it will be essential to access the innermost jet regions, with typical 230 GHz flux densities ranging from a fraction to a few tens of mJy. For this reason, a sub-mJy sensitivity will be an essential requirement for the ngEHT, which is achievable thanks to next-generation receivers and the addition of new telescopes to the array. More details about the sensitivity requirements will be presented in M. Johnson et al. in preparation. In order to unambiguously connect MWL flaring and neutrino emission events to changes in the jets, monitoring capabilities are required for the ngEHT at a cadence shorter than the lifetime of the above-mentioned multimessenger events. This could be possible during 1–2 week campaigns in the fortunate case that we are indeed observing a source with such an event. Either year-long or multiple weeks-long sessions per year would produce a high chance of success. Such multi-epoch ngEHT observations will need to be coordinated with other existing facilities, as well as with the new generation of instruments, such as the Cherenkov Telescope Array (CTA), LISA, and Athena. Moreover, the Russian-led next generation millimeter-band Space VLBI mission Millimetron, when launched, will provide a space arm to the ngEHT. This will enhance the resolution of the ngEHT to provide both space VLBI imaging, although with somewhat limited uv-coverage, and uv-domain studies of photon rings.

**Author Contributions:** Conceptualization and writing—original draft preparation, R.L., S.G.J. and A.P.M.; writing—review and editing, all authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** R.L., J.L.G., and R.D. acknowledge financial support from the State Agency for Research of the Spanish MCIU through the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709), from the Spanish Ministerio de Economía y Competitividad and Ministerio de Ciencia e Innovación (grants AYA2016-80889-P, PID2019-108995GB-C21), the Consejería de Economía, Conocimiento, Empresas y Universidad of the Junta de Andalucía (grant P18-FR-1769), the Consejo Superior de Investigaciones Científicas (grant 2019AEP112). A.P.M. and S.G.J. acknowledge financial support from NASA Fermi Guest Investigator grants 80NSSC20K1567 and 80NSSC22K1571.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Communication* **Probing Neutrino Production in Blazars by Millimeter VLBI**

**Yuri Y. Kovalev 1,2,\*, Alexander V. Plavin 3, Alexander B. Pushkarev 3,4 and Sergey V. Troitsky 5,6**

	- 4 Crimean Astrophysical Observatory, Crimea, 298409 Nauchny, Russia
	- 5 Institute for Nuclear Research of the Russian Academy of Sciences, 60th October Anniversary Prospect 7a, 117312 Moscow, Russia<EMAIL_ADDRESS>	- 6 Faculty of Physics, M.V. Lomonosov Moscow State University, 1-2 Leninskie Gory, 119991 Moscow, Russia
	- **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The advancement of neutrino observatories has sparked a surge in multi-messenger astronomy. Multiple neutrino associations among blazars are reported while neutrino production sites are located within their central (sub)parsecs. Yet, many questions remain on the nature of those processes. The next generation Event Horizon Telescope (ngEHT) is uniquely positioned for these studies, as its high frequency and resolution can probe both the accretion disk region and the parsec-scale jet. This opens up new opportunities for connecting the two regions and unraveling the proton acceleration and neutrino production in blazars. We outline observational strategies for ngEHT and highlight what it can contribute to the multi-messenger study of blazars.

**Keywords:** neutrinos; active galaxies; galaxy jets; quasars; radio continuum; interferometric techniques

#### **1. Introduction: Current Status of High-Energy Neutrino Studies, Blazar–Neutrino Connections**

Neutrino observatories, such as IceCube, ANTARES (Astronomy with a Neutrino Telescope and Abyss environmental RESearch project), and Baikal-GVD (Gigaton Volume Detector) have been convincingly detecting astrophysical neutrinos at TeV to PeV energies [1–4]. Despite these observations, little was known about the origin of energetic astrophysical neutrinos until recently.

Blazars, a class of active galactic nuclei (AGN), have been considered as potential neutrino sources since the very early days of multi-messenger astronomy [5]. Observational evidence for a blazar–neutrino connection has started to emerge in recent years. First, the blazar TXS 0506+056 was associated with a high-energy neutrino, which coincided with a gamma-ray flare in 2017 [6]. This association was in contrast with a lack of systematic connection between gamma-ray-loud blazars and neutrinos (see, e.g., [7,8]). Then, numerous radio-bright blazars were shown to emit neutrinos with energies from TeVs to PeVs [9–16]. The detection of this correlation is driven by the unique capabilities of very-long-baseline interferometry (VLBI): the only technique able to directly probe and resolve central (sub)parsecs in AGNs at cosmological distances. Blazars emit neutrinos preferentially at the times of their flares (Figure 1), visible in radio bands [10,12,15,17]. Still, the neutrino production mechanism and the physical regions where it occurs remain unclear. The observed connection of neutrinos with radio emission from compact jet regions emphasizes the importance of high-resolution studies in answering these questions. VLBI is the best direct visual evidence we can obtain in astronomy.

For a general discussion of multi-wavelength and multi-messenger studies with the ngEHT, see Lico et al. [18]. In this paper, we present the progress in multi-messenger astronomy studies of cosmic neutrinos, their probable association with blazars, challenges

**Citation:** Kovalev, Y.Y.; Plavin, A.V.; Pushkarev, A.B.; Troitsky, S.V. Probing Neutrino Production in Blazars by Millimeter VLBI. *Galaxies* **2023**, *11*, 84. https://doi.org/ 10.3390/galaxies11040084

Academic Editor: Luigina Feretti

Received: 17 April 2023 Revised: 1 July 2023 Accepted: 5 July 2023 Published: 10 July 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and a critical role to be played by ngEHT [19–21] in addressing exciting open questions of high-energy neutrino production.

**Figure 1.** Radio and gamma-ray light curves of TXS 0506+056. Top: *Fermi* LAT weekly binning light curve [22] of the gamma-ray source 4FGLJ0509.4+0542 positionally associated with the quasar TXS 0506+056. The gray horizontal line denotes the median gamma-ray flux. Bottom: RATAN-600 light curve at 11 GHz. The radio light curve is decomposed of three radio flares depicted by the blue dashed lines, the sum of which is represented by a blue thick line [23]. The vertical lines denote the times of neutrino arrival [6,23].

#### **2. Neutrino Production in Blazars: Open Questions**

Assuming no particle physics beyond the standard model, astrophysical neutrinos with energies above TeV can only be produced in interactions of relativistic hadrons protons or nuclei—with ambient matter or radiation, see, e.g., [24] for a recent review. This fits the observational evidence discussed in Section 1 because the non-thermal radiation of blazars gives a clear signal that particles are accelerated there. However, both the amount of relativistic hadrons in AGN, and the degree to which these hadrons contribute to the observed electromagnetic radiation, are uncertain. Population studies suggest [24–26] that their contribution is small, and neutrino luminosities of blazars are orders of magnitude lower than photon luminosities. Consequently, one may imagine neutrino production in various places in a blazar and by means of different mechanisms.

The main challenge is to explain the production of neutrinos of very different energies, from a few TeV [11,27] to sub-PeV [6,10], in sources of the same class. For the *pγ* mechanism, expected to dominate in blazars [28], the wide neutrino energy range requires the presence of target photons with a very broad distribution of energies. Conventional models of high-energy neutrino production in AGN, known for decades, e.g., [29–31], as well as their modern versions, e.g., [32–34], often experience problems in explaining the lowerenergy part of the observed neutrino flux, particularly because the target photons from the accretion disk are expected to have energies ∼(10 ... 100) eV, while ∼10 keV are required for the intense production of ∼10 TeV neutrinos.

While neutrinos have already been associated with VLBI-bright blazars [11,15] and with their radio flares [10,12], these results were based on observations at centimeter wavelengths. There, synchrotron self-absorption prevents one from obtaining detailed spatio-temporal studies of the AGN central sub-parsec parts, e.g., [35]. To summarize, the open questions of the blazar–neutrino astrophysics are the following: (i) how are protons accelerated; (ii) what is the neutrino production process, *pγ* or *pp*; (iii) from where do seed (X-ray) photons originate from in case of *pγ*; (iv) where are neutrinos produced? Note that (ii) and (iv) can be different, and multi-zone models may be required to explain all observations consistently.

#### **3. Neutrino Astronomy in the ngEHT Era**

Currently, studies of high-energy astrophysical neutrinos and their sources are limited by the sensitivity and resolution of neutrino observatories. The situation is rapidly changing, as their capabilities are increasingly improving. The next-generation IceCube-Gen2 will grow the telescope volume ten-fold, from 1 to 10 km3, aiming at a corresponding increase in detection rates by 2033 [36]. The Baikal-GVD detector has already reached the effective volume of 0.5 km3 and continues to grow and improve event reconstruction algorithms [37]. KM3Net (Cubic Kilometer Neutrino Telescope), a neutrino observatory in the Mediterranean, is being constructed and has already started yielding its first results [38]. Together, these instruments will provide a qualitative leap in both the number of detected astrophysical neutrinos and their precise localization.

An increasing number of well-localized neutrinos will lead to reliable identification of individual blazars as neutrino sources. Moreover, it should be possible to highlight specific time periods with more prominent neutrino emission. This brings new challenges and opportunities to the EM counterpart of such multi-messenger studies.

The planned ngEHT array [19–21] will provide superior angular resolution, dynamic range and sensitivity in Stokes I and polarization at 3, 1.3, and 0.7 mm. This will allow scientists to observe and monitor the central (sub)parsecs of neutrino-emitting blazars at the highest resolution and frequency possible, significantly alleviating the synchrotron opacity problem of the current centimeter-wavelength VLBI. The ngEHT will be able to probe both the accretion disk region and the parsec-scale jet base, opening new opportunities for connecting the two regions and unraveling the proton acceleration and neutrino production in blazars.

#### **4. Planning ngEHT Experiments**

Below, we discuss several approaches to study and understand the physics behind the connection between neutrino production and EM activity from the jet upstream of the central engine—a possibility which will be realized by ngEHT. Before elaborating on observational campaigns, we note the following important complications of neutrino astrophysics that affect the suggested scenarios below. A typical probability of a neutrino with an energy above 100 TeV to be of an astrophysical origin is around 50%, and it drops significantly for lower energies [39,40]. A typical 90% error region of a highly probable high-energy neutrino is several square degrees [40,41]. Some neutrinos might arrive from nearby non-jetted AGNs [42] or even from our Galaxy and its relativistic objects [43–48]. On top of this, we know very little about the mechanisms of neutrino production in blazars; therefore, there is no streetlight under which we can plan our search.

We expect that a variety of blazars could be associated with neutrinos, allowing us to select optimal ngEHT targets by accounting for both their physical properties and technical or observational limitations. Within our current understanding and the experience accumulated from observational searches for high-energy neutrino counterparts, the following three scenarios for monitoring observations are suggested.

**Scenario 1:** Observation of blazars associated with selected new high-energy neutrino alerts immediately after neutrino arrival. Several blazar-associated high-energy alerts per year are expected. When two or three neutrino telescopes become fully operational, one might conservatively require two alerts for a given target to arrive within several days.

Pros: The most efficient strategy since it is linked to a specific event.

Cons: It will only be able to probe the state of an associated object after neutrino arrival.

**Scenario 2:** Observation of a sample of selected blazars reliably identified previously as neutrino sources. See Table 1 for the current most probable neutrino candidates.

**Scenario 3:** Observations of a complete VLBI-flux-density limited sample of 50–100 of the brightest blazars with a 3 mm VLBI flux density above 1 Jy [49,50].

Pros: Offers full temporal coverage of the expected events, with the possibility to compare neutrino-emitting and neutrino-non-emitting blazars to calculate the robust significance of a coincidence [15,51]. Furthermore, the strategy provides the option to combine such observations with other ngEHT cases [20].

Cons: Observationally expensive.

**Table 1.** Most probable neutrino candidates among the VLBI-selected bright blazars.

Notes: Publications that selected each blazar as a highly probable neutrino-emitter and measured their flux densities are referenced in the last column. The 0506+056 blazar was the first and only blazar distinguished by the IceCube, while the others were found by statistical analysis of complete VLBI-selected samples. The dates for high-energy neutrinos are shown in the format YYYY-MM-DD. † Estimated from nearby VLBI measurements at 15 and 43 GHz of MOJAVE and Boston University VLBA programs.

Tracing changes in the compact structure of blazars during and around periods of increased neutrino emission requires multi-epoch monitoring at the high resolution provided by ngEHT. To roughly estimate the required observation time, we expect that one imaging epoch per target will take 4–8 h. The observations should happen with a cadence between two weeks and one month (an estimate based on experience gained by the 7 mm blazar VLBA-monitoring program [53]) and produce polarization images with Stokes I dynamic range or better than 1000:1, preferably multi-frequency with a possibility for Faraday rotation measurements (RMs) and spectral analysis. From this, we will be able to constrain the following source properties.

flares in disks and in jets, e.g., [64,65] related to neutrino production if the resolution, sensitivity, and opacity permit. Observing in this regime, we will be able to overcome significant delays related to synchrotron self-absorption at lower radio frequencies (see Figure 1 and [60]).

We underline that studying a complete sample of AGN with understandable properties will allow us to not only relate the observed changes to detected neutrinos but also set a robust significance on that association, following the approach suggested by Plavin et al. [15].

#### **5. Synergy with Other Facilities**

The Square Kilometer Array, SKA [66] and especially the next-generation Very Large Array, ngVLA [67,68] going as high as 100 GHz will allow the monitoring of much larger samples of VLBI-selected AGNs as well as the faster imaging of neutrino arrival fields, and pre-selecting probable neutrino candidates for targeted ngEHT studies. Wide-field telescopes such as the optical Legacy Survey of Space and Time, LSST [69] will allow scientists to better associate blazars with neutrinos in cases where flaring activity is confirmed as a valid indicator, e.g., [6,15,70]. Moreover, optical and UV telescopes can separate flares occurring in jets and accretion disks, analyzing the optical color and polarization. Seed photons are expected from X-rays [11,64,71], where current and new-generation space X-ray telescopes will be very helpful. High energies, e.g., the Cherenkov Telescope Array, CTA [72,73], will continue to support gamma-ray/TeV–neutrino analyses and identify whether neutrino production zones are opaque to gamma-rays.

#### **6. Summary**

The ngEHT will revolutionize VLBI-imaging capabilities by bringing together the power of very high resolution, advanced dynamic range, and sensitive polarization data. What makes it unique, however, is its remarkable immunity to synchrotron absorption. It will allow for regions to be extensively probed from the accretion disk to the parsec-scale jet, e.g., [74], and study the most probable sources of high-energy neutrinos.

When ngEHT is fully operational, three large high-energy neutrino telescopes will be fully functional: IceCube, KM3NeT, and Baikal-GVD. This paper formulates the case, presents eight most probable associations to date, and suggests observational strategies to address very exciting and open questions concerning proton acceleration and neutrino production.

**Author Contributions:** Conceptualization and methodology, Y.Y.K.; writing—original draft preparation as well as review and editing, Y.Y.K., A.V.P., A.B.P. and S.V.T.; visualization, A.B.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This publication was funded in part by the Gordon and Betty Moore Foundation. It was also made possible through the support of a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the author(s) and do not necessarily reflect the views of these foundations.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank the ngEHT team for discussions, Eduardo Ros as well as anonymous referees for useful comments on the manuscript, and Elena Bazanova for language editing. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Emission Modeling in the EHT–ngEHT Age**

**Richard Anantua 1,2,3,\*, Joaquín Dúran 1,†, Nathan Ngata 4,†, Lani Oramas 1,\*,†, Jan Röder 5,†, Razieh Emami 3, Angelo Ricarte 2,3, Brandon Curd 1,2,3, Avery E. Broderick 6,7,8, Jeremy Wayland 9,10, George N. Wong 11,12, Sean Ressler 13, Nitya Nigam 2,3,14 and Emmanuel Durodola 15**

**Abstract:** This work proposes a methodology for testing phenomenologically motivated emission processes that account for the flux and polarization distribution and global structure of the 230 GHz sources imaged by the Event Horizon Telescope (EHT): Messier (M)87\* and Sagittarius (Sgr) A\*. We introduce into general relativistic magnetohydrodynamic (GRMHD) simulations some novel models to bridge the largely uncertain mechanisms by which high-energy particles in jet/accretion flow/black hole (JAB) system plasmas attain billion-degree temperatures and emit synchrotron radiation. The "Observing" JAB Systems methodology then partitions the simulation to apply different parametric models to regions governed by different plasma physics—an advance over methods in which one parametrization is used over simulation regions spanning thousands of gravitational radii from the central supermassive black hole. We present several classes of viewing-angle-dependent morphologies and highlight signatures of piecewise modeling and positron effects, including a MAD/SANE dichotomy in which polarized maps appear dominated by intrinsic polarization in the MAD case and by Faraday effects in the SANE case. The library of images thus produced spans a wide range of morphologies awaiting discovery by the groundbreaking EHT instrument and its yet more sensitive, higher-resolution next-generation counterpart, ngEHT.

**Keywords:** accretion disk; relativistic jet; GRMHD

### **1. Introduction**

With some of the highest-resolution images ever obtained in astronomy, the Event Horizon Telescope has probed the horizon scale of the supermassive black holes M87\* [1] and Sagittarius A\* [2]. Incidentally, both sources possess ringlike morphologies with diameters of 42 and 52 μas, respectively [1,2]. Based on a seminal 2017 data collection campaign, the ring imaged around M87\* was seen in 2019 to possess brightness asymmetry

**Citation:** Anantua, R.; Dúran, J.; Ngata, N.; Oramas, L.; Röder, J.; Emami, R.; Ricarte, A.; Curd, B.; Broderick, A.E.; Wayland, J.; et al. Emission Modeling in the EHT–ngEHT Age. *Galaxies* **2023**, *11*, 4. https://doi.org/10.3390/ galaxies11010004

Academic Editor: Dimitris M. Christodoulou

Received: 11 November 2022 Revised: 9 December 2022 Accepted: 15 December 2022 Published: 23 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

dominated by a Doppler-boosted Southern bright spot, indicating a black hole spin direction pointing away from the Earth [3]. After the observations in 2017 and publication of the initial total intensity in 2019, the M87 image in linearly polarized light was published in 2021, revealing a dynamically important poloidal B-field threading a plasma with a polarization pattern that spiraled azimuthally into the hole and electrons with inferred synchrotron temperatures from 10 to 120 billion K [4]. Following this milestone, this year, the EHT published an image of the supermassive black hole at our Galactic Center, likewise observed in 2017, revealing a ring with a face-on inclination of *i* < 50◦ and azimuthal hotspots in most reconstructions [5], affirming work done by the very large infrared telescopes of GRAVITY four years prior [6].

 The M87 image libraries produced by the EHT Collaboration's Theory and Simulations Team [3] strongly favored non-zero black hole spin values in order to exceed lower limits for the relativistic jet power in concordance with the Blandford–Znajek mechanism [8]. Some of the greatest uncertainties, however, lie in the interpretation of the heating mechanisms required to produce up-to-billion-degree bright features and the overall flux distribution of not only the emitting rings of M87 and Sgr A\*, but of jet/accretion flow/black hole, or JAB, systems in general. To this end, we have developed "Observing" JAB Systems to bridge state-of-the-art simulations and cutting-edge observations.

#### **2. Methodology**

*2.1. "Observing" JAB Systems*

The "Observing" JAB Systems pipeline can be summarized as follows:

Repeated applications of "Observing" JAB Systems to broad classes of phenomenological processes can naturally lead to model feature libraries with significant clusters in parameter space, as shown in our application to Sgr A\* [9]. Note the provisos in (2) and (3) are often overlooked without adherence to this methodology.1 We illustrate the importance of local piecewise modeling and the inclusion of significant radiating particle species, such as positrons, for our M87 application.

#### *2.2. GRHMD*

The first step in "Observing" JAB Systems makes use of the powerful simulated plasma physics laboratories produced by general relativistic magnetohydrodynamic (GRMHD) numerical methods. GRMHD methods are (typically) conservative in mass (*μ*(*ρuμ*) = 0) and stress–energy–momentum (*μTμ ν* = 0); 2 however, recent advances have enabled the inclusion of dissipative effects, such as viscous heating and heat conduction [11], as well as simulations of both thick torii and thin disks [12] and different implementations of the interaction between radiation and the background plasma fluid [13].

We use implementations of the HARM method [14,15] as a testbed for emission models. The simulations presented in this work were generated by using descendants of the harm code [14,15], a conservative second-order explicit shock-capturing finite-volume method for

solving the equations of ideal GRMHD in arbitrary stationary spacetimes. On a coordinate basis, the governing equations are

$$\partial\_t \left( \sqrt{-g} \rho u^t \right) = -\partial\_i \left( \sqrt{-g} \rho u^i \right), \tag{1}$$

$$\partial\_t \left( \sqrt{-\mathfrak{J}} T^t{}\_{\nu} \right) = -\partial\_i \left( \sqrt{-\mathfrak{J}} T^i{}\_{\nu} \right) + \sqrt{-\mathfrak{J}} T^{\mathbf{x}}{}\_{\lambda} \Gamma^{\lambda}{}\_{\nu\kappa\nu} \tag{2}$$

$$
\partial\_t \left( \sqrt{-g} \mathcal{B}^i \right) = -\partial\_{\bar{j}} \left[ \sqrt{-g} \left( b^j u^i - b^i u^j \right) \right], \tag{3}
$$

along with the constraint

$$
\partial\_i \left( \sqrt{-g} \mathcal{B}^i \right) = 0,\tag{4}
$$

where the plasma is defined by its rest mass density *ρ* and its four-velocity *uμ*, and *bμ* is the magnetic field four-vector. More details about the simulation process can be found in [16].

All simulations were initialized from a Fishbone–Moncrief (FM) torus [17] with a spin of *a* = 0.5 (Sgr A\*) or *a* = −0.5 (M87) (the latter is loosely based on a well-performing model with respect to the M87 polarization constraints in [4]). The Sgr A\* simulations were performed by using the iharm3d code [18]. The M87 simulations were produced by using iharm3d's kokkos/GPU-based descendent, kharma.3 Other physically motivated flow geometries are possible. The Bondi spherical accretion flow solution [19] is analytically tractable, but requires an exquisite degree of symmetry and non-magnetized flows to be realistic. Magnetized Bondi flows in numerical calculations display similar properties to those of torus simulations in terms of horizon-scale emission [20,21]. Stellar winds from a small population of ∼ 30 Wolf–Rayet stars were theorized, e.g., in [22–25] to source the disk and inflow of Sgr A\* in a model that reasonably accounted for the diffuse X-ray emissions. However, the precise knowledge of the position and orbits of this population is currently only possible in our own Galactic Center, while the FM torus model can be generalized to other AGNs. This generalization includes the possibility that the FM torus can be tilted with respect to the black hole spin (a situation that could arise often in a low-luminosity AGN), but we neglect this possibility in our work, which is focused on emission physics.

The initial magnetic field geometry is computed as the curl of a prescribed axisymmetric electromagnetic vector potential *Aφ*(*r*, *φ*), which is computed at simulation zone corners. In the SANE case, the vector potential is linear in the plasma density *Aφ* ∼ *ρ*; in the MAD case, *Aφ* ∼ *ρr*3*er* . The field strength is normalized so that the ratio of the maximum gas pressure to the maximum magnetic pressure over the domain is 100. The simulations are performed in a cylindrified, modified spherical Kerr–Schild coordinate system that concentrates the resolution near the midplane and close to the event horizon. The domain extends to at least *r*max = 1000 *GM*/*c*2 with a resolution of *Nr* × *Nθ* × *Nφ* = 288 × 128 × 128. More details about the simulation procedure and the initial conditions can be found in [16].

Magnetic flux is a key distinguishing factor among accreting plasmas. The timeaveraged magnetic flux Φ = (1/2) *θ ϕ* |*Br*|*dAθϕ*/ *mr*˙ *gc* determines two distinct regimes.

For Φ - 50, the disk is magnetically arrested (MAD) by its own magnetic pressure as it plunges into the hole. A much smaller Φ governs standard and normal evolution (SANE). Analyses of EHT simulation libraries tend to prefer MAD models over SANE, e.g., with MAD *a* = −0.5 earning the highest average image score for M87 with the parametric likelihood estimation procedure THEMIS [26].

#### *2.3. Plasma-Heating-Based Emission Models*

### 2.3.1. *R* − *β* Model

In general, GRMHD simulations evolve only the bulk of a fluid, i.e., the dynamically important ions.4 Therefore, in radiative post-processing, we seek to bridge the ions to the radiating electrons by the ratio of their temperatures *T*i/*T*e. In a hot, low-density, collision-less plasma, electrons can radiate and, therefore, cool efficiently, while the ion

cooling through Coulomb collisions is suppressed. The heating processes believed to be in action, e.g., viscous, compressional, and turbulent heating, can have similarly asymmetric effects on the temperatures and are more poorly understood [27]. Based on the tendency of plasma turbulence to preferentially heat electrons at low gas-to-magnetic pressure ratios *β* = *Pg*/*PB* and ions at high *β* [28,29], the *R* − *β* model is well motivated (e.g., [10,30–33]):

$$R = \frac{T\_i}{T\_c} = \frac{\beta^2}{1 + \beta^2} R\_{\text{high}} + \frac{1}{1 + \beta^2} R\_{\text{low}}.\tag{5}$$

The *R*-*β* model is the primary model used by the EHT Collaboration [3,34]. The electron-to-ion temperature ratio asymptotically reaches 1/*R*low in the low *β* regime and 1/*R*high in the high *β* regime. In turn, this means that *R*high controls the electron temperature in the disk or torus, and *R*low governs *T*e in the jet or wind outflow. Often, *R*low is fixed to 1 and only *R*high is varied, since by normalizing the flux to fit observations, the jet appears comparatively brighter upon an increase in *R*high, i.e., upon a decrease in *T*e in the disk (see also Figure 3 in [10]). The *R* − *β* model has been extensively compared to GRMHD simulations that readily employ heating models [35,36], concluding that on scales probed by the EHT, the *R* − *β* model approximates the influences of electron heating physics reasonably well.

#### 2.3.2. Critical *β* Model

This alternative turbulent heating model has an exponential parameter *βc* that controls the transition between electron- and ion-dominated heating [9]

$$\frac{T\_c}{T\_c + T\_i} = f e^{-\beta/\beta\_c}.\tag{6}$$

The exponential parametric control over the maximum *β* contributing to the emitting region is the basis for distinctly different near-horizon electron heating behavior compared to the R-*β* model.

We compare the R-*β* and Critical Beta models in Figure 1. For the same range of electron-to-proton temperature ratios, the Critical Beta model can have a sharper or smoother decline in the emission contributions from the highest *β* regions, smoothly transitioning from funnel/outflow- to near-horizon/inflow-dominant heating profiles as *βc* increases.

At low *β*, the electron-to-ion temperature ratio approaches parametrically determined maximum values (*f* and 1/*R*low) for the Critical Beta and *R*-Beta models, respectively. However, at high *β*, the electron temperature always asymptotes to 0 in the former model and is adjustable (through *R*high) only for the latter model. The exponential rate of electron temperature fall-off in the Critical Beta parametrization should have testably different spectral properties, such as the lowering of the bremsstrahlung contribution to the spectral energy distribution.

**Figure 1.** Electron-to-proton temperature ratio variation between *Te*/*Ti* = 10 and ∼0 for the R-Beta and Critical Beta Models as a function of *β*. The electron-to-ion temperature ratio varies as *Te*/*Ti* = 1 (1+*β*2)−1*R*low+(1+*β*2)−1*β*2*R*high for the R-*β* Model and *f e*−*β*/*βc* 1−*f e*−*β*/*βc* for the Critical Beta Model. For the R-Beta parameters (*R*low, *R*high)=(0.1, 100) and Critical Beta parameters (*f* , *βc*)=(0.918, 0.1) (**Left**), the models nearly coincide. For the R-beta parameters (*R*low, *R*high)=(0.1, 100) and Critical Beta parameters (*f* , *βc*)=(0.91, 1) (**Right**), the Critical Beta model has a softer transition.

#### *2.4. Sub-Equipartition-Based Models*

#### 2.4.1. Constant *βe*

One of the simplest yet most powerful models for understanding jet emission in JAB systems is the Constant Electron Beta Model [37,38]:

$$P\_{\mathfrak{c}} = \beta\_{\mathfrak{c}0} P\_{\mathfrak{B}}.\tag{7}$$

This model exploits our analytic knowledge of relativistic jets as being commonly generated by the Blandford–Znajek mechanism [8], which converts energy from accreting magnetized plasmas and spinning black holes into Poynting flux relativistic electromagnetic outflows with a power of *PBZ* ∝ *a*2Φ2 *B*. We now assume that a fixed fraction *βe*0 of jet magnetic energy is available for the emission observed at radio VLBI frequencies, where *βe*0 is related to the efficiency of conversion.

#### 2.4.2. Magnetic Bias

This class of model generalizes the Constant *βe* model by relating the pressure *Pe* of relativistic emitters (electrons or positrons) to the conversion of magnetic energy into particle energy through the powers *n* of the magnetic pressure:

$$P\_{\mathcal{C}} = \mathcal{K}\_{\mathcal{W}} P\_{\mathcal{B}}^{\mathcal{W}} \sim B^{\mathcal{Z}n}. \tag{8}$$

The constant *Kn* ensures that the right-hand side has units of pressure, and it can be estimated from a simulation by taking the average value of *BN*/*B*2 (where *N* = 2*n*) over a simulation surface enclosing the flux threading the black hole. The Bias Model parameter *n* helps modulate jet collimation in circumstances in which azimuthal magnetic fields are expected to scale simply with the cylindrical radius [37,39], e.g., the Blandford–Königl [40] model, where *B* ∼ *Bφ* ∼ *r* for radio jets.

#### *2.5. Hybrid Models*

We recognize the broad diversity of plasma regimes represented in a typical GRMHD field of view. Disk plasmas at large radii (a few tens or hundreds of *M*) are often not in inflow/outflow equilibrium [9], and thus, many of the thermal emission modeling approaches fail. As a Keplerian or sub-Keplerian disk approaches the horizon, its behavior near the plunging region beyond the innermost stable circular orbit (ISCO) is heavily dependent on SANE and MAD. For the MAD case, the disk can barely trickle into the supermassive black hole. For SANE flows, plasma can continuously flow. In both cases, the black hole spin can interact magnetohydrodynamically with disk plasma to form stable relativistic jets in simulations [41]. These jets, which are of underdense material themselves, are subjected to instabilities that are often distinct from those found in the disk, such as *m* = 0 pinch, *m* = 1 kink, and magnetic Kelvin–Helmholtz instabilities.

The inflow/outflow division in JAB systems gives an impetus for generating piecewise models. Jet plasmas are characterized by their low density and high energy, making the magnetization

$$
\sigma = \frac{b^2}{\rho},
\tag{9}
$$

which is a natural demarcation for the transition from jet-dominated emission to accretionflow-dominated emission. When jets are well collimated, we may also use parabolic or other geometric cuts to isolate the jet region [37,39,42,43].

#### *2.6. Phenomenological Models*

We postulate other mechanisms here (detailing their full functional form in the model compendium in Appendix A.1) by noting the synchrotron emissivity *jν* ∝ *Pe*, where the pressure of relativistic emitters can be written as *Pe* ∝ *Wt*cool in terms of the dissipation rate per unit time *W* and the cooling time *t*. In the following, we proceed to specify more potential models in our arsenal, systematically carrying out the second step of "Observing" JAB Simulations by relating the pressure of relativistic emitters to energetic processes in an AGN.

The Current Density Model relates the dissipation of energy into emitting particles to the current density *W* ∝ *j* 2. This is seen to trace a co-axial current morphology with a central outgoing current and a return current, thus creating a boundary layer in jets [39].

The Jet Alpha Model parametrizes the efficiency of linear momentum transport in jets by *αJ*, i.e., *W* ∝ *αJ*, in a manner analogous to that of the Shakura–Sunyaev model for angular momentum transport in accretion disks [44]. This dissipation rate is also linear in the shear stress, as seen in Appendix A.1, thus enabling the disk–jet interface to be visible [39].

Lastly, the Shear Model adopts a Newtonian framework for velocity shear *τ* = *μS*, where *μ* is the dynamical viscosity and *S* is the shear stress (|*dvz*/*ds*| in cylindrical coordinates). Then, *W* ∝ *τS* ∝ *S*2. This quadratic dependence on shear stress gives us an edge-brightened model relative to the Jet Alpha Model.

These jet models may naturally be glued to disk and corona models (especially those based on turbulent heating), as they capture the behavior of the inflow/outflow interface. Together, these phenomenological emission prescriptions may form the building blocks for detailed hybrid models that account for different plasma flow physics throughout the JAB system.

#### *2.7. Electron Distribution Functions*

A fixed temperature for electron thermodynamics models is often an idealization when particles are found in nearly collisionless plasmas (such as in JAB systems with mean free paths of O(105*M*)). The influence of a non-thermal electron population on horizon- and jet-base-scale emission remains an important subject of investigation in the EHT Collaboration [3,7]. These non-thermal distributions need not be applied globally in the GRMHD domain. Rather, determining the proper region's home for charged particle acceleration processes is an integral part of these studies. We explore some possibilities for particle acceleration and the concomitant energy distributions below.

#### 2.7.1. Power Law

One of the simplest phenomenologically viable assumptions for the energy distribution of a population of *ne* emitting particles in an astrophysical plasma is a power-law energy decay:

$$\frac{d\eta\_{\varepsilon}}{d\gamma} = K\gamma^{-p} \tag{10}$$

The normalization factor *K* depends on the synchrotron pressure. Power-law particle distributions are naturally produced in shock waves through diffusive shock acceleration, whereby particles are energized by repeated interactions with magnetic inhomogeneities as they propagate alongside the shock. This is a first-order Fermi process, as the energy gain of the particles is linear in the shock velocity. Astrophysical shocks may occur at the interface of fluids with differing velocities, such as the interface of a jet with its ambient medium.

#### 2.7.2. The Kappa Model

It is often preferable to be able to model the full SED with a single distribution function. This can be achieved with the kappa electron energy distribution function (e.g., [45,46]), which has its theoretical foundation in non-extensive Tsallis statistics [47,48]. Looking like a thermal distribution at low energies, it smoothly transitions into a non-thermal power-law tail with index *s* towards high energies, so that *κ* = 1 + *s* (see, e.g., Figure 4 in [10]). The thermal [49] and kappa [50] electron energy distribution functions read

$$\frac{dn\_{\mathbf{c}}}{d\gamma\_{\mathbf{e}}d\cos\overline{\xi}\_{\mathbf{e}}d\phi} = \begin{cases} \frac{n\_{\mathbf{e}}}{4\pi\Theta\_{\mathbf{e}}} \frac{\gamma\_{\mathbf{e}}\left(\gamma\_{\mathbf{e}}^{2}-1\right)^{1/2}}{K\_{2}\left(1/\Theta\_{\mathbf{e}}\right)} \exp\left(-\frac{\gamma\_{\mathbf{e}}}{\Theta\_{\mathbf{e}}}\right), & \text{thermal} \\\\ \frac{N}{4\pi}\gamma\_{\mathbf{e}}\left(\gamma\_{\mathbf{e}}^{2}-1\right)^{1/2} \left(1+\frac{\gamma\_{\mathbf{e}}-1}{\kappa w}\right)^{-(\kappa+1)}, & \text{kappa} \end{cases}$$

where *n*e is the electron number density, *φ* is the gyrophase, *γ*e is the electron Lorentz factor, *ξ* is the electron pitch angle, and *K*2 is the modified Bessel function of the second kind. From these distribution functions, emission and absorption coefficients are determined by using fit functions [51]. *w* is the "width" of the distribution and describes the energy in the system.

From particle-in-cell (PIC) simulations of magnetized current sheets, it is evident that the kappa index is not constant in all sub-regions of the system, requiring kappa to be variable and dependent on plasma quantities [52]. Further, since many GRMHD codes do not provide accurate values close to the jet spine due to boundary conditions, the inner spine is usually excluded from the emission by imposing a maximum in the magnetization *σ*cut (e.g., [10,31,32]).

Additionally, the distribution function can be modified to account for a thermal and a magnetic contribution to the total energy [10,31,32,53]. That way, it is possible to control the amount of magnetically accelerated electrons and the distance of their point of injection into the jet from the central engine.

It is important to distinguish the influences of the individual dials of this non-thermal model on the image morphology. This requires extensive parameter surveys, which would be too computationally expensive in a two-temperature GRMHD simulation. Therefore, an implementation in radiative post-processing is currently the only feasible option for these surveys [10]. Increasing *ε* or *σ*cut raises the SED at energies past the synchrotron turnover, with a growing influence towards the highest energies (e.g., Figure 12 in [10]).

#### *2.8. Emission Modeling in Non-Kerr Spacetimes*

For over a century, Einstein's theory of general relativity has been thoroughly tested in many different ways (e.g., [54–63]). In the era of the EHT and ngEHT, a test of GR in an imaging-based approach is on the horizon for the first time. Specifically, the shape and size of the black hole shadow and photon ring are crucial characteristic properties of horizonscale images. Photons rings in the Kerr metric are predicted to be nested in subrings with exponentially decreasing separation, the lowest orders of which may be observed with the aid of longer baselines in the ngEHT [64]. In addition to this observational test of the Kerr solution for the spacetime around astrophysical black holes, many studies on alternative theories to general relativity make use of semi-analytical models for both plasma accretion and emission processes to compare models to observations (e.g., [65–67]). Alternative spacetime geometries are rarely investigated in full GRMHD and GRRT, and if so, the emission physics are based on a constant proton-to-electron temperature ratio and purely thermal radiation [68,69]. Only recently have advances been made to study the influence of emission models more akin to reality in an alternative spacetime [33,53]. The fundamental difficulties in finding deviations from the Kerr metric arise from the presence of greater astrophysical uncertainties. For example, the magnetic field configuration in a GRMHD simulation appears to have a much larger influence on the source morphology on horizon scales than the background spacetime [33,53]. Moreover, the differences in image features, such as the shadow size or photon sub-ring spacing caused by a deviation from GR, are often small and subject to degeneracy with accretion and emission models. Still, past EHT observations of black hole shadows have helped constrain alternative theories to GR [65,70]. With the help of the ngEHT, we aim to resolve the degeneracy between the effects of plasma physics and GR and put even more robust constraints onto alternative spacetime geometries.

#### **3. Commencing the Computing: Emission Models in Numerical Codes**

The emission models discussed above have been implemented in a variety of numerical codes (e.g., [71–77]; for an extensive comparison, see [78]). Since GRMHD simulations are already computationally expensive, radiation is commonly modeled in post-processing. While there are a handful of radiative GRMHD (GRrMHD) codes (e.g., [36]), when evolving a two-temperature plasma while accounting for ions and electrons, the computational cost only increases. In the case of pure jet simulations, however, it is possible to bring the cost down by moving to special relativity (e.g., [79–83]).

On event horizon scales, i.e., in the strong-gravity regime, we are required to take the full effects of GR into account. Usually, radiative transfer codes first calculate light rays by integrating the geodesic equation and, subsequently, solve the radiative transfer equation along those rays. The use of post-processing enables us to freely investigate a variety of emission models, with many going beyond thermal synchrotron radiation. For example, Compton/inverse Compton scattering, bremsstrahlung, and non-thermal emission processes all have their own imprints on the SED and the image morphology in both total intensity and polarized light [43,76].

#### *Positrons' Effects on Radiative Transfer*

Positrons in JAB systems can be produced through photon–photon (Breit–Wheeler) interactions in jet funnel walls [84] and spark gap processes near the magnetospheric poles of supermassive Kerr black holes [85]. Each of these processes may contribute electron– positron pair densities exceeding the Goldrich–Julien value [86], which is required in order to screen the large-scale electric fields responsible for high-energy lepton cascades. In sources such as M87 that have significant pair production through these channels, positron effects abound. These include the increase in linear and vanishing of circular polarization, as well as the higher energy fall-off of the circular polarization spectrum, all of which have been modeled as potential discriminators of JAB systems that are rich in ionic versus pair plasma [43,83].

#### **4. Results: Adding an Observer**

We present a suite of parametric emission models and illustrate them in GRMHD simulations that are scaled to the first horizons observed: M87\* and Sgr A\*. We start by painting a single semi-MAD Sgr A\* simulation with several phenomenologically moti-

vated prescriptions of electron thermodynamics. Thus, by using the second step of the "Observing" JAB Simulations pipeline, we illustrate the sensitivity of the interpretations of near-horizon structure in real-world observations to our models of plasma physics. We then use a SANE and MAD M87 simulation to reveal widely divergent—and observationally distinguishable—polarized emission signatures governed by Faraday and positron effects. The latter effects enable us to better constrain a fundamental, though largely uncertain, property of relativistic jets: ion- versus lepton-dominated composition. In what follows, we use perceptually uniform color maps to identify image morphology with the underlying emission mechanisms that energize relativistic particles to produce synchrotron radiation.

#### *4.1. Sgr A\**

#### 4.1.1. Parametric Model Comparison

We explore the Critical Beta model parameter space from *f* ∈ {0.1, 0.5} and *βc* ∈ {0.01, 0.1, 1}, starting with an edge in view to highlight the lensing effects on image morphology in Figure 2. As the *f* parameter increases, the overall electron temperature increases, and changes in the image morphology at fixed flux are effected through the *M*unit used in codes to scale the relative importance of inertial plasma properties, such as the density and mass accretion rate, relative to the plasma's electromagnetic properties. As the critical beta parameter increases, the locus of the greatest electron contributions moves from the low *β* outflow to higher *β* regions that are lensed around a compact crescent near the horizon due to higher values of *β* being found in the inflow around the black hole. Thus, *βc* serves as a dial in JAB emission modeling for compactifying emitting regions and asymmetrizing them from an edge-on orientation, as shown in Figure 2.

**Figure 2.** *Cont*.

**Figure 2.** A parameter scan of the Critical Beta Model with (**Left**) *f* = .1, (**Right**) *f* = 0.5, *βc* = 0.01 (**Top**), *βc* = 0.1 (**Middle**), and *βc* = 1 (**Bottom**) from an edge-on view. For Sgr A\* models, the cgs conversion into Jy is found by multiplying each cgs-intensity-colored pixel value by 57.9 to get its flux density in Jy [9]. The scale is *M* ≡ *GM*BH/*c*2 when used as a length and ≡ *GM*BH/*c*3 when used as time. The notation 1e2 is a compact form of the scientific notation 1 × 102.

In Figure 3, we see different morphologies associated with our equipartition-inspired Constant *βe* and Magnetic Bias models. The locus of emission for the Constant Electron Beta Model approaches the funnel for low *βe*0 and broadens into a thick, lensed torus for higher *βe*0. When the bias parameter goes to 0, the pressure from relativistic electrons goes to a constant and does not sharply decline with the radius, leaving extended outflow signatures.

**Figure 3.** *Cont*.

**Figure 3.** A parameter scan of the Constant Electron Beta Model (**Left**) with *f* = 0.1, (**Right**) *f* = 0.5, *βe*0 = 0.01 (**Top**), *βe*0 = 0.01 = 0.1 (**Middle**), and *βe*0 = 1 (**Bottom**). Magnetic Bias Model with *βe*0 = 1 and *N* = 0 (**Bottom Right**).

We now change the orientation and consider face-on models in Figure 4. An Sgr A\* 20◦ spin axis orientation was observationally preferred by GRAVITY [6] , and this was corroborated by the EHT, disfavoring inclination angles above 50◦. Here, tuning up the critical value of beta, *βc*, still leads to more compact images, but now they maintain a ring symmetry even near the gravitational lensing profile of a Kerr black hole.

**Figure 4.** *Cont*.

**Figure 4.** A parameter scan of the face-on Critical Beta Model with (**Left**) *f* = 0.1, (**Right**) *f* = 0.5, *βc* = 0.01 (**Top**), *βc* = 0.1 (**Middle**), and *βc* = 1 (**Bottom**).

#### 4.1.2. Morphological Classification

The ring morphology has dominated the focus of emission modeling since the advent of the first two ring-like images of the horizon scale were released by the EHT in 2019 and 2022. However, we showed examples of the strong inclination dependence of 230 GHz images above. It is worth taking stock of how fortuitous the near-face-on (∼20◦) spin axis viewing angles of Sgr A\* and M87\* are given that their main selection criterion is the exceptionally large angular width of their horizon gravitational radii as seen from Earth, which is independent of their spin inclination angle.

We preview the edge-on morphologies that the ngEHT will see based on morphological clusters in the parameter space of the Critical Beta, Constant Electron Beta, and Magnetic Bias Models in Figure 5. These types are:

In Figure 6, we plot the morphological types on a common intensity scale to emulate observing with a single instrument. The EHT's dynamical range of ∼2 orders of magnitude in intensity is used to set this scale. Morphological degeneracies, e.g., between Types I (large *f* , *βc*) and III (large *βe*0) and between Types II (small *f* , *βc*) and IV (small *n*), are more likely to emerge as more regions fall below the flux threshold of the observing instrument. The ngEHT's dynamical range will span a few orders of magnitude. The ngEHT's improved sensitivity down to ∼5 mJy and increased frequency range to -345 GHz [87] will enable us to resolve the currently excluded low-flux density regions in the field of view and at larger radii. This will lead to more accurate determination of the morphological type, thus breaking degeneracies. Thus, as the ngEHT increases the dynamic range of JAB system observations to lower fluxes, we are able to better classify images and pinpoint the underlying emission physics related to their morphologies.

**Figure 5.** Semi-MAD simulation [9] models that were ray traced at 230 GHz at *T* = 10,000 *M*: (**Top Left**) best-fit Critical Beta model (**Top Right**) R Beta with (*f* , *βc*)=(0.5, 1); (**Bottom Left**) Constant Electron Beta model with *βe*0 = 1; (**Bottom Right**) Magnetic Bias with *N* = 0 jet.

**Figure 6.** Semi-MAD simulation [9] model classes that were ray traced at 230 GHz at *T* = 10,000 *M*, as shown in Figure 5, but with a common intensity scale.

#### *4.2. M87*

4.2.1. Parametric Model Comparison

We implement hybrid models that include turbulent heating and sub-equipartition emitting regions for M87. We make a natural choice to partition the simulation region according to the magnetization (i.e., the magnetic energy density to enthalpy density ratio *σ* = *b*2 *ρ* ). The value *σ*transition = 1/2 that determines the jet regions is near the inflow/outflow interface and is less than the *σ*cut of 2 for these simulations.

A SANE/MAD dichotomy emerges in Figures 7 and 8. The SANE intensity maps are more ring symmetric, whereas the MAD case has a prominent flux loop (lower left of all panels in Figure 8). The SANE circular polarization varies on smaller spatial scales than the MAD does. The sign of the circular polarization changes on sub-*M* scales for SANE. The SANE and MAD hybrid models that have added a Constant Electron Beta jet for *σ* > *σ*transition have broader flux distributions over the field of view (all of the right panels in Figures 7 and 8).

**Figure 7.** Synthetic intensity with the electric vector polarization angle (EVPA) and circular polarization maps for models that include positron effects and piecewise modeling. For the *a* = −0.5 SANE at 230 GHz and at *T* = 25,000 *M*: (**Top Left**) R-Beta with a *βe*0 = 0.01 jet model; (**Top Right**) R-Beta with a *βe*0 = 0.01 jet; (**Bottom Left**) Critical Beta Model; (**Bottom Right**) Critical Beta with a *βe*0 = 0.01 jet. For each case, the intensity is overplotted with the electric vector polarization angle on the left panel, and the circular polarization degree is mapped on the right panel.

#### 4.2.2. Positron Effects

We see the effect of the addition of positrons in Figure 9. There, we start with an ionic plasma with a number density of *n*0 (see the top panels), and then add enough positron pairs *n*+ so that the pair number density equals that of the original electrons (and the unit for *M* is adjusted to match the normalized flux (see the bottom panels)). In Figure 9, we see that in the MAD case (left), the degree of circular polarization is proportional to the unpaired emitter fraction (1/3 in the case of *f*pos = 1) [43,88]. The SANE case (right), which has Faraday rotation depths that are thousands of times greater than the less dense, more highly magnetized MAD, does not have a simple linear relationship between the pair content and *V*/*I*. In the linear polarization, we also see greater positron effects for SANE than for MAD, as the addition of positrons scrambles the EVPA pattern only in the SANE case.

**Figure 8.** For the *a* = −0.5 MAD at 230 GHz and at *T* = 25,000 *M*: (**Top Left**) R-Beta Model; (**Top Right**) R-Beta with a *βe*0 = 0.01 jet; (**Bottom Left**) Critical Beta Model; (**Bottom Right**) Critical Beta with a *βe*0 = 0.01 jet. For each case, the intensity is overplotted with the electric vector polarization angle on the left panel, and the circular polarization degree is mapped on the right panel.

**Figure 9.** For *a* = −0.5 at 230 GHz and at *T* = 25,000 *M*: (**Top**) Critical Beta Model without positrons; (**Bottom**) Critical Beta with *n*pairs/*n*0 = 1.We also compare MAD (**Left**) and SANE (**Right**).

#### **5. Discussion and Conclusions**

We have seen from a SANE and a MAD M87 simulation and one semi-MAD simulation for Sgr A\* that the thermodynamic properties of emitting particle populations can be as important as those of the underlying plasma for predicting the properties of intensity and polarization maps near the horizon. Moreover, special relativistic effects, such as beaming and general relativistic lensing, can produce nonlinear modifications of the flux received by observers away from the JAB system. Understanding the relative contributions of these effects can be facilitated by idealized cases, such as in semi-analytic jet models, e.g., those in [42]. We considered the impacts of some of the competing effects listed here:

By using a variety of turbulent heating and magnetic-to-particle-energy-ratio-based emission models, we found a number of distinct image morphologies that we expect to encounter in future observational campaigns with the ngEHT. In particular, the ngEHT will have the resolution to probe images with areas a few times smaller than those probed with the EHT, allowing the imaging of more compact ringlike structures that are viewed face-on (corresponding to high values of Critical Beta, should this turbulent heating mechanism predominate). This will enable us to better determine which processes dominate the emissions. Furthermore, we showed the strong dependence on inclination of the relationship between electron thermodynamics and image morphology. We identified at least four distinct image signatures of turbulent heating and magnetic-to-particle energy ratios, mainly visible in the edge-on view. We also showed the EHT's/ngEHT's ability to discriminate between MAD and SANE accretion states based on ordered vs. scrambled polarization signatures. This dichotomy is accentuated by the presence of positrons, allowing us to probe the plasma composition of JAB systems, which remains a long-standing question.

#### *5.1. Limits of Instrumentation*

The 20 μas angular resolution of the EHT (enough to resolve an atom at arm's length or an apple on the Moon) is exquisite according to the VLBI observational standards. More remarkably still, the ngEHT may double this angular resolution. However, the angular gravitational radius of even the largest black hole shadow (*M* = 5 μas for Sgr A\*) is still smaller than the size of the synthesized beam of the ngEHT (20 μas). This leads to blurring in the image plane, washing out intricate features. Due to this instrumental limitation, numerical simulations are crucial for our understanding of the properties of accretion flows. Though simulations themselves produce grid-resolution-dependent models (cf. Appendix A.2), they remain our most reliable probe below the scale set by the gravitational radius, and will continue to motivate observing missions for the ngEHT and beyond to confirm their predictions.

#### *5.2. Universality of Select Measures*

The existence of flaring emissions in Sgr A\* and extended synchrotron emissions in M87's jet provide strong evidence for the presence of a high-energy tail to the lepton distribution, i.e., non-thermal lepton distribution functions. This is generally expected as a consequence of the microphysical processes responsible for the dissipation of turbulence and the attendant injection of energy explored above. However, the detailed shape and magnitude of this additional tail depend sensitively on the specific acceleration mechanism considered. Nevertheless, for the physical parameters that are relevant for Sgr A\* and M87, below 690 GHz, the synchrotron coefficients are only weakly sensitive to the particular choice of the extension of the lepton distribution function [89].

The lepton distributions are expected to be rapidly isotropized on the cyclotron scale by plasma instabilities [14,15]. For such distributions, the synchrotron emission and absorption coefficients for all Stokes parameters may be expressed as convolutions with broad kernels in frequency space [90]. As a consequence, all of these transfer coefficients are well approximated by a universal expression that is dependent only on the local spectral index, which conceptually corresponds to a measure of the relative number of "hot" non-thermal and "cold" thermal leptons [89].

At a single observation frequency, this approximation is as good as to better than 2%. Thus, for images at 230 or 345 GHz, for example, the specific nature of the acceleration mechanism may be effectively and efficiently parameterized by a single lepton distribution model (e.g., thermal, *κ*, power law, etc.), thus eliminating a key systematic degeneracy between the location of and microphysical processes responsible for turbulence dissipation. Across large frequency ranges, e.g., from 230 to 345 GHz, the accuracy of this approximation falls to ∼10%, or from 230 to 480 GHz to ∼40%. Therefore, multi-frequency image reconstructions remain a powerful discriminant between different acceleration mechanisms.

**Author Contributions:** Conceptualization, R.A.; methodology, R.A.; software, A.R., R.E.; validation, J.D., L.O. and N.N. (Nathan Ngata); formal analysis, J.W.; investigation, R.A.; resources, A.R.; data curation, B.C. and S.R.; writing—original draft preparation, R.A., J.R., A.E.B. and G.N.W.; writing review and editing, J.R., J.D., L.O. and N.N. (Nathan Ngata).; visualization, J.D., L.O., N.N. (Nathan Ngata), N.N. (Nitya Nigam) and E.D.; supervision, R.A.; project administration, R.A.; funding acquisition, R.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** Jan Röder received financial support for this research from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. Razieh Emami acknowledges the support from the Institute for Theory and Computation at the Center for Astrophysics, as well as grant numbers 21-atp21-0077, NSF AST-1816420, and HST-GO-16173.001-A, for their very generous support.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A**

*Appendix A.1. List of Emission Models*

We summarize the emission models employed above (and more) in Table A1.

**Table A1.** JAB Emission Model List. This is a compendium of the phonologically motivated emission models mentioned in this work. The shear stress *S* = *γ*2|*dvz*/*ds*| and dimensional parameters *Lj* and *LS* are set by the width of jet system.

**Table A1.** *Cont.*

*Appendix A.2. Resolution Dependence*

In Figure A1, we briefly consider the robustness of image properties to resolution effects. Changing the resolution can shift the local and global maxima, as seen in these 50 *M* × 50 *M* fields of view, where going from 100 pixels per side to 500 pixels per side makes the maximum intensity (red) pixels more concentrated in a more compact equatorial region. Increasing the resolution also increases the overall maximum intensity, as one may expect due to the breaking up of the high-intensity regions into more pixels over which to perform the maximization. In both cases shown, as discussed in Section 5.1, the sub-gravitational-radius pixel spacing in the synthetic images exceeds current and planned observational resolution limits. However, for even the modestly sized supermassive black hole Sgr A\*, *M* = 1010 cm, so the role of emission modeling to fill in the sizable gap left by subgrid physics will remain paramount for years to come.

#### **Notes**

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## **Accretion Flow Morphology in Numerical Simulations of Black Holes from the ngEHT Model Library: The Impact of Radiation Physics**

**Koushik Chatterjee 1,2,\*, Andrew Chael 3, Paul Tiede 1,2, Yosuke Mizuno 4,5,6, Razieh Emami 2, Christian Fromm 6,7,8, Angelo Ricarte 1,2, Lindy Blackburn 1,2, Freek Roelofs 1,2, Michael D. Johnson 1,2, Sheperd S. Doeleman 1,2, Philipp Arras 9,10, Antonio Fuentes 11, Jakob Knollmüller 10,12, Nikita Kosogorov 13,14, Greg Lindahl 2, Hendrik Müller 8, Nimesh Patel 2, Alexander Raymond 1,2, Efthalia Traianou 11 and Justin Vega 2,15**

**Abstract:** In the past few years, the Event Horizon Telescope (EHT) has provided the first-ever event horizon-scale images of the supermassive black holes (BHs) M87∗ and Sagittarius A∗ (Sgr A∗). The next-generation EHT project is an extension of the EHT array that promises larger angular resolution and higher sensitivity to the dim, extended flux around the central ring-like structure, possibly connecting the accretion flow and the jet. The ngEHT Analysis Challenges aim to understand the science extractability from synthetic images and movies to inform the ngEHT array design and analysis algorithm development. In this work, we compare the accretion flow structure and dynamics in numerical fluid simulations that specifically target M87∗ and Sgr A∗, and were used to construct the source models in the challenge set. We consider (1) a steady-state axisymmetric radiatively inefficient accretion flow model with a time-dependent shearing hotspot, (2) two time-dependent single fluid general relativistic magnetohydrodynamic (GRMHD) simulations from the H-AMR code, (3) a two-temperature GRMHD simulation from the BHAC code, and (4) a two-temperature radiative GRMHD simulation from the KORAL code. We find that the different models exhibit remarkably similar temporal and spatial properties, except for the electron temperature, since radiative losses substantially cool down electrons near the BH and the jet sheath, signaling the importance of radiative cooling even for slowly accreting BHs such as M87∗. We restrict ourselves to standard torus accretion flows, and leave larger explorations of alternate accretion models to future work.

**Keywords:** black holes; general relativity; accretion; relativistic jets; very-long-baseline interferometry

**Citation:** Chatterjee, K.; Chael, A.; Tiede, P.; Mizuno, Y.; Emami, R.; Fromm, C.; Ricarte, A.; Blackburn, L.; Roelofs, F.; Johnson, M.D.; et al. Accretion Flow Morphology in Numerical Simulations of Black Holes from the ngEHT Model Library: The Impact of Radiation Physics. *Galaxies* **2023**, *11*, 38. https://doi.org/10.3390/ galaxies11020038

Academic Editor: Dimitris M. Christodoulou

Received: 4 December 2022 Revised: 14 February 2023 Accepted: 20 February 2023 Published: 23 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

With the advent of the Event Horizon Telescope ([EHT; [1,2]), imaging the near-horizon structure of supermassive black holes (BHs) is now a reality. The primary targets of the EHT and the future next-generation EHT (or ngEHT1) are M87∗ (the supermassive BH in the elliptical galaxy M87; [1]) and Sgr A∗ in the Galactic Center [2], which are two of the most well-studied low-luminosity active galactic nuclei. Extracting information about the event horizon-scale accretion flows in these two sources using the EHT's enormous resolving power is an active area of research. With the ngEHT, we will achieve unprecedented levels of angular resolution and sensitivity to low-flux regions, with the dynamic range in flux expected to increase to ∼1000 compared to the EHT's current dynamic range of ∼10 (e.g., [3]). This would enable us to investigate the BH shadow shape with higher precision as well as provide a crucial connection between the accretion flow and the jet launching region. The expected advances in sensitivity require deeper investigations of feature extraction from simulated synthetic reconstructions of BH systems. Hence, we designed the ngEHT analysis challenges 2 [4] to test our ability to capture the complex dynamics of gas and magnetic fields around M87∗ and Sgr A∗ using the ngEHT reference array (e.g., [5]) with various analysis methods.

Black hole accretion and jet physics have been intensively studied over the past few decades (e.g., [6–15]). In the context of M87∗ and Sgr A∗, we expect the accretion flow to be highly sub-Eddington, radiatively inefficient, and geometrically thick, popularly known as radiatively inefficient accretion flows (RIAFs). This accretion flow solution has been used to successfully model the multiwavelength spectrum of Sgr A∗ (e.g., [16]). On the other hand, semi-analytical models of jets are preferred to explain the spectrum of M87∗ (e.g., [17]). Thus, these two sources already provide a means to probe two different components of BH accretion, namely, the inner accretion flow structure and turbulence in Sgr A∗ and the prominent jet feature in M87∗. The first three EHT numerical simulation papers [18–20] already provide us with important clues about the horizon-scale conditions of these BH systems based on numerical simulations: (1) these BHs probably have non-zero spin, (2) the accretion disk is expected to have colder electrons than the jet sheath, and (3) the observations favor the presence of dynamically important magnetic fields close to the BH. All of these results point us toward the magnetically arrested disk (MAD; [10,21]) state, an accretion mode where the BH magnetosphere becomes over-saturated with magnetic flux and exhibits quasi-periodic explosions of vertical magnetic field bundles. MAD flows around spinning BHs also have powerful relativistic jets, where the jet power can exceed the input accretion power [13], which is a definite signature of BH spin energy extraction via the Blandford and Znajek [22] process.

Building on the semi-analytical RIAF models, the time-dependent general relativistic magneto-hydrodynamic (GRMHD) simulations have become important tools for deciphering BH accretion physics in a variety of astrophysical systems (e.g., [11,13,23–27]). Indeed, the EHT regularly uses the libraries of GRMHD simulations to model the observed horizonscale BH images of M87∗ and Sgr A∗ as well as larger-scale jet images (such as for Centaurus A; [28]) in order to constrain the time-variable plasma properties. In designing the ngEHT reference array, it is therefore crucial to use GRMHD simulations for understanding the attainability of specific science goals, such as resolving the photon ring and the disk-jet connection region as well as tracing out time-variable features via the ngEHT analysis challenges.

In this work, we discuss the numerical fluid simulations that were used as source models for the ngEHT analysis challenges. In particular, our objective is to compare between the models that incorporate increasingly complicated levels of accretion and electron physics, focusing on M87∗ and Sgr A∗. Our model set consists of a time-dependent shearing hotspot stitched to a steady-state RIAF solution, two standard GRMHD simulations of MAD accretion flows, a GRMHD MAD simulation with electron heating via incorporating twotemperature physics, and a fully radiative, two-temperature, GRMHD MAD simulation. We describe the equations and setup of our numerical models in Section 2, show our comparison results in Section 3, and, finally, conclude in Section 4.

#### **2. Numerical Simulations**

In this section, we provide a brief description of the semi-analytical stationary RIAF and shearing hotspot model as well as the (two-temperature/radiative) GRMHD simulations used for the ngEHT analysis challenges.

#### *2.1. RIAF + Hotspot Solutions*

RIAF models attempt to describe the time and azimuthally averaged appearance of accretion flows. This is performed using a set of building blocks. The specific RIAF models used in the challenges are based on Yuan et al. [16], Broderick and Loeb [29], Broderick and McKinney [30] approach. We decompose the accretion flow into a set of phenomenological models that describe the electron density, temperature, magnetic field, and velocity profile. We have a cylindrical coordinate system *xμ* ≡ (*R*cyl, *ϕ*, *z*). The electron density profile is defined in terms of the cylindrical radius *R*cyl = *r*| sin(*θ*)| and vertical displacement *z* = *r* cos(*θ*), and is provided by

$$n\_{\mathfrak{e},\mathfrak{X}}(R\_{\text{cyl}},z) = n\_{\mathfrak{e},\mathfrak{X},0} R\_{\text{cyl}}^{p\chi} \exp\left(-\frac{z^2}{2h^2 R\_{\text{cyl}}^2}\right) \tag{1}$$

where *X* denotes the population of electrons. The disk height *h* is set to unity for this work. For the challenge dataset, we included both thermal synchrotron emitting (*X* ≡ th), and non-thermal synchrotron (*X* ≡ nth) emitting electrons. The thermal electrons have *ne*,th,0 = 1.3 × 108 and *p*th = −1.1, while the non-thermal electrons are provided by *ne*,nth,0 = 1.3 × 105, *p*nth = −2.02. These numbers are from Tiede et al. [31] and are set to match the best-fit parameters for the Sgr A∗ spectrum from Broderick et al. [32] and Broderick et al. [33].

The temperature profile of the thermal electrons is also provided by a radial power law with a Gaussian envelope describing the height:

$$T\_{\varepsilon}(t, r, \theta, \varphi) = T\_{\varepsilon, 0} R\_{\text{cyl}}^{-0.84} \exp\left(-\frac{z^2}{2h^2 R\_{\text{cyl}}^2}\right),\tag{2}$$

where, for the challenge data, we set *Te*,0 = 6.3 × 1010 K.

Following Pu et al. [34], we define the gas pressure by assuming that the ultra-relavistic protons are in roughly virial equilibrium with the gravitational force providing:

$$p\_{\rm gas}(\mathbf{x}^{\mu}) = \frac{m\_p n\_{c, tl}(\mathbf{x}^{\mu})}{6} \frac{GM\_{\rm BH}}{r} \,\mathrm{},\tag{3}$$

where *mp*, *G*, and *M*BH are the proton mass, gravitational constant, and the black hole mass.

For the local magnetic field, we then assume a constant *β* = *p*gas/*p*mag = 10 plasma, which, combined with Equation (3), provide us the magnetic field strength. The orientation of the magnetic field is then provided by a purely toroidal configuration, relative to the plasma observer. Finally, we extract the emission and absorption coefficients from the synchrotron self-absorption model in Broderick and Blandford [35] for both the thermal and non-thermal synchrotron electrons. For the non-thermal electrons, we use a power law prescription with radiation coefficients from Jones and O'Dell [36], and a photon spectral power law index of 1.25 (see [31] for more details).

We follow the velocity field prescription from Pu et al. [34] to describe the accretion flow dynamics. Using the notation from Tiede et al. [31], the velocity field *uμ* = *ut vμ* is provided by