2.1. Plasmoid-dominated reconnection regime

The main focus of this paper is on NTPA. Efficient high-energy NTPA requires strong electric fields coherent over substantial distances, at least comparable to or larger than the Larmor radii of energetic particles and therefore much larger than the plasma kinetic microscales. On such scales the electric field should be the motional (${\boldsymbol u}\times \boldsymbol {B}$) ideal-MHD electric field (cf. e.g. Drake et al. Reference Drake, Arnold, Swisdak and Dahlin2019). This simple reasoning underlies the need to understand the origin and structure of bulk plasma motions. For the resulting electric fields to be strong, these motions need to be fast, e.g. Alfvénic. Often, including in the plasmoid-dominated magnetic reconnection regime, the relevant fast motions arise as a result of the nonlinear development of various instabilities. Thus, it is important to identify what instabilities operate in a reconnecting current layer (or a reconnecting plasmoid chain) and what motions they drive. We will argue that 2-D plasmoid-mediated reconnection layers are subject to two important instabilities: the secondary tearing (aka plasmoid) instability, leading to the growth of plasmoids, and the coalescence instability, causing the plasmoids to move towards or away from each other along the global layer. These two instabilities lead to bulk plasma motions along the layer, and the associated electric fields lead to two channels for NTPA, as we will explain below.

In this section, we outline the basic physical picture of plasmoid-dominated reconnection that we believe adequately captures the physics most relevant for NTPA. This underlying picture is qualitatively the same for any plasma-physical framework of 2-D reconnection: resistive or collisionless, relativistic or non-relativistic.

2.1.1. Two-dimensional plasmoid-dominated reconnection layer: general picture

For simplicity, we will consider only 2-D reconnection, thus ignoring 3-D effects. In addition to the simplicity considerations, this approximation is motivated by the observation that recent 3-D PIC studies (e.g. Werner & Uzdensky (Reference Werner and Uzdensky2017); see also Sironi & Spitkovsky Reference Sironi and Spitkovsky2014) indicate that 3-D collisionless magnetic reconnection produces non-thermal particle spectra that are quite similar to the spectra produced by its 2-D counterpart, at least for relativistic reconnection in pair plasmas (see, however, Dahlin et al. (Reference Dahlin, Drake and Swisdak2015, Reference Dahlin, Drake and Swisdak2017) for 2-D/3-D comparison studies of NTPA in non-relativistic reconnection with a finite guide field, and Werner & Uzdensky (Reference Werner and Uzdensky2021), for trans-relativistic reconnection), even though the layer's morphology may be quite different. In addition, most of the numerical PIC studies of magnetic reconnection have so far been done only in two dimensions and hence most of what we know about reconnection-driven relativistic NTPA is limited to the 2-D case; it therefore makes sense to focus on the 2-D case first, before tackling the general 3-D situation. We acknowledge, however, that the role of the system dimensionality remains an open issue and should be investigated further.Footnote ³

For convenience, we introduce the following system of coordinates (see figure 1): $x$ is the direction of the reversing reconnecting magnetic field $\boldsymbol {B}_{\boldsymbol {0}}$ (the outflow direction); $y$ is the direction perpendicular to the current layer (the inflow direction); and $z$ is the direction of the electric current and of the main reconnection electric field $E_{\rm rec}$ (sometimes called the out-of-plane, or ignorable direction). In general, there may also be a guide magnetic field $B_g$ in the $z$ direction. Together, $x$ and $y$ form what is often called the ‘reconnection plane’, and $x$ and $z$ form the reconnection-layer midplane (at $y=0$).

Figure 1. The coordinate system (blue) used in this article. The black lines show magnetic field lines in a plasmoid-dominated reconnection layer. The $x$-direction is the direction of the reversing reconnecting magnetic field, the $y$-direction is across the current layer (the direction of the field reversal) and the $z$-direction is that of the main electric field and current in the sheet (the ignorable direction).

We envision a vigorous reconnection process taking place in the large-system, plasmoid-dominated regime (e.g. Shibata & Tanuma Reference Shibata and Tanuma2001; Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Daughton et al. Reference Daughton, Roytershteyn, Albright, Karimabadi, Yin and Bowers2009; Huang & Bhattacharjee Reference Huang and Bhattacharjee2010; Uzdensky, Loureiro & Schekochihin Reference Uzdensky, Loureiro and Schekochihin2010; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012; Petropoulou, Giannios & Sironi Reference Petropoulou, Giannios and Sironi2016; Sironi et al. Reference Sironi, Giannios and Petropoulou2016; Werner & Uzdensky Reference Werner and Uzdensky2017; Petropoulou et al. Reference Petropoulou, Christie, Sironi and Giannios2018; Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018; Schoeffler et al. Reference Schoeffler, Grismayer, Uzdensky, Fonseca and Silva2019). The reconnection layer is broken up into a highly dynamic hierarchical plasmoid chain as a result of the secondary tearing, aka plasmoid, instability (Loureiro, Schekochihin & Cowley Reference Loureiro, Schekochihin and Cowley2007; Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Daughton et al. Reference Daughton, Roytershteyn, Albright, Karimabadi, Yin and Bowers2009; Samtaney et al. Reference Samtaney, Loureiro, Uzdensky, Schekochihin and Cowley2009; Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012; Loureiro, Schekochihin & Uzdensky Reference Loureiro, Schekochihin and Uzdensky2013); (see Loureiro & Uzdensky (Reference Loureiro and Uzdensky2016) for a recent review). The chain is essentially one-dimensional, consisting of a broad distribution of plasmoids (aka magnetic islands, centred around magnetic O-points) of different sizes, strung on a single line $y=0$ (which is a projection of the reconnection midplane) and connected to each other by reconnecting inter-plasmoid current layers containing magnetic X-points (see figure 1). The plasmoids grow (i.e. accumulate magnetic flux and mass) continuously via reconnection taking place in these inter-plasmoid layers. The $E_z$ electric field associated with these inter-plasmoid reconnection processes would exist even if the plasmoids were themselves stationary; this field is responsible for one of the channels of particle acceleration and thus will play an important role in our analysis, as discussed in § 2.2.

In reality, however, the plasmoids are not stationary – they move about in the reconnection midplane (i.e. in the $x$ direction) with different velocities (generally, of order the Alfvén speed $V_A$) in a complicated, chaotic fashion, a kind of 1-D plasmoid turbulence. Sometimes they merge with each other, overall maintaining a broad statistical distribution of sizes (see §§ 3.5–3.6 for a more detailed description). The flow dynamics that controls plasmoid motions in the $x$-direction is complex and is governed by the interplay of two factors: the coalescence instability, developing on a broad range of scales, and the large-scale inhomogeneities along the global layer, which drive divergent large-scale flows (Hubble flow, $u_x \sim x$) out of the layer. Overall, one may distinguish between two somewhat different situations, discussed next, which we may call (i) the nonlinear evolution of a tearing-unstable plasmoid chain and (ii) reconnection proper; this distinction is somewhat similar in spirit to that made in turbulence studies between, respectively, decaying and driven turbulence.

2.1.3. Self-similar electric and magnetic fields

In any case, we envision that the hierarchical plasmoid chain during the active reconnection stage has a self-similar, fractal structure (see figure 2), as was first proposed by Shibata & Tanuma (Reference Shibata and Tanuma2001) and then confirmed numerically in both the resistive-MHD regime (Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Huang & Bhattacharjee Reference Huang and Bhattacharjee2010; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012) and collisionless regime (Daughton et al. Reference Daughton, Roytershteyn, Albright, Karimabadi, Yin and Bowers2009, Reference Daughton, Roytershteyn, Karimabadi, Yin, Albright, Bergen and Bowers2011; Sironi et al. Reference Sironi, Giannios and Petropoulou2016; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016, Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018). This means that, when one looks closely at an inter-plasmoid current layer between two neighbouring plasmoids of similar size $w$ at some given level in the middle of the hierarchy, one again finds a plasmoid-dominated reconnection region with essentially the same, universal characteristic values of three key electromagnetic (EM) field components as found at all other levels in the hierarchy. These field components are: the reconnecting ($x$-direction) upstream magnetic field $B_0$ (and hence the corresponding Alfvén speed $V_A$, which sets the characteristic scale for the plasma motions along the layer), the inter-plasmoid reconnected (i.e. in the $y$ direction) magnetic field $B_1$ (e.g. averaged over one half of the inter-plasmoid layer under consideration) and the effective reconnection rate $E_{\rm rec}$ (i.e. the electric field in the $z$ direction). This self-similarity extends all the way from the global layer as a whole (of size $L$) at the top of the hierarchy down to the smallest elementary current layers (which are marginally stable to tearing and thus essentially laminar) at the very bottom of the hierarchy, of characteristic thickness $\delta$ and length $\ell$. In the case of collisionless reconnection without a strong guide field, $\delta \sim \bar {\rho }$ (the typical Larmor radius of particles in the layer in the upstream reconnecting magnetic field $B_0$) and $\ell \sim 10\text {--}30 \delta$. At these smallest scales, the self-similarity of the EM field structure breaks down; in particular, in the case of electron–ion plasma reconnection, the Hall effect becomes important at these scales, leading to the emergence of a quadrupole out-of-plane ($z$-direction) magnetic field and an in-plane bipolar electrostatic electric field (Sonnerup Reference Sonnerup1979; Terasawa Reference Terasawa1983; Shay et al. Reference Shay, Drake, Denton and Biskamp1998; Uzdensky & Kulsrud Reference Uzdensky and Kulsrud2006; Melzani et al. Reference Melzani, Walder, Folini, Winisdoerffer and Favre2014a; Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018). We will ignore these fields in the present analysis because they only affect the typical, average-energy particles, but not the highly energetic, non-thermal particles of interest to us here.

Figure 2. Characteristic Larmor radius of a given energetic particle (red circle) appears different relative to the sizes of the neighbouring plasmoids in a self-similar hierarchical (Shibata–Tanuma) plasmoid chain (black), when viewed at different levels in the hierarchy. It is smaller than its flanking plasmoid sizes at the top level, comparable to them at the middle level and bigger at the bottom level.

All three above-mentioned field quantities – $B_0$, $B_1$ and $E_{\rm rec}$ – will play important roles in our analysis. Since we are ultimately interested in the effects that these fields have on NTPA, it is important to recognize that what matters is how the structure of these fields appears to a given energetic (i.e. with energy $\gamma m_e c^2$ well above the average particle energy in the layer, $\bar {\gamma } m_e c^2$) particle under consideration. Such a particle will have a large Larmor radius $\rho (\gamma )$ and so its motion will be blind, nearly insensitive to small-scale EM structures of modest field strength. It will interact most strongly with EM structures that exist on scales of order its Larmor radius or larger. It is thus important to look at the EM fields as a function of scale $\lambda$ (in the $x$-direction). This philosophy is similar to the one used in analysing the self-similar dynamics in the inertial range of turbulence (on scales much smaller than the large driving scale but much larger than the dissipative scale) and, more specifically, turbulent NTPA, where one is interested primarily in the resonant interaction of a particle with turbulent eddies or waves of a given size.

While $B_0$ is relatively straightforward, the other two quantities ($B_1$ and $E_{\rm rec}$) are less trivial and deserve a careful discussion. We first discuss the nature of the reconnected magnetic field $B_1$ (in the $y$-direction) in the plasmoid-dominated reconnection regime. Let us consider it somewhere in the middle of the plasmoid hierarchy, i.e. on the scale of an inter-plasmoid reconnecting layer (or, more precisely, inter-plasmoid reconnecting chain) between two plasmoids of some intermediate size $w$ such that $\delta \ll w \ll L$. The length of this inter-plasmoid current layer, $\lambda$, is then much smaller than the global length of the layer $L$, but larger than the length $\ell$ of the shortest elementary inter-plasmoid current sheets; in simulations it is typically $\sim$10 times greater that the width $w$ of the two plasmoids flanking it. The reconnected magnetic field $B_1$ can be defined, e.g. in terms of the net reconnected ($B_y$) magnetic flux between the layer's X point and the edge of adjacent plasmoid, divided by $\lambda /2$. While there may be many smaller-size (i.e. belonging to the next levels of the plasmoid hierarchy) magnetic islands inhabiting the inter-plasmoid layer under consideration, they consist of closed magnetic flux surfaces and thus do not contribute to the $\lambda$-scale $B_1$. Instead, this field is build up from patches of ‘semi-open’ reconnected flux – the field lines that intersect the reconnection midplane only once in the region between the two big ($w$-scale) islands under consideration; these field lines may then envelope these plasmoids and close (i.e. intersect the reconnection midplane again) somewhere outside the inter-plasmoid layer in question (see Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010). The self-similarity of the plasmoid chain then dictates that the typical value of $B_1$ is independent of the scale within the hierarchy. It is thus comparable to the typical reconnected magnetic field values in the outflow exhaust regions of elementary current layers at the very bottom of the hierarchy, i.e. approximately $\epsilon B_0 \simeq 0.1\,B_0$ for collisionless reconnection (and approximately 0.01 $B_0$ for collisional, resistive MHD reconnection; but in this study for concreteness we will adopt the fiducial collisionless value of 0.1 $B_0$). This point will be important in our analysis when we consider the motion of energetic particles in such a layer.

Likewise, the self-similar reconnecting plasmoid chain is characterized by a universal (scale-independent) typical average value of the out-of-plane electric field $|E_z| = E_{\rm rec}$ in the inter-plasmoid current sheet at all levels of the hierarchy. This effective reconnection rate is thus equal to the characteristic microscopic reconnection rate at individual inter-plasmoid X-points in the elementary current layers (Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010). This electric field is associated with the nonlinear development of the secondary tearing instability. It can be conveniently parametrized in terms of the upstream asymptotic magnetic field $B_0$ as

(2.1)$$ E_{\rm rec} = v_{\rm rec} B_0/c = \beta_{\rm rec} B_0 = \epsilon B_0 V_A/c, $$

where the dimensionless coefficients $\beta _{\rm rec} \equiv v_{\rm rec}/c$ and $\epsilon = v_{\rm rec}/V_A$ represent the reconnection inflow speed $v_{\rm rec}$ normalized, respectively, to the speed of light $c$ and to the Alfvén speed $V_A$ defined with the upstream reconnecting field $B_0$ and discussed in more detail below. First-principles PIC simulations (e.g. Hesse et al. Reference Hesse, Schindler, Birn and Kuznetsova1999; Birn et al. Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001; Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018) and analytical theories (Comisso & Bhattacharjee Reference Comisso and Bhattacharjee2016; Cassak, Liu & Shay Reference Cassak, Liu and Shay2017; Liu et al. Reference Liu, Hesse, Guo, Daughton, Li, Cassak and Shay2017) indicate that the dimensionless reconnection rate $\epsilon$ for collisionless reconnection is typically $\epsilon \simeq 0.1$ in both relativistic and non-relativistic regimes, the value that we will adopt in this paper.Footnote ⁴ This is true even for electron–positron pair plasmas Bessho & Bhattacharjee (Reference Bessho and Bhattacharjee2005, Reference Bessho and Bhattacharjee2007), despite the absence of the Hall effect that has been linked to fast ($\epsilon \simeq 0.1$) collisionless reconnection in electron–ion plasmas in previous studies (e.g. Birn et al. Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001).

In a laminar and stationary 2-D reconnection layer, the electric field $E_z$ would be steady and uniform, as dictated by Faraday's law. In contrast, however, a realistic stochastic reconnecting plasmoid chain is highly dynamic and $E_z$ can vary dramatically in both space and time due to rapid (Alfvénic) chaotic plasmoid motions up and down the chain, e.g. driven by the plasmoid coalescence instability as discussed above. The electric field at any given point can then have strong fluctuations on top of the average value of 0.1 $V_A B_0/c$. These fluctuations may have the form of rapid, intense alternating-sign spikes of the motional electric field of amplitude as high as $E_{\rm pl} = V_A B_0/c$ (i.e. 10 times the mean), e.g. when a circularized plasmoid with $B_y \sim B_0$ passes by with $v_x \sim V_A$ (e.g. Zenitani & Hoshino Reference Zenitani and Hoshino2001, Reference Zenitani and Hoshino2005; Sironi et al. Reference Sironi, Giannios and Petropoulou2016; Philippov et al. Reference Philippov, Uzdensky, Spitkovsky and Cerutti2019). It is this strong motional electric field that is responsible for short intense particle acceleration events associated with an energetic particle bouncing off a rapidly moving plasmoid – a key element in what is described as stochastic Fermi acceleration.Footnote ⁵

As will be described in more detail in § 2.2, the model presented in this paper will incorporate both acceleration mechanisms, viewed as separate channels: acceleration by $E_{\rm rec}$, as a particle traverses an inter-plasmoid current layer (the motional component of this electric field is associated with the reconnection outflows out of the individual inter-plasmoid current sheets); and acceleration by $E_{\rm pl}$ as it bounces off a large rapidly moving plasmoid. An essential assumption that is important to us here is that the statistical properties of the electric field, e.g. the distribution of the electric field strengths and the average value, be the same at all levels of the hierarchy, i.e. for reconnecting sublayers of each scale.

2.1.4. Hot and cold magnetization parameters and relativistic reconnection

The relativistic Alfvén speed $V_A$ that appears in (2.1) is associated with the reconnecting magnetic field $B_0$ and is defined in terms of the ambient upstream plasma conditions. In particular, it is convenient to express it in a dimensionless form

(2.2)$$ V_A \equiv \beta_A c = c \sqrt{\frac{\sigma_h}{{1+\sigma_h}}}, $$

in terms of the so-called ‘hot’ upstream magnetization parameter $\sigma _h$, defined as the ratio of the enthalpy density of the reconnecting magnetic field $B_0$ to the relativistic (including rest-mass) enthalpy density $h$ of the upstream plasma (Melzani et al. Reference Melzani, Walder, Folini, Winisdoerffer and Favre2014a; Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018),

(2.3)$$ \sigma_h \equiv {\frac{B_0^2}{4 {\rm \pi}h}}. $$

For example, in the case of a pair plasma that is relativistically cold, i.e. has an upstream background temperature $T_b = \theta _e m_e c^2 \ll m_e c^2$, or in the case of a pure electron–ion plasma (with $n_{e,b} = n_{i,b}$) that is non-relativistic ($T_b \ll m_e c^2$) or semi-relativistic ($m_e c^2 \ll T_b \ll m_i c^2$), the enthalpy is dominated by the rest mass of the dominant particle species (electrons and positrons in the pair-plasma case and ions in the electron–ion plasma case). In these cases, the ‘hot’ upstream magnetization becomes the same as the ‘cold’ plasma magnetization for this dominant species

(2.4)$$ \sigma_{c} \equiv {\frac{B_0^2}{4 {\rm \pi}n_{b} m_e c^2}}, $$

for the pair case with $n_b = 2 n_{b,e}$ being the total (electron + positron) particle number density, and

(2.5)$$ \sigma_{c} \simeq \sigma_{ci} \equiv {\frac{B_0^2}{4 {\rm \pi}n_{ib} m_i c^2}}, $$

for the electron–ion case.

However, in the opposite case of an upstream plasma that is ultra-relativistically hot, i.e. $T_b = \theta _e m_e c^2 \gg m_e c^2$ for the pair-plasma case (as found in, e.g. PWN) or $T_b \gg m_i c^2$ for the electron–ion plasma case, the rest-mass contribution to the upstream enthalpy is negligible and so the enthalpy density becomes simply $h \simeq 4 P_b = 4 n_b T_b = 4 n_b \theta _e m_e c^2$ (where $n_b$ is the total background particle density). In this case, the ‘hot’ magnetization simply becomes (apart from a factor of 1/2) the inverse of the upstream plasma-$\beta$ parameter ($\beta _b \equiv 8{\rm \pi} P_b/B_0^2$),

(2.6)$$ \sigma_h(T_b \gg m_e c^2, m_i c^2) = {\frac{B_0^2}{16 {\rm \pi}n_b T_b}} = {\frac{1}{2\beta_b}}, $$

and thus differs dramatically from the ‘cold’ magnetization, e.g. $\sigma _h = \sigma _c/4 \theta _e \ll \sigma _c$ for the pair-plasma case and $\sigma _h = \sigma _{ci}/8 \theta _i \ll \sigma _{ci}$ for the electron–ion case with $n_{i,b}=n_{e,b}$.

Both $\sigma _h$ and $\sigma _c$ are useful quantities that will play important roles in our analysis. In particular, $\sigma _c$ reflects (up to a factor or order unity) the upstream magnetic energy per background particle and thus sets the basic characteristic energy scale for particles energized by the reconnection process (normalized by the particle rest mass). And the hot magnetization $\sigma _h$ controls how relativistic the upstream Alfvén velocity – and hence the bulk fluid motions in the layer – are. Namely, it allows us to distinguish two limiting cases:

1. (i) the relativistic reconnection regime: $\sigma _h \gg 1$ and hence $\beta _A \equiv V_A/c = [\sigma /(1+\sigma )]^{1/2} \rightarrow 1$ and $E_{\rm rec} \simeq \epsilon B_0 \simeq 0.1 B_0$;

2. (ii) the non-relativistic reconnection regime: $\sigma _h \ll 1$ and hence $\beta _A \simeq \sigma _h^{1/2} \ll 1$ and $E_{\rm rec} \simeq \epsilon \sigma _h^{1/2} B_0$.

For the sake of completeness, we mention here how some of the above key relationships are modified in the presence of a guide ($\hat {z}$-component) magnetic field $B_g$, even though this paper focusses on the zero-$B_g$ case. Following Werner & Uzdensky (Reference Werner and Uzdensky2017), in the case of a finite guide field, the plasma outflows that control the reconnection electric field correspond not to the full Alfvén velocity but to its projection onto the outflow ($x$) direction, i.e. to the in-plane Alfvén velocity (see also Liu et al. Reference Liu, Guo, Daughton, Li and Hesse2015)

(2.7)$$ {\frac{V_{A,x}^2}{c^2}} = \frac{B_0^2}{B_{\rm tot}^2} {\frac{B_{\rm tot}^2}{B_{\rm tot}^2 + 4 {\rm \pi}h}} = {\frac{B_0^2}{B_{\rm tot}^2 + 4 {\rm \pi}h}}, $$

where $B_{\rm tot}^2 = B_0^2 + B_g^2$. This expression can be rewritten as

(2.8)$$ {\frac{V_{A,x}^2}{c^2}} = {\frac{B_0^2}{B_0^2 + (B_g^2+ 4 {\rm \pi}h)}} = {\frac{B_0^2}{B_0^2 + 4 {\rm \pi}h_{\rm eff}}}, $$

where we have introduced the effective total enthalpy, $h_{\rm eff} \equiv h + B_g^2/4{\rm \pi}$, that includes the contribution of the guide magnetic field, $B_g^2/4{\rm \pi}$, in addition to the relativistic plasma enthalpy $h$. Physically, this additional contribution needs to be included because the guide field is advected together with the plasma out of the layer by the (Alfvénic) reconnection outflows, and this its inertia (and, more specifically, its enthalpy) has to be taken into account (Werner & Uzdensky Reference Werner and Uzdensky2017). With this in mind, we can write the relevant in-plane Alfvén velocity in the standard form

(2.9)$$ {\frac{V_{A,x}^2}{c^2}} = \frac{\sigma_{\rm eff}}{1+\sigma_{\rm eff}}, $$

where we have defined the effective hot magnetization (Werner & Uzdensky Reference Werner and Uzdensky2017)

(2.10)$$ \sigma_{\rm eff} \equiv {\frac{B_0^2}{4 {\rm \pi}h_{\rm eff}}} = {\frac{B_0^2}{(B_g^2+ 4 {\rm \pi}h)}}. $$

Thus, in the case of a relativistically strong guide field, $B_g^2 \gg 4 {\rm \pi}h$, when the plasma enthalpy is negligible and we are dealing with a relativistic force-free field, we get $\sigma _{\rm eff} \rightarrow B_0^2/B_g^2$, and so $V_{A,x}/c \rightarrow B_0/B_{\rm tot}$. And in the opposite, non-relativistic case, $B_g^2 \ll 4 {\rm \pi}h$, (note that this does not imply that the guide field is weak compared with the reconnecting field $B_0$), the effect of the guide field on $\sigma _{\rm eff}$ and $V_{A,x}$ can be ignored and we recover the standard expression $V_{A,x} = c[\sigma _h/(1+\sigma _h)]^{1/2}$.

2.2. Acceleration of relativistic particles in a reconnecting plasmoid chain

We shall now discuss the motion and acceleration of energetic relativistic particles in the reconnecting plasmoid chain. We will focus here on relativistic particles with energies $\gamma m_e c^2$ in the high-energy non-thermal tail of the distribution function, far above the average particle energy $\bar {\gamma } m_e c^2$. We will call such particles simply ‘energetic particles’ or ‘high-energy particles’. The Larmor radii of such particles (corresponding to the upstream magnetic field $B_0$), given by $\rho (\gamma ) = \gamma \rho _0$, where $\rho _0 \equiv m_e c^2/e B_0$, are much greater than the average electron Larmor radius, $\bar {\rho } \equiv \bar {\gamma } \rho _0$. Since the thickness $\delta$ of the smallest elementary current layers in collisionless reconnection is usually of order $\bar {\rho }$, this means that the energetic particles under consideration will have Larmor radii greater than $\delta$.

Because of this, the questions of whether the acceleration of such highly energetic particles is done by a non-ideal electric field or by an ideal (motional) electric field, or whether this accelerating electric field is mostly parallel or perpendicular to the local magnetic field, may not be very relevant or even well posed. A given highly energetic particle does not know or care whether the electric field accelerating it is ideal or non-ideal. What matters is just the electric field component along the particle's direction of motion, smoothed on the appropriate scale of the particle's motion. For a particle with $\gamma \gg \bar {\gamma }$, this scale is generally larger than the microscopic plasma scales (like $\delta$ or $\bar {\rho }$), on which one can determine whether the electric field is ideal or not. Thus, the electric field responsible for energetic particle acceleration may in general have both ideal and non-ideal components, but in our view this question is not particularly relevant for understanding NTPA of highly energetic particles.Footnote ⁶ Moreover, the distinction between ideal and non-ideal electric fields may be even less relevant in the kinetic picture; this concept requires identifying the plasma bulk velocity ${\boldsymbol u}$, i.e. the average particle velocity, but individual particles will have velocities that may in general be very different from ${\boldsymbol u}$. Furthermore, in relativistic plasmas even defining the plasma's bulk velocity is not a completely trivial task, as there are two different ways of doing this (the Landau & Lifshitz (Reference Landau and Lifshitz1959) and the Eckart (Reference Eckart1940) frames). In this case, it may be sensible to define the non-ideal parallel electric field in the de Hoffmann–Teller frame, which corresponds to the $\boldsymbol {E}\times \boldsymbol {B}$ drift and in which the perpendicular (to the magnetic field) electric field component vanishes.Footnote ⁷

Let us now discuss the main characteristics of an accelerating energetic particle's motion. The first two key features of this motion are: (i) confinement to the reconnection layer in the $y$ direction (i.e. across the layer) by the reversing upstream magnetic field, which always deflects the particle back towards the layer; and (ii) the particle's continuous and relatively steady acceleration along the layer in the $z$ direction by the main reconnection electric field $E_{\rm rec}$; along with the electric field $E_{\rm pl}$ involved in the Fermi acceleration by particle reflection off of moving large plasmoids (discussed below), it is this electric field that is ultimately responsible for primary particle acceleration in our model. Thus, the predominant motion of the particle in the $yz$ plane during its acceleration stage (i.e. until it gets magnetized by the reconnecting ($B_y$) field and eventually trapped by a large plasmoid, see below) can be described by the relativistic version of a Speiser orbit (Speiser Reference Speiser1965; Zenitani & Hoshino Reference Zenitani and Hoshino2001; Uzdensky et al. Reference Uzdensky, Cerutti and Begelman2011; Cerutti et al. Reference Cerutti, Uzdensky and Begelman2012a), where the particle wiggles in and out of the thin current layer into the region of stronger ($B_0$) upstream magnetic field, crossing the reconnection midplane $y=0$ multiple times as it continuously gains energy. Since the particle's Larmor radius in the upstream magnetic field is greater than $\delta$, as we discussed above, the particle may initially spend a substantial fraction of time in one of the two upstream regions on either side of the current sheet. However, as the particle moves along this trajectory and is accelerated by the reconnection electric field, over time its relativistic Speiser trajectory focusses closer and closer to the midplane, with the $y$-meandering getting progressively smaller (Kirk Reference Kirk2004; Contopoulos Reference Contopoulos2007; Uzdensky et al. Reference Uzdensky, Cerutti and Begelman2011; Cerutti et al. Reference Cerutti, Uzdensky and Begelman2012a), and hence the particle's confinement to the layer becoming tighter.

2.2.2. Magnetization of particles by the reconnected magnetic field

This process of regular acceleration by $E_{\rm rec}$ proceeds for a while, but then there is a certain probability per unit time, or unit path length, for the particle to encounter a patch of reconnected field $B_y$ that is strong and extended enough to magnetize it, i.e.

(2.11)$$ \rho(\gamma, B_y) = {\frac{\gamma m_e c^2}{eB}} =\gamma \rho_0 \frac{B_0}{B_y} \lesssim \Delta x, $$

where $\Delta x$ is the extent of the patch in the direction perpendicular to the reconnected field ${\boldsymbol B}_y$. A more general and precise (but basically equivalent) formulation of this condition is

(2.12)$$ \epsilon \equiv \gamma m_e c^2 < e \Delta \varPsi, $$

where $\Delta \varPsi \equiv \Delta x B_y$ is the magnetic flux (i.e. the $z$-component of the magnetic vector potential) of the patch. Interestingly, this condition is similar to that for electrostatic trapping in an electrostatic potential well of depth $\Delta \varphi$, with the electrostatic potential $\Delta \varphi$ replaced by the drop in the vector potential $|\Delta A_z| = \Delta \varPsi$. That is, particle trapping in plasmoids can be viewed as a magnetic analogue of particle trapping in electrostatic potential wells of nonlinear electrostatic waves (electron holes, etc.).

When this happens, the particle may with some probability be reflected back into the acceleration region and continue accelerating, or it may be taken out of the active acceleration process, at least temporarily. In a plasmoid-dominated reconnection layer, the reconnected field $B_y$ can be thought of as having a more-or-less bimodal distribution:

1. (i) the relatively weak, distributed reconnected field of typical strength $B_1 \sim \epsilon B_0 \sim 0.1 B_0$ present almost everywhere in inter-plasmoid current layers outside of circularized plasmoids (see § 2.1);

2. (ii) the $B_{\rm pl} \sim B_0$ field inside fully formed, circularized plasmoids with aspect ratios of order 1 (Sironi et al. Reference Sironi, Giannios and Petropoulou2016) (this field can be even stronger deep inside plasmoid cores).

Because of this bimodal nature of the $B_y$-distribution, in our analysis we will treat the magnetization by $B_1$ (see § 3.3) and particle interaction with large plasmoids as separate processes, with the latter further subdivided into particle trapping inside plasmoids (§ 3.4) and particle reflection by moving plasmoids (§ 3.7). Here, we present a qualitative physical discussion of these processes.

First, if a particle becomes magnetized by the $B_1$ field in an inter-plasmoid reconnection layer, it starts performing electric drift associated with $E_{\rm rec}$ and $B_1$, thus moving with the general bulk plasma outflow in the $\pm x$-direction, away from that layer's X-point. Its acceleration then slows down dramatically. However, since the reconnected magnetic field $B_y$ generally strengthens in the outflow direction, the particle's energy still increases gradually, responding adiabatically to the magnetic field compression while obeying the conservation of the particle's magnetic moment (the first adiabatic invariant), and also to the field-line shortening while preserving the second adiabatic invariant. Eventually, this stage ends when the particle reaches a big plasmoid at the end of the inter-plasmoid layer under consideration (see below); the particle then is either trapped by the plasmoid or kicked out back into the active acceleration zone where it can resume rapid acceleration. One may, therefore, examine whether it is possible for a particle that got caught in a patch of cross-layer reconnected magnetic field extended enough to magnetize it and thus inhibit its acceleration, to ever be ‘released back into the wild’ again. To investigate this question, let us suppose that the patch initially has a field $\sim B_1$, representing the typical field at the edge of a reconnection-layer outflow, and is subsequently pulled and absorbed into an adjacent big magnetic island (plasmoid). This phase corresponds to the contraction and circularization of this newly added part of the island, with the field strength in the patch increasing from $B_1$ to $B_0$ while preserving its flux $\Delta \varPsi$. Will the particle continue to be magnetized within this patch of reconnected flux? If this were indeed so, then during this contraction phase the particle's perpendicular energy would rise in betatron acceleration due to the conservation of the relativistic particle's magnetic moment, $\mu \sim \epsilon ^2 /B = {\rm const}$. Thus, the particle's energy would increase by a factor of $(B_0/B_1)^{1/2} \sim 3$. Therefore, if the magnetization condition (2.12) was initially satisfied only marginally, i.e. $\epsilon \simeq e \Delta \varPsi$, then it may no longer be satisfied by the end of this process; the particle may then escape from the patch under consideration, perhaps receiving a modest (order-unity) energy boost. It may then again be able to enter another acceleration region, e.g. another inter-plasmoid reconnection layer, and be accelerated by $E_{\rm rec}$ once again. If, however, the plasmoid is sufficiently large then even if the particle is no longer magnetized to the patch in question, it would still continue to be trapped by this plasmoid.

2.2.4. Particle reflection off large plasmoids and Fermi acceleration

The third possible outcome of an energetic particle's encounter with a large plasmoid is that it can bounce off the plasmoid and get back into the main acceleration region; during this bounce, however, the particle interacts with the motional electric field associated with the plasmoid's motion as a whole. This interaction is brief, of order the gyro-period of the particle in the $B_y$ magnetic field of the plasmoid, which for circularized plasmoids is of order $B_0$; however, this interaction can be quite intense because the electric field here is $\sim V_A B_0/c$. As a result, the particle may gain or lose a substantial amount of energy. Namely, assuming that an individual reflection is elastic in a moving plasmoid's frame, when considered in the laboratory frame, the particle can either gain or lose energy to the plasmoid depending on whether the plasmoid moves (in the $x$-direction) head-on or tail-on relative to the particle. From the microscopic point of view, this interaction can be understood as follows (see figure 3). Let us for definiteness consider a positively charged particle approaching a large plasmoid from the left, i.e. emerging from the inter-plasmoid layer located to the left of the plasmoid and moving to the right (positive $x$, say), i.e. in the direction of the general reconnection outflow. The plasmoid, however, can be moving either to the left, opposite to the particle's direction of motion (a head-on collision), or to the right, in the same direction as the particle (a tail-on collision). If the interaction is head-on (as in figure 3), the plasmoid moves in the direction opposite to the overall large-scale reconnection outflow, and, because $B_y$ in the left half of the plasmoid (facing the approaching particle) is in the same direction as in the adjacent inter-plasmoid current sheet, the motional electric field ($E_z$) in this part of the plasmoid is opposite to the main reconnection electric field $E_{\rm rec}$. Then, when the particle approaching from the left encounters this plasmoid, it gets partially magnetized by it for a short period of time: the plasmoid's $B_y$ field deflects the particle so that it performs an incomplete gyro-orbit (more than a half but less than full) inside this plasmoid, thereby reversing its $x$-motion, and then escapes from the plasmoid back into the inter-plasmoid current layer from which it came. Effectively, the particle is reflected by the plasmoid's magnetic field. Importantly, while the particle is covering this incomplete gyro-orbit inside the plasmoid, it is moving backwards in $z$ relative to its $z$-motion prior to the encounter. Thus, both the particle's $z$ direction of motion and the electric field are reversed relative to what they were in the main current-layer acceleration zone; therefore, they are again aligned and so that the particle gains energy from the plasmoid (see § 3.7 for quantitative details). This process is a special version of gyro-resonance acceleration over one gyro-orbit: the electric field and the particle's velocity ($z$ components) stay aligned during one gyro-orbit.

Figure 3. Trajectory of an energetic positively charged particle (shown in blue) initially moving to the right (positive-$x$) direction in the $xz$-plane (the current-layer midplane) as it is being reflected by a large plasmoid moving to the left (negative-$x$ direction). The vertical dashed line shows the plasmoid's left boundary. As the particle's $z$-component of motion is reversed by the strong reconnected magnetic field $B_y\sim B_0$ (out of the page) inside the plasmoid (to the right of the vertical dashed line), it maintains its alignment with the plasmoid's motional electric field, and thus continuously gains energy.

In the opposite case (a tail-on collision), when the plasmoid in moving to the right, the reconnected magnetic field on its left side (interacting with the particle) is still in the same direction as in a head-on collision, and so the particle gets temporarily magnetized for $\gtrsim$ half of a gyro-orbit and reflected by the plasmoid leftward, back into the layer, in a way similar to a head-on collision. However, the direction of the motional $E_z$ electric field in this part of the plasmoid is now the same as, not opposite to, that in the inter-plasmoid layer on the left. Hence, this electric field and the backward $z$-motion of the particle in its half gyro-orbit inside the plasmoid are now counter-aligned and thus the particle loses energy in the encounter.

The net effective result of many such encounters is Fermi (Reference Fermi1949) acceleration. If plasmoid motions are random, uncorrelated, then one gets second-order diffusive Fermi acceleration (considered in § 3.7). However, if a particle bounces repeatedly between two approaching large plasmoids, the energy kicks are strongly correlated and the particle rapidly gains energy in a first-order Fermi acceleration.

In summary, we would like to reiterate that here we try to avoid concepts like parallel or perpendicular electric field acceleration. The particles are primarily accelerated by the out-of-plane reconnection electric field, ${\boldsymbol E}_z$. In the absence of a guide magnetic field, as is the main focus of the paper, and of the quadrupolar out-of-plane Hall magnetic field that arises in collisionless electron–ion plasma reconnection, the accelerating electric field is naturally perpendicular to the magnetic field almost everywhere. But this fact may not be very important; the exact orientation of the electric field relative to the actual local, microscopic magnetic field (e.g. whether it is parallel or perpendicular) is not particularly relevant to high-energy particles. A given energetic particle is not sensitive to EM fields with scales much smaller than its Larmor radius. In particular, even though the electric field is perpendicular to the local magnetic field, this does not necessarily mean that its acceleration is due to a drift motion (such as, e.g. the curvature drift) along the perpendicular electric field, since, in order to describe a particle's motion as a drift, the particle needs to be magnetized in the first place.

3.1. Kinetic equation for NTPA in a plasmoid-dominated reconnection layer

We describe relativistic particle acceleration by reconnection quantitatively in terms of a kinetic equation governing the energy distribution function $f(\gamma )$ of the particles in the acceleration zone

(3.1)$$ \partial_t f(\gamma,t) = S(\gamma) - \partial_\gamma (\dot{\gamma}_{\rm acc} f) - \frac{f(\gamma)}{\tau(\gamma)} + \partial_\gamma (D_\gamma \partial_\gamma f). $$

Here, the first term on the right-hand side, $S(\gamma )$, is the source term describing the injection of the particles into the reconnection layer from the upstream region; $\dot {\gamma }_{\rm acc}$ in the next term is the rate of regular, continuous acceleration by the reconnection electric field $E_{\rm rec}$; the third term, $-f(\gamma )/\tau$, represents the effective escape of particles from the acceleration process by their capture onto cyclotron orbits around the reconnected field $B_y$; and the last term represents diffusive Fermi acceleration due to the particle bouncing off moving plasmoids.

We shall now make several simplifications. First, we will focus on the quasi-steady-state solution, which is appropriate for the subpopulation of energetic particles undergoing active, rapid acceleration before they are taken out of this process by being magnetized by the reconnected magnetic field and trapped inside plasmoids. In contrast, the subpopulation of energetic particles trapped in plasmoids grows continuously and is thus, of course, not stationary; it is being constantly fed by the escape (the $f/\tau$ term) of particles from the actively accelerating subpopulation under consideration here; however, we relegate the analysis of this trapped subpopulation to future studies (see, e.g. Hakobyan et al. Reference Hakobyan, Petropoulou, Spitkovsky and Sironi2021). As for the actively accelerating particles, it is reasonable to expect that they may develop a quasi-stationary spectrum governed by the balance between the injection $S(\gamma )$ at small energies, a flow up in energy due to the regular reconnection acceleration $\dot {\gamma }_{\rm acc}$ and Fermi acceleration by moving plasmoids (the second and last terms in (3.1)), and escape, $-f/\tau$.

Our second simplification is based on the assumption that the background upstream plasma is cold (magnetically dominated, $\beta _b \ll 1$). Then, the source term is concentrated at small energies characteristic of the upstream conditions, $\gamma _{\rm inj}$, much less than the energies $\bar {\gamma } \sim \sigma _c$ of accelerated particles of interest to us here. We will therefore ignore the injection term in our analysis of acceleration of high-energy non-thermal particles; this is similar in spirit to ignoring the details of energy injection at the large driving scale when analysing the inertial range of turbulence.

The steady-state kinetic equation (3.1) at $\gamma \gg \gamma _{\rm inj}$ then becomes

(3.2)$$ {-}\partial_\gamma (\dot{\gamma}_{\rm acc} f) - \frac{f(\gamma)}{\tau(\gamma)} + \partial_\gamma (D_\gamma \partial_\gamma f)= 0. $$

Next, as discussed in § 2, the reconnected magnetic field $B_y$ has a roughly bimodal distribution in a reconnecting plasmoid chain, corresponding to two types of regions: inter-plasmoid current layers with a relatively weak distributed reconnected field of order $B_1 \sim \epsilon B_0$, and circularized plasmoids where the field is of order $B_0$. To reflect this dichotomy, it is convenient to split the escape term in the kinetic equation into two separate terms encapsulating the two effective escape channels that check the continuous particle acceleration: magnetization in the general reconnected field $B_1$ and trapping in large plasmoids (see § 2)

(3.3)$$ \frac{1}{\tau(\gamma)} = \frac{1}{\tau_{\rm magn}(\gamma)} + \frac{1}{\tau_{\rm trap}(\gamma)}. $$

We will now discuss the individual terms in (3.2) one by one to better understand their underlying physics and the role they play in shaping up the particle distribution. For most of this discussion, however, for the sake of simplicity and analytical tractability, we will ignore the diffusive acceleration term. We will come back to it only in § 3.7, where we will evaluate its importance in different regimes and will obtain an explicit analytical solution in one special case. We will then discuss the general kinetic equation in § 3.8.

3.2. Regular acceleration by the reconnection electric field

Consider a particle that is moving in the layer before it gets magnetized by the reconnected field $B_1$ or trapped in a big plasmoid. The particle then undergoes a regular, steady acceleration by the main reconnection electric field $E_{\rm rec} = \epsilon \beta _A B_0$, with the acceleration rate given by

(3.4)$$ \dot{\gamma}_{\rm acc} = \tau_{\rm reg}^{{-}1} = {\frac{e E_{\rm rec} v_z}{m_e c^2}} \sim \beta_z \varOmega_0 \frac{E_{\rm rec}}{B_0} = \beta \varOmega_0 \epsilon \frac{V_A}{c}= \epsilon \beta_z \beta_A \varOmega_0, $$

where $\varOmega _0 \equiv e B_0/m_e c$ is the nominal non-relativistic electron cyclotron frequency and $\beta _z \equiv v_z/c$. As long as the particle is not yet strongly deflected in the $x$ direction by the Lorentz force due to the reconnected field $B_y$ (which is related to the condition that it is not magnetized by this field), it should have a finite, sizable velocity component in the $z$ direction, so that $\beta _z \simeq 1$ (assuming the particle is ultra-relativistic). Importantly, we then see that the acceleration rate (3.4) is independent of the particle energy,

(3.5)$$ \dot{\gamma}_{\rm acc} \simeq \epsilon \beta_A \varOmega_0 = {\rm const}. $$

Substituting this expression into the kinetic equation (3.1) without the diffusion term, we see that a steady-state distribution is governed by

(3.6)$$ \dot\gamma_{\rm acc} {\frac{{\rm d} f}{{\rm d} \gamma}} ={-} {\frac{f(\gamma)}{\tau(\gamma)}} \quad \Rightarrow \quad {\frac{{\rm d} \ln f}{{\rm d} \gamma}} ={-} \frac{1}{\dot\gamma_{\rm acc} \tau(\gamma)} \simeq{-} \frac{1}{\epsilon \beta_A \varOmega_0 \tau(\gamma)}. $$

Integrating (3.6), we obtain the stationary distribution function as

(3.7)$$ f(\gamma) = C \exp \left(-\frac{1}{\dot\gamma_{\rm acc}} \int^{\gamma} \frac{{\rm d}\gamma'}{\tau(\gamma')}\right) \simeq C\exp \left(-\frac{1}{\epsilon \beta_A \varOmega_0} \int^{\gamma} \frac{{\rm d}\gamma'}{\tau(\gamma')}\right). $$

Recalling (3.3), this distribution function can be represented as a product of two factors, one due to magnetization and one due to trapping

(3.8)$$ f(\gamma) = C\exp \left(-\frac{1}{\dot\gamma_{\rm acc}} \int^{\gamma} \frac{{\rm d}\gamma'}{\tau_{\rm magn}(\gamma')}\right) \times \exp \left(-\frac{1}{\dot\gamma_{\rm acc}} \int^{\gamma} \frac{{\rm d}\gamma'}{\tau_{\rm trap}(\gamma')}\right). $$

We will now proceed to discuss the two escape terms on the right-hand side of (3.3). As we will see, these two terms affect the mathematical shape of the resulting distribution function in different ways, because their characteristic time and length scales have different scalings with the particle energy $\gamma m_e c^2$. Thus, $\tau _{\rm magn}$ is basically directly proportional to $\gamma$, while $\tau _{\rm trap}$ is tied to the plasmoid distribution function. As we will show below, the first term controls the power-law slope and the second one controls the high-energy cutoff.

3.3. Particle magnetization by inter-plasmoid reconnected magnetic field $B_1$: governing the power-law index

As a particle moves through a reconnection region and experiences acceleration by the $E_{\rm rec}$ field, it also interacts with the reconnected magnetic field $B_y \sim B_1 \simeq \epsilon B_0$, which continuously deflects it more and more towards $\pm x$ direction, out of the accelerating layer. At some point, after the particle travels a certain distance $l_{\rm magn}$ in the $x$-direction, this deflection may become so large that the particle becomes effectively magnetized by this field. Its subsequent motion in the $xz$ plane then becomes dominated by the cyclotron motion associated with the field $B_1$, coupled with ${\boldsymbol E}\times \boldsymbol {B}$ drift in the $x$-direction, which corresponds to the particle moving with the general (fluid-level) reconnection outflow. The resulting inability of the particle to move unimpeded along the reconnection electric field (in the $z$ direction) means that its rapid energy gain becomes greatly diminished and its acceleration slows down.Footnote ¹¹ This phase continues until the particle encounters a large plasmoid, at which point it will have a finite chance to get absorbed into it (see § 3.4 below).

The distance $\ell _{\rm magn}(\gamma )$ that a particle of a given energy $\gamma m_e c^2$ has to travel in the $x$-direction before becoming effectively magnetized can be estimated simply as the Larmor radius corresponding to $B_1$ (e.g. Zenitani & Hoshino Reference Zenitani and Hoshino2001): $\ell _{\rm magn} (\gamma ) \sim \rho _L(\gamma, B_1) \sim (B_0/B_1) \rho _L(\gamma, B_0) \sim \epsilon ^{-1} \rho _L(\gamma, B_0)$ $= \epsilon ^{-1} \gamma \rho _0 \sim 10 \gamma \rho _0$. Since a particle in the acceleration region typically moves with a finite angle with respect to the $z$ axis, so that $v_x \sim v_z \sim c$, we can estimate the time for a particle to get magnetized as

(3.9)$$ \tau_{\rm magn}(\gamma) \sim \ell_{\rm magn} /c \sim \epsilon^{{-}1} \gamma \varOmega_0^{{-}1}. $$

Furthermore, the typical amount of energy that the particle gains by regular acceleration while it crosses the distance $\ell _{\rm magn}$ before becoming magnetized is of order $\delta \gamma (\ell _{\rm magn}) \sim e E_{\rm rec} \ell _{\rm magn}/m_e c^2 $ $= (E_{\rm rec}/B_0) \ell _{\rm magn}/\rho _0 \sim \gamma (E_{\rm rec}/\epsilon B_0) \sim \gamma (E_{\rm rec}/B_1) \sim \gamma \beta _A$. Thus, the fractional energy gain, $\delta \gamma (\ell _{\rm mag})/\gamma \sim \beta _A$, is of order unity in the case of relativistic reconnection ($\sigma _h > 1$ and hence $\beta _A\simeq 1$, see (2.2)), but becomes small, of order $\beta _A \sim \sigma _h^{1/2}$ in the non-relativistic reconnection case.

Importantly, the magnetization time (3.9) is directly proportional to the particle's energy, and this enables a power-law distribution to develop. Indeed the corresponding factor in (3.8) is

(3.10)$$ \exp \left(-\frac{1}{\dot\gamma_{\rm acc}}\int \frac{{\rm d}}{\tau_{\rm magn}(\gamma)}\right) = \exp \left(-\frac{1}{\epsilon \beta_A \varOmega_0} \int \frac{\epsilon \varOmega_0 {\rm d}\gamma}{\gamma}\right) = \exp \left(-\frac{1}{\beta_A}\int^{\gamma} \frac{{\rm d}\gamma'}{\gamma'}\right) \sim \gamma^{{-}p}, $$

where the power-law index is given by

(3.11)$$ p = p_{\rm magn}(\sigma_h) \sim \frac{1}{\beta_A} = \sqrt{\frac{1+\sigma_h}{\sigma_h}}. $$

Thus we see that the balance between regular acceleration by the reconnection electric field $E_{\rm rec} \sim \epsilon \beta _A B_0$ and particle magnetization by the reconnected magnetic field $B_1\sim \epsilon B_0$ produces NTPA with a power-law index $p_{\rm magn}$ that exhibits the same dependence on the hot magnetization $\sigma _h$ as was observed in recent numerical PIC studies (e.g. Guo et al. Reference Guo, Li, Daughton and Liu2014; Sironi & Spitkovsky Reference Sironi and Spitkovsky2014; Guo et al. Reference Guo, Liu, Daughton and Li2015; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016; Ball et al. Reference Ball, Sironi and Özel2018; Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018), namely:

1. (i) $p \rightarrow {\rm const}$ of order unity in the ultra-relativistic reconnection limit, $\sigma _h \gg 1$, $V_A \simeq c$;

2. (ii) $p \sim \sigma _h^{-1/2}$ in the non-relativistic case, $\sigma _h \ll 1$ and hence $V_A \simeq c\sigma _h^{-1/2} \ll c$.

We also note that this physical picture and the theoretical arguments are similar to those presented by Zenitani & Hoshino (Reference Zenitani and Hoshino2001) for a simple laminar (without secondary plasmoids) reconnection layer in the ultra-relativistic limit.

3.4. Particle trapping in plasmoids and the high-energy cutoff: general discussion

We shall now discuss the effects of the second factor in (3.8) – particle trapping in large plasmoids – and will argue that, in large systems, where reconnection proceeds in the plasmoid-mediated regime, this process controls the extent of the power-law segment of the particle energy distribution; in particular, it induces an exponential-like high-energy cutoff $\gamma _{c}$. Thus, in our model the particle spectrum at highest energies is shaped by the distribution of plasmoids. We will postpone the discussion of the other important aspect of particle–plasmoid interaction, i.e. the Fermi acceleration by particle reflections off rapidly moving plasmoids, until § 3.7.

As discussed in § 2.2, we use the term ‘large plasmoid’ as a technical term with a specific meaning: a plasmoid that is large enough to trap a particle of a given energy $\gamma m_e c^2$, i.e. a plasmoid with the size comparable to, or larger than, the particle's Larmor radius in the plasmoid's magnetic field. Also as discussed in § 2.2, the trapping condition is actually most accurately cast in terms of the plasmoid's magnetic flux $\psi$, i.e. the absolute value of the difference in the out-of-plane ($z$) component of the EM vector potential $A_z$ between the plasmoid's centre (i.e. the O-point) and its edge: $\epsilon = \gamma m_e c^2 < e \psi$. Thus, strictly speaking, we should be dealing with the plasmoid distribution function with respect to their fluxes. However, while the formulation of the theory in terms of plasmoid fluxes is more rigorous, its formulation in terms of sizes $w$ is arguably more intuitive and easier to visualize, and so this is the language we will use in this paper.

Once fully formed, plasmoids tend to be roughly circularized, with aspect ratios of order 1, and with characteristic magnetic fields comparable to $B_0$.Footnote ¹² Then, as discussed in § 2.2, we postulate a 1-to-1 correspondence between plasmoid fluxes and sizes, $\psi = w B_0$, and then the trapping condition can be written as

(3.12)$$ w \geq w_{\rm tr} (\gamma) \equiv \rho_L(\gamma, B_0) = \gamma \rho_0. $$

Correspondingly, there exists a characteristic distance that a given particle is able to travel in the $x$-direction before it is trapped by a large plasmoid: this is the characteristic separation $\lambda _{\rm pl}[w_{\rm tr}(\gamma )]$ between plasmoids of this size $w_{\rm tr}(\gamma )$. This separation, in turn, is controlled by the plasmoid-size distribution function in the reconnecting plasmoid chain, i.e.

(3.13)$$ \lambda_{\rm pl}(w) = \frac{L}{N(w)}. $$

Here, $L$ is the global length of the layer, and $N(w)$ is the number of plasmoids with size equal or greater than $w$, i.e. the cumulative plasmoid-size distribution function. It is related to the plasmoid distribution density $F(w)$ as

(3.14)$$ F(w) ={-} {\rm d} N/{\rm d} w, $$

i.e.

(3.15)$$ N(w) = \int_w F(w')\,{\rm d} w'+ {\rm const}. $$

Using (3.13), the characteristic trapping time $\tau _{\rm trap}(\gamma )$ that enters the kinetic equation (3.2) can be estimated as

(3.16)$$ \tau_{\rm trap}(\gamma) \sim \frac{\lambda_{\rm pl}[w_{\rm tr}(\gamma)]}{c} = \frac{L}{c} \frac{1}{N[w_{\rm tr}(\gamma)]}, $$

where we have assumed that a relativistic particle's motion through the acceleration region before it is trapped or magnetized is generally at a finite angle with respect to the $z$-axis, so that its typical velocity in the $x$-direction (i.e. along the layer) is relativistic, $v_x \sim c$. From (3.16) we see that higher-energy particles take longer to get trapped because they need to encounter a large enough plasmoid, which is rare.

The corresponding term in the kinetic equation is

(3.17)$$ {-}\frac{f(\gamma)}{\tau_{\rm trap}(\gamma)} \sim{-}f\frac{c}{\lambda_{\rm pl}[w_{\rm tr}(\gamma)]}\sim{-} f \frac{c}{L} N[w_{\rm tr}(\gamma)], $$

and the resulting suppression factor in the particle energy distribution function (3.8) is

(3.18)$$\begin{align} \exp \left(-\frac{1}{\dot\gamma_{\rm acc}} \int^{\gamma} \frac{{\rm d}\gamma'}{\tau_{\rm trap}(\gamma')}\right) & = \exp \left(-\frac{1}{\epsilon \beta_A \varOmega_0}\frac{c}{L} \int^{\gamma} N[w_{\rm tr}(\gamma')]\,{\rm d}\gamma' \right) \nonumber\\ & = \exp \left(-\frac{1}{\epsilon \beta_A}\frac{\rho_0}{L} \int^{\gamma} N[w_{\rm tr}(\gamma')]\,{\rm d}\gamma' \right). \end{align}$$

These expressions explicitly embody the relationship between particle acceleration and plasmoid distribution $N(w)$.

In order to evaluate the effect of this factor on particle acceleration in more concrete terms, we need a model for the plasmoid-size distribution function. For illustration, let us first consider a chain of equidistant identical plasmoids of a single size $w_*$ separated by a distance $\lambda _*$; thus, the number of such plasmoids in the chain is $N_* = L/\lambda _*$. The plasmoid distribution density then is $F(w) = N_* \delta (w-w_*)$ and the cumulative distribution function is given by the Heaviside step function: $N(w) = N_* [1- \varTheta (w-w_*)]$. These plasmoids can trap all particles with energies $\gamma \leq \gamma _* \equiv w_*/\rho _0$, and cannot trap particles with higher energies. Then the resulting suppression factor for particles with $\gamma \leq \gamma _*$ is

(3.19)$$\begin{align} \exp \left(-\frac{1}{\epsilon \beta_A} \frac{\rho_0}{L} \int^{\gamma} N[w_{\rm tr}(\gamma')]\,{\rm d}\gamma' \right) & = \exp \left(-\frac{1}{\epsilon \beta_A}\frac{N_* \rho_0}{L} \int^{\gamma} {\rm d}\gamma' \right) \nonumber\\ & \sim \exp \left(-\frac{1}{\epsilon \beta_A}\frac{\rho_0}{\lambda_*}\gamma \right) = \exp \left(-\frac{\gamma}{\gamma_c}\right), \end{align}$$

where the exponential cutoff $\gamma _c$ is given by

(3.20)$$ \gamma_c \equiv \epsilon \beta_A \frac{\lambda_*}{\rho_0} = \epsilon \beta_A \frac{\lambda_*}{w_*}\gamma_*. $$

This simple example illustrates the idea that (at least in two dimensions) large plasmoids can serve as particle traps, effectively taking the particles out of the rapid acceleration process, at least temporarily (see Dahlin et al. Reference Dahlin, Drake and Swisdak2015; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016; Dahlin et al. Reference Dahlin, Drake and Swisdak2017; Kagan et al. Reference Kagan, Nakar and Piran2018). The origin of the exponential cutoff can then be understood as follows. To get to a certain high energy, a particle has to spend proportionally longer time in the acceleration region, while all the while being subject to a certain, constant in this case, probability per unit time of being trapped and sequestered in a plasmoid, $\tau _{\rm trap}^{-1}(\gamma <\gamma _*) \sim c/\lambda _*$.

In the next two subsections we shall discuss two more realistic examples of $N(w)$, characterized by broad, power-law distributions.

3.5. Hierarchical plasmoid chain i: a single power-law model

Let us now consider the case when the plasmoid distribution function in the chain is given by a single power law,

(3.21)$$ F(w) \sim w^{-\alpha}. $$

We normally expect $\alpha \in [1,2]$ in realistic reconnecting plasmoid chains (Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010; Huang & Bhattacharjee Reference Huang and Bhattacharjee2012; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012; Sironi et al. Reference Sironi, Giannios and Petropoulou2016; Petropoulou et al. Reference Petropoulou, Christie, Sironi and Giannios2018), but here we shall keep the discussion general.

The cumulative distribution $N(w)$ is given by

(3.22)$$\begin{gather} N(w) \sim {\rm const} - \ln w \quad {\rm if}\ \alpha=1, \end{gather}$$

(3.23)$$\begin{gather}N(w) \sim w^{1-\alpha} \quad {\rm if}\ \alpha> 1. \end{gather}$$

Let us now further imagine that this power-law distribution extends up to some maximum plasmoid size $w_{\rm max}$, which we will tentatively associate with the size of the so-called ‘monster plasmoids’ (Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012), typically of order $\epsilon L \sim 0.1 L$ (Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012; Sironi et al. Reference Sironi, Giannios and Petropoulou2016; Petropoulou et al. Reference Petropoulou, Christie, Sironi and Giannios2018). Then, we can establish the normalization of the plasmoid distribution function by demanding that the number of these largest plasmoids of maximum size $w_{\rm max}$ found in the chain at any given time is approximately 1. In addition, we should account for the finite length $L$ of the layer, which imposes a strict limit on a particle's free accelerating motion set by the ejection time when a given particle traverses the entire layer, $\tau _{\rm ej} \sim c/L$. This can be effectively taken into account by requiring that $N$ cannot be less than 1, so that $\lambda _{\rm pl} \leq L$. Thus, we shall write

(3.24)$$\begin{gather} N(w \leq w_{\rm max}) \simeq 1 - \ln\left(\frac{w}{w_{\rm max}}\right),\quad \alpha = 1, \end{gather}$$

(3.25)$$\begin{gather}N(w \leq w_{\rm max}) \simeq \left(\frac{w}{w_{\rm max}}\right)^{1-\alpha},\quad \alpha > 1. \end{gather}$$

Now, let us consider the interaction of energetic relativistic particles with such a plasmoid chain, using the trapping condition $w \gtrsim w_{\rm tr}(\gamma ) = \rho _L(\gamma, B_0) = \gamma \rho _0$ as discussed above. It is convenient to define the maximum particle energy that can be confined by the largest plasmoids of size $w_{\rm max}$, i.e.

(3.26)$$ \gamma_{\rm max} \equiv w_{\rm max}/\rho_0. $$

We can then express the inter-plasmoid separation $\lambda _{\rm pl}(w_{\rm tr})$, and hence the typical path length that a particle can travel before it is captured, as follows. First, in the special case $\alpha =1$ (3.24) yields

(3.27)$$ \lambda_{\rm pl}^{\alpha=1}[w_{\rm tr}(\gamma)] = \frac{L}{N[w_{\rm tr}(\gamma)]} \simeq \frac{L}{1+ \ln [w_{\rm max}/w_{\rm tr}(\gamma)]} = \frac{L}{1+ \ln (\gamma_{\rm max}/\gamma)}, $$

which for $\gamma \ll \gamma _{\rm max}$ (and hence for $w_{\rm tr} \ll w_{\rm max}$) can be approximated as $\lambda _{\rm pl}^{\alpha =1}[w_{\rm tr}(\gamma \ll \gamma _{\rm max})] \simeq L /\ln (\gamma _{\rm max}/\gamma )$.

Next, for $\alpha > 1$, we have

(3.28)$$ \lambda_{\rm pl}[w_{\rm tr}(\gamma)] = \frac{L}{N[w_{\rm tr}(\gamma)]} \simeq L \left(\frac{w_{\rm tr}(\gamma)}{w_{\rm max}}\right)^{\alpha-1} \simeq L \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{\alpha-1}, \quad \alpha > 1. $$

In particular, in the important special case $\alpha =2$, we have $N(w) \sim w_{\rm max}/w$ and consequently

(3.29)$$ \lambda_{\rm pl}^{\alpha=2}[w_{\rm tr}(\gamma)] \simeq L \frac{w_{\rm tr}(\gamma)}{w_{\rm max}} = L \frac{\gamma}{\gamma_{\rm max}} = \rho_L(\gamma)\frac{L}{w_{\rm max}}. $$

For example, adopting $w_{\rm max} \simeq 0.1 L$, we can estimate this as

(3.30)$$ \lambda_{\rm pl}^{\alpha=2}[w_{\rm tr}(\gamma)] \simeq~10 w_{\rm tr}(\gamma) \simeq 10\rho_L(\gamma). $$

These expressions allow us to estimate the corresponding trapping rates (for $\gamma <\gamma _{\rm max}$)

(3.31)$$ \tau_{\rm trap}^{{-}1}(\gamma) \sim \frac{c}{L} \left[1 - \ln\left(\frac{\gamma}{\gamma_{\rm max}}\right)\right],\quad \alpha = 1, $$

approaching $(c/L) \ln (\gamma _{\rm max}/\gamma )$ for $\gamma \ll \gamma _{\rm max}$; and

(3.32)$$ \tau_{\rm trap}^{{-}1}(\gamma) \sim \frac{c}{L} N[w(\gamma)] \sim \frac{c}{L}\left(\frac{\gamma}{\gamma_{\rm max}}\right)^{1-\alpha},\quad \alpha > 1, $$

and, in particular,

(3.33)$$ \tau_{\rm trap}^{{-}1}(\gamma) \sim \frac{c}{L} \frac{\gamma_{\rm max}}{\gamma},\quad \alpha= 2. $$

Let us now investigate the implications of these estimates for our model of particle acceleration. Substituting them into (3.18), we find the corresponding expressions for the suppression factors (for $\gamma <\gamma _{\rm max}$):

(i) $\alpha =1$:

(3.34)$$ \exp \left( -\frac{\rho_0}{L\beta_A \epsilon} \gamma [2-\ln(\gamma/\gamma_{\rm max})] \right). $$

(ii) $1< \alpha < 2$:

(3.35)$$ \exp \left[-\frac{w_{\rm max}}{\epsilon \beta_A L(2-\alpha)} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{2-\alpha}\right] = \exp \left[-\left(\frac{\gamma}{\gamma_c}\right)^{2-\alpha}\right], $$

where

(3.36)$$ \gamma_c = \gamma_{\rm max} \left[ \beta_A \frac{\epsilon L}{w_{\rm max}} (2-\alpha) \right]^{1/{(2-\alpha)}}. $$

Interestingly, for $\alpha \rightarrow 1$ this factor approaches a simple exponential cutoff

(3.37)$$ \exp \left[-\frac{\rho_0}{L} \frac{1}{\beta_A \epsilon}\gamma \right] = \exp (-\gamma/\gamma_{c}), $$

where the cutoff Lorentz factor is

(3.38)$$ \gamma_{c} \equiv \beta_A \epsilon\frac{L}{\rho_0}, $$

corresponding to $\rho (\gamma _c) = \gamma _c \rho _0 \simeq \beta _A \epsilon L$ and formally independent of $w_{\rm max}$. This case corresponds to a situation where the plasmoid distribution function is so shallow that small plasmoids are not sufficiently numerous to have a strong effect on particle acceleration. The cutoff $\gamma _c$ then simply corresponds to the voltage drop due to the electric field $E_{\rm rec} = \epsilon \beta _A B_0$ over the entire layer's length $L$. This limit is thus equivalent to the natural high-energy cutoff expected in small laminar reconnection layers not succumbing to secondary tearing and plasmoid formation (Larrabee et al. Reference Larrabee, Lovelace and Romanova2003), as has been confirmed in numerical PIC simulations (e.g. Lyubarsky & Liverts Reference Lyubarsky and Liverts2008; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016). It is interesting to note that, if $w_{\rm max} \simeq \epsilon L \simeq 0.1L$, then we can express the cutoff as

(3.39)$$ \gamma_c \simeq \beta_A \frac{w_{\rm max}}{\rho_0} = \beta_A \gamma_{\rm max}. $$

Thus, in the case of ultra-relativistic reconnection, $\sigma _h \gg 1$ and $\beta _A \rightarrow 1$, this cutoff energy $\gamma _c$ is simply equal to $\gamma _{\rm max}$, the maximum energy of particles that can be confined in the largest plasmoid (consistent, e.g. with numerical observations by Sironi et al. Reference Sironi, Giannios and Petropoulou2016). However, in the case of non-relativistic reconnection, $\beta _A \sim \sigma _h^{1/2} \ll 1$, the cutoff may in principle be smaller than $\gamma _{\rm max}$. Thus, in this case there may be a substantial, measurable range of particle energies where the energy distribution is exponentially suppressed.

(iii) $\alpha =2$:

This case is indeed special. There is no exponential cutoff in this case; instead, particle trapping in plasmoids leads to an additional power-law factor, $\gamma ^{-p_{\rm tr}}$, similar to the effect of magnetization by $B_1$ (see § 3.3). The corresponding power-law index is given by

(3.40)$$ p_{\rm tr} \simeq \frac{1}{\epsilon \beta_A} \frac{w_{\rm max}}{L}. $$

This expression is similar to that for the power-law index $p_{\rm magn}$ due to magnetization, see (3.11). The combined effect of the two processes leads to a power law $\gamma ^{-p_{\rm tot}^{\alpha =2}}$ with an index given by the sum of the two

(3.41)$$ p_{\rm tot}^{\alpha=2} = p_{\rm magn} + p_{\rm tr} \sim \beta_A^{{-}1}\left(1 + \frac{w_{\rm max}}{\epsilon L}\right). $$

That is, trapping in plasmoids in the $\alpha =2$ case leads to a steepening of the non-thermal power law compared with expected from magnetization alone, while preserving (in contrast to the $\alpha < 2$ case) the overall power-law shape of the distribution function.

The relative importance of the magnetization and trapping processes for $\alpha =2$ does not depend on $\gamma$ and is instead controlled by the ratio $w_{\rm max}/\epsilon L$. In particular, if the chain is truncated at relatively small sizes, i.e. if $w_{\rm max} \ll \epsilon L$, then $p_{\rm tr} \ll p_{\rm magn}$ and the trapping correction to the total power-law index is relatively minor. If, however, the plasmoid chain extends all the way up to nearly the monster plasmoid scale, $w_{\rm max} \sim \epsilon L \sim 0.1 L$, then we get

(3.42)$$ p_{\rm tr} \simeq \beta_A^{{-}1} \sim p_{\rm magn}, $$

and hence the two processes generally play comparable roles for $\gamma < \gamma _{\rm max}$. The total power-law index is then greater by a factor of order unity (e.g. double) than that due to magnetization alone, i.e. the resulting steepening is significant. We will see, however, in § 3.7 that Fermi acceleration by moving plasmoids tends to compensate this steepening tendency to some degree and may even cancel it completely.

3.6. More realistic double power-law plasmoid chain

Numerical simulations of 2-D reconnection in the large-system, plasmoid-dominated regime, done both in the non-relativistic resistive-MHD case (Huang & Bhattacharjee Reference Huang and Bhattacharjee2012; Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012) and in the relativistic collisionless pair-plasma case (Sironi et al. Reference Sironi, Giannios and Petropoulou2016) (see also Petropoulou et al. Reference Petropoulou, Christie, Sironi and Giannios2018 for an associated semi-analytical model), show that realistic reconnecting systems develop plasmoid-size distributions that are more complex than a single power law discussed in the previous section. Namely, the resulting plasmoid distributions are best described by a double power law, i.e. a relatively shallow power law with some index $\alpha _1$ below a certain break plasmoid size $w_{\rm br}$ and another, steeper power law with an index $\alpha _2>\alpha _1$ above $w_{\rm br}$

(3.43)$$\begin{gather} F \sim w^{-\alpha_1},\quad w \leq w_{\rm br}, \end{gather}$$

(3.44)$$\begin{gather}F \sim w^{-\alpha_2},\quad w_{\rm br} < w \leq w_{\rm max}, \end{gather}$$

with a sharp, perhaps exponential, cutoff above $w_{\rm max}$ (see figure 4). The large-$w$ cutoff $w_{\rm max}$ of the second power law may correspond to the monster plasmoid size, typically $0.1 L$, as in § 3.5. Also, based on the above-mentioned simulation studies and analytical theory (Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010), we generally expect $\alpha _1 \sim 1$ and $\alpha _2 \sim 2$, but here will consider these indices as free variable parameters.

Figure 4. A more realistic, broken power-law plasmoid-size distribution with a shallow small-$w$ power law $F\propto w^{-\alpha _1}$ transitioning at $w=w_{\rm br}$ to a steeper large-$w$ power law $F\propto w^{-\alpha _2}$, with $\alpha _2>\alpha _1$.

The intermediate break scale $w_{\rm br}$ that marks the transition between the two power laws will play an important role in our analysis. Unfortunately, what governs it and how it depends on the system parameters is not yet understood; in particular, it is not known whether this size scales with (i.e. just some finite fraction of) the global size $L$, or with (is a fixed multiple of) the thickness of the smallest elementary inter-plasmoid layers $\delta \sim \bar {\rho } = \bar {\gamma } \rho _0 \sim \sigma _c \rho _0$ (where the latter estimate assumed that the upstream region is magnetically dominated, so that the average particle energy of the heated plasma in the layer is just a finite fraction of $\sigma _c m_e c^2$). Answering these important questions is beyond the scope of this paper and should be obtained by more careful numerical studies, e.g. similar to Sironi et al. (Reference Sironi, Giannios and Petropoulou2016), in conjunction with analytical and semi-analytical studies (e.g. Petropoulou et al. Reference Petropoulou, Christie, Sironi and Giannios2018). The analysis presented in this paper provides a strong motivation for such studies as it underscores the strong connection between the statistical properties of the plasmoid chain and NTPA. In the meantime, in this study we will consider $w_{\rm br}$ to be a variable parameter.

For $w \in (w_{\rm br}, w_{\rm max})$, the cumulative plasmoid distribution function is

(3.45)$$ N (w_{\rm br} < w \leq w_{\rm max}) \sim w^{1-\alpha_1}. $$

Using the normalization condition $N(w_{\rm max}) \simeq 1$, we can write

(3.46)$$ N (w_{\rm br} < w \leq w_{\rm max}) \simeq \left(\frac{w}{w_{\rm max}}\right)^{1-\alpha_2}, $$

thus also fixing the normalization for $F(w)$ in this range

(3.47)$$ F (w_{\rm br} < w \leq w_{\rm max}) ={-}N'(w_{\rm br} < w \leq w_{\rm max}) \simeq {\frac{\alpha_2-1}{w_{\rm max}}} \left(\frac{w}{w_{\rm max}}\right)^{-\alpha_2}. $$

We will also introduce

(3.48)$$\begin{gather} N_{\rm br} \equiv N (w=w_{\rm br}) \simeq \left(\frac{w_{\rm max}}{w_{\rm br}}\right)^{\alpha_2-1}, \end{gather}$$

(3.49)$$\begin{gather}\lambda_{\rm br} \equiv \lambda_{\rm pl}(w_{\rm br}) \equiv \frac{L}{N_{\rm br}} \simeq L \left(\frac{w_{\rm br}}{w_{\rm max}}\right)^{\alpha_2-1}, \end{gather}$$

and

(3.50)$$ \gamma_{\rm br} \equiv \gamma(w_{\rm br}) = \frac{w_{\rm br}}{\rho_0}, $$

to denote, respectively, the total number of plasmoids larger than the critical break size $w_{\rm br}$, the characteristic separation between them, and the corresponding particle energy.

We shall now consider the plasmoid distribution below $w_{\rm br}$, which will allow us to estimate the particle-trapping factors for $\gamma \leq \gamma _{\rm br}$. The normalization of the plasmoid distribution function $F(w)$ in this range is determined by the requirement that $F(w)$ be continuous at $w_{\rm br}$

(3.51)$$ F (w< w_{\rm br}) \simeq {\frac{\alpha_2-1}{w_{\rm br}}}N_{\rm br} \left(\frac{w}{w_{\rm br}}\right)^{-\alpha_1} \simeq {\frac{\alpha_2-1}{w_{\rm max}}}\left(\frac{w_{\rm br}}{w_{\rm max}}\right)^{-\alpha_2} \left(\frac{w}{w_{\rm br}}\right)^{-\alpha_1}. $$

The cumulative plasmoid distribution function below $w_{\rm br}$ is then determined by a straightforward integration of $F (w< w_{\rm br})$ with the condition that $N(w)$ be continuous at $w_{\rm br}$.

In particular, in the case $\alpha _1=1$ we get

(3.52)$$\begin{gather} F(w< w_{\rm br}) \simeq {\frac{\alpha_2-1}{w}} N_{\rm br}; \end{gather}$$

(3.53)$$\begin{gather}N(w\leq w_{\rm br}) = \int_w^{w_{\rm br}} F(w')\,{\rm d} w' + N(w_{\rm br}) \simeq N_{\rm br}\left[1 - (\alpha_2-1) \ln\frac{w}{w_{\rm br}}\right]; \end{gather}$$

(3.54)$$\begin{gather}\int^{\gamma<\gamma_{br}} N[w_{\rm tr}(\gamma')]\,{\rm d} \gamma' \simeq N_{\rm br} \gamma\left[\alpha_2 - (\alpha_2-1) \ln\frac{\gamma}{\gamma_{\rm br}}\right]. \end{gather}$$

Then, using the formalism developed in the previous two subsections (see (3.18)), we obtain the exponential factor describing the particle trapping in plasmoids, for $\gamma \leq \gamma _{\rm br}$

(3.55)$$\begin{align} & \exp \left(-\frac{1}{\dot\gamma_{\rm acc}} \int^{\gamma} \frac{{\rm d}\gamma'}{\tau_{\rm trap}(\gamma')}\right) = \exp \left[-\frac{1}{\beta_A}\frac{w_{\rm br}}{\epsilon L} \int^{\gamma} N[w_{\rm tr}(\gamma')]\,{\rm d} \left(\frac{\gamma'}{\gamma_{\rm br}}\right)\right] \nonumber\\ & \quad \simeq \exp \left[-\frac{1}{\beta_A} \frac{w_{\rm br}}{\epsilon \lambda_{\rm br}} \left(\frac{\gamma}{\gamma_{\rm br}}\right) \left[\alpha_2 - (\alpha_2-1) \ln\frac{\gamma}{\gamma_{\rm br}}\right]\right], \end{align}$$

where we used $\lambda _{\rm br} = L/N_{\rm br}$. Thus, ignoring the logarithmic correction, we see that there is simple exponential cutoff, $\exp (-\gamma /\gamma _c)$, with a cutoff energy that is of order

(3.56)$$ \gamma_c \sim \gamma_{\rm br} \beta_A \left({\frac{\epsilon \lambda_{\rm br}}{w_{\rm br}}}\right) = \epsilon \beta_A {\frac{\lambda_{\rm br}}{\rho_0}}. $$

Next, for the case $\alpha _1>1$, integration of (3.51) yields

(3.57)$$\begin{align} N(w\leq w_{\rm br}) = \int_w^{w_{\rm br}} F(w')\,{\rm d} w' + N(w_{\rm br}) \simeq N_{\rm br} \left({\frac{\alpha_2-1}{\alpha_1-1}}\right) \left[\left(\frac{w}{w_{\rm br}}\right)^{1-\alpha_1} - {\frac{\alpha_2-\alpha_1}{\alpha_2-1}}\right]. \end{align}$$

Then, using (3.57), we can calculate, for $\gamma \leq \gamma _{\rm br}$,

(3.58)$$\begin{align} \int^{\gamma} N[w_{\rm tr}(\gamma')]\,{\rm d} \gamma' & = \frac{1}{\rho_0} \int^{w_{\rm tr}(\gamma)} N(w') \,{\rm d} w' \nonumber\\ & \simeq \gamma_{\rm br} N_{\rm br} \left({\frac{\alpha_2-1}{\alpha_1-1}}\right) \left[\frac{1}{2-\alpha_1}\left(\frac{w}{w_{\rm br}}\right)^{2-\alpha_1} - {\frac{\alpha_2-\alpha_1}{\alpha_2-1}}\frac{w}{w_{\rm br}}\right]_{w=w_{\rm tr}(\gamma)} \nonumber\\ & = \gamma N_{\rm br} \left({\frac{\alpha_2-1}{\alpha_1-1}}\right) \left[\frac{1}{2-\alpha_1}\left(\frac{\gamma}{\gamma_{\rm br}}\right)^{1-\alpha_1} - {\frac{\alpha_2-\alpha_1}{\alpha_2-1}}\right] \nonumber\\ & = \gamma N_{\rm br} \left[{\frac{\alpha_2-1}{(\alpha_1-1)(2-\alpha_1)}} \left(\frac{\gamma}{\gamma_{\rm br}}\right)^{1-\alpha_1} - {\frac{\alpha_2-\alpha_1}{\alpha_1-1}}\right], \end{align}$$

and thus obtain the particle-trapping factor for $\gamma \leq \gamma _{\rm br}$

(3.59)$$\begin{align} & \exp \left(-\frac{1}{\dot\gamma_{\rm acc}} \int^{\gamma} \frac{{\rm d}\gamma'}{\tau_{\rm trap}(\gamma')}\right) = \exp \left(-\frac{1}{\epsilon \beta_A \varOmega_0}\frac{c}{L} \int^{\gamma} N[w_{\rm tr}(\gamma')] {\rm d}\gamma'\right) \nonumber\\ & \quad \simeq \exp \left(-\frac{1}{\beta_A} \frac{w_{\rm br}}{\epsilon \lambda_{\rm br}} \frac{\gamma}{\gamma_{\rm br}}\left[{\frac{\alpha_2-1}{(\alpha_1-1)(2-\alpha_1)}} \left(\frac{\gamma}{\gamma_{\rm br}}\right)^{1-\alpha_1} - {\frac{\alpha_2-\alpha_1}{\alpha_1-1}}\right]\right). \end{align}$$

We see that there are two cutoff factors here: one pure exponential, $\sim \exp (-\gamma /\gamma _{c1})$, with

(3.60)$$ \gamma_{c1} \simeq \gamma_{\rm br} \beta_A {\frac{\alpha_1-1}{\alpha_2-\alpha_1}} \left({\frac{\epsilon \lambda_{\rm br}}{w_{\rm br}}}\right) = \epsilon \beta_A {\frac{\alpha_1-1}{\alpha_2-\alpha_1}} {\frac{\lambda_{\rm br}}{\rho_0}} $$

(similar to what we have obtained above for the $\alpha _1=1$ case, apart from the logarithmic correction and factors of order unity), and the other sub-exponential, $\sim \exp [-(\gamma /\gamma _{c2})^{2-\alpha _1}]$, with

(3.61)$$ \gamma_{c2} \simeq \gamma_{\rm br}\left[\beta_A {\frac{(\alpha_1-1)(2-\alpha_1)}{\alpha_2-1}} \left({\frac{\epsilon \lambda_{\rm br}}{w_{\rm br}}}\right)\right]^{1/{(2-\alpha_1)}}. $$

As we can see, there is an important factor $\epsilon \lambda _{\rm br}/w_{\rm br}$ that governs the ratio of the cutoff to $\gamma _{\rm br}$ in all these expressions for any $\alpha _1$. Assuming first that $\alpha _2<2$, this factor can be evaluated, using (3.49), as

(3.62)$$ {\frac{\epsilon \lambda_{\rm br}}{w_{\rm br}}} \simeq {\frac{\epsilon L}{w_{\rm max}}}\left(\frac{w_{\rm max}}{w_{\rm br}}\right)^{2-\alpha_2} = {\frac{\epsilon L}{w_{\rm max}}}N_{\rm br}^{{(2-\alpha_2)}/{(\alpha_2-1)}}, $$

which is large if $w_{\rm max}\gg w_{\rm br}$ (and hence $N_{\rm br}\gg 1$). Substituting this into (3.60)–(3.61), we get

(3.63)$$ \gamma_{c1} \simeq \gamma_{\rm br} \beta_A {\frac{\alpha_1-1}{\alpha_2-\alpha_1}} \left({\frac{\epsilon L}{w_{\rm max}}}\right) \left(\frac{w_{\rm max}}{w_{\rm br}}\right)^{2-\alpha_2}, $$

and

(3.64)$$ \gamma_{c2} \simeq \gamma_{\rm br}\left[\beta_A {\frac{(\alpha_1-1)(2-\alpha_1)}{\alpha_2-1}} \left({\frac{\epsilon L}{w_{\rm max}}}\right) \left(\frac{w_{\rm max}}{w_{\rm br}}\right)^{2-\alpha_2}\right]^{1/{(2-\alpha_1)}}. $$

Thus we can see that, in the case of relativistic reconnection, $\beta _A = \mathcal {O}(1)$, the cutoffs are large compared with $\gamma _{\rm br}$ if $\alpha _2<2$; this means that particle trapping in plasmoids is not important for $\gamma < \gamma _{\rm br}$. For non-relativistic ($\beta _A \sim \sigma _h^{1/2} \ll 1$) reconnection, however, particle trapping may become important below $\gamma _{\rm br}$ (i.e. one of the cutoffs may drop below $\gamma _{\rm br}$) if the scale separation between $w_{\rm max}$ and $w_{\rm br}$ is not too large. Provided that the numerical coefficients $(\alpha _1-1)$, $(2-\alpha _1)$, $(\alpha _2-\alpha _1)$, etc., are not close to zero, we can regard them as finite constants of order unity; ignoring them for simplicity, we can then write the condition $\gamma _{c1,2} < \gamma _{\rm br}$ as

(3.65)$$ \beta_A {\frac{\epsilon L}{w_{\rm max}}} \left(\frac{w_{\rm max}}{w_{\rm br}}\right)^{2-\alpha_2} < 1. $$

If this condition is satisfied, we have $\gamma _{c2} < \gamma _{c1} < \gamma _{\rm br}$, i.e. the cutoff is dominated by $\gamma _{c2}$.

Next, we will investigate the relative importance of particle trapping above the break energy $\gamma _{\rm br}$. To do this, we need to evaluate the integral in (3.18) for $\gamma >\gamma _{\rm br}$. We will do this by splitting the whole region of integration into two sub-regions: up to $\gamma _{\rm br}$ and from $\gamma _{\rm br}$ to $\gamma$, and correspondingly obtain two multiplicative factors.

The first factor is a constant (independent of $\gamma$) that can be evaluated as

(3.66)$$\begin{gather} \exp \left(-\frac{c}{\epsilon \beta_A \varOmega_0 L} \int^{\gamma_{\rm br}} N[w_{\rm tr}(\gamma')]\,{\rm d}\gamma' \right) \simeq \exp \left[- \frac{\alpha_2 }{\beta_A}\frac{w_{\rm br}}{\epsilon \lambda_{\rm br}}\right] \quad {\rm for}\ \alpha_1=1, \end{gather}$$

(3.67)$$\begin{gather}\exp \left(-\frac{c}{\epsilon \beta_A \varOmega_0 L} \int^{\gamma_{\rm br}} N[w_{\rm tr}(\gamma')]\,{\rm d}\gamma'\right) \simeq \exp \left(-\frac{1}{\beta_A}\frac{w_{\rm br}}{\epsilon \lambda_{\rm br}} {\frac{1+\alpha_2-\alpha_1}{2-\alpha_1}}\right) \quad {\rm for}\ \alpha_1>1. \end{gather}$$

Once again, in light of (3.62), and adopting the assumptions $w_{\rm max}/\epsilon L = \mathcal {O}(1)$, $w_{\rm max}\gg w_{\rm br}$ and $\alpha _2<2$, we see that in both cases this exponential factor is unimportant for relativistic reconnection ($\beta _A \sim 1$). But it does become important for extremely non-relativistic reconnection, namely, when $\beta _A \simeq \sigma _h^{1/2}$ becomes comparable to $(w_{\rm br}/w_{\rm max})^{2-\alpha _2}$, see (3.65).

To evaluate the second factor, involving the integral of $\tau _{\rm tr}^{-1}(\gamma ')$ from $\gamma '=\gamma _{\rm br}$ to $\gamma >\gamma _{\rm br}$, we combine (3.16) and (3.46) to get

(3.68)$$\begin{align} \tau_{\rm tr}^{{-}1}(\gamma_{\rm br} < \gamma \leq \gamma_{\rm max}) & \simeq \frac{c}{L} N(w_{\rm br} < w \leq w_{\rm max})|_{w=w_{\rm tr}(\gamma)} \simeq \frac{c}{L} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{1-\alpha_2} \nonumber\\ & \simeq \frac{c}{L} N_{\rm br} \left(\frac{\gamma}{\gamma_{\rm br}}\right)^{1-\alpha_2}, \end{align}$$

where $\gamma _{\rm max} = w_{\rm max}/\rho _0$.

The resulting particle-trapping factor is then naturally independent of $\alpha _1$ and, for $\alpha _2<2$, is given by

(3.69)$$ \exp \left(-\frac{1}{\epsilon \beta_A \varOmega_0}\frac{c}{L} \int_{\gamma_{\rm br}}^{\gamma} N[w_{\rm tr}(\gamma')]\,{\rm d}\gamma' \right) \simeq \exp \left[-\frac{w_{\rm max}}{\beta_A \epsilon L} \frac{1}{2-\alpha_2} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{2-\alpha_2} \right], $$

which is, of course, identical to (3.35) obtained in the single power-law plasmoid distribution case. The resulting cutoff energy

(3.70)$$ \gamma_c = \gamma_{\rm max} \left[\beta_A {\frac{\epsilon L}{w_{\rm max}}}(2-\alpha_2)\right]^{1/{(2-\alpha_2)}} $$

can also be expressed in terms of $\gamma _{\rm br}$ as

(3.71)$$ \gamma_c =\gamma_{\rm br} \left[ \beta_A {\frac{\epsilon L}{w_{\rm max}}}(2-\alpha_2) \left(\frac{w_{\rm max}}{w_{\rm br}}\right)^{2-\alpha_2} \right]^{1/{(2-\alpha_2)}}. $$

We then see that $\gamma _c$ is above $\gamma _{\rm br}$ unless the inequality (3.65) is satisfied. If, however, the extreme non-relativistic condition (3.65) is satisfied, then $\gamma _c$ formally drops below $\gamma _{\rm br}$ (and even below $\gamma _{c1}$ and $\gamma _{c2}$, as one can see by by comparing (3.71) with (3.63)–(3.64) and ignoring order-unity factors like $(2-\alpha _2)$). This means that the particle distribution suffers even more severe suppression above $\gamma _{\rm br}$ in this case.

The case $\alpha _2=2$ is, again, special. Examining this case is motivated by analytical arguments (Uzdensky et al. Reference Uzdensky, Loureiro and Schekochihin2010) and numerical evidence (Loureiro et al. Reference Loureiro, Samtaney, Schekochihin and Uzdensky2012; Sironi et al. Reference Sironi, Giannios and Petropoulou2016) that this situation is indeed realized in practice. In this case, $\lambda _{\rm br} \simeq w_{\rm br} L/w_{\rm max}$ (see (3.49)), and hence

(3.72)$$ {\frac{\epsilon \lambda_{\rm br}}{w_{\rm br}}} \simeq {\frac{\epsilon L}{w_{\rm max}}}. $$

Thus, since we generally expect $w_{\rm max} \sim \epsilon L$, the cutoff of the small-$\gamma$ ($\gamma < \gamma _{\rm br}$) segment of the particle distribution (see (3.56) and (3.60)–(3.61)) can be estimated, ignoring factors of order unity and the logarithmic correction in the $\alpha _1=1$ case, as

(3.73)$$ \gamma_c \sim \gamma_{\rm br} \beta_A {\frac{\epsilon \lambda_{\rm br}}{w_{\rm br}}} \sim \gamma_{\rm br} \beta_A {\frac{\epsilon L}{w_{\rm max}}} \sim \gamma_{\rm br} \beta_A. $$

That is, the cutoff is comparable to $\gamma _{\rm br}$ for relativistic reconnection, and small compared with $\gamma _{\rm br}$ for non-relativistic reconnection. This indicates that the trapping of particles by plasmoids is important in this case even for moderate-energy ($\gamma \lesssim \gamma _{\rm br}$) relativistic particles. In particular, in the ultra-relativistic reconnection case, $\beta _A \simeq 1$, we get $\gamma _c \sim \gamma _{\rm br} = w_{\rm br}/\rho _0$, which can be recast as $\gamma _c \sim \sigma _c w_{\rm br}/\rho _{\sigma }$, where $\rho _{\sigma } \equiv \sigma _c \rho _0$. We then see that if $w_{\rm br}$ scales as a finite multiple of $\delta \sim \rho _{\sigma }$, then we get a result consistent with Werner et al. (Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016), i.e. the existence of an exponential cutoff proportional to $\sigma _c$, e.g. $\gamma _{c1} \simeq 4 \sigma _c$.

Furthermore, for particles above the break energy $\gamma _{\rm br}$, the trapping by large plasmoids manifests itself not via an exponential or quasi-exponential suppression factor but, just as in the case of a single power-law plasmoid chain with $\alpha =2$ (see (3.40)), via a power-law steepening due to an additional power-law factor

(3.74)$$ \exp \left(-\frac{1}{\epsilon \beta_A} \frac{\rho_0}{L} \int_{\gamma_{\rm br}}^{\gamma} N[w_{\rm tr}(\gamma')]\,{\rm d}\gamma' \right) \simeq \exp \left[-\frac{w_{\rm max}}{\beta_A\epsilon L} \ln\left(\frac{\gamma}{\gamma_{\rm br}}\right) \right] \simeq \left(\frac{\gamma}{\gamma_{\rm br}}\right)^{{-}p_{\rm tr}}, $$

with the same power-law index as in (3.40)

(3.75)$$ p_{\rm tr} \simeq \beta_A^{{-}1}\frac{w_{\rm max}}{\epsilon L}. $$

Once again, assuming that $w_{\rm max} \sim \epsilon L$, we see that this additional power-law index resulting from particle trapping in plasmoids scales as $p_{\rm tr} \sim \beta _A^{-1} \sim (1+1/\sigma _h)^{1/2}$.

To sum up, for $\alpha _2 < 2$ and relativistic reconnection ($\beta _A \sim 1$) we find that the existence of the power-law break in the plasmoid distribution at $w_{\rm br} \ll w_{\rm max}$ does not result in any significant changes to the particle energy distribution relative to the case of a single power law with $\alpha = \alpha _2$. However, in the case of non-relativistic reconnection ($\beta _A \ll 1$), the presence of the break in the plasmoid spectrum $F(w)$ may lead to a modification in the particle distribution if the dynamic range $w_{\rm max}/w_{\rm br}$ of the second ($\alpha = \alpha _2$) power-law segment of $F(w)$ is not too large, namely if the condition (3.65) is satisfied. In this case, the high-energy cutoff of the particle distribution falls below $\gamma _{\rm br}$ and is given by (3.64).

In the special case $\alpha _2 =2$ we expect a power-law distribution, with the usual index $p_{\rm magn}$ scaling inversely with $\beta _A$, at moderately suprathermal energies ($\gamma \lesssim \gamma _{\rm br}$), checked by an exponential or quasi-exponential cutoff (3.73). At higher energies, above $\gamma _{\rm br}$, the distribution transitions to a second, steeper power law with index $p_{\rm tot} = p_{\rm magn} + p_{\rm tr}$. The cutoff $\gamma _c$ connecting the two power laws is of order $\gamma _{\rm br} \beta _A$, and is thus fundamentally controlled by the break $w_{\rm br} = \gamma _{\rm br} \rho _0$ in the plasmoid-size distribution; in particular, it can be much smaller than the extreme-acceleration Hillas limit of $\sim L/\rho _0$, and may even just scale with $\sigma _c \beta _A$ if $w_{\rm br} \sim \sigma _c \rho _0$. Moreover, in the case of strongly non-relativistic reconnection, $\beta _A \sim \sigma _h^{1/2} \ll 1$, the cutoff $\gamma _c$ becomes much smaller than $\gamma _{\rm br}$, thus greatly reducing the dynamic range of the first power-law segment of the particle distribution. This has a strong effect on the normalization of the second, high-energy power-law segment above $\gamma _{\rm br}$: its normalization is suppressed by the exponential factors (3.66)–(3.67); the suppression is just by a finite factor of order unity in the ultra-relativistic case $\beta _A \simeq 1$ but could be exponentially strong, viz. $\exp (-1/\beta _A)$, in the non-relativistic case $\beta _A \ll 1$. This may explain the difficulties in obtaining robust extended non-thermal relativistic power laws in numerical PIC studies of 2-D non-relativistic magnetic reconnection.

3.7. Fermi acceleration by reflection off moving plasmoids

In this subsection we discuss the diffusive Fermi acceleration due to a particle's reflection off moving plasmoids that fly chaotically up and down along the plasmoid chain (see also Guo et al. Reference Guo, Liu, Daughton and Li2015; Nalewajko et al. Reference Nalewajko, Uzdensky, Cerutti, Werner and Begelman2015). We stress again that this acceleration avenue is distinct from the above-considered acceleration by the main (regular) reconnection electric field $E_{\rm rec}$ responsible for the growth of plasmoids (see § 3.2). Indeed, since the Fermi acceleration process under consideration here is related to the relative motions of plasmoids along the layer (i.e. in the $x$-direction), it is fundamentally driven by the nonlinear development of the coalescence instability, as opposed to the tearing instability associated with $E_{\rm rec}$. Building on the qualitative discussion of the underlying physics of this process in § 2.2, here we develop a more quantitative description.

The detailed microphysical underpinning of this acceleration process can be described as follows (see figure 3). When a particle in its motion in the $(xz)$-plane encounters a large plasmoid moving (relative to the laboratory frame) with a speed $v_{\rm pl}$, it interacts with the plasmoid's motional electric field of magnitude $|E_z| \sim E_{\rm pl} = v_{\rm pl} B_{\rm pl} /c$; the sign of this electric field can be positive or negative depending on the plasmoid's direction of motion. During this interaction, as discussed in § 2.2, the particle completes a $\gtrsim 1/2$ fraction of its cyclotron orbit in the plasmoid's magnetic field, thus moving backward in the $z$-direction by about the orbit's diameter, $|\Delta z| \sim 2 \gamma m_e c^2/eB_{\rm pl}$. Thus, the magnitude of the resulting particle's energy change due to the work done on it by the electric field $E_{\rm pl}$ during this interaction can be estimated as

(3.76)$$ |\Delta \epsilon| \simeq e E_{\rm pl} |\Delta z| \simeq (v_{\rm pl}/c) \gamma m_e c^2 = (v_{\rm pl}/c) \epsilon. $$

For simplicity, we shall take the characteristic plasmoid velocities to be independent of plasmoid size, namely, of order $v_{\rm pl} \sim V_A = \beta _A c$ for all plasmoids.Footnote ¹³ Then, the particle's fractional energy change is of order

(3.77)$$ \frac{|\Delta \epsilon|}{\epsilon} = \frac{|\Delta \gamma|}{\gamma} \simeq \frac{v_{\rm pl}}{c} \sim \beta_A, $$

which is consistent with the result of an elastic collision in the plasmoid frame.

As discussed in § 2.2, the sign of this energy change can be positive or negative depending on whether the encounter is head-on or tail-on. If we ignore correlations in the neighbouring large plasmoids’ directions of motionFootnote ¹⁴ and thus regard the scatterings as random, the particle's energy undergoes a random walk, resulting in diffusive second-order Fermi acceleration.

For a high-energy particle with a Lorentz factor $\gamma \gg \bar {\gamma }$, the characteristic time step of this random walk, i.e. the interval between successive encounters with large plasmoids, is approximately the particle's flight time over the typical inter-plasmoid distance between such large plasmoids, i.e. $\Delta t (\gamma ) \sim \lambda _{\rm pl}[w_{\rm tr} (\gamma )]/c$. Then, the corresponding energy diffusion coefficient can be estimated, ignoring factors of order unity, as

(3.78)$$ D_\gamma(\gamma) = {\frac{\Delta \gamma^2}{\Delta t}} \sim \beta_A^2 \gamma^2 \frac{c}{\lambda_{\rm pl}[w_{\rm tr}(\gamma)]} = \beta_A^2 \gamma^2 \frac{c}{L}N_{\rm pl}[w_{\rm tr}(\gamma)] = \beta_A^2 \gamma^2 \tau_{\rm trap}^{{-}1}(\gamma), $$

and the corresponding term in the kinetic equation (3.1) or (3.2) becomes

(3.79)$$ \partial_\gamma (D_\gamma \partial_\gamma f) \sim \beta_A^2 \frac{c}{L}\partial_\gamma \left( \gamma^2 N_{\rm pl}[w_{\rm tr}(\gamma)] \partial_\gamma f \right) = \beta_A^2 \partial_\gamma \left(\gamma^2 \tau_{\rm trap}^{{-}1}(\gamma) \partial_\gamma f \right). $$

We can get further insight by defining the nominal characteristic diffusive Fermi acceleration time

(3.80)$$ \tau_{\rm dif}(\gamma) \equiv \frac{\gamma^2}{D_\gamma(\gamma)} \sim \beta_A^{{-}2} \frac{\lambda_{\rm pl}[w_{\rm tr}(\gamma)]}{c} = \beta_A^{{-}2} \frac{L}{c} \frac{1}{N_{\rm pl}[w_{\rm tr}(\gamma)]}, $$

and comparing it with the characteristic time scales corresponding to the other terms in the kinetic equation (3.2).

First, the ratio of $\tau _{\rm dif}$ to the particle escape time due to trapping by plasmoids is particularly simple and is independent of the particle energy

(3.81)$$ \frac{\tau_{\rm dif}(\gamma)}{\tau_{\rm trap}(\gamma)} \sim \beta_A^{{-}2} = 1+1/\sigma_h. $$

Thus, the two time scales are automatically comparable in the case of relativistic reconnection, $\sigma _h \gtrsim 1$, $\beta _A \simeq 1$; but for non-relativistic reconnection, $\sigma _h \ll 1$, we find that $\tau _{\rm dif} \gg \tau _{\rm trap}$.

Next, by comparing (3.80) with the magnetization time scale $\tau _{\rm magn}$ (see (3.9)) and the regular acceleration time scale $\tau _{\rm reg}(\gamma ) \equiv \gamma /\dot {\gamma }_{\rm acc} \sim \rho (\gamma )/c \epsilon \beta _A$ (see (3.5)), we find

(3.82)$$ \frac{\tau_{\rm dif}(\gamma)}{\tau_{\rm magn}(\gamma)} \sim \beta_A^{{-}2} \frac{\lambda_{\rm pl}[w_{\rm tr}(\gamma)]}{c} \frac{\epsilon\varOmega_0}{\gamma} = \beta_A^{{-}2} \frac{\epsilon\lambda_{\rm pl}[w_{\rm tr}(\gamma)]}{\rho(\gamma)} \sim \beta_A^{{-}2} \left.\frac{\epsilon\lambda_{\rm pl}(w)}{w}\right|_{w=w_{\rm tr}(\gamma)}, $$

and

(3.83)$$ \frac{\tau_{\rm dif}(\gamma)}{\tau_{\rm reg}(\gamma)} \sim \beta_A^{{-}2} \frac{\lambda_{\rm pl}[w_{\rm tr}(\gamma)]}{c} \frac{c \epsilon\beta_A}{\rho(\gamma)} \simeq \beta_A^{{-}1} \frac{\epsilon\lambda_{\rm pl}[w_{\rm tr}(\gamma)]}{\rho(\gamma)} \sim \beta_A^{{-}1} \left.\frac{\epsilon\lambda_{\rm pl}(w)}{w}\right|_{w=w_{\rm tr}(\gamma)}. $$

We can make two observations from these two expressions. First, the last factor in these expressions, $\epsilon \lambda _{\rm pl}(w)/w = \epsilon L/[wN(w)]$, already familiar from the discussion in § 3.6, manifestly underscores that the relative importance of Fermi acceleration (compared with particle magnetization by $B_1$ and regular acceleration by $E_{\rm rec}$) depends on the plasmoid distribution function $N(w)$. Thus, further analysis requires adopting a specific choice for $N(w)$. For concreteness, let us consider the case of a single power-law plasmoid distribution (3.25) with $\alpha >1$ (see § 3.5): $N\simeq (w/w_{\rm max})^{1-\alpha }$. We then have

(3.84)$$\begin{gather} D_\gamma(\gamma) \sim \beta_A^2 \gamma_{\rm max}^2 \frac{c}{L} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{3-\alpha}, \end{gather}$$

(3.85)$$\begin{gather}\tau_{\rm dif}(\gamma)\sim \beta_A^{{-}2} \frac{L}{c} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{\alpha-1}, \end{gather}$$

and the factor $\epsilon \lambda _{\rm pl}(w)/w$, controlling the ratios of $\tau _{\rm dif}$ to $\tau _{\rm magn}$ and $\tau _{\rm reg}$, can be written as

(3.86)$$ {\frac{\epsilon\lambda_{\rm pl}(w)}{w}} \simeq \frac{\epsilon L}{w} \left(\frac{w}{w_{\rm max}}\right)^{\alpha-1} = \frac{\epsilon L}{w_{\rm max}} \left(\frac{w}{w_{\rm max}}\right)^{\alpha-2}. $$

Substituting this estimate into (3.82)–(3.83), we obtain

(3.87)$$ \frac{\tau_{\rm dif}(\gamma)}{\tau_{\rm magn}(\gamma)} \sim \beta_A^{{-}2}\frac{\epsilon L}{w_{\rm max}} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{\alpha-2}, $$

and

(3.88)$$ \frac{\tau_{\rm dif}(\gamma)}{\tau_{\rm reg}(\gamma)} \sim \beta_A^{{-}1}\frac{\epsilon L}{w_{\rm max}} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{\alpha-2}. $$

Considering first the case $\alpha <2$ and adopting a reasonable assumption that $w_{\rm max}\sim \epsilon L \sim L/10$, we see that the ratio $\epsilon \lambda _{\rm pl}(w)/w$ given by (3.86) can be expected to be of order unity only for the biggest plasmoids, $w\sim w_{\rm max}$. For the majority of plasmoids, however, $w \ll w_{\rm max}$, and this ratio becomes large. Then, since $\beta _A \leq 1$, (3.87)–(3.88) imply that, for the majority of the non-thermal particles, namely those with $\gamma \ll \gamma _{\rm max}$, the diffusive Fermi acceleration by moving plasmoids is relatively unimportant compared with the regular acceleration by $E_{\rm rec}$ and magnetization by $B_1$, and, furthermore, that its relative importance decreases as the particle energy is lowered. In other words, while smaller plasmoids, which are capable of scattering lower-energy particles, are more numerous, they are not numerous enough as long as the plasmoid spectrum is sufficiently shallow, $\alpha < 2$. On the other hand, however, for the highest-energy particles with $\gamma \sim \gamma _{\rm max}$ (corresponding to $w_{\rm tr}(\gamma ) \sim w_{\rm max}$) the diffusive Fermi acceleration by moving plasmoids may still be important, provided that the $\beta _A^{-1}$ and $\beta _A^{-2}$ factors can be circumvented (see below).

Next, let us consider the practically important case $\alpha =2$, which is, again, special. In this case, (3.86) tells us that $\epsilon \lambda _{\rm pl}/w$ is independent of $w$, and hence the ratios $\tau _{\rm dif}/\tau _{\rm magn}$ and $\tau _{\rm dif}/\tau _{\rm reg}$ are the same for all particles in the non-thermal power law, and are just governed by $\beta _A$. For reference, the diffusion coefficient and the nominal diffusive acceleration time in this case are given by

(3.89)$$\begin{gather} D_\gamma^{\alpha=2}(\gamma)\sim \beta_A^2 \gamma_{\rm max}^2 \frac{c}{L}\frac{\gamma}{\gamma_{\rm max}}, \end{gather}$$

(3.90)$$\begin{gather}\tau_{\rm dif}^{\alpha=2}(\gamma)\sim \beta_A^{{-}2} \frac{L}{c} \frac{\gamma}{\gamma_{\rm max}}. \end{gather}$$

The second observation we can make by looking at (3.82)–(3.83), as well as (3.81), is that the relative strength of the Fermi acceleration appears to be suppressed in the non-relativistic ($\beta _A\ll 1$) reconnection case relative to the relativistic ($\beta _A \simeq 1$) case by the factors of $\beta _A$ or $\beta _A^2$. We point out, however, that it is not quite fair to judge the relative importance of Fermi acceleration just by the ratios of $\tau _{\rm dif}$ to the other characteristic time scales in the kinetic equation, especially in the non-relativistic case. The reason for this is that the diffusion term involves higher-order (namely, second-order) derivatives of $f(\gamma )$ than the other terms. Thus, if $f(\gamma )$ has a strong dependence on $\gamma$ (e.g. declines rapidly), then this term may become important (although not dominant). In particular, for a power-law spectrum $f(\gamma ) \sim \gamma ^{-p}$, corresponding to a power-law plasmoid distribution with an index $\alpha$, we can estimate the diffusion term (3.79), using (3.84), as

(3.91)$$ \partial_\gamma (D_\gamma \partial_\gamma f) \sim p (\alpha+p-2)\beta_A^2 \frac{c}{L} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{1-\alpha} f(\gamma) \sim p (\alpha+p-2)\frac{f(\gamma)}{\tau_{\rm dif}(\gamma)}, $$

which thus can be significantly greater than just the naive estimate $f/\tau _{\rm dif}$ if $p \gg 1$. This is in fact expected in the non-relativistic reconnection case, where $p \sim \beta _A^{-1} \gg 1$ (see (3.11) and (3.41)); in this case, the diffusion term is by a factor $\beta _A^{-2} \simeq \sigma _h^{-1} \gg 1$ greater than $f/\tau _{\rm dif}$ if $p \gg 1$, and becomes automatically comparable to the plasmoid-trapping term.

In the special case $\alpha =2$ expression (3.91), for any $\beta _A$, simplifies to

(3.92)$$ \partial_\gamma (D_\gamma \partial_\gamma f)|_{\alpha=2} \sim p^2 \beta_A^2 \frac{c}{L} \frac{\gamma_{\rm max}}{\gamma} f(\gamma) \sim p^2 \frac{f(\gamma)}{\tau_{\rm dif}(\gamma)}, $$

which is, again, by a factor $p^2 \sim \beta _A^{-2} \sim \sigma _h^{-1}$ greater than $f/\tau _{\rm dif}$ in the non-relativistic limit.

We can then see that the Fermi diffusion term may affect our calculation of the particle spectrum power-law index $p$. To remind the reader, in our theory this index is governed, in the absence of Fermi acceleration, by the interplay between the regular acceleration $\dot {\gamma }_{\rm acc}$ by the main reconnection electric field $E_{\rm rec}$ and particle magnetization by the reconnected field $B_1$ (see §§ 3.2 and 3.3). And in the special case $\alpha =2$, particle trapping in plasmoids also affects the power-law index, resulting in a moderate steepening (see § 3.5 and, in particular, (3.41)). We can now examine the role of diffusive Fermi acceleration in the balance between all these terms, which controls $p$.

To do this, let us consider the steady-state kinetic equation (3.2) and focus on the power-law part of the particle distribution function, i.e. $f\sim \gamma ^{-p}$, thus ignoring a possible high-energy cutoff due to trapping by large plasmoids. Then, considering first the case $1< \alpha <2$, and hence excluding the plasmoid-trapping term, we can write this equation as

(3.93)$$ \epsilon \beta_A \varOmega_0 p \frac{f}{\gamma} - \epsilon \varOmega_0 \frac{f}{\gamma}+ p (\alpha+p-2)\beta_A^2 \frac{c}{L}\left(\frac{\gamma}{\gamma_{\rm max}}\right)^{1-\alpha} f(\gamma) = 0, $$

where we have used (3.5) and (3.9) to evaluate the first two terms (regular acceleration by $E_{\rm rec}$ and magnetization by $B_1$). Cancelling $f$, multiplying through by $\gamma /\epsilon \varOmega _0$, and using $w_{\rm max} = \gamma _{\rm max} \rho _0 = \gamma _{\rm max} c \varOmega _0^{-1}$, we obtain

(3.94)$$ \beta_A p - 1 + p (\alpha+p-2)\beta_A^2 \frac{w_{\rm max}}{\epsilon L} \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{2-\alpha} = 0. $$

Since we here assume $\alpha < 2$, we again see that, with the exception of the most energetic particles with $\gamma \sim \gamma _{\rm max}$, the last term, corresponding to Fermi acceleration by moving large plasmoids, is small – consistent with our discussion after (3.88) above. We then get a familiar result $p \simeq \beta _A^{-1}$.

Next let us consider the special case $\alpha =2$. In this case, the problem can be solved exactly. The particle trapping by plasmoids now affects the power-law index instead of introducing a high-energy cutoff (see § 3.5), and needs to be included in the above kinetic equation (3.93). Using (3.33) to evaluate this term, (3.93) is modified as

(3.95)$$ \epsilon \beta_A \varOmega_0 p \frac{f}{\gamma} -\epsilon \varOmega_0 \frac{f}{\gamma}- \frac{c}{L}\left(\frac{\gamma_{\rm max}}{\gamma}\right) f(\gamma) + p^2 \beta_A^2 \frac{c}{L}\left(\frac{\gamma_{\rm max}}{\gamma}\right) f(\gamma) = 0, $$

which simplifies to the following quadratic equation for $p\beta _A$, involving a single dimensionless parameter $w_{\rm max}/{\epsilon L}$:

(3.96)$$ p \beta_A - \left(1 + \frac{w_{\rm max}}{\epsilon L}\right) + p^2 \beta_A^2 \frac{w_{\rm max}}{\epsilon L} = 0. $$

The positive solution of this equation, for any $w_{\rm max}/{\epsilon L}$, is simply

(3.97)$$ p^{\alpha=2} = \beta_A^{{-}1}, $$

coinciding with (3.11). Thus, in the $\alpha =2$ case, diffusive Fermi acceleration by rapidly moving plasmoids effectively cancels, or negates, the steepening of the particle distribution due to particle trapping in plasmoids described by (3.41). Note, however, that in reality the diffusion term may not lead to an exact complete cancellation of the trapping effect because, in a more precise and realistic model, both of these terms would in general come in with some numerical coefficients of order unity, which don't have to be equal. We discuss this more general situation in the next subsection.

3.8. General form of the kinetic equation

Finally, for completeness and to set the stage for future, more detailed studies of this problem, in this subsection we present, for reference, the general second-order linear ordinary differential equation for the stationary particle distribution function including all the terms, for an arbitrary plasmoid distribution. To obtain it, we plug in the estimates (3.5), (3.9) and (3.79) into the steady-state kinetic equation (3.2). Note however, that, generally speaking, all these estimates are uncertain up to factors of order unity; we will treat these uncertainties by assigning some unknown positive constant coefficients, $c_1$, $c_2$, $c_3$, to the last three terms when we enter these expressions into (3.2). We then obtain

(3.98)$$ f'(\gamma) + c_1 \frac{1}{\beta_A} \frac{f}{\gamma} + c_2 \frac{1}{\beta_A} \frac{f}{\epsilon\varOmega_0 \tau_{\rm trap}(\gamma)} - c_3 \beta_A \left(\gamma^2 \frac{1}{\epsilon\varOmega_0 \tau_{\rm trap}(\gamma)} f'\right)'=0, $$

where prime denotes derivative with respect to $\gamma$.

Using $\tau _{\rm trap}(\gamma ) \simeq \lambda _{\rm pl}[w_{\rm tr}(\gamma )]/c$, and $c/\varOmega _0 = \rho _0 = w_{\rm tr}(\gamma ) /\gamma$, it is convenient to define a new function encapsulating the information about the plasmoid distribution

(3.99)$$ q(\gamma) \equiv \left.\frac{w}{\epsilon \lambda_{\rm pl}(w)} \right\vert_{w=w_{\rm tr}(\gamma) = \gamma \rho_0} = \frac{\gamma}{\epsilon\varOmega_0 \tau_{\rm trap}(\gamma)}. $$

Then, (3.98) becomes

(3.100)$$ f'(\gamma) + \frac{1}{\beta_A} \frac{f}{\gamma} [c_1 + c_2 q(\gamma)] - c_3 \beta_A \left( \gamma q(\gamma) f'\right)'=0. $$

Using the substitution $\gamma = {\rm e}^{\beta _A t}$ and $\tilde {q}(t) = q[\gamma (t)]$, we can rewrite this as

(3.101)$$ \dot{f} (t) + f(t)[c_1 + c_2 \tilde{q}(t)] - c_3 \frac{{\rm d}}{{\rm d} t}\left( \tilde{q}(t) \dot{f} \right) =0, $$

where $\dot {f} \equiv {\rm d} f/{\rm d} t$.

In particular, if the plasmoid distribution is a single power law with $\alpha >1$, then, using (3.32)

(3.102)$$ q(\gamma) = q_0 \left(\frac{\gamma}{\gamma_{\rm max}}\right)^{2-\alpha} = q_0 \exp({(2-\alpha)\beta_A (t-t_{\rm max})}), $$

where we defined $q_0 = w_{\rm max}/\epsilon L$, and $t_{\rm max} \equiv \beta _A^{-1}\ln \gamma _{\rm max}$. Equation (3.101) then transforms into

(3.103)$$\begin{align} & \dot{f} (t) + f(t)[c_1 + c_2 q_0 \exp({(2-\alpha)\beta_A (t - t_{\rm max})})] \nonumber\\ & \quad - c_3\frac{{\rm d}}{{\rm d} t}\left( q_0 \exp({(2-\alpha)\beta_A (t - t_{\rm max})})\dot{f} \right) =0. \end{align}$$

For example, in the special case $\alpha = 2$, we get $q(\gamma ) = \tilde {q}(t) = q_0 = {\rm const}$, (3.100) becomes an Euler equation

(3.104)$$ f'(\gamma) + \frac{1}{\beta_A} \frac{f}{\gamma} [c_1 + c_2 q_0 ] - c_3 \beta_A \left(q_0\gamma f'\right)'=0, $$

while the corresponding equation (3.103) simplifies to a linear homogeneous ordinary differential equation with constant coefficients

(3.105)$$ \dot{f} (t) + f(t) [c_1 + c_2 q_0] - c_3 q_0 \ddot{f} =0. $$

This can be readily solved as

(3.106)$$ f(t) = C_1 \,{\rm e}^{\lambda_1 t} + C_2 \,{\rm e}^{\lambda_2 t}, $$

where $C_1$, $C_2$ are arbitrary constants and $\lambda _1$ and $\lambda _2$ are the roots of the characteristic polynomial

(3.107)$$ P(\lambda) = \lambda^2 - \frac{1}{c_3 q_0} \lambda - {\frac{c_1+c_2 q_0}{c_3 q_0}}. $$

For illustration, for the case $c_1 = c_2 = c_3 = 1$ considered in this paper, these roots are $\lambda _1 =-1$, $\lambda _2 = 1+1/q_0$, and the corresponding solution is

(3.108)$$ f = C_1 \,{\rm e}^{\lambda_1 t} + C_2 \,{\rm e}^{\lambda_2 t} = C_1 \gamma^{\lambda_1/\beta_A} + C_2 \gamma^{\lambda_2/\beta_A} = C_1 \gamma^{-({1/{\beta_A}})} + C_2 \gamma^{{(1+q_0)}/{q_0 \beta_A}}. $$

Imposing the condition that $f(\gamma )$ be a declining function of $\gamma$, we have to discard the second solution and thus recover the result $f(\gamma ) \sim \gamma ^{-1/\beta _A}$ we obtained previously, see (3.97). We note that the solution $\lambda _1 =-1$, $f(\gamma ) \sim \gamma ^{-1/\beta _A}$ remains valid even more generally, as long as $c_1 = 1$ and $c_2=c_3$.