The Properties of Reconnection Current Sheets in GRMHD Simulations of Radiatively Inefficient Accretion Flows

General relativistic magnetohydrodynamic (GRMHD) simulations are often used to study the physics of accretion systems around compact objects and to explain their observed properties. Such simulations employ the assumption of "ideal" MHD, which enforces that the plasma is infinitely conductive. This assumption leads to a couple of important properties: any electric fields in the fluid frame are shorted out and, as a result, the magnetic fields are frozen into the fluid (see, e.g., Kulsrud 2005).

For many systems, the ideal approximation is adequate: astrophysical plasmas often have extremely low resistivities. However, even in systems where the approximation may globally seem appropriate, there can arise regions that violate the underlying assumptions. One such example is when two regions of plasma with opposing magnetic flux encounter each other. This may be a site of reconnection, where magnetic fields can change topology and dissipate their energy into the plasma. Such a configuration will have a rapidly changing magnetic field in space, resulting in a high current density. Ohmic dissipation, which scales as $\eta {J}^{2}$, where η is the resistivity and J is the current density, can then lead to significant dissipation where the current is high enough. This shows that dissipative terms can change the energetics of the flow and, given their localization, can lead to time-dependent phenomena that are not properly captured in global ideal MHD simulations.

In the case of magnetic reconnection, even the inclusion of non-ideal terms may not be sufficient to capture the entire behavior of the plasma. For instance, reconnection has been shown to be an efficient source of non-thermal particle acceleration under certain conditions (Drake et al. 2013; Melzani et al. 2014; Sironi & Spitkovsky 2014; Guo et al. 2015; Li et al. 2015; Sironi et al. 2016; Werner et al. 2016; see Kagan et al. 2015 for review), with the most efficient electron acceleration occurring in highly magnetized and cold plasmas. MHD, however, cannot capture the physics of non-thermal particle acceleration in reconnection. In these cases, we must turn to computational methods of solving the Vlasov equation that do not make assumptions about the particle distribution, such as particle-in-cell (PIC) simulations.

PIC simulations are applicable to collisionless systems, where particle acceleration is possible due to the long timescale of Coloumb collisions. The assumption of a collisionless plasma employed in PIC simulations may appear at odds with the use of MHD for the global flow, which is based on a fluid (i.e., collisional) description of the plasma. In principle, both assumptions can be satisfied at the same time if the collisional mean free path is much larger than the characteristic length of the reconnection region, but smaller than the scale of the global flow. Reconnection occurs in collisionless systems through kinetic effects acting as an "effective resistivity," whereby electromagnetic fields scatter particles, allowing for reconnection. The precise origin of the resistivity operating in collisionless reconnection is not yet understood.

The supermassive black hole at the center of the Milky Way, Sgr A*, has an accretion flow that falls into a category broadly referred to as radiatively inefficient accretion flows (see Yuan & Narayan 2014 for a recent review). These flows are characterized by geometrically thick, optically thin disks, low accretion rates, and low luminosities. In some accretion rate regimes, such as that of Sgr A*, these accretion flows are not fully collisional. Nevertheless, because of current computational capabilities, GRMHD has been employed in a number of studies to infer properties of the accretion flow around Sgr A* (e.g., Dexter et al. 2012; Narayan et al. 2012; Drappeau et al. 2013; Mościbrodzka et al. 2014; Chan et al. 2015b). Even more recently, studies have begun to evolve the electron entropy equation, accounting for electron heating and anisotropic conduction (Ressler et al. 2015). Additionally, Chael et al. (2017) developed a scheme for coevolving a population of non-thermal electrons, including effects of adiabatic compression and expansion as well as radiative cooling. While these simulations successfully match a number of broadband steady-state properties, they show very little X-ray variability, contrary to observations (Eckart et al. 2004; Eckart et al. 2006; Neilsen et al. 2013). In Ball et al. (2016), we showed that a population of non-thermal electrons in highly magnetized regions of a radiatively inefficient accretion flow studied with GRMHD, where they are likely to be accelerated via reconnection, can result in X-ray variability with properties that are roughly consistent with observations.

In this paper, we use representative GRMHD simulations to assess whether reconnection regions frequently occur in global simulations. We consider simulations with standard and normal evolution (SANE) and magnetically arrested disk (MAD) initial magnetic field configurations (see Narayan et al. 2012 & Sa̧dowski et al. 2013 for a more technical discussion of these initial field geometries). In the SANE case, the magnetic field is initialized with alternating poloidal loops, while the MAD initial field consists of a single poloidal loop. These initial field geometries play a large role in the evolution of the accretion flow. In MAD simulations, because the initial field is uniformly aligned, it quickly builds up large-scale toroidal flux, which inhibits turbulence in the innermost accretion disk. In the SANE case, because the initial field alternates in polarity, there is less net flux through the disk when the poloidal loops are sheared, hence more turbulence develops.

We devise criteria to locate regions of field reversal and characterize the properties of the plasma in these regions. We focus on the plasma-β and magnetization parameter σ, which have been shown to play an important role in particle acceleration. We also identify field components that are orthogonal to the reversing field, often referred to as guide fields, and quantify their strengths. Our results will guide future PIC studies of low-luminosity accretion flows. Finally, we compute the time-dependent magnetic energy available in reconnection regions to assess whether this is a plausible mechanism to generate the observed X-ray variability of Sgr A*.

Magnetic reconnection takes place in regions where there is a reversal of magnetic field over a short characteristic length scale in which the current density becomes large. In typical simulations of the local dynamics of reconnection, the initial condition is specified in terms of a Harris sheet, which has the magnetic field profile

$$\begin{eqnarray}&&{\boldsymbol{B}}={B}_{0}\tanh \displaystyle \frac{x}{L}\hat{y}.\end{eqnarray} \tag{ 1 }$$

In this geometry, the y-component of the magnetic field reverses direction over a characteristic length L in the x-direction. This field reversal has a high curl associated with it, leading to a sudden peak in the current density, which scales as

$$\begin{eqnarray}&&{\boldsymbol{J}}=\displaystyle \frac{{B}_{0}}{L}{{\rm{sech}} }^{2}\displaystyle \frac{x}{L}\hat{z}.\end{eqnarray} \tag{ 2 }$$

There are only a small number of parameters that determine the particle heating and acceleration that results from reconnection events. These are the magnetization parameter

$$\begin{eqnarray}&&\sigma \equiv \displaystyle \frac{{B}^{2}}{4\pi \rho {c}^{2}},\end{eqnarray} \tag{ 3 }$$

which is the ratio of magnetic energy density to rest mass energy density, and the plasma-β parameter

$$\begin{eqnarray}&&\beta \equiv \displaystyle \frac{{P}_{\mathrm{gas}}}{{P}_{\mathrm{magnetic}}}=\displaystyle \frac{8\pi {nkT}}{{B}^{2}},\end{eqnarray} \tag{ 4 }$$

which specifies the ratio of gas pressure to magnetic pressure. Here, ρ, n, and T are the mass density, number density, and temperature of the plasma particles, respectively.

Another important quantity to consider for magnetic reconnection is the magnitude and direction, if present, of the so-called guide field. This is the component of the magnetic field in the sheet perpendicular to the reconnecting field. The effect of such a guide field on particle acceleration has been studied in certain regimes (Dahlin et al. 2016; Stanier et al. 2016; Wang et al. 2016) and, in some cases, can have an effect on the resulting electron energy distribution.

Our first goal is to devise an algorithm that will allow us to identify the location and relevant properties of potential reconnection regions, i.e., Harris sheets, in global GRMHD simulations, which we describe in the following section.

As an illustrative example, we use two 60 hr (about 11,000 GM/c³) long GRMHD simulations of a radiatively inefficient accretion flow onto a black hole (Narayan et al.2012; Sa̧dowski et al. 2013), which were performed using the HARM code (Gammie et al. 2003). These simulations were employed in a large study of the broadband, time-dependent emission from Sgr A* (Chan et al. 2015a, 2015b) where we coupled HARM to the radiative transfer algorithm GRay (Chan et al. 2013) and varied the black hole spin, density normalization, observer inclination, initial magnetic field configuration, as well as the electron thermodynamic prescription. These simulations resolve the magnetorotational instability and are run for sufficient time to reach a steady accretion rate. From these investigations, we identified five models that best fit the steady-state broadband spectrum, as well as the previously observed 1.3 mm image size of Sgr A*, and also characterized their variability properties. In the present study, we use two representative models from Chan et al. (2015b): a SANE model with a black hole spin a = 0.7 and a MAD model with a black hole spin a = 0.9. In general, the thermal SANE models tend to show short-lived, high-amplitude variability in their IR and mm flux, while the thermal MAD models tend to show light curves dominated by smooth and long-timescale flux changes (Chan et al. 2015a). The stark differences in variability properties between SANE and MAD simulations can be explained by how the initial field configuration evolves and affects the dynamics in the disk. In MAD simulations, the initial field is uniformly aligned (i.e., no changes in polarity). As a result, when the field is advected inward and sheared, the disk quickly builds up magnetic flux, resulting in strong toroidal fields that play a dominant role in the dynamics of the disk. The large-scale fields inhibit turbulence, resulting in less variability in MAD simulations (see, e.g., Chan et al. 2015a; Medeiros et al. 2016). The SANE models, however, are initialized with alternating poloidal loops, which results in many regions of opposing toroidal flux in the disk; as the loops are sheared and advected inward, there will be many regions of opposing toroidal flux in the disk, hence a smaller net flux throughout the disk than in the MAD case. Therefore, the magnetic fields are not strong enough to suppress turbulence as effectively as in MAD configurations, leading to more variability in SANE disks.

Our goal is to identify in each snapshot from these simulations potential regions of reconnection. Because of the large shear in the accretion flow, the magnetic fields are primarily toroidal and the alternating components occur primarily in the azimuthal direction. For this reason, we search through the simulation volume for cells that have both a high current relative to the mean value in the snapshot, as well as very low values of the azimuthal component of the magnetic field B_ϕ to pick out the sheets where reconnection may occur.

Specifically, we consider 2D slices of the simulation volume at each azimuthal angle ϕ at each snapshot and identify the points that (i) have current magnitudes $\sqrt{{J}_{\mu }{J}^{\mu }}$ that are higher than four times the mean current of that snapshot and (ii) ϕ-components of the magnetic field smaller than a fiducial value, characterized by the usual magnetization parameter ${\sigma }_{\phi }\,={B}_{\phi }^{2}/(4\pi {{mnc}}^{2})$. We use a ${\sigma }_{\phi }$ threshold of 10⁻⁶. We then apply an algorithm similar to the one described in Zhdankin et al. (2013) for identifying and analyzing the statistics of current sheets in shearing box simulations of MHD turbulence. For every point with grid indices ($i,j$) on an azimuthal slice of our domain picked out by the above criteria with current magnitude ${J}_{\mathrm{peak}}$, we consider all four adjacent points in the grid. If the current at an adjacent point is above ${J}_{\mathrm{peak}}/2$, while also satisfying ${\sigma }_{\phi }\lt {10}^{-6}$, we consider it as part of the same current sheet. We continue this process of scanning every point in the sheet, considering all neighboring points, and applying these criteria to them until no more points are being added to the sheet.

In the top panel of Figure 1, we show the result of applying this algorithm on a snapshot of the SANE simulation, with the B_ϕ configuration shown in the bottom panel of the same figure. It is evident that the regions between flux tubes of opposing azimuthal magnetic flux are effectively picked out. When we repeat this procedure on all adjacent azimuthal slices, we find that the current sheets show large azimuthal extents throughout the flow, providing ample surface area for neighboring flux tubes to reconnect over. We show in Figure 2, a representative snapshot of the MAD simulation for comparison. We see that the magnetic field strengths about the current sheet in the MAD simulation are higher than that in SANE. Additionally, we find that there is typically less fine structure to the MAD current sheets; the MAD current sheets tend to consist of one or two relatively flat sheets, while the SANE current sheets are often highly curved and twisted into complicated geometries, as shown in the top panel of Figure 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Sampling The Plasma Properties Associated with Current Sheets

Once we identify the current sheets in each snapshot, we characterize the plasma parameters of these sheets that are relevant to magnetic reconnection. The location where we want to measure these parameters is not in the sheet itself, but where the magnetic field reaches its asymptotic value some distance away from the sheet in a direction perpendicular to it. This breaks down into two problems: finding the direction perpendicular to the sheet and determining how far to go along this direction until the magnetic field reaches its appropriate asymptotic value.

To approximate the direction perpendicular to the current sheet at a point ($i,j$) that has been flagged as belonging to the sheet, we first find the local slope of the sheet about this point. To do this, we consider a box around each point ($i,j$) in the sheet, with width S + 1, whose corners are at ($i\pm S/2,j\pm S/2$). We use a value of S = 10 pixels, which is generally smaller than the radius of curvature of a current sheet. We then calculate the slope from the point in the center of this box ($i,j$) to every other point (${i}^{{\prime} },{j}^{{\prime} }$) in the box which is flagged as being part of the current sheet. Taking the inverse tangent of this slope gives the angle with respect to the horizontal of the line that passes through ($i,j$) and (${i}^{{\prime} }$, ${j}^{{\prime} }$). We calculate the average of these angles, approximating the angle of the current sheet about point ($i,j$) as

$$\begin{eqnarray}&&{\theta }_{\mathrm{mean}}=\displaystyle \frac{1}{N}\displaystyle \sum _{{({i}^{{\prime} },{j}^{{\prime} })}_{n=1}}^{{({i}^{{\prime} },{j}^{{\prime} })}_{N}}\arctan \left[\displaystyle \frac{z({j}_{n}^{{\prime} })-z(j)}{r({i}_{n}^{{\prime} })-r(i)}\right].\end{eqnarray} \tag{ 5 }$$

We then calculate the mean slope,

$$\begin{eqnarray}&&{m}_{\mathrm{mean}}=\tan {\theta }_{\mathrm{mean}},\end{eqnarray} \tag{ 6 }$$

and take the direction perpendicular to this slope as

$$\begin{eqnarray}&&{m}_{\perp }=-\displaystyle \frac{1}{{m}_{\mathrm{mean}}}.\end{eqnarray} \tag{ 7 }$$

We sample the plasma properties at some distance along the normal where the toroidal magnetic field has reached its asymptotic value. We approximate this location by scanning along the normal direction, given by Equation (7), until the field profile flattens out. We consider the field sufficiently flat when the fractional change in magnetic field from one computational cell along the normal to the next is less than three percent, averaged over two adjacent cells.⁵

As an illustrative example, Figure 3 shows the magnetic field and current density profile along the normal of a current sheet picked out by our algorithm in the SANE simulation (the typical shapes of these profiles generalize to the MAD simulations). We indeed see a Harris-sheet-like structure, with a magnetic field profile that passes through 0 and asymptotes to a fixed value at a distance ≈0.2–0.3 GMc⁻² away from the center; the current density has the expected maximum associated with the steep gradient in magnetic field. To illustrate the variety of current sheets and show their typical length scales and field profiles, Figure 4 shows a sample of current sheets identified in different snapshots and locations of the SANE simulation. The magnetic field profiles again follow structures reminiscent of Harris sheets, with the magnetic field passing in a linear fashion through 0 and reaching an asymptotic value at a distance that is typically 0.2 to 0.6 GMc⁻² away from the center of the sheet. We find that while the current sheets are often approximately symmetric (approaching similar magnetic field strengths to either side of the current sheet), there are also cases of non-symmetric current sheets, where the asymptotic magnetic field strength differs between the two sides of the sheet. Studies of magnetic reconnection almost always employ symmetric current sheets, but the asymmetry in magnetic field profile could influence the outcome of reconnection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Having established the frequent occurrence and geometry of potential reconnection regions, our second goal is to investigate the properties of current sheets in time-dependent simulations of accretion flows and characterize the parameters relevant to non-thermal particle acceleration to inform further PIC studies. We ultimately wish to determine the role of magnetic reconnection in contributing to the multiwavelength variability of low-luminosity accretion flows.

Iterating through timesteps in our simulations, we find the current sheets and, at every point in each sheet, determine the asymptotic values of σ and β as well as the guide field strength at the center of the sheet for both our SANE and MAD simulations, as described below.

4.1. Properties of SANE Current Sheets

For the SANE simulation, in the regions where reconnection may occur, the magnetization σ ranges from 10⁻⁴ to 1, while the plasma-β ranges from 0.1 to 10³, as shown in Figure 5.

Zoom In Zoom Out Reset image size

The anticorrelation evident in Figure 5 (see also Figure 7) occurs because the magnetization parameter σ scales as ${B}^{2}/n$, while the plasma-β scales as ${({B}^{2}/n)}^{-1}$. The spread arises because the plasma-β also depends on the plasma temperature. The temperatures in these regions are typical of temperatures throughout the accretion flow, with values of around ${kT}/{m}_{i}{c}^{2}\approx 0.05$.

The most promising subspace of this region for particle acceleration to be efficient is the high-σ, low-β (bottom right) regime, where there is maximal magnetic energy to dissipate into the particles and fairly little gas pressure relative to the magnetic pressure, such that the plasma is magnetically dominated. The inferred ranges of σ and β are interesting for a number of reasons. Studies have only recently begun to investigate this transrelativistic ($\sigma \approx 1$) regime (Werner et al. 2016; Rowan et al. 2017), and the physics of particle acceleration and heating in these conditions are not yet fully understood. While the ions in this regime remain non-relativistic (because σ is of order one), the electrons will likely be accelerated (or heated on average) to highly relativistic speeds, as ${\sigma }_{e}\equiv {B}^{2}/(4\pi {\rho }_{e}{c}^{2})=\sigma {m}_{i}/{m}_{e}\approx {10}^{3}$, which is an estimate of the characteristic electron Lorentz factors expected from reconnection.

Finally, Figure 6 shows a histogram of the relative guide field strengths in the SANE simulation. It is evident that both cases of weak (${B}_{r}/{B}_{\phi }\lt 0.5$) and of no guide fields are of interest for the purposes of these simulations. Even weak guide fields may play an important and potentially adverse role in determining the outcome of particle acceleration in magnetic reconnection and must be explored via PIC simulations in the transrelativistic regime.

Zoom In Zoom Out Reset image size

4.2. Properties of MAD Current Sheets

For the MAD simulation, we find that in the regions of potential reconnection, σ ranges from 10⁻³ to 10, while β ranges from 0.03 to 10³, as shown in Figure 7. This is roughly an order of magnitude higher (lower) than the σ (β) values in the SANE simulation, hinting that particle acceleration may be more efficient in these systems. The higher values of σ in the MAD current sheets reflect the stronger, large-scale torroidal fields expected in MAD simulations, as previously discussed.

Zoom In Zoom Out Reset image size

We show the guide fields in the MAD simulation in Figure 8. In stark contrast with the SANE guide fields, which are weak relative to the reconnecting field, the MAD guide fields are stronger and can be comparable with the reconnecting ones. Previous studies (e.g., Narayan et al. 2012) have shown larger magnetic shear stresses in MAD simulations, which is consistent with the stronger guide fields (${B}_{\theta }$, B_r) we identify. While the typical values of the magnetization σ and the plasma-β are more favorable in terms of particle acceleration in the MAD simulations, the stronger guide fields may alter the outcome of the reconnection event in regard to the particle acceleration.

Zoom In Zoom Out Reset image size

We finally examine the time variability of energy available to reconnection throughout the accretion flow. One motivation for this is to assess whether reconnection events can contribute substantially to the high-energy variability of low-luminosity accretion flows, as has been extensively observed in the case of Sgr A*. We integrate the magnetic energy density, ${B}^{2}/8\pi$, over the reconnecting volume bounded by the surfaces defined by the asymptotic magnetic field location and obtain in this way the total magnetic energy in the reconnection regions throughout the flow. We plot the results of this in Figure 9 for both the SANE and MAD simulations.

Zoom In Zoom Out Reset image size

We see that the turbulent nature of the accretion flow very often leads to the formation of transient current sheets that result in a highly time-varying magnetic energy being available to reconnection. This indicates that magnetic reconnection likely is a significant contributor to the variability of such systems. The SANE model produces persistent variability due to the high levels of turbulence in the disk. The MAD system has less turbulence, hence fewer variations, but the higher degree of magnetization means when current sheets develop, they typically have more magnetic energy associated with them. For this reason, we find that the MAD simulation is characterized by fewer but stronger variations in the magnetic energy available to reconnection.

Figure 10 shows histograms of σ for the current sheets of both SANE (left) and MAD (right) simulations, shown at a quiescent and flaring timestep. We see that during times of increased magnetic energy, there tends to be an associated population of current sheets with higher-than-typical magnetizations (around $\sigma \approx 1$ for the SANE case, and $\sigma \approx 10$ for MAD). This has significant implications for electron acceleration, because the typical energy of non-thermal electrons, as well as the efficiency of electron acceleration both scale with σ.

Zoom In Zoom Out Reset image size

Including the variability properties of non-thermal electrons that are accelerated in these current sheets will likely significantly alter the earlier finding of Chan et al. (2015a), who used models that assumed a purely thermal electron distribution. In that early work, the variability of MAD simulations was characterized by smooth long timescale variations in the flux. It is clear, however, from the present analysis of the MAD simulation that there is the potential for having a sudden injection of non-thermal electrons associated with the spikes in Figure 9, which can then result in corresponding flares in the light curve.

Based on previous studies, we expect some fraction of this magnetic energy to go toward accelerating particles. This acceleration efficiency will, in principle, depend on the flow conditions, such as σ and β, and can be found through PIC simulations. To go from a picture such as the one shown in Figure 9 to the non-thermal particle energy as a function of time, the magnetic energy as a function of time must be combined with the acceleration efficiency as a function of flow parameters, which will likely result in even more dramatic variation of energy on short timescales.

To estimate whether the energy available for reconnection is a plausible explanation for flares of these magnitudes, we calculate the total energy from an average X-ray flare from Sgr A*. Observations of X-ray flares from Sgr A* show typical luminosities from $\sim {10}^{34}$ erg s⁻¹ to $2\times {10}^{35}$ erg s⁻¹ and typical timescales from hundreds of s to 8 ks (Neilsen et al. 2013). With a luminosity of $5\times {10}^{34}$ erg s⁻¹ and a duration of 1000 s, about $5\times {10}^{37}$ erg is being released in a typical flare.

Considering Figure 9, we see for the SANE model that the energy available to reconnection peaks at typical values around 10³⁸ erg, while typical MAD energies are an order of magnitude higher than this. This shows that there is enough energy available to reconnection in these simulations to plausibly account for the observed energy released during these flares. Moreover, the efficiency ξ that determines the fraction of magnetic energy that goes into particle acceleration must be quite high in the case of SANE models, which have typically lower magnetizations and hence less magnetic energy associated with their current sheets.

Ball et al. (2016) characterized the non-thermal particle distribution using η, the fraction of non-thermal to thermal energy densities in the fluid and power-law index, p. We found that significant X-ray flares can occur while satisfying the observed quiescent X-ray constraints for values of $\eta =0.1$ and a conservative power-law index of $p=-3.5$. To connect our present results to these earlier findings, we express η in terms of β and ξ as

$$\begin{eqnarray}&&\eta \equiv \displaystyle \frac{{E}_{\mathrm{nt}}}{{E}_{\mathrm{th}}}=\displaystyle \frac{\xi {B}^{2}}{{nkT}8\pi }=\xi {\beta }^{-1}.\end{eqnarray} \tag{ 8 }$$

We can rewrite this as a constraint on the plasma-β using the η found in our previous study to result in significant X-ray variability, i.e.,

$$\begin{eqnarray}&&0.1\left(\displaystyle \frac{\eta }{0.1}\right)\beta \leqslant \xi .\end{eqnarray} \tag{ 9 }$$

This places a constraint on the plasma-β, given a local ξ, which must be found as a function of flow parameters via PIC simulations. Note that ξ, by definition, cannot be greater than 1, placing a strict upper limit of $\beta =10$ for regions where there is sufficient magnetic energy to accelerate particle to the energies required to generate the flux excursions demonstrated in Ball et al. (2016). More realistically, ξ is likely to be of order 0.1, resulting in an upper limit of $\beta \approx 1$. As shown in Figure 5, we find that we indeed identify many current sheets satisfying this condition. In the MAD simulation (Figure 7), we see a larger number of current sheets extending well below this β threshold (and to higher-σ) indicating that particle acceleration via reconnection may be more efficient in MAD than SANE configurations.

In this paper, we investigated the detailed structure of current sheets and their plasma properties in GRMHD simulations of radiatively inefficient accretion flows. We found that the regimes of plasma parameters relevant to magnetic reconnection have been relatively unexplored in terms of non-thermal particle acceleration. Specifically, we found that the magnetization σ in the vicinity of current sheets in the SANE simulation is of order 10⁻⁴ to 1 , while the plasma-β is of order 0.1 to 10³. Current sheets in the MAD simulation have magnetization σ ranging from 10⁻³ to 10 and plasma-β from 0.03 to 10³. Additionally, we find that in these regions there is a relatively small spread in temperature, leading to a tight correlation between the parameters σ and β. We also characterized the guide fields found in current sheets, which can play a role in governing the details of particle acceleration, and found that the ratio of guide field to reconnecting field strength is typically 0–0.5 for SANE simulations, but can be of order unity in MAD simulations.

GRMHD simulations need to use subgrid models to account for physical effects that cannot be resolved or incorporated in MHD. To employ correctly subgrid models of reconnection, we must improve our understanding of particle acceleration and heating in the parameter space we lay out here.

In addition to characterizing the plasma properties of current sheets, we also calculated the magnetic energy available to reconnection throughout the simulations. We found that the turbulent nature of the accretion flow leads to current sheets of varying characteristics continuously forming and dissipating in the flow. This leads to a highly variable amount of energy available to reconnect and dissipate into heating and particle acceleration and makes magnetic reconnection a promising candidate for contributing to the X-ray variability of Sgr A* and other black holes with similar accretion characteristics. Additionally, we found that there is indeed enough energy available to reconnection around current sheets to account for typical flares observed from Sgr A*. We conclude that if this mechanism is responsible for the X-ray flares, then the acceleration efficiency must be reasonably high for SANE disks and can be lower for the MAD model.

We gratefully acknowledge support for this work from Chandra Award No. TM6-17006X and from NASA TCAN award NNX14AB48G. D.P. acknowledges support from the Radcliffe Institute for Advanced Study at Harvard University. F.O. gratefully acknowledges a fellowship from the John Simon Guggenheim Memorial Foundation in support of this work. L.S. acknowledges support from DoE DE-SC0016542, NASA Fermi NNX-16AR75G, NASA ATP NNX-17AG21G, NSF ACI-1657507, and NSF AST-1716567.

To find regions in our GRMHD simulation where magnetic reconnection would occur, if it was explicitly included, we need to identify regions of high current density. However, GRMHD simulations typically evolve only four quantities, i.e., the magnetic field, density, fluid velocity, and internal energy, and the current density is typically not explicitly computed in the simulation. To understand the structure of the current density throughout the flow, we calculate it from the electromagnetic tensor,

$$\begin{eqnarray}{F}^{\mu \nu }=\left(\begin{array}{cccc}0 & -{E}_{1} & -{E}_{2} & -{E}_{3}\\ {E}_{1} & 0 & {B}_{3} & -{B}_{2}\\ {E}_{2} & -{B}_{3} & 0 & {B}_{1}\\ {E}_{3} & {B}_{2} & -{B}_{1} & 0\end{array}\right)\end{eqnarray} \tag{ 10 }$$

via

$$\begin{eqnarray}&&{J}^{\nu }={{\rm{\nabla }}}_{\mu }{F}^{\nu \mu },\end{eqnarray} \tag{ 11 }$$

where ${{\rm{\nabla }}}_{\mu }$ represents the covariant derivative.

Breaking up the four current into its zeroth and ith components (henceforth, Greek indices run from 0–3, while Latin indices go from 1–3), we can rewrite Equation (11) as

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{j}{F}^{{ij}}-{{\rm{\nabla }}}_{0}{F}^{0i}={J}^{i},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{i}{F}^{0i}={J}^{0}.\end{eqnarray} \tag{ 13 }$$

In the comoving frame, ${F}^{0i}=0$, Equation (13) reads ${J}^{0}=0$. Additionally, the displacement current is 0 (${{\rm{\nabla }}}_{0}{F}^{0i}=0$). We can then write the three current as

$$\begin{eqnarray}&&{J}^{i}={\partial }_{j}{F}^{{ij}}+{{\rm{\Gamma }}}_{j\lambda }^{i}{F}^{\lambda j}+{{\rm{\Gamma }}}_{j\lambda }^{j}{F}^{i\lambda }.\end{eqnarray} \tag{ 14 }$$

Considering the second term on the right hand side, ${{\rm{\Gamma }}}_{j\lambda }^{i}{F}^{\lambda j}$, we see that for every term in the sum, the symmetry of ${{\rm{\Gamma }}}_{{jk}}^{i}$ about its lower two indices and the anti-symmetry of ${F}^{\mu \nu }$ implies perfect cancellation. Hence, Equation (14) reduces to

$$\begin{eqnarray}&&{J}^{i}={\partial }_{j}{F}^{{ij}}+{{\rm{\Gamma }}}_{j\lambda }^{j}{F}^{i\lambda }.\end{eqnarray} \tag{ 15 }$$

We use the Christoffel symbols for the metric of an uncharged Kerr black hole, which in Boyer–Lindquist coordinates is given by

$$\begin{eqnarray}{{ds}}^{2} & = & -\displaystyle \frac{{\rm{\Delta }}}{{\rm{\Sigma }}}{({dt}-a{\sin }^{2}\theta d\phi )}^{2}+\displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}}{{dr}}^{2}+{\rm{\Sigma }}d{\theta }^{2}\\ & & +\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Sigma }}}{(({r}^{2}+{a}^{2})d\phi +{adt})}^{2}\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}{\rm{\Delta }} & = & {r}^{2}-2{Mr}+{a}^{2}\\ {\rm{\Sigma }} & = & {r}^{2}+{a}^{2}{\cos }^{2}\theta .\end{eqnarray} \tag{ 17 }$$

Here, M and a are the mass and angular momentum per unit mass of the black hole, respectively.

For this metric, the relevant nonzero Christoffel terms in Equation (15) are ${{\rm{\Gamma }}}_{11}^{1},{{\rm{\Gamma }}}_{12}^{1},{{\rm{\Gamma }}}_{31}^{3},{{\rm{\Gamma }}}_{32}^{3},{{\rm{\Gamma }}}_{21}^{2},{{\rm{\Gamma }}}_{22}^{2}.$ With these, we can write out the individual components of Jⁱ as

$$\begin{eqnarray}&&{J}^{1}={\partial }_{j}{F}^{1j}+({{\rm{\Gamma }}}_{12}^{1}+{{\rm{\Gamma }}}_{32}^{3}+{{\rm{\Gamma }}}_{22}^{2}){F}^{12},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{J}^{2}={\partial }_{j}{F}^{2j}+({{\rm{\Gamma }}}_{11}^{1}+{{\rm{\Gamma }}}_{31}^{3}+{{\rm{\Gamma }}}_{21}^{2}){F}^{21},\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}{J}^{3} & = & {\partial }_{j}{F}^{3j}+({{\rm{\Gamma }}}_{11}^{1}+{{\rm{\Gamma }}}_{21}^{2}+{{\rm{\Gamma }}}_{31}^{3}){F}^{31}\\ & & +({{\rm{\Gamma }}}_{12}^{1}+{{\rm{\Gamma }}}_{32}^{3}+{{\rm{\Gamma }}}_{22}^{2}){F}^{32}.\end{eqnarray} \tag{ 20 }$$

Using these equations, we solve for the current density in GRMHD simulations at every timestep and identify current sheets where reconnection may take place.