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

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 plasmoids could therefore observationally manifest themselves as hot spots that 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 plasmoids, we determine the 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. Interestingly, 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 behaviour is similar. Plasmoid occurrence 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 a decrease in black-hole-horizon-penetrating magnetic flux. These results demonstrate the efficacy of the newly developed algorithm which has enabled an extensive quantitative analysis of plasmoid formation for black hole accretion simulations.


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 center 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  ≈ 4 × 10 6  ⊙ and lies at a distance of  ≈ 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. 2009aGillessen et al. ,b, 2017;;Gravity Collaboration et al. 2018a, 2019;Do et al. 2019).At sub-mm / mm wavelengths, Sgr A* is known to be a stochastically (O (10%) over hours) variable source that 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-wavelenghts 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).Previously, this ★ E-mail: jt.vos@astro.ru.nl was also established in the NIR band (Gravity Collaboration et al. 2018b).The physical mechanism that causes these flares is currently not well-understood, but a number of working theories associate them with strongly magnetized anisotropies in the accretion flow (Broderick & Loeb 2005, 2006;Gravity Collaboration et al. 2020b;Dexter et al. 2020;Porth et al. 2021;Vos et al. 2022Vos et al. , 2023;;Ripperda et al. 2022).
One such scenario that may explain these flares and the formation of hot spots is the formation of plasmoids as part of a magnetic reconnection event (e.g., Ripperda et al. 2020;Ripperda et al. 2022;El Mellah et al. 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 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 with plasmoid formation.To model the plasmoid-unstable regime, one has to adopt a numerical approach via particle-in-cell (PIC) or magnetohydrodynamical (MHD) simulations.Both methods will be outlined in detail in the following paragraphs.
General relativistic 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).Almost exclusively, one assumes an ideal GRMHD description that is an inadequate framework for capturing magnetic reconnection and the formation of plasmoids (Ripperda et al. 2020), as the plasmoid-instability is triggered due to numerical limits rather than consistently resolving the underlying current sheet.Resistive GRMHD 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 () in the simulations.While the physical resistivity is likely spatially and temporally variable, 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) and is more akin to the methods employed in PIC studies by Sironi et al. (2016); Hakobyan et al. (2019Hakobyan et al. ( , 2021)).However, as in a fluid MHD description one does not have the luxury of individual particle trajectories, we apply our analysis fully in post-processing which gives it more flexibility.Using 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 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 paper is structured as follows.An in-depth description of the methods we use to simulate and identify these features are outlined in Section 2. The results and their interpretation are presented in Section 3. The discussion and conclusion can be found in Sections 4 and 5.

METHODS
In the following sections, we will describe the algorithm that identifies the plasmoids and outline the two setups we investigate.

Relativistic MHD primer: ideal and resistive
The plasma flows of 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), which solves the (resistive; Ripperda et al. 2019a) MHD equations in stationary spacetimes.These equations are defined as; where ∇  denotes 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.BHAC is a versatile code that sets the speed of light  to unity and utilises Lorentz-Heaviside units, which effectively incorporates the √ 4 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 is denoted by Note that the resistivity is denoted by  = 1/  where   is the conductivity.While in resistive MHD the electric field (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×B, also known as the "frozen-in condition").Effectively, one assumes the plasma to be perfectly conducting (  → ∞ ⇒  → 0) in the ideal MHD limit, which is a useful and macroscopically valid approximation in large parts of the accretion disc domain but not when it comes to the formation of plasmoids and other non-ideal effects.More specifically, the resistivity  is not exactly zero in the ideal case (except for infinite resolution), but rather determined numerically by the underlying resolution (or cell size Δ) which implies that  ide ∝ Δ  with  ≈ 2 depending on the accuracy of the fluid evolution scheme (Ripperda et al. 2022).The physical interpretation of the resistivity  is that it acts as a proxy for kinetic effects within the plasma.
We investigate plasmoid formation from fast relativistic plasmoiddominated reconnection.Whether the plasma becomes plasmoid unstable is determined by the Lundquist number  =  ′   /, with typical length of the current sheet  ′ and the Alfvén velocity   (see section 2.3 for definition).In order to trigger the fast reconnection and tearing-or plasmoid-unstable regime, the Lundquist number needs to satisfy  >  crit where  crit ∼ 10 4 (Loureiro et al. 2007;Bhattacharjee et al. 2009;Uzdensky et al. 2010).Note that the Lundquist number is largely determined by the underlying resistivity ( = 5 • 10 −5 ) which is set as a constant and uniform quantity in our resistive simulations.Then, if we estimate probable values of  ′ ≈ 1 and   ≈  = 1, we find  = 2 × 10 4 which lies above the threshold.At first glance, for the ideal simulations, one might think that as  ide is very small it reaches a sufficiently high Lundquist number.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.

Plasmoid identification
The starting point of our plasmoid identification routine lies in finding a quantity that lays bare the intrinsically circular magnetic field geometry.A natural choice for this identification quantity would then fall to the magnetic flux function; == where √ − is the metric determinant.Note that √ −  corresponds to the magnetic field in the Eulerian frame and that the magnetic flux function Ψ B corresponds to   , except for a minus sign discrepancy (Sironi et al. 2016).The magnetic flux function Ψ B is a good choice as its isocontours will follow the inplane magnetic field lines (i.e., B • ∇Ψ B = 0).More specifically, as plasmoids are characterized by their circular magnetic field configuration, the plasmoid center will correspond to the local maxima or minima in Ψ B ("O-points").Due to our methodology, the base Ψ B structure is not the ideal starting-point of the pipeline.We, therefore, work with the following quantity; where ΨB scalar denotes the (image) averaged flux function at a given time.The removal of the averaged flux function allows for clearer identification of plasmoids in Ψ B .
As we now have a suitable medium from which we can start to identifying the plasmoids, we will need a method that is able to classify the magnetic island structure reliably.For this purpose, we have developed an algorithm that consists of four steps: (i) All simulations contain a lot of fine-structure in the magnetic flux function.This makes it hard to differentiate between (magnetic) turbulence and more global features that correspond to a presence of a plasmoid.Therefore, to make certain we filter out much of the turbulence, we apply a blurring (Gaussian or flat) kernel to the flux function ( Ψ B ).This also gives us control over the size of features we want 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 from highly turbulent primary needs the most blurring relatively, while the GRMHD simulation are well-served with a fairly light blurring method.
(ii) Following the blurring step, we identify the local minima or maxima that will correspond to the plasmoid's center.
(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) maxima, so that we have control over what is being fitted while simultaneously improving the quality of the fit.
(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 center along x for the Harris sheet and inner and outer radii along 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 will be summarized below.(i) As the Harris sheet setups have periodic boundaries, one needs to be careful to catch plasmoids that are on the boundary.(ii) Additionally, capturing both "big" and "small" plasmoids in the Harris sheet setup required two different approaches, mainly concerning the blurring kernel.For the big features, one has to apply to a relatively small kernel many times (several hundred times works well in our experience) to not flatten out the global structure too much.To capture the smallest features, 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 following the ULS criterion.(iii) For GRMHD, one has to take into account that resolution is concentrated near the black hole and in the equatorial plane and therefore has non-uniform cell-sizes.(iv) Due to this nonuniform grid layout (for the GRMHD simulations), 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, then we interpolate the data to an uniform grid structure.
Lastly, we would like to note more explicitly how plasmoids are identified using the magnetic flux function in other works.In essence, one identifies plasmoids via the so-called "O"-and "X"-points.Opoints 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;Servidio et al. 2010;Zhdankin et al. 2013;Kadowaki et al. 2018;Zhou et al. 2020) via; Then, one calculates the matrix determinant of the Hessian (| Ψ B |) to find the critical points that correspond to | Ψ B (x)| = 0 at a given coordinate x.The eigenvalues of the Hessian then determined if we have an O-point 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 picked up 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.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 employed 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).

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 et al. 2001;Birn & Hesse 2001, also in Goedbloed et al. 2010).We start with a (wide) rectangular box with periodic boundary conditions on all sides and initialize two sheets of matter on top of an uniform background density that is scaled with  0 ; where   ,   ,  bg , and  are the box half-size in x and ŷ, the background factor, and the layer half-thickness, respectively.The initial values that were used for these parameters (and others) are denoted in Table 1.
We assume an uniform resistivity;  = 5 • 10 −5 , and an initialized magnetic and electric field according to Here,   denotes a (1%) 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) for PIC simulations.Note that we do not apply the typical guiding magnetic field perturbation that guides the initial plasmoids to the edges and creates a well-controlled reconnection region in the middle of the simulation domain (as perscribed for the GEM challenge, also in Keppens et al. 2013).To acquire pressure equilibrium at initialization we define the fluid pressure to be Additionally, we define the length and time scale as a function of system length  = 2   , so that (, ) ∈ [−0.5, 0.5] × [−0.125, 0.125] for Hb and (, ) ∈ [−0.5, 0.5] × [−0.25, 0.25] for Hs with a typical time unit of   = / (see Table 1).For completeness, we note that the computational length unit is  = 1 with corresponding time-scale / = 1, which both reduce to unity due to the geometrical unit assumption ( =  = 1).
If one were interested in relating the initial layer half-thickness  (see Table 1) to the resistivity , then one finds that / = 1000 (/ = 2000) for Hb (Hs).Nevertheless, we will connect it to a more intrinsic plasma-physical timescale in our unit set in the following paragraph.This is typically determined by the upstream Alfvén velocity   , which is defined as where ℎ = 1+ γ /( γ−1)  is the specific enthalpy with adiabatic index γ = 13/9 and  = √  2 = √︁     denotes the magnetic field strength.Additionally, the ("hot") magnetization is defined as  =  2 /ℎ.While we will primarily use the light-crossing time, it is worthwhile to connect it to the Alfvén and (resistive) diffusion timescales of the system, which then become   ≈  ′ /  and   ≈  ′2 / with  ′ being the current sheet's length (Ripperda et al. 2019b).Figure 2 gives an overview of the evolution of the Harris sheet (for the Hb case).From the magnetization () panels, we find that  ∼ 5 near the sheet, which indicates an upstream Alfvén velocity   ∼ .Then, one can determine the Lundquist number via  =   /  , but it becomes clear   is very large and   ∼  ′ which indicates that  will be similarly large.
Lastly, we would like to note that all boundaries are fully periodic (similar to Keppens et al. 2013;Takamoto 2013;Cerutti et al. 2013Cerutti et al. , 2014, and some quasi-periodic works in Sironi & Spitkovsky 2014;Petropoulou & Sironi 2018).This implies that no matter is lost so that evolution eventually saturates after having formed several 'monster' plasmoids that effectively act as a reservoir spanning a large part of 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 plasmoids become too large the sheets are influenced by one another.Another approach has 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).This tends to give less chaotic current sheets and allows for longer evolution times as, for periodic boundaries, the large plasmoids will eventually affect the opposing current sheet.The periodic Harris sheet simulations are primarily meant to have another verification case for the identification algorithm, but tend to display more complex behavior than what is found for the outflowing variety, especially combined with a global magnetic field perturbation (so that sign() •   ≳ 0; Loureiro et al. 2012;Sironi et al. 2016).Nevertheless, we did make sure that the magnetization was comparable to the GRMHD simulations.

GRMHD configuration
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 et al. 1973), we will only comment on the modification from the standard KS coordinates (, , , ), which is done via; Here,  and  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  and .All our GRMHD simulation use user-defined parameters ℎ = 0.25 and  0 = 0, which implies that the resolution of the underlying grid will be more concentrated in the equatorial plane.
Before continuing, we would like to outline a few specifics about the 3 + 1 split that is employed in BHAC.The line element is described as follows; with , ,  denoting the lapse, shift, and geometric part of the metric (  ), where Roman characters ,  ∈ {1, 2, 3} and Greek characters ,  ∈ {0, 1, 2, 3}.The metric determinant is then defined as √ − =  √ .Consistent with the conventions introduced in Porth et al. ( 2017), we denote electromagnetic quantities in the Eulerian frame with capitalized letters while lower-case 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   = {−, 0, 0, 0} (or contravariantly;   = {1/,   /}).
In this work, we will only consider Magnetically Arrested Disc (MAD, Igumenshchev et al. 2003;Narayan et al. 2003) models which are initialized via the vector potential The simulations are initialize with a torus that is in hydrodynamic equilibrium (Fishbone & Moncrief 1976), except for a perturbation to the fluid pressure , and is threaded by a single poloidal magnetic field loop (that is initialized via B = ∇ × A with A = (0, 0,   )).The inner and pressure maximum radii of the torus that determine the size and available matter are set to  in = 20  and  max = 41  for a black hole spin of  * = 0.9375.The scaling of the vector potential is set so that  = / mag = 100, with  being the gas pressure and  mag the magnetic pressure.Other user-defined parameters of the evaluated configurations can be found in Table 2.For completeness, we will note that the less magnetized accretion scenario is known as the Standard And Normal Evolution model (hereafter SANE, De Villiers et al. 2003;Narayan et al. 2012;Sadowski et al. 2013).

Energetics
An important objective in this work is to quantify if plasmoids 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 correspond- ing components of the stress-energy tensor   ; Here, the hereto unexplained quantities are , , and   , which are the specific interal energy, the fluid pressure, and the fully antisymmetric symbol, respectively. em denotes the electro-magnetic energy density (Qian et al. 2017),  kin the kinetic energy density, and  en the thermal energy density (McKinney et al. 2012;Ripperda et al. 2019a).The subscripts "EM", "PAKE" and "EN" correspond to the electro-magnetic, free particle, and enthalpy terms of the stressenergy tensor   (primarily following McKinney et al. 2012).The free thermokinetic energy (denoted as "MAKE" in McKinney et al. 2012) is the sum of  kin ("PAKE") and  th ("EN").This is important to note because  kin is predominantly negative in our GRMHD simulation, which can be interpreted from the geometric Bernoulli criterion (  ⩽ −1) corresponding to unbound matter.The term (  + 1) will therefore be negative (positive) when the fluid element is unbound (bound) and as   is positive we will end up with a negative  kin for bound matter that is typically found within the accretion disc.Lastly, note that the minus-sign in front of    is due to the metric signature (−, +, +, +) and is needed to get positive values.
Next, we define the covariant surface average (denoted by a bar, Q, over a given fluid variable) by with the surface , in an arbitrary coordinate system, denoted as The  corresponds to the geometric part of the metric as explained in section 2.4.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 quantities are calculated in the Eulerian (or laboratory) frame.

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 plasmoids that are known as "primary" plasmoids (see, e.g., Loureiro et al. 2007;Uzdensky & Loureiro 2016;Comisso et al. 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, they have a higher density 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.Overall, this results in a lower overall magnetization, but also a relatively lower magnetic field strength in relation to the surface.Their composition is primarily determined by the initial conditions.Following the initial break-up of the layer (at ∼1.56   for Hb and ∼5.27   for Hs), a continuous and steady creation of "secondary" plasmoids has started in the reconnection layer between the primary islands that remains active till the very end of the simulation window.These plasmoids do probe the underlying plasma characteristics and are relatively unaffected by the initial conditions.Two animations are attached to Fig. 2 which show both a window correponding to the figure and the entire simulation domain over time.
Following the formalism by Uzdensky, Loureiro & Schekochihin (2010) (hereafter ULS) that when a plasmoid coalesces with a larger plasmoid, then the smaller one is considered to be part of the larger body, and is therefore no longer taken into account.In practice, however, the small plasmoid will retain it's structure for some time depending on the size (ranging several 0.05   ) before conforming to 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  = 0.135  and is roughly 0.02  in width initially) takes approximately O (0.1   ) from the moment of impact to being fully absorbed by the primary plasmoid.When two plasmoids of similar size coalesce, then this timescale tend to be even longer and significant perturbation is needed before one of the two loses it's structure.
Generally, it is not simple to enforce the ULS criterion, which is reflected by the two-step approach outlined in section 2.2.Starting with secondary plasmoids, the minimum size for which we identify this population is set to ∼10 −4  (0.005 ), 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 plasmoids are identified with a very high fidelity.The primary plasmoids 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 plasmoids keep colliding and merging into it.These continuous perturbations also give rise to some magnetic reconnection within the primary plasmoid structure.As described in section 2.2, we need to apply an aggressive blurring kernel to identify the global primary plasmoid structure, but we still want to pick up on the distinct plasmoid structure if they have not fully merged.This implies that two plasmoids that have a similar magnetic flux signature (an example is seen at  = 0.34) 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 ULS criterion.At the interface of these two plasmoids one often finds new plasmoids forming.Naturally, all previously mentioned points become less pronounced at lower resolutions as one is resolving the current sheets less well which results in less formed plasmoids and less fine-structure.
The end of the evaluated window (at 4.1   for Hb and 8.79 Hs) is determined by the amount of interference the current sheets have on one another.Beyond these times, the few primary plasmoids 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 occurring at interfaces between the primary plasmoids that now have lost their elliptical shape and have become more hexagonal in shape.The simulation has a closer resembles to a turbulent box simulation than to the initial double current sheet configuration.This is beyond the scope of this work and therefore we chose the evaluated time windows to correspond to a clear current sheet structure.

Plasmoid statistics
Figure 3 displays two-dimensional 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 ( ρ) panels, one finds a main triangular distribution that spans −1.25 < log 10 ρ < −0.25.In addition to the main distribution, there is a secondary channel corresponding to −0.25 < log 10 ρ < 0.25 that corresponds to the densest plasmoids which also seem to occur over the entire width range.These plasmoids are, to summarize, in part due to minor misclassifications and due to the simulation conditions quickly after break-up of the initial layer.For the former, we find plasmoids that often correspond to fluctuation in the magnetic flux function within a large primary plasmoid.This population corresponds with a small plasmoid half-width.For the latter scenario, there are a number of high density reservoirs of matter that will eventually contract into the primary plasmoid population and will generally correspond to a large plasmoid half-width.Returning to the "true" plasmoid population, spanned by −1.25 < log 10 ρ < −0.25, we find that the smallest detected plasmoids have a half-width  ≈ 2 • 10 −4  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 Δ ≈ 0.02 ) for the evaluated data.For the surface-averaged magnetization ( σ), we find that the main population spans −2.5 < log 10 σ < 0.5.As is also seen in Fig. 2, the secondary plasmoids have a remarkably similar  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 σ rises with half-width, up to  ≈ 6.3 • 10 −3 .For log 10 / > −2, the σ mean plateaus and even seems to decrease slightly for the largest plasmoids.After the growth phase (in log 10 / ≲ −2), it seems that the increase in density and magnetic field strength is roughly matched.Lastly, for β, we find a similar but inverse trends to what we described for σ.The part of the distribution with the largest plasmoids ( ∼ 0.1 ) seems to deviate significantly from the main population and possesses a relatively high β ≳ 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  values that in turn affects the surface-averaged quantity ( β).
For the energies ( εem , εkin , and εth ), we find that the thermal energy ( th ) is the leading term in the total energy budget of the plasmoids with a mean (denoted by the green line) that remains fairly constant (0.0 < log 10 εth < 0.25) as a function of half-width ().At smallest , it appears the second term is the electro-magnetic energy (at εem ≈ 10 −1.5 ) that steadily becomes more significant for increasing width.As the kinetic energy ( kin ) is closely tied to the velocity of the plasmoid, we find that it can actually become a competing term for the electro-magnetic energy, especially in the active reconnection regions and merging (or colliding) plasmoids (see Fig. 2).The distribution of εth and εkin are wide indicating significant variance, while εem closely follows the distribution of σ and seems to show a more consistent trend.This trend is explained by secondary plasmoids becoming more magnetized with time, until they grow up to a size of  ∼ 0.01, after which they generally encounter a primary plasmoid and are absorbed after which the growth in magnetization ( εem ) stagnates.The high variance in εkin is explained by the fact that acceleration of plasmoids only happens in very localized regions -predominantly in active reconnection regions and just before plasmoids coalesce.As soon as the secondary are absorbed by the primary plasmoids, εth will be the leading term by a significant factor.Even though εth is still most dominant in the secondary plasmoid, both εkin and, especially, ε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 talked about the Hb, in the left-most panels of Fig. 3. Nevertheless, we find that all findings 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 () that is twice as wide and resolution that is lower by a factor two.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, simulation box length (in x, long side) is halved and it contains more matter due to the thicker initial layers when compared to the Hb case.As these are only minor differences, we find that the evolution is similar, which is also reflected by the results here, except that the primary plasmoids seem to span a greater part 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, we find there is reasonable variation between the probability density function over time, but a consistent image emerges as well.Generally speaking, at the smallest plasmoid half-widths (up to / ≈ 10 −3 ), we find a plateau followed by a steady decrease in occurrence frequency as the plasmoids become larger, up to the largest plasmoids that span a tenth of the simulation domain ( ∼ 0.1).Mainly via the slope  = −d log  /d log(/), we will be able to quantity the growth rate of plasmoids in the system.

Plasmoid distribution functions
The scaling laws have been studies in detail in the past (Uzdensky et al. 2010;Loureiro et al. 2012;Huang & Bhattacharjee 2012;Sironi et al. 2016).The density function of plasmoid width was predicted and verified to scale according to  () ∼   1) for both the Hb (left panels) and Hs (right 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.
perturbation (or outflowing boundaries), relative velocities between plasmoids are stochastically determined and relatively low, so we expect a greater similarity with Huang & Bhattacharjee (2012).Overall, we find that Ψ B and  do not scale with the same , which is contradictory with 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.
For half-width, we find  = 1.81 ± 0.05 for the Hb case and  = 1.48±0.06for the Hs case.Overall, for , we find that that we recover a scaling that is close to  () ∼  −2 corresponding to  = 2.For the mean trend in magnetic flux (in dark grey), we find  = 0.64 ± 0.10 for the Hb case and  = 0.59 ± 0.06 for the Hs case.However, the trend described by the smallest values per bin (in light grey) is  ≈ 1, which indicates agreement with Huang & Bhattacharjee (2012).The evolution of the distributions is characterized by a relative overrepresentation of large plasmoids, with |Ψ B / 0 | ∈ [5•10 −2 , 10 −1 ], that expands itself both to the left (lower |Ψ B |, smaller plasmoids) and right (higher |Ψ B |, larger plasmoids) over time.The smallest plasmoids have the lowest magnetic fluxes (as is also verified in Fig. 2) and the largest plasmoid will increase in |Ψ B | over time.This evolution also creates the sizable one-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  = 1 over time (especially for a low |Ψ B |), but this relation is affected by a high |Ψ 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.
It is important to note that our simulation setup differs substantially on at least two fronts from the previously mentioned scaling law studies, namely that it is relatively unperturbed and that it has no outflowing boundaries.With unperturbed, we mean that there is no guiding magnetic field perturbation present, which recreates 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.3).In practice, this implies that; (i) coalescence of plasmoids is a relatively more 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%) error on the values of .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.

General evolution
Figures 6 and 7 display the typical structure of the (two-dimensional) MAD (Tchekhovskoy et al. 2011;McKinney et al. 2012) simulations.After having evolved sufficiently, they will saturate in magnetic flux that penetrates the event horizon (see section 3.2.4).Following such a saturation event, the accretion flow is completely halted in axisymmetric (2.5D) simulation, while in 3D a so-called "flux tube" forms (Dexter et al. 2020;Porth et al. 2021).There, instead of halting the accretion flow completely, a localized less dense, more magnetized cavity moves outward from the black hole.These outbursts occur semi-periodically and seems to be even more prevalent in the relatively more confining 2D simulations.
Another feature is that the Magneto-Rotational Instability (MRI, Balbus & Hawley 1991), responsible for angular momentum transport, is suppressed as the main magnetic field component component is strongly poloidal in MAD simulations (Porth et al. 2021).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) that causes incursions into the disc structure (Marshall et al. 2018, and references therein).The Kelvin-Helmholtz Instability (KHI) becomes important in regions with strong shear flows as are the conditions at the jetdisc interface and are characterized by swirl-like vortices (see, e.g., Begelman et al. 1984;Hillier 2019).All these instabilities are perturbatative channels that are able to set off magnetic reconnection in the accretion disc.Therefore, 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 plasmoids' ability to produce flares, which are known to originate close to the central black hole, we apply our algorithm only within the inner 25  g .In Figs. 6 and 7, we find the plasma quantities and energies (similar to Fig. 2) for the iM5 and rM5 cases.The magenta and green colors denote found plasmoid corresponding to local maximum and local minimum in the magnetic flux function, respectively.Both figures show typical phases of MAD evolution that happen in all the GRMHD simulations in this work.The panels ( − ℎ) of Fig. 6 correspond to a flux eruption where we find the accretion flow is entirely halted.The panels ( − ℎ) 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 plasmoids are the collisional (non-pair-production plasma) equivalent to what has been seen in GRPIC simulation of diffuse collisionless magnetospheres around BHs (Crinquand et al. 2021;Bransgrove et al. 2021).
The overall structure and location of the plasmoid chains indicate that at the disc-jet boundary one finds plasmoids that correspond to local maxima (magenta) while when plasmoids occur within the disc they correspond to local minima (green).The magenta contours seem to have a lower density () and higher magnetization () than the ones in the disc.They also seem to be smaller when compared to the green contours.Their location and smaller size indicates that they are likely created by the shear-induced KHI.The purple contours also tend to leave the identification domain ( ⩽ 25  g ) on short timescales (5-10  g /) as they rapidly move outwards with turbulent jet-disc layer (also referred to as jet sheath).The green contours are tied to the bulk motion of the disc's fluid giving them more time and matter to interact with which explains their larger size.The energy and plasma parameter distributions will be explained in more detail in the next section (3.2.2).However, before we continue, we would like to point out that the values (visible in the , | kin |, and  th maps) near the vertical axis ( = 0  g ) are due to floor violations, which happen sufficiently far from our areas of interest and will therefore not interfere with the analysis.

Plasmoid statistics
Figure 8 shows two-dimensional histograms with various plasmoid quantities as a function of width for the two GRMHD simulations (iM5 and rM5).Before we comment on the general findings from the histograms, we would like to point out that we find a significantly lower plasmoid count for the iM5 case when compared to the the rM5 case.The ideal distributions are therefore more sparsely sampled.We will address this point in more detail in section 3.2.4.Overall, however, we do find that the distributions of iM5 and rM5 are consistent with one another, except for the aforementioned difference in occurrence rate.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 this does not inherently take into account the spacetime curvature.Therefore, we have chosen to display the distribution as a function of "circular" radius  S = √︁ / as the surface  calculation is taking into account the curvature.As plasmoids are generally elliptical, we loose information about the plasmoid shape as the ratio between width and length is no longer defined.
Starting with the distribution of ρ (panels . & .), we find that the surface-averaged density is highest for the smallest plasmoids at log 10 ρ ≈ 0.75 and then plateaus at log 10 ρ ≈ −0.5 from −0.5 < log 10  S < 1.0.For σ (. & .), we find a roughly constant mean value of log 10 σ ≈ −1, but a wide spread in values is also present.For β (. & .), one finds a very elongated distributions centered around a mean of roughly log 10 β ≈ 0.5 which has a complicated origin.This behavior is largely explained by the 'green' (local minima in Ψ B ) and 'purple' (local maxima in Ψ B ) plasmoid populations.For the purple contours, we find the origin of the elongated β distribution as the plasmoids detected in the jet sheath correspond to a distribution centered on a relatively low log 10 β ≈ −1.The typical distribution of β is fairly uniformly distributed with −2 ≲ log 10 β ≲ 2 centered on the mean value of log 10 β ≈ 0 − 0.5.Similar arguments apply to the distributions for ρ and σ where they both near identical means but a larger variance is present of the purple contours.For the green contours, we find more uniform and compact distributions overall that are located around the means for the entire (both green and purple) distribution as shown in Fig. 8.
For the energies ( em , | kin |, and  th ), there are only minor difference between the green and purple distributions, so we will just discuss the combined distributions for the energies in panels (. -.ℎ) and (. -.ℎ).Interestingly, the mean for all energy distributions describe an almost identical path -starting at log 10 ε ≈ 0 to ending at log 10 ε ≈ −2 for increasing  S .After a rapid decline up to log 10  S ≈ −0.5, we find that the log 10 ε means plateau, especially for εkin and εth .Additionally, the distributions indicate that the various surface-averaged energies are of similar strength.Nevertheless, εem does stand out with respect to the other energies as it has 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 εkin , we would like to note that εkin is negative, except in the jet-sheath where  kin ∼ O (1).This is explained in detail in section 2.5.Here, we will look at the absolute value | εkin | (. & .).The dashed purple and green lines in these panels correspond to the means of the distributions containing only the positive or negative values of εkin , respectively.So, it becomes clear that the vast majority of plasmoids has a negative εkin values as the global mean (in solid green) lies close to the dashed green line.Lastly, for εth , we find similar behavior as for the other energies combined with a relatively more considerable contribution at the lowest plasmoid sizes.
We define the magnetic Bernoulli factor as   = −(ℎ + /2)  , which incorporates the contribution of the magnetic pressure (/2) and therefore deviates slightly from the standard relativistic Bernoulli  = −ℎ  (Rezzolla & Zanotti 2013).The Bernoulli criterion states that the fluid is unbound when   > 1.Note that we have taken the liberty to incorporate a minus sign within the Bernoulli factor.
Returning to the distributions in panels (. & .), we find the majority of surface-averaged plasmoids is unbound as they pass the criterion, but there is still a significant number that lies under and close to the critical value of B  = 1 and are therefore bound.The mean of the function does however indicate B  ≈ 1 with a small number going up to relatively high values of B  ≈ 2. In the panels next to B  , we find the distributions of Ψ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 (Δ = |3000 − 4000|  g /).At the beginning ( = 3000  g /), we find a mean of ΨB ≈ 6.5, while at the end ( = 4000  g /) we find a mean of ΨB ≈ 11.5.
The last unexplained panels of Fig. 8 are two variations on the orbital velocity Ω =   /  .First, in panel (.& .), we investigate the ratio between the surface-average within the plasmoid contour ( Ωin ) with the surface-average for a shell directly outside the plasmoid contour ( Ωshell ).The outer edge of the shell corresponds to one-anda-half times the distance to the central O-point.From this quantity we can gauge if the plasmoid moves with its surroundings ( Ωin / Ωshell = 1) or disconnected from it ( Ωin / Ωshell ≠ 1).From the distributions, we find that the mean is consistent with Ωin / Ωshell = 1, but their While the other parameters have been outlined before, the magnetic Bernoulli factor is defined as B  = − (ℎ + /2)  and the orbital velocity Ω =   /  with the surface-averaged quantity inside the plasmoid being denoted as Ωin .We stack the distributions as a function of time from 3000  g / to 4000  g / with a 1  g / cadence and divide by the maximum.The green line denotes the mean per width bin that has more than 20 counts total.The left panels denotes iM5 case, while the right panels denotes the rM5 case.
is also a significant variance indicating that the plasmoid can move twice as fast or slow with respect to its direct environment.This can potentially be interpreted in a number of ways, which includes, e.g., that the plasmoid seems to be dynamically disconnected from the accretion disc in its direct surroundings (potentially driven by a plasma-instability).
Second, in panels (. & .), we evaluate the ratio of Ωin divided by the Keplerian circular orbital velocity in the equatorial plane which is defined as Ω K = ( 3/2 +  * ) −1 with  the cylindrical radius (corresponding the horizontal axis in Figs. 6 and 7) and  * = 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 Ωin /Ω K ≈ 0.8.Nevertheless, the broad distribution with 0.1 ≲ Ωin /Ω K ≲ 1.3 indicates the potential for plasmoids 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 plasmoids have non-zero   or   velocities that break both the circular and equatorial assumption for Ω 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.The quantitative differences in plasmoid identication rate ( P ) will, however, be outlined explicitly for all cases in section 3.2.4.

Plasmoid distribution functions
Figure 9 displays the probability density function (  ) of plasmoid radius  S , while Fig. 10 displays the probability density function of plasmoid half-width .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 3.1.3)and reflects the plasmoids shown in Figs. 6 and 7. Interestingly, we recover the power law indices of  = 1.88 ± 0.06 ( = 1.90 ± 0.05) and  = 2.15 ± 0.11 ( = 2.09 ± 0.09) for iM5 and rM5 in Fig. 9 (10), respectively.These are similar to the results described in section 3.1.3,which gives an indication that plasmoid formation is driven by the same principles even while taking into account the curvature of the spacetime.So, even with the additional perturbations by the plasma instabilities outlined in section 3.2.1,we still find scaling that is consistent with  ≈ 2. 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.
In a way, the identification strategy we employ for the GRMHD simulations is more consistent with the aforementioned works that have outflowing boundary conditions as we stop identifying plasmoids when  > 25  g .Another "outflowing" boundary lies at the horizon but the vast majority of plasmoids moves outwards (in r) 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 yields that plasmoid radius goes all the way up to  S ≈ 10 g .The Cartesian projection equivalent in Fig. 10 yields a radius of  ≈ 2  g .These largest plasmoids are visible in Fig. 6.The smallest detected plasmoid radii correspond to  S ≈ 10 −1  g (and  ≈ 10 −2  g ).Especially the largest plasmoids seem to be comparable in size to the 'hot spots' that were used to interpret flares around Sgr A * (Gravity Collaboration et al. 2020b;Wielgus et al. 2022b;Vos et al. 2022).From our simulations, we find that the plasmoids are of sufficient size to give a physical origin to these hot spots.However, currently, we do not explicitly interpret their emission potential, but as plasmoids are typically hot (/ ≳ 1) and magnetized (⟨ σ⟩ ≳ 0.1, as per Fig. 8) they are likely to create a 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 rM5, plasmoids with radii  S > 2.5  g occur at least once and three times on average for all evaluated time instances (corresponding to 8.2% of all identified plasmoids), while plasmoids with (Cartesian-projected) widths  > 1  g are much less common as they occur in only half (51.4%) of the evaluated snapshots (corresponding to 1.8% of all identified plasmoids).This perceived discrepancy is partially due to the space-time curvature (not taken into account for ) and the mixing of plasmoid length and width for the  S quantity.For iM5, the occurrence rates of at least one plasmoid passing the  S and  criteria are 57.7% and 16.5% (with 6.6% and 2.3% for all identified plasmoids over the entire time window), respectively.Overall, if we take into account the much lower plasmoid counts for iM5, we find that the percentage between the two cases are comparable, except for having at least one  > 1 plasmoid per evaluated time.This is well-explained, however, in section 3.2.4.
Lastly, we note that the power law gradient  = −d log  /d log( S / 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 colors 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 3.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 power law.Even though close to  ≈ 2, there is a minor break visible and the gradient becomes shallower at  S / g ≈ 4. As especially the larger plasmoid size bins contain more counts, this naturally pushes  to slightly lower values for iM5.Nevertheless, it is interesting that rM5 indicates a somewhat steeper gradient with  = 2.15 ± 0.11.However, combined with the points raised at the end of section 3.1.3,we conclude that iM5 and rM5 are consistent with a power law with  ≈ 2 as more robust claims can not be made without further investigation.

Timeseries of plasmoids and horizon penetrating fluxes
Plasmoids form within the accretion disc and are mostly governed by their local plasma conditions.Nevertheless, it would be interesting to see how these accretion disc probes connect dynamically to the central black hole.Quantities that are typically calculated to assess this are the mass accretion rate ( ) and surface-penetrating magnetic flux (Φ B ) that are defined as (in Porth et al. 2019); iM5 rM5  4 is also applicable here.Note that the quantities here do not correctly take into account the spacetime curvature, which is the case for Fig. 9.
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 two-dimensional simulations, the accretion flow will be stopped completely due to the constraining nature of the setup.The parameter that is used to quan-tify this behavior is the so-called MAD parameter  BH = Φ B / √ , which corresponds to the normalized magnetic flux.The MAD parameter saturates (in 3D) at  BH ≈ 151 (cf.Yuan & Narayan 2014).In our simulations, as shown in Fig. 11, we find that  occasionally rises to  BH ∼ 120.This is due to the confining nature of the 2D simulation, which allows for a greater accumulation of magnetic flux before an eruption.It is consistent with behaviour found for simulations in Ripperda et al. (2020).As we used a different adiabatic index γ = 13/9 (vs.γ = 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 (- ) of Fig. 11 display the number of identified plasmoids  P per simulation.While not shown explicitly in the figure, we confirm that plasmoids for either polarity (i.e., purple and green contours in Figs. 6 and 7) are equally abundant.As we already indicated (in section 3.2.2), a significantly lower number of plasmoids is detected for the ideal simulations than for the resistive ones where a factor 2 − 10 difference (in  P ) is common.The mechanism that triggers plasmoid formation, via the tearing instability, is not well-defined in ideal simulation and, more specifically, resolution dependent (∝ Δ 2 , with Δ being the cell-size).This implies that numerical resistivity ( ide =  num ) is lower close to the black hole then further away due to the MKS coordinate system and is significantly smaller than  ris =  ( num << ).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  ≈ 3700  g /, 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 plasmoids are created is somewhat increased (except for iM5).
The resistive simulations possess a surprisingly constant number of plasmoids ( 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  Δ ≡  Δ / Δ (see EHTC et al. 2022b, and description therein), with  Δ and  Δ being the standard deviation and mean, respectively, over a time-interval Δ = |3000 − 4000|  g /.We calculated the modulation index for both the accretion rate  and magnetic flux Φ B that penetrate the spherical shell at  = 2.5  g .The modulation indices for our simulation are listed in Table 3.There seems to be little difference between   or  Ψ B for the ideal and resistive simulations.This is surprising as  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 ⟨ P ⟩ ∼ 75.While for the higher resolution cases rM4 and rM5, we find ⟨ P ⟩ ∼ 100.As we find converging plasmoid numbers for both resolution cases, we conclude that the current sheet width (set by  = 5 × 10 −5 ), within the 25  g domain, is fully resolved starting from a resolution of 4096 2 .
In the last panel (), we cross-correlate the plasmoid number ( P ) with the negative gradient of the magnetic shell-penetrating flux  (−∇Ψ B ) and find a positive relation for most cases.Except for rM3, which is the uncorrelated component on the background (in lightest purple), we find a clear correlation between the most prominent peak in  P and a decrease in Ψ 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 rM3 results as there is no clear decrease in Φ B present.For iM5 at  ≈ 3780  g /, the flux eruption is rather large as is indicated by the decrease in Φ B , which has pushed the maximum corr(−∇Φ B ,  P ) further to the right.Just before the flux eruption, we find that an increase (of several factors) in Φ B after which it will start to drop.The positive correlation is naturally explained by the fact that the flux eruption, that 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 rM4 case displays different behavior with a drop in  P directly after the flux eruption.This is in part explained by our identification strategy which only identifies plasmoids within 25  g and as the disc has receded during the flux eruption the domain in which plasmoids can form also shrinks and effectively delays the peak in  P .Additionally, the shell-penetrating magnetic flux (Φ B ) only shows a relatively minor depression which indicates a relatively minor flux eruption and subsequent perturbation of the disc structure.So, in short, one can expect a reaction in the plasmoid formation rate following a flux eruption, which tends to increase the plasmoid count as it perturb the disc triggering reconnection.

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 three-dimensional (3D) results, (iv) effects of resistivity, and (v) a discussion of the flaring potential of plasmoids.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× more plasmoids.Both aforementioned works utilize the Bernoulli factor ( = −ℎ  ) 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 attempt with this proposed method but we did not reach the desired efficacy or fidelity, which started the development of algorithm outlined in this work.Overall, we typically find 5 − 10× more plasmoids than the previously mentioned works, which are not all attributed to the detection method difference.Other po-tential 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.
Starting with the identification medium, which we took to be the magnetic flux function Ψ B as it naturally identifies places with circular magnetic field structure.When we compare this with using the Bernoulli factor , then it is 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 advantageous to using , because one can apply wellestablished image-recognition algorithms if one is able to increase the contrast (i.e., only show a limited color-range) to which the  lends itself well.Nevertheless, this comes with the cost that one can only identify a subset of the plasmoid population.
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) utilizes an initial single dipolar loop up to intricate multi-polar initial magnetic field configurations with an evolution that can be described as SANE-like (with low  BH ∼ 2).Especially, the multi-polar 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 multi-polar 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 simulation were using ideal MHD, so we only compare it to the iM3, iM4, and iM5 cases.The differences between resistive and ideal GRMHD will be discussed in detail in section 5.
3D.How applicable are 2D results to a 3D reality?A number of arguments come into play here.First, 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 φ), one significantly complicates the situation.First, the plasmoid morphology would change and gain the resemblance of a flux rope.Second, the chance for interaction would decrease significantly as it is simply less likely to come across another flux rope.Third, the definition of flux ropes coalesce is difficult as they likely merge in a single place but not in it's entirety.These points are clearly demonstrated for the 3D equivalent of the Harris sheet as presented in, e.g., Sironi & Spitkovsky (2014); Cerutti et al. (2014).There, one finds complex behavior 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) 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.
Resistivity.In essence, setting a resistivity () 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 clearly outlined in 3.2.4,there is a clear discrepancy between the resistive and ideal simulations.While the former has a relatively consistent plasmoids number  P ∼ 100, the latter has a non-flaring count comparable to  P ∼ 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 the light of recent finding by the Even 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 this is further worsened by the use of resistive MHD.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. 2022).A rigorous (GRMHD) study including several resisitivities 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.
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 plasmoids that are in close vicinity to each other are grouped as they have very similar Ψ B signatures.Except that this diminished the detected plasmoid count somewhat, 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 lengthscales while not picking up on erroneous features.This is largely determined by the blurring layer, which dictates the minimal sizescale to which one is sensitive and gives a handle on how much finestructure one want to include.As the large plasmoid 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 to run, even despite parallelization attempts which should be intensified in the future.At present, we do not give an exact number of misclassifications, but one is able to find a few in most snapshots while the vast majority (of O (100)) is classified correctly.The number of plasmoids that were not classified is also of 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.
Flaring potential.While we started this paper by talking about plasmoids as a potential connection to flares, it is nevertheless difficult to make direct emission interpretations.The main reason for this is that the emission properties of plasmoids in the BH accretion environment are still very unknown, especially as one would expect a significant non-thermal contribution.The utilization of a thermal synchrotron proxy (as 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 access the flaring potential of plasmoids including a non-thermal electron population or reconnection-dedicated description (Rowan et al. 2017).This is an interesting avenue to pursue in the future, as it is possible to pin-point the plasmoid's location with the algorithm.

CONCLUSIONS
We have been able to identify plasmoids in highly turbulent accretion disc surrounding SMBHs with a higher fidelity than has been achieved before, which allows for creating complete time-series and distributions with sufficient counts to assess the statistics.Additionally, we have also verified our methodology with a set of previously well-investigated Harris current sheet simulations and found they are consistent with finding from previous studies (Uzdensky et al. 2010;Loureiro et al. 2012;Huang & Bhattacharjee 2012;Sironi et al. 2016).Interestingly, the scaling laws (outlined in sections 3.1.3and 3.2.3)for both the Harris sheet and the GRMHD simulation are very similar, which indicates that plasmoid formation in the more complex accretion disc environment does not differ fundamentally from the Harris sheet picture.Using this 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 MHD that warrant further investigation with a more systemic study that includes other accretion scenario (e.g., SANEs).
The typical plasmoid in a MAD GRMHD simulation is equally dense and somewhat under-magnetized with respect to their surrounding, moves with its surroundings, and is likely to be unbound according to the Bernoulli criterion.Nevertheless, this behavior describes the averages of distributions and does not describe the deviations which occur frequently.Especially for the orbital velocities and boundedness of the plasmoids, one finds large spreads in the distributions.This indicates that plasmoids can both occur as superor 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 a violent event that (partially) pushes 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 plasmoids is, therefore, still a strong candidate for explaining both Keplerian (Gravity Collaboration et al. 2018b, 2020a) and super-Keplerian (Matsumoto et al. 2020) nearinfra-red observation of flares around Sgr A * .More specifically, we regularly recover plasmoid sizes that are comparable to the hot spots that were used to interpret flares at both NIR-and mm-wavelengths (Gravity Collaboration et al. 2020a,b;Wielgus et al. 2022a;Vos et al. 2022).Also, as we outlined in section 3.2.4,flux eruptions (corresponding with a decrease in horizon-penetrating magnetic flux Φ B ) and plasmoid formation are likely strongly correlated with one another, indicating that flux eruptions 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 wellcharacterised as a 'closed contour-detector around local extrema'.So, in the future, we are planning to apply our methodology to mapping 3D structure of plasmoids and/or flux tubes for accretion onto SMBHs.It would also lend itself well to other MHD or PIC identification applications, such as shearing or turbulent box simulations.

APPENDIX A: RESOLUTION CONVERGENCE FOR GRMHD SIMULATIONS
We have performed our simulations at three resolution levels, from lowest to highest; 2048 × 2048 (iM3 & rM3), 4096 × 4096 (iM4 & rM4), and 8192 × 8192 (iM5 & rM5).These correspond to the third, fourth, and fifth AMR level, which we will use for referance.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  g /, have been run at the fourth and fifth level.These simulation have been started from the third level snapshot at 2900  g /, 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  ∈ [3000, 4000]  g /.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 Ψ B , density , and magnetization  for all cases where we are especially interested in the difference with respect to the highest resolution case.Figures A1 and A2 display the magnetic flux function Ψ 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 ( > 0  g ) region, while the bottom jet-sheath shows 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  ≈ 45  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 ≳ 6%, which is quite small and reasonable when compared to the other quantities.This once again confirm that the flux function is a very suitable choice in identification medium as it is not very variable, which makes identification difficult.
That leaves the time-averaged density () results in Figs.A3 and  A4.Before we start the discussion, we would like to point out that the ARDs go up to 20%, which signifies that the differences are significantly 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 Ψ 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 in the past that the ideal case comes across as atypically quiet, which may have been the result of an unfortunate coincidence 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 simulation's variability.So, in the future, it would be interesting to undertake a more systematic study of resistive GRMHD simulations to see if the equatorial plasmoid is a common occurrence.Nevertheless, it is likely related to the confining nature of 2D simulations, as we commented before.

Figure 1 .
Figure1.A schematic decomposition of the plasmoid identification algorithm.In the top panels (a-d), we display a snapshot of a GRMHD simulation (rM3 at  = 3000  g /) 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  = 2.93   ).In the left column (panels a&e), one finds the base magnetic flux function Ψ B -the starting point.To apply the watershed (panels c&g), one needs to make sure that the plasmoid corresponds to a local minimum which is done with the quantity -Ψ B (panels b&g).The last column (panels d&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 circles.Other O-points in the displayed simulation domain are denoted by the open grey circles.

Figure 2 .
Figure 2. Representative state for the evolution of the Harris sheet for the Hb case corresponding to  = 2.46   .Rows (a) till (g) show the density , 'hot' magnetization  =  2 /ℎ, plasma  = /(  2 /2), electro-magnetic energy density  em , kinetic energy density  kin , thermal energy density  th , and magnetic flux function Ψ B .The purple contours denote plasmoid detections corresponding to local maxima in the flux function (Ψ 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.

Figures 4
Figures 4 and 5 display the probability density function (  ) of plasmoid half-width and the absolute surface averaged magnetic flux function (| ΨB |), respectively.A distribution is calculated at each time starting at the beginning of the evalutated window at  = 1.76   for Hb ( = 5.47   for Hs) in dark blue up to  = 4.1   ( = 8.79   ) in bright yellow.Starting with Fig.4, we find there is reasonable variation between the probability density function over time, but a consistent image emerges as well.Generally speaking, at the smallest plasmoid half-widths (up to / ≈ 10 −3 ), we find a plateau followed by a steady decrease in occurrence frequency as the plasmoids become larger, up to the largest plasmoids that span a tenth of the simulation domain ( ∼ 0.1).Mainly via the slope  = −d log  /d log(/), we will be able to quantity the growth rate of plasmoids in the system.The scaling laws have been studies in detail in the past(Uzdensky et al. 2010;Loureiro et al. 2012;Huang & Bhattacharjee 2012;Sironi et al. 2016).The density function of plasmoid width was predicted and verified to scale according to  () ∼  −2(Uzdensky et al. 2010;Loureiro et al. 2012), while for magnetic flux both  (Ψ B ) ∼ Ψ −2

Figure 3 .
Figures 4 and 5 display the probability density function (  ) of plasmoid half-width and the absolute surface averaged magnetic flux function (| ΨB |), respectively.A distribution is calculated at each time starting at the beginning of the evalutated window at  = 1.76   for Hb ( = 5.47   for Hs) in dark blue up to  = 4.1   ( = 8.79   ) in bright yellow.Starting with Fig.4, we find there is reasonable variation between the probability density function over time, but a consistent image emerges as well.Generally speaking, at the smallest plasmoid half-widths (up to / ≈ 10 −3 ), we find a plateau followed by a steady decrease in occurrence frequency as the plasmoids become larger, up to the largest plasmoids that span a tenth of the simulation domain ( ∼ 0.1).Mainly via the slope  = −d log  /d log(/), we will be able to quantity the growth rate of plasmoids in the system.The scaling laws have been studies in detail in the past(Uzdensky et al. 2010;Loureiro et al. 2012;Huang & Bhattacharjee 2012;Sironi et al. 2016).The density function of plasmoid width was predicted and verified to scale according to  () ∼  −2(Uzdensky et al. 2010;Loureiro et al. 2012), while for magnetic flux both  (Ψ B ) ∼ Ψ −2 B (following the same works) or  (Ψ B ) ∼ Ψ −1 B (Huang & Bhattacharjee 2012) were established.The main difference between scaling found byUzdensky et al. (2010) andHuang & Bhattacharjee (2012) lies in how they treat the relative velocity between plasmoids.WhileUzdensky et al. (2010) assumed it to be ∼  ,Huang & Bhattacharjee (2012) evaluates a size-dependent relative velocity (see alsoSironi et al. 2016).As our simulations have no guiding magnetic field

Figure 4 .
Figure 4. Probability density function  (/) of plasmoid half-width  (along x) that is scaled according to the total simulation box width ( = 2  ) for both the Harris sheet cases; Hb and Hs.From scaling arguments, proposed in (Uzdensky et al. 2010; Loureiro et al. 2012), it has been shown that the distribution is expected to scale ∝  −2 .The various probability density function profiles are colored according the time at which they occur over a range of  ∈ [1.76, 4.10]   for Hb and  ∈ [5.47, 8.79]   for Hs.The mean density profile and one-sigma error over time are denoted by the dashed dark grey line.The power law slope is determined via  = − d log  / d log(/).

Figure 5 .
Figure 5. Probability density function  ( | ΨB |/ 0 ) of plasmoid half-width | ΨB |/ 0  that is scaled according to the total simulation box width ( = 2  ) and initial magnetic field strength ( 0 ) for both the Harris sheet cases; Hb and Hs.The light grey dashed lines denote the mean of the lowest (1%) values per bin with a corresponding linear fit shown in the same color.For the rest, the description of Fig. 4 also applies here.

Figure 6 .
Figure 6.Overview of the iM5 simulation at  = 3840  g /.Here, we are in the middle of a flux eruption event with pushed back the accretion disc.The purple contours corresponds to local maxima and the green contours correspond to local maxima in the magnetic flux function (Ψ B ).In panels ()-(), one finds the density (), the magnetization (), the ideal to magnetic pressure ratio (), and the magnetic flux function (Ψ B ).In the panels ()-(ℎ), we find the electro-magnetic energy ( em ), kinetic energy ( kin ), thermal energy ( th ), and the magnetic Bernoulli factor (  = − (ℎ + /2)  ).The corresponding animations can be found in the following repository; https://doi.org/10.5281/zenodo.8318522.

Figure 7 .
Figure 7. Overview of the rM5 simulation at  = 3500   /.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.

Figure 8 .
Figure 8. Two dimensional distributions  (,  ) of the plasma quantities  ∈ { ρ, σ, β, ΨB , Ωin / Ωshell , εem , εkin , εth , B  , Ωin /Ω K } as a function of "circular" radius  S =√︁ /  of the plasmoid for both the iM5 and rM5 cases.While the other parameters have been outlined before, the magnetic Bernoulli factor is defined as B  = − (ℎ + /2)  and the orbital velocity Ω =   /  with the surface-averaged quantity inside the plasmoid being denoted as Ωin .We stack the distributions as a function of time from 3000  g / to 4000  g / with a 1  g / cadence and divide by the maximum.The green line denotes the mean per width bin that has more than 20 counts total.The left panels denotes iM5 case, while the right panels denotes the rM5 case.

Figure 9 .Figure 10 .
Figure 9. Probability density function  ( S /  ) of "circular" plasmoid radius  S for both the high-resolution cases iM5 (left) and rM5 (right).All identification takes place within a circle of radius  = 25  and we evaluate a time-window of  ∈ [3000, 3001, . . ., 3999, 4000]   /.The rest of the description for Fig. 4 is also applicable here, except now we utilize  S .

Figure A1 .
Figure A1.Time-averaged magnetic flux function Ψ B over time interval  ∈ [3000, 4000]  g / 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.

Figure A2 .
Figure A2.The description of Fig. A1 applies here as well, except here we show the resistive cases; rM5, rM4, and rM3.

Figure A3 .
Figure A3.The description of Fig. A1 applies here as well, except here we show the density .

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  0 , layer half-thickness , background factor  bg , and resolution (with corresponding AMR level).The total dimensions of the box are denoted by  ∈ [ −  ,   ],  ∈ [ −  ,   ].

Table 3 .
The modulation index  Q ≡  Q / Q with  Q and  Q denoting the standard deviation and mean of quantity Q ∈ { , Φ B }.This index gives a measure of the variability in the simulations' timeseries.