<figure><img src="image_1.png" /></figure>

HAL
open science

# Plasmoid identification and statistics in two-dimensional Harris sheet and GRMHD simulations

Jesse Vos, Hector Olivares, Benoit Cerutti, Monika Moscibrodzka

<figure><img src="image_2.png" /></figure>

## To cite this version:

Jesse Vos, Hector Olivares, Benoit Cerutti, Monika Moscibrodzka. Plasmoid identification and statistics in two-dimensional Harris sheet and GRMHD simulations. Monthly Notices of the Royal Astronomical Society, 2024, 531 (1), pp.1554-1577. 10.1093/mnras/stae1046. hal-04229996

HAL Id: hal-04229996

https://hal.science/hal-04229996v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# Plasmoid identification and statistics in two-dimensional Harris sheet and GRMHD simulations

J. T. Vos<sup>1,1★</sup> H. Olivares<sup>1</sup> B. Cerutti<sup>1,2</sup> and M. Mościbrodzka<sup>1,1</sup>

Accepted 2024 April 10. Received 2024 March 6; in original form 2023 September 7

## ABSTRACT

Magnetic reconnection is a ubiquitous phenomenon for magnetized plasma and leads to the rapid reconfiguration of magnetic field lines. During reconnection events, plasma is heated and accelerated until the magnetic field lines enclose and capture the plasma within a circular configuration. These so-called plasmoids could therefore observationally manifest themselves as hotspots, which are associated with flaring behaviour in supermassive black hole systems, such as Sagittarius A*. We have developed a novel algorithm for identifying plasmoid structures, which incorporates watershed and custom closed contouring steps. From the identified structures, we determine the plasmoids’ plasma characteristics and energetics in magnetohydrodynamical simulations. The algorithm’s performance is showcased for a high-resolution suite of axisymmetric ideal and resistive magnetohydrodynamical simulations of turbulent accretion discs surrounding a supermassive black hole. For validation purposes, we also evaluate several Harris current sheets that are well-investigated in the literature. We recover the characteristic power-law distribution of plasmoid sizes for both the black hole and Harris sheet simulations. This indicates that while the dynamics are vastly different, with different dominant plasma instabilities, the plasmoid creation characteristics are similar. Plasmoid formation rates for resistive general relativistic magnetohydrodynamical simulations are significantly higher than for their ideal counterpart. Moreover, the largest identified plasmoids are consistent with sizes typically assumed for semi-analytical interpretation of observations. We recover a positive correlation between the plasmoid formation rate and decreases in black-hole-horizon-penetrating magnetic flux, during which the accretion flow is temporarily halted. These results demonstrate the efficacy of the newly developed algorithm which has enabled an extensive quantitative analysis of plasmoid formation in black hole accretion simulations.

**Key words:** accretion, accretion discs – black hole physics – magnetic reconnection – MHD – methods: numerical.

## 1 INTRODUCTION

Flaring events, at the X-ray and infrared wavelengths, are known to occur on a daily basis for the supermassive black hole (SMBH) at the centre of the Milky Way, Sagittarius A* (hereafter Sgr A*, Baganoff et al. 2001; Genzel et al. 2003; Eckart et al. 2004; Witzel et al. 2021). The SMBH has an estimated mass of $M \approx 4 \times 10^6$ M$_\odot$ and lies at a distance of $D \approx 8$ kpc as was established by long-term monitoring programs of the source and dynamics of orbiting stars (Ghez et al. 2008; Gillessen et al. 2009a, b, 2017; Gravity Collaboration 2018a, 2019; Do et al. 2019). The daily flares are most apparent at near-infrared (NIR; Gravity Collaboration 2018b) and X-ray wavelengths (e.g. Baganoff et al. 2001; Porquet et al. 2003) as the flux density increases by several orders of magnitude. At sub-mm/mm wavelengths, Sgr A* is known to be a stochastically ($O(10$ per cent) over hours) variable, which is associated with the stereotypical *quiescent* accretion state. While flares at NIR/X-ray wavelengths correspond to significant increases in flux, ‘flaring’ events at mm-wavelengths are typically hard to disentangle from the background variability (EHTC et al. 2022a; Wielgus et al. 2022b). Recently, it was shown

that mm-wavelength light curves observed with the Atacama Large Millimeter/submillimeter Array suggest orbital motion of a hotspot quickly after an X-ray flare (Wielgus et al. 2022a). The physical mechanism that causes these flares is currently not well-understood, but a number of working theories associate them with strongly magnetized anisotropies and locations of particle acceleration in the accretion flow (Broderick & Loeb 2005, 2006; Dexter et al. 2020; Gravity Collaboration 2020b; Porth et al. 2021; Ripperda et al. 2022; Vos, Mościbrodzka & Wielgus 2022; Vos et al. 2023).

One such scenario that may explain these flares and the creation of hotspots is the formation of plasmoids as part of a magnetic reconnection event (e.g. Ripperda, Bacchini & Philippov 2020; Ripperda et al. 2022; El Mellah, Cerutti & Crinquand 2023). This is a phenomenon that occurs in a vast number of astrophysical sources, including pulsar wind nebulae, magnetars, black hole and neutron star magnetospheres, or relativistic jets of active galactic nuclei (Kagan et al. 2015). Magnetic reconnection (Uzdensky 2022, for a review) can broadly be thought of as a rapid reconfiguration of the magnetic field geometry at the interface of opposite polarity magnetic fields that results in the formation of a magnetic island with a typical circular magnetic field morphology. After the closing of the magnetic field lines, plasma is trapped within the magnetic field structure, creating what is known as a plasmoid. The reconfiguration is often

<footnote n="1"><sup>1</sup>Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands</footnote>

<footnote n="2"><sup>2</sup>Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote>* E-mail: jt.vos@astro.ru.nl</footnote>
accompanied by particle acceleration to high (non-thermal) energies (Werner et al. 2017) – effectively converting electromagnetic energy into particle kinetic energy (thermal and non-thermal). A theoretical description for the large-scale dynamics of magnetic reconnection in idealized configurations was established by Sweet (1958) and Parker (1957). This picture is, however, too simplistic for our purposes as it does not deal directly with the formation of plasmoid. To model the plasmoid-unstable regime, one has to adopt a numerical approach via particle-in-cell (PIC) or magnetohydrodynamical (MHD) simulations. The fact that accretion flows near Sgr A*, and in low-luminosity active galactic nuclei in general, are thought to be low in density with corresponding low accretion rates (e.g. EHTC et al. 2022a, b), indicates that the plasma in these environments behaves in a largely collisionless manner. General relativistic (GR) PIC methods are the most accurate numerical model for plasmas in this regime.

Fully kinetic PIC methods have been used to model collisionless ion-electron, electron-positron (pair), or ion-pair plasmas (Kagan et al. 2015, for a review). These methods are considered first principle as they naturally impose both a spatial (skin depth $c/\omega_p$) and temporal ($\omega_p^{-1}$) scale via the plasma oscillation frequency $\omega_p = \sqrt{4\pi nq^2/m}$, where n, q, m are the particle number density, charge, and mass, respectively. While MHD methods only describe the plasma’s bulk motion and characteristics, PIC methods track the velocities, trajectories, and energies of individual particles. Collisionless plasma studies have been conducted to investigate various physical scenarios; isolated (Harris) current sheets in 2D (Zenitani & Hoshino 2001; Cerutti et al. 2012; Nalewajko et al. 2015; Kagan, Nakar & Piran 2016; Sironi, Giannios & Petropoulou 2016; Petropoulou & Sironi 2018) and 3D (Zenitani & Hoshino 2008; Cerutti et al. 2014; Sironi & Spitkovsky 2014; Guo et al. 2016; Werner & Uzdensky 2017), configurations investigating magnetic turbulence (Comisso & Sironi 2019; Bacchini et al. 2022; Borgogno et al. 2022), and (general-relativistic) accretion simulations describing plasma within the magnetosphere of compact objects (for black holes; Parfrey, Philippov & Cerutti 2019; Crinquand et al. 2020, 2021; El Mellah et al. 2022, 2023, or neutron stars; Chen & Beloborodov 2014; Cerutti et al. 2015; Philippov & Spitkovsky 2018; Guépin, Cerutti & Kotera 2020). Although PIC methods are instrumental in, e.g. understanding the origin of non-thermal emission, they remain confined to microscopic plasma scales, which makes interpretation at astrophysically large scales difficult.

General relativistic (ideal) magnetohydrodynamical (GRMHD) methods have been extensively and successfully used to describe the macroscopic picture of accretion onto SMBHs (for M87*; EHTC et al. 2019a, b, 2021, for Sgr A*; EHTC et al. 2022a, b). Resistive (GR)MHD does give a scale to the current sheet and makes it resolvable (Ripperda et al. 2019a, and reference therein) by means of imposing a constant resistivity ($\eta$) in the simulations. While the physical resistivity is likely spatially and temporally variable (and would ideally be constrained with a dedicated PIC study), a uniform scalar resistivity already helps to consistently capture the dynamics associated with magnetic reconnection in the accretion flow. Even though not physically or numerically well-constrained, we point out that magnetic reconnection and plasmoid formation does occur in ideal GRMHD, where numerical limits effectively impose the minimally achievable resistivity.

In this work, we investigate plasmoid formation from fast relativistic reconnection for plasmoid-forming astrophysical plasma in both ideal and resistive GRMHD. To be able to assess the plasmoid formation dynamics, we need to address another, equally important aspect which is that plasmoid structures are difficult to isolate from their surroundings. Therefore, we have developed a novel analysis

algorithm for detecting them. It deviates significantly from plasmoid-finding methods employed previously for GRMHD simulations (Nathanail et al. 2020). Plasmoid-identification studies have been amply applied in more idealized plasma-physical configurations both for MHD (Huang & Bhattacharjee 2012; Loureiro et al. 2012) and PIC studies (Sironi et al. 2016; Petropoulou & Sironi 2018; Hakobyan, Philippov & Spitkovsky 2019; Banesh et al. 2020; Hakobyan et al. 2021; Winarto & Kunz 2022). For a MHD description, one does not have the luxury of individual particle trajectories as typically acquired in PIC simulations. We therefore apply our analysis fully in post-processing which gives it considerable flexibility. With our methodology, we investigate the differences in occurrence rate, morphology, size, and typical plasma parameters of plasmoids in both ideal and resistive GRMHD for a newly created suite of axisymmetric (2.5D) simulations with exquisite resolution. To showcase the validity and high fidelity of the algorithm, we also apply it to a set of Harris current sheet simulations that are equally well-resolved. The application of our methodology to the GRMHD accretion disc simulations also enables a study of the connection between plasmoids formation and the occurrence of flares in the Sgr A* system.

The paper is structured as follows. An in-depth description of the two types of simulations (i.e. Harris current sheet and GRMHD accretion disc setups) is outlined in Section 2. The plasmoid identification procedure is described in detail in Section 3. The results and their interpretation are presented in Section 4. The discussion and conclusion can be found in Sections 5 and 6.

## 2 METHODS

In the following sections, we outline the two simulation classes we investigate, which are the Harris current sheet (Section 2.2) and GRMHD (Section 2.3) simulations. We also outline how we acquire the energetics and surface-averaged quantities (Section 2.4), which are used extensively throughout this work.

### 2.1 Relativistic MHD primer: ideal and resistive

Both the Harris sheet and BH accretion disc are simulated within the framework of the Black Hole Accretion Code (BHAC, Porth et al. 2017; Olivares et al. 2019). BHAC solves the ideal or resistive (Ripperda et al. 2019a) MHD equations in stationary and arbitrary spacetimes. The covariant ideal MHD equations are written as

$$ \nabla_{\mu}(\rho u^{\mu}) = 0, \qquad (1) $$

$$ \nabla_{\mu} T^{\mu\nu} = 0, \qquad (2) $$

$$ \nabla_{\mu} * F^{\mu\nu} = 0, \qquad (3) $$

where $\nabla_\mu$, $\rho$, $u^\mu$, $T^{\mu\nu}$, and $*F^{\mu\nu}$ are the covariant derivative, the rest-mass density, the fluid four-velocity, the energy-momentum tensor (containing both ideal fluid and electromagnetic fields), and the (Hodge) dual of the Faraday tensor, respectively. BHAC is a versatile code that sets the speed of light $c$ to unity and utilizes Lorentz-Heaviside units, which effectively incorporates the $\sqrt{4\pi}$ factors into the electromagnetic quantities.

In this work, we utilize both ideal and resistive MHD. The main difference between both these approaches is the way they handle the evolution of the electric field, which defined as

$$ \mathbf{E} = -\mathbf{v} \times \mathbf{B} + \eta \mathbf{J}. \qquad (4) $$

The resistivity is denoted by $\eta = 1/\sigma_c$ where $\sigma_c$ is the conductivity. While in resistive MHD the electric field ($\mathbf{E}$) includes an explicit
calculation of the resistive Ohm's law to get an expression for J, in ideal MHD it is inferred directly from the magnetic field (via $E = -v \times B$, also known as the frozen-in condition). In the ideal MHD limit, one effectively assumes the plasma to be perfectly conducting (i.e. $\sigma_c \to \infty \Rightarrow \eta \to 0$), which gives a macroscopically valid approximation in large parts of the accretion disc domain. At the basis of plasmoid formation lies the tearing of an induced current sheet via the similarly named tearing- or plasmoid-instability, which is ill-defined in ideal MHD. More specifically, the resistivity $\eta$ is not exactly zero in the ideal case (except for infinite resolution), but rather determined numerically by the underlying resolution (or cell size $\Delta x$) which implies that $\eta_{ide} \propto \Delta x^k$ with $k \approx 2$ depending on the accuracy of the fluid evolution scheme (Ripperda et al. 2022). The physical interpretation of the resistivity $\eta$ is that it acts as a proxy for kinetic effects within the plasma (and has similar properties as a diffusivity).

We investigate plasmoid formation in the fast relativistic plasmoid-dominated reconnection regime. Whether the plasma becomes plasmoid unstable is determined by the Lundquist number $S = L'v_a/\eta$, with typical length of the current sheet $L'$ and the Alfvén velocity $v_a$ (see Section 2.2 for definition). In order to trigger the fast reconnection and tearing- or plasmoid-unstable regime, the Lundquist number needs to satisfy $S > S_{crit}$ where $S_{crit} \sim 10^4$ (Loureiro, Schekochihin & Cowley 2007; Bhattacharjee et al. 2009; Uzdensky, Loureiro & Schekochihin 2010). Note that the Lundquist number is largely determined by the underlying resistivity, which is set as a constant and uniform quantity. For the resistive simulations, which are described in the following sections, we set $\eta = 5 \times 10^{-5} \tau$, where $\tau$ is the arbitrary time unit used in the simulations.<fnref n="1" /> Then, if we estimate probable values of $L' \approx 1$ and $v_a \approx c = 1$, we find $S = 2 \times 10^4$ which lies above the fast reconnection threshold. For the ideal simulations, one may think that as $\eta_{ide}$ is very small, Lundquist number will be sufficiently large. Even though this is the case, the resulting current sheet will always be under-resolved (as it is determined by the underlying resolution) and typically has a width comparable to a singular grid cell (Ripperda et al. 2020). This indicates that the tearing instability is not triggered in the same way as for the resistive simulation and will likely result in differences in plasmoid formation statistics.

## 2.2 Harris sheet configuration

To validate the methodology for a well-known case, we investigate a relativistic 2D Harris sheet in resistive MHD. The implementation is broadly based on what was prescribed for the Geospace Environmental Modeling (GEM) challenge (Birn & Hesse 2001; Birn et al. 2001, also in Goedbloed, Poedts & Keppens 2010). We start with a (wide) rectangular box with periodic boundary conditions on all sides and initialize two current sheets on top of a uniform background density that is scaled with $\rho_0$;

$$ \rho = \rho_0 \left[ \cosh^{-2} \left( \frac{y+L_y/2}{\delta} \right) + \cosh^{-2} \left( \frac{y-L_y/2}{\delta} \right) + f_{bg} \right], \quad (5) $$

where $L_x$, $L_y$, $f_{bg}$, and $\delta$ are the box half-sizes (in the $\hat{x}$ and $\hat{y}$ directions), the background factor, and the sheet half-thickness, respectively. The values of these parameters (and others) are summarized in Table 1.

<table><caption>**Table 1.** The user-defined initial parameters for the Harris sheet simulations that include the model names (acronyms are derived from *Harris-small* and *Harris-big*), density scaling $\rho_0$, layer half-thickness $\delta$, background factor $f_{bg}$, and resolution (with corresponding AMR level). The total dimensions of the box are denoted by $x \in [-L_x, L_x]$, $y \in [-L_y, L_y]$. In addition to the listed parameters, there are several parameters that are constant between (all) the simulations; magnetic field scaling $B_0 = 1$, density scaling $\rho_0 = 1$, resistivity $\eta = 5 \times 10^{-5}$, and adiabatic index $\dot{\gamma} = 13/9$.</caption><thead><tr><th>Name</th><th>δ<br>[l]</th><th>L<sub>x</sub><br>[l]</th><th>L<sub>y</sub><br>[l]</th><th>f<sub>bg</sub></th><th>Effective Res.<br>N<sub>x</sub> × N<sub>y</sub></th><th>AMR<br>levels</th></tr></thead><tbody><tr><td>Hs</td><td>0.1</td><td>25.6</td><td>12.8</td><td>0.2</td><td>6144 × 3072</td><td>6</td></tr><tr><td>Hb</td><td>0.05</td><td>51.2</td><td>12.8</td><td>0.2</td><td>24576 × 6144</td><td>5</td></tr></tbody></table>

We assume a uniform resistivity; $\eta = 5 \times 10^{-5}$, and an initialized magnetic and electric fields that are defined as

$$ B^x = \begin{cases} B_0 \tanh \left( \frac{y-L_y/2}{\delta} \right) + B_0 \epsilon_p & \text{for } y > 0 \\ -B_0 \tanh \left( \frac{y+L_y/2}{\delta} \right) + B_0 \epsilon_p & \text{for } y < 0, \end{cases} \quad (6) $$

$$ B^y = B_0 \epsilon_p, \quad (7) $$

$$ B^z = 0, \quad (8) $$

$$ E^x = E^y = E^z = 0. \quad (9) $$

Here, $\epsilon_p$ denotes a uniform (1 per cent) white noise perturbation to the magnetic field that varies between -0.01 and 0.01. This perturbation is similar to what is introduced (more naturally) in PIC simulations. Notice that we do not apply the typical magnetic field perturbation that guides the initial plasmons to the edges and creates a well-controlled reconnection region in the middle of the simulation domain (as prescribed for the GEM challenge, also in Keppens et al. 2013). To acquire pressure equilibrium at initialization, we define the fluid pressure to be (also outlined in Porth et al. 2014):

$$ p = \frac{B_0^2 \rho}{2}. \quad (10) $$

This implies that the initial plasma temperature $T = p/\rho = B_0^2/2 = 0.5$ and is therefore relativistically hot (cf. Ripperda et al. 2020). Additionally, we define the length and time-scale as a function of total system length $L = 2 L_x$, so that $(x, y) \in [-0.5L, 0.5L] \times [-0.125L, 0.125L]$ for Hb and $(x, y) \in [-0.5L, 0.5L] \times [-0.25L, 0.25L]$ for Hs with a typical time unit of $t_c = L/c$. Table 1 outlines the simulation length-scales and the Adaptive Mesh Refinement (AMR) resolutions, which are dynamically updated and applicable in all regions with sufficient density. It becomes apparent that Hb is effectively better resolved (i.e. $2L_x/N_x \approx 4 \times 10^{-3} l$ for Hb versus $8 \times 10^{-3} l$ for Hs).

For completeness, we note that the computational length unit is $l=1$ with corresponding time-scale $l/c=1$. If one were interested in relating the initial layer half-thickness $\delta$ (see Table 1) to the resistivity $\eta$, then one finds that $\delta/\eta = 1000$ ($\delta/\eta = 2000$) for Hb (Hs).

Nevertheless, we will connect it to a more intrinsic plasma-physical time-scale in our unit set in the following paragraph. This is typically determined by the upstream Alfvén velocity $v_a$, which is defined as

$$ v_a = \frac{B}{\sqrt{\rho h + B^2}} = \frac{\sqrt{\sigma}}{\sqrt{1+\sigma}}, \quad (11) $$

where $h = 1 + \hat{\gamma} p/(\hat{\gamma} - 1)\rho$ is the specific enthalpy with adiabatic index $\hat{\gamma} = 13/9$ and $B = \sqrt{B^2} = \sqrt{B^i B_i}$ denotes the magnetic field strength. Additionally, the ('hot') magnetization is defined as $\sigma =$

<footnote n="1">¹In our geometrized unit system, this is numerically equal to the magnetic diffusivity. However, conversion to other unit systems requires another factor in addition to the value of the time unit employed, namely $c^2$ for Lorentz-Heaviside units and $c^2/4\pi$ for Gaussian units.</footnote>
$B^2/\rho h$ and is set to an upstream value of $\sigma \approx 2$ at initialization. The Harris sheet magnetization $\sigma$ was set so that it is comparable with regions (i.e. disc-jet interface layer as seen in Fig. 6) in the GRMHD simulations. While we will primarily use the light-crossing time, it is worthwhile to connect it to the Alfvén and (resistive) diffusion time-scales of the system, which then become $\tau_a \approx L'/v_a$ and $\tau_d \approx L'^2/\eta$ with $L'$ being a representative but arbitrary current sheet length (Ripperda et al. 2019b). Fig. 2 gives an overview of the evolution of the Harris sheet (for the Hb case). From the magnetization ($\sigma$) panels, we find that $\sigma \sim 5$ near the sheet, which indicates an upstream Alfvén velocity $v_a \sim c$. Then, one can determine the Lundquist number via $S = \tau_d/\tau_a$, but it becomes clear that $\tau_d$ is very large and $\tau_a \sim L'\epsilon(0.5, 5)$ (in the relevant computational unit, i.e. $l$ or $r_g$), which indicates that $S = L'/\eta \in (10^4, 10^5)$ will be similarly large and we (generally) operate in the fast reconnection regime.

As mentioned, we set all boundaries conditions to fully periodic (similar to Cerutti et al. 2013, 2014; Keppens et al. 2013; Takamoto 2013, and some quasi-periodic works in Sironi & Spitkovsky 2014; Petropoulou & Sironi 2018). This implies that no matter is lost from the simulation domain so that the evolution eventually saturates after having formed several ‘monster’ plasmons that effectively act as matter reservoirs spanning a considerable part of the simulation domain. Up to a point, each sheet will evolve independently and uniquely due to the minor non-uniform perturbation to the initialized magnetic field, but when the primary plasmons become too large the sheets are influenced by one another. A different approach utilizes outflowing boundaries at the short sides of the box corresponding to the y-boundaries in our simulation (Loureiro et al. 2012; Sironi et al. 2016). Such boundary conditions tend to give less chaotic current sheets and allows for longer evolution times as, for periodic boundaries, the large plasmons will eventually affect the opposing current sheet. The periodic Harris sheet simulations are primarily meant as an additional, well-investigated verification case for the identification algorithm. Nevertheless, the periodic Harris sheet simulations tend to display more erratic behaviour than what is found for the outflowing variety, especially when combined with a global magnetic field perturbation (so that $\text{sign}(x) \cdot u_x \gtrsim 0$; Loureiro et al. 2012; Sironi et al. 2016).

## 2.3 GRMHD accretion disc configuration

We first outline a few specifics about the 3 + 1 split that is employed in the BHAC code. The line element is described as follows:

$$ds^2 = -\alpha^2 dt^2 + \gamma_{ij} (dx^i + \beta^i dt)(dx^j + \beta^j dt), \quad (12)$$

with $\alpha, \beta, \gamma$ denoting the lapse, shift, and geometric part of the metric ($g^{\mu\nu}$), where Roman characters $i, j \in \{1, 2, 3\}$ and Greek characters $\mu, \nu \in \{0, 1, 2, 3\}$. The metric determinant is then defined as $\sqrt{-g} = \alpha\sqrt{\gamma}$. Consistent with the conventions introduced in Porth et al. (2017), we denote electromagnetic quantities in the Eulerian frame with capitalized letters while lowercase letters denote quantities in the co-moving fluid (or plasma) frame. With Eulerian frame, we imply an Eulerian observer that is moving with four-velocity $n_\mu = \{-\alpha, 0, 0, 0\}$ (or contravariantly; $n^\mu = \{1/\alpha, \beta^i/\alpha\}$).

The initialized torus, that is in hydrodynamic equilibrium (Fishbone & Moncrief 1976, except for a small perturbation in the fluid pressure $p$), is threaded by a single poloidal magnetic field loop (that is initialized via $\mathbf{B} = \nabla \times \mathbf{A}$ with $\mathbf{A} = (0, 0, A_\phi)$), which is set by

$$A_\phi \propto \max \left( \frac{\rho}{\rho_{\max}} \left( \frac{r}{r_{\text{in}}} \right)^3 \sin^3 \theta \exp \left( -\frac{r}{400 r_g} \right) - 0.2, 0 \right). \quad (13)$$

<table><caption>**Table 2.** The model identifiers (IDs) and corresponding resolutions of the GRMHD accretion disc simulations. These simulations are all run with a dimensionless black hole spin $a_* = 0.9375$, adiabatic index $\hat{\gamma} = 13/9$, and simulation domain $r \in [1.185r_g, 1500r_g]$, $\theta \in [0, \pi]$. The density floor and magnetization ceiling are set to $\rho_{\min} = 10^{-4}$ and $\sigma_{\max} = 10^3$, respectively.</caption><thead><tr><th>ID<sub>s</sub></th><th>Type<br>GRMHD</th><th>η</th><th>Effective resolution<br>N<sub>r</sub>, N<sub>θ</sub></th><th>AMR<br>levels</th></tr></thead><tbody><tr><td>iM3</td><td>Ideal</td><td>–</td><td>2048 × 2048</td><td>3</td></tr><tr><td>iM4</td><td>Ideal</td><td>–</td><td>4096 × 4096</td><td>4</td></tr><tr><td>iM5</td><td>Ideal</td><td>–</td><td>8192 × 8192</td><td>5</td></tr><tr><td>rM3</td><td>Resistive</td><td>5 × 10<sup>-5</sup></td><td>2048 × 2048</td><td>3</td></tr><tr><td>rM4</td><td>Resistive</td><td>5 × 10<sup>-5</sup></td><td>4096 × 4096</td><td>4</td></tr><tr><td>rM5</td><td>Resistive</td><td>5 × 10<sup>-5</sup></td><td>8192 × 8192</td><td>5</td></tr></tbody></table>

The magnetic vector potential combined with the size of the disc at initialization eventually bring about the Magnetically Arrested Disc (MAD; Igumenshchev, Narayan & Abramowicz 2003; Narayan, Igumenshchev & Abramowicz 2003) state. The inner and pressure maximum radii of the torus that determine the size and available matter are set to $r_{\text{in}} = 20r_g$ and $r_{\max} = 41r_g$ for a black hole spin of $a_* = 0.9375$. The initialized magnetic field strength is scaled so that plasma-$\beta_{\max} = \max(p)/\max(p_{\text{mag}}) = \max(p)/\max(B^2/2) = 100$ is satisfied within the disc interior, with $p$ being the gas pressure and $p_{\text{mag}}$ the magnetic pressure. Other user-defined parameters of the evaluated configurations can be found in Table 2. For completeness, we note that the less magnetized accretion scenario is known as the Standard And Normal Evolution model (hereafter SANE; De Villiers, Hawley & Krolik 2003; Narayan et al. 2012; Sadowski et al. 2013), which is initialized with a different vector potential and tends to display a more (small-scale) turbulent accretion state.

To evolve the accretion disc surrounding the BH, we utilize the Modified Kerr-Schild (MKS) coordinate system (that is clearly described in McKinney & Gammie 2004; Porth et al. 2017). As the Kerr-Schild (KS) metric is well-documented (Misner, Thorne & Wheeler 1973), we will only comment on the modification from the standard KS coordinates $(t, r, \theta, \phi)$, which is done via;

$$r = R_0 + e^s, \quad (14)$$

$$\theta = \vartheta + \frac{h}{2} \sin(2\vartheta). \quad (15)$$

Here, $s$ and $\vartheta$ are the code's internally used coordinates, which can be converted to KS coordinates with the listed relations. We will exclusively show results in KS coordinates $r$ and $\theta$. All our GRMHD simulation use user-defined parameters $h = 0.25$ and $R_0 = 0$, which implies that the resolution of the underlying grid will be more concentrated in the equatorial plane.

BHAC inherently utilizes an AMR grid structure, but we only employ static refinement criteria that do not change over time for the GRMHD simulations. Without explicitly showing the grid structure, we note that all AMR levels are applied for $r < 200 r_g$ according to pre-defined criteria. We study three ideal and three resistive GRMHD models, which are run with different maximum AMR levels (see Table 2). Throughout this work, however, we will mostly comment on the highest resolution cases, which are: iM5 and rM5. Appendix A outlines a resolution convergence study between all AMR levels for both ideal and resistive GRMHD. There, we find increased activity in the equatorial plane for the resistive simulations, which is especially noteworthy as it is less pronounced in their ideal counterparts.
## 2.4 Energetics and surface averaged quantities

An important objective of this work is to quantify if plasmons are able to produce flaring events or create hot spots that would stand out with respect to the background. Therefore, we associate the electromagnetic, kinetic, and thermal fluid energies with their corresponding components of the stress–energy tensor $T^{\mu\nu}$ according to:

$$ \epsilon_{em} = -T_{EM}^{t} = -(b^2 + e^2)(u'u_t + \frac{1}{2}g't) + b'b_t + e'e_t \quad (16) \\ + \frac{u_{\lambda}e_{\beta}b_{\kappa}}{\sqrt{\gamma}}(u'\eta_t^{\lambda\beta\kappa} + u_t\eta^{\lambda\beta\kappa}), $$

$$ \epsilon_{\text{kin}} = -T_{\text{PAKE}}^{t} = -(u_t + 1)\rho u^t, \quad (17) $$

$$ \epsilon_{\text{th}} = -T_{\text{EN}}^{t} = -(\epsilon + p)u^t u_t - p. \quad (18) $$

Here, the hereto unexplained quantities are $\epsilon$, $p$, and $\eta^{\nu\lambda\beta\kappa}$, which are the specific internal energy, the fluid pressure, and the fully antisymmetric symbol, respectively. $\epsilon_{em}$ denotes the electromagnetic energy density (Qian et al. 2017), $\epsilon_{kin}$ the kinetic energy density, and $\epsilon_{en}$ the thermal energy density (McKinney, Tchekhovskoy & Blandford 2012; Ripperda et al. 2019a). The subscripts ‘EM’, ‘PAKE’, and ‘EN’ correspond to the electromagnetic, free particle, and enthalpy terms of the stress–energy tensor $T^{\mu\nu}$ (primarily following McKinney et al. 2012). The free thermokinetic energy (denoted as MAKE in McKinney et al. 2012) is the sum of $\epsilon_{kin}$ (PAKE) and $\epsilon_{th}$ (EN). This is important to note because $\epsilon_{kin}$ is predominantly negative in our GRMHD simulation, which can be interpreted from the geometric Bernoulli criterion ($u_t \le -1$) corresponding to unbound matter. The term ($u_t + 1$) will therefore be negative (positive) when the fluid element is unbound (bound) and as $\rho u^t$ is positive we will end up with a negative $\epsilon_{kin}$ for bound matter that is typically associated with the interior of the accretion disc. Lastly, note that the minus-sign in front of $T^{t}$ is due to the metric signature ($-+, +, +, +$) and is needed to get positive values.

Next, we define the covariant surface average (denoted by a bar, $\bar{Q}$, over a given fluid variable) by

$$ \bar{Q} = \frac{\int Q \sqrt{\gamma} dx^1 dx^2}{S} \quad (19) $$

with the surface $S$, in an arbitrary coordinate system, denoted as

$$ S = \int \sqrt{\gamma} dx^1 dx^2. \quad (20) $$

The $\gamma$ corresponds to the geometric part of the metric as explained in Section 2.3. Note that by surface average we imply that we take the average of a given quantity that is enclosed by a plasmoid-describing contour found by the algorithm. All averages are calculated in the Eulerian (or laboratory) frame.

## 3 PLASMOID IDENTIFICATION

The starting point of our plasmoid identification routine lies in finding a quantity that reveals the intrinsically circular magnetic field geometry. A natural choice for this quantity would then fall to the magnetic flux function, which is defined as

$$ \Psi_B \stackrel{\text{KS}}{=} \int \sqrt{-g} B^r d\theta, \quad (21) $$

$$ \stackrel{\text{Cart.}}{=} \int B^x dy - \int B^y dx, \quad (22) $$

where $\sqrt{-g}$ is the metric determinant. Part of the simulations utilize a special variety of Kerr–Schild (KS) coordinates $(r, \theta)$, which were introduced in Section 2.3. Note that $\sqrt{-g} B^r$ corresponds to the

magnetic field in the Eulerian frame and that the magnetic flux function $\Psi_B$ corresponds to the out of plane component of the magnetic vector potential (i.e. $A_\phi$, except for a minus sign discrepancy, Sironi et al. 2016). The magnetic flux function $\Psi_B$ is a good choice as its isocontours will follow the inplane magnetic field lines (i.e. $\mathbf{B} \cdot \nabla \Psi_B = 0$). More specifically, as plasmons are characterized by their circular magnetic field configuration, the plasmoid centre will correspond to a local maxima or minima in $\Psi_B$ (‘O-points’).

For our methodology, we work with a modified quantity with respect to the base $\Psi_B$ structure, which is defined as

$$ \tilde{\Psi}_B = \bar{\Psi}'_B - \Psi_B, \quad (23) $$

where the $\bar{\Psi}'_B$ scalar denotes the spatially averaged flux function. The removal of the averaged flux function yields a higher accuracy for the identification of plasmons. We note that this choice is purely numerically motivated – a particular step (i.e. iii) in the algorithm is only applicable for local minima and achieves better accuracy with the average subtracted.

After having selected a suitable plasmoid identification quantity (i.e. $\tilde{\Psi}'_B$), we need a method that can reliably classify the magnetic island structure and size. For this purpose, we developed an algorithm that consists of the following four steps:

(i) All simulations contain a significant fine-structure in the magnetic flux function. This can make it hard to differentiate between (magnetic) turbulence and more global features that correspond to a presence of a plasmoid. Therefore, to guarantee we filter out much of the turbulence, we apply a blurring (Gaussian or flat) kernel to the flux function ($\tilde{\Psi}'_B$). This also gives us control over the size of the features we want our algorithm to be sensitive to. The blurring step, however, requires (manual) fine-tuning depending on resolution and nature of the setup. Interestingly, to extract the global structure of the highly turbulent primary plasmoid in the isolated Harris sheet simulation(s), one needs the strongest blurring procedure, while the GRMHD simulation are well-served with a fairly light blurring procedure.

(ii) Following the blurring step, we identify the local minima or maxima that will correspond to the plasmoid’s centre.

(iii) Then, we apply a watershed algorithm (well-described in, e.g. Beucher & Meyer 2018) to isolate the domain of interest around the local minimum. We have chosen an implementation that is based on Vincent & Soille (1991). The watershed segmentation is then used to make an informed cut-out of the domain that will contain a single (local) extrema, so that we have control over what is being fitted while simultaneously improving the quality of the fit. We note that as the watershed is only able to identify local minima one has to run the pipeline twice to identify all local extrema–once for $\tilde{\Psi}'_B$ and once for $-\tilde{\Psi}'_B$.

(iv) Lastly, we draw the maximally possible contour within the isolated segment. Utilizing the inherent symmetry in the systems, we sample the space efficiently by means of a binary search from opposite sides (i.e. left and right from centre along $\hat{x}$ for the Harris sheet and inner and outer radii along $\hat{r}$ for the GRMHD setups). The resulting contour enables us to gauge the plasmoid’s size and orientation, and enables calculations of the plasma quantities associated with the plasmoid and its direct vicinity.

In Fig. 1, one finds a schematic summary of the points discussed above. Additionally, it becomes clear that both setups differ fundamentally from one another and therefore warrant a different configuration of the algorithm. The main differences are summarized below. (i) As the Harris sheet setups have periodic boundaries, one needs to be careful to catch plasmons that are on the boundary. (ii) Additionally, capturing both ‘big’ and ‘small’ plasmons in the Harris
<figure><img src="image_5.png" /><figcaption>**Figure 1.** A schematic decomposition of the plasmoid identification algorithm. In the top panels (a–d), we display a snapshot of a GRMHD simulation (±M3 at T = 3000 $r_g/c$) at various points in the pipeline. In the bottom panels (e–h), we find the same but for one of the Harris sheet cases (Hb at T = 2.93 $t_c$). In the left column (panels a and e), one finds the base magnetic flux function $\Psi_B$ – the starting point. To apply the watershed (panels c and g), one needs to make sure that the plasmoid corresponds to a local minimum which is done with the quantity $-\tilde{\Psi}_B$ (panels b and f). The last column (panels d and h) showcases how the maximal contour is found for the watershed segment and how the plasmoid’s width and height are determined (between the orange diamonds). The evaluated O-point is denoted by the black circle. Other O-points in the displayed simulation domain are denoted by the open grey circles.</figcaption></figure>

sheet setup requires two different approaches, mainly concerning the blurring kernel. For the big features, one has to apply a relatively small kernel many times (several hundred times works well in our experience) to retain as much global structure as possible. To capture the smallest plasmons, one only has to apply the small blurring kernel a few times. One acquires the master set by combining the output from both described configurations where the largest plasmoid structure is leading (following the later defined Uzdensky et al. 2010 criterion). (iii) For GRMHD, one has to take into account that the resolution is concentrated near the black hole and in the equatorial plane and therefore has non-uniform cell-sizes. (iv) Due to the non-uniform GRMHD grid layout, applying a kernel blur manifests itself differently in various regions of the simulation domain. When applying a relatively small blurring kernel, this effect is minor and manageable. If this is not sufficient, we interpolate $\tilde{\Psi}_B$ to a uniform grid structure.

Lastly, we would like to note more explicitly how plasmons are identified using the magnetic flux function in other works. In essence, one identifies plasmons via the so-called ‘O’- and ‘X’-points. O-points corresponds to the local minima and maxima of the magnetic flux function and denote the center of a plasmoid. X-points are saddle points and lie in between O-points. Along a current sheet one therefore expects these points to succeed one another. One typically finds the extrema by calculating the Hessian matrix of the magnetic flux function (Servidio et al. 2009, 2010; Zhdankin et al. 2013; Kadowaki, Pino & Stone 2018; Zhou, Loureiro & Uzdensky 2020) via;

$$H_{ij}^{\Psi_B}(\mathbf{x}) = \frac{\partial^2 \Psi_B(\mathbf{x})}{\partial x_i \partial x_j}. \qquad (24)$$

Then, one calculates the matrix determinant of the Hessian ($|H^{\Psi_B}|$) to find the critical points that correspond to $|H^{\Psi_B}(\mathbf{x})| = 0$ at a

given coordinate $\mathbf{x}$. For an O-point, the eigenvalues of the Hessian then determine if it is a local minima (positive definite Hessian) or maxima (negative definite Hessian). For an X-point, one finds both positive and negative eigenvalues of the Hessian (Servidio et al. 2010). However, in our methodology, there is no need to explicitly calculate the Hessian to identify the O- and X-points as these are naturally isolated by the watershed algorithm. The X-points, which are typically harder to identify (Zhdankin et al. 2013), will lie on the border of a watershed segment (which is a property that was utilized explicitly in Winarto & Kunz 2022). For the O-points, we straight-forwardly calculate the local extrema in a segment. Finding the critical points in these turbulent maps is a complicated endeavour, as is also illustrated by the computationally intensive mitigation techniques (as introduced in Servidio et al. 2010). Our methodology works around this problem in a relatively natural manner, but this implies that we do not know the exact orientation of the current sheet as the X-point locations are not calculated (other than being on the watershed segment’s border). Additionally, one can end up with two O-points per watershed segment, but this is straightforwardly mitigated by the contour-finding algorithm as it only selects the contour enclosing the O-point in question. Even though we may sacrifice some accuracy, our methodology saves us from having to employ (relatively) computationally and memory intensive mitigation strategies and will therefore provide a significant speed-up with respect to (e.g. Servidio et al. 2010).

## 4 RESULTS

In this section, we outline the findings of the plasmoid detection routine for both classes of simulation separately. Section 4.1 summarizes
the results associated with Harris sheet simulations, while Section 4.2 summarizes the results associated with the GRMHD simulations.

## 4.1 Harris sheet

### 4.1.1 General evolution

In Fig. 2, a well-developed and representative state of the Hb case is shown. Before this state is reached, the current sheet needs to evolve for some time before it becomes (plasmoid- or tearing-) unstable enough, as the sheet becomes thinner, to break up and form the first magnetic islands. This first tearing mode creates the first plasmons that are known as primary plasmons (see e.g. Loureiro et al. 2007; Comisso et al. 2016; Uzdensky & Loureiro 2016; Petropoulou & Sironi 2018) and have significantly different plasma characteristics than the ones that are created at later times in the secondary tearing-unstable regions of the sheets. First, the primary plasmons have higher densities and, second, they possess a characteristic magnetic field profile with a lower magnetic field strength at the center than in the rings further on the outside. This results in a lower overall magnetization, but also a relatively lower surface-averaged magnetic field strength. Their plasma properties are primarily determined by the initial conditions. Following the initial break-up of the layer (at ~1.56 $t_c$ for Hb and ~5.27 $t_c$ for Hs), a continuous and steady creation of secondary plasmons commences in the reconnection layers between the primary islands. The layer remains active till the very end of the evaluated time window. The secondary plasmons do probe the underlying plasma characteristics and are relatively unaffected by the initial conditions. Two animations are attached to Fig. 2 which show both the spatial layout corresponding to the figure as well as the entire simulation domain over time.

Following the criterion outlined by Uzdensky et al. (2010) that when a plasmoid coalesces with a larger plasmoid, then the smaller one is considered to be incorporated into the larger body. From coalescence onwards, the smaller plasmoid is therefore no longer considered as an separate entity. In practice, however, the small plasmoid will retain its structure for some time (ranging several 0.05 $t_c$ depending on its size) before dissipating into the global structure of the primary plasmoid. This is clearly illustrated in Fig. 2 and accompanying animations, the coalescence of the plasmoid on the left-hand side (at $X = 0.135 L$ and is roughly 0.02 $L$ in width initially) takes approximately 0.1 $t_c$ from the moment of impact to being fully absorbed by the primary plasmoid. When two plasmons of similar size coalesce, then this time-scale tends to be longer and significant perturbation is needed before one of the two loses its structure.

Generally, it is not simple to enforce the Uzdensky et al. (2010) criterion, which is reflected by the two-step approach outlined in Section 3. Starting with secondary plasmons, the minimum size for which we identify this population is set to ~$10^{-4} L$ (0.005 $l$), but in practice the algorithm tends to detect a plasmoid when it starts to deviate from the straight current sheet configuration (i.e. gain some width). Overall, we find that the secondary plasmons are identified with a very high fidelity. The primary plasmons are typically much harder to identify as they are the end point of the inverse cascade (or plasmoid coalescence) and, therefore, act as highly turbulent plasma reservoirs that will never relax as smaller plasmons keep colliding and merging into them. These continuous perturbations also give rise to magnetic reconnection and plasmoid formation within the primary plasmoid structure. As described in Section 3, we need to apply an aggressive blurring kernel to identify the global primary plasmoid structure, but we still want to remain sensitive to the distinct plasmoid

structure if they have not fully merged yet. This implies that two plasmons that have a similar magnetic flux signature (an example is seen at $X = 0.34L$ in Fig. 2g) are still picked up as two separate entities even though one can argue that they are actually part of one global body, especially when following the Uzdensky et al. (2010) criterion. At the interface of these two plasmons one often finds new plasmons forming. Naturally, all previously mentioned points become less pronounced at lower resolutions as one is resolving the current sheets less well which results in a lower number of formed plasmons and less fine-structure.

The end of the evaluated time window (at 4.1 $t_c$ for Hb and 8.79 $t_c$ for Hs) is determined by the amount of interference the current sheets have on one another. Beyond these times (for which no representative state is shown in this work), the few primary plasmons become of sufficient size that they start to incorporate the opposing current sheet. This brings about an interesting new turbulent mode that is similar to the ABC structure described in Lyutikov et al. (2016). Magnetic reconnection is then no longer confined to the current sheets but occurs at interfaces between the (opposite polarity) primary plasmons that now have lost their elliptical shape and are more hexagonal in shape. As the turbulent ABC mode is beyond the scope of this work, we selected the evaluated time windows based on the presence of a clear current sheet structure.

### 4.1.2 Plasmoid statistics

Fig. 3 displays 2D histograms with various plasmoid quantities as a function of width for both Harris sheet (Hs and Hb) cases. First, we would like to point out that the distributions show the same general trends. Starting with the surface-averaged density ($\bar{\rho}$) panels, one finds a main triangular distribution that spans $-1.25 < \log_{10} \bar{\rho} < -0.25$. In addition to the main distribution, there is a secondary channel corresponding to $-0.25 < \log_{10} \bar{\rho} < 0.25$ that corresponds to the densest plasmons which also seem to occur over the entire width range. This dense plasmoid population is partly explained by misclassification of the global structure and partly by the simulation conditions quickly after the initial break-up of the layer. For the former, we find that the primary plasmoid is occasionally not identified properly, which is often tied to strong fluctuation in the magnetic flux function associated with the primary plasmoid. These fluctuations, within the primary plasmoid's interior, will then be identified and are responsible for the 'dense' population with a small plasmoid half-width. For the latter scenario, there are a number of high density matter reservoirs that will eventually become the primary plasmoid population and will generally correspond to a large plasmoid half-width. Overall, misclassifications occur rarely, as is demonstrated by the otherwise compact distributions.

Returning to the 'true' plasmoid population, spanned by $-1.25 < \log_{10} \bar{\rho} < -0.25$, we find that the smallest detected plasmons have a half-width $w \approx 2 \times 10^{-4} L$ for both the Hb and Hs cases. This lower limit is partially set by an identification requirement that either the width or height of the contour spans at least 5 cells (which equates to a minimal width or height of $\Delta x \approx 0.02 l$) for the evaluated data. For the surface-averaged magnetization ($\bar{\sigma}$), we find that the main population spans $-2.5 < \log_{10} \bar{\sigma} < 0.5$. As is also seen in Fig. 2, the secondary plasmons have a remarkably similar $\sigma$ profile with the outer shells being more magnetized than the interior (similar to findings in Petropoulou & Sironi 2018). Nevertheless, we do find a trend where the $\bar{\sigma}$ rises with half-width, up to $w \approx 6.3 \times 10^{-3} L$. For $\log_{10} w/L > -2$, the $\bar{\sigma}$ mean plateaus and even seems to
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Representative state for the evolution of the Harris sheet for the Hb case corresponding to $T = 2.46 t_c$. Rows (a) till (g) show the density $\rho$, ‘hot’ magnetization $\sigma = B^2/\rho h$, plasma $\beta = p/(B^2/2)$, electromagnetic energy density $\epsilon_{em}$, kinetic energy density $\epsilon_{kin}$, thermal energy density $\epsilon_{th}$, and magnetic flux function $\Psi_B$. The magenta contours denote plasmoid detections corresponding to local maxima in the flux function ($\Psi_B$), while green contours correspond to local minima. The evolution over time is displayed in two animations; one for the zoom-in corresponding to this figure and another displaying the entire simulation domain, which can be found in the following repository; https://doi.org/10.5281/zenodo.8318522.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** 2D distributions $N(w, x)$ of the plasma quantities $x \in \{\bar{\rho}, \bar{\sigma}, \bar{\beta}, \bar{\epsilon}_{em}, \bar{\epsilon}_{kin}, \bar{\epsilon}_{th}\}$ as a function of plasmoid half-width $w$ (with $L = 2L_x$, as per Table 1) for both the Hb (left-hand panels) and Hs (right-hand panels) cases. We stack the distributions as a function of time and divide by the maximum. The green line denotes the mean per width bin that has more than ten counts total.</figcaption></figure>

decrease slightly for the largest plasmons. After the growth phase (in $\log_{10}w/L \lesssim -2$), it seems that the increase in density and magnetic field strength is roughly matched. Lastly, for $\bar{\beta}$, we find a similar but inverse trends to what we described for $\bar{\sigma}$. The part of the distribution with the largest plasmons ($w \sim 0.1 L$) seems to deviate significantly from the main population and possesses a relatively high $\bar{\beta} \gtrsim 10^3$. This happens because at the center of the plasmoid the magnetic field strength becomes very small due to the circular configuration. This generates some very high $\beta$ values that in turn affects the surface-averaged quantities ($\bar{\beta}$).

For the energies ($\bar{\epsilon}_{em}$, $\bar{\epsilon}_{kin}$, and $\bar{\epsilon}_{th}$), we find that the thermal energy ($\epsilon_{th}$) is the leading term in the total energy budget of the plasmons with a mean (denoted by the green line) that remains fairly constant ($0.0 < \log_{10}\bar{\epsilon}_{th} < 0.25$) as a function of half-width ($w$). At smallest $w$, it appears the second term is the electromagnetic energy (at $\bar{\epsilon}_{em} \approx 10^{-1.5}$) that steadily becomes more significant for increasing width. As the kinetic energy ($\epsilon_{kin}$) is closely tied to the velocity of the plasmoid, we find that it can actually become a competing term for the electromagnetic energy, especially in the active reconnection regions and merging (or colliding) plasmons (see Fig. 2). The distribution of $\bar{\epsilon}_{th}$ and $\bar{\epsilon}_{kin}$ is wide indicating significant variance, while $\bar{\epsilon}_{em}$ closely follows the distribution of $\bar{\sigma}$ and seems to show a more consistent trend. This trend is explained by secondary plasmons becoming

more magnetized with time until they grow up to a size of $w \sim 0.01L$. Then, the secondary plasmoid generally encounters a primary plasmoid and is subsequently absorbed, after which the growth in magnetization ($\bar{\epsilon}_{em}$) stagnates. The high variance in $\bar{\epsilon}_{kin}$ is explained by the localized nature of plasmoid acceleration – predominantly in active reconnection regions and just before plasmoid coalesce. As soon as the secondary are absorbed by the primary plasmons, $\bar{\epsilon}_{th}$ will be the leading term by a significant factor. Even though $\bar{\epsilon}_{th}$ is still most dominant in the secondary plasmoid, both $\bar{\epsilon}_{kin}$ and, especially, $\bar{\epsilon}_{em}$ can become close in significance.

Lastly, we would like to briefly comment on the differences between the two cases; Hb and Hs. So far, we have mainly described the Hb case shown in the left-most panels of Fig. 3. Nevertheless, we find that all findings based on Hb are also applicable to Hs. The description of both simulations is outlined in Table 1, where we find that the main differences lie in the initial layer (half-)thickness ($\delta$) that is twice as wide and that the resolution is twice lower. This also explains why the evolution starts later for Hs; it takes longer for the perturbations to create a sufficiently thin current sheet to activate the tearing instability. Additionally, the simulation box length (in $\hat{x}$, long side) is halved and contains more matter due to the thicker initial layers when compared to the Hb case. As there are only relatively minor differences, we find that their evolution is similar, which is
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Probability density function $f(w/L)$ of plasmoid half-width $w$ (along $\hat{y}$) that is scaled according to the total simulation box width ($L = 2L_x$) for both the Harris sheet cases; Hb and Hs. From scaling arguments, proposed in Uzdensky et al. (2010) and Loureiro et al. (2012), it has been shown that the distribution is expected to scale $\propto w^{-2}$. The various probability density function profiles are coloured according the time at which they occur over a range of $T \in [1.76, 4.10]t_c$ for Hb and $T \in [5.47, 8.79]t_c$ for Hs. The mean density profile and $1\sigma$ error over time are denoted by the dashed dark grey line. The power-law index is determined via $p = -d \log f/d \log(w/L)$.</figcaption></figure>

also reflected by the results here. Most notably, it seems that the primary plasmons spans a greater section of the simulation domain for Hs. To gain insight into the dependence of the plasmoid dynamics on starting conditions, a more detailed study is needed, but that lies beyond the scope of this work.

### 4.1.3 Plasmoid distribution functions

Figs 4 and 5 display a probability density function ($f$) of plasmoid half-width and the absolute plasmoid surface averaged magnetic flux function ($|\tilde{\Psi}_B|$), respectively. The distributions are calculated at a time cadence of $10^{-2}t_c$ starting at the beginning of the evaluated window at $T = 1.76t_c$ for Hb ($T = 5.47t_c$ for Hs) up to $T = 4.1t_c$ ($T = 8.79t_c$) shown by the colour range of dark blue up to bright yellow. Starting with Fig. 4, we find variation of the probability density function as a function of time, but a consistent image emerges as well. Generally speaking, at the smallest plasmoid half-widths (up to $w/L \approx 10^{-3}$), we find a plateau followed by a steady decrease in occurrence frequency as the plasmons become larger, up to the largest plasmons that span a tenth of the simulation domain ($w \sim 0.1L$). We primarily quantity the growth rate of plasmons in the system by evaluating the power-law index $p = -d \log f/d \log(w/L)$ of the density function.

The plasmoid scaling laws have been studied in detail in the past (Uzdensky et al. 2010; Huang & Bhattacharjee 2012; Loureiro et al. 2012; Sironi et al. 2016). The density function of plasmoid width was predicted and verified to scale according to $f(w) \sim w^{-2}$ (Uzdensky et al. 2010; Loureiro et al. 2012), while for magnetic flux both $f(\tilde{\Psi}_B) \sim \tilde{\Psi}_B^{-2}$ (following the same works) or $f(\tilde{\Psi}_B) \sim \tilde{\Psi}_B^{-1}$ (Huang & Bhattacharjee 2012) were established. The main difference between scaling found by Uzdensky et al. (2010) and Huang & Bhattacharjee (2012) lies in how they treat the relative velocity between plasmons. While Uzdensky et al. (2010) assumed it to be $\sim v_a$, Huang & Bhattacharjee (2012) evaluate a size-dependent relative velocity (see also Sironi et al. 2016). As our simulations have no large-scale magnetic field perturbation (or outflowing boundaries), relative velocities between plasmons are stochastically determined and relatively low, so we expect a greater similarity with Huang & Bhattacharjee (2012). Overall, we find that $\tilde{\Psi}_B$ and $w$ do not scale with the same power-law index $p$, which is contradictory to earlier works (Loureiro et al. 2012; Sironi et al. 2016). However, there are clear explanations for this perceived discrepancy that will be outlined in the next paragraphs. We note that we display the $\tilde{\Psi}_B$ quantity in Fig. 5, which is the surface-average of $\Psi_B$ as enclosed by the identified plasmoid contour. This quantity is comparable to the scaling relations for $\Psi$ listed in Loureiro et al. (2012); Sironi et al. (2016).

For the half-width ($w$) scaling relations in Fig. 4, we find a power-law index $p = 1.81 \pm 0.05$ for the Hb case and $p = 1.48 \pm 0.06$ for the Hs case. Overall, we find a scaling relation that is close to $f(w) \sim w^{-2}$ corresponding to $p = 2$ (indicating a surprising similarity with non-relativistic results by Uzdensky et al. 2010; Loureiro et al. 2012). For the mean trend in magnetic flux (in dark grey), we find $p = 0.64 \pm 0.10$ for the Hb case and $p = 0.59 \pm 0.06$ for the Hs case. However, the trend described by the smallest values per bin (in light grey) is $p \approx 1$, which indicates agreement with Huang & Bhattacharjee (2012). The evolution of the distributions is characterized by a relative over-representation of large plasmons, with $|\tilde{\Psi}_B/B_0L| \in [5 \times 10^{-2}, 10^{-1}]$, that expands itself both to the left (lower $|\tilde{\Psi}_B|$, smaller plasmons) and right (higher $|\tilde{\Psi}_B|$, larger plasmons) over time. The smallest plasmons have the lowest magnetic fluxes (as is also verified in Fig. 2) and the largest plasmoid will increase in $|\tilde{\Psi}_B|$ over time. This evolution also creates
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Probability density function $f(|\tilde{\Psi}_B|/B_0L)$ of plasmoid surface-averaged flux $|\tilde{\Psi}_B|/B_0L$ that is scaled according to the total simulation box width ($L = 2L_\alpha$) and initial magnetic field strength ($B_0$) for both the Harris sheet cases; Hb and Hs. The light grey dashed lines denote the mean of the lowest (1 per cent) values per bin with a corresponding linear fit shown in the same colour. For the rest, the description of Fig. 4 also applies here.</figcaption></figure>

the sizable $1\sigma$ error (visually made worse by the log-scale) as the density function evolves significantly over time. So, in short, the magnetic flux distributions evolve with $p = 1$ over time (especially for a low $|\tilde{\Psi}_B|$), but this relation is affected by a high $|\tilde{\Psi}_B|$ population (that is present from the start). This population is there because of the periodic boundary conditions and would not be over-represented when utilizing outflowing boundary conditions, which was done by the comparative studies.

Next to differences in the physical regime (i.e. non-relativistic, incompressible, and higher Lundquist numbers for Huang & Bhat-tacharjee 2012; Loureiro et al. 2012), we note that our simulation configuration differs substantially on at least two fronts from the previously mentioned scaling law studies, namely that it is rela-tively unperturbed and that it has no outflowing boundaries. With unperturbed, we mean that there is no large-scale magnetic field perturbation present. This type of perturbation would recreate a clean reconnection layer in the middle of the box and guides the primary plasmoids to the edge of the simulation domain (also discussed in detail in Section 2.2). In practice, this implies that (i) coalescence of plasmoids is a relatively prominent growth channel in our simulations and (ii) large plasmoids could disproportionately affect the distribution. The latter point is two-fold; as the primary plasmoids become larger they effectively shrink the domain where the (secondary) current sheets can form and they will eventually start interfering with the opposing current sheet. Especially for the Hs case, these points are influential, which is also accentuated by the larger deviations. All these effects are likely to play a role in explaining the differences in scaling found in this work with respect to previous works. Additionally, the informed (but arbitrary) choice regarding which bins to include for the fit combined with the imperfect sampling of the distribution by the bins also introduces a $O(5$ per cent) error on the values of $p$. Although, even despite the differences in simulation configuration (and the

numerical uncertainties), we still reach a remarkable consistency with previous studies that employed more idealized configurations for finding plasmoid scaling.

## 4.2 GRMHD

### 4.2.1 General evolution

Figs 6 and 7 display the typical structure of the axisymmetric MAD (cf. Tchekhovskoy, Narayan & McKinney 2011; McKinney et al. 2012) simulations. After having evolved sufficiently, they will saturate in magnetic flux that penetrates the event horizon (see Section 4.2.4). Following such a saturation event, the accretion flow is completely halted in axisymmetric simulation, while in 3D a so-called ‘flux tube’ forms (Dexter et al. 2020; Porth et al. 2021). Instead of halting the accretion flow completely, a localized less dense, more magnetized cavity (i.e. the flux tube) moves outward from the black hole and eventually dissipates into the accretion flow (after up to several orbits) in 3D. These outbursts occur semi-periodically and seem to be even more prevalent in the relatively more confining 2D simulations.

Another feature of MAD models is that the Magneto-Rotational Instability (MRI; Balbus & Hawley 1991), responsible for angular momentum transport, is suppressed or only marginally influential as the main magnetic field component is strongly poloidal (Porth et al. 2021; Begelman, Scepi & Dexter 2022, and reference therein). The MRI does play a role in the early developing phase of the simulation, when it is less magnetized, but then one of the leading causes of turbulence (close to the BH) is the Rayleigh-Taylor Instability (RTI; Avara, McKinney & Reynolds 2016; Marshall, Avara & McKinney 2018, and references therein). The RTI drives the break-up of coherent (magnetic) structures such as flux tubes towards the central BH. The turbulence associated with the Kelvin-
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Overview of the iM5 simulation at $T = 3840 r_g/c$. Here, we are in the middle of a flux eruption event with pushed back the accretion disc. The magenta contours corresponds to local maxima and the green contours correspond to local maxima in the magnetic flux function ($\Psi_B$). In panels (a)-(d), one finds the density ($\rho$), the magnetization ($\sigma$), the ideal to magnetic pressure ratio ($\beta$), and the magnetic flux function ($\Psi_B$). In the panels (e)-(h), we find the electromagnetic energy ($\epsilon_{em}$), kinetic energy ($\epsilon_{kin}$), thermal energy ($\epsilon_{th}$), and the magnetic Bernoulli factor ($B_m = -(h + \sigma/2)u_t$). The corresponding animations can be found in the following repository; https://doi.org/10.5281/zenodo.8318522.</figcaption></figure>

Helmholtz Instability (KHI) becomes important in regions with strong shear flows and is characterized by swirl-like vortices (see e.g. Begelman, Blandford & Rees 1984; Hillier 2019). In the black hole accretion environment, such shear layers are naturally associated with the jet-disc interface. Both of these instabilities are perturbative channels that are able to set off magnetic reconnection in the accretion disc. Therefore, for MAD accretion discs, we find a much more turbulent environment than for the Harris current sheet for which reconnection is only determined by the tearing instability (Ripperda et al. 2017) that is triggered in a relatively controlled scenario.

As we are mainly interested in the plasmons' ability to produce flares, which are known to originate close to the central black hole, we apply our algorithm only within the inner 25 $r_g$. In Figs 6 and 7, we display maps of the plasma quantities and energies (similar to Fig. 2) for the iM5 and xM5 cases. The magenta and green colours denote plasmons found with a local maximum and local minimum in the magnetic flux function, respectively. Both figures show typical phases of MAD evolution that occur in all the GRMHD simulations in this work. Both the ideal (i.e. iM5) and resistive (i.e. xM5) GRMHD simulations have a similar global evolution, except that xM5 produces significantly more plasmoid and does therefore have a more variable
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Overview of the rM5 simulation at $T = 3500 r_g/c$. Here, we find an accretion state that is standard for MAD simulations with a turbulent but fairly steady flow. The rest of the description is analogous to Fig. 6. The corresponding animations can be found in the following repository; https://doi.org/10.5281/zenodo.8318522.</figcaption></figure>

disc structure as is discussed in detail in Section 4.2.4. The panels (a–h) of Fig. 6 correspond to a flux eruption where we find the accretion flow is entirely halted. The panels (a–h) of Fig. 7 show a fairly generic accretion state with the turbulent accretion flow extending up to the horizon. Even though the density is low near the BH, one does find a reconnection layer along the equatorial plane (denoted by the magenta contours). These plasmons are the collisional (non-pair-production plasma) equivalent to what has been seen in GRPIC simulation of diffuse collisionless magnetospheres around BHs (Bransgrove, Ripperda & Philippov 2021; Crinquand et al. 2021).

The overall structure and location of the plasmoid chains indicate that at the disc-jet boundary, one finds plasmons that correspond to

local maxima (magenta), while when plasmons occur within the disc they correspond to local minima (green). The magenta contours seem to have a lower density ($\rho$) and higher magnetization ($\sigma$) than those that lie within the disc (denoted in green). The magenta contours also seem smaller when compared to their green counterparts. Their location and smaller size indicates that they are likely created by the shear-induced KHI. The magenta contours also tend to leave the identification domain ($r \le 25 r_g$) on short time-scales (5–10 $r_g/c$) as they rapidly move outwards with turbulent jet-disc layer (often referred to as the jet sheath). The green contours are tied to the bulk motion of the disc giving these plasmons more time and matter to interact with, which explains their larger sizes. Lastly, we would like to point out that the quantities (visible in the $\rho$, $|\epsilon_{\text{kin}}|$, and $\epsilon_{\text{th}}$ maps)
near the vertical axis ($x = 0$ $r_g$) are due to floor violations, which happen far from our areas of interest and will therefore not interfere with the analysis.

## 4.2.2 Plasmoid statistics

Fig. 8 shows 2D histograms with various plasmoid quantities as a function of width for the two GRMHD simulations (iM5 and rM5). We find a significantly lower plasmoid count for the iM5 case compared to the rM5 case. Hence, the ideal distributions are more sparsely sampled. We will address this point in more detail in Section 4.2.4. Overall, however, we do find that the distributions of iM5 and rM5 are consistent with one another. Before we start describing the distributions, we would like to note that we can no longer use the Euclidean width for the GRMHD cases, as it does not inherently take into account the space-time curvature. Instead, we have chosen to display the distribution of plasmoid sizes as a function of ‘circular’ radius $R_S = \sqrt{S/\pi}$ as the surface S calculation is taking into account the curvature. As plasmons are generally elliptical, we loose information about the shape because the distinction between width and length is no longer made.

Starting with the distribution of $\bar{\rho}$ (panels i.a and r.a of Fig. 8), we find that the surface-averaged density is highest for the smallest plasmons at $log_{10} \bar{\rho} \approx 0.75$ and then plateaus at $log_{10} \bar{\rho} \approx -0.5$ from $-0.5 < log_{10} R_S < 1.0$. For $\bar{\sigma}$ (i.b and r.b), we find a roughly constant mean value of $log_{10} \bar{\sigma} \approx -1$, but a wide spread in values is also present. For $\bar{\beta}$ (i.c and r.c), one finds a very elongated distributions centred around a mean of roughly $log_{10} \bar{\beta} \approx 0.5$ which has a complicated origin. This behaviour is largely explained by the ‘green’ (local minima in $\Psi_B$) and ‘magenta’ (local maxima in $\Psi_B$) plasmoid populations (hereafter simply referred to as ‘green’ population or distributions, analogously for ‘magenta’). For the magenta population, we find the origin of the elongated $\bar{\beta}$ distribution as the plasmons detected in the jet sheath correspond to a distribution centred on a relatively low $log_{10} \bar{\beta} \approx -1$. The typical plasmoid values for $\bar{\beta}$ are fairly uniformly distributed with $-2 \lesssim log_{10} \bar{\beta} \lesssim 2$ centred around a mean of $log_{10} \bar{\beta} \approx 0 - 0.5$. Similar behaviour is also seen for the distributions of $\bar{\rho}$ and $\bar{\sigma}$ where they both display near-identical means but a larger variance is present for the magenta contours. For the plasmons corresponding to green contours, we find more uniform and compact distributions overall that are located around the mean trends for the entire (combined green and magenta) distribution as shown in Fig. 8. We note that we do not explicitly show the results for both green and magenta populations separately, but rather only comment on it in the text for brevity’s sake.

For the energies ($\epsilon_{em}$, $|\epsilon_{kin}|$, and $\epsilon_{th}$ shown in Fig. 8), there are only minor differences between the green and magenta distributions, so we will just discuss the combined distributions for the energies in panels (i.f-i.h and (r.f-r.h). Interestingly, the means for all energy distributions follow an almost identical trend – starting at $log_{10} \bar{\epsilon} \approx 0$ to ending at $log_{10} \bar{\epsilon} \approx -2$ for increasing $R_S$. After a rapid decline up to $log_{10} R_S \approx -0.5$, we find that the means of $log_{10} \bar{\epsilon}$ plateau, especially for $\bar{\epsilon}_{kin}$ and $\bar{\epsilon}_{th}$. Additionally, the distributions indicate that the various surface-averaged energies are of similar strength. $\bar{\epsilon}_{em}$ does stand out, however, with respect to the other energies as it displays a more compact distribution with a clear, gradually declining trend. Generally speaking, we find that all energy densities are of similar strength independent of the plasmoid size. Continuing with $\bar{\epsilon}_{kin}$, $\bar{\epsilon}_{kin}$ is mostly negative, except in the jet-sheath where $\epsilon_{kin} \sim O(1)$. This is

explained in detail in Section 2.4. Here, we only consider the absolute value $|\bar{\epsilon}_{kin}|$ (i.g and r.g). The dashed magenta and green lines in these panels correspond to the means of the distributions containing only the positive or negative values of $\bar{\epsilon}_{kin}$, respectively. So, it becomes clear that the vast majority of plasmons has a negative $\bar{\epsilon}_{kin}$ value as the global mean (in solid green) lies close to the dashed green line. Lastly, for $\bar{\epsilon}_{th}$, we find similar behaviour as for the other energies, except it has the strongest contribution at the lowest plasmoid sizes.

To assess how freely a plasmoid moves through the accretion disc, one can evaluate the boundedness of the plasmoid to the global accretion flow, which is described by the Bernoulli factor or parameter (shown in panels i.i and r.i of Fig. 8). We define the magnetic Bernoulli factor as $Be_m = -(h + \sigma/2)u_t$, which incorporates the contribution of the magnetic pressure ($\sigma/2$) and therefore deviates slightly from the standard relativistic Bernoulli $Be = -hu_t$ (Rezzolla & Zanotti 2013). The Bernoulli criterion states that the fluid is unbound when $Be_m > 1$. Note that we have taken the liberty to incorporate a minus sign within the Bernoulli factor. Returning to the distributions in panels (i.i and r.i), we find the majority of surface-averaged plasmons is unbound as they pass the criterion, but there is still a significant number that lies under and close to the critical value of $\bar{Be}_m = 1$ and are therefore bound. The mean of the function does, however, indicate $\bar{Be}_m \approx 1$ with a small number going up to relatively high values of $\bar{Be}_m \approx 2$. In the panels next to $\bar{Be}_m$, we find the distributions of $\bar{\Psi}_B$ which seem elongated and somewhat non-uniform. However, they are easily explained as the accretion disc is still undergoing a global evolution over the duration of the evaluated time-window ($\Delta T = [3000-4000] r_g/c$). At the beginning ($T = 3000 r_g/c$), we find a mean of $\bar{\Psi}_B \approx 6.5$, while at the end ($T = 4000 r_g/c$) we find a mean of $\bar{\Psi}_B \approx 11.5$.

The last unexplained panels of Fig. 8 are two variations on the orbital velocity $\Omega = u^\phi/u^r$. First, in panel (i.e and r.e), we investigate the ratio between the surface-average within the plasmoid contour ($\bar{\Omega}_{in}$) with the surface-average for a shell directly outside the plasmoid contour ($\bar{\Omega}_{shell}$). The outer edge of the shell corresponds to one-and-a-half times the distance to the central O-point. From this quantity, we can gauge if the plasmoid moves with its surroundings ($\bar{\Omega}_{in}/\bar{\Omega}_{shell} = 1$) or disconnected from it ($\bar{\Omega}_{in}/\bar{\Omega}_{shell} \neq 1$). From the distributions, we find that the mean is consistent with $\bar{\Omega}_{in}/\bar{\Omega}_{shell} = 1$, but there is also significant variance indicating that the plasmoid can move twice as fast or slow with respect to its direct environment. This can be interpreted in a number of ways, which includes, e.g. that the plasmoid is dynamically disconnected from its direct surroundings in the accretion disc.

Second, in panels (i.j and r.j), we evaluate the ratio of $\bar{\Omega}_{in}$ divided by the Keplerian circular orbital velocity in the equatorial plane which is defined as $\Omega_K = (x^{3/2} + a_*)^{-1}$ with $x$ the cylindrical radius (corresponding the horizontal axis in Figs 6 and 7) and $a_* = 0.9375$ the black hole spin parameter. It has been established that MAD discs are sub-Keplerian (Igumenshchev 2008; Porth et al. 2021) which explains the mean of $\bar{\Omega}_{in}/\bar{\Omega}_K \approx 0.8$. Nevertheless, the broad distribution with $0.1 \lesssim \bar{\Omega}_{in}/\bar{\Omega}_K \lesssim 1.3$ indicates the potential for plasmons to be super- or sub-Keplerian, which has interesting observational implications. However, one still has to take into account that our estimate of the Keplerian orbital velocity is somewhat crude as the plasmons have non-zero $u^\theta$ or $u^r$ velocities that break both the circular and equatorial assumption for $\Omega_K$.

Even though the distributions of iM3, rM3, iM4, and rM4 are not explicitly shown, we have confirmed that the general trends described for iM5 and rM5 are consistent with the lower resolution simulations.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** 2D distributions $N(R_S, x)$ of the plasma quantities $x \in \{\bar{\rho}, \bar{\sigma}, \bar{\beta}, \bar{\Psi}_B, \bar{\Omega}_{\text{in}}/\bar{\Omega}_{\text{shell}}, \bar{\epsilon}_{\text{em}}, \bar{\epsilon}_{\text{kin}}, \bar{\epsilon}_{\text{th}}, \bar{B}_{\text{em}}, \bar{\Omega}_{\text{in}}/\Omega_{\text{K}}\}$ as a function of ‘circular’ radius $R_S = \sqrt{S}/\pi$ of the plasmoid for both the iM5 and xM5 cases. While the other parameters have been outlined before, the magnetic Bernoulli factor is defined as $\bar{B}_{\text{em}} = -(h + \sigma/2)u_t$ and the orbital velocity $\Omega = u^\phi/u^t$ with the surface-averaged quantity inside the plasmoid being denoted as $\bar{\Omega}_{\text{in}}$. We stack the distributions as a function of time from 3000 $r_g/c$ to 4000 $r_g/c$ with a 1 $r_g/c$ cadence and divide by the maximum. The green line denotes the mean per width bin that has more than 20 counts total. The left-hand panels denotes iM5 case, while the right-hand panels denotes the xM5 case.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Probability density function $f(R_S/r_g)$ of 'circular' plasmoid radius $R_S$ for both the high-resolution cases $\text{iM5}$ (left) and $\text{rM5}$ (right). All identification takes place within a circle of radius $R = 25r_g$ and we evaluate a time-window of $T \in [3000, 3001, \dots, 3999, 4000] r_g/c$. The rest of the description for Fig. 4 is also applicable here, except now we utilize $R_S$.</figcaption></figure>

The quantitative differences in plasmoid identification rate ($N_P$) will, however, be outlined explicitly for all cases in Section 4.2.4.

### 4.2.3 Plasmoid distribution functions

Fig. 9 displays the probability density function ($f$) of plasmoid radius $R_S$, while Fig. 10 displays the probability density function of plasmoid half-width $w$. We show both distributions to illustrate the general relativistic effects in Fig. 9, while Fig. 10 is straightforwardly compared with the Harris sheet's density function (in Section 4.1.3) and reflects the plasmoids shown in Figs 6 and 7. Interestingly, we recover the power-law indices of $p = 1.88 \pm 0.06$ ($p = 1.90 \pm 0.05$) and $p = 2.15 \pm 0.11$ ($p = 2.09 \pm 0.09$) for $\text{iM5}$ and $\text{rM5}$ in Fig. 9 (10), respectively. These are similar to the results described in Section 4.1.3, which indicates that plasmoid formation is driven by the same principles, even when taking into account the underlying curvature of the space-time. Even though more plasma instabilities and subsequent perturbation are activated (as outlined in Section 4.2.1), we still find scaling laws that is consistent with $p \approx 2$ for the GRMHD simulations. While the onset of magnetic reconnection in the isolated Harris sheet simulations occurs somewhat spontaneously, in GRMHD it is subjected to global dynamics (such as the RTI and KHI) that trigger magnetic reconnection. Although one clearly sees Harris-sheet-like structures forming in GRMHD, they also rapidly fall apart which interestingly does not affect the trends in the density functions. One therefore concludes that the width distributions are robust features of reconnection, no matter how it is triggered.

The identification strategy we effectively employ for the GRMHD simulations is more consistent with studies (i.e. Huang & Bhat-tacharjee 2012; Loureiro et al. 2012; Sironi et al. 2016) that have outflowing boundary conditions as we stop identifying plasmoids for $r > 25 r_g$. Another 'outflowing' boundary lies at the

horizon but the vast majority of plasmoids moves outwards (in the $\hat{r}$ direction) in the jet-disc region. Some plasmoids, typically associated with green contours, even move into the identification domain along the equatorial plane to then exit via the upper or lower identification boundaries. Only a relatively minor fraction of plasmoids is accreted onto the BH and the majority of those are created in close proximity to the BH in the equatorial current sheet.

Close examination of Fig. 9 reveals that plasmoid radius goes all the way up to $R_S \approx 10r_g$. The Cartesian projection equivalent in Fig. 10 shows a maximum radius of $w \approx 2r_g$. These largest plasmoids are visible in Fig. 6. The smallest detected plasmoids have radii corresponding to $R_S \approx 10^{-1}r_g$ (and $w \approx 10^{-2}r_g$). The largest plasmoids seem comparable in size to the 'hotspots' that were used to interpret flares around Sgr A* (Gravity Collaboration 2020b; Vos et al. 2022; Wielgus et al. 2022b). From our simulations, we find that the plasmoids are of sufficient size to give a physical origin to these hotspots. However, currently, we do not explicitly interpret their emission potential, but as plasmoids are typically hot ($p/\rho \gtrsim 1$) and magnetized ($\langle \bar{\sigma} \rangle \gtrsim 0.1$, as per Fig. 8) they are likely to create an emission feature, albeit undetermined if predominantly thermal or non-thermal (Werner et al. 2017; Petropoulou & Sironi 2018). Nevertheless, the occurrence rate of these large, and potentially bright, plasmoids is still quite low. More specifically, for $\text{rM5}$, plasmoids with radii $R_S > 2.5 r_g$ occur at least once (and three times on average) per individually evaluated snapshot (corresponding to 8.2 per cent of all identified plasmoids), while plasmoids with (Cartesian-projected) widths $w > 1 r_g$ are much less common as they occur in only half (51.4 per cent) of the evaluated snapshots (corresponding to 1.8 per cent of all identified plasmoids). This perceived discrepancy is partially due to the space-time curvature (not taken into account for $w$) and the mixing of plasmoid length and width for the $R_S$ quantity. For $\text{iM5}$, the occurrence rates of at least one plasmoid passing the $R_S$ or $w$ criteria are 57.7 and 16.5 per cent (with 6.6 and 2.3 per cent for
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Probability density function $f(w/r_g)$ of plasmoid half-width $w$ for both the high-resolution cases iM5 (left) and rM5 (right). The rest of the description for Fig. 4 is also applicable here. Note that the quantities here do not correctly take into account the space-time curvature, which is the case for Fig. 9.</figcaption></figure>

all identified plasmoids over the entire time window), respectively. Overall, if we take into account the much lower plasmoid count for iM5, we find that the percentages of the two cases are comparable, except for having at least one $w > 1$ plasmoid per evaluated snapshot. This is well-explained, however, in Section 4.2.4.

Lastly, we note that the power-law index $p = -\log f/\log(R_S/r_g)$ is less steep for iM5 than for rM5. We believe this is largely explained by the lower plasmoid number, but we also note that the colours indicate that at later times (more yellow) the plasmoid density function spans more radius (or width) bins and therefore lies slightly lower than at earlier times (dark purple to black). This indicates there is some evolution in the density function as is confirmed in Section 4.2.4. For rM5, we find a relatively consistent density function over time. Next to a potential difference in evolution, we find that a singular linear relation (in log-log space) is not the best description of the downwards trend. Even though close to $p \approx 2$, there is a minor break visible and the slope becomes shallower at $R_S/r_g \approx 4$. As especially the larger plasmoid size bins contain more counts, this naturally pushes $p$ to slightly lower values for iM5. Nevertheless, it is interesting that rM5 indicates a steeper power-law index with $p = 2.15 \pm 0.11$. However, combined with the points raised at the end of Section 4.1.3, we conclude that iM5 and rM5 are consistent with a power law with $p \approx 2$ as more robust claims can not be made without further investigation.

### 4.2.4 Timeseries of plasmoids and horizon penetrating fluxes

Plasmoids form within the accretion disc and are mostly governed by their local plasma conditions. Nevertheless, it is interesting to see how plasmoids, as probes of the accretion disc, connect to the global fluxes on the central black hole. Quantities that are customarily calculated to assess this are the mass accretion rate ($\dot{M}$) and surface-

penetrating magnetic flux ($\Phi_B$), which are defined as (Porth et al. 2019):

$$ \dot{M} = - \int_{0}^{2\pi} \int_{0}^{\pi} \rho u^r \sqrt{-g} \, d\theta \, d\phi, \qquad (25) $$

$$ \Phi_B = \frac{1}{2} \int_0^{2\pi} \int_0^{\pi} |F^{rt}| \sqrt{-g} \, d\theta \, d\phi. \qquad (26) $$

MAD models are known to saturate in horizon-penetrating magnetic flux. This implies that magnetic energy will be building up and will eventually be released in a sudden flux eruption that partly and temporarily halts the accretion flow onto the BH. In 2D simulations, the accretion flow will be stopped completely due to the constraining nature of the setup. The parameter that is used to quantify this behaviour is the so-called MAD parameter $\phi_{BH} = \Phi_B/\sqrt{\dot{M}}$, which corresponds to the normalized magnetic flux. The MAD parameter saturates (in 3D) at $\phi_{BH} \approx 15^2$ (cf. Yuan & Narayan 2014). In our simulations, as shown in Fig. 11, we find that $\phi$ occasionally rises to $\phi_{BH} \sim 120$. This is partly due to the confining nature of the 2D simulation, which allows for a greater accumulation of magnetic flux before an eruption. These findings are consistent with behaviour found for simulations presented in Ripperda et al. (2020). As we used a different adiabatic index $\dot{\hat{\gamma}} = 13/9$ (versus $\dot{\hat{\gamma}} = 4/3$ for Ripperda et al. 2020), we have a thicker disc at initialization which allows for greater accumulation of magnetic flux.

The middle to lower panels (d-f) of Fig. 11 display the number of identified plasmoids ($N_P$) per simulation. While not shown explicitly in the figure, we confirm that plasmoids for either polarity (i.e. magenta and green contours in Figs 6 and 7) are equally abundant. As we already indicated (in Section 4.2.2), a significantly lower number of plasmoids is detected for the ideal simulations than for

<footnote n="2">²Note that this differs by a factor $\sqrt{4\pi}$ from works like Tchekhovskoy et al. (2011), McKinney et al. (2012), and Tchekhovskoy & McKinney (2012). There, $\phi_{BH}$ saturates at $\sim 50$.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Timeseries of the mass accretion rate $\dot{M}$ (panel a), magnetic flux $\Phi_B$ (b), normalized magnetic flux $\phi = \Phi_B/\sqrt{|\dot{M}|}$ (c), number of identified plasmons $N_{\text{plasmons}}$ per simulation (panels d, e, and f), and normalized cross-correlation function (xcorr) between $-\nabla\Phi_B$ and $N_P$ (g). The fluxes are calculated at $2.5 r_g$. We display both the ideal (iM3, iM4, and iM5 in shades of orange) and resistive (rM3, rM4, and rM5 in shades of purple) GRMHD simulations.</figcaption></figure>

<table><caption>**Table 3.** The modulation index $M_Q \equiv \sigma_Q/\mu_Q$ with $\sigma_Q$ and $\mu_Q$ denoting the standard deviation and mean of quantity $Q \in \{M, \Phi_B\}$. This index gives a measure of the variability in the simulations' timeseries.</caption><thead><tr><th>Name</th><th>μ<sub>M</sub></th><th>σ<sub>M</sub></th><th>M<sub>M</sub></th><th>μ<sub>Φ<sub>B</sub></sub></th><th>σ<sub>Φ<sub>B</sub></sub></th><th>M<sub>Φ<sub>B</sub></sub></th></tr></thead><tbody><tr><td>iM3</td><td>6.86</td><td>7.05</td><td>1.03</td><td>55.99</td><td>8.87</td><td>0.16</td></tr><tr><td>iM4</td><td>7.36</td><td>11.46</td><td>1.56</td><td>54.19</td><td>8.93</td><td>0.16</td></tr><tr><td>iM5</td><td>7.29</td><td>17.35</td><td>2.38</td><td>50.45</td><td>10.17</td><td>0.2</td></tr><tr><td>rM3</td><td>6.34</td><td>8.16</td><td>1.29</td><td>59.38</td><td>6.27</td><td>0.11</td></tr><tr><td>rM4</td><td>7.29</td><td>9.23</td><td>1.27</td><td>59.37</td><td>6.55</td><td>0.11</td></tr><tr><td>rM5</td><td>10.78</td><td>17.28</td><td>1.6</td><td>57.79</td><td>7.53</td><td>0.13</td></tr></tbody></table>

the resistive ones where a factor 2–10 difference (in $N_P$) is common. The mechanism that triggers plasmoid formation, via the tearing instability, is not well-defined in ideal MHD simulations, as the ideal resistivity ($\eta_{\text{ide}}$) is resolution dependent (Zhang et al. 2003). This implies that numerical resistivity ($\eta_{\text{ide}} = \eta_{\text{num}}$) is lower close to the BH due to the MKS coordinate system and is significantly smaller than $\eta_{\text{ris}} = \eta$ ($\eta_{\text{num}} < <\eta$). Overall, the tearing instability is triggered less often, due to the relatively lower resistivity, and less reliably as it is determined by (stochastic) numerical effects. Visually, the ideal simulations are significantly calmer, which is explained by the suppression of the MRI after the initial few thousand time-steps. Starting from $T \approx 3700 r_g/c$, however, a sudden increase in plasmoid formation rate is visible, which roughly corresponds to the state shown in Fig. 6 for iM5. After this ‘flaring’ event, the rate at which plasmons are created is somewhat increased (except for iM5).

The resistive simulations possess a strikingly constant number of plasmons ($N_P$) indicating a steady rate of plasmoid formation. As the MRI is also suppressed for the resistive simulations, we can assume that the tearing instability is a sufficient perturbation in itself to keep plasmoid formation up. To get an indication of how the flux eruptions could contribute to this process, we verified if there are significant changes in the modulation index $M_Q \equiv \sigma_Q/\mu_Q$ (see EHTC et al. 2022b, and description therein), with $\sigma_Q$ and $\mu_Q$ being the standard deviation and mean, respectively, over time-interval $\Delta T = |3000-4000| r_g/c$ for a given quantity $Q$. The modulation index can be universally applied to a given timeseries to assess the level of variability. We calculate the modulation index for both the accretion rate $\dot{M}$ and magnetic flux $\Phi_B$ that penetrate the spherical shell at $r = 2.5 r_g$. An overview of the simulation modulation indices is listed in Table 3. There seems to be little difference between $M_M$ or $M_{\Phi_B}$ for the ideal and resistive simulations. This is surprising as $N_P$ indicates a more turbulent disc for the resistive cases, as this would give rise to the greater plasmoid count. Nevertheless, one is not able to ascertain this directly from the shell-penetrating fluxes. Another consequence of setting a fixed resistivity is that there is a fixed length scale (i.e. width of the current sheet) that determines when the tearing instability is triggered. When this length-scale is sufficiently resolved, one finds consistent results starting from a certain critical resolution and upwards. It is therefore interesting that we see this being verified in the panels (d)-(f). For rM3, the lowest resolution case, we find that the mean plasmoid count $\langle N_P \rangle \sim 75$. While for the higher resolution cases rM4 and rM5, we find $\langle N_P \rangle \sim 100$. As we find converging plasmoid numbers for both resolution cases, we conclude that the current sheet thickness (set by $\eta = 5 \times 10^{-5}$), within the 25 $r_g$ domain, is fully resolved starting from a resolution of 4096².

For the last panel (g), we investigate the relation between a (MAD) flux eruption (i.e. characterized by a drop in horizon-penetrating magnetic flux $\Phi_B$) and the detected plasmoid count ($N_P$). We cross-
correlate the plasmoid number ($N_p$) with the negative gradient of the magnetic shell-penetrating flux ($-\nabla\Phi_B$) and find a positive relation for most cases. This indicates that the perturbing effect of the flux eruption brings about a period of increased plasmoid formation after a characteristic time-delay. Except for $\tau M_3$, which is the uncorrelated component on the background (in lightest purple), we find a clear correlation between the most prominent peak in $N_p$ and a decrease in $\Psi_B$. The maxima of the correlation function coincides with beginning of a drop in the magnetic flux function and are denoted by vertical dashed line in their corresponding panels. This is a consistent trend as long as one has a clear flux eruption, which also explains the uncorrelated $\tau M_3$ results as there is no clear decrease in $\Phi_B$ present. For $iM_5$ at $T \approx 3780 \ r_g/c$, the flux eruption is rather large as is indicated by the decrease in $\Phi_B$, which has pushed the maximum xcorr($-\nabla\Phi_B, N_p$) further to the right. Just before the flux eruption, we find that an increase (of several factors) in $\Phi_B$ after which it drops quickly. The positive correlation is naturally explained by the fact that the flux eruption, which is accompanied by the temporary halting of the accretion flow, is a significant perturbation to the accretion disc that is able to initiate reconnection in numerous places. Even though this general picture applies, we find that the dynamics are likely also stochastic in nature as the $\tau M_4$ case displays different behaviour with a drop in $N_p$ directly after the flux eruption. This is in part explained by our identification strategy which only identifies plasmoids within 25 $r_g$ and as the disc has receded during the flux eruption the domain in which plasmoids are detected also shrinks and effectively delays the peak in $N_p$. Additionally, the shell-penetrating magnetic flux ($\Phi_B$) only shows a relatively minor depression which indicates a relatively minor flux eruption and subsequent perturbation of the disc structure. So, to summarise, one expects a reaction on the plasmoid formation rate following a flux eruption, which tends to increase the plasmoid count as it perturbs the accretion disc and subsequently triggers reconnection.

# 5 DISCUSSION

In this section, we will discuss our results following in the context of earlier works following four main points; (i) direct comparison to GRMHD-related plasmoid detection methods, (ii) specifics from our simulation library, (iii) implication for 3D results, (iv) effects of resistivity, and (v) a discussion of the flaring potential of plasmoids.

## 5.1 GRMHD plasmoid detection

Comparison with earlier works that have identified plasmoid structures in GRMHD (Nathanail et al. 2020; Jiang et al. 2023) suggests that the approach outlined in this work finds $5-10 \times$ more plasmoids. Both aforementioned works utilize the Bernoulli factor ($Be = -hu_t$) as underlying identification medium and use a canny-edge detection algorithm on a Gaussian blurred segment (as provided by the scikit Python package). We have made initial attempts with this proposed method. However, as we did not reach the desired efficacy or fidelity, we started the development of algorithm outlined in this work. Overall, we typically find $5-10 \times$ more plasmoids than the previously mentioned works, which are not all attributed to the detection method difference. Other potential causes for the discrepancy can be the identification medium, resolution, simulation configuration, and the inherent differences between resistive and ideal GRMHD. A number of these points will be discussed in detail in the following paragraphs.

Regarding the choice of identification medium, we opted for the magnetic flux function $\Psi_B$ as it naturally identifies circular magnetic

field structures. When we compare this criterion with the utilization of the Bernoulli factor $Be$, then it becomes clear from the results in this work that not all plasmoids are unbound as demonstrated in Fig. 8. One is likely to miss the plasmoids created in the equatorial plane as those tend to be bound (as was also pointed out in Jiang et al. 2023). There are also clear advantages to using $Be$, because one can apply well-established image-recognition algorithms if one is able to increase the contrast (i.e. only show a limited colour range) to which the $Be$ factor lends itself well. Nevertheless, this comes with the cost that one can only identify a relatively small subset of the entire plasmoid population.

## 5.2 Simulation library

When visually comparing our simulation to those of Ripperda et al. (2020), with a highest resolution of 6144×3072 with respect to our 8192×8192, then we infer that the number of plasmoids does not differ significantly based on the presented figures, except perhaps at the smallest scales. More importantly, one may even draw the conclusion that SANE simulations produce clearer and more abundant plasmoid structures. Nathanail et al. (2020) utilize an initial single dipolar loop up to intricate multipolar initial magnetic field configurations with an evolution that can be described as SANE-like (with low $\phi_{BH} \sim 2$). Especially, the multipolar configurations are expected to produce a lot of plasmoids, as is confirmed in their Fig. 6. However, they do not show any statistics. This is done, however, in Jiang et al. (2023) using the same methodology, but their configuration has a multipolar initial magnetic field and evolves to be heavily magnetized (i.e. MAD-like). The evolution is very chaotic and consistent with MAD but only relatively few plasmoids are visible indicating that the lower resolution (up to 4096 × 2048) and identification technique are likely to play a role. It is important to note that those simulations all utilized an ideal MHD description, so we only compare them to the $iM_3, iM_4$, and $iM_5$ cases. The differences between resistive and ideal GRMHD will be discussed in detail in a following paragraph (i.e. ‘resistivity’).

## 5.3 3D

How applicable are 2D results to a 3D reality? A number of arguments come into play here. The plasmoids in our simulations describe predominantly elliptical (close to circular) structures and have long merging chains. This is in part explained by the confined nature of the 2D simulations. As, due to this confining nature, plasmoids have a greater probability to interact and merge, they are likely to become larger. If one were to add an additional dimension (in $\phi$), one significantly complicates the situation. First, the plasmoid morphology would change and gain the resemblance of a flux rope. Secondly, the chance for interaction would decrease significantly as it is simply less likely to come across another flux rope. Thirdly, the definition of flux ropes coalescence is difficult as they likely merge in a single place but not in its entirety. These points are clearly demonstrated for the 3D equivalent of the (Harris) current sheet (see e.g. Kagan, Milosavljević & Spitkovsky 2013; Cerutti et al. 2014; Sironi & Spitkovsky 2014; Zhou et al. 2020). There, one finds complex behaviour of and interaction between flux ropes that is partially due to the presence of the kink instability (e.g. Bromberg et al. 2019; Davelaar et al. 2020; Werner & Uzdensky 2021), which is absent in axisymmetric simulations. For high-resolution 3D GRMHD simulations, some evidence for the presence of plasmoids, or flux ropes, was presented in Ripperda et al. (2022). Nevertheless,
the typical appearance and how much it stands out with respect to its
environment is relatively unknown in 3D.

## 5.4 Resistivity

In essence, setting an explicitly resistivity ($\eta$) allows for consistently resolving the underlying current sheets in the simulation, which in ideal (GR)MHD is ill-defined as it is numerically determined and therefore has a stochastic (and coordinate-dependent) component. As is outlined in Section 4.2.4, there is a clear discrepancy between the resistive and ideal simulations. While the former has a relatively consistent plasmons number $N_p \sim 100$, the latter has a non-flaring count comparable to $N_p \sim 10$. So, even though these discrepancies were expected, they were not verified in regard to plasmoid count till now. In part it can be a selection effect as the ideal simulation(s) entered a 'quiet' phase with few perturbations to the disc structure, but it is interesting this does not happen for the resistive case. However, in light of recent finding by the Event Horizon Telescope Collaboration (EHTC et al. 2022b), where was pointed out that the (ideal) GRMHD simulations produce too variable emission signatures, one can draw the tentative conclusion that the variability could be even more significant for resistive GRMHD simulations. Additionally, the physical interpretation of resistivity is that it is a proxy for kinetic effects, which are simulated self-consistently with PIC methods, but to assess what is the 'correct' resistivity for our physical scenario is a non-trivial question (Selvi et al. 2023). A rigorous (GRMHD) study including several resistivities is therefore needed to make more robust claims, but this is rather computationally expensive as one needs to assure that the current sheets are well-resolved.

## 5.5 Misidentification

For the approach outlined in this work, we are indiscriminate as to what properties the plasmoid should contain, except that it should correspond to a circular magnetic field geometry. Even though this allows us to get a rather complete distribution, it is slightly sensitive to misclassifications, which happens mainly for overly dense region. This is explained by the sensitivity of both the local extrema finder and the watershed algorithm – even though it is only a minor deviation from the background, it is treated as if it is a plasmoid. Overall, this happens only rarely. What occurs more often is that plasmons that are in close vicinity to each other are grouped as they have very similar $\Psi_B$ signatures. Except that this diminishes the detected plasmoid count, it does not influence the surface-averaged quantities (and distributions) as they still probe the plasmoid structure. As with all identification problems, the difficulty lies in finding a strategy that is able to bridge the various length-scales while not picking up on erroneous features. This is largely determined by the blurring layer, which dictates the minimal size-scale to which one is sensitive and gives a handle on how much fine-structure one wants to include. As the large plasmons tend to have a lot of fine-structure, one should apply a more aggressively blurring strategy. Even though our algorithm is accurate, it is by no means computationally fast, even despite parallelization attempts that will be intensified in the future. At present, we do not give an exact number of misclassifications, but we are able to find a few in most snapshots while the vast majority (of $\mathcal{O}(100)$) is classified correctly. The number of plasmons that were not classified is also of order $\mathcal{O}(1)$ and are predominantly caused by numerical instabilities in the contour-finding step of the algorithm that typically occur for relatively unclear 'plasmoid' structures.

## 5.6 Flaring potential

While we started this paper by talking about plasmons as a potential source for flares, it is nevertheless difficult to make direct emission interpretations. The main reason for this is that the emission properties of plasmons in the BH accretion environment are still unknown, especially as one would expect a significant non-thermal contribution. The utilization of a thermal synchrotron proxy (as utilized in, e.g. Porth et al. 2019) would therefore likely give an unrealistic picture. Ripperda et al. (2020) gave estimates of the synchrotron emission and its potential to explain flares and our estimates are of the same order. Nevertheless, it would be beneficial to conduct a full radiative transfer study to accurately assess the flaring potential of plasmons including a non-thermal electron population or reconnection-dedicated electron temperature description (Rowan, Sironi & Narayan 2017). This is an interesting avenue to pursue in the future, as it is possible to track the plasmoid's location with the algorithm.

# 6 CONCLUSIONS

We have been able to identify plasmons in highly turbulent accretion discs surrounding SMBHs with a higher fidelity than has been achieved before, which allows for creating complete time-series and distributions with sufficiently large numbers to assess the statistics. Additionally, we have verified our methodology with a set of previously well-investigated Harris current sheet simulations and found they are consistent with findings from previous studies (Uzdensky et al. 2010; Huang & Bhattacharjee 2012; Loureiro et al. 2012; Sironi et al. 2016). Interestingly, the scaling laws (outlined in Sections 4.1.3 and 4.2.3) for both the Harris sheet and the GRMHD simulations are very similar, which indicates that plasmoid formation in the more complex accretion disc environment does not differ fundamentally from the Harris sheet picture. Utilization of the newly developed algorithm has enabled us to better study the plasmoid population within MAD accretion discs, and has clearly laid bare discrepancies in plasmoid occurrence rates between ideal and resistive (GR)MHD that warrant further investigation with a more systemic study that includes other accretion scenarios (e.g. SANEs).

The typical plasmoid in a MAD GRMHD simulation is equally dense and somewhat undermagnetized with respect to their surrounding, moves with its surroundings, and is likely to be unbound according to the Bernoulli criterion. Nevertheless, this behaviour describes the averages of distributions and does not describe the deviations which occur frequently. Especially for the orbital velocities and boundedness of the plasmons, one finds large spreads in the distributions. This indicates that plasmons can both occur as superior sub-Keplerian features, which is currently still an active point of investigation within the community. Magnetic saturation at the BH event horizon produces flux tubes in violent events that (partially) push back the accretion flow for MAD simulations. Even though this is one of the leading theories to explain flares around SMBHs (Dexter et al. 2020; Porth et al. 2021), they are established to orbit with strongly sub-Keplerian velocities, which is at odds with some observations. The formation of plasmons is, therefore, still a candidate for explaining both Keplerian (Gravity Collaboration 2018b, 2020a) and super-Keplerian (Matsumoto, Chan & Piran 2020) near-infrared observation of flares around Sgr A*. More specifically, we regularly recover plasmoid sizes that are comparable to the hotspots that were used to interpret flares at both NIR- and mm-wavelengths (Gravity Collaboration 2020a, b; Vos et al. 2022; Wielgus et al. 2022a). Also, as outlined in Section 4.2.4, flux eruptions (corresponding
to a decrease in horizon-penetrating magnetic flux $\Phi_B$) and plasmoid formation are likely strongly correlated with one another, indicating that flux eruptions can act as an instigator of magnetic reconnection. Both the flux tube and plasmoid (or flux rope) pictures do therefore not have to be mutually exclusive, but rather have a complementary co-existence.

Lastly, we would like to point out that the identification algorithm is much more universally applicable as its function can be well-characterized as a 'closed contour-detector around local extrema'. So, in the future, we are planning to apply our methodology to investigate the 3D structure of plasmoids, flux ropes, and flux tubes in SMBH accretion scenarios. It would also lend itself well to other MHD or PIC identification applications, such as shearing or turbulent box simulations.

## ACKNOWLEDGEMENTS

We thank Bart Ripperda, Jordy Davelaar, Fiorenze Stoppa, Alejandra Jimenez-Rosales, and Aristomenis Yfantis for the helpful discussions and comments on the manuscript. We extend particular thanks to the anonymous referee, who provided exceptionally thorough and thoughtful comments that substantially improved this work. JV acknowledges support from the Dutch Research Council (NWO) supercomputing grant No. 2021.013. MM acknowledges support by the NWO grant No. OCENW.KLEIN.113 and support by the NWO Science Athena Award. BC acknowledges the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 863412). HO acknowledges funding from Radboud University through a Virtual Institute of Accretion (VIA) post-doctoral fellowship from the Netherlands Research School for Astronomy (NOVA).

Software used to create the results in this work; BHAC (Porth et al. 2017; Olivares et al. 2019), PYTHON (Van Rossum & Drake 2009), NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), GNU PARALLEL (Tange 2022).

## DATA AVAILABILITY

The data used for this work will be shared following a reasonable request to the authors. The plasmoid identification algorithm will be made publicly available in following repository; https://github.com/JesseVos/Plasmoid_Finder, in due time, but will also be shared following a reasonable request before that time.

## REFERENCES

* Avara M. J., McKinney J. C., Reynolds C. S., 2016, MNRAS, 462, 636

* Bacchini F., Arzamasskiy L., Zhdankin V., Werner G. R., Begelman M. C., Uzdensky D. A., 2022, ApJ, 938, 86

* Baganoff F. K. et al., 2001, Nature, 413, 45

* Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214

* Banesh D., Lo L.-T., Kilian P., Guo F., Hamann B., 2020, preprint (arXiv:2010.01959)

* Begelman M. C., Blandford R. D., Rees M. J., 1984, Rev. Mod. Phys., 56, 255

* Begelman M. C., Scepi N., Dexter J., 2022, MNRAS, 511, 2040

* Beucher S., Meyer F., 2018, Mathematical Morphology in Image Processing. CRC Press, p. 433

* Bhattacharjee A., Huang Y.-M., Yang H., Rogers B., 2009, Phys. Plasmas, 16, 112102

* Birn J., Hesse M., 2001, J. Geophys. Res., 106, 3737

* Birn J. et al., 2001, J. Geophys. Res., 106, 3715

* Borgogno D., Grasso D., Achilli B., Romé M., Comisso L., 2022, ApJ, 929, 62

* Bransgrove A., Ripperda B., Philippov A., 2021, Phys. Rev. Lett., 127, 055101

* Broderick A. E., Loeb A., 2005, MNRAS, 363, 353

* Broderick A. E., Loeb A., 2006, MNRAS, 367, 905

* Bromberg O., Singh C. B., Davelaar J., Philippov A. A., 2019, ApJ, 884, 39

* Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2012, ApJ, 754, L33

* Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2013, ApJ, 770, 147

* Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2014, ApJ, 782, 104

* Cerutti B., Philippov A., Parfrey K., Spitkovsky A., 2015, MNRAS, 448, 606

* Chen A. Y., Beloborodov A. M., 2014, ApJ, 795, L22

* Comisso L., Sironi L., 2019, ApJ, 886, 122

* Comisso L., Lingam M., Huang Y.-M., Bhattacharjee A., 2016, Phys. Plasmas, 23

* Crinquand B., Cerutti B., Philippov A., Parfrey K., Dubus G., 2020, Phys. Rev. Lett., 124, 145101

* Crinquand B., Cerutti B., Dubus G., Parfrey K., Philippov A., 2021, A&A, 650, A163

* Davelaar J., Philippov A. A., Bromberg O., Singh C. B., 2020, ApJ, 896, L31

* De Villiers J.-P., Hawley J. F., Krolik J. H., 2003, ApJ, 599, 1238

* Dexter J. et al., 2020, MNRAS, 497, 4999

* Do T. et al., 2019, ApJ, 882, L27

* EHTC et al., 2019a, ApJ, 875, L1

* EHTC et al., 2019b, ApJ, 875, L5

* EHTC et al., 2021, ApJ, 910, L13

* EHTC et al., 2022a, ApJ, 930, L12

* EHTC et al., 2022b, ApJ, 930, L16

* Eckart A. et al., 2004, A&A, 427, 1

* El Mellah I., Cerutti B., Crinquand B., Parfrey K., 2022, A&A, 663, A169

* El Mellah I., Cerutti B., Crinquand B., 2023, A&A, 677, A67

* Fishbone L. G., Moncrief V., 1976, ApJ, 207, 962

* Genzel R., Schödel R., Ott T., Eckart A., Alexander T., Lacombe F., Rouan D., Aschenbach B., 2003, Nature, 425, 934

* Ghez A. M. et al., 2008, ApJ, 689, 1044

* Gillessen S., Eisenhauer F., Trippe S., Alexander T., Genzel R., Martins F., Ott T., 2009a, ApJ, 692, 1075

* Gillessen S., Eisenhauer F., Fritz T. K., Bartko H., Dodds-Eden K., Pfuhl O., Ott T., Genzel R., 2009b, ApJ, 707, L114

* Gillessen S. et al., 2017, ApJ, 837, 30

* Goedbloed J. P., Poedts S., Keppens R., 2010, Advanced Magnetohydrodynamics: With Applications to Laboratory and Astrophysical Plasmas, Cambridge Univ. Press, Cambridge

* Gravity Collaboration, 2018a, A&A, 615, L15

* Gravity Collaboration, 2018b, A&A, 618, L10

* Gravity Collaboration, 2019, A&A, 625, L10

* Gravity Collaboration, 2020a, A&A, 635, A143

* Gravity Collaboration, 2020b, A&A, 643, A56

* Guépin C., Cerutti B., Kotera K., 2020, A&A, 635, A138

* Guo F., Li H., Daughton W., Li X., Liu Y.-H., 2016, Phys. Plasmas, 23, 055708

* Hakobyan H., Philippov A., Spitkovsky A., 2019, ApJ, 877, 53

* Hakobyan H., Petropoulou M., Spitkovsky A., Sironi L., 2021, ApJ, 912, 48

* Harris C. R. et al., 2020, Nature, 585, 357

* Hillier A., 2019, Phys. Plasmas, 26, 082902

* Huang Y.-M., Bhattacharjee A., 2012, Phys. Rev. Lett., 109, 265002

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Iğumenshchev I. V., 2008, ApJ, 677, 317

* Iğumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042

* Jiang H.-X., Mizuno Y., Fromm C. M., Nathanail A., 2023, MNRAS, 522, 2307

* Kadowaki L. H. S., Pino E. M. D. G. D., Stone J. M., 2018, ApJ, 864, 52

* Kagan D., Milosavljević M., Spitkovsky A., 2013, ApJ, 774, 41

* Kagan D., Sironi L., Cerutti B., Giannios D., 2015, Space Sci. Rev., 191, 545

* Kagan D., Nakar E., Piran T., 2016, ApJ, 826, 221
Keppens R., Porth O., Galsgaard K., Frederiksen J. T., Restante A. L., Lapenta G., Parnell C., 2013, *Phys. Plasmas*, 20, 092109

Loureiro N. F., Schekochihin A. A., Cowley S. C., 2007, *Phys. Plasmas*, 14, 100703

Loureiro N. F., Samtaney R., Schekochihin A. A., Uzdensky D. A., 2012, *Phys. Plasmas*, 19, 042303

Lyutikov M., Sironi L., Komissarov S., Porth O., 2016, preprint (arXiv:1603.05731)

Marshall M. D., Avara M. J., McKinney J. C., 2018, *MNRAS*, 478, 1837

Matsumoto T., Chan C.-H., Piran T., 2020, *MNRAS*, 497, 2385

McKinney J. C., Gammie C. F., 2004, *ApJ*, 611, 977

McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, *MNRAS*, 423, 3083

Misner C. W., Thorne K. S., Wheeler J. A., 1973, Gravitation. W. H. Freeman and Company, United States

Nalewajko K., Uzdensky D. A., Cerutti B., Werner G. R., Begelman M. C., 2015, *ApJ*, 815, 101

Narayan R., Iğumenshchev I. V., Abramowicz M. A., 2003, *PASJ*, 55, L69

Narayan R., Sá dowski A., Penna R. F., Kulkarni A. K., 2012, *MNRAS*, 426, 3241

Nathanil A., Fromm C. M., Porth O., Olivares H., Younsi Z., Mizuno Y., Rezzolla L., 2020, *MNRAS*, 495, 1549

Olivares H., Porth O., Davelaar J., Most E. R., Fromm C. M., Mizuno Y., Younsi Z., Rezzolla L., 2019, *A&A*, 629, A61

Parfrey K., Philippov A., Cerutti B., 2019, *Phys. Rev. Lett.*, 122, 035101

Parker E. N., 1957, *J. Geophys. Res.*, 62, 509

Petropoulou M., Sironi L., 2018, *MNRAS*, 481, 5687

Philippov A. A., Spitkovsky A., 2018, *ApJ*, 855, 94

Porquet D., Predehl P., Aschenbach B., Grosso N., Goldwurm A., Goldoni P., Warwick R. S., Decourchelle A., 2003, *A&A*, 407, L17

Porth O., Xia C., Hendrix T., Moschou S. P., Keppens R., 2014, *ApJS*, 214, 4

Porth O., Olivares H., Mizuno Y., Younsi Z., Rezzolla L., Moscibrodzka M., Falcke H., Kramer M., 2017, *Computat. Astrophys. Cosmol.*, 4, 1

Porth O. et al., 2019, *ApJS*, 243, 26

Porth O., Mizuno Y., Younsi Z., Fromm C. M., 2021, *MNRAS*, 502, 2023

Qian Q., Fendt C., Noble S., Bugli M., 2017, *ApJ*, 834, 29

Rezzolla L., Zanotti O., 2013, Relativistic Hydrodynamics. Oxford Univ. Press, Oxford

Ripperda B., Porth O., Xia C., Keppens R., 2017, *MNRAS*, 467, 3279

Ripperda B. et al., 2019a, *ApJS*, 244, 10

Ripperda B., Porth O., Sironi L., Keppens R., 2019b, *MNRAS*, 485, 299

Ripperda B., Bacchini F., Philippov A. A., 2020, *ApJ*, 900, 100

Ripperda B., Liska M., Chatterjee K., Musoke G., Philippov A. A., Markoff S. B., Tchekhovskoy A., Younsi Z., 2022, *ApJ*, 924, L32

Rowan M. E., Sironi L., Narayan R., 2017, *ApJ*, 850, 29

Sadowski A., Narayan R., Penna R., Zhu Y., 2013, *MNRAS*, 436, 3856

Selvi S., Porth O., Ripperda B., Bacchini F., Sironi L., Keppens R., 2023, *ApJ*, 950, 169

Servidio S., Matthaeus W. H., Shay M. A., Cassak P. A., Dmitruk P., 2009, *Phys. Rev. Lett.*, 102, 115003

Servidio S., Matthaeus W. H., Shay M. A., Dmitruk P., Cassak P. A., Wan M., 2010, *Phys. Plasmas*, 17, 032315

Sironi L., Spitkovsky A., 2014, *ApJ*, 783, L21

Sironi L., Giannios D., Petropoulou M., 2016, *MNRAS*, 462, 48

Sweet P. A., 1958, in Lehnert B. ed., Proc. IAU Symp. 6, Electromagnetic Phenomena in Cosmical Physics. Kluwer, Dordrecht, p. 123

Takamoto M., 2013, *ApJ*, 775, 50

Tange O., 2022, GNU Parallel 20220722 ('Roe vs Wade'). https://doi.org/10.5281/zenodo.6891516

Tchekhovskoy A., McKinney J. C., 2012, *MNRAS*, 423, L55

Tchekhovskoy A., Narayan R., McKinney J. C., 2011, *MNRAS*, 418, L79

Uzdensky D. A., 2022, *J. Plasma Phys.*, 88, 905880114

Uzdensky D. A., Loureiro N. F., 2016, *Phys. Rev. Lett.*, 116, 105003

Uzdensky D. A., Loureiro N. F., Schekochihin A. A., 2010, *Phys. Rev. Lett.*, 105, 235002

Van Rossum G., Drake F. L., 2009, Python 3 Reference Manual. CreateSpace, Scotts Valley, CA

Vincent L., Soille P., 1991, *IEEE Trans. Pattern Anal. Mach. Intell.*, 13, 583

Vos J., Močibrodzka M. A., Wielgus M., 2022, *A&A*, 668, A185

Vos J., Davelaar J., Olivares H., Brinkerink C., Falcke H., 2023, preprint (arXiv:2310.16938)

Werner G. R., Uzdensky D. A., 2017, *ApJ*, 843, L27

Werner G. R., Uzdensky D. A., 2021, *J. Plasma Phys.*, 87, 905870613

Werner G. R., Uzdensky D. A., Begelman M. C., Cerutti B., Nalewajko K., 2017, *MNRAS*, 473, 4840

Wielgus M. et al., 2022a, *A&A*, 665, L6

Wielgus M. et al., 2022b, *ApJ*, 930, L19

Winarto H. W., Kunz M. W., 2022, *J. Plasma Phys.*, 88, 905880210

Witzel G. et al., 2021, *ApJ*, 917, 73

Yuan F., Narayan R., 2014, *ARA&A*, 52, 529

Zenitani S., Hoshino M., 2001, *ApJ*, 562, L63

Zenitani S., Hoshino M., 2008, *ApJ*, 677, 530

Zhang Y.-T., Shi J., Shu C.-W., Zhou Y., 2003, *Phys. Rev. E*, 68, 046709

Zhdankin V., Uzdensky D. A., Perez J. C., Boldyrev S., 2013, *ApJ*, 771, 124

Zhou M., Loureiro N. F., Uzdensky D. A., 2020, *J. Plasma Phys.*, 86, 535860401

## APPENDIX A: RESOLUTION CONVERGENCE FOR GRMHD SIMULATIONS

We have performed our simulations at three effective resolution levels, from lowest to highest; 2048 × 2048 (iM3 and rM3), 4096 × 4096 (iM4 and rM4), and 8192 × 8192 (iM5 and rM5). These correspond to the third, fourth, and fifth AMR level, which we will use for reference. In principle, the current sheet are well-resolved starting from the fourth level, which is consistent with the plasmoid number findings in Fig. 11. Nevertheless, it is important to note that only relatively short periods, of 1000 $r_g/c$, have been run at the fourth and fifth level. These simulations have been started from the third level snapshot at 2900 $r_g/c$, but then the resolution is increased up to the desired level. After a period where the simulation adapts to the new resolution level, we start evaluating the window $T \in [3000, 4000] r_g/c$. Next to the analysis described in the main text, it would be interesting to see how the structure changes as a function of resolution level. Therefore, we have calculated the time-averaged profiles over the aforementioned time-window for $\Psi_B$, density $\rho$, and magnetization $\sigma$ for all cases where we are especially interested in the difference with respect to the highest resolution case.

Figs A1 and A2 display the magnetic flux function $\Psi_B$ and the absolute relative difference (ARD) between the various resolution levels of the ideal and resistive simulations. Starting with the ideal simulations, we predominantly find differences in the jet regions. The inner jet region, near the axis, is dominated by numerical floor violations and does therefore not have a physical origin. In the jet sheath, the transition regions between disc and jet, we do find most of the activity and difference between the various resolution levels. Interestingly, for the ideal cases, most of the variability occurs in the upper ($x > 0 r_g$) region, while the bottom jet-sheath acshows little variability. This once again confirms that the evaluated ideal case is atypically quiet. For the resistive cases, we find a similar structure with the highest values within the disc equatorial plane that then drops off the further you move away. However, the jet sheath does have relatively higher flux function values. The most striking difference is the gigantic plasmoid that lies at $x \approx 45 r_g$ on the equatorial plane. This is likely a remnant of the initial poloidal loop at initialization. Overall, we find much variability and differences between the various resolution cases, which indicates more activity overall (as was established throughout the main text). It is also important to note that the maximal differences are $\gtrsim 6$ per cent, which is quite small and reasonable when compared to the other quantities. This once again confirm that the flux function is a very
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Time-averaged magnetic flux function $\Psi_B$ over time interval $T \in [3000, 4000]$ $r_g/c$ for the iM5 (top left), iM4 (top middle), and iM3 (top right) cases. The bottom panels show the absolute relative difference (ARD) between the iM5 and the cases in the panels above.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure A2.** The description of Fig. A1 applies here as well, except here we show the resistive cases; rM5, rM4, and rM3.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure A3.** The description of Fig. A1 applies here as well, except here we show the density $\rho$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure A4.** The description of Fig. A1 applies here as well, except here we show the density $\rho$ for the resistive cases; rM5, rM4, and rM3.</figcaption></figure>
suitable choice in identification medium as it is not very variable,
which makes identification difficult.

That leaves the time-averaged density ($\rho$) results in Figs A3 and A4. We point out that the ARDs go up to 20 per cent, which signifies that the differences are considerably larger. This is in part due to the nature of the density itself as it tends to be small. Nevertheless, we find that the results are consistent with what is shown for the $\Psi_B$ maps. For the resistive cases, we find a lot of activity in both the equatorial plane (up to and concentrated around the giant plasmoid) and the jet sheath. When compared to the ideal cases, we find that especially the equatorial current sheet activity is low. This further outlines the clear differences between the ideal and resistive cases. Nevertheless, we should be wary to take this as a

general result, as it could very well be subject to selection effects.

We already noted that the ideal cases come across as atypically quiet, which may have been the result of an unfortunate choice in time window. Additionally, it is likely coincidental that a large equatorial plasmoid was created for the resistive cases, which has probably enhanced the resistive simulations' variability. So, in the future, it would be worthwhile to undertake a more systematic study of resistive GRMHD simulations to see if the equatorial plasmoid is a common occurrence. Nevertheless, it is likely (partially) related to the confining nature of 2D simulations, as we commented before.

This paper has been typeset from a TeX/LATEX file prepared by the author.