© The Authors 2024

1. Introduction

Relativistic jets have been established as a prominent feature of some active galactic nuclei (AGN). Even though prominent jets are present for M87 (Owen et al. 1989; de Gasperin et al. 2012) and Centaurus A (Israel 1998; Hardcastle et al. 2003; Janssen et al. 2021), it is not a common feature of AGN in general (Hada 2020). The compact radio source Sagittarius A* (Sgr A*) at the center of the Milky Way does not boast a direct jet detection, even though claims of jet-like features spanning several parsecs have been made (Li et al. 2013). The exact origin of Sgr A^*’s emission signature is still highly debated (as described in, e.g., Issaoun et al. 2019; EHTC 2022a), but we can infer the presence of a (compact or unresolved) jet from the flat or inverted radio spectrum that is difficult to create without an extended emission structure.

Supermassive black holes (SMBHs) have been pivotal in understanding the origin and dynamics of jets (Blandford et al. 2019). Collimated outflows from SMBH systems are expected as a means of angular momentum and energy dissipation either via the Blandford-Znajek (BZ; Blandford & Znajek 1977) or Blandford-Payne (BP; Blandford & Payne 1982) processes. While the former invokes the rotation of the central BH as the main driver, the latter relies on the rotational energy of the disk itself to produce outflows or winds. If sufficient poloidal magnetic field is present and the BH has a non-negligible spin, then the BZ mechanism is assumed to be dominant (e.g., Tchekhovskoy et al. 2012).

Brinkerink et al. (2015, 2021) (hereafter BR15; BR21) observed flares at various radio frequencies (ν = 19 − 47 GHz) and quantified the differences in the arrival times of light curve features as a function of frequency with a cross-correlation analysis. The acquired “time lags” can be used to infer properties of the emitting plasma, mostly with regard to the flow velocity. Previously, time lags have been interpreted with simple in- or outflow models, either via a jet-dominated (Falcke et al. 2009) or an expanding cloud model (van der Laan 1966; Yusef-Zadeh et al. 2008). Until now, however, no rigorous study based on general relativistic magnetohydrodynamical (GRMHD) simulations of BH accretion had been undertaken in 3D (a preliminary work was done in 2D by Okuda et al. 2023). Such GRMHD simulations describe the global accretion dynamics around the central compact object and are known to self-consistently bring about a coupled disk-jet system with the associated in- and outflows.

In this paper, we investigate whether GRMHD models of Sgr A^* can self-consistently capture the time lags, as seen by BR15. Our modeling effort starts from a new set of GRMHD simulations with various black hole spin parameters that are post-processed to acquire synchrotron emission images and light curves. Due to its low Eddington luminosity, accretion onto Sgr A^* is best described by an advection-dominated and radiatively inefficient accretion flow (ADAF and RIAF; Narayan & Yi 1994; Narayan et al. 1995; Abramowicz et al. 1995). Coupled jet-ADAF models were subsequently developed to better recover the spectra of Sgr A^*, and low-luminosity AGN in general (Yuan et al. 2002). Dynamically evolved GRMHD models, which naturally display all previously mentioned characteristics, are typically split in two subclasses: SANE (from Standard And Normal Evolution; e.g., De Villiers et al. 2003) or MAD (from Magnetically Arrested Disk; e.g., Igumenshchev et al. 2003; Narayan et al. 2003). While the former is relatively weakly magnetized and turbulence-driven, the latter has a larger, more highly magnetized disk with a coherent large-scale magnetic field.

For this work, we limit ourselves to the investigation of MAD models since they are currently favored for Sgr A^* (EHTC 2022a,b). A MAD model characteristic is that they display magnetic flux eruptions, associated with a critical accumulation of magnetic field, that produce under-dense and over-magnetized regions in the accretion disk called flux tubes (McKinney et al. 2012; Dexter et al. 2020; Porth et al. 2021). These flux tubes are the end product of magnetic reconnection between the upper and lower hemispheres that culminates in a vertical, poloidal field structure (Ripperda et al. 2020, 2022; Davelaar et al. 2023). Flux tubes turn out to create clear synchrotron emission features (also shown in Najafi-Ziyazi et al. 2024; Davelaar et al. 2023) and are prominently represented in the best-corresponding simulation time lag windows at the evaluated radio frequencies (ν = 19 − 47 GHz). Moreover, we find evidence for interaction between the flux tube exhaust and jet-disk boundary, which is a prominent source of variability in the light curves.

The paper is structured as follows. The methods include a brief overview of our GRMHD simulation configuration, the post-processing procedure, and how one calculates the time lags. This is all outlined in Sect. 2. The results and their interpretation are presented in Sect. 3. The discussion and conclusion can be found in Sect. 4.

2. Methods

In the following sections, we will outline the methods employed for our GRMHD simulations (Sect. 2.1), the radiative transfer procedure (Sect. 2.2), the cross-correlation procedure and observational goodness-of-fit estimation (Sect. 2.3), and an overview of our slow light implementation within the radiative transfer method (Sect. 2.5).

2.1. General relativistic magnetohydrodynamics

We simulate the accretion flow surrounding a rotating Kerr black hole with the Black Hole Accretion Code (BHAC, Porth et al. 2017; Olivares et al. 2019), which solves the MHD equations in a general-relativistic context. These equations are defined as:

$$\begin{array}{cc}& {\mathrm{\nabla }}_{\mu }(\rho u^{\mu })=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\mathrm{\nabla }}_{\mu }{T}^{\mu \nu }=0,\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\mathrm{\nabla }}_{\mu }{}^{\star }F{}^{\mu \nu }=0,\hfill \end{array}$$(3)

where ∇_μ denotes the covariant derivative, ρ the rest-mass density, u^μ the fluid four-velocity, T^μν the energy-momentum tensor (containing both ideal fluid and electromagnetic components), and ^⋆F^μν the (Hodge) dual of the Faraday tensor. We assume an ideal MHD description, which implies that the electric field is acquired via the “frozen-in” condition (i.e., E = −v × B).

Our 3D simulations investigate three black hole spins, a_* ∈ { − 0.9375, 0.0, 0.9375}, which will be referred to as MADM, MAD0, and MADP. The simulations were performed using Modified Kerr-Schild coordinates (MKS, McKinney et al. 2012), which are logarithmically spaced in the radial direction and concentrates resolution in the equatorial plane (with h = 0.25 and R = 0, see Eqs. (49) and (50) in Porth et al. 2017). Nevertheless, we will list all following quantities in Kerr-Schild (KS) coordinates (t, r, θ, ϕ).

One of BHAC’s main features is the inherent adaptive-mesh-refinement (AMR) construction of the underlying grid. Thus, dynamical refine- and derefinement is possible based on user-defined criteria. This allows for a more efficient resolution allocation and, therefore, increases computational performance. While the resolution is typically concentrated in the disk region to ensure a well-resolved magneto-rotational instability (MRI), our simulations have an additional static high refinement block on the jet-disk boundary. As is shown explicitly in Appendix C, these blocks are at the highest AMR level corresponding to the resolution level that is listed in Table 1, alongside other simulation parameters. Due to the logarithmic radial spacing of the grid, this is still at a lower effective resolution level than close to be BH. Another feature of our simulation suite is the high output cadence of 1r_g/c, which is a factor of 5–10 larger than comparable works. This makes it well-suited to conduct variability and “slow light” studies (see both Sects. 2.5 and 3.4).

The simulations are initialized from a hydrostatic equilibrium solution (FM; Fishbone & Moncrief 1976) with the main user-defined parameters and naming conventions being introduced in Table 1. The disk is seeded with a (perturbing) poloidal magnetic field loop that is initialized via the following vector potential (see, e.g., Wong et al. 2021);

$$\begin{array}{c}\hfill A_{\varphi }\propto max(\frac{\rho }{{\rho }_{\max }}{(\frac{r}{r_{\mathrm{in}}}sin\theta )}^3{e}^{-r/400}-0.2,0).\end{array}$$(4)

Here, ρ and r_in denote the density and inner radius of the torus at initialization. The resulting, evolved magnetic field brings about the MAD state which is characterized by magnetic flux eruption that locally and temporarily halt the accretion flow onto the black hole as demonstrated in Fig. 1. These flux eruption come about when the black hole saturates in horizon-penetrating magnetic flux after which the flux eruption generates a vertical magnetic field that creates a barrier and halts the accretion flow. The resulting under-dense, over-magnetized regions are dubbed flux tubes (McKinney et al. 2012; Dexter et al. 2020; Porth et al. 2021), as outlined more explicitly in Sect. 1.

Fig. 1. Flux tubes in the 3D GRMHD simulations. They are visualized with meridial (panels a, b, c) and equatorial (d, e, f) slices and correspond to low-density (ρ), high-magnetization (σ in Fig. C.1) regions that come about after a magnetically driven flux eruption starting from the central BH. The chosen slices correspond to promising time lag windows, as listed in Table D.1. Most jet-related emission is expected from low-density (ρ ∼ 10−2) transition region known as the jet sheath (see panels a, b, c).

2.2. Radiative transfer

To obtain synthetic images, light curves, and spectra, we post-processed our GRMHD simulations with the general-relativistic radiative transfer code RAPTOR (Bronzwaer et al. 2018, 2020), which computes null-geodesics along which the radiative transfer equations are solved assuming thermal synchrotron emission (according to Leung et al. 2011). RAPTOR’s radiative transfer prescription takes into account both the synchrotron emissivity and absorptivity and is shown to be a robust numerical solver (Gold et al. 2020; Prather et al. 2023). It is currently the only code that supports a direct readout of the native non-uniform (AMR) data structure as generated by BHAC (Davelaar et al. 2019). After the null-geodesics are calculated, we infer the electron temperature via a coupling relation as we are only evolving a (proton) one-fluid. We utilize the so-called Rβ description (Mościbrodzka et al. 2016) that allows for a parametrization of the ion-to-electron temperature ratio (T_ratio) as a function of the local plasma-β (ratio of gas to magnetic pressure) from which we then acquire the dimensionless electron temperature (Θ_e) according to:

$$\begin{array}{cc}& T_{}\text{ratio}=\frac{T_{\mathrm{p}}}{T_{\mathrm{e}}}=R_{}\text{low}\frac{1}{1+{\beta }^2}+R_{}\text{high}\frac{{\beta }^2}{1+{\beta }^2},\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\mathrm{\Theta }}_{\mathrm{e}}=\frac{U(\widehat{\gamma }-1)m_{\mathrm{p}}/m_{\mathrm{e}}}{\rho (T_{}\text{ratio}+1)}.\hfill \end{array}$$(6)

Hitherto undefined quantities are the plasma-β = 2p/B², internal energy density (U), proton mass (m_p), and electron mass (m_e). In this work, we exclusively image our models with R_high = 100 and R_low = 1, which tends to produce more extended (i.e., jet-like) emission structures. Generally speaking, a high R_high produces relatively hot jet-sheath electrons and relatively cool disk electrons, while for low R_high, the electron temperature is more uniform. We look at the accretion system under three inclinations angles (i = 10 ° ,30 ° ,50°). Our choices in R_high, as well as the inclinations, are consistent with promising regions of the parameter-space as explored by EHTC (2022a). We evaluate the light curves and images at frequencies of ν ∈ {19, 21, 23, 25, 30, 32, 34, 36, 41, 43, 45, 47} GHz (following BR21). To appropriately cover the large-scale emission features present at these radio-wavelengths, we utilized a spherical polar camera grid, where the pixels are spaced logarithmically in the radial direction (outlined in Davelaar et al. 2018) with a resolution of 256 × 256 pixels, which has been shown to converge with results at higher resolutions (Davelaar et al. 2018), and a radial size of 500r_g. As regions in the inner jet-funnel tend to be dominated by emission from floor violations, it is necessary to exclude these from the radiative transfer procedure. We employed both a geometrical and a standard (cold) magnetization criterion (σ = B²/ρ > 1.0), for which we set the synchrotron emissivity to zero. For the geometrical cut, we assumed a conical (7.5°) cutout up to a radius r = 50r_g; afterward, we utilized a staged cylindrical (with cylindrical radius r_c) cutout (i.e., [i] r_c < 10 for r ≥ 50r_g, [ii] r_c < 25r_g for r ≥ 300r_g).

Lastly, we introduced the relevant scaling relations. Overall, GRMHD simulations are naturally scale-invariant and it is thus possible to scale it to any (RIAF) system. To scale our simulations, we need to define a black hole mass (M_BH), which sets the time and length units, along with a user-defined parameter called the mass unit (ℳ) that sets the overall energy budget of the accretion flow. The total produced emission is also effectively set with ℳ and this parameter is therefore varied until the desired integrated emission level is achieved. We will predominantly be working with a theoretically advantageous unit-set that is based on gravitational radius r_g = GM/c² and its time-unit t_g = r_g/c, which both reduce simply to M in geometrized units G = c = 1. The black hole mass is set to M_BH = 4.1 × 10⁶ M_⊙ in combination with a distance to the central black hole of D = 8.125 kpc (GRAVITY Collaboration 2018). Expressed in cgs units, we find that t_g ≈ 20.2 s and the angular size corresponding to a single r_g is θ_g = 5 μas (see also, e.g., EHTC 2022b).

2.3. Local cross-correlation function and goodness-of-fit estimation

To compute the time lags from our synthetic light curves, we used the local cross-correlation function (Welsh 1999; Max-Moerbeck et al. 2014). The LCCF enables cross-correlation between unevenly sampled timeseries and will be used exclusively to estimate the time lags in this work. Even though the simulated light curves are evenly sampled by construction, we still utilized the LCCF as this is consistent with the methodology outlined in BR15; BR21. We cross-correlated all the light curves with a reference of ν₁ = 19 GHz light curve, which implies that the auto-correlation will serve as the zeroth time lag point. Overall, the choice of reference frequency does not significantly alter the found time lag relations, but some minor variations are perceivable, as is outlined in detail in Appendix E.

To determine the similarity between the simulated (E_sim) and observed time lag (O_obs), we evaluated the χ² statistic for the six listed time lag values (excluding 100 GHz) in BR15 to see which light curve sections are most similar to the observational equivalent. However, as the acquired time lags can be arbitrarily shifted (ξ) in time, we minimized the following function to assess goodness of fit:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i^N}{\left[\frac{O_{\mathrm{obs},i}-E_{\mathrm{sim},i}-\xi }{{\sigma }_{\mathrm{obs},i}}\right]}^2,\end{array}$$(7)

where σ_obs denotes the observational standard deviation and the sum is over the number of frequencies, N = 12. In addition to the goodness-of-fit estimation with the observations directly, we also calculated it with respect to the linear fit from BR15 (with a slope of A = 41 ± 14 cm min⁻¹) to gauge the similarity to the preferred linear relation, as expected from theory (Falcke et al. 2009). We note that we calculate the χ² estimates in a sliding window fashion, which is determined by the starting time and window length. A tabulated overview of the best-fitting windows is listed in Appendix D. These windows are also displayed with horizontal (grey and black) bars in Figs. 4–6.

2.4. Favored linear fit

With the analysis explained here, we have a two-fold objective, namely to assess (i) how linear the acquired time lag relation is and (ii) how well it corresponds to the linear fit acquired in BR15 (with a slope of A = 41 ± 14 cm min⁻¹). As we know from isotropic outflow arguments, a linear relation is expected in the case of a relatively simplistic conical jet structure (Falcke et al. 2009) and this serves as an interesting test to find the deviations from this picture. Thus, more specifically, we will fit a linear relation (y = Ax + B) to all time lags acquired from the simulations and assess the root-mean-square error (RMSE) with a linear fit to this relation to evaluate goodness-of-fit according to:

$$\begin{array}{c}\hfill \mathrm{RMSE}={\displaystyle \sqrt{\frac{1}{N}\sum _{i=1}^N{(T_i-T_{\mathrm{f},i})}^2}.}\end{array}$$(8)

Here, T is the simulation time lag and T_f is the time lag relation prescribed by the fitted linear relation for the (N=)12 frequencies evaluated in this work. Typically, we would find RMSE ≲ 5 for time lag windows that show a reasonably linear relation.

2.5. Slow light implementation

The fast- versus slow-light paradigm has been quite widely investigated for mm-wavelength emission originating close to the BH (see, e.g., Dexter et al. 2010; Bronzwaer et al. 2018; Mościbrodzka et al. 2021). For ray-tracing GRMHD simulations, the standard in the field is to use a “fast light” prescription, which implies that the plasma description is not evolved while the radiative transfer equations are evaluated along the null-geodesics. This effectively renders the speed of light to be infinite. For a more realistic description, denoted as “slow light”, we needs to evolve the plasma and light simultaneously. Overall, we can identify two scenarios where a slow-light prescription could make a significant difference. First, this effect can be important when a strong space-time curvature present, which predominantly indicates emission originating from within the inner 5 − 10r_g. Even though matter moves fast (v ≅ 0.5c) for this scenario, relatively few photons escape from the innermost accretion regions and we typically recover only minor (or even negligable) time lags when only near-horizon emission features are evaluated, which is partly explained by the confined region in which they are created. Second, when matter is moving (relatively) fast over long distances, which we associate with emission from the jet(-sheath). Contrary to all previously mentioned studies at 230 GHz (1.3 mm), our scenario falls within the second category, which is the more complex case as we will outline in the following paragraph.

The complexity of this latter case lies predominantly with the large-scale and diffuse emission structure that presents itself, for example, in the 19 GHz images in this work. So, when the emission structure is compact (e.g., at 230 GHz), then we would only have to read in a relatively low number of GRMHD snapshots (𝒪(50)) to get a consistent picture; typically, the number of required snapshots is comparable to the size of the mission structure. To determine the required number of snapshots, we performed a convergence study that indicates that slow-light images at 19 GHz converge for ∼600 read-in GRMHD snapshots, where we used a 1r_g/c time cadence. This convergence test was conducted by cross-comparing images with an increasing amount of read-in snapshots, where the maximal case was calculated with 1200 read-in GRMHD snapshots. Increasing the snapshot number higher than ∼600 results in minor 𝒪(0.5%) differences in the images and the light curves. For the remainder of the slow light analysis, we therefore used 600 snapshots. Naturally, the convergence number is (partially) determined by when one reaches the limiting optical depth with τ_ν > 1, after which the medium becomes (too) optically thick and the integrated flux density no longer increases. This is the case for all radiative transfer studies, however, it is only for slow light that this limiting surface is dependent on the plasma conditions over a given time period.

In practice, we employed a hybrid slow and fast light approach, where we calibrated the slow-light integration region along the geodesic so that (almost) all emission is described in the slow light window. As mentioned, for this slow light implementation, we read-in 600 GRMHD files (∼350 GB in RAM) that will give the plasma description at different times along the geodesic and it is therefore a memory-intensive procedure. The implementation relies on two dependent parameters, r_slow and N_snap, while the latter denotes how many GRMHD snapshots are read in, the former determines from which distance (to origin) the slow light description is assumed. As the geodesics are integrated backwards, the GRMHD description starts to move back in time until it has reached T_END = T_START − N_snap. All integration before T_START and after T_END assumes fast light. For this approach, we have to calibrate for N_snap and r_slow to find the number after which differences in emission structure are at an acceptably minimal level. In this work, we chose to include a study of the most prominent fast light time window, as we are investigating a scenario where a slow light description could have a significant impact (Sect. 3.4).

3. Results

In the following sections, we discuss the fitted spectra (Sect. 3.1), the qualitative (Sect. 3.2) and quantitative (Sect. 3.3) aspects of the time lag analysis, and the implications of adopting a slow light radiative transfer approach (Sect. 3.4).

3.1. Spectral energy distributions, accretion rates, and jet power estimates

Figure 2 displays the spectral energy distributions (SEDs) for all our ray-traced models. The spectra were fitted with a binary search algorithm (also utilized and outlined in Bronzwaer et al. 2021), which launches the ray-tracing code in an iterative manner and fits the time-averaged (over T ∈ [6300, 15 300]r_g/c) 230 GHz flux density F_ν to 2.5 ± 0.03 Jy. With this relatively simple procedure and an exclusively thermal electron population, we obtained good fits to the quiescent Sgr A^* spectrum at inclinations of i = 10°, 30°, and 50°. Especially in the cases of i = 30° and i = 50°, we obtained SEDs that correspond well to the binned observational data that are gathered from Falcke et al. (1998), Schödel et al. (2011), Bower et al. (2015, 2019), Witzel et al. (2018), von Fellenberg et al. (2018), GRAVITY Collaboration (2020). At i = 10°, we tend to either under- or overproduce the longer wavelength radio emission (i.e., ν ≲ 100 GHz).

Fig. 2. SEDs for three inclinations. From top to bottom, we find i = 10° (panel a), i = 30° (b), and i = 50° (c). Each panel contains the fit corresponding to the models MADP (a* = 0.9375 in dark blue), MAD0 (a* = 0.0 in cyan), and MADM (a* = −0.9375 in red), where the shaded region denotes the one-sigma standard deviation for the 1rg/c spaced light curves per frequency. The origin of the observational points are listed in Sect. 3.1.

For our fits, we find a corresponding accretion rate via Ṁ_cgs = ℳ · ⟨Ṁ_sim⟩/t_g(¹), where the simulation accretion rate (Ṁ_sim) is the averaged over the evaluated time-domain (i.e., T ∈ [6300, 15 300]r_g/c). For inclination i = 30°, we find accretion rates of Ṁ_cgs = 5.69 × 10⁻⁹, 1.52 × 10⁻⁸, 1.62 × 10⁻⁸ M_⊙ yr⁻¹ for MADP, MAD0, and MADM, respectively. As ℳ is calculated for every inclination (and simulation) separately, the accretion rate varies by several factors, but it remains low, as is outlined in Appendix B. Here, we list the i = 30° accretion rates explicitly, as they correspond to the favored models in this work.

As some of the classical time lag interpretation finds its origin in jetted outflow models, it could be informative to gauge the jet power of our GRMHD models. Following the methodology outlined in (detail in Appendix A of) EHTC (2019), we can estimate the jet power (P_jet = ⟨P_jet, sim⟩⋅ℳc²) (²) of our GRMHD simulations (averaged over the evaluated time-domain, T), which were calculated to be P_jet = ∼10³⁹, ∼10³⁷, ∼10³⁹ erg s⁻¹ at a radius of r = 100r_g for MADP, MAD0, and MADM, respectively. We found only minor variations in P_jet based on inclination (i). These estimates are consistent with what was presented previously for Sgr A^* (EHTC 2022a). As no robust jet feature has been observed for Sgr A^* (yet), it is hard to gauge if these values are realistic (see discussion in EHTC 2022a, and references therein). We do note that our interpretation of time lags does not hinge on the presence of a clear or strong jet signature; rather, it focuses on the complex magnetically driven dynamical interactions, as we outline below.

3.2. Qualitative analysis of time lags

Figure 3 displays composite ray-traced images of thermal synchrotron emission at 19 (blue-green), 32 (pink-purple), and 47 GHz (red-orange) for selected windows of the MAD0 case. This allows for a more intuitive interpretation of how emission features move through frequency space, where dark regions indicate coincident emission at all frequencies. As time lags generally move from high to low frequencies, it will start out as a red hue and gradually move to a dark purple hue indicating coincidence of multiple (lower) frequencies as is illustrated in Appendix A. Window 2 (W2, in panels d, e, f) gives a prime example of this behavior, especially if one looks at the flux tube that starts out with a red color and gradually becomes darker as it cools (which is illustrated explicitly in Appendix C).

Fig. 3. Ray-traced images for inclination i = 30°. The images consist of three combined colormaps corresponding to the thermal synchrotron emission at ν = 19, 32, and 47 GHz in shades of green-blue, red-purple, and orange-red, respectively. When emission from all frequencies is present, one creates a black-purple hue in the composite image. All base image maps describe the flux density (Fν) between 0.0 and 0.025 Jy per pixel. The base images are combined with equal weight to acquire the composite images. More details on the procedure as well as animations for all i = 30° cases are provided in Appendix A.

Figures 4–6 display the corresponding light curves and LCCF coefficients of three sections of the light curves (W1, W2, and W3) for the MADP, MAD0, and MADM cases, respectively. These windows are arbitrarily chosen to represent a variety of time lag relations as found for our simulations (detailed statistics are calculated in Sect. 3.3). The classical interpretation of the time lags rests on the frequency-dependent synchrotron emission signature from a (mildly) relativistic, conical jet (see, e.g., Blandford & Königl 1979; Falcke et al. 1993, 2009; Falcke & Biermann 1999), or even an expanding blob of plasma (van der Laan 1966). Even though the jet picture is broadly correct, it is not sufficient to explain the best-fitting time lag windows from our simulations. As these are relatively short observing windows (of up to 400r_g/c ≈ 2.25 h), the variable emission component is arguably more important than the relatively static (jet) emission structure, which is confirmed by our findings. More specifically, we note that a number of windows (in Figs. 4–6) contain good time lag fits and follow closely after a flux eruption which is denoted by the drop in Φ_B (magenta line), as will be explained in more detail in next paragraph(s).

Fig. 4. LCCF coefficients and light curves for the MADP case at i = 30°. For the LCCF calculations (top), we chose one light curve (ν1 = 19 GHz) that we cross-correlated with all other light curves (ν2). The black crosses (“x”) denotes the maximum of the LCCF within a given window. The cyan diamants with errors denote the fit obtained in BR15. Both light curves, ranging from 19 to 47 GHz, and horizon-penetrating magnetic flux ΦB (right-hand axis) are shown in the bottom panel. The three windows W1, W2, and W3 are 400rg/c in length and correspond to the LCCF panels shown in the top. The grey and black horizontal bars denote the best-fitting time lag windows (as prescribed by criteria (iii) and (iv), respectively, in Appendix D).

Fig. 5. LCCF coefficients and light curves for the MAD0 case at i = 30°. The rest of the description of Fig. 4 is also applicable here.

Fig. 6. LCCF coefficients and light curves for the MADM case at i = 30°. The rest of the description of Fig. 4 is also applicable here.

We postulate that the main driver behind the variable component in our MAD simulations is the formation of flux tubes after a magnetic flux eruption. During these flux eruptions, a low-density, high-magnetization region is created after the BH saturates in horizon-penetrating magnetic flux (Φ_B; Tchekhovskoy et al. 2011). When this occurs, a part of the magnetic flux Φ_B (³) is dissipated by means of a magnetic reconnection event that generates a strong vertical magnetic field component that partly halts and pushes back the accretion flow effectively creating the flux tube. After this event, the horizon-penetrating magnetic flux Φ_B (magenta line) drops and seems to strongly correspond with the overall variation of the light curves, especially for the MADP case and, to a lesser extent, for the MADM case. The flux tube is aligned with the jet and orbits the BH at sub-Keplerian velocities before expanding and eventually dissipating into the disk. This can span several orbital periods (∼1500r_g/c for, e.g., W2 in Fig. 5). At its birth, it is predominantly visible in the higher frequency emission (red) before it expands and cools so it will start emitting at lower frequencies also (creating the purple-black hue).

Interestingly, as it matures, the tail-end of the flux tube (furthest away from the BH) produces the strongest emission feature as the flux tube is pushing against the accretion flow which compresses (creating an over-density) and heats the plasma as is clearly shown in Appendix C (also commented on in Najafi-Ziyazi et al. 2024). This emission feature resembles a “hot spot” (Broderick & Loeb 2006; Vos et al. 2022), especially when the flux tube is receding with respect to the observer as one then directly looks at (the back of) the compression region. We note that this picture deviates in interpretation from the standard hot spot model (cf. Vos et al. 2022). There, the point of maximal emission typically occurs when the spot approaches the observer and is maximally (Doppler) beamed. For i = 30 ° /50°, however, we can clearly see a vertical tubular emission structure, which is quite different from the typical spherical hot spot picture. Another channel by which variable emission features are introduced is the non-homogeneous nature of the jet sheath (or wall), which produces the majority of the emission overall. Over-dense outflows will occur sporadically and are sheared and subsequently “smeared out” over the jet sheath cone. Nevertheless, we find that the jet-sheath-related bursts in emission are often closely preceded by a flux eruption event; so, even though they originate is different parts of the simulation domain, they are triggered by the same physical mechanism. While it is typically hard to pinpoint the main emission contributor to well-corresponding time lags, it is striking that flux tubes are prominently featured in the best windows (as can be confirmed with the animations provided in Appendix A together with the horizontal markers in Figs. 4–6).

3.3. Quantitative analysis of time lags

To quantify how well our simulations agree with the findings of BR15; BR21, we calculated the time lags over sections of the synthetic light curves for time windows of 100, 200, 300, and 400r_g/c. These sliding windows are shifted with steps of 10r_g/c until the entire light curve is covered. The acquired time lags are then scored with two types of χ² estimate. Following the method outlined in Sect. 2.3, we calculated ${\chi }_{\text{o}}^2$ directly with the observations. Then, we calculated ${\chi }_{\text{f}}^2$ between the linear fits listed in BR15 and our simulated ones to get an indication of whether we are able to recover this slope. The overview in Table D.1 is listed in Appendix D. We note that the 25 GHz measurement (by BR15) deviates significantly from the otherwise linear trend in the data and is therefore not well represented in the linear fit. We also investigate the preferred linear relations based solely on our simulated data. This will be outlined later in this section and displayed in Fig. 7. Overall, ${\chi }_{\text{o}}^2$ is significantly less constraining than ${\chi }_{\text{f}}^2$ as a result of the large errors on the observational time lags. Nevertheless, the combined constrains (both χ² < 7) are restricting, but still identify a number of passing windows with a preference for the MAD0 case, especially with i = 30 ° /50° and longer correlation windows. The choice for χ² < 7 is arbitrary and quite stringent, but it also clearly highlights windows where the simulation behavior is highly consistent with BR15, as we find a reduced ${\chi }_{\nu }^2$ ≲ 1.2 (for the six observational points in the case of ${\chi }_{\text{o}}^2$ and 12 points for ${\chi }_{\text{f}}^2$ as shown in, e.g., Fig. 4).

Fig. 7. RMSE between the time lag acquired from the simulations and their linear fit. The points are color-coded according to the time lag window length; 100 (magenta), 200 (red), 300 (green), and 400 (blue) rg/c. The color-dashed lines denote the mean slopes of the best-fitting (RMSE < 5) time lags for the vertical lines, while the horizontal lines denote the mean RMSE of the entire population. The grey dashed line, with standard deviation area, corresponds to the linear fit to the VLA-only data from BR15, which is A = 42 ± 14 cm min−1. The contours denote the 25%, 50%, and 75% levels of the normalized kernel density estimations for each of the window length results (in the corresponding colors).

The physical manifestation of the best-fitting time lags corresponds to a flux tube that starts directly in front of the main jet structure and continues its counter-clockwise (as seen from i = 0°) trajectory until it moves behind the jet structure. Even behind the jet, the flux tubes still contributes to the observed flux. Both these points are clearly demonstrated in the provided animations in Appendix A, especially when combined with an evaluation of the light curves displayed in Figs. 4–6. As clearly seen in W2 and W3 of Fig. 3, the flux eruption that creates the flux tube also give an increase in the high frequency (red) emission that will gradually move to lower frequencies in the jet sheath (dark purple). The time lag displays quite an characteristic kink in the low frequency bands to accommodate for the outlier (at 25 GHz). The largest passing window is broadly denoted by W3 for the MAD0 case, as is shown in both Figs. 3 and 5. For completeness, we note that the data listed in BR21 is more complex in nature than what was used for BR15, but is still consistent with the linear relation that we investigate (i.e., ${\chi }_{\text{f}}^2$ in Table D.1).

Figure 7 answers two main questions; how linear are the simulation time lag relations and which slope would best fit them? While all distributions display considerable variance, we note that one finds quite a considerable number of points in the area corresponding to the fit of BR15 in the grey shaded area (for A = 41 ± 14 cm min⁻¹). The mean slopes of the best-fitting time lag distributions (denoted by the black-color dashed vertical lines) are significantly lower than the observational slope. This indicates that the most linear time lags from the simulations favor a more shallow slope than was found for BR15 or that the steeper (better fitting) time lags do not have a (very) linear profile. Therefore, based on the results presented here, we find that the relatively steep slope is a rare occurrence according to our simulations. A last point to consider is that a significant portion of observational time lags do not show a clear relation (BR21). Interestingly, this behavior is also recovered for our idealized simulations as they do regularly not show a clear time lag relation (i.e., there are breaks or multiple competing features in the LCCF coeff.), as is indicated by the black dashed horizontal lines that denotes the mean RMSE of all points. Overall, based on the kernel density contours shown in Fig. 7, we find that most of the simulated time lags show a predominantly linear (RMSE < 5) trend that is typically lower than A = 41 ± 14 cm min⁻¹. The linear trend that is preferred from simulations is A ∈ [7, 23] cm min⁻¹ based on the 400 window length calculations, with A ≈ 20 cm min⁻¹ for MAD0 and lower values for MADP and MADM, which increase with inclination from A ≈ 10 cm min⁻¹ at i = 10° to A ≈ 20 cm min⁻¹ at i = 50°. These listed values are denoted by the vertical lines in Fig. 7.

3.4. Slow light comparison

After one has acquired the slow light curve, it is beneficial to find the optimal offset between both the fast and slow light curves. This offset can then be applied to one of the time series to allow for intuitive comparison. However, finding the offset can be a somewhat non-trivial procedure as the timing between various features differs depending on the specific light description. While there are several options available to find the offset (e.g., via cross-correlation of either the images or the light curves), we have have opted for the minimization of the Euclidean distance between segments of the slow and fast light curves (i.e., 1D equivalent of the down-hill-simplex method; Press et al. 2002). Taking ν = 19 GHz as the reference frequency, we find an offset in time of Δt = 365r_g/c, while for ν = 47 GHz, we find Δt = 402r_g/c. Nevertheless, we would like to note that for the outcome of the (LCCF) time lag calculations, the alignment shift is inconsequential as we are calculating the time lags between the various frequency light curves within their corresponding groups.

Figure 8 displays the fast and slow light curves corresponding to the promising window centered around T = 10 550r_g/c for the MAD0, i = 30° case. What becomes clear is that the slow light description does have a non-negligible effect on resulting time lags. We only calculated the alignment shift for ν = 19 GHz and shifted all slow light curves accordingly. As we expected, timing between various light curve features has changed significantly and results in a relatively worse correspondence with the observational time lags (especially for ${\chi }_{\text{f}}^2$, as listed in Table D.1). This is also reflected in the slow light time lag panels in the top of Fig. 8. Nevertheless, what has not changed is that a relatively large number of time lags passes the ${\chi }_{\text{o}}^2$ criterion. For the ${\chi }_{\text{f}}^2$ constraint, however, we find that almost no time lag passes, which mostly explains the somewhat diminished size in promising window size (grey horizontal bars) when compared to the fast light equivalent in Fig. 5. This indicates that while there is a relatively good consistency when the time lag is compared to the observational data directly, it is not consistent with the proposed linear fit. This is mostly because the slow light time lags display somewhat shallower slopes. We note that W2 in Fig. 8 corresponds almost exactly to W3 shown in Fig. 5. While we find a similar ${\chi }_{\text{o}}^2$ ≈ 7 for both the slow and fast time lags, we find that ${\chi }_{\text{f}}^2$ differs significantly, which highlights the sensitivity of this particular diagnostic. Nevertheless, if we focus on the more agnostic direct observational comparison (${\chi }_{\text{o}}^2$), then we find good agreement overall with the observations.

Fig. 8. Fast and slow light curves for the MAD0 case at i = 30°. The slow light window was chosen to correspond to the most promising time lag window for the fast light curves, as indicated by Table D.1. The color-schemes corresponding to the various frequencies for the fast (solid) and slow (dashed lines) light curves are denoted by “f” and “s”, respectively. The three windows W1, W2, and W3 are 400rg/c in length and correspond to the LCCF panels (calculated for the slow light curves) shown at the top, analogously to Fig. 5 for fast light. The horizon-penetrating magnetic flux (ΦB) is denoted in green (correspond to the right axis) and the horizontal (grey; only for criterion (iii) as outlined in Appendix D) bars denote the best-fitting time lag regions.

To conclude, we comment on why the slow light time lags can differ substantially from their fast light counterpart. When employing the simplest interpretation, there are two main ways in which the light curves change when comparing fast with slow light results, namely: the timing between (light curve or image) features differs or the features themselves are different. While the former point is largely determined by relative emission structure sizes and the associated light travel time differences, the latter point is directly associated with the radiative transfer calculations. The preferred shallower time lag slope does indicate a change in timing and from the light curves, we can already ascertain the more peaky nature of the slow light curves. For the slow light approach, we would expect to find differences in regions where the plasma moves with or against the (integration) direction of the light rays, which either results in increased or reduced emission, respectively. The higher the velocity of plasma, the greater this effect will be. Another clear effect is related to the size of the emission structure, which is smaller for the higher frequencies indicating shorter light travel times through it. Therefore, we would expect the differences at higher frequencies to be smaller than at low frequencies. From this preliminary study, it becomes clear that the adaption of a slow light approach is able to significantly alter the light curves and underlying images at the evaluated wavelengths. Also, it does not necessarily result in better correspondence to the observational time lags (for the evaluated window).

4. Discussion and conclusions

We have demonstrated that flux eruption events and the emergence of flux tubes are consequential for the emission features at radio frequencies (ν = 19 − 47 GHz). Selected windows from our simulations are fully able to reproduce the characteristic observational time lags as presented in BR15; BR21. Interestingly, flux tube emission is prominently featured in the best-fitting time lag windows, but that does not mean that the classical picture with an expanding and cooling jet sheath is no longer applicable. We advocate a more complex picture where the flux eruption, subsequent flux tube creation, and jet sheath emission structure are all affecting and perturbing one another, which is clearly seen in the most promising windows (in Fig. 3). What happens to these features at larger distances still remains an open question. As our simulations display relatively strong outflows extending several hundreds of r_g, we consider the exact origin of the emission in VLBI jet (Lobanov & Zensus 1999; Kim et al. 2023) observations, which could perhaps be explained with a strong outflow (variability) starting at the central black hole rather than the (re)collimation shock picture. However, we also note that our modeling only deals with the inner 500r_g and that these VLBI jet observations span thousands of r_g.

A number of effects that are not within the scope of this work are likely to influence our results. For the sake of improvement, we may consider the addition of non-thermal emission (Özel et al. 2000; Chan et al. 2009; Chael et al. 2017; Davelaar et al. 2018), along with a GRMHD simulation focusing on better resolving the jet sheath (as seen in, e.g., Ripperda et al. 2022) and more elaborate electron heating prescriptions (see, e.g., Howes 2010; Rowan et al. 2017, 2019; Kawazura et al. 2019); alternatively, even a more in-depth Rβ study could be useful. Another possibility would be the inclusion of non-ideal currents in the vicinity of the jet-disk interface with a resistive GRMHD approach (Ripperda et al. 2020; Vos et al. 2024). Additionally, polarized light curves (see, e.g., Mościbrodzka et al. 2021; Najafi-Ziyazi et al. 2024; Davelaar et al. 2023) are likely to shed more light on the preferred magnetic field orientation – if it diverges significantly from the commonly applied singular poloidal loop, as outlined in Sect. 2.1.

In this work, we have also undertaken an exploratory study to establish the sensitivity of our results to the fast-light assumption and found that it does indeed have a significant effect on the recovered time lags (see Sect. 3.4). At face value, the slow light time lags did not exhibit a better correspondence with the observational findings, which highlights the fact that the fast versus slow light paradigm is not intuitively understood, especially for a timing-sensitive study (as is the case here). Nevertheless, we have demonstrated that it is possible to explain the observed time lags originating from Sgr A^* utilizing the modeling techniques outlined in this work.

Lastly, we note that the time lag slope of A = 41 ± 14 cm min⁻¹ (BR15; BR21) is only recovered on relatively rare occasions (see Appendix D) and seems to correspond to a particular configuration of the emission structure. Here, the flux tube starts its orbit on the front-right side of the main emission structure and recedes with respect to the observer over the evaluated time window. Nevertheless, as noted explicitly in BR21, observational time lags tend to span across a range of 20 ≲ A ≲ 40 cm min⁻¹ and a considerable part of the simulation time lag distributions is consistent with this range as shown in Fig. 7. Our analysis favors a strong flux tube emission component to be present to explain the high observational time lag slope of A ≈ 40 cm min⁻¹. Overall, it seems that the a_* = 0 (MAD0) case outperforms the others two spin cases (a_* = 0.9375 and a_* = −0.9375 for MADP and MADM, respectively) and has a more variable time lag signature as shown in the confidence intervals in Fig. 7. This is likely tied to the comparative flux tube to base emission structure strength, where we expect the base emission structure of MAD0 case to be less prominent due to the suppression of the BZ mechanism. In a future study, we aim to investigate the implications of our findings at higher frequencies and including polarization.

Acknowledgments

We thank Monika Mościbrodzka for the insightful discussions and comments that have improved the quality of the manuscript. J.V. acknowledges support from the Dutch Research Council (NWO) supercomputing grant No. 2021.013. J.D. is supported by a Joint Columbia University and Flatiron Institute Postdoctoral Fellowship. Research at the Flatiron Institute is supported by the Simons Foundation. J.D. acknowledge support from NSF AST-2108201. During part of this project, funding for H.O. came from Radboud University Nijmegen through a Virtual Institute of Accretion (VIA) postdoctoral fellowship from the Netherlands Research School for Astronomy (NOVA). H.O. is supported by the Individual CEEC program - 5th edition funded by the Portuguese Foundation for Science and Technology (FCT). This work is supported by the Center for Research and Development in Mathematics and Applications (CIDMA) through the FCT, references UIDB/04106/2020, UIDP/04106/2020 (https://doi.org/10.54499/UIDB/04106/2020 and https://doi.org/10.54499/UIDP/04106/2020). H.O. acknowledges support from the projects PTDC/FIS-AST/3041/2020, CERN/FIS-PAR/0024/2021 and 2022.04560.PTDC. This work has further been supported by the European Union’s Horizon 2020 research and innovation (RISE) programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740 and by the European Horizon Europe staff exchange (SE) programme HORIZON-MSCA-2021-SE-01 Grant No. NewFunFiCO-10108625. Software: python (Oliphant 2007; Van Rossum & Drake 2009), scipy (Virtanen et al. 2020), numpy (Harris et al. 2020), matplotlib (Hunter 2007), RAPTOR (Bronzwaer et al. 2018, 2020), BHAC (Porth et al. 2017; Olivares et al. 2019).

References

Appendix A: Animations

To display the time lag in a visually intuitive manner, we combined the images at multiple frequencies into a single multi-frequency image, from which one can interpret regions of simultaneous emission or when a feature is predominantly present at a single frequency. The colors per frequency are as follows: green-blue for ν = 19 GHz, pink-purple for ν = 32 GHz, orange-red for ν = 47 GHz that are combined into a red to dark-purple and black master color scheme. All images saturate at a flux density (F_ν) per pixel of 0.025, which allows for interpreting which image feature (at which frequency) is most dominant. As we modify the RGBA codes of the base figures to acquire the top figure, there is a potential loss of some physical interpretability, but the visual interpretation of the time lag becomes more intuitive. The interpretation of the strongest emission features is naturally still robust. An example of how the combination of base images occurs is shown in Fig. A.1. As mentioned, all i = 30° simulations have accompanying animations that are available in this online repository⁴.

Fig. A.1. Schematic outlining the components of the multi-frequency graphics as shown in Fig. 3. The panels (a,b,c) display the flux density maps (Fν) at frequencies ν = 19, 32, 47 at T = 9230 of the MADM, i = 30° case, which is the same instance as shown in panels (c) and (f) of Fig. 1. The image’s RGBA codes are combined into the final result in the bottom right (panel d).

Appendix B: Accretion rate for all simulations

Table B.1 lists the acquired accretion rates as Ṁ_cgs = ℳ · ⟨Ṁ_sim⟩/t_g as a function of the time-averaged simulation accretion rate (⟨Ṁ_sim⟩) and the user-defined mass unit, ℳ, as outlined previously in Sects. 2.2 and 3.1. The flux calibration procedure, which we outline in more detail here, is based on incrementally changing ℳ based on the average flux acquired from a 100r_g/c spaced light curve. After updating ℳ, RAPTOR is launched 900 times as the total length of the light curve is 9000r_g/c. The ℳ space is sampled by means of a binary search and typically converges on the desired flux density (F_ν = 2.5 Jy at 230 GHz) in ten cycles or less.

Appendix C: Grid layout of the GRMHD simulations and flux tube temperature profiles

Figure C.1 displays the AMR grid layout of our simulations, with a high refinement block covering the jet-disk boundary, and gives insight into the temperature profile of the flux tubes. The high refinement blocks allow for capturing structural variability in the jet sheath better than is standardly applied. As it is typically advantageous to keep the blocks at lowest level along the jet axis, we opted for the user-defined and intricate refinement scheme. Nevertheless, to sufficiently resolve plasma-driven (i.e., Kelvin-Helmholtz instabilities; Sironi et al. 2021) or even current-driven (i.e., magnetic reconnection; Ripperda et al. 2020; Vos et al. 2024) instabilities one needs much higher resolution levels to resolve these properly. Additionally, as discussed in Sect. 3.2, we find that tail-end of the flux tube is both heated and over-dense which gives it a clear emission feature, while the very center of the tube (having low density) is typically excluded from the GRRT analysis (as σ > 1).

Fig. C.1. Flux tubes in the 3D GRMHD simulations. They visualized with meridial (panels a, b, c) and equatorial (panels d till i) slices. Here, in addition to Fig. 1, we display the (cold) magnetization σ = B2/ρ (panels a till f) with the AMR structure (in grey) and temperature T = p/ρ (panels g till i). Each AMR block contains 8 × 8 × 8 cells. For intuitive comparison with the GRRT images, we display T regions with σ > 1 in magenta as the emissivity is set to zero there. As mentioned, the chosen slices correspond to promising time lag windows and are listed in Table D.1.

Appendix D: Best-fit time windows

Table D.1 lists the time-windows of the simulated time lags that fit best to the relation found in BR15. These windows are identified by means of ${\chi }_{\text{f}}^2$ (related to the linear fit with slope A = 40 ± 14 cm/min) or ${\chi }_{\text{o}}^2$ (direct comparison with the six observational time lags). The main identification criteria (${\chi }_{\text{x}}^2$ < 7, x ∈ {o, f}) are stringent, which ensure that we only identify the very best-fitting time-windows. We note that the reduced chi-squared, ${\chi }_{\nu }^2$, corresponding to the aforementioned criteria are ${\chi }_{\nu ,\text{o}}^2=1.17$ and ${\chi }_{\nu ,\text{f}}^2=0.58$, which clearly outlines that the ${\chi }_{\text{f}}^2$ criterion is the main decider for the best-fitting windows. This also explains the main trends shown in the table, which we briefly discuss now.

Generally, from all evaluated cases, MAD0 is best able to recover the desired time lag, especially for inclination of i = 30° and i = 50°. Overall, column (i) outlines that numerous light curve sections have passable correspondence to the observations, where the i = 10° cases are slightly disfavored (especially for MADP and MADM). Column (ii) is where the clear preference arises for the MAD0 case, as it seems that the other cases are less able to recover the needed slope which is also clearly seen for Fig. 7. Columns (iii) and (iv) further accentuate the previously discussed point, but are nevertheless useful for selecting the very best time lag windows. Nevertheless, we note that the ${\chi }_{\text{o}}^2$ < 7 criterion is leading as it is (most) free from any predefined assumption (i.e., linearity with a certain slope). From this result, we conclude that there are numerous windows that are consistent, but there are indications for a preferred medium inclination (i = 30 ° /50°) and low BH-spin (MAD0).

Appendix E: Different reference frequency

We have exclusively reported our results using the reference frequency of ν₁ = 19 GHz. Even though we typically find very consistent time lags when choosing a different reference frequency, it is still possible that the time lag behavior changes. Here, we show the differences in results when using a reference frequencies of ν₁ = 32 GHz. As mentioned before, only the slope of the relation is important for our analysis, as we would naturally get a different offset when changing the reference frequency. In essence, we will therefore be repeating the analysis outlined in Sects. 2.4 and 3.3 to estimate how much the choice in reference frequency affects the final time lag interpretation.

Based on the comparison of Fig. E.1 (ν₁ = 32 GHz) and Fig. 7 (ν₁ = 19 GHz), we conclude that there are differences, which are best pointed out by the density contours, between both time lag distributions, but overall they are in good agreement. More specifically, the variation in the distribution itself is lower (i.e., density contours are more confined) and the best-fitting linear slopes (for RMSE < 5) are slightly lower than what was found for ν₁ = 19 GHz. Nevertheless, we still find a good number of time lag slopes that lie within the shaded region. In our experience, when the reference light curve is smooth or devoid of small-scale features, it results in a more “tolerant” time lag profile, which often contains competing features. When we correlate the light curves of neighboring frequencies, we are able to acquire the clearest outcomes. Although the reference light of ν₁ = 19 potentially results in a more variable time lag profile, we find that the results and the interpretation are still robust, even when a different reference frequency is chosen.

Fig. E.1. RMSE between the time lags acquired from the simulations and their linear fits for a reference frequency of ν1 = 32 GHz. The points are color-coded according to the time lag window length; 100 (magenta), 200 (red), 300 (green), and 400 (blue) rg/c. The color-dashed lines denote the mean slope of the best-fitting (RMSE < 5) time lags for the vertical lines, while the horizontal lines denote the mean RMSE of the entire population. The grey dashed line, with standard deviation area, corresponds to the linear fit to the VLA-only data from BR15, which is A = 42 ± 14 cm/min. The contours denote the 25%, 50%, and 75% levels of the normalized kernel density estimations for each of the window length results (in the corresponding colors). Figure 7 displays the results for a reference frequency of ν1 = 19 GHz, which was the base assumption throughout this work.

All Tables

All Figures

Fig. 1. Flux tubes in the 3D GRMHD simulations. They are visualized with meridial (panels a, b, c) and equatorial (d, e, f) slices and correspond to low-density (ρ), high-magnetization (σ in Fig. C.1) regions that come about after a magnetically driven flux eruption starting from the central BH. The chosen slices correspond to promising time lag windows, as listed in Table D.1. Most jet-related emission is expected from low-density (ρ ∼ 10−2) transition region known as the jet sheath (see panels a, b, c).

In the text

Fig. 2. SEDs for three inclinations. From top to bottom, we find i = 10° (panel a), i = 30° (b), and i = 50° (c). Each panel contains the fit corresponding to the models MADP (a* = 0.9375 in dark blue), MAD0 (a* = 0.0 in cyan), and MADM (a* = −0.9375 in red), where the shaded region denotes the one-sigma standard deviation for the 1rg/c spaced light curves per frequency. The origin of the observational points are listed in Sect. 3.1.

In the text

Fig. 3. Ray-traced images for inclination i = 30°. The images consist of three combined colormaps corresponding to the thermal synchrotron emission at ν = 19, 32, and 47 GHz in shades of green-blue, red-purple, and orange-red, respectively. When emission from all frequencies is present, one creates a black-purple hue in the composite image. All base image maps describe the flux density (Fν) between 0.0 and 0.025 Jy per pixel. The base images are combined with equal weight to acquire the composite images. More details on the procedure as well as animations for all i = 30° cases are provided in Appendix A.

In the text

Fig. 4. LCCF coefficients and light curves for the MADP case at i = 30°. For the LCCF calculations (top), we chose one light curve (ν1 = 19 GHz) that we cross-correlated with all other light curves (ν2). The black crosses (“x”) denotes the maximum of the LCCF within a given window. The cyan diamants with errors denote the fit obtained in BR15. Both light curves, ranging from 19 to 47 GHz, and horizon-penetrating magnetic flux ΦB (right-hand axis) are shown in the bottom panel. The three windows W1, W2, and W3 are 400rg/c in length and correspond to the LCCF panels shown in the top. The grey and black horizontal bars denote the best-fitting time lag windows (as prescribed by criteria (iii) and (iv), respectively, in Appendix D).

In the text

	Fig. 5. LCCF coefficients and light curves for the MAD0 case at i = 30°. The rest of the description of Fig. 4 is also applicable here.
In the text

	Fig. 6. LCCF coefficients and light curves for the MADM case at i = 30°. The rest of the description of Fig. 4 is also applicable here.
In the text

Fig. 7. RMSE between the time lag acquired from the simulations and their linear fit. The points are color-coded according to the time lag window length; 100 (magenta), 200 (red), 300 (green), and 400 (blue) rg/c. The color-dashed lines denote the mean slopes of the best-fitting (RMSE < 5) time lags for the vertical lines, while the horizontal lines denote the mean RMSE of the entire population. The grey dashed line, with standard deviation area, corresponds to the linear fit to the VLA-only data from BR15, which is A = 42 ± 14 cm min−1. The contours denote the 25%, 50%, and 75% levels of the normalized kernel density estimations for each of the window length results (in the corresponding colors).

In the text

Fig. 8. Fast and slow light curves for the MAD0 case at i = 30°. The slow light window was chosen to correspond to the most promising time lag window for the fast light curves, as indicated by Table D.1. The color-schemes corresponding to the various frequencies for the fast (solid) and slow (dashed lines) light curves are denoted by “f” and “s”, respectively. The three windows W1, W2, and W3 are 400rg/c in length and correspond to the LCCF panels (calculated for the slow light curves) shown at the top, analogously to Fig. 5 for fast light. The horizon-penetrating magnetic flux (ΦB) is denoted in green (correspond to the right axis) and the horizontal (grey; only for criterion (iii) as outlined in Appendix D) bars denote the best-fitting time lag regions.

In the text

	Fig. A.1. Schematic outlining the components of the multi-frequency graphics as shown in Fig. 3. The panels (a,b,c) display the flux density maps (Fν) at frequencies ν = 19, 32, 47 at T = 9230 of the MADM, i = 30° case, which is the same instance as shown in panels (c) and (f) of Fig. 1. The image’s RGBA codes are combined into the final result in the bottom right (panel d).
In the text

Fig. C.1. Flux tubes in the 3D GRMHD simulations. They visualized with meridial (panels a, b, c) and equatorial (panels d till i) slices. Here, in addition to Fig. 1, we display the (cold) magnetization σ = B2/ρ (panels a till f) with the AMR structure (in grey) and temperature T = p/ρ (panels g till i). Each AMR block contains 8 × 8 × 8 cells. For intuitive comparison with the GRRT images, we display T regions with σ > 1 in magenta as the emissivity is set to zero there. As mentioned, the chosen slices correspond to promising time lag windows and are listed in Table D.1.

In the text

Fig. E.1. RMSE between the time lags acquired from the simulations and their linear fits for a reference frequency of ν1 = 32 GHz. The points are color-coded according to the time lag window length; 100 (magenta), 200 (red), 300 (green), and 400 (blue) rg/c. The color-dashed lines denote the mean slope of the best-fitting (RMSE < 5) time lags for the vertical lines, while the horizontal lines denote the mean RMSE of the entire population. The grey dashed line, with standard deviation area, corresponds to the linear fit to the VLA-only data from BR15, which is A = 42 ± 14 cm/min. The contours denote the 25%, 50%, and 75% levels of the normalized kernel density estimations for each of the window length results (in the corresponding colors). Figure 7 displays the results for a reference frequency of ν1 = 19 GHz, which was the base assumption throughout this work.

In the text