Magnetic Reconnection and Hot Spot Formation in Black Hole Accretion Disks

Here we assume a flat-spacetime Minkowski metric and consider a relativistic ideal gas with an adiabatic index $\hat{\gamma }=4/3$ with an initial uniform pressure p = 10 and rest-mass density ρ = 1. The magnetic field ${\boldsymbol{B}}={\rm{\nabla }}\times {\boldsymbol{A}}$ is obtained from a vector potential ${\boldsymbol{A}}=(0,0,{A}_{z})$ on a 2.5 dimensional Cartesian grid $(x,y,z)$:

$$\begin{eqnarray}&&{A}_{z}=\displaystyle \frac{1}{2}\cos (2x)+\cos (y),\end{eqnarray} \tag{ 1 }$$

and the vortex imposes an initial velocity field ${\boldsymbol{v}}=({v}_{x},{v}_{y},0)$ given by

$$\begin{eqnarray}&&{v}_{{\rm{x}}}=-{v}_{\max }\sin (y),\quad {v}_{{\rm{y}}}={v}_{\max }\sin (x),\end{eqnarray} \tag{ 2 }$$

where ${v}_{\max }=0.99c/\sqrt{2}$ is set to ensure that the initial speed v is limited by the speed of light, where the maximum initial Lorentz factor is ${\rm{\Gamma }}={(1-{v}^{2}/{c}^{2})}^{-1/2}=10$. The electric field is initialized as ${\boldsymbol{E}}=-{\boldsymbol{v}}\times {\boldsymbol{B}}/c$. These settings result in a minimum value of the gas-to-magnetic-pressure ratio $\beta \,=2p/{b}^{2}=10$ and a maximum magnetization $\sigma ={b}^{2}/(\rho {{hc}}^{2})\,\approx 0.05$, where $h=1+4p/(\rho {c}^{2})$ is the specific enthalpy for an ideal gas with adiabatic index $\hat{\gamma }=4/3$ and b² is the magnetic field strength comoving with the fluid. We use a uniform base grid with 128² cells spanning $x,y\in [0,2\pi ]$ with periodic boundary conditions. We then apply subsequently increasing AMR levels to study the convergence of the numerical solution, with each level quadrupling the effective resolution (which is the total number of cells if the highest AMR levels were fully utilized, whereas the actual number of cells is smaller). The mesh refinement is based on second derivatives in the quantities B_z, E_z, and ρ according to the Löhner (1987) scheme. All simulations run until final time t = 10t_c, with ${t}_{{\rm{c}}}=L/c$ being the light-crossing time, where we assume L = 1 as the typical length scale of the system.

We consider a set of uniform and constant resistivities, $\eta \in [0,{10}^{-5},2.5\times {10}^{-5}$, $5\times {10}^{-5},{10}^{-4}$, $5\times {10}^{-4},{10}^{-3},5\times {10}^{-3}]$, chosen such that the corresponding Lundquist numbers $S={Lc}/\eta \in [\infty ,{10}^{5},4\times {10}^{4},2\times {10}^{4}$, ${10}^{4},2\,\times {10}^{3},{10}^{3},2\times {10}^{2}]$ range from infinity (ideal GRMHD) to well above and below the threshold for plasmoid formation ${S}_{\mathrm{th}}\approx {10}^{4}$ (Bhattacharjee et al. 2009; Uzdensky et al. 2010).

The initial magnetic field topology of the vortex is characterized by four alternating X-points and magnetic nulls, resulting in two magnetic islands along horizontal lines y = 0 and y = π. The vortex motion immediately displaces the left island at y = π diagonally upwards and the right island diagonally downwards, resulting in diagonal current sheets. The sheets get compressed, and depending on the Lundquist number, they can become unstable to tearing, demonstrating the break-up in a series of plasmoids (see, e.g., also van der Holst et al. 2008).

For $S\geqslant {10}^{4}$, the current sheets shrink until they become plasmoid-unstable. In Figure 1, we show the rest-mass density ρ for S = 10⁵ at t = 10t_c, where plasmoids are recognized as overdense blobs of plasma in the thin current sheets. The current sheet is characterized by an antiparallel magnetic field configuration resulting in a large out-of-plane component of the current density in the fluid frame,

$$\begin{eqnarray}&&{\boldsymbol{J}}={\eta }^{-1}{\rm{\Gamma }}[{\boldsymbol{E}}+{\boldsymbol{v}}\times {\boldsymbol{B}}/c-({\boldsymbol{E}}\cdot {\boldsymbol{v}}){\boldsymbol{v}}/{c}^{2}],\end{eqnarray} \tag{ 3 }$$

as can be seen in the middle and right panels of Figure 2 for S = 2 × 10⁴ and resolution 8192². For lower Lundquist numbers $S\lt {10}^{4}$, the current sheet does not become thin enough for the plasmoid instability to grow (see, e.g., the left panel of Figure 2 for S = 2 × 10³ and resolution 8192²).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We claim converged results if the evolution of the domain-averaged magnetic energy density ${\bar{B}}^{2}\equiv {\iint }_{V}{B}^{2}{dxdy}/{\iint }_{V}{dxdy}$ no longer changes between successive doubling of the number of grid cells per direction over the full domain V (see Figure 3). The plasmoid-unstable cases ($S\geqslant {10}^{4}$) are the most demanding, and effective resolutions of 8192² cells are needed to resolve the current sheets by at least 10 cells over their width. For our fiducial case of S = 2 × 10⁴, there are no visual differences between runs with resolutions of 8192² and 16,384² grid cells (see middle panel of Figure 3). For lower resolutions $\lt {8192}^{2}$, the (average) magnetic energy density is underestimated (particularly after $t=2{t}_{{\rm{c}}}$) and the dissipation is governed by numerical resistivity that is larger than the explicit resistivity $\eta =5\times {10}^{-5}$, affecting the energetics, heating, and plasmoid statistics. We also ensure that the effect of the AMR on the evolution of B² is negligible compared to a (significantly more expensive) run with uniform resolution of 8192² grid cells (see the dashed black lines in the left and middle panels of Figure 3).

Zoom In Zoom Out Reset image size

In Figure 4, we show the fraction of the domain, calculated as the number of $\mathrm{cells}\times {4}^{{n}_{\max }-n}$, with n the AMR level and ${n}_{\max }$ the maximum number of levels, that is covered by the different AMR levels in time for $S=2\times {10}^{4}$. It is clear that the lowest levels, 1 and 2, corresponding to an effective resolution of 128² and 256², respectively, are only important in the first part of the simulation $t\lt 3{t}_{{\rm{c}}}$ (i.e., when the plasmoid instability is not yet activated). In the turbulent final stage of the simulation, less than 25% of the domain is covered by the highest level corresponding to an effective resolution of 8192², and less than 40% by level 6, corresponding to an effective resolution of 4096², resulting in a major speed-up compared to a uniform resolution of 8192². The highest level is only activated around current sheets, as can also be seen in the right panel of Figure 2.

Zoom In Zoom Out Reset image size

We performed an ideal GRMHD (i.e., η = 0) run (see the right panel of Figure 3), showing that in this case, the magnetic energy density evolution, and hence the (numerical) dissipation and heating of the plasma, does not converge even for the highest resolutions considered. Due to the absence of a resistive dissipation scale set by $\eta \gt 0$, the current sheet shrinks to the grid scale and is captured by one grid cell over its width for all resolutions considered here. Notice that the evolution of B² for 4096² and 16,384² grid cells seems comparable for η = 0, whereas the magnetic energy density is clearly lower for the case with 8192² cells. The lowest resolution (1024²) considered here is generally high enough to capture the essential (nonresistive) effects of the magnetorotational instability (MRI; Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991), the accretion dynamics, and the Blandford–Znajek jet launching process (Blandford & Znajek 1977) in global disk simulations according to recent studies in ideal GRMHD (Porth et al. 2019; White et al. 2019). One can see that at resolutions $\lt {1024}^{2}$, there is hardly a difference between the ideal and S = 2 × 10⁴ resistive results (cyan line in middle and right panels of Figure 3), and the linear growth phase of the magnetic energy density (at $t\lesssim 3{t}_{{\rm{c}}}$, i.e., where no plasmoids have formed yet) is accurately captured, confirming that numerical resistivity due to the finite grid dominates the magnetic energy dissipation and that higher resolutions are needed to resolve magnetic reconnection and plasmoid formation. For S = 2 × 10³, the evolution of the current sheet is easier to capture because the violent tearing instability is not triggered and the explicit resistivity is always larger than the numerical resistivity in this case. The resistive scale is already fully resolved for an effective resolution of 1024² cells (see the left panel of Figure 3), and the sheet does not shrink below the threshold for plasmoid formation (see the left panel of Figure 2).

We analyze the reconnection rate for all considered Lundquist numbers in the converged runs by measuring the inflow velocity into the current sheet. In order to determine the inflow speed, we use the upstream ${\boldsymbol{E}}\times {\boldsymbol{B}}$ velocity such that ${v}_{\mathrm{up}}/c={E}_{z,\mathrm{up}}/{B}_{\mathrm{up}}$. We take five slices across the current sheet (ensuring that the slice does not cut across a plasmoid) at t = 10t_c for all cases with resolutions of 16,384² such that even the thinnest sheet for S = 10⁵ is resolved. We then obtain a profile of ${B}_{\mathrm{up}}$ by projecting the magnetic field along the current sheet and find the location where the profile becomes flat (see the bottom panel of Figure 5 for S = 2 × 10⁴). We measure both the upstream electric and magnetic fields by averaging over five locations in the upstream on the slices on both sides of the sheet. We account for the bulk flow velocity (${v}_{\mathrm{bulk}}\ll c$ such that ${{\rm{\Gamma }}}_{\mathrm{bulk}}=1$) of the vortex by defining the speed on the left of the sheet as ${v}_{\mathrm{up},\mathrm{left}}/c=({v}_{\mathrm{bulk}}+{v}_{\mathrm{in}}$)/c and on the right of the sheet as ${v}_{\mathrm{up},\mathrm{right}}/c=({v}_{\mathrm{bulk}}-{v}_{\mathrm{in}})/c$, such that ${v}_{\mathrm{rec}}/c=({v}_{\mathrm{up},\mathrm{left}}-{v}_{\mathrm{up},\mathrm{right}})/2c$. Assuming that the outflow ${v}_{\mathrm{out}}={v}_{A}\approx c$, this directly yields the reconnection rate ${v}_{\mathrm{rec}}/c={v}_{\mathrm{in}}/{v}_{\mathrm{out}}={v}_{\mathrm{in}}/c$. We confirmed that locally, in the upstream region of the current sheet, $\sigma ={b}^{2}/(\rho {{hc}}^{2})\approx 8$ such that the Alfvén speed around the sheet is ${v}_{{\rm{A}}}\,=c{(\sigma /(\sigma +1))}^{1/2}\approx c$ and that the reconnection is relativistic. In the top panel of Figure 5, we observe a Sweet–Parker scaling ${v}_{\mathrm{rec}}/c\sim {S}^{-1/2}$ for $S\lt {10}^{-4}$ (indicated by the blue circles and the dashed black line) and plasmoid-dominated "fast" reconnection with a rate independent of the Lundquist number of ${v}_{\mathrm{rec}}/c\approx 0.01$ for $S\geqslant {10}^{4}$ (indicated by the red circles and the dotted black line).

Zoom In Zoom Out Reset image size

We consider both MAD and SANE magnetic field configurations around a Kerr black hole, for three representative values of the dimensionless spin parameter $a\in [-0.9375;0;0.9375]$, corresponding to near-extremal retrograde Kerr, Schwarzschild, and near-extremal prograde Kerr black holes, respectively. We use geometrized units with gravitational constant, black hole mass, and speed of light G = M = c = 1, such that length scales are normalized to the gravitational radius ${r}_{g}={GM}/{c}^{2}$ and times are given in units of r_g/c. We employ spherical Kerr–Schild coordinates, where r is the radial coordinate, θ and ϕ are the poloidal and toroidal angular coordinates, respectively, and t is the temporal coordinate. We start our simulations from a torus in hydrodynamic equilibrium (Fishbone & Moncrief 1976) threaded by a single weak poloidal magnetic field loop, defined by the vector potential

$$\begin{eqnarray}&&{A}_{\phi }\propto \max (q,0),\end{eqnarray} \tag{ 4 }$$

where the function q is set to obtain a large torus resulting in a MAD state:

$$\begin{eqnarray}&&q=\displaystyle \frac{\rho }{{\rho }_{\max }}{\left(\displaystyle \frac{r}{{r}_{\mathrm{in}}}\right)}^{3}{\sin }^{3}\theta \exp \left(-\displaystyle \frac{r}{400}\right)-0.2,\end{eqnarray} \tag{ 5 }$$

with an inner radius ${r}_{\mathrm{in}}=20{r}_{g}$ and the density maximum ρ_max is located at ${r}_{\max }=41{r}_{g}$ for a = 0.9375, ${r}_{\max }=42{r}_{g}$ for $a=-0.9375$, and ${r}_{\max }=41.4{r}_{g}$ for a = 0 to start from an initial torus of similar size. For a SANE state, we set a smaller torus with

$$\begin{eqnarray}&&q=\displaystyle \frac{\rho }{{\rho }_{\max }}-0.2,\end{eqnarray} \tag{ 6 }$$

where the inner radius ${r}_{\mathrm{in}}=6{r}_{g}$ and the density maximum ρ_max is located at ${r}_{\max }=12{r}_{g}$ for a = 0.9375, ${r}_{\max }=22{r}_{g}$ for a = −0.9375, and ${r}_{\max }=15{r}_{g}$ for a = 0 to start from an initial torus of similar size. In both cases, the magnetic field strength is set such that $\beta =2{p}_{\max }/{b}_{\max }^{2}=100$. The plasma β and the magnetization $\sigma ={b}^{2}/(\rho {c}^{2})$ for a cold (i.e., $p\ll \rho$) plasma are defined using the magnetic field strength b² comoving with the fluid (see Ripperda et al. 2019a for a definition). We set an atmospheric rest-mass density and pressure as ${\rho }_{\mathrm{atm}}={\rho }_{\min }{r}^{-3/2}$ and ${p}_{\mathrm{atm}}={p}_{\min }{r}^{-5/2}$, where ${\rho }_{\min }={10}^{-4}$ and ${p}_{\min }={10}^{-6}/3$. We apply floors on rest-mass density, pressure, and Lorentz factor ${\rm{\Gamma }}\lt {{\rm{\Gamma }}}_{\max }=20$ such that the magnetization $\sigma ={b}^{2}/(\rho {c}^{2})\lt {\sigma }_{\max }=100$ and ${\beta }^{-1}={b}^{2}/2p\lt {\beta }_{\max }^{-1}\,={(10{\sigma }_{\max })}^{(\hat{\gamma }-1)}$. We adopt an equation of state for a relativistic ideal gas with an adiabatic index of $\hat{\gamma }=4/3$. The equilibrium fluid pressure is perturbed to trigger the MRI as $p={p}_{\mathrm{eq}}(1+{X}_{p})$ with a random variable uniformly distributed between ${X}_{p}\in [-0.02;0.02]$.

We model dissipation in the accretion flow by assuming a small, uniform, and constant resistivity $\eta =5\times {10}^{-5}$ in the set of GRRMHD equations (Ripperda et al. 2019a), resulting in a large Lundquist number $S={Lc}/\eta \sim 2\times {10}^{5}$ that is well above the plasmoid threshold, where we assume the typical length scale to be the length of a current sheet $L\,={L}_{\mathrm{sheet}}\sim { \mathcal O }(10{r}_{g})$, and the speed of light c as the typical velocity. A small resistivity allows us to capture both the global near-ideal accretion dynamics (Ripperda et al. 2019a) and the fast reconnection resulting in plasmoid formation in localized thin current sheets (Ripperda et al. 2019c).

We show in Figure 6 that both MAD and SANE configurations in ideal and resistive runs reach a quasi-steady state after approximately $2500{r}_{g}/c$ and $500{r}_{g}/c$, respectively, as indicated by the magnetic flux through the horizon:

$$\begin{eqnarray}&&\dot{{\rm{\Phi }}}:= \displaystyle \frac{1}{2}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }| {B}^{r}| \sqrt{-\gamma }d\theta d\phi ,\end{eqnarray} \tag{ 7 }$$

where γ is the determinant of the spatial part of the metric.

Zoom In Zoom Out Reset image size

Figure 6 shows the magnetic flux for the highest resolution considered, 6144 × 3072 cells, compared to a resolution of 3072 × 1536, confirming that the quasi-steady-state phase of accretion is independent of the resolution in our simulations. We confirm convergence of the thinning process of the current sheets by restarting from the quasi-steady state and increasing the resolution up to 12,288 × 6144, and measure the cells per thickness of the sheets. In all considered cases, the sheets are captured by eight or more cells over their widths, which is discussed in Section 3.3.

Current sheets are expected to form above the disk or in the jet's sheath where magnetic flux tubes are twisted by global shearing motion. These structures can inflate while they accrete onto the black hole and get thinner, after which their magnetic energy is dissipated close to the event horizon through reconnection. The magnetic energy is released into heat and bulk motion of the plasmoids that can either fall into the black hole or get ejected. This process has been studied in the force-free paradigm, assuming the plasma to be infinitely magnetized (Parfrey et al. 2014; Yuan et al. 2019a, 2019b; Mahlmann et al. 2020). Additionally, when the net magnetic flux in the accretion disk is relatively small, turbulence resulting from the MRI can produce magnetic fields with alternating polarities prone to reconnection (Davis et al. 2010; Zhu & Stone 2018). MHD turbulence is known to intermittently form large plasmoid-unstable current sheets (Zhdankin et al. 2013, 2017; Dong et al. 2018), and the plasmoid instability can significantly modify the turbulent MHD cascade at relatively small scales and high Lundquist numbers $S\sim {10}^{6}$ (Boldyrev & Loureiro 2017; Comisso et al. 2018; Dong et al. 2018).

In Figure 7, we observe both processes in a SANE configuration and detect current sheets in the disk and along the jet's sheath, indicated by a small ${\beta }^{-1}$ and antiparallel field lines. In the leftmost two panels at $t=1460{r}_{g}/c$ and $t=1480{r}_{g}/c$, a large flux tube falls onto the black hole in the bottom-left half at approximately $x={r}_{\mathrm{KS}}\cos \theta \approx 8{r}_{g}$, $y={r}_{\mathrm{KS}}\sin \theta \approx -8{r}_{g}$. In the third panel, at $t=1540{r}_{g}/c$, the flux tube has both its footpoints attached to the black hole, it opened up after it inflated and became thin enough to be tearing unstable and form multiple plasmoids (Parfrey et al. 2014 observe a similar process). In the fourth panel, at $t=1580{r}_{g}/c$, the plasmoids that are advected away from the black hole along the jet's sheath have formed a large structure at $x\approx 6{r}_{g},y\approx -10{r}_{g}$ through coalescence. At $x\approx 9{r}_{g},y\,\approx 17{r}_{g}$, a similar process occurs in the third panel and a large plasmoid has formed through mergers of multiple smaller plasmoids that can be seen in the left panels. In the first and second panels, one can detect a flux tube with one footpoint connected to the black hole, which is then twisted by the shear flow, becomes thinner and inflates, at $x\approx 3{r}_{g},y\approx 2{r}_{g}$, after which it ejects plasmoids into the accretion disk in the fourth panel at $x\approx 6{r}_{g},y\approx 2{r}_{g}$. A similar process occurs in the first panel for a loop with one footpoint on the black hole at $x\approx 7{r}_{g},y\approx -2{r}_{g}$. The large-scale plasmoids form on a timescale of $\sim { \mathcal O }(100{r}_{g}/c)$, growing to circular objects with a radius $\gtrsim { \mathcal O }({r}_{g})$, indicating a reconnection rate of $\approx 0.01c$. These structures eventually break up and lose coherence due to interaction with the ambient flow.

Zoom In Zoom Out Reset image size

In Figure 8 we show the magnetization $\sigma ={b}^{2}/(\rho {c}^{2})$, temperature $T=p/\rho$, magnetic dissipation $| {E}^{i}{J}_{i}|$, and the current density magnitude $J=\sqrt{{J}^{i}{J}_{i}}$, where Eⁱ and Jⁱ are the Eulerian electric field and current density (see Ripperda et al. 2019a for definitions) at $t=1540{r}_{g}/c$ (the third panel in Figure 7). We mask the regions with $\sigma \geqslant 5$ where the flow dynamics is dominated by the density and pressure floors. In the top-right panel, one can detect ubiquitous current sheets along the jet's sheath and inside the disk within $10{r}_{g}$, indicated by a strong current density. The current sheets are tearing unstable, and plasmoids form close to the event horizon after which they either fall into the black hole or are advected along the jet's sheath, merge with other plasmoids, and grow into larger objects with a size of the order of a few r_g. The magnetization in the top-left panel shows that current sheets along the jet's sheath reconnect in a highly relativistic regime (i.e., $\sigma \gg 1$), whereas in the disk reconnection occurs in the transrelativistic regime $\sigma \lesssim 1$. Reconnection in the highly magnetized plasma surrounding the jet's sheath and the current sheets in the equatorial plane efficiently heats plasmoids to relativistic temperatures $T=p/\rho \sim 1$ as can be seen in the bottom-right panel. These hot plasmoids are advected along the jet's sheath or into the accretion disk. Heating of plasmoids that form inside the disk due to reconnection induced by MRI turbulence is less efficient due to the low magnetization in the disk. Instead, hot plasmoids observed at $x\approx 6{r}_{g},y\approx \pm 2{r}_{g}$ were energized in the inner region close to the event horizon and ejected into the disk. The plasmoids are mainly heated close to the event horizon by ohmic heating as can be seen from the bottom-left panel. The strong parallel electric field indicated by $| {E}^{i}{J}_{i}|$ in the bottom-left panel can potentially accelerate particles to nonthermal energies.

Zoom In Zoom Out Reset image size

Figure 9 shows a zoom into the current sheet along the jet's sheath with several interacting plasmoids, as visible in Figure 8. The current sheets are captured by typically 10 cells along their widths, as can be seen from the overplotted grid-block structure (white rectangles), where each block consists of 64 × 32 cells.

Zoom In Zoom Out Reset image size

In a MAD configuration, the accretion disk is threaded by a strong magnetic flux suppressing the MRI and forming a more powerful magnetized jet than in the SANE case. In Figure 10, we show the typical evolution of a current sheet indicated by a small ${\beta }^{-1}$ in the highly magnetized equatorial region close to the black hole at representative times $t=[2941,2971,2988,3009]{r}_{g}/c$ (from left to right). In the left panel, a magnetic flux tube falls onto the black hole at $x\approx 1{r}_{g},y\approx 3{r}_{g}$, and in the second panel, the tube connects its two footpoints to the black hole. A near-equatorial current sheet forms and produces multiple small plasmoids within $5{r}_{g}$ from the event horizon. The plasmoids that escape the gravitational pull merge into a large structure in the third panel at $x\approx 6{r}_{g},y\approx 2.5{r}_{g}$. The large-scale plasmoid grows to a circular structure with a radius $\gtrsim 1{r}_{g}$ and escapes along the jet's sheath at $x\approx 9{r}_{g},y\approx 6{r}_{g}$ in the fourth panel. The process of formation to ejection takes place on a timescale of $\sim 70{r}_{g}/c$. In the fourth panel, the process restarts with an infalling flux tube at $x\approx 9{r}_{g},y\approx 0{r}_{g}$.

Zoom In Zoom Out Reset image size

In Figure 11, we show the magnetization, current density, ohmic heating, and temperature for the MAD state, where we again mask the regions with $\sigma \geqslant 5$. Due to the higher magnetization (see the top-left panel) surrounding the equatorial current sheets (see the top-right panel) in the MAD state, the plasmoids are heated to relativistic temperatures $T\sim 10$ (bottom-right panel), an order of magnitude higher than in the SANE case. The current density in the sheets is significantly higher than in the SANE case. The plasmoids are heated through ohmic heating close to the event horizon, as can be seen from the bottom-left panel.

Zoom In Zoom Out Reset image size

We calculate the reconnection rate in a similar way to the Orszag–Tang vortex for both MAD and SANE configurations. We first transform the Eulerian electric and magnetic fields into a local inertial frame (see, e.g., White et al. 2016) to apply the standard reconnection analysis. We project the fields in the flat frame along the direction parallel to the current layer to determine the upstream geometry, and a typical Harris-type sheet structure is found in Figure 12 both for the magnetic field and the current density magnitude J. All three magnetic field components switch sign in the current sheets, indicating that zero-guide-field reconnection occurs in both MAD and SANE cases.

Zoom In Zoom Out Reset image size

In the local inertial frame, we determine the inflow speed from the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ velocity that we project along the direction perpendicular to the current sheet, and then calculate the reconnection rate as ${v}_{\mathrm{rec}}/c=({v}_{\mathrm{up},\mathrm{left}}-{v}_{\mathrm{up},\mathrm{right}})/2c$. In both MAD and SANE configurations, we select 10 current sheets at different times during the quasi-steady-state phase of accretion and consistently find a reconnection rate between 0.01c and 0.03c. This finding is in accordance with analytic resistive MHD predictions for plasmoid-dominated reconnection in isolated current sheets (Bhattacharjee et al. 2009; Uzdensky et al. 2010). Note that the actual Lundquist number is approximately $S={L}_{\mathrm{sheet}}c/\eta \gtrsim { \mathcal O }({10}^{5})$, as all current sheets have a typical length scale of $\sim { \mathcal O }(5\mbox{--}20{r}_{g})$, confirming that reconnection occurs in the plasmoid-dominated regime as $S\gg {10}^{4}$.