In the millimeter, we observe Sgr A\* and M87\* near the peak of emission from synchrotron radiation, where the flow transitions from optically thick to optically thin. Here, the emissivity scales approximately as *jν* ∝ *nB*2Θ5/2 *e* , where *n* is the electron number density, *B* is the magnetic field strength, and Θ*e* is the electron temperature2 (in units of the electron

rest mass energy). Each of these quantities can vary by orders of magnitude among different models, and thus even the total flux is informative for jointly constraining these parameters. To match the total flux of a given system, the fluid in ideal GRMHD simulations can be rescaled via *n* → M*n*, *B* → √M*B*, and *u* → M*u*, where *u* is the internal energy and M is a scalar. After doing so, both MAD and SANE simulations are capable of matching the total flux of EHT sources at a single frequency, as well as broad image characteristics such as the image size [5,13]. However, this rescaling causes SANE simulations to typically have orders of magnitude larger number density than MADs, due to their intrinsically weaker magnetic fields and lower temperatures [8,13]. Consequently, any additional observables sensitive to these variables immediately help break degeneracies and distinguish between models.

For example, the degeneracies between *n*, *B*, and Θ*e* can be partially resolved with the spectral index, *α* ≡ *d* log *I*/*d* log *ν*. Spectral index is mainly sensitive to the optical depth *τν* as well as the temperature and magnetic field strength in the combination *B*Θ2 *e* (on which the critical synchrotron frequency is dependent), e.g., [45]. Ricarte et al. [41] show that GRMHD models span a wide range of spectral indices, and that SANE models typically exhibit more negative spectral indices than MADs at a fixed optical depth due to their lower temperatures. Although lacking by construction in ideal GRMHD simulations, spectral index can also provide insight into synchrotron cooling processes that affect the temperature.

In the ngEHT era, multi-frequency VLBI will enable spatially resolved spectral index maps. Since the most important parameters (Θ*e*, *B*, and *τν*) all decline with radius, spectral index maps should generically grow more negative as radius increases. Equivalently, the image becomes smaller as the frequency grows larger, e.g., [46]. One example from Ricarte et al. [41] is shown in Figure 2: a MAD simulation of Sgr A\* with *a*• = 0, *R*high = 40, and a non-thermal "kappa" electron distribution function3 with *κ* = 5 [47,48], inclined at 50◦. The true spectral index map across 214 to 228 GHz is shown in the top central panel, while a one-zone analytic prediction is shown in the top right, using a *κ* = 5 eDF combined with the plasma variables computed in the bottom row. To obtain the analytic prediction, each pixel is treated as a one-zone model using plasma properties computed via an emissivity-weighted average along the geodesic. In this simulation, both *τν* and *B* decline with radius, but Θ*e* stays in a relatively narrow range. Thus, the decline of *α* with radius can be attributed to a decline in *τν* and *B*. Models also exhibit a generic spike in *α* in the photon ring, whose geodesics plunge into regions with strong magnetic fields and acquire a larger optical depth due to their longer path lengths in the emitting region.

**Figure 2.** Intensity and spectral index map of a MAD model of Sgr A\* adapted from Figure 4 of [41]. The top left panel plots total intensity in log scale averaged between 214 and 228 GHz, the top center panel plots the spectral index across this bandwidth calculated by ray tracing the image at two different frequencies, and the top right panel plots an analytic prediction of the spectral index in each pixel obtained by combining the three quantities in the bottom panel: electron temperature, optical depth, and magnetic field strength, each computed by performing an emissivity-weighted average long each geodesic. The excellent agreement between the true spectral index map and the analytic prediction illustrates the power of spectral index maps to jointly constrain these plasma quantities.

#### **3. Linear Polarization**

In flat space, a uniform parcel of synchrotron emitting optically thin plasma intrinsically produces a linear polarization fraction of ≈ 70%, with an orientation perpendicular to its magnetic field projected onto the plane perpendicular to the photon wave-vector, e.g., [49]. The observed orientation of this linear polarization gets modified by two effects: achromatic effects of propagation through the curved space-time, and chromatic "Faraday" effects of propagation through a magnetized plasma (recent theoretical investigations in the vicinity of a black hole include [50–53]). Both rotate and potentially scramble the electric vector position angle (EVPA, or *χ*). Thus, resolved images of linear polarization can inform us about the magnetic field geometry, the magnitude of Faraday effects, and potentially the space-time itself. In Figure 3, we plot the linear polarization pattern of a ring of emission moving clockwise on the sky around a Schwarzschild black hole using the model of Narayan et al. [51]. Here, only the direct (*n* = 0) sub-image is included. The ring is located at a radius of 4.5 *GM*•/*c*2 and has a Gaussian width of 2 *GM*•/*c*2. The angular momentum vector of the ring projected onto the sky has a position angle of 288◦ East of North and is viewed at an inclination of 17◦, consistent with that of M87\* [54]. In the toroidal and radial field cases, the resulting linear polarization pattern is simply perpendicular to the

magnetic fields projected onto the line-of-sight. In the idealized vertical case, the EVPA pattern becomes more sensitive to the direction of the fluid's motion [8].

**Figure 3.** Polarization pattern of a ring of emission around a Schwarzschild black hole threaded with magnetic fields of different geometries: toroidal, radial, and vertical adapted from Figure 3 of [8]. The toroidal and radial magnetic field cases clearly illustrate the fact that synchrotron emission is polarized perpendicular to the magnetic field projected onto the sky. The orientation of the ticks in the vertical field case encodes the direction of the fluid's motion [8], chosen here to be clockwise on the sky. These maps were computed using the analytic ring model of Narayan et al. [51]. Here, the color map encodes the total intensity, and unlike in Figure 1, the linear polarization ticks do not scale with the polarized intensity.

EHT studies have identified the linear polarization fraction (on both resolved and unresolved scales) as well as the morphology of polarization ticks as important observables for theoretical interpretation. This "twistiness" can be quantified by the complex number *β*2, the rotationally symmetric component of a Fourier decomposition of the polarization pattern [55]. M87\* and Sgr A\* both exhibit percent level spatially-unresolved linear polarization fractions at 1.3 mm and ∼10–20% level resolved polarization fractions, e.g., [7,56], much lower than the ideal value of 70% for a uniform parcel of emitting plasma, e.g., [57–59], as do other AGN observed with cm-VLBI, e.g., [60,61]. This could be due to both Faraday depolarization and disordered magnetic field structures. For the EHT observations of M87\*, Event Horizon Telescope Collaboration et al. [8] found that the low polarization fraction was attributable mainly to Faraday rotation in GRMHD models. Combining resolved linear polarization information with an upper limit on circular polarization, Event Horizon Telescope Collaboration et al. [8] found that MAD models were favored over their SANE counterparts for M87\*, which could not be concluded based on total intensity alone. Fundamentally, this can be attributed to linear polarization's sensitivity to the geometry of the magnetic field, as well as Faraday rotation's sensitivity to cooler electrons that may otherwise be invisible. MAD models tend to have more ordered fields with stronger poloidal components, which produces twistier polarization patterns. Meanwhile, SANE models tend to require orders of magnitude larger mass density to compensate for their intrinsically weaker magnetic fields and lower temperatures, resulting in much greater Faraday depths. In retrograde systems, images can exhibit flips in the handedness of the polarization spiral with radius, quantified as a sign flip in the imaginary component of *β*2. This results from a flip in the angular velocity of inflowing streams due to frame dragging [62].

Sgr A\* also exhibits interesting time variability in linear polarization, especially during flares, which are accompanied by large polarization fractions, swings in polarization angle, and "*Q*-*U* loops" on the timescale of hours [63–67]. These can be interpreted as the motion of hotspots or other structures as they light up different parts of the magnetic field structure during their orbit [52,68–72]. The hotspots themselves may originate from "flux eruption events" and magnetic reconnections that occur naturally in MAD accretion flows [73–75]. Thus, time variability of linear polarization offers unique insights into the magnetic field structure and direction of orbital motion that could potentially be linked to the inclination and spin of Sgr A\*. The GRAVITY Collaboration has detected centroid motion coincident with a flare [66,76]. Spatially resolved movies created by the ngEHT would help test the hotspot interpretation, motivating high-cadence monitoring of this source.

#### **4. Rotation Measure**

In an ionic plasma, circularly polarized waves of opposite handedness propagate at different speeds, resulting in a circular birefringence effect known as Faraday rotation. The EVPA of propagating emission rotates an amount sensitive to the density, temperature, and line-of-sight magnetic field. As examined in several studies, internal Faraday rotation is important for depolarizing and scrambling images of GRMHD models of black hole accretion flows [8,42–44]. The magnitude of Faraday rotation has a wavelength-squared dependence, thus it is useful observationally to define the rotation measure RM = *dχ*/*dλ*2, which offers insights into physical parameters of the Faraday rotating plasma. For a linearly polarized emitter entirely behind a uniform Faraday screen, the RM is related to the properties of the screen via

$$\text{RM} = 8.1 \times 10^5 \text{ rad m}^{-2} \int\_{\text{source}}^{\text{observer}} f\_{\text{rel}}(\Theta\_{\varepsilon}) \frac{n\_{\varepsilon}}{1 \text{ cm}^{-3}} \frac{B\_{||}}{G} \frac{ds}{\text{pc}'} \tag{1}$$

where *ne* is the electron number density, *B*|| is the component of the magnetic field parallel to the photon wave-vector, and *f*rel is a correction term suppressing Faraday rotation at relativistic temperatures [77]. For relativistic plasmas, *f*rel(Θ*e*) ≈ log(Θ*e*)/(2Θ2 *e*), while *f*rel asymptotes to 1 as Θ*e* → 0. Here Θ*e* ≡ *kBTe*/*mec*2, *kB* is the Boltzmann constant, *Te* is the electron temperature, *me* is the electron rest mass, and *c* is the speed of light [78]. In GRMHD models, the plasma responsible for synchrotron emission is sometimes completely separate from the plasma responsible for Faraday rotation. For example, some large *R*high SANE models exhibit a cold Faraday rotating midplane sandwiched between emission from their hot jet sheaths [44]. Moreover, while the EVPA probes the magnetic field as projected onto the sky, Faraday rotation is sensitive to the magnetic field along the lineof-sight. Hence, rotation measure and linear polarization can offer a view into electron populations that may otherwise be undetectable from total intensity alone.

At the time of writing, narrow observing bandwidths inhibit our ability to create spatially resolved rotation measure maps with currently published EHT data, but spatially unresolved RM measurements at millimeter wavelengths exist for Sgr A\*, the core of M87\*, and several other LLAGN, e.g., [79–84]. Note also that rotation measure from AGN generally increases with increasing frequency, e.g., [85] and can reach values of the order of 107rad m−2 [86] due to the opacity effect probing regions close to the central engine at the ngEHT frequencies. Without spatial resolution, unresolved rotation measure measurements are difficult to interpret because the assumptions underlying Equation (1) are not believed to generally hold. In GRMHD models, Faraday rotation occurs co-spatially with the plasma, can vary by orders of magnitude in different locations, and can also flip sign across the image due to turbulence [44]. As a result, unresolved EVPA measurements may exhibit significant temporal variation and not strictly follow a *λ*2 law. Figure 4 plots a rotation measure map of a MAD GRMHD model of M87\*, with the spatially unresolved RM written at the bottom of each panel. This turbulence can explain the intra-week time variability of the RM observed for M87\* [59]. On the other hand, Sgr A\* has exhibited a constant sign of RM for decades, suggesting the existence of a more stable (but still variable) foreground Faraday screen [84]. Spatially resolved rotation measure maps could help disentangle the Faraday screen and give insights into both the turbulence of the accretion flow and the magnetic field structure of jets. This may be of increased importance, since EHT observations of Sgr A\* indicate that GRMHD models are too variable [13].

**Figure 4.** Rotation measure map of a MAD simulation of M87\* adapted from Figure 13 of Ricarte et al. [44]. Both positive and negative RM regions are simultaneously present, reflecting flips in the line-of-sight magnetic field direction due to turbulence in the accretion flow. The motion of these structures produces a time variable spatially unresolved RM, written at the bottom of each panel.

#### **5. Circular Polarization**

Circular polarization, Stokes V, can be generated both intrinsically through synchrotron emission or through Faraday conversion, which exchanges linear and circular polarization states, e.g., [35,78]. Circular polarization fractions are much lower than their linear counterparts for both Sgr A\* (*V*/*I* ≈ −1%) and M87\* (|*V*/*I*| 0.8%), making it more challenging to study than linear polarization. In addition, the circular feed basis used for EHT sites makes it more challenging to construct circularly polarized images. However, Stokes V has the potential to encode not only the magnetic field direction and geometry, but also the plasma composition.

Unlike the near unity linear polarization fractions produced by a uniform parcel of plasma, intrinsically emitted Stokes V is typically produced at the ∼1% level for plasma parameters appropriate for M87\* or Sgr A\*. The sign of intrinsically emitted circular polarization generated encodes the sign of the magnetic field along the line of sight, following the right hand rule. Stokes V can also be generated via Faraday conversion, which exchanges linear and circular polarization states. The sign of circular polarization generated by conversion depends on the relative orientation of the EVPA with respect to the local magnetic field. Specifically, Stokes V generated by conversion inherits the sign of Stokes U, defined in the local plasma frame. Thus, Stokes V from Faraday conversion is sensitive to the line of sight "twist" in the magnetic field as well as any Faraday rotation affecting the EVPA of the linear polarization that gets converted. Interestingly, Faraday conversion has no effect in a unidirectional magnetic field lacking Faraday rotation.

Both intrinsic synchrotron emission and Faraday conversion are believed to be important for generating circular polarization on event horizon scales [36–38]. For Faraday conversion, both the line of sight twist in the magnetic field and Faraday rotation are important for determining the orientation of the linear polarization that is converted into circular. Small amounts of Faraday rotation can affect the relative alignment between linear polarization and the converting magnetic field. Large amounts of Faraday rotation can lead to depolarization by randomizing the sign of Stokes U that is converted into Stokes V.

Figure 5 depicts a cartoon of a typical helical field geometry as well as the circular polarization produced by a model that exhibits this structure reproduced from [37]. In the left panel, the background colors depict the line of sight direction of the magnetic field, viewed edge-on. In the time averaged image of *V*/*I* viewed at 90◦, this structure produces a "four quadrants" pattern. Outside the immediate vicinity of the event horizon, Stokes V originates from intrinsic emission in this model, which reflects the line-of-sight direction of the magnetic field see also [38]. Another interesting feature that arises due to a generic helical field geometry is the successive sign flipping of sub-images in the photon ring also discussed in [36], which can be explained by Faraday conversion and parallel transport in a generic helical field geometry viewed face-on [37]. In this particular model, the spatially unresolved Stokes V is surprisingly dominated by this sign-flipped photon ring.

**Figure 5.** Maps of circular polarization encode properties of the geometry of the magnetic field, both its line-of-sight direction and twist. A cartoon of a generic helical field geometry is depicted on the left. On the right, we plot the circular polarization of a MAD model of M87\* at two inclinations. Both are reproduced from Ricarte et al. [37]. The top row depicts a 5◦ viewing angle, and the bottom row depicts a 90◦ viewing angle. The first column shows the time averaged circularly polarized image, the second column shows the same at a single snapshot, and the third column shows fractional circular polarization. For face-on viewing angles, the photon ring exhibits an interesting sign flip due to Faraday conversion and the sourcing of photons from the opposite side of the disk. For edge-on viewing angles, circular polarization exhibits a "four quadrants" pattern that reflects the line-of-sight magnetic field direction.

Finally, circular polarization is strongly affected by plasma composition, and can potentially be used to distinguish pair plasmas from ionic plasmas, e.g., [39]. Most models used to study EHT images have contained only ionic plasma. Electron-positron pairs can also be naturally produced on event horizon scales, but their abundance is highly theoretically uncertain. In a pair plasma, intrinsic circularly polarized emission and Faraday rotation both vanish, but Faraday conversion persists. Intuitively, a pair plasma with equal parts positively and negatively charged particles should not gyrate in a preferred direction. This can cause dramatic differences in images of circular polarization, and potentially also those of linear polarization [29,30].

#### **6. Scattering**

A major challenge for studies of Sgr A\* with the EHT and ngEHT is interstellar scattering by dilute plasma in the ionized interstellar medium. This scattering arises from variations in the index of refraction of the plasma because of fluctuations in the electron number density see, e.g., [87]. In particular, the line of sight to Sgr A\* is heavily scattered by plasma in the spiral arms of the Milky Way [88], resulting in angular broadening that is approximately three orders of magnitude larger than median values for lines of sight at higher galactic latitudes. The effects of scattering are two-fold: (1) small-scale modes in the scattering material result in diffractive "blurring," described by a convolution with an anisotropic kernel, and (2) large-scale modes in the scattering material result in refractive "substructure," described by additive image noise with a slowly falling power spectrum. For detailed discussion of the scattering of Sgr A\*, see Psaltis et al. [89], Johnson et al. [90,91].

While the scattering severely affects images of Sgr A\*, many polarimetric properties of the images are comparatively immune because the scattering is not significantly birefringent. For example, in the case of purely diffractive scattering, the image-integrated fractional polarization is independent of scattering. More generally, the interferometric fractional polarization *m*˘(**u**) ≡ *P*˜(**u**)/ ˜*I*(**u**) is independent of diffractive scattering because convolution is multiplicative in the visibility domain and is identical for all Stokes parameters, thereby canceling in the quotient. Here, **u** is the dimensionless baseline vector, projected orthogonal to the line of sight and measured in units of the observing wavelength, which corresponds to the sampled image wavevector in Fourier space; *F*˜ denotes the Fourier transform of some quantity *F*.

Figure 6 shows an example GRMHD snapshot before and after scattering, while the images look substantially different, key polarimetric observables such as the *β*2 mode, which is highly constraining for GRMHD models [8] and carries information about black hole spin [55], are almost unaffected by scattering. Likewise, certain interferometric observables, such as the interferometric fractional polarization, are only mildly affected by scattering (see Figure 7). For additional discussion of how the deterministic frequency dependence of scattering can be used for scattering mitigation on images, see [56]; for discussion of how the lack of birefringence can be used to study the relative power spectra in different polarization modes, see Ni et al. (in prep).

**Figure 6.** GRMHD model of Sgr A\* at 230 GHz before (**left**) and after (**right**) including the effects of interstellar scattering. This simulation is a MAD with *a*• = 0.7, *R*high = 20, and *i* = 30◦. The background image shows total intensity with respect to the image peak, while the ticks show the polarization magnitude and direction, colored by fractional polarization, while scattering severely affects the image, key polarimetric measures are nearly immune to scattering. For example, the unresolved fractional polarization is 10.5% before scattering and is 10.6% after scattering. Likewise, the *β*2 mode in polarization [55] has |*β*2| = 0.40 and arg(*β*2) = 52.1◦ before scattering, and |*β*2| = 0.37 and arg(*β*2) = 51.0◦ after scattering.

**Figure 7.** Interferometric properties of the GRMHD model shown in Figure 6. Solid lines show the normalized intensity ˜*I*(**u**)/ ˜*I*(**0**) before (blue) and after (red) scattering, and the dashed lines show the interferometric fractional polarization magnitude |*m*˘(**u**)|. For these curves, baselines are oriented along the East–West direction: **u** = (*u*, 0). Over the full range of baseline lengths accessible from the ground, the fractional polarization is largely immune to scattering, while diffractive scattering causes a substantial reduction in the flux on long baselines.

#### **7.  In practice, since the measured visibilities are actually samples of the Fourier transform, image reconstruction can introduce significant systematic uncertainties. Reconstruction methods must find images consistent with incomplete and noisy information in Fourier space, to which there can be multiple families of solutions. Images cannot be constructed at all without sufficient *uv* coverage (or strict image priors). Thus, it can be useful to study signals in their native visibility space.

The visibility-domain response of polarimetric observables has some key differences from that of total intensity. For instance, the visibility amplitude for total intensity is guaranteed to be maximal for the zero baseline because the image is positive. However, because the Stokes parameters *Q*, *U*, and *V* are not constrained to be positive, their visibility functions may not be maximal on the zero baseline. This simple property can be used to make powerful inferences from sparse measurements (e.g., from the EHT or the space VLBI project *RadioAstron*). For instance, a single measurement of |*V*˜(**u**)| > |*V*˜(**0**)| would demonstrate that the image does not have uniform *sign* of circular polarization. Likewise, if *P*˜(**u**) = *P*˜(−**u**) , where *P*˜(**u**) ≡ *Q*˜(**u**) + *iU*˜ (**u**), then the linear polarization field *must* have variations in direction. This test can be performed with as few as two stations, but each must have dual polarization receivers. In addition, because polarization breaks the trivial baseline reflection degeneracy that is present in total intensity, the polarimetric structure can be compared on baselines **u** and −**u**, effectively doubling the angular resolution relative to total intensity alone!

Another significant difference between the total intensity and polarization is that the linear polarization images can have changes in both amplitude and *direction*, allowing it to have significant image substructure relative to the intensity image. In the visibility domain, this substructure translates to a relative increase in the power on long baselines in polarization versus total intensity. In the limit of a heavily resolved source, the polarimetric signal may exceed that of the total intensity, even for a source with a low fractional polarization!

This can be quantified using the interferometric fractional polarization, *m*˘(**u**) ≡ # *Q*˜(**u**) + *iU*˜ (**u**) \$ / ˜*I*(**u**). On a zero baseline, this complex quantity corresponds to the unresolved fractional polarization, *m*˘(**0**)=(*Q*tot + *iU*tot)/*I*tot, where the "tot" subscript denotes an image-integrated quantity. However, unlike the image fractional polarization, |*m*˘(**u**)| can exceed unity on long baselines. This was found out from ground-based observations of for Sgr A\* [92] as well as within the *RadioAstron* Space VLBI survey of AGN [93]. These observations had one common feature: a very high angular resolution corresponding to several G*λ* spatial scales. In general, for a heavily resolved source with polarized substructure, we expect |*m*(**u**)| to generically grow with increasing baseline length.

From a calibration perspective, the interferometric fractional polarization has the benefit of properties analogous to VLBI closure quantities, since the source of rapidly varying gains at VLBI sites (e.g., changing atmospheric delay and reference frequency errors) are equivalent for both polarization feeds. In addition, the interferometric fractional polarization is resilient to the effects of interstellar scattering, which is likewise not significantly birefringent (see Figure 7). Finally, |*m*(**u**)| is a useful observable to measure the relative coherence of the polarization field. For a perfectly uniform polarization field, |*m*(**u**)| will be independent of baseline length. However, for a polarization field that varies on spatial scales much smaller than the the beam size of the interferometer, |*m*(**u**)| will grow roughly as 1/|˜*I*(**u**)|, as the observations resolve the structure in total intensity without resolving the structure in polarization.

In Figure 8, we plot total intensity and linear polarization maps of two models of M87\* with very different polarization characteristics. The top model is a MAD simulation with *a*• = 0.9 and *R*high = 1, which exhibits an ordered polarization pattern due to ordered magnetic fields and little Faraday rotation. Meanwhile, the bottom model is a SANE simulation with *a*• = −0.3 and *R*high = 40. As mentioned in Section 4, SANE simulations tend to have much larger Faraday depths than MADs, causing this model's linear polarization to be much more disordered. The pixel-to-pixel intensity weighted Faraday depth *τF* is written at the bottom of each panel.

**Figure 8.** Two GRMHD models imaged at 228 GHz and corresponding maps of linear polarization in visibility space. The top row corresponds to a MAD model of M87\* with *a*• = 0.9 and *R*high = 1, while the bottom row corresponds to a SANE model with *a*• = −0.3 and *R*high = 40. Due to a much larger Faraday depth, written at the bottom of the images, the SANE model exhibits a much more disordered linear polarization pattern. In the ordered model, measures of the linear polarization rise dramatically with radius in the Fourier domain, while the disordered model is characterized by blobs with a coherence length corresponding to the size of the image.

These characteristics are reflected in their Fourier space maps of |*P*˜| and |*m*˘ |, shown in the second and third columns, which can be directly sampled using an interferometer. Here, *P*˜ ≡ *Q*˜ + *iU*˜ and *m*˘ ≡ *P*˜/ ˜*I*, where ∼ denotes a Fourier transform. For the MAD model, *P*˜ and *m*˘ both rise dramatically with radius in Fourier space in the central region. *P*˜ rises because the linear polarization is higher on resolved scales than a spatially unresolved measurement would suggest. The rotational symmetry of the polarization pattern causes substantial cancellation of polarization without spatial resolution. *m*˘ also rises for this reason, and also because ˜*I* exhibits nulls in Fourier space that do not necessarily coincide with the nulls in *P*˜. Meanwhile, the disordered SANE simulation exhibits a mottled pattern in *P*˜ with a characteristic length scale corresponding to the size of the image. These phenomena should not change much qualitatively as a function of wavelength in the sub-millimeter.

The previous discussion has focused on the relationship between expected image features and their appearance in the (Fourier-conjugate) interferometric visibility domain. However, a crucial consideration is how to study frequency-dependent effects, such as spectral index and rotation measure, using interferometry. Because the fringe spacing *u* ∝ *ν*, interferometric measurements across multiple frequencies necessarily mix the effects of a changing dimensionless baseline with those of a changing image.

Specifically, the polarized visibility *P*˜ on a physical vector baseline **b** at an observing frequency *ν* is

$$\tilde{P}(\mathbf{b},\boldsymbol{\nu}) = \int d^2\theta \, P(\theta,\boldsymbol{\nu}) e^{-2\pi i \frac{\boldsymbol{\nu} \cdot \mathbf{b}}{\boldsymbol{\ell}} \cdot \mathbf{b} \cdot \mathbf{b}} \tag{2}$$

$$\Rightarrow \nu \partial\_{\boldsymbol{\nu}} \mathcal{P}(\mathbf{b}, \boldsymbol{\nu}) = \int d^2 \theta \left[ \nu \partial\_{\boldsymbol{\nu}} P(\theta, \boldsymbol{\nu}) - 2\pi i \theta \cdot \mathbf{u} P(\theta, \boldsymbol{\nu}) \right] e^{-2\pi i \theta \cdot \mathbf{u}},\tag{3}$$

where *∂ν* denotes a partial derivative with respect to *ν*. The first term in the square brackets of Equation (3) accounts for the frequency dependence of the image, while the second accounts for the changing dimensionless baseline with frequency. For the first term, an image with spectral index *α* has *P* ∝ *να*, so |*ν∂νP*|∼|*αP*|; likewise, an image with finite rotation measure has *P* ∝ *e*2*iλ*2RM, so |*ν∂νP*| ∼ 4 RM × *λ*2*P* . Roughly speaking, we expect that the relative dominance or subdominance of spectral index versus RM are independent of baseline length, so the relative effects on long baselines are likely similar to those for unresolved measurements of a source. For instance, the effects of rotation measure for observations of Sgr A\* at millimeter wavelengths are likely to vastly dominate the effects of spectral index. Sgr A\* has *α* = 0.0 ± 0.1 [94] but 4 × RM × *λ*2 ≈ −2.7 [59]. The second term gives a relative contribution that increases as the image is increasingly resolved. It becomes significant when the spanned frequencies change the baseline length by the inverse field-of-view, *F*.

#### **8. Discussion and Conclusions**

The EHT and upcoming ngEHT enable us to probe accreting supermassive black holes via a variety of multi-frequency polarimetric observables. In this contribution, we have discussed the many ways in which the physical properties of underlying accretion flow are mapped onto these observables. Total intensity and spectral index encode the density, temperature, and magnetic field strength of emitting plasma in different regions. Linear polarization encodes the geometry of the magnetic field, and its depolarization via Faraday rotation offers an observational probe into what may be otherwise invisible cool electrons, and a different component of the magnetic field vectors. Rotation measure maps probe this cooler Faraday rotating electron population directly, and can probe the magnetic field direction, which can reveal turbulent structures. Finally, circular polarization encodes both overall geometry and direction of the magnetic field via emission, Faraday rotation, and Faraday conversion. We have discussed that even if imaging proves prohibitively challenging for some datasets, constraining information exists already in visibility space. For some models, low polarization fractions in spatially unresolved measurements hide large polarization fractions in spatially resolved measurements. In the ngEHT era, we will have access to not only single snapshots, but also movies, with a much larger dynamic range in intensity than is presently possible with the EHT. Multiple snapshots will also enable cleaner theoretical connections via time averaging, e.g., [95,96]. This will enable an unprecedented deluge of data about the nearest SMBHs that will help us understand their accretion and feedback processes.

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

**Funding:** We thank the National Science Foundation (AST-1716536, AST-1935980, AST-2034306, AST-1816420, and OISE-1743747) for financial support of this work. YYK was supported by the Russian Science Foundation grant 21-12-00241. RE acknowledges the support from NASA via grant HST-GO-16173.001-A. This work was supported in part by the Black Hole Initiative, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University. The opinions expressed in this publication are those of the author(s) and do not necessarily reflect the views of the Moore or Templeton Foundations.

**Data Availability Statement:** The data underlying the figures of this article can be provided upon reasonable request to the corresponding author.

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

#### **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.

## *Article* **Modeling Reconstructed Images of Jets Launched by SANE Super-Eddington Accretion Flows around SMBHs with the ngEHT**

**Brandon Curd 1,2,\*, Razieh Emami 2, Freek Roelofs 1,2 and Richard Anantua 3**

**Abstract:** Tidal disruption events (TDEs) around supermassive black holes (SMBHs) are a potential laboratory to study super-Eddington accretion disks and sometimes result in powerful jets or outflows which may shine in the radio and sub-millimeter bands. In this work, we modeled the thermal synchrotron emission of jets by general relativistic radiation magneto-hydrodynamics (GRRMHD) simulations of a BH accretion disk/jet system which assumed the TDE resulted in a magnetized accretion disk around a BH accreting at ∼12–25 times the Eddington accretion rate. Through synthetic observations with the Next Generation Event Horizon Telescope (ngEHT) and an image reconstruction analysis, we demonstrate that TDE jets may provide compelling targets within the context of the models explored in this work. In particular, we found that jets launched by a SANE super-Eddington disk around a spin *a*∗ = 0.9 reach the ngEHT detection threshold at large distances (up to 100 Mpc in this work). A two-temperature plasma in the jet or weaker jets, such as a spin *a*∗ = 0 model, requires a much closer distance, as we demonstrate detection at 10 Mpc for limiting cases of *a*∗ = 0, R = 1 or *a*∗ = 0.9, R = 20. We also demonstrate that TDE jets may appear as superluminal sources if the BH is rapidly rotating and the jet is viewed nearly face on.

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

### **1. Introduction**

Tidal disruptions of stars by supermassive black holes (SMBHs), or tidal disruption events (TDEs), have recently become a regularly observed transient phenomenon. Stars which enter the tidal radius

$$R\_t \simeq 7 \times 10^{12} \left(\frac{M\_{\rm BH}}{10^6 M\_\odot}\right)^{1/3} \left(\frac{M\_\*}{M\_\odot}\right)^{-1/3} \left(\frac{R\_\*}{R\_\odot}\right) [\text{cm}] \tag{1}$$

of the central SMBH in their host galaxy will be disrupted [1,2], either partially or fully, depending on the orbit and equation of state of the star [3,4]. The bound stream of gas returns towards the BH delivering mass at the fallback rate (*M*˙ fb). Apsidal precession of the returning stream leads to self-intersection with material that has yet to pass through pericenter and leads to dissipation and disk formation. Dissipation due to self-intersection may also be a source of early emission in a TDE [5]. After the initial rise to peak, the fallback rate follows a power law behavior which can be approximated as

$$
\dot{M}\_{\rm fb} = \dot{M}\_{\rm fb,peak} \left(\frac{t}{t\_{\rm fb}}\right)^{-5/3}.\tag{2}
$$

**Citation:** Curd, B.; Emami, R.; Roelofs, F.; Anantua, R. Modeling Reconstructed Images of Jets Launched by SANE Super-Eddington Accretion Flows around SMBHs with the ngEHT. *Galaxies* **2022**, *10*, 117. https://doi.org/10.3390/ galaxies10060117

Academic Editor: Bidzina Kapanadze

Received: 14 November 2022 Accepted: 7 December 2022 Published: 13 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/).

Here, *M*˙ fb,peak is the peak mass fallback rate:

$$\frac{\dot{M}\_{\text{fb,peak}}}{\dot{M}\_{\text{Edd}}} \approx 133 \left( \frac{M\_{\text{BH}}}{10^6 M\_{\odot}} \right)^{-3/2} \left( \frac{M\_{\ast}}{M\_{\odot}} \right)^2 \left( \frac{R\_{\ast}}{R\_{\odot}} \right)^{-3/2} \tag{3}$$

making the "frozen in" approximation as in Stone et al. [6], and

$$t\_{\rm fb} = 3.5 \times 10^6 \,\mathrm{s} \left(\frac{M\_{\rm BH}}{10^6 M\_{\odot}}\right)^{1/2} \left(\frac{M\_\*}{M\_{\odot}}\right)^{-1} \left(\frac{R\_\*}{R\_{\odot}}\right)^{3/2},\tag{4}$$

is the fallback time, which is the orbital time of the most bound part of the stream. Of note is the fact that the mass fallback rate can greatly exceed the Eddington mass accretion rate *M*˙ Edd. Since the mass accretion rate is expected to be similar in magnitude to the mass fallback rate, it is possible that TDEs result in a super-Eddington accretion disk, which are geometrically thick, radiatively inefficient accretion disks [7,8]. The exact power law behavior varies with the properties and orbit of the star [9].

TDEs are typically seen as optical/X-ray transients [10,11], but several TDEs have resulted in outflows or jets which shine in the radio bands [12]. In the most common case in which no relativistic jet is launched (commonly referred to as "non-jetted" TDEs), the X-ray and optical/UV luminosity follows a roughly *t* −5/3 decline, similar to the fallback rate. If the TDE leads to prompt disk formation, the X-rays are thought to arise from an accretion disk, and the optical/UV emission arises from a large scale reprocessing layer [13,14]. TDEs have also been observed to launch relativistic X-ray jets in a few cases. These jetted TDEs have been argued to arise due to a magnetically arrested disk (MAD) [15] forming around the BH during the TDE [13,16,17], which leads to jet production via the Blandford–Znajek (BZ) mechanism [18] extracting spin energy from the BH. Alternatively, these powerful relativistic jets may be produced thanks to radiative acceleration of gas through a narrow funnel region [19,20].

It is important to note that the BZ mechanism is the strong field limit of the magnetic Penrose process (MPP), but the MPP can operate with weaker magnetic field strengths [21,22]. Note as well that the magnetic Penrose process, as its name suggests, is itself the magnetic flavor of the Penrose process (PP) [23]. While the PP relies on the existence of negative energy orbits in the vicinity of the BH, the MPP relies on negative energy orbits and quadrupole electric fields produced by twisting of the magnetic field lines threading the BH horizon. A key distinction, and why the BZ mechanism is often cited while interpreting astrophysical jets, is that the energy extraction from the MPP can exceed 100 percent, whereas the maximum energy extraction from a PP is ∼20% [21,22].

A handful of non-jetted TDEs have been observed to produce radiowaves peaking at tens of GHz with *L*radio ∼ 1037−39erg s−1. This emission is thought to arise from an outflow launched by the TDE with velocity *v* ∼ 0.1*c* shocking on the gas surrounding the BH. Meanwhile, jetted TDEs produce bright radio emission peaking at *L*radio ∼ 1040−42erg s−1. The appearance of radio emission is often delayed by several weeks from the initial appearance of the optical/UV/X-ray emission in non-jetted TDEs, which hints at some connection to the disk formation process to the occurrence of outflows.

When *M*˙ fb rises to peak, it has previously been assumed that by this stage a circularized accretion disk has formed [13,17]. The first direct demonstration of circularization near the peak fallback rate was recently demonstrated in a numerical simulation by Steinberg and Stone [5]. Multiple authors have argued in favor of a picture in which an inner accretion flow is surrounded by a quasi-spherical reprocessing layer, since this naturally explains the sometimes delayed appearance of X-ray emission in optical/UV-discovered TDEs [13,14,24]. This picture naturally arises if the accretion flow is actually super-Eddington [13], which has motivated multiple studies of GRRMHD simulations magnetized of super-Eddington disks. Motivated by this fact and the demonstration that a "standard and normal evolution" (SANE) [15] super-Eddington accretion disk also leads to viewing angle effects which may explain the behavior in non-jetted TDEs [17], we conducted a study of the outflows

launched by SANE models in Curd et al. [14] and studied their radio-submm emission. We focus on these SANE models in this work as well. The BZ mechanism also acts in the SANE model [14,25]; however, the disk and BH have a much lower magnetic flux and thus a substantially lower jet power. Although BZ driven outflows in SANE disks can be relativistic, they are weakly comptonized and do not appear as powerful, beamed X-ray sources like the presently known jetted TDEs [14].

The Next Generation Event Horizon Telescope (ngEHT) will provide more baseline coverage and faster response times than the previous mission. An estimate for ngEHT is to double the antenna sites [26] of its 20 μas predecessor EHT, and thus the number of possible baseline pairs and triads of sites available for imaging jet/accretion flow/black hole systems will scale combinatorially (the number of baselines grows with the number of antennae as (N(N − 1)/2). This could allow for interesting sources, such as jets from nearby tidal disruption events, to be imaged directly. In our previous work [14], we provided the first demonstration that SANE super-Eddington accretion flows can produce radio emission which is bright enough at 230 GHz to be detected and resolved. Here we take things a step further and produce reconstructed images assuming such jets happen in the nearby universe.

The detection rate of TDEs in the optical/UV/X-ray is expected to grow rapidly once the Large Synoptic Survey Telescope comes online [27,28]. Assuming rapid follow-up of TDEs in radio-submm bands finds detectable emission, this could provide a large number of targets for the ngEHT. As we demonstrated in Curd et al. [14], some models produced detectable emission, even at ∼180 Mpc. At this distance, a conservative estimate of the volume integrated TDE rate suggests more than 200 TDEs per year assuming volumetric TDE rates based on Stone and Metzger [29]. Even at < 40 Mpc, we estimate that several TDEs should occur per year (see Figure 2 in Curd et al. [14]), which suggests some nearby TDEs may become targets of opportunity during the ngEHT mission.

We stress that jets such as those in our first work on the subject of jets from SANE models of TDE accretion disks [14] do not resemble any previously detected radio TDEs. This may suggest most, or even all, TDEs do not form accretion disks which resemble SANE models to begin with. However, the number of TDEs that have appeared in the radio-submm in the first place is extremely small as of this writing; fewer than twenty radio TDEs have been reported. Furthermore, magnetic fields were certainly present in the forming disk, albeit dynamically subdominant to hydrodynamic effects early in the disk formation [30,31]. Nevertheless, it is possible that after the disk circularizes, which Steinberg and Stone [5] suggests may take tens of days, the magnetic field builds up in a dynamo effect similar to Sadowski et al. [30]. In this case, one would almost certainly expect the magnetic field to become dynamically important, in which case a magnetized outflow may be launched as in Curd et al. [14]. In addition, if the currently known jetted TDEs are indeed MAD disks around rapidly rotating BHs, it stands to reason that a sub-class of less strongly magnetized disks which could be described by a SANE model should exist. TDEs continue to surprise observers in terms of the range of behavior, so such a jet formation channel may yet be discovered.

#### **2. Numerical Methods**

#### *2.1. GRRHMD Simulations*

Throughout this work, we often use gravitational units to define length and time. In particular, we use the gravitational radius

$$r\_{\mathcal{S}} = \frac{GM\_{\rm BH}}{c^2} \tag{5}$$

and the gravitational time

$$t\_{\rm g} = \frac{GM\_{\rm BH}}{c^3} \,, \tag{6}$$

where *M*BH is the mass of the black hole (BH). We also adopt the following definition for the Eddington mass accretion rate:

$$
\dot{M}\_{\text{Edd}} = \frac{L\_{\text{Edd}}}{\eta\_{\text{NT}}c^2} \tag{7}
$$

where *L*Edd = 1.25 × 1038 (*M*BH/*M* ) erg s−1 is the Eddington luminosity and *η*NT is the radiative efficiency of a Novikov–Thorne thin disk around a BH with spin parameter *a*∗ [32].

We conducted an imaging analysis of GRRMHD simulations presented in Curd et al. [14]. In particular, we analyzed the most massive BH models m7a0.0-HR and m7a0.9-HR, which are *M*BH = 107 *M* BHs of spin *a*∗ = 0 and 0.9 BHs. We specify the simulation diagnostics relevant for this work in Table 1. The simulations were conducted in 2D (*r*, *ϑ*) coordinates on a *Nr* × *Nϑ* = 640 × 256 grid with added resolution near the poles to adequately resolve both the disk and jet. The radial grid cells were logarithmically spaced with a maximum domain radius of *R*max = 105 *rg* to capture the large-scale features of the jet.

**Table 1.** We tabulate the mass accretion rate *M*˙ , jet efficiency *η*jet ≡ *L*jet/*Mc* ˙ 2 (where *L*jet is the jet power as defined in [14]), and total simulation duration *t*sim for each KORAL simulation. Note that *M*˙ and *η*jet are time averaged over the final 50,000 *tg* of each simulation.

On horizon scales, gas is flowing across the BH horizon in an accretion disk due to angular momentum transport driven by the magneto-rotational instability. The disk is optically thick and turbulent, and gas inside of the disk is advected with the gas across the BH horizon. However, an optically thin funnel above and below the disks exists. Here, radiation can escape freely and pushes on gas, accelerating a significant outflow. In addition, the funnel is magnetized and sometimes exhibits magnetization parameter *σ* = *b*2/*ρc*2 > 1, where *b* is the magnetic field strength and *ρ* is the mass density of the plasma. In the jet, magnetic energy is partially converted into kinetic energy, as it contributes to accelerating gas into an outflow.

In both simulations, radiative and Poynting acceleration drove fast outflows. The jet reached relativistic speeds with Lorentz factor *γ* > 5 for model m7a0.9-HR, which was likely due to the BZ effect extracting spin energy from the BH, which may produce roughly ∼1 percent of the jet efficiency, even though the magnetic flux threading the black hole was well below the MAD limit [14,25]. The primary sites of dissipation were the jet head and internal shocks inside of the jet. Internal shocks were due to fast- and slow-moving gas interacting downstream of the jet head, in addition to recollimation shocks. As we show in Figure 1, this resulted in a hot, magnetized jet which reached large scales (*r* > 30,000 *rg*) by the end of the simulation. The *a*∗ = 0.9 model had a significantly more magnetized jet and also produced a more powerful jet by a significant fraction (see Table 1).

The simulations were conducted in 2D to overcome the substantial numerical requirements of resolving the accretion disk and jet in a global GRRMHD simulation. While SANE super-Eddington disks can be described by an axisymmetric flow quite well [33], the jet may undergo 3D instabilities which cannot be captured in our numerical setup. In particular, the jet may become kink-unstable as it propagates [34]. This will certainly change the morphological features, and additional dissipation of magnetic energy along the jet may occur. If the dissipation is extreme, especially at large scales such as those we consider here, the synchrotron emission which we describe in the next section may differ from our analysis of 2D jets.

**Figure 1.** Here we show snapshots of the GRRMHD KORAL simulations that we post-processed with ipole. All data are shown for *t* = 78, 000 *tg* in both m7a0.0-HR (**top**) and m7a0.9-HR (**bottom**). The colors indicate the gas density *ρ* (**left**), gas temperature *T*gas (**middle**), and magnetic field strength |*B*| (**right**); the yellow contours indicate the *σ* = 1 boundary, in which we set *ρ* = 0 in the ray tracing step to prevent emission.

#### *2.2. 230 GHz Emission*

We post-processed the KORAL simulation data with the general relativistic ray tracing (GRRT) code ipole [35–37], which includes synchrotron and Bremsstrahlung emission and absorption. Since the KORAL simulation was conducted in 2D, we first copied the 2D data onto a full 3D (*r*, *ϑ*, *φ*) grid with 32 cells in azimuth which spanned *φ* = 0 − 2*π* before ray tracing. We did this by simply copying the original 2D data at each *φ* cell such that the new 3D data maintained the assumed axisymmetry in *φ*. The electron distribution function was assumed to be thermal. Ohmura et al. [38,39] demonstrated that large-scale active galactic nuclei (AGN) jets can produce a two-temperature plasma. Motivated by their findings and the possibility that a two-temperature plasma will be produced due to shocks at the jet head and within the jet itself, we tested a simple two-temperature jet model by scaling the electron temperature relative to the ion temperature via the plasma temperature ratio:

$$\mathcal{R} = \frac{T\_i}{T\_\varepsilon},\tag{8}$$

where *Ti* is the temperature of the ions and *Te* is the temperature of the electrons. Note that *Ti* was obtained directly from the KORAL simulation by setting *Ti* = *T*gas.

The peak of the radio-submm spectra in m7a0.0-HR is lower than that of m7a0.9-HR, so increasing R has a much more significant impact on the 230 GHz emission and can make the jet undetectable even at 10 Mpc for values of R > 2 [14]. It is possible that a non-thermal electron distribution will have greater high-energy emission even as R increases, but we saved an exploration of non-thermal electron models for a future analysis.

Each model was imaged at 230 GHz. For both models, we imaged the simulation at times *t* = 38,000 *tg* and *t* = 78,000 *tg* for a difference in observing times of ∼23 days. We chose a distance *D* = 10 Mpc, R = 1, observing angles relative to the jet axis (*z* in Figure 1) of *θ* = 10, 45, and 90◦, respectively. Note that we use *θ* for the observer angle, and *ϑ* is the polar angle in the KORAL grid coordinates. For model m7a0.9-HR, we also tested limiting cases *D* = 100 Mpc, R = 1, and *D* = 10 Mpc, R = 20, imaged at *θ* = 90◦. The total 230 GHz flux of each ray-traced model is tabulated in Table 2.

We show a full library of each of the ipole images convolved with a Gaussian beam with a 20 μas full width at half maximum (FWHM) in Figures A1 and A3.

**Table 2.** Here we tabulate the 230 GHz flux density for each model given a specific time, viewing angle *θ*, distance *D*, and temperature ratio R.

#### *2.3. Synthetic ngEHT Observations and Image Reconstruction*

In order to test to what extent the jet features in our models can be observed, we simulated observations with a potential ngEHT array, consisting of the 2022 EHT stations and 11 additional stations, selected from Raymond et al. [40] and similar to the ngEHT reference array used in the ngEHT Analysis Challenges [41]. The new dishes were assumed to have a diameter of 10 m and a receiver temperature of 50 K, with the array operating at a bandwidth of 8 GHz. For each image, we simulated a 24-hour observation with a 50% duty cycle with this array, using the ngehtsim 1 library, which makes use of eht-imaging [42,43]; see also [44]. The atmospheric opacity was set to reflect a good day in April, using the top 1*σ* quantile from the MERRA-2 data interpolated and integrated for each site ona3h cadence for a 10-year period [45,46]. Thermal noise was added to

the complex visibilities, visibility phases were randomized, and no systematic visibility amplitude errors were added to the data.

We subsequently used the regularized maximum likelihood framework in eht-imaging to produce image reconstructions, with maximum entropy and (squared) total variation regularizers, fitting to visibility amplitudes and closure phases; see [42,43,47]. After establishing a set of well-performing imaging parameters on the m7a0.9-HR, R = 1 model at *tg* = 78,000 and a distance of 10 Mpc (Figure 2), we applied the same script to all other simulated datasets.

#### **3. Results**

In this section, we comment on the detectability of our models and then compare the ipole images with the reconstructed images. We comment on features which may be of interest in terms of the broader study of astrophysical jets.

#### *3.1. Reconstructed Images*

In the ray-traced image (i.e., see the left panel in Figure 2), the jet head produces bright emission as it shocks on the circumnuclear medium (CNM). In addition, various shocks occur within the jet due to both slow- and fast-moving components colliding radially and due to recollimation shocks. This leads to dissipation within the jet and bright "bubbles" of emission at 230 GHz. As we show in the right panel of Figure 2, the jet head and the structures in the jet are faithfully reproduced in the reconstruction for favorable viewing angles (*θ* = 45◦ and 90◦) as long as the source is nearby (*D* = 10 Mpc). Jets viewed near *θ* = 10◦ are dominated by emission from the jet head, and distinguishing internal jet features would be unlikely. This can be seen by comparing the base images with the full library of reconstructed images for models m7a0.0-HR and m7a0.9-HR (Figures A1–A4).

For distant sources (*D* = 100 Mpc), distinguishing internal features is impossible, and only the jet head can be fully distinguished in the reconstruction (right panel in Figure 3). This would still allow for the jet motion to be tracked, but detailed information is lost.

**Figure 2.** Model m7a0.9-HR at *t* = 78,000 *tg* imaged at *D* = 10 Mpc with *θ* = 90◦ and R = 1. We show the base ipole image with no blurring (**left**), the base ipole image blurred via convolution with a 20 μas FWHM Gaussian beam (**middle**), and the reconstructed image blurred using the same Gaussian beam (**right**).

**Figure 3.** The same as Figure 2 but for model m7a0.9-HR at *t* = 78,000 *tg* imaged at *D* = 100 Mpc with *θ* = 90◦ and R = 1.

#### *3.2. Tracking Jet Motion*

In this subsection, we demonstrate that the original ray-traced images and the reconstructed images allow for the jet motion to be tracked and yield similar results for the time evolution of the jet. The jet features are approximately Lorentzian, so we fit Lorentzian profiles to the image to find the position of the top and bottom jet in both the base images and the reconstructed images. We detail the peak finding algorithm in Appendix B. Since we cannot properly center the jet in the reconstructed images (there is no bright, central radiation from the near BH), we only measured the distance between the two jet peaks *y*1 and *y*2, respectively. We define the apparent jet length as:

$$l\_{\text{jet}} = |y\_2 - y\_1|. \tag{9}$$

Note that we have not differentiated the "top" or "bottom" jet here, as we are only concerned with the total distance between the jet heads. We obtained errors on the jet length from the error estimates of the jet head locations using standard error propagation analysis:

$$
\delta l\_{\rm jet} = \sqrt{\delta y\_1^2 + \delta y\_2^2}. \tag{10}
$$

We computed the relative difference between the jet lengths for the ray-traced (*l*jet) and reconstructed (*l*jet,rec) images in order to quantify the extent to which measurements of the jet length and velocity agree. We found that, generally, the ray-traced and reconstructed images yield similar jet lengths within 3 standard deviations (Figure 4). In general, there are much larger errors on the fit for the Lorentzian profile's center at steep angles (i.e., see the relative difference for *θ* = 10◦). The agreement is also effected by how bright the source is, as illustrated by the shift to the right for model m7a0.9-HR when *D* = 100 Mpc. We compiled estimates of the jet lengths for each image in Table 3.

**Figure 4.** Here we show the relative difference (open circles) and the error within 3 standard deviations (horizontal bars) between the jet lengths obtained from the ray-traced images (*l*jet) and the reconstructed images (*l*jet,rec). Models m7a0.0-HR (**left panel**) and m7a0.9-HR (**right panel**) are shown. For each choice of viewing angle *θ*, distance to the source *D*, and R, we show the data at *t* = 38,000 *tg* (blue) and *t* = 78,000 *rg*. We indicate the viewing angle next to each pairing of error bars. For each model, except where explicitly indicated to be otherwise, we show *D* = 10 Mpc and R = 1. In general, there is a very small relative difference when the jet is viewed edge on. As the viewing angle approaches face on or the jet is placed at a larger distance, the relative difference increases along with the error. In all but one image, even when the relative difference shifts away from zero, the zero relative difference line is within 3 standard deviations, which suggests excellent agreement between jet lengths derived from the raw ipole and reconstructed images.

Since we lack centering information in the reconstructed images, we chose to estimate the jet velocity perpendicular to the line of sight by assuming both the top and bottom jet have the same speed. Then, the velocity of the jet in the source's frame is:

$$v = \frac{1}{2} \frac{|l\_{\text{jet}}(t\_2) - l\_{\text{jet}}(t\_1)|}{t\_2 - t\_1}. \tag{11}$$

Note that we used the same expression to derive the velocity in the reconstructed images (*v*rec) but replaced *l*jet with *l*jet,rec in Equation (11). Very long baseline interferometry (VLBI) observations can provide the apparent motion of the jet, which may be superluminal due to relativistic effects. To account for this, we assumed the jet has a velocity of *v*rec and then estimated the apparent velocity (*v*app,rec) via the time in the observer's frame (*t* 1,*t* 2):

$$\upsilon\_{\rm app,rec} = \frac{1}{2} \frac{|l\_{\rm jet,rec}(t\_2') - l\_{\rm jet,rec}(t\_1')|}{t\_2' - t\_1'} = \upsilon\_{\rm rec} \left(1 - \frac{\upsilon\_{\rm rec}}{c} \frac{\cos \theta}{\sin \theta} \right)^{-1}. \tag{12}$$

We used the relationship

$$t\_2' - t\_1' = (t\_2 - t\_1) \left( 1 - \frac{v\_{\text{rec}}}{c} \frac{\cos \theta}{\sin \theta} \right) \tag{13}$$

in the last expression. Note that the division by sin *θ* was to account for the fact the *v*rec measures the velocity parallel to the line of sight with no time delay effects while we require an estimate of the velocity along the jet axis (which we can obtain since the geometry is fully known).

**Table 3.** Here we tabulate the estimated jet length for each model at each time for various choices of the distance *D*, viewing angle *θ*, and plasma temperature ratio R. We compare the jet length as computed from the base ipole image (*l*jet) and the reconstructed image (*l*jet,rec).

We used the data at *t*1 = 38,000 *tg* and *t*2 = 78,000 *tg* for each model to estimate the jet velocity. We tabulated the estimated velocity and apparent velocity for each model in Table 4. We found excellent agreement between the velocities derived from the raytraced and reconstructed images. It is interesting to note the apparently faster jet for R = 20. We suspect the increased speed is due to the emitting material being dominated by material near the jet axis rather than some of the slower moving material around the jet head, which does not produce much emission at 230 GHz as the electron temperature is reduced. The more powerful jet model m7a0.9-HR demonstrates that such jets may appear as superluminal sources, as we found a maximum *v*app,rec ≈ 1.526*c* at *θ* = 10◦.

**Table 4.** Here we tabulate the estimated jet velocity for each model for different choices of the distance *D*, viewing angle *θ*, and plasma temperature ratio R. The velocities shown were calculated using the base ipole images (*v*), the reconstructed images (*v*rec), and the reonstructed images, while accounting for the possibility of superluminal motion (*v*app,rec).

#### **4. Discussion**

#### *4.1. Extracting Jet Physics from VLBI Images*

A key feature of the jets in our models is the bright "bubbles" (or knots) of 230 GHz emission, which appear to correlate with recollimation shocks. Such structures have been seen in VLBI images of various AGN jets [48–50]. Previous simulations of jets in various astrophysical contexts have demonstrated that recollimation occurs when there is a pressure mismatch between the jet and the surrounding medium, which could be a static atmosphere or a slower-moving jet sheath [51–54]. The number of recollimation shocks along the jet axis is dependent on the properties of the jet and medium. It is therefore possible that direct VLBI of TDE jets will allow in-depth modeling of jet launching and could also aid in constraining the properties of the surrounding medium. For instance, one work successfully applied simulations of MHD jets to constrain properties of BL Lacartae, which is a blazar jet with recollimation features [55].

We suggest that a similar approach may be applied in TDE jets. With a suitable exploration of the parameter space, it is conceivable that an analysis similar to that of [55] could be applied to TDE jets in cases where VLBI is possible. A broader exploration of TDE jets through various simulation methodologies is strongly suggested. We plan to explore the effects of the ambient medium, magnetic field strength, and disk accretion rate on the jet properties in the case of a SANE, super-Eddington disk in a future work.

#### *4.2. Proposed Observational Methodology*

Our synthetic ngEHT observations demonstrate that a 24-hour observation may be sufficient to study both the structure and/or motion of newly born TDE jets. This is much shorter than the fallback time, which is on the order of month(s), and the duration of radio emissions in several TDEs, which can sometimes be visible for years [12]. Our suggested observational strategy is conducting rapid followup of newly discovered optical/X-ray TDEs when they are near the peak of their emission in order to study both the early- and late-time properties of their jets (if present). Observations with a single telescope at 230 GHz could be conducted to search for TDEs emitting in the radio. If emission at 230 GHz is detected, we suggest that the ngEHT conduct VLBI follow-up of targets within no more than a month. Our imaging simulations were done assuming a full ngEHT array consisting of the 2022 EHT stations plus 11 additional sites, but depending on the target, not all sites may need to be available in order to obtain a high-fidelity image reconstruction.

Unlike many other EHT/ngEHT targets, TDEs will appear randomly across the sky, and the ngEHT will need to be capable of follow-up observations in the order of a week to weeks. Our current modeling of jets and outflows from super-Eddington disks is too sparse to let us make predictions regarding how long the jets will be visible at 230 GHz. However, if radio TDE observations are any indicator, emissions may persist for many months [12].

TDEs provide an excellent laboratory for studying jet/accretion/black hole systems across a wide range of accretion states over a relatively short period of time (∼1-a few years). In several cases, TDEs have shown state transitions after several hundred days in the X-ray range, which are likely associated with the evolution of the disk as the mass accretion rate declines. Stone and Metzger [29] argue, for instance, that the transition from a thick, super-Eddington disk to a thin disk can explain the jet shut-off in jetted TDEs, such as *Swift* J1644+57 and *Swift* J2058+05 [56,57], but recent simulations [58,59] demonstrated that MAD is possible even for thinner accretion disks. As such, long-term VLBI monitoring is strongly suggested. as this would allow for (1) the radio-submm emission of the outflows to be characterized and compared to the behavior of the accretion flow, and (2) the direct study of how the jet evolves morphologically as the disk state changes.

Another attractive potential target which we have yet to attribute a self-contained study to is jetted TDEs. These TDEs are extremely rare, and current observations suggest only about 1% of all TDEs will produce powerful relativistic jets. These jets will produce extremely bright radiation in the X-ray and the radio-submm regions. However, most have been distant due to the lower probability of their occurrence. Should a jetted TDE occur nearby enough for VLBI to resolve the jet, we strongly suggest such jets be treated as targets for the ngEHT.

Lastly, a recent TDE AT2018hyz showed a late outflow (∼3 years after the initial outburst) and brightened in the radio over several hundred days [60]. Unlike many other radio TDEs, a relatively bright 240 GHz emission was detected. If AT2018hyz is in fact a jet instead of a spherical outflow, Cendes et al. [60] estimates that the velocity could reach 0.6c. AT2018hyz is a relatively nearby TDE at ∼204 Mpc, but the flux density at the time of detection (∼0.2 mJy at 240 GHz) makes it too dim for ngEHT follow-up. However, placing AT2018hyz at ∼50 Mpc would shift the flux density to ∼5 mJy, which is the minimum estimated flux density required for an ngEHT VLBI detection. Future TDEs will likely be monitored across the radio-submm, so nearby targets of opportunity such as late radio TDEs like AT2018hyz should be considered, should they show significant radio emission.

While the EHT operated using dedicated time, the ngEHT will include more arrays and some telescopes may be sequestered for targeting of opportune VLBI observations (Doeleman et al., in prep). Development and proposals for the ngEHT are still underway, but we expect that fast (within 1-several weeks) response coverage of sources of interest in order to conduct VLBI of TDEs (and other transient sources) should be possible (Shepard Doeleman, private communication). The Atacama Large Millimeter Array (or ALMA), the most sensitive antenna, will maintain sustained coverage [40]. As such, we anticipate that at least unresolved radio data can be captured, which would prove equally invaluable to the study of radio TDEs.

#### **5. Conclusions**

In this work, we have demonstrated through a synthetic imaging analysis that TDE jets resembling the GRRMHD models presented in Curd et al. [14] are compelling ngEHT targets. We also confirmed that the detection limits considered in Curd et al. [14] are roughly applicable, as m7a0.0-HR did not produce detectable emission at a distance of 100 Mpc in our imaging analysis.

Various shock features in the jet are visible for the 10 Mpc images we consider, and studying the jet morphology in these cases could aid in characterizing the environment of the BH. Most TDEs that occur during the ngEHT mission will be farther away, but the apparent motion, which may be superluminal, can be extracted in such cases.

We suggest that the ngEHT be utilized for radio follow-up of TDEs. Our models study the birth of a TDE jet in the first ∼ 48 days after the disk forms under the assumption that the disk is SANE, super-Eddington, and threaded by a dynamically important magnetic field. However, TDEs which occur nearby and have jet properties similar to those of jetted TDEs, such as *Swift* J1644+57 or AT2018hyz, may provide interesting targets.

**Author Contributions:** Conceptualization, R.E. and R.A.; Formal analysis, B.C. and F.R.; Investigation, B.C.; Writing—original draft, B.C., F.R. and R.A.; Writing—review & editing, B.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** Brandon Curd was supported by NSF grant AST-1816420, and made use of computational support from NSF via XSEDE/ACCESS resources (grant TG-AST080026N). Razieh Emami acknowledges the support by the Institute for Theory and Computation at the Center for Astrophysics and grants 21-atp21-0077, NSF AST-1816420, and HST-GO-16173.001-A for very generous support. Freek Roelofs was supported by NSF grants AST-1935980 and AST-2034306. This work was supported by the Black Hole Initiative at Harvard University, made possible through the support of grants from the Gordon and Betty Moore Foundation and the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the Moore or Templeton Foundations.

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

**Acknowledgments:** We graciously acknowledge Koushik Chatterjee, Lani Oramas, Joaquin Duran, and Hayley West for fruitful discussions that were helpful in the preparation of this work.

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

#### **Appendix A. Full Image Library**

In Figures A1–A4, we show the full library of images analyzed in this work. All of the base images were ray traced at *ν* = 230 GHz and then convolved with a Gaussian beam with a FWHM of 20 μas. Similarly, we blurred the reconstructed image using the same beam for comparison. Note that for Figures A2 and A4, we also shifted the reconstructed image to be approximately centered for comparison with the base image. In general, bright features are represented quite well in the reconstruction; however, some noise is introduced in dimmer sources (i.e., see the reconstruction of m7a0.9-HR at *D* = 100 Mpc).

**Figure A1.** Here we show the full library of base ipole images for model m7a0.0-HR. The time of each column is indicated at the top; the distance *D* and plasma temperature ratio R for each row are indicated on the right. The angle of the observer relative to the jet axis is indicated above each set of two rows. Each image spans 550 × 770 μas2 and is blurred by convolving the base image with a Gaussian beam with a FWHM of 20 μas. The color scale is logarithmic, spanning three orders of magnitude, and each image uses the same maximum for the intensity scale.

**Figure A2.** The same as Figure A1 but showing the reconstructed images. Note that the intensity scale for each panel is the same as the corresponding panel in Figure A1 for comparison.

**Figure A3.** Here we show the full library of base ipole images for model m7a0.9-HR. Each image in the top and middle rows spans 1000 × 1400 μas2, and the bottom row spans 170 × 238 μas2. Each image was blurred by convolving the base image with a Gaussian beam with a FWHM of 20 μas (indicated by the white circle in the bottom right panel). The time of each column is indicated at the top; the distance *D* and plasma temperature ratio R for each row are indicated on the right. The angle of the observer relative to the jet axis is indicated above each set of two rows. The color scale is logarithmic, spanning three orders of magnitude in each image. We used the same color scale for the *D* = 10 Mpc images, but reduced the maximum by an order of magnitude in the *D* = 100 Mpc images to better show the image features.

**Figure A4.** The same as Figure A3 but showing the reconstructed images. Note that the intensity scale for each panel is the same as the corresponding panel in Figure A3 for comparison.

#### **Appendix B. Fitting Procedure for Jet Head Position**

Here we describe the algorithm implemented to estimate the positions of the top and bottom jet heads in each image. In order to track the jet's motion, we first smoothed the data (either the base image or the reconstructed image) with a Gaussian beam which assumed an angular resolution for the VLBI observations of Δ*θ* = 20 μas. We then binned the data along the symmetry axis of the jet (*y*) by summing the flux along each row (*x*).

Since the images appear to be roughly Lorentzian, we first attempted to fit a double Lorentzian function of the form:

$$f\_d(y) = a\_1 \frac{w\_1^2}{(y - y\_{\text{cen},2})^2 + w\_1^2} + a\_2 \frac{w\_2^2}{(y - y\_{\text{cen},2})^2 + w\_2^2} \tag{A1}$$

where *a* is the amplitude, *w* is the width, and *y*cen is the center defining the curve. We took *y*cen as a measure of the jet head location. If this fitting procedure did not produce a good fit, we found the peaks by performing a two-step fitting procedure in which we fit a single Lorentzian:

$$f\_1(y) = a\_1 \frac{w\_1^2}{(y - y\_{\text{cen},1})^2 + w\_1^2} \tag{A2}$$

and then subtracted the fit *f*1(*y*) from the data and then fit the second peak with:

$$f\_2(y) = a\_2 \frac{w\_2^2}{(y - y\_{\text{cer},2})^2 + w\_2^2}. \tag{A3}$$

We implemented the Python package SciPy [61] to optimize the curve(s) and estimate the jet head positions and errors. We show an example of the data and the Lorentzian fit for both the base image and the reconstruction in Figure A5. The reconstruction tends to be a bit broader, but the fitting procedure works equally well for all images and reconstructed images.

**Figure A5.** We demonstrate the fit performance using a ray-traced image of model m7a0.9-HR at *t* = 78, 000 *tg*, *θ* = 90◦, *D* = 10 Mpc, and R = 20. We show the *x*-binned data and a double Lorentzian (Equation (A1)) fit for (**a**) the blurred base image and (**b**) the blurred reconstruction.

#### **Notes**

1 https://github.com/Smithsonian/ngehtsim (accessed on 1 June 2022).

#### **References**

### *Article* **Spectral Line VLBI Studies Using the ngEHT**

**Dong-Jin Kim 1,2,\*,† and Vincent Fish 1,†**

**Abstract:** Spectroscopy in the mm/sub-mm wavelength range is a powerful tool to study the gaseous medium in various astrophysical environments. The next generation Event Horizon Telescope (ngEHT) equipped with a wide-bandwidth backend system has great potential for science using high angular resolution spectroscopy. Spectral line VLBI studies using the ngEHT will enable us to scrutinize compact astrophysical objects obscured by an optically thick medium on unprecedented angular scales. However, the capabilities of ngEHT for spectroscopy and specific scientific applications have not been properly envisioned. In this white paper, we briefly address science cases newly achievable via spectral line VLBI observations in the mm/sub-mm wavelength ranges, and suggest technical requirements to facilitate spectral line VLBI studies in the ngEHT era.

**Keywords:** very long baseline interferometry (1769); radio astronomy (1338); millimeter astronomy (1061); submillimeter astronomy (1647); radio telescopes (1360); high angular resolution (2167)

### **1. Introduction**

The next generation Event Horizon Telescope (ngEHT) is primarily designed to image the continuum emission from the vicinity of the photon ring of supermassive black holes (SMBHs) with a better imaging performance than that of the current EHT array in terms of sensitivity and sampling in the spatial-frequency domain. A major upgrade of the ngEHT would be a new front and back end system, supporting a bandwidth of 8 GHz per sideband, dual polarization, and simultaneous dual band 230/345 GHz capability with a recording rate of 256 Gbps [1]. This new capability provokes interest for spectral line studies using the ngEHT for two reasons: (1) VLBI observations with such a wide bandwidth enable spectral line studies without additional frequency tunings at each station, and (2) the ngEHT offers several orders of magnitude higher angular resolutions than those available with single-dish telescopes or connected arrays, such as SMA, NOEMA, and ALMA. Microarcsecond scale angular resolution in spectral line studies would open new doorways to explore intriguing physical phenomena in compact astrophysical systems, such as the core of star-forming regions (SFRs), circumstellar envelope (CSE) of evolved stars, and circumnuclear disk (CND) of active galactic nuclei (AGNs).

Spectral line studies with a mm/sub-mm VLBI array and their possible scientific impacts have attracted the attention of researchers in recent years, especially due to the significant increase in sensitivity provided by the Atacama Large Millimeter/submillimeter Array (ALMA) Phasing System [2–5]. ALMA has already offered a prototype spectral line VLBI capability in Band 3 (3.5 mm) in conjunction with the Global mm-VLBI Array (GMVA) and anticipates offering more flexible spectral line VLBI modes in additional bands in the near future. However, much work is still needed to utilize its unique capability for VLBI studies. One of the most pressing tasks is to categorize possible candidates and types of spectral lines suited for mm/sub-mm VLBI experiments. In this context, here we highlight the types of spectral line sources available for VLBI studies and their potential use in various scientific applications.

**Citation:** Kim, D.-J.; Fish, V. Spectral Line VLBI Studies Using the ngEHT. *Galaxies* **2023**, *11*, 10. https:// doi.org/10.3390/galaxies11010010

Academic Editor: Luigina Feretti

Received: 18 November 2022 Revised: 30 December 2022 Accepted: 1 January 2023 Published: 6 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/).

#### **2. Molecular Maser Lines**

Microwave amplification by stimulated emission of radiation (maser) is a non-thermal process in the interstellar medium (ISM) that can result in remarkably bright spectral line emission. Molecular maser lines were discovered in both galactic and extragalactic sources [6–8]. The physical sizes of masing clouds measured by interferometric observations range from sub-AU to a few AU for galactic masers and sub-pc to pc scales for extragalactic masers [9,10]. Considering the distance and physical size of maser sources, their angular scale is typically a milli-arcsecond (mas). For instance, a galactic maser source with a physical scale of 1 AU corresponds to an angular scale of 1 mas at a source distance of ∼1000 pc. In this case, a proper motion of 10 km s−1 (∼2.1 AU yr−1) corresponds to 2.1 mas yr−1 with an additional annual parallax of up to 1 mas. Monitoring VLBI observations at cm/mm wavelengths over a few months allow us to measure the proper motion or annual parallax of nearby galactic maser sources. VLBI observations at shorter wavelengths (mm/sub-mm) achieve higher angular resolutions, which can enlarge the number of spectral line source for astrometry, such as more distant galactic maser sources (>a few kpc) or sporadic maser features that last only a few weeks or less. In addition, mm/sub-mm maser transitions are theorized to have distinct excitation conditions, and hence they can unveil different part of the gaseous medium that cannot be traced by conventional maser transitions available at cm wavelengths [11,12].

#### *2.1. Circumstellar Envelope of AGB Stars*

In the late stages of stellar evolution, the outer layers of intermediate mass stars (1–8 M ) expand to a few AU in radius, pulsating with a period of up to several years [13]. They are called asymptotic giant branch (AGB) stars, and their strong stellar winds (10–30 km s−1) form a thick circumstellar envelope (CSE) [14]. AGB stars show morphological transitions of CSEs on relatively short time scales (105–106 yr), evolving into planetary nebulae (PNs) [14]. Despite such drastic changes, CSEs are known as the cradle of complex molecular species owing to their dense and cool environments, which are well-suited for molecular synthesis [15]. The characteristics of CSEs vary depending on their chemical and physical environments. For instance, the chemical composition of AGB stars determines the strength of stellar winds driven by heated dust [16]. Meanwhile, the presence of a companion star or giant planet likely plays an important role in the morphological evolution of CSEs, such as spiral patterns, collimated jets, and multi-polar structures [17,18]. To investigate the evolution of CSEs, multiple different types of AGB stars need to be studied separately, tracing stellar winds from the inner to the outer part of the CSEs.

In oxygen-rich (O-rich) AGB stars, SiO masers show a symmetric ring-like structure within a few stellar radii (*R*∗), while H2O masers are located farther out up to a few tens of stellar radii with either symmetric or asymmetric spatial distributions [19–21]. To understand the development of asymmetric structure, it is necessary to probe the region between symmetric SiO masers and asymmetric H2O masers clouds. According to radiative transfer models, the sub-mm H2O maser emission at 321 GHz traces denser and warmer regions (2–5 *R*∗) compared to physical conditions required for masing at 22 GHz (2–10 *R*∗) [11,22]. VLBI observations of 321 GHz H2O masers using the 345 GHz ngEHT receivers can test this theoretical prediction, constraining the size of masing clouds down to sub-AU scales (e.g., 15 μas ∼ 0.015 AU at 1 kpc distance) and also fill in the gap in our understanding of the development of stellar winds in CSEs. In contrast, carbon-rich (C-rich) AGB stars mainly form carbon-based molecules in their CSEs, and thus maser lines from oxygen-bearing molecules (e.g., H2O, SiO, and OH) are barely detected in their CSEs. Spectral line surveys have detected HCN and SiS maser lines in the sub-mm wavelength range, and these lines seem widespread in C-rich AGB stars [23,24]. The CSEs of C-rich AGB stars are supposed to generate stronger stellar winds than those in O-rich AGB stars, but much is unknown about their kinematics due to the lack of high-resolution observations [16]. VLBI observations of HCN masers in C-rich AGB stars would shed light on this matter by spatially resolving maser spots and measuring their proper motions.

#### *2.2. Star-Forming Regions*

Gravitational collapse of clouds in the ISM triggered by internal or external physical conditions forms dense and compact clumps suited for star formation. Since star-forming regions (SFRs) are optically thick, much is still unknown about the process of star formation and differences in regions of high- and low-mass star formation [25] (and references therein). A variety of maser lines observed in massive SFRs (e.g., OH, H2O, class I and II CH3OH masers) are unique tools to classify different stages or classes of star formation. For instance, two distinct types of methanol masers excited collisionally (class I) or radiatively (class II) trace outer and inner regions of SFRs, respectively [26,27]. VLBI observations of class II methanol masers have shown complex inner structure in SFRs, including a ring-like structure [28–30]. On the contrary, class I methanol masers have revealed outflow features at a large distance from the central object [31]. Proper motion measurements via multiepoch VLBI campaigns show diverse kinematics within scales of a few hundred AU, such as infall, outflow, and rotation. On the other hand, maser lines are rarely detected in low- and intermediate-mass SFRs compared to those in high-mass SFRs [32,33]. H2O and methanol masers detected in some low-mass SFRs likely trace the inner parts of outflows [34]. Observations in the sub-mm wavelength range have detected new types of H2O and CH3OH maser lines in SFRs, but they have not been spatially resolved due to insufficient angular resolutions. [35–38]. Their excitation conditions are different from those of cm maser lines, suggesting the possibility to probe different parts of SFRs via VLBI observations [27]. Kinematics of clumpy molecular clouds revealed by VLBI observations will advance our understanding on the process of star-formations.

#### *2.3. Circumnuclear Gas of AGNs*

Extraordinarily bright maser lines (e.g., H2O and OH) have been detected in a number of active galactic nuclei (AGNs) at cm wavelengths. VLBI observations show that they originate in either the accretion disk or outflow in the circumnuclear region of AGNs [10]. Measurements of the angular scale and rotational curve of the disk masers in AGNs provide an accurate measurement of the mass of the SMBHs [39,40]. The host galaxies of H2O megamaser sources are commonly Seyfert galaxies, having relatively small SMBH masses (106–108 M ) [41]. H2O masers are very rare in elliptical galaxies. NGC 1052 and NGC 4261 are two such cases where the line profiles and VLBI observations indicate that the 22 GHz H2O masers are likely associated with the ambient ISM interacting with radio jets [42,43]. Sub-mm H2O masers have been detected in a number of AGNs, but their fine structure has not been spatially resolved [44–46]. VLBI observations of those sub-mm masers have potential to explore nuclear region of AGNs, such as molecular outflows, inner structure of molecular accretion disks, or jet-ISM interactions. We note that VLBI observations of 22 GHz H2O masers towards the nearby AGN NGC 3079 show clumpy clouds with size in the range of 0.002 to 0.02 pc [9]. The ngEHT observations at 345 GHz offer an angular resolution of ∼15 μas (15 μas = ∼0.0015 pc at 20 Mpc distance), which can resolve a compact circumnuclear gas cloud in a nearby AGN.

#### **3. Atomic Maser Lines**

Ionized H II regions are known as the site of radio recombination lines (RRLs). In ionized regions, hydrogen RRLs sometimes show maser features in the sub-mm wavelength range [47]. Thanks to their high brightness temperatures, atomic maser lines are promising candidates for spectral line VLBI observations. The ionized region is commonly identified by optical emission lines (e.g., H*α*), but a thick ambient medium often absorbs or reflects the optical emission from the ionized H II regions. Atomic maser lines in the radio wavelength range are optically thin, and thus help scrutinize gas kinematics in the obscured ionized region.

#### *3.1. Post-AGB Stars and Pre-Planetary Nebulae*

At the end of the AGB phase, the central star starts to ionize its surrounding materials. The ionized core region is the key to understand the development of asymmetric multipolar jets or winds seen in PNs, but a dense dusty envelope heavily obscures the central region, impeding investigations of the core region of CSEs during the post-AGB and preplanetary nebulae (pPNe) phases. *η* Carinae is the first RRL maser source detected in AGB stars [48]. ALMA observations of H30*α* emission show that a slightly blueshifted narrow line feature in the H30*α* spectrum is extended, but broad line features are compact [49,50]. More RRL maser sources have been detected in post-AGB stars, and they are characterized by different line profiles [51]. Follow-up high-resolution observations towards those RRL masers are necessary to clarify their origins and possible connections with bipolar/multipolar outflows in the early phase of post-AGB stars.

#### *3.2. Ultra-Compact H II Regions of Massive Stars*

Strong radiation from newborn massive stars forms ionized regions. In the process of massive star formation, the ambient medium is thick and dense enough to cover the ionized region. This type of object is known as an ultra-compact (UC) H II region. Unveiling the kinematics of UC H II region is one of the most important topics in the study of starformation. Theoretical studies point out that feedback from young stars prevents further matter accretion in high-mass star formations. However, somehow high-mass stars do exist in various star-forming regions. Studying hyper-compact (HC) H II regions might hint how matter keeps accreting during the early stage of high-mass star formation [52,53].

RRLs appear in ionized gas, such as UC H II regions. Since ionized gas also radiates continuum emission, line/continuum ratio is an important parameter to detect RRL. At submm wavelengths, the line/continuum ratio increases remarkably and maser emission can arise in certain local conditions. MWC 349 is the high-mass star-forming region, showing RLL maser at sub-mm wavelengths [54]. The H30*α* maser emission in MWC 349 traces the innermost regions of the cocoon. Submillimeter Array (SMA) observations of RRL maser reveal an ionized accretion disk in MWC 349 [55]. More RRL masers sources have been detected in UC H II region. The double-peaked RRL maser spectrum with very broad line widths (up to a few hundreds km s−1) suggests that RRL masers are associated with a biconical or collimated radio jet [56,57]. High-angular resolution VLBI observations towards RRL maser lines are needed to unveil the nuclear region of very young massive stars.

#### **4. Molecular Absorption Lines**

Molecular clouds become visible in absorption against the featureless continuum emission of the AGN jet. Absorption studies allow us to detect even low abundance molecular species regardless of their cosmological distance if the continuum jet is bright enough, effectively using the jet as a strong flashlight to illuminate the material. The absorption technique can be used to access any intervening molecular clouds lying on the line of sight. Spectral-line VLBI observations towards molecular absorption lines spatially resolve the structure and kinematics of the obscuring molecular clouds, and therefore can clarify the origin of obscuring molecular clouds. Gravitationally lensed blazars (GLB) and circumnuclear gas in radio AGNs are well-known such systems, showing molecular absorption lines from the intervening molecular clouds. Previous absorption studies focused on cm wavelengths hint interaction between radio jets and ambient medium [58,59], but observable molecular species and angular resolutions are limited. In the sub-mm range, various molecular species in different chemical environments are detectable in absorption with improved angular resolutions. Therefore, mm/sub-mm absorption line studies are desirable to advance our understanding on the chemistry of extragalactic ISM and kinematics of circumnuclear molecular gas in radio AGNs.

#### *4.1. Gravitationally Lensed Blazars*

PKS 1830-211 is one of the most well-known gravitationally lensed blazar systems, consisting of a distant blazar at a redshift of 2.5 and a face-on spiral galaxy at a redshift of 0.89 [60]. The blazar serves as a background continuum source for the intervening galaxy, and thus various molecular lines have been detected in absorption. With bright continuum, absorption line observations allow us to detect rare molecular species compared to observations using emission lines. Indeed, over 50 molecular species have been detected in the intervening spiral galaxy via absorption lines [60–64]. However, observations using single-dish telescopes or connected arrays tend to underestimate absorption depths due to low filling factors, since the beam size is larger than the size of absorbed area [65]. Absorption line observations using VLBI can spatially resolve individual continuum components related to absorption features, hence providing a better physical constraint on the obscuring molecular gas. Furthermore, VLBI monitoring of molecular absorption lines can provide a unique tool to study a structural variation of the background blazar through absorption line variability at multiple frequency bands [66].

#### *4.2. Circumnuclear Gas in AGNs*

AGNs are among the most powerful energy sources in the Universe. The radio jet outflows triggered by AGN activity are believed to play an important role in galaxy formation and evolution [67–69]. However, the nature of different types of radio jets and their influence on their host galaxies (i.e., radio jet feedback) are still poorly understood. As a matter reservoir in the vicinity of the central engine, the circumnuclear gas is the key to study different accretion environments and radio jet feedback in AGNs harboring radio jets (i.e., radio AGNs). Although radio AGNs tend to be gas-poor and distant, molecular absorption line observations using bright continuum radio jets enable us to study the circumnuclear gas. Sparse molecular species can be detected in absorption, providing a unique tool to study the small-scale (<pc) structure of the circumnuclear gas via spectralline VLBI observations. Such high resolution observations are particularly important to reveal the onset of radio jet feedback or the fate of a fueling flow in the vicinity of the central engine. Despite the advantages of absorption line study, only a handful of radio AGNs have been identified as molecular absorption systems. In the sub-mm wavelength range, plenty of molecular lines are available for absorption line observations. Indeed, recent ALMA observations shows molecular absorption lines in radio AGNs [70–72]. VLBI observations towards sub-mm molecular absorption lines are the next step to clarify the origin of obscuring gas in radio AGNs and their roles in AGN fueling and feedback.

#### **5. Technical Overview for Spectroscopy**

#### *5.1. Spectral Line Sensitivity*

New ngEHT stations may have an aperture size of 10 m or less, resulting in a lower spectral line sensitivity on new baselines than for the current EHT array. Assuming a conventional quantization loss in the data recording process (*ηs* = 0.8), the expected spectral line sensitivities of ngEHT6 *m*-ALMA and ngEHT6 *m*-ngEHT6 *m* baselines are 88 mJy and 1040 mJy, respectively with an integration time of 15 min and a spectral resolution of 1 MHz at 345 GHz (see Table 1). An order of magnitude higher spectral line sensitivity can be achievable by adding ALMA as an anchor station. Therefore, the optimal configuration of the ngEHT for spectral line VLBI observations would be including NOEMA and ALMA as anchor stations. This will provide exquisite spectral line sensitivity for sources in both the northern and southern hemispheres. Despite its relatively small aperture size, additional ngEHT stations will eventually improve the total spectral line sensitivity of ngEHT array by increasing the number of baselines (N*baseline* = N*ant*(N*ant*−1)/2). Considering the expected spectral line sensitivity, a minimum peak flux density of a narrow maser emission line (FWHM < 10 km s−1) would be >10 Jy and >20 Jy at 230 GHz and 345 GHz, respectively for imaging the spectral line with the standalone ngEHT array. Table 2 shows a representative sub-mm spectral line sources mentioned in Sections 2–4 with line parameters.

**Table 1.** Single-dish spectral line sensitivity of ngEHT stations with an integration time of 15 min at 230 GHz and 345 GHz, respectively.

Note: Aperture efficiencies are 0.8, 0.67, 0.8 at 230 GHz and 0.8, 0.62, 0.7 at 345 GHz for ngEHT, NOEMA, and ALMA stations, respectively. The same amount of precipitable water vapor (PWV) was assumed at all stations (3 mm at 230 GHz and 1 mm at 345 GHz). Reference: https://www.iram.fr/GENERAL/NOEMA-Phase-A.pdf. Accessed on 30 December 2022.

**Table 2.** Different types of spectral line sources at sub-mm wavelengths.

Reference: a [73–75]. b [76] c [46]. d [49,50]. e [54,55]. f [77] g [71].

Adding an anchor station, such as ALMA and NOEMA will lower those limits by a factor of 12. Observations of a broad (FWHM > 10 km s−1) absorption or emission line with a moderate spectral resolution will enlarge the number of detectable spectral line sources with the ngEHT. In previous H2O maser surveys in the cm/mm wavelength ranges, about 2300 H2O maser sources have been detected (SFR: ∼1400, AGB: ∼700, and extragalactic: ∼180) [38,78]. Assuming 5% detection rate of bright (>20 Jy) sub-mm H2O maser among known maser sources, about 115 sources will be feasible targets for the ngEHT. Observations of other sub-mm maser species, such as CH3OH and SiO, will increase the number of available targets to a few hundreds. Meanwhile, an increased coherence time by frequency phase transfer (FPT) technique will further improve the spectral line sensitivity of ngEHT [79,80], but a detailed investigation is needed to evaluate the capability of FPT for spectral line VLBI observations.

#### *5.2. Frequency Coverage*

A simultaneous triple-band (3 mm, 1 mm, and 0.8 mm) receiver system has been considered for the ngEHT to apply the FPT technique. Simultaneous multi-band VLBI observations can yield several atomic and molecular lines in a single run. This unique capability is beneficial for spectral line VLBI studies in various aspects. For instance, a combined study using multiple VLBI images of different molecular maser lines would provide a comprehensive view of a masing region [21,81]. In case of molecular absorption line studies, observing multiple rotational lines provides solid constraints on the physical quantities of the obscuring gas, such as temperature and column density. We note that the observing frequency of 1 mm and 0.8 mm receivers should have integer frequency ratios of the 3 mm receiver for FPT. This technical requirement would help simultaneous observations of several rotational lines at different bands, such as HCN, HCO+, SiO, and CO. However, simultaneous multi-band observations will result in narrower bandwidth at each band than single-band observations due to a limited recording rate of the backend system. Therefore, we need to consider several spectral line observation modes for individual science projects. For instance, non-integer and wide bandwidth observations are preferred

for RRL maser line observations. On the contrary, integer and multi-band observations are preferred for molecular maser or absorption line observations. In practice, offering a flexible frequency setup would not be feasible in the early phase of array operation. Therefore, we suggest two optimal frequency tunings for 8 GHz and 16 GHz bandwidth observations, respectively considering continuum observations with FPT. Specific frequency tunings and distribution of spectral lines are present in Table 3 and Figure 1.

**Figure 1.** Distribution of spectral lines over the frequency range of 1.3 mm and 0.8 mm. Each color indicates the type of spectral lines. Frequency coverage of 8 GHz and 16 GHz observations is marked with two different frequency tunings (low and high).

**Table 3.** Suggested frequency tunings with different bandwidths (8 GHz and 16 GHz) at each band for spectral line VLBI observations.

#### **6. Conclusions**

The ngEHT will be equipped with multi-band receivers supporting simultaneous wideband observations. It will not only improve sensitivity for continuum imaging, but also open new opportunities for high angular resolution spectroscopy. In this white paper, we have outlined the types of spectral lines suited for VLBI studies and briefly described several scientific applications newly achievable with the ngEHT. Spectral lines observable with the 230 and 345 GHz receivers are listed in Table 4. Masers appear in various astrophysical environments, and they trace certain chemical and physical environments depending on atomic/molecular species and transitions. Atomic and molecular maser lines observable in the sub-mm wavelength range possibly unveil regions unprobed by previous cm-maser observations. The ngEHT can be used to image sub-mm masers, investigating their peculiar physical environments and kinematics.

**Table 4.** Selected spectral lines near 230 GHz and 345 GHz and their energy level above ground state. Three different types of spectral lines available for VLBI observations are listed with physical origins.

Note: Rest frequencies were obtained from the Splatalogue database available at www.splatalogue.net. Accessed on 30 December 2022.

However, VLBI experiments with maser lines could be challenging due to the limited spectral line sensitivity of the ngEHT and the rarity of sub-mm masers compared to cm/mm masers. Absorption line observations using a bright background continuum source (e.g., radio jet) are an alternative approach to enable spectral line VLBI studies towards common molecular clouds regardless of their distance. Molecular absorption from circumnuclear gas in AGNs or the intervening galaxy of a lensed blazar are promising candidates for spectral

line VLBI studies. Lastly, we have highlighted the capabilities of ngEHT for spectroscopy. Considering the rest frequencies of atomic and molecular transitions, we suggest frequency tunings to facilitate the operation of ngEHT in spectral line mode.

**Author Contributions:** Conceptualization, D.-J.K.; methodology, D.-J.K.; writing—original draft, D.-J.K.; writing—review and editing, V.F.; supervision, V.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Science Foundation through grant numbers AST-1935980 and AST-2034306.

**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.

MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Galaxies* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/galaxies

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.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-9355-5