© The Authors 2025

1. Introduction

In the center of the Milky Way, at a distance of D ≃ 8.3 kpc, lies Sagittarius A* (Sgr A*), a 4.3 × 10⁶ M_⊙ (GRAVITY Collaboration 2022) supermassive black hole. The accretion disk around it is optically thin, with collisionless, high-temperature, low-density plasma. Its bolometric luminosity is L ∼ 10³⁶ erg s⁻¹, which is in the sub-Eddington range (Genzel et al. 2010). The emission from Sgr A* is systematically monitored over a broad range of frequencies, allowing us to investigate the dynamics of the disk in the vicinity of a supermassive black hole.

Observations from the GRAVITY Collaboration have shown a hot spot rotating around Sgr A* with a period of ∼60 minutes at a distance of a few gravitational radii (hereafter r_g) from the black hole (GRAVITY Collaboration 2018b). Sgr A* is known to exhibit variability in the near-infrared, which manifests itself in the form of several flares over a single day (Ghez et al. 2005b; Do et al. 2019). For the first time, GRAVITY observations revealed that a flare in Sgr A* coincided with a hot spot moving around the central black hole. Several investigations suggest that flares originate in highly magnetized structures that are formed in the innermost region of the black hole and produce synchrotron radiation (Dodds-Eden et al. 2009; Boehle et al. 2016; Ponti et al. 2017; Chatterjee et al. 2021; Scepi et al. 2022).

Analytical models have been applied to study the trajectory of the hot spots in the resulting light curve and show the importance of the effect of gravitational lensing (GRAVITY Collaboration 2005; Younsi & Wu 2015; Bauböck et al. 2020; Ball et al. 2021). As shown by Matsumoto et al. (2020), assuming that its motion is in the equatorial plane, the GRAVITY hot spot follows a circular path at super-Keplerian velocity. However, an outflowing conical sub-Keplerian orbit can fit these observations equally well (Antonopoulou & Nathanail 2024).

The theoretical astrophysics community has increasingly turned to numerical simulations, specifically general relativistic magnetohydrodynamic (GRMHD) simulations, in order to gain a comprehensive understanding of the dynamics involved in accreting black holes. Observations of flares and moving hot spots further emphasize the importance of considering both magnetic fields and fluid dynamics within the extreme gravitational environment near a black hole (Dexter et al. 2020; Porth et al. 2021; Chatterjee et al. 2021; Čemeljić et al. 2022; Scepi et al. 2022; Mellah et al. 2023; Lin & Yuan 2024).

Numerical investigations explore the formation of current sheets in the vicinity of a black hole and the subsequent production of plasmoids and plasmoid chains (Nathanail et al. 2020, 2022b; Ripperda et al. 2022). Current sheets are responsible for particle acceleration and the generation of variable nonthermal radiation. These models need to be supported by numerical studies of turbulence and magnetic reconnection in collisionless plasmas that self-consistently capture the dynamic interplay between particles and fields on kinetic plasma scales. Such investigations can be performed with particle-in-cell (PIC) simulations (Drake et al. 2012; Guo et al. 2014, 2015; Sironi & Spitkovsky 2014; Shay et al. 2014; Dahlin et al. 2014; Li et al. 2015, 2023; Petropoulou & Sironi 2018; Werner et al. 2018; Ball et al. 2018; Comisso & Sironi 2019). In particular, these investigations help us study the acceleration processes in detail by tracking individual particles, as well as characterize the properties of the particle distribution function, for example the shape of the distribution (power law or log-parabolic), the fraction of energy carried by nonthermal particles, the maximum particle energy, and others (e.g., Werner et al. 2018; Ball et al. 2018; Petropoulou et al. 2019; Zhang et al. 2023). Some of the particle properties, like the slope of the power-law distribution, can be related, through empirical expressions, to the local plasma properties, such as magnetization, σ ≔ B²/ρ (where B is the magnetic field strength and ρ is the density) and plasma β ≔ 2P/B² (where P is the gas pressure). These prescriptions can then be combined with GRMHD simulations to anticipate the generation of nonthermal, high-energy electron distributions with properties that depend on the local physical conditions (Chatterjee et al. 2021; Scepi et al. 2022; Aimar et al. 2023; Lin et al. 2023; Mellah et al. 2023). The way to connect GRMHD simulations and PIC results is the following: After the plasma accretion process has reached an inflow equilibrium (Dexter et al. 2020), the thermal radiation is attributed to the hot electrons with a parameterized temperature derived from the simulation ion temperature (Mościbrodzka et al. 2016). To include nonthermal, radiating electrons, formulas from PIC simulations are employed. Nonthermal electron distributions are assumed in each cell, and they are assumed to depend on the magnetization (for σ < 5) and plasma β (Davelaar et al. 2019; Fromm et al. 2022; Cruz-Osorio et al. 2022). This study proposes a novel approach for incorporating nonthermal particles into simulation results: Current sheets are first identified (see also Vos et al. 2024), and then their environment is characterized in terms of averaged quantities of magnetization (σ) and plasma β. Finally, nonthermal particles are sourced from the plasma within the current sheet. In this study we utilized the 2D GRMHD simulations of Nathanail et al. (2020) to investigate the ability of the current sheet and plasmoids to produce the observed nonthermal flaring activity of Sgr A*. We applied a radiation model to produce light curves at 230 GHz (thermal radiation from the disk) and 2.2 μm (nonthermal radiation from current sheets; see also Scepi et al. 2022).

Section 2 describes the methodology of our investigation and is divided into three subsections: we introduce the GRMHD simulation setup (Sect. 2.1), the modeling of the radiation (Sect. 2.2), and the method with which we determined current sheets and their environment from the GRMHD simulation (Sect. 2.3). In Sect. 3 we present our results for the thermal (Sect. 3.1) and nonthermal radiation (Sect. 3.2). Finally, in Sect. 4 we present our conclusions and discuss prospects for future work.

2. Numerical setup

The numerical methods employed for the GRMHD simulation in this study closely mirror those used in previous research (Nathanail et al. 2020). We initialized the simulation with magnetic field configurations designed to generate multiple current sheets as the accretion system evolves. The formation of current sheets occurs in the vicinity of the black hole, and their subsequent reconnection gives rise to the development of multiple plasmoids and plasmoid chains. To assess the radiation characteristics of the current sheets formed and the subsequent generation of plasmoids, we employed a radiation proxy model. This model enables the calculation of thermal and nonthermal radiation at 230 GHz and 2.2 μm, respectively.

2.1. GRMHD simulation

The accretion disk surrounding a black hole can be conceptualized as a hydrodynamic system in the context of curved space-time, where the magnetic field significantly influences its dynamics. A GRMHD code integrates these physical properties into a unified numerical simulation and allows a robust evolution in time, where the accretion and accumulation of magnetic flux produce energetic phenomena. We employed the BHAC code (Porth et al. 2017), which solves the ideal magnetohydrodynamic equations in general relativity, namely

$$\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 }{}^{\ast }{F}^{\mu \nu }=0,\hfill \end{array}$$(3)

where T^μν and ^*​F^μν are the energy momentum tensor and the dual of the Faraday tensor, respectively. Here, we denote with ρ the rest-mass density and with u^μ the fluid four-velocity. Some important characteristics of the code are the following: it implements second-order high-resolution shock-capturing finite-volume methods and adaptive mesh refinement wherever needed. It also implements constrained-transport (Del Zanna et al. 2007) to guarantee a divergence-free magnetic field (Olivares et al. 2019).

Furthermore, the model is axisymmetric (2D), and the coordinates of the initial torus-like plasma distribution are spherical (r × θ × ϕ). We used a logarithmic radial grid and the domain extends out to 2500 r_g. The resolution of the model is 4096 × 2048 × 1 cells. The torus at the initial equilibrium state has constant specific angular momentum l = 4.28 (Fishbone & Moncrief 1976). The inner radius was set to r_in = 6 r_g and the pressure maximum radius to r_max = 12 r_g. All quantities are calculated in geometrized units (G = c = 1) in which the gravitational radius is equal to r_g = M. We considered a Kerr black hole with dimensionless spin a = J/M² = 0.93, where J and M are the angular momentum and mass of the black hole, respectively. The radius of the event horizon is r_h = 1.341 r_g, and there are 29 × θ_res (θ_res = 2048 resolution in θ direction) grid cells inside the horizon.

The initial conditions for the magnetic field inside the torus consist of several loops of interchangeable clockwise–counterclockwise orientation between neighboring loops (see Fig. 1). To generate a magnetic field topology as described above, we used a vector potential, A_ϕ:

$$\begin{array}{c}\hfill A_{\varphi }\propto A\times B,\end{array}$$(4)

Fig. 1. Initial magnetic field (lines) and density (color) configuration in the initial torus. Axis coordinates are calculated in rg. Red and blue denote a clockwise and counterclockwise direction of the magnetic field.

with

$$\begin{array}{c}\hfill A\equiv \max (\rho /{\rho }_{\max }-0.2,0),\end{array}$$(5)

$$\begin{array}{c}\hfill B\equiv cos[(N-1)\theta ]sin[2\pi (r-r_{\mathit{in}})/{\lambda }_r],\end{array}$$(6)

where N is the number of poloidal loops and λ_r their characteristic length scale, and ρ_max is the maximum rest mass density of the disk.

State-of-the-art GRMHD codes have provided significant insights into the role of initialized magnetic fields with nested loops in jet formation (Narayan et al. 2022; Wong et al. 2021; Porth et al. 2019). The strength of the resulting jet is highly dependent on the initial conditions, particularly the size of the initial torus. A larger torus results in greater magnetic flux. For example, a standard and normal evolution (SANE) configuration with a small torus produces a weaker jet, while a magnetically arrested disk (MAD) configuration with a larger torus generates a more powerful jet.

In our case (SANE multi-loop), the initial setup causes the inverse polarity of the magnetic field to reconnect with the first polarity as the loops accrete. This prevents the black hole from accumulating sufficient magnetic flux to sustain a stable jet. We note that our simulations are based on a small torus configuration. For studies involving multi-loop and large torus setups, we refer readers to Jiang et al. (2023). It is important to note that the behavior is highly dependent on the initial conditions. To provide a physical basis for the initial multi-loop magnetic field, we refer to magnetorotational instability (MRI) simulations, which demonstrate a cycle of B_ϕ changing polarity over time. This behavior is precisely what the initial multi-loop magnetic field aims to replicate (Nathanail et al. 2022a).

A key debate within the scientific community is whether the supermassive black hole at the center of our galaxy is capable of producing a stable jet. We have not yet observed a jet from Sgr A*, and it is an open question in which direction our models should go for the interpretation of phenomena near the black hole. As discussed in Nathanail et al. (2022a), images at 43 and 86 GHz will produce an extended stable jet base in the case of the MAD model, whereas in the case of the SANE multi-loop model, the images will possibly have an extended jet-like structure only in the flaring state. We expect that new images of the Sgr A* black hole will further clarify the landscape of pertinent models.

Furthermore, it is important to emphasize that while 2D GRMHD simulations are valuable for exploring the dynamics of accretion flows and jet formation, they inherently impose certain limitations compared to 3D simulations. In 2D, the axisymmetric assumption restricts the development of non-axisymmetric instabilities, such as the MRI turbulence, which plays a crucial role in angular momentum transport and magnetic field amplification. This simplification can lead to an over-simplified representation of the physical processes within the accretion disk and jet, potentially underestimating or misrepresenting phenomena like turbulent magnetic field structures and dynamo effects. Additionally, 2D simulations cannot fully capture the three-dimensional nature of jet collimation, twisting, and asymmetries caused by disk-jet interactions. As a result, while 2D models offer computational efficiency and valuable insights, they may not accurately represent the complexity and richness of the physics inherent to astrophysical systems, underscoring the importance of 3D simulations for comprehensive studies.

2.2. Radiation model

In this subsection we discuss and present in detail the model that we used for the thermal and nonthermal radiation expected from plasma in the accretion flow. We employed an approximate “proxy” model for radiation that integrates the thermal emissivity across the entire disk. The validity of this method has been tested in the past against general relativistic radiative transfer (GRRT) calculations (Porth et al. 2019). These comparisons demonstrated a close correspondence between the proxy method’s results, the detailed GRRT calculations, and the observed variability in the light curve. Scepi et al. (2022) employed a similar strategy to investigate flares originating from current sheets and plasmoids for MAD models.

The total luminosity L_ν emitted by the fluid at a specific frequency ν (in units of ergs per second per Hz), is calculated based on the following equation:

$$\begin{array}{cc}\hfill L_{\nu }& ={\int }_0^{2\pi }{\int }_0^{\pi }{\int }_{r_h}^{r_{\max }}{e}^{-{\tau }_{\nu }}\times j_{\nu }(B,p,\rho )\times g_{\mathrm{redshift}}^3\hfill \\ \hfill & \sqrt{-g}\mathrm{d}\varphi \mathrm{d}\theta \mathrm{d}r,\hfill \end{array}$$(7)

where j_ν is the emissivity of thermal or nonthermal processes at a specific frequency¹ and τ_ν is the opacity of the fluid to a given process at frequency ν. The terms $\sqrt{-g}\mathrm{d}\varphi \mathrm{d}\theta \mathrm{d}r$ represents the infinitesimal volume element dV for each grid point. Close to the black hole radiation in specific is suppressed by the gravitational redshift term (g_redshift³; Viergutz 1993), where $g_{\mathrm{redshift}}=\sqrt{\mathrm{\Sigma }\mathrm{\Delta }/A}$. Here, A = (r² + a²)² − a²Δsin²θ, Σ = r² + a²cos²θ and Δ = r² − 2rM + a². The parameter a corresponds to the spin of the black hole, while M represents its mass, which is fixed at 1 in these forms (an estimation of g_redshift for r = 2r_g, θ = π/4 and a = 0.93 is approximately on the order of 10⁻¹). The observable radiation flux is given by F_ν = L_ν/4πD², where the distance from the black hole is D = 8.277 kpc (GRAVITY Collaboration 2022). The integration of the luminosity works best when the viewing angle is nearly face on, as we are not taking the effects of Doppler beaming and gravitational lensing into account. The difference of such a method compared with a full GRRT calculation is of order of unity (Scepi et al. 2022).

At 230 GHz, the thermal electron synchrotron radiation comes predominantly from the disk region, namely π/3 ≲ θ ≲ 2π/3. To calculate the synchrotron emission of a thermal distribution of plasma, we utilized the formula developed by Leung et al. (2011, see also Pandya et al. 2016). At 230 GHz, the emissivity of a fluid element in the disk is given by the equation

$$\begin{array}{c}\hfill j_{\nu }\propto n_e{e}^2{\nu }_s{\mathrm{\Theta }}_e^{-2}{c}^{-1}X{e}^{-X^{1/3}}.\end{array}$$(8)

Here n_e is the electron number density and X ≡ ν/ν_s, where ν_s ≡ (2/9)ν_LΘ_e²sin λ (λ is the particle pitch angle with respect to the direction of the magnetic field), and ν_L ≡ eB/(2πm_ec) is the Larmor frequency. The dimensionless proton temperature Θ_p, is defined as Θ_p = k_βT_p/(m_pc²), where k_β is the Boltzmann constant and T_p is the proton temperature. This then calculated from the GRMHD simulation as Θ_p = P/ρ (Porth et al. 2017). We use Θ_e to denote the dimensionless temperature of the electrons, defined as Θ_e = k_βT_e/(m_ec²). Using the above relations, the formula is transformed to

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_e=\frac{m_p}{m_e}\frac{T_e}{T_p}\frac{P}{\rho },\end{array}$$(9)

where m_p and m_e are the masses of the proton and electron, respectively. It is important to mention that the composition of the fluid was assumed to be electron-proton plasma. The ratio T_p/T_e was calculated from the prescription of Mościbrodzka et al. (2016), which describes the electron-to-proton coupling in low (disk) and high magnetized regions (funnel), namely

$$\begin{array}{c}\hfill \frac{T_p}{T_e}=R_{\mathrm{high}}\frac{{\beta }^2}{1+{\beta }^2}+R_{\mathrm{low}}\frac{1}{1+{\beta }^2},\end{array}$$(10)

where R_high and R_low are free parameters that determine the heating ratio of electrons. In our calculations we set R_high = 20 and R_low = 1, since we did not focus on a large parametric investigation of flares from such configurations. In Dihingia et al. (2023a) and Dihingia et al. (2023b), was made a comparison between two-temperature GRMHD model with cooling or no cooling effect and various prescriptions of temperature ratios proposed by various studies (Mościbrodzka et al. 2016; Dihingia et al. 2023a; Meringolo et al. 2023). Some results suggest that, within the disk region, the model proposed by Mościbrodzka et al. (2016) provides a closer approximation to the GRMHD two-temperature cooling model compared to other prescriptions. The electron temperature in a GRMHD simulation can be determined through two temperature simulations that independently evolve this quantity (Jiang et al. 2023). Such studies highlight the importance of self-consistently evolving the electron temperature to accurately explore the thermal synchrotron radiation (Jiang et al. 2024).

The thermal radiation absorption term e^−τ_ν is important only inside the disk. The optical depth is defined as τ_ν = ∫_π/3^θ₀a_ν, thds where $\mathrm{d}s=\sqrt{g_{\theta \theta }}\mathrm{d}\theta$ (π/3 ≤ θ₀ ≤ 2π/3 are the latitudinal limits of the disk). In the case of local thermal equilibrium, the source function of the emitting material is given by the Planck function B_ν. Then, the absorption coefficient of the material can be expressed as

$$\begin{array}{c}\hfill a_{\nu ,th}=j_{\nu }/B_{\nu },\end{array}$$(11)

where

$$\begin{array}{c}\hfill B_{\nu }=(2h{\nu }^3/c^2){[exp(h\nu /k{T}_e)-1]}^{-1}\end{array}$$(12)

is the Planck function. So, the total luminosity for thermal radiation at each grid point (r₀, θ₀) is calculated as follows: First, we determined the emissivity at the point. Next, we computed the absorption term at all grid points (r, θ), where r = r₀ and π/3  ≤  θ  ≤  θ₀. The sum of absorption coefficients (a_ν, th) at these points provides us with the value of the optical depth τ_ν at (r₀, θ₀). In this way, it is like integrating along light rays perpendicular to the equatorial plane.

The nonthermal emissivity at 2.2 μm, which is attributed to the synchrotron radiation of relativistic electrons, is given by (Rybicki & Lightman 1986; Ghisellini 2013)

$$\begin{array}{c}\hfill j_{\nu }\propto \frac{{\sigma }_{T}cK{U}_B}{{\nu }_L}{\left(\frac{\nu }{{\nu }_L}\right)}^{-\frac{p-1}{2}},\end{array}$$(13)

where U_B = B²/(8π) is the energy density of the magnetic field. One needs also to determine the nonthermal electron distribution n(γ) = Kγ^−p, so K is the electron number density multiplied by the efficiency ϵ of the radiation mechanism, and p is the slope of the distribution.

In the case of nonthermal radiation the integration of Eq. (7) is performed along the radial direction. In the funnel region (see, e.g., the colored region in Fig. 3), this is like integrating along rays perpendicular to the equatorial plane. The optical depth in this case is τ_ν = ∫_r^r_maxa_ν, nthds, where $\mathrm{d}s=\sqrt{g_{\mathit{rr}}}\mathrm{d}r$. a_ν, nth is the dimensionless absorptivity; for the nonthermal process, it is defined as

$$\begin{array}{c}\hfill a_{\nu ,nth}=\frac{\sqrt{\pi }{e}^{2}K}{8m_{e}c}{\nu }_L^{\frac{p+2}{2}}{\nu }^{-\frac{p+4}{2}}{f}_a(p)\end{array}$$(14)

(Ghisellini 2013; Pandya et al. 2016), where the function f_a(p) is a product of Γ functions, approximated as

$$\begin{array}{c}\hfill f_a(p)3^{\frac{p+1}{2}}(\frac{1.8}{p^{0.7}}+\frac{p^2}{40}),\end{array}$$(15)

also in this case, the total luminosity for nonthermal radiation at each grid point is calculated as follows: First, we determined the emissivity at the point (r₀, θ₀). Next, we computed the absorption term at all grid points (r, θ), where r₀ ≤ r ≤ r_max and θ = θ₀. The sum of absorption coefficients (a_ν, nth) at these points provides us with the value of the optical depth, τ_ν, at (r₀, θ₀).

A GRMHD simulation does not have the necessary micro-physics to follow the acceleration of particles and extract the electron distribution function or the efficiency for microscopic processes occurring within the plasma. To incorporate how energy is transferred from the magnetic field to the plasma, we used two different post-processing formulas from results of PIC simulations. The first approach is from a local investigation of idealized current sheets (Ball et al. 2018). The second one is from an investigation of microphysical properties of special-relativistic turbulence (Meringolo et al. 2023). The first approach was applied to the current sheets in the funnel region above the disk boundary, and the second one in the vicinity of the disk boundary. The results of such local investigations give us the opportunity to approximate the properties of the accelerated particles that will eventually radiate through the local plasma characteristics. We first defined the nonthermal acceleration efficiency, ϵ, as the fraction of the kinetic energy carried by the nonthermal particles to the kinetic energy of the total electron population,

$$\begin{array}{c}\hfill ϵ=\frac{{\int }_{{\gamma }_{\mathit{pc}}}^{\infty }(\gamma -1)[\frac{\mathrm{d}N}{\mathrm{d}\gamma }-f_{\mathit{MB}}(\gamma ,{\mathrm{\Theta }}_e)]\mathrm{d}\gamma }{{\int }_{{\gamma }_{\mathit{pc}}}^{inf}(\gamma -1)\frac{\mathrm{d}N}{\mathrm{d}\gamma }\mathrm{d}\gamma },\end{array}$$(16)

where γ_pc denotes the peak of the spectrum and f_MB is a relativistic Maxwellian fitting function.

According to Ball et al. (2018), Eq. (17) gives us the slope of the spectrum with respect to the magnetization, σ, and plasma β:

$$\begin{array}{c}\hfill p=A_p+B_{p}tanh(C_p\beta ),\end{array}$$(17)

where $A_p=1.8+0.7/\sqrt{\sigma }$, B_p = 3.7σ^−0.19, and C_p = 23.4σ^0.26. Moreover, the electron nonthermal efficiency with respect to the magnetization (σ) and plasma β is

$$\begin{array}{c}\hfill ϵ=A_{ϵ}+B_{ϵ}tanh(C_{ϵ}\beta ),\end{array}$$(18)

where A_ϵ = 1 − 1/(4.2σ^0.55 + 1), B_ϵ = 0.64σ^0.07, and C_ϵ = −68σ^0.13. On the other hand, Eq. (19) approximates the slope, p (for Eq. (19) it holds that p = k − 1) with respect to the same parameters based on Meringolo et al. (2023):

$$\begin{array}{c}\hfill k=k_0+\frac{k_1}{\sqrt{\sigma }}+k_2{\sigma }^{-6/10}tanh(k_3\beta {\sigma }^{1/3}),\end{array}$$(19)

where k₀ = 2.8, k₁ = 0.2, k₂ = 1.6, and k₃ = 2.25. The efficiency of particle acceleration in turbulent plasmas is equal to

$$\begin{array}{c}\hfill ϵ=e_0+\frac{e_1}{\sqrt{\sigma }}+e_2{\sigma }^{0.1}tanh[e_3\beta {\sigma }^{0.1}],\end{array}$$(20)

where e₀ = 1, e₁ = −0.23, e₂ = 0.5 and e₃ = −10.18.

In Sect. 2.3 we elaborate on how we determined current sheets and their associated parameters, such as the magnetization and plasma β, in their vicinity. The reason we are interested in the latter is because dissipative current sheets are fed with plasma from their environment. These conditions determine the strength of the nonideal electric fields that accelerate particles and the rate of energy dissipation. All empirical relations that describe properties of the accelerated particle populations refer to the σ and plasma β parameters in the upstream regions of the current sheets. Once these areas of interest are identified, these formulas can then be applied to determine the slope of the distribution of accelerated electrons within the current sheets or turbulent regions at the disk boundary.

2.3. Determination of current sheets and their environment

The PIC simulations that investigate current sheets start with a specific magnetization, a specific plasma β and an initial current layer in the middle. As the magnetic field evolves magnetic reconnection events take place, magnetic energy is dissipated, and particles are accelerated from the thermal pool to higher energies. Equations (17) and (19) are fitting models of the results of such PIC simulations and give us the slope of the spectrum of the nonthermal electrons as functions of GRMHD quantities. It is important to emphasize that the parameter values that must be used in these formulas are roughly those set as initial conditions in a PIC simulation. As it turns out, these values are maintained during the calculation only in the environment surrounding the current sheets.

So, in order to use Eqs. (17) and (19), it is necessary to constrain the current sheet and its surrounding environment. First, it is useful to distinguish the regions where current sheets form from the polarity reversal of the magnetic fields (Fig. 2, reversal of B_z). We utilized five variables to determine the current sheets, which are as follows: toroidal current J_ϕ := (∇×B)_ϕ, density ρ, magnetization σ, plasma β and dimensionless temperature Θ_p (normalized to mc², which is equal to P/ρ). The remaining plots in Fig. 2 provide limits for each parameter that defines the exact region where a current sheet occurs. As we see in Fig. 2, a current sheet is characterized by high values of toroidal current, high values of density, low values of magnetization, high values of plasma β and high values of dimensionless temperature. Thus, we set the following limits for each parameter: |J_ϕ|> 10⁻⁴, a density floor (cutoff) at ρ_cut = 2 × 10⁻⁵, a magnetization ceiling cutoff at σ_cut = 10, a plasma β floor at β_cut = 10⁻³ and finally a temperature floor at Θ_p = 10⁻³. By combining these constraints, we obtained the red regions in Fig. 3, which are clearly associated with the current sheets we identified. Additionally, the plasma in these red-colored regions has undergone the reconnection process. Efficient particle acceleration occurs in regions with high magnetization and low plasma β (Ball et al. 2018). Such regions develop only in the funnel region, so we applied the above methodology exclusively to the plasma within this region.

Fig. 2. Upper panel (from left to right): Magnetic field Bz component, toroidal current, and density. Lower panel (from left to right): Magnetization, plasma β, and dimensionless temperature Θp. By setting limits on these parameters, we define the red regions in Fig. 3.

Fig. 3. Determination of plasma (red) in the reconnection region and places of high magnetic turbulence. They are selected after applying limits on the toroidal current, magnetization, plasma β, density, and dimensionless temperature. The green grid corresponds to the discretization for the purpose of parameter averaging. The snapshot is the same as that shown in Fig. 2.

Following this procedure, we identified the current sheets developed in the simulation. The next step involved characterizing the current sheet environment in terms of averaged quantities in order to be able to use Eqs. (17) and (19) to determine the plasma within the current sheet. To achieve this, a discretization method was employed, wherein the region was divided into cells with a radial length of 1 r_g and a meridional extent of 10 degrees in the θ direction. In each cell that contains part of a current sheet (red region), we calculated the mean values of magnetization, σ, and plasma β, of their environment (white region) over two cells to the right and two to the left in the θ direction, namely

$$\begin{array}{c}\hfill \overline{\sigma }=\frac{1}{4}\sum _{\mathrm{cell}=-2}^{\mathrm{cell}=+2}\sigma \end{array}$$(21)

$$\begin{array}{c}\hfill \overline{\beta }=\frac{1}{4}\sum _{\mathrm{cell}=-2}^{\mathrm{cell}=+2}\beta .\end{array}$$(22)

Subsequently, the slope of the nonthermal electrons inside the current sheet or in the turbulent plasma in each cell is computed using Eqs. (17) and (19) with parameters values as $\sigma =\overline{\sigma }$ and $\beta =\overline{\beta }$. The plasma assumed to be accelerated corresponds to the region depicted in red in Fig. 3.

3. Results

In the next two subsections, we present our results on the thermal and nonthermal radiation. The first part covers calculations at 230 GHz (thermal component), which help calibrate the physical parameters for Sgr A*. Further calculations provide information about the model’s variability compared to observations. Subsequently, we present results from the nonthermal radiation calculations at 2.2 μm, including the resulting light curve, an analysis of flares, and a comparison with the GRMHD simulation results.

3.1. Radiation at 230 GHz

Observations of Sgr A* at 230 GHz have been conducted by various telescopes, including the Submillimeter Array (SMA), the Atacama Large Millimeter/submillimeter Array (ALMA), and notably, the Event Horizon Telescope (EHT). The EHT in particular has provided a wealth of high-resolution images from the innermost region of the accretion disk onto Sgr A* itself, presenting an invaluable resource for the scientific community (EHT Collaboration 2022a,b,c,d,e,f). Thus, there exists a plethora of data on the properties of the light curve at that particular frequency. Therefore, utilizing much of these data, our main goals are to reproduce the main features of the light curve at 230 GHz (amplitude and variability) by calibrating the various physical parameters from dimensionless code units to actual physical units. When we combine how the flux changes with distance and angle (as illustrated in Fig. 4), it becomes clear that the majority of the radiation is emanating from a distance of approximately 10 r_g. This underscores that the events we are studying occur in extremely close proximity to the black hole, and consequently, the critical data points we are interested in are highly concentrated in this area.

Fig. 4. Cumulative 230 GHz flux from all angles θ for radii less than or equal to r.

As we know from observations (Wielgus et al. 2022), the radiation flux of Sgr A* at 230 GHz is approximately stable. To turn the parameters of the GRMHD simulation from code units to physical cgs units, a conversion and calibration of variables is necessary. We defined a scaling parameter (S) that appropriately multiplies physical quantities such as the magnetic field (B), pressure (P), and density (ρ) to convert values of physical quantities from geometric (code) units to cgs units. Obviously, the conversion leaves dimensionless quantities such as Θ_e, σ, and β unchanged:

$$\begin{array}{cc}\hfill B_{\mathit{cgs}}& =\sqrt{4\pi }c{S}^{1/2}{B}_{\mathrm{geometric}}\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill {\rho }_{\mathit{cgs}}& =S{\rho }_{\mathrm{geometric}}\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill P_{\mathit{cgs}}& =c^{2}S{P}_{\mathrm{geometric}}.\hfill \end{array}$$(25)

Through an iterative procedure applied in various snapshots of the GRMHD calculation, the model limited the radiation flux in a range around 4 Jy and finally produced the value of the scaling parameter S. After calculating the parameter S, we could calculate the radiation flux at all snapshots and finally produce the light curve at 230 GHz. The upper panel of Fig. 5 presents the light curve at 230 GHz (red line). Noticeably, the radiation flux remains almost constant throughout the calculation period. The higher flux values observed in the beginning of the light curve gradually decrease toward the end of the evolution.

Fig. 5. Top panel: Simulated light curve at 230 GHz (red line), the moving average of the simulated light curve (m3, in a time window of 3 hours; black line), and the 1s3 standard deviation (gray area). Lower panel: Measure of the variability of the simulated light curve (M3), defined as M3 = s3/m3.

To investigate the variability of the light curve, we defined the moving average radiation flux over a 3-hour period (m₃, black line in Fig. 5) and the standard deviation over the same duration (s₃, gray area in Fig. 5). The lower panel of Fig. 5 shows the spread of M₃ where M₃ = s₃/m₃. M₃ is a measure of the variability of the light curve and can be compared with historical observations from Sgr A* (Fig. 6).

Fig. 6. Distribution of the M3 of the simulated light curve at 230 GHz (blue) compared with historical data from observations of Sgr A* (gray). We have divided the simulated light curve and its respective M3 into three time intervals: t ∈ [0, 370] min (green), t ∈ [370, 746] min (red), and t ∈ [746, 1143] min (yellow).

The Sgr A* light curve shifts and changes over time due to a combination of stellar winds and turbulence on various scales. Longer changes come from variations in the stellar winds, especially near the S stars. Shorter changes happen because of turbulence closer to the center. The spread of M₃ in Fig. 5 exhibits a downward trend, suggesting that if we keep running the same simulation for further time, the flare variability will eventually stabilize at lower levels. Figure 6 is a histogram of M₃ that provides a broader view derived from observing Sgr A* (see Wielgus et al. 2022); the other distributions are derived from the simulation results for different time windows. Our model appears to align well with the observations. As previously mentioned, our GRMHD simulation corresponds to the SANE configuration (SANE multi-loop). According to the EHT publication (EHT Collaboration 2022e), these models demonstrate better consistency with the variability constraints compared to MAD models. It is important to keep in mind about Fig. 6 that the yellow line distribution represents the last part of the M₃ curve inside the time interval t ∈ [746, 1143] min. This closely aligns with the observations and shows that a new calculation of the same model for longer times will probably yield the desired variability in the light curve. However, for a more accurate study of variability, it would be necessary to incorporate phenomena such as Doppler beaming and lensing effects into the calculations. These aspects were deliberately excluded, as their inclusion would deviate from the primary objective of this work.

3.2. Results at 2.2 μm (near-infrared flares)

3.2.1. Light curve and statistics

In Fig. 3 we carefully mark the specific areas where nonthermal radiation processes occur. These dynamic areas reside within the funnel region and in the limit of the disk. There is an environment characterized by high magnetization and low plasma (β) values, which are needed to activate the radiation mechanisms with sufficient efficiency. In Sect. 2.3 we discuss in detail that it is this local environment that characterizes the process of acceleration for the plasma residing within current sheets.

Using data from these regions, we calculated the light curve at 2.2 μm, and our results are depicted in Fig. 7. Notably, during quiescent states, the average flux is remarkably low, consistent with observations. Within the simulation, several flares occur, each typically exceeding a flux of > 1.5 mJy (Witzel et al. 2021; von Fellenberg et al. 2023). Noteworthy is that fact that there are four flares during the simulation time that surpass this threshold, with an additional three flares marginally meeting the criteria for acceptance as flares. One particularly significant radiation flare, with a flux of approximately 7 mJy and a longer duration than typical flares, stands out as a potential candidate within the observed limits for both flux and duration, resembling the unprecedented flare reported in Do et al. (2019).

Fig. 7. Simulated light curve at 2.2 μm. Our model can produce one very bright flare of around 7 mJy and another three that pass the 1.5 mJy threshold.

In a theoretical comparison with observations, we do not have a flux limit that allows us to define what is considered a flare and what not, but it is obvious that the values of the flux in Fig. 7 increase by several orders of magnitude between a quiescent and a flare state (F ≳ 1.5 mJy). Sgr A* produces flares with lower flux but also some very bright ones (F > 6 mJy; Do et al. 2019) like those produced in our model. In Fig. 8 we see a comparison of the historical flux distribution in histogram for Sgr A* from Witzel et al. (2018) as presented in Do et al. (2019) compared to the high flux flare observed in 2019. The probability of observing such a flare (> 6 mJy) was computed to be less than 0.05% (Witzel et al. 2018). Our model is capable of reproducing high-flux flares that cover the full range of observations.

Fig. 8. Histogram comparing our model with observations from Do et al. (2019). The lower panel is a zoomed-in view.

According to Fig. 8 our model appears to have fairly good statistics in terms of the population of flares. This is not the case for Sgr A* flares generated from MAD models, which seem to show a very large population of flares and tend to overproduce the quiescent state flux (Dexter et al. 2020; Scepi et al. 2022; White & Quataert 2022). This is probably due to the very strong magnetic field they possess, and the fact that the main source of the flares is the current sheet in the equator, which feeds high values of absorptivity in the density due to the ̇M (as seen characteristically in the paper of Scepi et al. 2022).

It is important to remark that in our method for computing the light curve, we have not taken the effect of cooling into account. For reference, we estimated the cooling time of a particle with The Lorentz factor, γ, due to synchrotron radiation was calculated as

$$\begin{array}{c}\hfill t_{\mathrm{syn}}=\frac{\gamma m_e{c}^2}{4/3{\sigma }_{T}c{U}_B{\gamma }^2{\beta }^2},\end{array}$$(26)

where the Lorentz factor, γ, is γ = (ν/ν_L)^1/2. Using typical values for magnetic field B from our model, approximately 60 G for the region of interest calculating the cooling time from relation (26) we obtain t_syn ∼ 3 min. In future research, cooling should be examined more thoroughly, as it can influence the light curve results at near-infrared and higher wavelengths.

3.2.2. Properties of flares

In the upper panel of Fig. 9 we zoom in on the first major flare within the time range t ∈ [86, 150], hereafter referred to as Flare A. In the lower panel, we zoom in on a second, smaller flare within the time range t ∈ [879, 926], hereafter referred to as Flare B. We examined how the properties of each flare change across three different stages: the pre-flare state, which occurs when the flux is several orders of magnitude smaller than 1 mJy; the dim state, which occurs when the flux begins to increase (as shown in Fig. 11, where the spectral index also increases) and finally, the bright state, which occurs when the flux exceeds 1 mJy. In both flares, we can see the specific points and their corresponding colors that represent the three different states. Each point is derived from a specific snapshot of the GRMHD simulation. The first flare is characterized by a distinct peak, while the second flare has two peaks. Similar phenomena have been observed (GRAVITY Collaboration 2018a). The classification of points into the different flare states was based on flux and spectral index. The rationale for this classification is clearly illustrated in Fig. 11.

Fig. 9. Part of the light curve at 2.2 μm. Top: Flare A during the period t ∈ [86, 150] min (the first bright flare with approximately 7 mJy radiation flux). Bottom: Flare B during the period t ∈ [879, 926] min (a bright flare with two peaks). The duration of the flares is about 40–60 minutes. We divided the light curve into three states: pre-flare, dim state, and bright state. This separation is based on the flux and spectral index

Following the same calculation methods used for the nonthermal radiation, we estimated the radiation flux at three additional frequencies: 6.66, 7.95, and 18.1 GHz, which correspond to the M, L, and H bands, respectively. The primary frequency of our calculation, 2.2 μm, corresponds to the K band for Sgr A*. For this calculation we used two snapshots at times t = 136 min and t = 920 min, one for each flare at the peak of the bright state. We thus extracted the spectrum of our model for this flare state (Fig. 10). The fitting curves in Fig. 10 can give us the spectral index (which is defined as the exponent of a power-law-like F ∝ ν^a) for each flare. The two values that we obtained (a = −0.46 and a = −0.44) are within the observed range (Ghez et al. 2005a,b; Bremer et al. 2011).

Fig. 10. Radiation flux around 2.2 μm for the frequencies 66.6 (M band), 79.5 (L band), 138 (K band), and 181 THz (H band) for Flare A (stars) and Flare B (circles). The fitting lines (black and orange) correspond to power laws where the exponential gives us the spectral index of the spectrum, in our case a = −0.44 for Flare A and a = −0.46 for Flare B.

During a flare, the flux and the spectral index increase simultaneously. We categorized the points into three groups, distinguished by the colors used to represent them. Figure 11 shows the evolution of the spectral index as the flares transition from the pre-flare state to the dim state and finally to the bright state. The fitting power-law curve (orange line) approximates how the spectral index evolves as the flare transitions from one state to the next, which is comparable to observations (Bremer et al. 2011). The same colors were used to distinguish the points in Fig. 9.

Fig. 11. Spectral index vs. flux for Flare A (top) and Flare B (bottom). The dots correspond to the snapshots that are marked in Fig. 9: the pre-flare (cyan), dim (green), and bright states (red) of each flare. The orange lines are a fitting power-law curve, which can be compared with observations.

3.2.3. Flare states compared with the GRMHD simulation

We wanted to see how the peaks in the light curve relate to the GRMHD calculation from which the phenomena originated. Therefore, we plotted the parameters ρ, σ, and plasma β, which play crucial roles in shaping the current sheet (Fig. 12). The plots refer to the three states of the second flare we investigated, specifically the snapshots at times t = 886 min (pre-flare), t = 900 min (dim state), and t = 910 min (bright state). The corresponding points of these snapshots can be seen in the light curve of Fig. 9.

Fig. 12. Top row: Density of Flare B (t ∈ [886, 963] min) in its three states: pre-flare, dim state, and bright state, at snapshots t = 886 min, t = 900 min, and t = 910 min, respectively. Middle row: Magnetization (σ) at the same snapshots. Bottom row: Plasma β at the same snapshots.

The plot of the density (ρ; top row in Fig. 12) shows a clear structure in the upper funnel region that appears as the flare transitions from one state to the next. This structure is also characterized by low magnetization and high plasma β. According to the definition given above, it constitutes a clear current sheet that, in the bright state, is surrounded by regions with high values of magnetization (σ; middle row in Fig. 12) and low plasma β (bottom row in Fig. 12). These are ideal conditions for radiating nonthermal radiation. This current sheet in the bright state can be characterized as a plasmoid chain, similar to those seen in localized PIC simulations (Ball et al. 2018; Sironi & Spitkovsky 2014; Petropoulou & Sironi 2018).

It is important to emphasize that this increase in magnetization, σ, and decrease in plasma β in the bright state of the flares (conditions ideal for the activation of flares) are also found at the boundaries of the disk. These are possible regions for flare generation through high magnetic turbulence in the plasma, as mentioned in Sect. 2.2. Figure 12 clearly shows that the creation of the current sheet inside the funnel region is a phenomenon perfectly connected to the organization of the magnetic field and results in the creation of conditions (magnetization and plasma β) that will activate nonthermal particles in the current sheet itself and along the disk boundary.

Indeed, the areas that radiate are the current sheet and the disk boundary as shown in Fig. 13. Figure 13 illustrates the radiation flux (in mJy), the slope p and the efficiency ϵ at time t = 910 min (the bright state of Flare B that we investigated). The visible boxes in each plot correspond to the identification of current sheets and their environment as shown in Fig. 3. Within each box, all plasma quantities are averaged in order to characterize this specific reconnection region. However, only the plasma within the current sheet layer – determined by the constraints and limits discussed in Sect. 2.3 – is considered to have the specific slope and efficiency, producing the reported flux. Two-thirds of the nonthermal radiation flux comes from the disk boundary, where the turbulent plasma dominates, while one-third comes from the current sheet.

Fig. 13. Radiation flux at 2.2 microns (in mJy), slope (p), and efficiency (ϵ) in each grid point of the GRMHD simulation shown in the left, middle, and right panels, respectively, for the bright state of Flare B.

Similar phenomena can be observed in Flare A at time t = 136 min. A clear current sheet has produced a plasmoid chain, and the local environment is characterized by high magnetization and low plasma β (see Fig. 15). Accordingly, the slope p and the efficiency ϵ (middle and right panel of Fig. 14) peak at the current sheet and produce the radiation flux (left panel of Fig. 14), which strengthens our conclusions about how the flares are generated in this particular calculation.

Fig. 14. Same as Fig. 13 but for the bright state of Flare A.

Fig. 15. Density (left panel), magnetization, σ (middle panel) and plasma β (right panel) of Flare A in the bright state.

4. Conclusions

We have analyzed data from the GRMHD simulation of Nathanail et al. (2020), specifically focusing on model D, which is characterized as a SANE multi-loop model. Unlike MAD simulations, this model describes the accretion disk of a black hole without the production of a stable jet.

Applying a thermal radiation model to this particular simulation enabled us to reproduce the light curve at 230 GHz, from which we were able to calibrate the simulation variables. Encouragingly, the model’s variability demonstrated good agreement with observations, prompting us to conduct larger simulations with extended time evolution.

To calculate the nonthermal radiation from the simulation, we introduced a novel method that identifies current sheets and places with high magnetic turbulence. This identification, based on the magnetic field polarity reversals, is done by constraining primarily the current density together with the micro-physical plasma parameters, such as magnetization and plasma β. After the identification, all quantities are averaged in the local environment. The constraints we placed on the parameters provided clear spatial boundaries within which to apply the model that calculates the nonthermal radiation. This association both advances our comprehension of the fundamental processes at play and provides a crucial framework for future investigations in similar contexts.

The results for the light curve at 2.2 μm are very encouraging as they produce both small flares up to 2 mJy and a couple of brighter ones that reach 7 mJy. The duration of such flares is also consistent with observations. Further analysis of the flare spectral index verified the success of the model in reproducing the observations. It is interesting that, during the evolution of the flare, the spectrum follows a power law similar to those given by the observations with spectral index a = −0.44.

In summary, we emphasize once again that a very important result of this work is the identification of the source of flaring nonthermal radiation in GRMHD simulations. We have shown that such flares most likely originate from current sheets and their associated plasmoid chains in the funnel region, as well as from the disk boundary due to magnetic turbulence.

5. Data availability

The data underlying this article will be shared upon reasonable request to the corresponding author (URL source: Zenodo Record

Acknowledgments

ID is supported by the Hellenic Foundation for Research and Innovation (HFRI) under the 4th Call for HFRI PhD Fellowships (Fellowship Number: 9239). Support comes from the ERC Advanced Grant “JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales” (Grant No. 884631). This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility – ARIS. Simulations were performed also on the GOETHE-HLR cluster at CSC in Frankfurt).

References

All Figures

	Fig. 1. Initial magnetic field (lines) and density (color) configuration in the initial torus. Axis coordinates are calculated in rg. Red and blue denote a clockwise and counterclockwise direction of the magnetic field.
In the text

	Fig. 2. Upper panel (from left to right): Magnetic field Bz component, toroidal current, and density. Lower panel (from left to right): Magnetization, plasma β, and dimensionless temperature Θp. By setting limits on these parameters, we define the red regions in Fig. 3.
In the text

	Fig. 3. Determination of plasma (red) in the reconnection region and places of high magnetic turbulence. They are selected after applying limits on the toroidal current, magnetization, plasma β, density, and dimensionless temperature. The green grid corresponds to the discretization for the purpose of parameter averaging. The snapshot is the same as that shown in Fig. 2.
In the text

	Fig. 4. Cumulative 230 GHz flux from all angles θ for radii less than or equal to r.
In the text

	Fig. 5. Top panel: Simulated light curve at 230 GHz (red line), the moving average of the simulated light curve (m3, in a time window of 3 hours; black line), and the 1s3 standard deviation (gray area). Lower panel: Measure of the variability of the simulated light curve (M3), defined as M3 = s3/m3.
In the text

	Fig. 6. Distribution of the M3 of the simulated light curve at 230 GHz (blue) compared with historical data from observations of Sgr A* (gray). We have divided the simulated light curve and its respective M3 into three time intervals: t ∈ [0, 370] min (green), t ∈ [370, 746] min (red), and t ∈ [746, 1143] min (yellow).
In the text

	Fig. 7. Simulated light curve at 2.2 μm. Our model can produce one very bright flare of around 7 mJy and another three that pass the 1.5 mJy threshold.
In the text

	Fig. 8. Histogram comparing our model with observations from Do et al. (2019). The lower panel is a zoomed-in view.
In the text

Fig. 9. Part of the light curve at 2.2 μm. Top: Flare A during the period t ∈ [86, 150] min (the first bright flare with approximately 7 mJy radiation flux). Bottom: Flare B during the period t ∈ [879, 926] min (a bright flare with two peaks). The duration of the flares is about 40–60 minutes. We divided the light curve into three states: pre-flare, dim state, and bright state. This separation is based on the flux and spectral index

In the text

	Fig. 10. Radiation flux around 2.2 μm for the frequencies 66.6 (M band), 79.5 (L band), 138 (K band), and 181 THz (H band) for Flare A (stars) and Flare B (circles). The fitting lines (black and orange) correspond to power laws where the exponential gives us the spectral index of the spectrum, in our case a = −0.44 for Flare A and a = −0.46 for Flare B.
In the text

	Fig. 11. Spectral index vs. flux for Flare A (top) and Flare B (bottom). The dots correspond to the snapshots that are marked in Fig. 9: the pre-flare (cyan), dim (green), and bright states (red) of each flare. The orange lines are a fitting power-law curve, which can be compared with observations.
In the text

	Fig. 12. Top row: Density of Flare B (t ∈ [886, 963] min) in its three states: pre-flare, dim state, and bright state, at snapshots t = 886 min, t = 900 min, and t = 910 min, respectively. Middle row: Magnetization (σ) at the same snapshots. Bottom row: Plasma β at the same snapshots.
In the text

	Fig. 13. Radiation flux at 2.2 microns (in mJy), slope (p), and efficiency (ϵ) in each grid point of the GRMHD simulation shown in the left, middle, and right panels, respectively, for the bright state of Flare B.
In the text

	Fig. 14. Same as Fig. 13 but for the bright state of Flare A.
In the text

	Fig. 15. Density (left panel), magnetization, σ (middle panel) and plasma β (right panel) of Flare A in the bright state.
In the text