Modeling Nearby Low-Luminosity Active-Galactic-Nucleus Jet Images at All VLBI Scales

: Relativistic jets from nearby low-luminosity active-galactic-nuclei (LLAGN) were observed by Very-Long Baseline Interferometry (VLBI) across many orders of magnitude in space, from milliparsec to sub-parsec scales, and from the jet base in the vicinity of black holes to the jet collimation and acceleration regions. With the improved resolution for VLBI observations, resolved VLBI jet morphologies provide valuable opportunities for testing and constraining black hole jet physics. In this review, we summarize and discuss the current progress of modeling nearby LLAGN jet images from horizon scales to large scales, including the construction of jet models and the assumed emission details. Illustrative examples for jet image modeling are also given to demonstrate how jet image features may vary with the underlying physics.

Observed spectra from black hole accretion jet systems include the contribution from both the accretion and the jet. The spectra of LLAGN are usually featured with a near millimeter bump (∼10 2 GHz) [15][16][17][18], which is usually referred to as the thermal synchrotron emission contributed by the innermost region of a radiatively inefficient accretion flow (RIAF), a type of black hole accretion at a low accretion rate (Ṁ < 0.01Ṁ Edd , wherė M Edd ≡ 10L Edd /c 2 is the Eddington accretion ratio and L Edd is the Eddington luminosity) (e.g., [18][19][20][21]). The radio emission from radio-loud LLAGNs at frequencies below ∼10 2 GHz is believed to be associated with the non-thermal synchrotron contributed by their jet. At different observational frequencies, the surface where the jet non-thermal synchrotron emission transits from optically thin (τ < 1, where τ is the optical depth) to optically thick (τ > 1) is responsible for the observed "core shift" (e.g., [22,23]); the shift in positions of the unresolved optically thick core (τ ∼ 1) [24].
The details of jet structure and physics can be constrained by the combination of the observed spectra and morphologies. The latter requires enough resolution for the Notes. * the size of black shadow is estimated by 10 R g (e.g., [49]). † the accretion rates are estimated by fitting the spectra of the sources with the RIAF model spectra [43], based on the semi-analytical RIAF model [18].

GRMHD Paradigm of Jet Formation for LLAGN
The accretion flows around LLAGNs are believed to be of the RIAF type. For a RIAF, the ion temperature (T i ∼ 10 12 K) is much higher than electron temperatures (T e ∼ 10 10−11 K) due to the inefficient Coulomb interaction between the ions and electrons at low accretion rates (e.g., [21,50,51]). In turn, the heat stored in the accretion flow results in a puffed-up, thick geometry of the flow.
With the rapid progress of GRMHD simulations since ∼2000 (e.g., [52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67]), details of a geometrically thick, negligible radiation feedback flow environment (such as LLAGN) have been revealed. An example GRMHD simulation result is shown in Figure 1. While large-scale magnetic field lines threaded onto the accretion flow are absent [53], large-scale magnetic field threaded onto the black hole event horizon is naturally developed in the funnel region along with the black hole's rotational axis. In such a funnel region, GRMHD extraction of black hole rotational energy along the large-scale magnetic fields can take place, and the jet is launched at the expense of black hole energy. In fact, relativistic jets powered by the rotation of a black hole via a large-scale magnetic field attached to the event horizon are common scenes in GRMHD simulations (e.g., [68,69]).

Figure 1.
Examples of magnetic field configuration from a 2D GRMHD simulation for accretion flow with negligible radiation feedback, with small (left panel) and large (right panel) fields of view. The black hole rotational axis is aligned with that of the accretion flow. Black hole-threaded, ordered, large-scale magnetic fields (indicated as white lines) are developed along the black hole rotational axis, enclosed by the contours of σ ≡ b 2 /ρc 2 = 1 (purple lines). Magnetic field configurations in the region where σ < 1 are not shown. The contours for u r = 0 (a radial component of the fluid's four velocities) are shown as cyan lines. Relativistic jets are launched from the funnel region at the expense of the rotational energy of black holes. The color shows the density on a logarithmic scale, with blue and red representing small and large densities, respectively.
The stationary GRMHD extraction process under ideal MHD conditions [70] can be viewed as an extension of the electromagnetic extraction process under force-free conditions, known as the Blandford-Znajek process [71]. The electromagnetic part is responsible for extracting the rotational energy from the rotating black hole. Therefore, efficient extraction is performed by a Poynting-dominated GRMHD flow (see also Section 3.1 for a discussion of the fluid and electromagnetic part of the GRMHD flow).
An important feature of the above-mentioned GRMHD paradigm for jet formation is the existence of the boundary for GRMHD flows streaming inward or outward along the same field line as the result of the competition of magneto-centrifugal force and gravity force acting on the plasma [70]. Such a boundary, usually called the stagnation surface or separation surface is shown as the cyan lines in Figure 1. In the inflow region, the GRMHD extraction process mentioned above takes place; in the outflow region, the flow is accelerated as a result of the conversion from electromagnetic energy to fluid energy of the GRMHD flow (e.g., [58,72]). The maximum possible jet velocity is obtained when all the electromagnetic energy converts to the fluid's kinetic energy (e.g., [73]). For both the inflow and outflow regions, starting from negligible poloidal velocities, the accelerated GRMHD flow will pass several characteristic surfaces where the poloidal fluid velocities equal the slowmagnetosonic, Alfvén, and fast-magnetosonic speed [58,70,[74][75][76].
The launching at the expense of black hole energy implies the relationship between jet power and the efficiency of extraction, which depends on black hole spin and magnetic field strength. GRMHD simulation has demonstrated that the jet power can be larger than the accretion power with strong enough magnetic fields [66]. As the magnetic field is supported by the accretion flow environment, there exists an upper limit for possible magnetic field strength, over which the magnetic field cannot be confined by the accretion flow. Accretion flow close to such limit is conventionally referred to as MAD (Magnetically Arrested Disk) [77][78][79][80] and otherwise referred to as SANE (Standard And Normal Evolution) (e.g., [81,82]). Detailed GRMHD simulation and image modelings have been extensively studied for the two target sources of EHT observation, M87 and Sgr A * (e.g., [28,30,35,[83][84][85][86][87][88][89][90][91]). In our following discussions, we do not further distinguish the jet emission features between SANE and MAD because the emission features are related to both the jet structure and assumed emission details, as will be discussed later.

Constructing Jet Model
In this section, we present an overview of constructing global jet structures and dynamics and a collection of flow properties along different large-scale magnetic field lines. As large-scale magnetic field lines are developed only in the funnel (see Figure 1), here we consider only large-scale magnetic field lines threaded onto the event horizon as the field configuration of the jet. A schematic magnetic configuration of a jet model is shown in Figure 2.

Semi-Analytical Jet Models
Semi-analytical jet models are usually assumed to be axisymmetric and stationary. MHD flow streaming in a magnetosphere consists of an ordered, large-scale magnetic field; it is convenient to decompose the equation of motion for MHD flow along the transfield direction and parallel-field direction, which results in the so-called Grad-Shafranov equation (GSE) and wind equation (or the Bernoulli equation) (e.g., [92]). The solution to GSE corresponds to the magnetic field configuration, and the solution to the wind equation corresponds to the flow acceleration along the field.
While solving the GSE is a challenging task (e.g., [93][94][95]), a common working proposal to describe the global magnetic field geometry of a black hole jet is simple: applying the approximated solution for a force-free magnetic field [61,65] where 0 ≤ ν ≤ 1.25 controls the geometry (ν = 0: split-monopole-like configuration; ν = 1: parabolic field line). With such a stream function, different poloidal magnetic field lines are determined by the contours of Ψ = constant. Even for an MHD jet, we can expect that the configuration described by the force-free solution, such as Equation (1), remains good given that the flow is Poynting dominated. By measuring the jet width, the observed jet collimation for M87 is consistent with ν = 0.75 [96][97][98]. Before we move onto the MHD jet models, it is illustrative to estimate the jet dynamics with the drift velocity under the force-free assumption (see [99] for the details of the model , and also [61,100] for the detailed analysis of the acceleration properties). Without consideration of the plasma effect, the velocity is determined solely by the assumed angular velocity of the field and the field configuration. In addition, with a prescribed poloidal magnetic field, the toroidal component of the magnetic field can be estimated under forcefree assumption far from the jet [99]: where R is the cylindrical radius and Ω F (Ψ) (a conserved quantity along the field line Ψ) is the angular velocity of the magnetic field. An overview of the magnetic configuration and velocity of the force-free jet model is given in Table 2. Note that B r and B φ change signs above and below the equatorial plane, while B θ retains the same sign (also see Figure 2). Along the field line, the continuity equation requires the number density n as follows where b is the strength of the magnetic field [99].
Notes. † region z > 0 and z > 0 can be referred to in Figure 2.
As the gravity effect is ignored, the velocity in the force-free jet model is always outward. The near horizon features suggested by the GRMHD framework (such as the existence of the stagnation surface and the inflow region) are, therefore, absent. Nevertheless, by applying a distance floor under which the emission near the horizon is excluded, the covariant form of the force-free jet model provides helpful insights on horizon-scale emission features (e.g., [99,[101][102][103]). It is worth noting that the predicted velocity of a force-free jet seems too large compared to the observed values (e.g., [10,98]). By reducing the toroidal magnetic component (to mimic the effect due to mass loading of magnetic fields), further modification of the force-free jet with a reasonable terminal Lorentz factor is considered in [101]. For practical applications to scales much larger than black hole size, the vector form of the force-free jet model can be applied [104,105].
Under the ideal MHD condition, there are four conserved quantities for an axisymmetric, stationary, cold (i.e., the pressure of the flow is ignored) GRMHD flow: the angular velocity of the magnetic field Ω F , the mass loading (particle number flux per unit electromagnetic flux ; see also Equation (6)), the total energy E, and the angular momentum L of the flow. The flow velocity can be obtained by solving the wind equation (along a field line Ψ), which has the form where u 2 p = u r u r + u θ u θ , µ is the relativistic enthalpy, and U is related to the conserved quantities and the assumed background spacetime [70,[106][107][108][109][110][111].
As an interesting feature of the solution to the wind equation, only the correct combination of the four conserved quantities can provide physical flow, which can successfully pass the fast magnetosonic point and reach infinity [70,72,112]). In addition, the toroidal field can not be prescribed as a prior but need to be obtained from the solution of the wind equa-tion. Furthermore, for Poynting-flux-dominated flow, the plasma loading onto a large-scale force-free magnetic field line can perturb the field configuration during its MHD acceleration (e.g., [113,114]. All these have made the construction of a semi-analytical GRMHD jet model challenging. By analysis of the wind equation, sophisticated relationships between the covariant components of toroidal and poloidal magnetic field [74,115] have been found to guarantee physical, trans-fast magnetosonic MHD outflow solutions. An increase in the ratio can lead to a lower terminal velocity [74,75]. Based on the relationships, a semi-analytical GRMHD jet model for collective flow properties have been recently constructed [75,116]. Further away from the central black hole (where gravitational effect is not essential), GRMHD jet solutions are in good agreement with relativistic MHD solutions [75]. An overview of the GRMHD flow properties for both inflow (u r < 0, the radial component of the four-velocity of the flow) and outflow (u r > 0) regions is presented in Table 3. As shown in the table, in the inflow region (u r < 0, the radial component of the fourvelocity of the flow), a successful extraction and positive outward energy flux E r = Eu r > 0 is possible only when E < 0, and, therefore, only for Poynting-flux-dominated GRMHD inflow (E EM E FL , the subscripts EM and FL correspond to the electromagnetic part and the fluid part of the GRMHD flow). A continuous propagation of the positive outward energy flux can be applied to match the inflow and outflow solutions [72]. Required by the continuity equation, written in a way similar to Equation (3), the number density n follows where subscript p represents the poloidal component.

GRMHD Simulation Models
The definition of the jet for GRMHD simulations varies with different studies and motivations. The boundary between the jet and the accretion flow may be defined by the ratio between magnetic energy and fluid rest energy (σ = b 2 /ρc 2 , where b is the field strength and ρ is the density), the Bernoulli parameter (B = −hu t , where h is the relativistic enthalpy and u t is the t-component of the four velocities of the fluid), or even the angle (see [117] for comparisons of some of these definitions). Among the above possible choices, here we refer to the σ ≥ 1 region as the jet region in GRMHD numerical simulations because the large-scale field lines are enclosed within the region σ > 1 (see also Figure 1). Note that the σ = 1 contours attach to the event horizon. In comparison, in the region where σ < 1, the magnetic field contours reveals chaotic magnetic configurations in the accretion flow, as a result of magneoto-rotational instability (MRI) [118].
Complementarily to the semi-analytical approach, the numerical simulation may suffer from the numerical dissipate process. To prevent the density from being too low in the computational domain (usually within the domain r < 10 R g and σ > 1) and from crashing the numerical simulation, an artificial density threshold is applied (e.g., [58,119]). Table 3. Poynting-flux-dominated GRMHD flow properties along large-scale magnetic field threading black hole † . The table is modified from [72].
Notes. † In the table, all subscript FL and EMs correspond to the fluid part and the electromagnetic part of the GRMHD flow, respectively. For a Poynting-flux-dominated flow, E EM > E FL . Region z > 0 and z > 0 can be referred to in Figure 2.
In principle, GRMHD simulations can provide sophisticated modeling, including both jet and the circum-jet environment, as well as their dynamical features. Large-scale simulations with a computation domain, e.g., 10 4 R g , can be relatively expensive (especially for 3D simulations) as a longer computational time is required to reach a stationary and meaningful physical state of the system.

General Features
For GRMHD flow along large-scale fields, there are two light surfaces respectively in the inflow and outflow region. Here we are interested in the light surface for the outflow, which provides an important reference for the structure and dynamics of the jet. In flat spacetime, the outer light surface has a cylinder profile (if Ω F is constant across field lines), with the cylindrical radius R L = c/Ω F . Beyond the light cylinder, the poloidal flow velocity must be larger than the toroidal velocity, as shown in Figure 3. For a Poyntingflux-dominated MHD flow, the Alfvén surface almost coincides with the light cylinder, implying that the toroidal magnetic field becomes larger than the poloidal field outside the light cylinder.
Noting that the location of the light cylinder is related to Ω F , the jet properties as a function of black hole spin can be understood by how Ω F = αΩ H is associated with the black hole angular velocity Ω H , where α 0.5 is usually considered (e.g., [58,71]). The horizon scale properties, such as the black hole spin, are then related to the large-scale jet properties through the jet boundary and the angular velocity of the large-scale magnetic field lines. For semi-analytical models, the jet boundary can be defined by the last largescale magnetic field line that is attached to the event horizon at the equatorial plane. For GRMHD simulation models, the jet boundary can instead be defined by the σ = 1 contour (see Figure 1).

Figure 3.
Ratio between the toroidal velocity v φ = u φ /u t and radial velocity v r = u r /u t of a force-free jet, as computed semi-analytically from the 4-velocity u α of the covariant force-free jet model with black hole spin a = 0.9. The field line geometry is ν = 1 (Equation (1)), and the angular velocity of all the field lines is half of the angular velocity of the hole. Located within the light cylinders (whose distance to the black hole rotational axis can be approximated by the location where R L = c/Ω F , as shown by the white vertical dashed lines), the toroidal velocity-dominated regions (v φ /v r > 1) are shown in red. Representative field lines are shown with black lines. The boundary of the jet region is defined by the last field line, which is attached to the event horizon at the equatorial plane. Note that the radial velocity is positive in all regions as the gravity effect is not included in the model.

Modeling Jet Emission
Synchrotron radiation is responsible for the observed radio emission from relativistic jets. The unpolarized and polarized emissivity and absorption coefficients of synchrotron radiation from different electron energy distribution have been extensively studied (e.g., [120][121][122]). The coefficients, in general, depend on the electron number density, the local magnetic field strength, and parameters that characterize the energy distribution. For example, for an ensemble of electrons, relativistic Maxwellian energy distribution is associated with the electron temperate, and power-law energy distribution is associated with the power-law index and minimum and maximum particle electrons. A hybrid distribution (thermal plus power-law) is of interest to reflect a more realistic energy distribution. In a hybrid case, the source function S ν for the radiative transfer is constructed by the sum of the contributed emissivity coefficients j ν and absorption coefficients α ν (e.g., [123]): Recently, the kappa distribution, which smoothly connects the thermal core to a power-law tail, has also been considered (e.g., [88][89][90]124]).
Due to the relativistic speed of the jet and the strong gravity in the vicinity of the black hole, the energy shift between the comoving frame and the observer's frame requires additional care in the radiative transfer computation. The covariant form for the energy shift (e.g., [125,126]) is where g αβ is the background spacetime metric, p α is the four-momentum of the photon, and u α is the four-velocity of the fluid, u α | local , or the distant observer, u α | ∞ = (1, 0, 0, 0). The background spacetime would affect the g αβ term in the above equation and also the geodesics. A number of general relativistic radiative transfer (GRRT) tools (e.g., see [127] and the references therein) have been developed to take care of the radiative transfer computation in curved spacetime by solving its Lorentz invariant form (e.g., [126]). When the scattering effect is not crucial, the observed flux of the received ray can be integrated along the geodesics backward in time. In flat spacetime (g αβ = η µν , the Minkowski metric), the photon geodesic is simply a straight line, and Equation (7) reduces to the familiar relativistic Doppler effect characterized by the Doppler factor D (e.g., [24]): where Γ is the Lorentz factor of the jet, β = v/c, and θ is the angle between the jet and the observer. The origin of non-thermal electrons responsible for the jet emission remains a subject of considerable debate. There are different approaches to model the jet emission, with different assumptions of electron energy distribution and spatial distribution. In fact, the injection of non-thermal electrons can rely on the microscope process and is outside the scope of the simulated parameters themselves. If it is the case, we may treat the nonthermal electron as a free parameter. Following such a philosophy, the electron number density does not need to satisfy the relationship derived from the continuity equation, Equations (3) and (6). A possible approach is to link the internal energy of the non-thermal electrons to the magnetic field energy with a fraction η and a possible modification function F [117,128] (see also [129] for further considerations) Interestingly, if the brightness temperature of the observed source is comparable to or smaller than the theoretical possible electron temperature, thermal synchrotron is also capable of being responsible for the observed emission (e.g., [28,35]). From a modeling point of view, thermal synchrotron emissions may originate from regions outside the funnel where large-scale magnetic field configuration appears. A funnel wall can be defined by the region between the funnel and the corona of the flow [59]. In one-fluid GRMHD simulations, only the species (i.e., ions) that dominate the dynamics are simulated. Therefore, the ratio R between T i and T e can be treated as a free parameter: While R is related to microscopic electron thermodynamics, the phenomenological relationship of R can be constructed with σ and plasma beta β P = P gas /P mag , defined by the ratio between the gas pressure and magnetic pressure. The electron thermodynamics are included in more sophisticated GRMHD simulations (see Section 5).
In the following, we consider four illustrative models. The models are constructed by combining the above considerations, and their properties are summarized in Table 4. As the synchrotron radiation becomes more optically thin at higher frequencies, low frequency observations are capable of observing the downstream, extended jet emission. We first demonstrate the modeling of sub-parsec scale jet emissions at relatively lower frequencies (43 and 86 GHz), then milliparsec jet emissions at a higher frequency (230 GHz), assuming that the mass of the central black hole is 5 × 10 9 M and the distance to the black hole is 10 Mpc (1 R g ∼ 4.9 µas for reference). With a black hole spin of a = 0.9, the viewing angle to the jet axis is assumed to be 135 • , and, therefore, the black hole spin vector is pointing away from the observer. The black hole mass, distance, and viewing angle are similar to the parameters for M87.
For the semi-analytical model, the covariant force-free jet model is applied. To mimic the funnel shape of the GMRHD jet model, the jet geometry ν = 1 is applied in Figure 3. The dynamics of the force-free model can be referred to in Figure 3. For GRMHD jet models, the numerical simulation data are from a 2D GRMHD simulation (as shown in Figure 1) performed by the public GRMHD code HARM [130,131], with initial and boundary conditions similar to the setup in [98] but with a larger magnetic flux reservoir. The post-processing for radiative transfer is performed by the public GRRT code Odyssey [132].  Notes. * M BH = 5 × 10 9 M , a = 0.9, D = 10 Mpc (1 R g ∼ 4.9 µas), and viewing angle i = 135 • . ignore all emissions from jet regions below the distance floor z = r sin θ < 10R g . ‡ only from region σ = b 2 /ρc 2 > 1 in the simulation data. † only from region σ < 1 in the simulation data.

Modeling Black Hole Jets at Sub-Parsec Scales
The 43 and 86 GHz model images for the GRMHD jet model A (non-thermal synchrotron from the region σ > 1) and model B (thermal synchrotron from the region σ < 1) are shown in Figure 4 on a logarithmic scale. The projected black hole spin axis is pointing to the left, and the forward jet is pointing to the right. For both of the models, the electron number density is normalized so that the total flux in the same field of view (500R g × 500R g ) is ∼0.4 Jy at 43 GHz. For all jet model images, the imprints of jet dynamics are shown: the forward jet is brighter than the counter jet, and the incoming (bottom) side of the jet is brighter than the receding (top) side. For each model, the 86 GHz image is less extended compared to the 43 GHz image, as the jet becomes more transparent. Such a feature is associated with the observed shift in the position of the unresolved VLBI core (i.e., core shift) of jetted radio sources.
With those common features, the resulting jet model images are related to the emission details considered in the radiative transfer. Non-thermal synchrotron emissions from the funnel region (σ > 1) are computed for GRMHD jet model A. Adopting η = 0.1 and F = 1 in Equation (9) (same assumption in [117]), the number density of non-thermal electrons can be determined by [117]: where p = 3.5 is the power law index of the electron energy distribution between the low-energy cutoff (γ min = 50) and the high-energy cutoff. Under the assumption, the electron number density is solely determined by the magnetic field distribution b of the jet (and independent of the number density given from the simulation data). GRMHD jet model B shows thermal synchrotron emission outside the jet funnel (regions where σ < 1). The number density of thermal electrons is considered by scaling its numerical value to fit the total target flux (∼ 0.4 Jy at 43 GHz), with the ratio R in Equation (10) following the physically motivated relationship with β P [133], with R ∼ R high = 80 at high β P regions (preferentially in the main flow body near the equatorial plane) and with R ∼ R low = 1 at low β P regions (preferentially near the jet funnel region). The limb-brightening feature observed for M87 (e.g., [9] can be reproduced with this approach by enhancing the R high value [133]. Other post-processing options for thermal synchrotron jet emissions have been applied in previous studies. For example, one can assume a constant T e (isothermal jet) [83,84], or constant R (constant ratio jet) [84] in the region where β P < β cut ∼ 0.2, with or without an excised σ < σ cut region. Applying constant R for all (both jet and accretion) regions has also been considered in [62,85]. In addition, hybrid synchrotron emissions with a kappa distribution (kappa jet) are also possible considerations [88][89][90]124]. As all the emissions within σ < 1 are included in GRMHD jet model B, the image morphology near the jet base (the regions near the white circles in the plots) has a much larger size compared to that of GRMHD jet model A. These emissions are mainly contributed by the accretion flow. Another striking jet feature between GRMHD jet models shown in Figure 4 is the relatively brighter stripe near the jet "spine" in the GRMHD jet model A images. Such a feature is also shown in the model images of the semi-analytical force-free jet model, as discussed below.
In Figures 5 and 6, a floor condition is applied to both GRMHD jet model A and the force-free jet model. The emission below a certain floor height z = 10 R g to the black hole is excluded. The purpose is twofold. First, to avoid the emission feature from the unphysical flow dynamics in the force-free jet model. Second, to mimic the case if nonthermal electrons are injected at a certain height. The non-thermal electron number density of the force-free model (with the floor) is assigned to fit the total flux of GRMHD jet model A (with the floor) at 43 GHz. In Figures 5 and 6, the relatively brighter stripe near the "spine" in the jet is shown in both models. This is due to the jet dynamics: across the jet (vertical direction of the figures), from the jet spine to the jet boundary, the flow velocity transits from toroidal-dominated to poloidal-dominated (see also Figure 3). In between, the flow direction would be swept through the observer's line of sight, resulting in a strong relativistic Doppler beaming effect (θ ∼ 0 in Equation (8)). A similar feature has also been found in [105].

Modeling Black Hole Jets Close to the Jet Base
The 230 GHz model images for GRMHD jet models A and B are shown in Figure 7, with the same parameters applied to Figure 4. For our setup, at 230 GHz, the black hole shadow can be observed. As can be seen in the images, the horizon scale image is sensitive to the ratio between the contribution of accretion flow and jet. For GRMHD jet model A, the black hole shadow (the dashed circle in the plot) is enclosed by the jet emission. Surrounding the black hole shadow, the ring-like structure with a brighter side at the bottom is actually mainly contributed by the emission of the counter jet (see also [117]). For GRMHD jet model B, in comparison, the horizon scale image further reveals the motion and emission of the accretion flow, resulting in a larger size of the emission region. For GRMHD jet model B, the black hole shadow is partially blocked by the accretion flow. The imprint of the funnel wall jet can also be seen outside the bright accretion flow emission (near the middle, at the bottom of the figure).
The models with floors are shown in Figure 8. The image morphology is similar for GRMHD jet model A with a floor and for the force-free jet model. As the injection location moves further away from the black hole, the black hole shadow image becomes less obvious. There appears a compact emission region at the incoming side of the jets, with its emission centroid being well outside the location of the black hole shadow. Such features have also been shown in [101].
As demonstrated, the horizon scale images can be largely affected by the uncertainties of emission regions and types of synchrotron radiation. Nevertheless, the background knowledge of jet dynamics and radiative transfer details can help to interpret the environment around the jet base.

Discussion and Outlook
The aim of the review is to summarize and discuss the current progress in modeling the rich features of LLAGN jets at VLBI scales. The modeling of jet images can be performed with jet models and emission details. Jet models can be constructed by semi-analytical approaches or from numerical simulations. The latter is a powerful tool for modeling the black hole accretion jet system, including its time-dependent features. The simulations with a large enough computational domain ( 10 3 R g ) can provide enough information to model jet properties across different spatial scales and frequencies and compare them with observations, e.g., spectra and core shift measurements (e.g., [76,89]), directly. The former, alternatively, provides a flexible and heuristic way of exploring jet physics and properties, given that the key physics can be properly included in the model. The emission details add different complications to the energy distribution of electrons and the emission region from the black hole accretion jet system. To explore the combined effects of the above considerations on the resulting images, jet image features at different frequencies, and different spatial scales are demonstrated with illustrative models.
As demonstrated by our illustrative models, exploring jet images down to the innermost region of the jet is important for several reasons. First, the jet image is jointly determined by the jet dynamics and mass loading. As the jet dynamics are associated with the central engine, the jet image is a possible diagnostic for the property of the central black hole and the required origin of energetic electrons [99,104,105]. Second, the ratio of the contributions between the jet and accretion flow near the sub-mm bump of LLAGN cannot be solely determined by the spectra (see, e.g., [117]) for examples). Based on the background knowledge of the emission features of jet and accretion, the observed morphology of the jetted LLAGN can provide valuable constraints to their relative importance [99,117].
Among the uncertainties in VLBI jet image modeling, electron heating microphysics or injection of non-thermal electrons are important directions to be explored [134]. Interestingly, these may link to possible correlations with other observational frequencies, e.g., X-rays [135]. Is the injection process stationary or intermittent? After injection, what are the effecd of the subsequent cooling process? Does the stagnation surface play any special role in the injection [90,[136][137][138]? All these details are expected to leave observable imprints on the resulting jet image. More advanced, state-of-the-art simulations, including electron thermodynamics, radiation, or beyond one-fluid approaches, have brought more sophisticated investigations into (part of) the above questions [86,87,91,[139][140][141][142][143]. It is also possible to construct dynamical features onto the stationary semi-analytical jet model [103].
Although not discussed in the paper, within the VLBI scales of our interest, there are also several important topics for jet modeling. For example, the jet polarization due to circum-jet materials [128], the jets from tilted disks [144][145][146][147][148], jet composition [149][150][151], and the possible mass-loading due to mixing between the jet and surroundings [76]. More discussion of GRMHD simulation modeling of jet formation can be seen in [152]. The magnetic configuration (Figure 2) may also be revealed by polarization observations at horizon scales [153].

Acknowledgments:
We thank Wen-Ping Lo for helping with the preparation of Table 1. We thank the anonymous referees for helpful comments and suggestions.

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

Abbreviations
The following abbreviations are used in this manuscript: