1 Introduction↩︎

Strong magnetic fields in astrophysical compact sources provide a reservoir of magnetic energy. This energy can be released to the plasma—resulting in particle acceleration and nonthermal emission—when anti-aligned field lines annihilate in a process called magnetic reconnection. In a number of astrophysical sources, reconnection occurs in the relativistic regime, where the magnetic energy exceeds the plasma rest mass energy [1]–[4]. Relativistic reconnection can power a variety of high-energy phenomena, such as emission from black hole coronae, magnetar flares, blazar jet flares, radio and gamma-ray emission from pulsar magnetospheres, fast radio bursts, and flares from supermassive black hole magnetospheres.

Magnetic reconnection refers to the breaking and reconnecting of oppositely directed field lines. This requires the “ideal” condition \[{\boldsymbol{E}}+\frac{\langle{\boldsymbol{v}}_s\rangle}{c}\times {\boldsymbol{B}}=0\] to be violated for each relevant plasma species \(s\). Here, \({\boldsymbol{E}}\) and \({\boldsymbol{B}}\) are the electromagnetic fields, while \(\langle {\boldsymbol{v}}_s\rangle\) is the mean three-velocity of species \(s\). In collisional plasmas, the ideal condition can be broken by resistive effects due to binary particle collisions—encoded by the resistivity appearing in Ohm’s law. When resistive effects are not important, the magnetic field is “frozen” into the fluid, as prescribed by Alfvén’s theorem ([5]; also known as the flux freezing theorem). In dilute astrophysical plasmas, binary collisions are rare, so the collisional resistivity is often insufficient to break flux freezing on interesting time and length scales.

Reconnection occurring in the collisionless regime requires a kinetic description. Since the typical separation between plasma scales and global scales is very large, kinetic descriptions, e.g., employing the particle-in-cell (PIC) method, are unaffordable at realistic scale separations. Fluid-type approaches such as magnetohydrodynamics (MHD), while suitable to model the global dynamics, are by construction collisional, and therefore unable to capture collisionless effects. In fact, single-fluid resistive MHD simulations with uniform resistivity yield reconnection rates in the plasmoid-dominated regime that are much lower (by nearly one order of magnitude) than equivalent kinetic calculations [6]–[9]. This discrepancy impacts the timescale of reconnection-powered flares, e.g., in black hole magnetospheres [10], [11].

A large body of work has focused on identifying the processes that can break the ideal condition in collisionless or weakly collisional plasmas—here, wave-particle interactions provide a form of effective collisionality. In pair plasmas, fast reconnection is mediated by the off-diagonal terms of the pressure tensor [12]–[16], which are also important for electron-ion plasmas in the small, electron-scale diffusion region [17]–[21].

By identifying the dominant contributors to the breaking of flux freezing in collisionless plasmas, it may be possible to write the corresponding nonideal electric field as \(\eta_{\rm eff} {\boldsymbol{J}}\)—here, \(\eta_{\rm eff}\) is some effective resistivity and \({\boldsymbol{J}}\) the electric current density—which could be incorporated in resistive, single-fluid MHD approaches as a kinetically motivated subgrid prescription [22]–[29]. In general, parameterizing kinetic effects as an effective resistivity is a nontrivial task (e.g., [30], [31]). By means of a statistical analysis based on PIC simulations, [32] identified the mechanisms driving the nonideal electric field in the generalized Ohm’s law, for the case of relativistic pair plasma reconnection. The effective resistivity proposed by [32] for the zero guide field case (and earlier suggested by [13], [33]) has been shown to successfully enhance the reconnection rate in resistive MHD simulations [34].

As we discuss below, the form of effective resistivity proposed by [32] suffers from a few limitations, which may hamper its applicability. In this work, rather than analyzing nonideal terms in the generalized Ohm’s law, we adopt an empirical approach. We perform a suite of PIC simulations of relativistic pair-plasma reconnection with varying guide field strength, and we formulate an empirical prescription for the effective resistivity \(\eta_{\rm eff}\), which is derived directly from our PIC runs through a data-driven parameterization. Our proposed model depends only on the electric current density and the plasma number density, both of which are readily available in resistive MHD codes. As compared to [32], the form of \(\eta_{\rm eff}\) that we obtain has four main advantages: it is written explicitly in single-fluid MHD quantities, does not depend on spatial derivatives, is coordinate-agnostic, and is valid for any guide field. We demonstrate that the formulation of \(\eta_{\rm eff}\) we propose successfully reproduces the spatial structure and strength of nonideal electric fields in our PIC simulations, thus providing a promising strategy for enhancing the reconnection rate in resistive MHD approaches.

2 PIC Simulation Setup↩︎

Our simulations are performed with the 3D PIC code TRISTAN-MP [35]. We use a 2D \(x-y\) domain, but we track all components of the particle velocity and of the electromagnetic fields. Although the physics of particle acceleration in relativistic reconnection is dramatically different between 2D and 3D [36], [37], the reconnection rate—which is defined as the plasma inflow velocity and is the focus of our work—is roughly the same [38], [39].

The in-plane magnetic field is initialized in a “Harris equilibrium” [40], \({\boldsymbol{B}_{\rm in}} = B_0{\boldsymbol{\hat{x}}}\tanh(2\pi y/\Delta)\), where the direction of the in-plane field reverses at \(y = 0\) over a thickness \(\Delta = 70\,c/\omega_{p}\). Here, \(c/\omega_p\) is the depth of the plasma skin, and \(\omega_p=\sqrt{4\pi n_0e^2/m}\) the plasma frequency, where \(n_0\) is the total number density of electron-positron pairs far from the layer, \(m\) the electron / positron mass and \(e\) the elementary charge. We parameterize the field strength in the plane \(B_0\) by the magnetization \[\sigma = \frac{B_0^2}{4\pi n_0 m c^2}~,\] which we take to be \(\sigma=50\). We consider guide fields of magnitude \(B_g/B_0 = 0.0, 0.3, 0.6, 1.0\).

The upstream region is initialized with \(n_0=64\) particles per cell (including both species). We resolve the plasma skin depth \(c/\omega_p\) with five cells and evolve the simulation up to \(4500 \, \omega_p^{-1}\). In 6, we choose \(c/\omega_p=20\) cells and demonstrate that our results are robust to spatial resolution. In 6 we also display cases that include strong synchrotron cooling losses. For our fiducial runs, the length of the domain in the \(x\)-direction of plasma outflows is \(L_x = 1920\,c/\omega_p\). We use open boundaries for fields and particles along the \(x\)-direction. The box grows in the \(y\)-direction as the simulation progresses, allowing for more plasma and magnetic flux to enter the domain. At the end of the simulations, the length of our box along the \(y\)-axis is comparable to \(L_{x}\).

3 Resistivity Formulation↩︎

We begin by considering the Ohm’s law for resistive relativistic single-fluid MHD [41]: \[\label{eq:32ohm95full} \Gamma \left[\mathbf{E}+\frac{\mathbf{v}}{c} \times \mathbf{B} - \frac{1}{c^2}(\mathbf{E}\cdot \mathbf{v})\mathbf{v} \right] = \eta \left(\mathbf{J}-\rho_e \mathbf{v} \right)\tag{1}\] where \(\Gamma\) is the bulk fluid Lorentz factor, \(\mathbf{v}\) the fluid three-velocity, \(\rho_e\) the electric charge density, and \(\eta\) the collisional resistivity. In a collisionless plasma, the replacement of \(\eta\) by \(\eta_{\mathrm{eff}}\) in 1 can be regarded as the definition of an effective resistivity that incorporates kinetic effects in single-fluid resistive MHD. As we justify in 7, we can further assume that \(|\rho_e {\boldsymbol{v}}|\ll |{\boldsymbol{J}}|\) and \(\Gamma\simeq 1\), and that the third term in the square bracket is negligible, which yields \[\label{Eq:eta} \mathbf{E}^*\equiv \mathbf{E}+\frac{\mathbf{v}}{c} \times \mathbf{B} = \eta_{\mathrm{eff}} \mathbf{J}\tag{2}\] where \(\mathbf{E}^*\) is the nonideal electric field. The spatial structure of the \(z\)-component of the nonideal electric field, \(E_z^*\), is shown in the bottom row of 1, at a representative time after the simulation has achieved a quasi-steady state (i.e., the reconnection rate attains a quasi-steady value). The figure emphasizes that nonideal regions are generally larger for increasing guide field. We also present the spatial structure of the total particle density \(n_t\) (top row; in units of \(n_0\)) and of the magnetic energy density (middle row; in units of \(B_0^2/8\pi\)), for both \(B_g/B_0=0\) (left column) and \(B_g/B_0=1\) (right column).

In order to determine the effective resistivity \(\eta_{\rm eff}\), we focus on the \(z\)-component of 2 , which dominates the nonideal field for the whole range of \(B_g/B_0\) we explore. The \(z\)-component \(E_z^*\) is the only significant component for zero guide field, being one to two orders of magnitude larger than other components; for non-zero guide fields, we still determine \(\eta_{\rm eff}\) from \(E_z^*=\eta_{\rm eff}J_z\), but we show that the same effective resistivity properly describes other components, specifically \(E_x^*=\eta_{\rm eff}J_x\) (see 8). Our new prescription for the effective resistivity is derived using a data-driven phenomenological model with two free parameters, which are benchmarked with PIC simulations. We compare the performance of our prescription to the resistivity model from [32]—based on a kinetic approach—and to its extension employing MHD quantities.

3.1 Kinetically motivated resistivity↩︎

[32] analyzed PIC simulations of relativistic reconnection in pair plasmas and identified the terms that dominate the nonideal electric field in the generalized Ohm’s law [14]. Their analysis was restricted to regions of electric dominance, defined as having \(E_z^2> B_x^2+B_y^2\) (which is nearly identical to the condition \(|\mathbf{E}| > |\mathbf{B}|\) in the case of zero guide field). They found that the \(z\)-component of the nonideal electric field could be written as \[\label{Eq:Selvi} E_z^* = \eta_{S23,\mathrm{kin}}J_z = \left[\frac{m}{n_t e^2}\frac{\langle u_{e z}\rangle}{\langle v_{e z}\rangle}\partial_y\langle v_{ey} \rangle\right]J_z\tag{3}\] where \(n_t\) is the total number density (including both electrons and positrons), \(\langle v_{e z}\rangle\) and \(\langle u_{e z}\rangle\) are respectively the mean electron three- and four-velocity in the \(z\)-direction,¹ and \(\langle v_{ey} \rangle\simeq v_{y}\) is the mean three-velocity along \(y\), which is roughly the same for both species (hereafter, we call \(v_{y}\) the single-fluid \(y\) velocity).

The effective resistivity proposed by [32] in 3 has a few limitations: (i) it provides a satisfactory description of the nonideal electric field only for \(B_g=0\); and (ii) it was derived considering regions of electric dominance, which are only a subset of the regions hosting nonideal fields [42], [43], where resistive effects are important. In order to derive 3 , [32] used the approximation \[\label{Eq:approx} \partial_y \left(n_e \langle v_{ey} \rangle \langle u_{ez} \rangle \right) \approx n_e \langle u_{ez} \rangle \partial_y \langle v_{ey}\rangle\tag{4}\] which is valid only in the vicinity of the center of the current sheet. In fact, as shown in 2, the effective resistivity in 3 (hereafter, \(\eta_{S23,\rm kin}\)) provides a reasonable description of the nonideal field near the center of the layer (\(|y|\omega_p/c\lesssim 1\)), where \(|\mathbf{E}| > |\mathbf{B}|\) (blue shaded area), but it significantly overestimates the ground truth (i.e., the direct measurement of \(E_z^*\) from PIC runs) farther away from the layer (\(|y|\omega_p/c\gtrsim 1\)).

For use in single-fluid MHD codes, 3 needs to be rewritten using fluid quantities. As we have already discussed above, the mean three-velocity along \(y\) is roughly the same for the two species, \(\langle v_{ey} \rangle\simeq\langle v_{y} \rangle\). The most reasonable approximation for the ratio between the mean four- and three-velocities of a given species is \(\langle u_{ez} \rangle/\langle v_{ez} \rangle\simeq \langle \gamma \rangle\), where the mean particle Lorentz factor (including both bulk and internal motions) can be derived from the \(T^{00}\) component of the stress energy tensor as \(\langle \gamma \rangle=T^{00}/n_t mc^2\). This leads to a form of 3 that can be implemented in MHD: \[\label{Eq:SelviMHD} E_z^*=\eta_{S23, \rm MHD}J_z =\left[\frac{m}{n_t e^2}\langle \gamma\rangle \;\partial_y v_y \right] J_z~.\tag{5}\] As shown in 3, 5 is an excellent approximation of the kinetic form in 3 (compare top right and bottom right panels). However, as anticipated in 2, the two forms overestimate the true resistivity (top left of 3), especially at the boundaries of the current layer. In an earlier version of [32], 3 was cast in an alternative form, approximating \[\label{eq:approx} \frac{\langle u_{ez} \rangle}{\langle v_{ez} \rangle}\simeq \frac{1}{\sqrt{1-(J_z/en_t c)^2}}~,\tag{6}\] which only holds if each species has negligible internal motions and moves in the \(z\) direction with dimensionless drift speed of \(|J_z|/e n_t c\). This was recently rewritten by [34] in the form \[\label{Eq:bugli} \eta_{B24} = \frac{1}{en_t c}\sqrt{\left(\frac{mc}{e}\partial_y v_y\right)^2+\left(\Gamma E^*_z\right)^2}~.\tag{7}\] While the approximation in 6 leading to 7 does not generally hold, as shown by the poor agreement between the top right and bottom left panels in 3, 7 appears to provide a remarkably good proxy for the ground truth (compare top left and bottom left). While 7 appears to improve upon the kinetically-motivated model by [32], it loses some of the physical motivation of 3 and 5 .

While useful, the forms of effective resistivity presented in this subsection have some undesirable properties: (i) they only apply to the case of zero guide field; (ii) they contain a spatial derivative, which makes them difficult to include in relativistic MHD codes while maintaining causality [44]; (iii) they only apply to the main layer, and not to the anti-reconnection layers in between merging plasmoids (which extend along \(y\), and for which the relevant velocity derivative is \(\partial_x v_x\)); (iv) they retain a dependence on the system geometry (e.g., via the \(z\) component \(E_z^*\)), which makes it hard to incorporate in global MHD simulations where current sheets will be curved, oscillating, and generally not aligned with the coordinate axes. In the next subsection we turn to a more agnostic approach that avoids some of these issues.

3.2 Prescriptive resistivity↩︎

To overcome the limitations of the model by [32], we propose an empirical approach. We expect that in regions of strong current—as defined below—the nonideal electric field should approach \(|{\boldsymbol{E}^*}|\rightarrow (v_{in}/c)B_0\), where \(v_{in}\) is the reconnection rate (i.e., the inflow velocity of plasma into the layer, see 4), which implies that the effective resistivity should be \[\eta_{\rm eff}\rightarrow \frac{v_{in}}{c} \frac{B_0}{|\mathbf{J}|} ~.\] A choice of \(\eta_{\rm eff}\propto |\mathbf{J}|^{-1}\) is inapplicable in regions of small electric current, where the resistivity should vanish. We therefore design a form such that \(\eta_{\rm eff}\propto |\mathbf{J}|^p\) for small current densities, where \(p>0\) is a free parameter. More precisely, this should occur where \(|\mathbf{J}|\ll e n_t c\). Adding a normalization factor \(\alpha\), this motivates choosing a form

\[\label{Eq:first95guess} \eta_{\mathrm{eff}} = \frac{ \alpha B_0 \mathbf{|J|}^p}{|\mathbf{J}|^{p+1}+\left(e n_t c \right)^{p+1}}=\frac{ \alpha B_0}{|\mathbf{J}|\left[1+\left({en_tc}/{|\mathbf{J}|}\right)^{p+1}\right]}.\tag{8}\] We will determine free parameters \(\alpha\) and \(p\) from PIC simulations. This scales as \(\eta_{\rm eff}\propto |\mathbf{J}|^p /(e n_t c)^{p+1}\) at small currents and approaches \(\eta_{\mathrm{eff}}=\alpha B_0/(2 |\mathbf{J}|)\) for \(|\mathbf{J}|\simeq e n_t c\). We therefore expect \(\alpha/2\simeq v_{in}/c\), as we indeed find below (see also 4). The condition \(|\mathbf{J}|\simeq e n_t c\) corresponds to the charge starvation regime, i.e., all charge carriers move at near the speed of light. This limit is indeed realized in the inner region of the current sheet: as the bottom panel of 2 shows, the 1D profiles of \(J_z\) and \(n_t\) have the same shape, suggesting a nearly constant drift velocity \(J_z/e n_t\simeq 0.9\,c\). In fact, if we define the drift velocity \({\boldsymbol{v}}_{dr}\equiv {\boldsymbol{J}}/e n_t\), our prescription can be written as \[\label{Eq:Guess} \eta_{\mathrm{eff}} = \frac{ \alpha B_0}{|\mathbf{J}|\left[1+\left({c}/{|{\boldsymbol{v}}_{dr}|}\right)^{p+1}\right]}~.\tag{9}\] In the inner region of the current sheet, where \(|{\boldsymbol{v}}_{dr}|\simeq c\) (green line in the bottom panel of 2), we obtain \(\eta_{\rm eff}\propto |{\boldsymbol{J}}|^{-1}\), which matches the double-peaked shape of the ground truth (i.e., \(E_z/J_z\)) in the middle panel of 2. We emphasize that the density dependence in \({\boldsymbol{v}}_{dr}\propto {\boldsymbol{J}}/n_t\) is a key ingredient of our resistivity model—in fact, the density in the middle of the sheet can be significantly larger than in the immediate upstream, see bottom panel of 2.

In Section 5 we provide an equivalent, more general version of Equation 8 suitable for implementation within resistive MHD codes.

4 Results↩︎

To determine the optimal values of \(\alpha\) and \(p\) in 9 we consider the \(z\) component of the nonideal field and define a loss, or data-fit metric \[\label{eq:loss} L(\alpha, p) = \sum_{x,y} |E_z^*-\eta_{\rm eff}(\alpha, p) J_z|^2|E_z^*|\tag{10}\] i.e., we minimize the L2 loss (the mean squared error) between \(\eta_{\mathrm{eff}} J_z\) and the measured \(E_z^*\). The loss is weighted by \(|E_z^*|\) to ensure that the large regions with negligible nonideal fields do not skew our findings. We calculate the optimal parameters \(\alpha\) and \(p\) by minimizing this loss through a simple grid search. We create a composite domain including several time snapshots of the PIC simulations. For each case with varying guide field, the snapshots (roughly 15 in each case) are equally spaced from the time when reconnection first attains a quasi-steady state up to the end of our simulations, \(\omega_p t=4500\). For each snapshot, we consider a region extending along the whole domain in \(x\) and with thickness \(64\,c/ \omega_p\) along \(y\) (sufficient to enclose the largest plasmoids), centered around the current sheet. As a representative case, the loss for \(B_g/B_0=1\) is shown in 5. The manifold shows that a valley of small loss, with values of \(L\) near the global minimum (shown by the red point), stretches across a wide range of \(p\).

In order to determine the optimal \(p\) and define a range of acceptable values we adopt the following procedure. We begin by minimizing the loss on many small, randomly selected regions (hereafter, “patches”) of the composite domain. In each small patch we find that there is a clearly preferred value of \(p\) (i.e. a sharp minimum of 10 , as opposed to the wide minimum we find on global scales) which we will use to define a range of acceptable values of \(p\). We continue adding regions until the results converge, meaning that repeatedly selecting the same number of random patches produces the same outcome, regardless of which regions are chosen. We vary the patch size depending on the guide field strength, such that the patch is twice larger than the typical extent of a region with significant nonideal fields (see bottom panels in 1). The number of patches and the patch size used in this step are indicated in 1. To ensure that the loss in a given patch is informative (which is not the case for patches with small \(E_z^*\)), we require \[\label{eq:patch} \mathrm{median}(E_{z}^*)_{\rm patch} > \mathrm{median}(E^*_{z})_{\rm global}+ \mathrm{threshold}\tag{11}\] where the median is computed in a given patch (left hand side) or over the whole composite domain (right hand side). The threshold indicated on the right hand side is the difference between the 55th and 45th percentiles of the distribution of \(E_{z}^*\) in the whole domain. We use these percentiles instead of the standard interquartile range to ensure a more robust analysis that includes a greater portion of the domain. We define \(p_{\rm best}\) as the value that minimizes the loss when considering the combined area of all patches.

We then create the distribution of values of \(p\) that minimize the loss in each patch. The difference \(\Delta p_{\rm low}\) between the 50th and 16th percentiles of this distribution gives the lower limit on allowed values of \(p\), \(p_{\text{low}}=\max[p_{\rm best}-\Delta p_{\rm low},0]\), while the difference \(\Delta p_{\rm high}\) between the 84th and 50th percentiles gives the upper limit \(p_{\text{high}}=p_{\rm best}+\Delta p_{\rm high}\). These values are indicated for each guide field in 1. We plot \(p_{\mathrm{best}}\) and \(p_{\mathrm{high}}\) as a function of guide field in 6. 3 demonstrates that our findings do not depend much on the patch size or the threshold percentiles employed in 11 . We find that the optimal value of \(p\) is robust to patch size, and varies by \(\lesssim 2\%\) when the threshold percentiles are altered. Similarly, the upper and lower limits on on \(p\) decrease by a modest amount when varying these parameters.

Within the range \([p_{\rm low},p_{\rm high}]\), we then consider 15 evenly spaced values of \(p\). For each value, we find the optimal \(\alpha\) using the loss function on the combined area of all patches. This reveals that the two parameters are correlated. We interpolate to find \(\alpha(p)\) and show the resulting linear fits in 2. The range of \(\alpha/2\) corresponding to the interval \([p_{\rm low},p_{\rm high}]\) is shown by the colored shaded bands in 4. As expected, \(\alpha/2\) matches well with the measured reconnection rate, for all guide field cases (solid lines of the same color).

We finally assess how well the reconstructed \(\eta_{\rm eff}J_z\) captures the nonideal field \(E_z^*\) obtained directly from our PIC simulations. [tab:loss] shows the L2 loss obtained for \(p=p_{\rm best}\) or \(p=p_{\rm high}\). For each of the two choices, the corresponding value of \(\alpha\) is obtained from the linear fit in 2. We find that, regardless of the weight adopted in the loss function (no weight, \(|E_z^*|\) or \((E_z^*)^2\)), the L2 loss increases by less than \(\sim 10\%\) when using \(p_{\rm high}\), as compared to choosing \(p_{\rm best}\) (and for the unweighted loss of \(B_g/B_0=0\), \(p_{\rm high}\) performs better than \(p_{\rm best}\)). We therefore regard all solutions within the range of \([p_{\rm best}, p_{\rm high}]\) as acceptable.

This is also confirmed by the 2D spatial profiles shown in 7. For all the guide fields we explore, we present the ground truth in the left column (i.e., the nonideal field \(E_z^*\) measured directly from our simulations), the reconstruction \(\eta_{\rm eff}(\alpha,p_{\rm best})J_z\) in the middle column (here, \(\alpha\) is the value corresponding to \(p_{\rm best}\) based on the linear fit in 2), and the reconstruction \(\eta_{\rm eff}(\alpha,p_{\rm high})J_z\) in the right column (here, \(\alpha\) is the value corresponding to \(p_{\rm high}\)). The plot shows that the two reconstructions are equally good for non-zero guide fields, while for \(B_g/B_0=0\) the case \(\eta_{\rm eff}(\alpha,p_{\rm high})J_z\) seems to capture better the longitudinal extent of nonideal regions. Most importantly, our prescriptive resistivity performs clearly better than the kinetically-motivated models presented in 3.

Although our prescriptive resistivity was benchmarked with the \(z\) component of the nonideal field, it can successfully model other non-trivial components that appear for non-zero guide fields. This is shown in 8. While \(E_y^*\) is negligible for all guide fields, there are distinct areas in which \(E_x^*\) is significant for non-zero guide field cases. We calculate \(\eta_{\text{eff}}J_x\) via 9 using the same \(\alpha\) and \(p\) from the analysis of the \(z\)-component described above. From the results in 8 we can conclude that the scalar resistivity in 9 provides a satisfactory description of all components of the nonideal electric field, across the whole range of guide fields that we explore.

5 Discussion↩︎

We have performed a suite of PIC simulations of relativistic pair-plasma reconnection with varying guide field strength, and we have formulated an empirical prescription for the effective resistivity \(\eta_{\rm eff}\) in 8 or equivalently 9 . Our prescription depends on two free parameters, \(\alpha\) and \(p\), which are derived directly from our PIC runs —with \(\alpha/2\) expected to be comparable to the dimensionless reconnection rate. As compared to the kinetically-motivated model proposed by [32], the form of \(\eta_{\rm eff}\) that we propose has four main advantages: it is explicitly written in single-fluid MHD quantities, does not depend on spatial derivatives, is coordinate-agnostic, and is valid for any guide field. It depends only on the electric current density and the particle number density (and the two free parameters discussed above). We have demonstrated that the scalar resistivity we propose successfully describes the spatial structure and strength of all components of the nonideal field. It thus provides a promising strategy for enhancing the reconnection rate in relativistic resistive MHD approaches.

To confirm the robustness of our findings, we demonstrate in 6 that the form in 9 (with \(\alpha\) and \(p\) determined from our reference runs) provides an excellent description of nonideal fields in independent simulations which either include synchrotron cooling or resolve the plasma skin depth with \(20\) cells (as compared to \(5\) cells for our reference runs).

We conclude with an important remark. Our prescription in 8 can be equivalently written as \[\label{eq:eta95eff95alt} \eta_{\text{eff}} = \frac{|\mathbf{E}^*|}{e n_t c}\left[\frac{\alpha B_0 - |\mathbf{E}^*|}{|\mathbf{E}^*|} \right]^{\frac{1}{p+1}}.\tag{12}\] In the limit of very high \(p\), the square bracket is elevated to a very small power, yielding a contribution of order unity. Furthermore, 5 suggests that, as long as \(p\) is large, our results do not significantly depend on its precise value. In the limit \(p\gg1\), the effective resistivity simplifies to \[\label{eq:best} \eta_{\text{eff}}\simeq \frac{|\mathbf{E}^*|}{e n_t c}\tag{13}\] which has several advantages: it is simple, coordinate-agnostic, and no longer depends on the free parameters \(\alpha\) and \(p\), i.e., it holds for any guide field strength. It retains the dependence on density which we already emphasized as being of key importance. The approximation \(p\gg1\) holds for all guide fields \(B_g/B_0\geq 0.3\), see 6. In 8, we demonstrate that solutions with \(p\gg1\) provide a satisfactory fit also for the case of zero guide field. We therefore regard 13 as the most promising form of effective resistivity to implement in resistive MHD simulations of relativistic reconnection, especially in global problems where it is non-trivial to determine the guide field strength.

We conclude with three caveats. First, our results are based on 2D simulations. While the physics of particle acceleration in relativistic reconnection is dramatically different between 2D and 3D [36], [37], [45], the nonideal physics of field dissipation—the focus of our work—is roughly the same [38], [39]. Yet, dedicated 3D simulations should be performed to confirm our findings. Second, we have employed an electron-positron composition, and future work is needed to confirm our results in the case of electron-proton and electron-positron-proton plasmas. Finally, the generalization of our prescription to the regime of trans- or non-relativistic reconnection is far from trivial. In fact, the importance of charge starvation and compressibility effects in our prescriptive model, as emphasized in 3.2, is likely to change in the case of low magnetization. There, the plasma beta becomes another important parameter governing the reconnection physics. We defer the investigation of the effective resistivity in trans- and non-relativistic reconnection (for different plasma beta) to future work.

6 Additional Validations↩︎

We validate our results on two additional sets of simulations having zero guide field: first, we increase the spatial resolution, and then we perform simulations with strong synchrotron cooling. In all the cases, we find that our prescription in 9 —using the same \(\alpha\) and \(p\) as determined in the main text, see 2—provides a successful reconstruction of the nonideal field.

We first confirm our findings with a higher resolution simulation (\(c/\omega_p\) = 20 cells) having zero guide field. The length of the domain in the \(x\) direction is \(L_x = 480\,c/\omega_p\). The results in 9 confirm the robustness of our conclusions, with \(\eta_{\rm eff}(\alpha,p_{\rm high}) J_z\) visually providing the best proxy for the nonideal field.

We also perform simulations with synchrotron cooling losses and the fiducial resolution of \(c/\omega_p=5\) cells. We quantify the cooling strength via the radiation reaction Lorentz factor \(\gamma_{\rm rad}\), also known as the classical “burnoff” limit, at which the radiation-reaction drag force balances the accelerating force of the reconnection electric field, yielding \[\gamma_{\rm rad}=\sqrt{\frac{e (v_{in}/c) B_{0}}{(4/3)\sigma_{\rm T}(B_0^2/8\pi)}}~.\] The results in 10 and 11 confirm the robustness of our conclusions, both for weak (\(\gamma_{\rm rad}=100\)) and strong (\(\gamma_{\rm rad}=25\)) cooling. In particular, our prescription \(\eta_{\rm eff}(\alpha,p_{\rm high}) J_z\) visually appears to provide the best proxy for the nonideal field. In summary, 9 —with \(\alpha\) and \(p\) determined from the fiducial simulations discussed in the main text—can be successfully applied to other runs, including the important case of strong cooling losses.

7 The Full Ohm’s Law↩︎

In the main text, we reduced the full Ohm’s law for resistive relativistic single-fluid MHD (1 ) to the simpler form in 2 . We verify in 12 that our results still hold when using the full relativistic Ohm’s law for resistive MHD, as given in 1 . Differences with respect to 7 are minor, especially in lower guide field cases. For simulations with stronger guide fields, we see that when the full Ohm’s law is used, the agreement between our model and the ground truth in plasmoid cores improves.

8 Extending the range of p for zero guide field↩︎

In 7, we have shown that higher values of \(p\) appear to reconstruct more accurately the nonideal electric fields in the case of zero guide field, despite yielding formally higher loss values. Motivated by this, we explore how the 2D spatial structure of \(\eta_{\rm eff}J_z\), with \(\eta_{\rm eff}\) in 9 , changes when using values of \(p\) higher than \(p_{\rm high}\) (for each \(p\), we use the \(\alpha\) value given by the function in 2). The results shown in 13 show that values of \(p\) greater than \(p_{\rm high}\) up to at least \(p=4\) provide an excellent reconstruction of the ground truth.

References↩︎