# Relativistic reconnection with effective resistivity: I. Dynamics and reconnection rate

M Bugli, E.F Lopresti, E Figueiredo, A Mignone, B Cerutti, G Mattia, L del Zanna, G Bodo, V Berta

► To cite this version:

M Bugli, E.F Lopresti, E Figueiredo, A Mignone, B Cerutti, et al. Relativistic reconnection with effective resistivity: I. Dynamics and reconnection rate. Astronomy & Astrophysics - A&A, 2025, 693, pp.A233. 10.1051/0004-6361/202452277. hal-04796130

HAL Id: hal-04796130

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

Submitted on 26 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Relativistic reconnection with effective resistivity

## I. Dynamics and reconnection rate

M. Bugli<sup>1,2,3,*</sup> (D), E. F. Lopresti<sup>1</sup> (D), E. Figueiredo<sup>4</sup> (D), A. Mignone<sup>1,3</sup> (D), B. Cerutti<sup>4</sup> (D), G. Mattia<sup>5,6</sup> (D), L. Del Zanna<sup>7,6,8</sup> (D), G. Bodo<sup>9</sup> (D), and V. Berta<sup>1</sup> (D)

Received 17 September 2024 / Accepted 8 December 2024

## ABSTRACT

**Context.** Relativistic magnetic reconnection is one of the most fundamental mechanisms that is considered responsible for the acceleration of relativistic particles in astrophysical jets and magnetospheres of compact objects. Understanding the properties of the dissipation of magnetic fields and the formation of non-ideal electric fields is of paramount importance to quantify the efficiency of reconnection at energizing charged particles.

**Aims.** Recent results from particle-in-cell (PIC) simulations suggest that the fundamental properties of how magnetic fields dissipate in a current sheet might be captured by an “effective resistivity” formulation, which would locally enhance the amount of magnetic energy dissipated and favor the onset of fast reconnection. Our goal is to assess this ansatz quantitatively by comparing fluid models of magnetic reconnection with a non-constant magnetic diffusivity and fully kinetic models.

**Methods.** We performed 2D resistive relativistic magnetohydrodynamic (ResRMHD) simulations of magnetic reconnection combined to PIC simulations using the same initial conditions (i.e., a Harris current sheet). We explored the impact of crucial parameters such as the plasma magnetization, its mass density, the grid resolution, and the characteristic plasma skin depth.

**Results.** Our ResRMHD models with effective resistivity are able to quantitatively reproduce the dynamics of fully kinetic models of relativistic magnetic reconnection. In particular, they lead to reconnection rates consistent with PIC simulations, whereas for constant-resistivity fluid models, the reconnection dynamics is generally ten times slower. Even at modest resolutions, adopting an effective resistivity can qualitatively capture the properties of kinetic reconnection models and produce reconnection rates compatible with collisionless models (i.e., on the order of ~$10^{-1}$).

**Key words.** acceleration of particles – magnetic reconnection – magnetohydrodynamics (MHD) – plasmas – relativistic processes – methods: numerical

## 1. Introduction

Relativistic magnetic reconnection is one of the main physical processes invoked to explain the high-energy emission produced by strongly magnetized astrophysical plasmas. Its core mechanism consists of the rearrangement of magnetic field lines, which rapidly transfers the magnetic energy stored in the plasma into heat, high-speed flows, and accelerated particles. This process is believed to power a large number of highly energetic peculiar transients, such as flares from black hole magnetospheres (Nathanail et al. 2022; Ripperda et al. 2022), pulsar magnetospheres (Lyubarskii 1996; Uzdensky & Spitkovsky 2014) and pulsar wind nebulae (Cerutti et al. 2012a; Olmi et al. 2016), fast radio bursts Mahlmann et al. (2022), Most & Philippov (2023), gamma-ray bursts Zhang & Yan (2011), McKinney & Uzdensky (2012), and flares from blazar jets (Giannios 2013; Petropoulou et al. 2016).

Particle in cell (PIC) numerical simulations are an essential numerical tool to model the fundamental physics behind the reconnection of current sheets in astrophysical plasmas (Zenitani & Hoshino 2001; Cerutti et al. 2012b; Sironi & Spitkovsky 2014). Their ability to capture the collisionless nature of astrophysical plasmas makes it possible to quantitatively model the dissipation of magnetic fields and the consequent acceleration of relativistic particles from first principles. However, their high computational cost poses a serious challenge when working to address the large separation between typical length scales of the plasma (e.g., particles gyro-radius, plasma skin depth) and the size of the astrophysical systems where reconnection takes place. On the other hand, resistive relativistic magnetohydrodynamics (ResRMHD; see, e.g., Komissarov 2007; Palenzuela et al. 2009) offers a powerful framework to study the impact of magnetic dissipation on the large-scale dynamics of accretion flows around compact objects (Qian et al. 2017; Ripperda et al. 2019a; Vos et al. 2024) and the launch of relativistic polar outflows (Mattia et al. 2023, 2024).

<footnote n="1">¹ Dipartimento di Fisica, Università di Torino, Via P. Giuria 1, I-10125 Torino, Italy</footnote>

<footnote n="2">² Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ INFN - Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy</footnote>

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

<footnote n="5">⁵ Max-Planck Institute for Astronomy (MPIA), Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="6">⁶ INFN, Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, (FI), Italy</footnote>

<footnote n="7">⁷ Dipartimento di Fisica e Astronomia, Università di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, (FI), Italy</footnote>

<footnote n="8">⁸ INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="9">⁹ INAF, Osservatorio Astrofisico di Torino, Strada Osservatorio 20, I-10025 Pino Torinese, (TO), Italy</footnote>

<footnote>* Corresponding author; matteo.bugli@unito.it</footnote>
By introducing an explicit physical magnetic diffusivity, such simulations can quantitatively address the formation of current sheets and the departure from an ideal relativistic magnetohydrodynamic (RMHD) regime (Del Zanna et al. 2016), while also decoupling the dissipation properties from the specific numerical implementation of the code adopted. Nonetheless, such resistivity is only a proxy for the actual magnetic dissipation occurring in the plasma, as its value is typically assumed to be constant across the domain and is significantly limited by the numerical diffusion imposed by the grid discretization. Moreover, reconnection rates produced by fluid models are generally one order of magnitude lower than those obtained with fully kinetic models (Sironi & Spitkovsky 2014), suggesting a fundamental shortcoming in the ability of ResRMHD simulations to capture the properties of collisionless reconnection.

One possible way to improve the description of magnetic dissipation in the ResRMHD framework is to introduce an “anomalous resistivity” in specific regions close to the current sheet, where the dissipation is significantly enhanced Watanabe & Yokoyama (2006), Zanotti & Dumbser (2011). Such an approach can significantly increase the rate at which magnetic field lines are reconnected, but it requires one to select a priori a fraction of the numerical domain where dissipation needs to be amplified. Some works have instead adopted the more practical strategy to increase resistivity in regions with high values of current density Zenitani et al. (2010), Ripperda et al. (2019b), which is more flexible and can deal with non-trivial configurations of the current sheet. This prescription can lead to faster reconnection as well, but the amplitude of the “effective” resistivity was still set by an arbitrary constant that was not constrained by a specific physical model of collisionless reconnection. Selvi et al. (2023) provided the first physically motivated model for effective resistivity informed by first-principles PIC simulations of 2D Harris sheets that is consistent with a classical collisionless inertial resistivity associated with a finite permanence of charged particles within the reconnecting region (Speiser 1970). By measuring the dominant component of the electric field within the reconnecting layer in the comoving frame and comparing it to the local current density, they constructed a closed form for an effective resistivity, $\eta_{\text{eff}}$, that depends on the particles’ mean velocity and density.

The present study is aimed at testing, for the first time, the prescription proposed by Selvi et al. (2023) by using 2D ResRMHD models of magnetic reconnection. We reformulated the expression for the effective resistivity in terms of fluid quantities available in MHD models and measured its effects on the dynamics of a standard 2D Harris current sheet, quantifying the impact of the different quantities that set its magnitude. By performing a set of PIC simulations sharing the same initial conditions as the fluid models, we quantified the accuracy of the effective resistivity prescription to capture the non-collisional properties of magnetic dissipation, obtaining a very good agreement between the two numerical frameworks.

Our study takes the following structure. In Sect. 2, we introduce the numerical methods we adopted and the derivation of the closed form of $\eta_{\text{eff}}$ we use in the rest of the paper. In Sect. 3, we present a detailed analysis of the reconnection dynamics introduced by the effective resistivity and the comparison with the constant resistivity and kinetic cases. In Sect. 4, we quantify the impact of numerical dissipation (through a convergence study) and physical parameters (e.g., magnetization, density, and plasma skin depth) in combination with PIC models. Finally, we share our conclusions in Sect. 5 and present numerical benchmarks validating our resistivity implementation in Appendix A.

We leave the analysis of the particle acceleration process and the impact of effective resistivity on their spectral distribution to an upcoming follow-up study. We use CGS Gaussian units throughout the paper.

## 2. Numerical methods

### 2.1. Resistive RMHD models

Our ResRMHD simulations were performed using the PLUTO code Mignone et al. (2007). We solved a set of conservation laws that take the form:

$$ \partial_t U = -\nabla \cdot F(U) + S(U), \quad (1) $$

with the vectors of conservative variables U and fluxes F given by

$$ U = \begin{pmatrix} D \\ m_i \\ \mathcal{E} \\ B_i \\ E_i \end{pmatrix}, \quad F_j = \begin{pmatrix} \rho u_j \\ T_{ij} \\ c^2 m_j \\ c \epsilon^{ijk} E_k \\ -c \epsilon^{ijk} B_k \end{pmatrix}. \quad (2) $$

Here, $D = \Gamma \rho$ is the mass density of the plasma, $m = \rho h \Gamma u + \frac{1}{4\pi c} E \times B$ is the total momentum density, $\mathcal{E} = \rho c^2 h \Gamma^2 - p + u_{\text{em}}$ is the total energy density, $T_{ij} = \rho h u_i u_j - \frac{1}{4\pi}(E_i E_j + B_i B_j) + p_{\text{tot}} \delta_{ij}$ is the total stress tensor, $B$ is the magnetic field, $E$ is the electric field (all measured in the laboratory frame), and $\epsilon^{ijk}$ is the 3D Levi-Civita symbol. In the previous expressions we introduced the Lorentz factor, $\Gamma$, the electromagnetic energy density, $u_{\text{em}} = \frac{1}{8\pi}(E^2 + B^2)$, the total pressure, $p_{\text{tot}} = p + u_{\text{em}}$, and the specific enthalpy, $h$, along with the so-called “primitive” variables, namely: the rest mass density, $\rho$, the fluid spatial velocity, $v$, or the related spatial part of the four-velocity, $u = \Gamma v$ (then $\Gamma = \sqrt{1 + (u/c)^2}$), and the thermal pressure, $p$.

The source term vector in Eq. (1) is instead given by

$$ S = \begin{pmatrix} 0_{\times 8} \\ -4\pi J_i \end{pmatrix}, \quad (3) $$

where the current density, $J$, is expressed in terms of primitive variables and the electromagnetic field using the relativistic Ohm’s law for a resistive plasma (Komissarov 2007):

$$ J = \rho_e v + \frac{1}{\eta} \Gamma \left[ E + \frac{1}{c} v \times B - \frac{1}{c^2} (E \cdot v) v \right], \quad (4) $$

where $\rho_e$ is the electric charge density and $\eta$ is the electric resistivity, taken here as the inverse of the fluid conductivity.

A well-known property of the evolution equation for the electric field, as first pointed out by Komissarov (2007), is that it can be potentially stiff for low values of the magnetic resistivity, since parts of its source term (i.e., those inversely proportional to $\eta$) change over characteristic timescales that may end up being much shorter than the rest of the contributions (Palenzuela et al. 2009; Bucciantini & Del Zanna 2013; Del Zanna & Bucciantini 2018). To ensure numerical stability to the time integration of Eq. (2), we used the IMEX-SSP3(3,3,2) scheme (Pareschi & Russo 2005), following a similar procedure to the one presented in Mignone et al. (2019) and Tomei et al. (2020).

The evolutionary laws described in Eqs. (1), (2), and (3) need to be coupled to the solenoidal constraint for the magnetic field and the conservation of electric charge, namely:

$$ \nabla \cdot B = 0, \quad (5) $$

$$ \nabla \cdot E = 4\pi\rho_e. \quad (6) $$
<table><caption>**Table 1.** ResRMHD and PIC models and their parameters: upstream “cold” magnetization, Alfvèn speed, resistivity, effective plasma skin depth, upstream density, and magnetic field.</caption><thead><tr><th>Model</th><th>σ<sub>0</sub></th><th>c<sub>A,0</sub></th><th>η&#773;</th><th>δ&#773;<sub>u</sub></th><th>ρ<sub>0</sub></th><th>B<sub>0</sub></th></tr></thead><tbody><tr><td>s10eEr1</td><td>10</td><td>0.945</td><td>η<sub>eff</sub></td><td>0.002</td><td>1</td><td>3.16</td></tr><tr><td>s10eEd0r1</td><td>10</td><td>0.945</td><td>η<sub>eff</sub></td><td>0.02</td><td>1</td><td>3.16</td></tr><tr><td>s10eC2r1</td><td>10</td><td>0.945</td><td>2.5 × 10<sup>-2</sup></td><td>—</td><td>1</td><td>3.16</td></tr><tr><td>s10eC3r1</td><td>10</td><td>0.945</td><td>2.5 × 10<sup>-3</sup></td><td>—</td><td>1</td><td>3.16</td></tr><tr><td>s10eC4r1</td><td>10</td><td>0.945</td><td>2.5 × 10<sup>-4</sup></td><td>—</td><td>1</td><td>3.16</td></tr><tr><td>s10eC5r1</td><td>10</td><td>0.945</td><td>2.5 × 10<sup>-5</sup></td><td>—</td><td>1</td><td>3.16</td></tr><tr><td>s10eEr01</td><td>10</td><td>0.945</td><td>η<sub>eff</sub></td><td>0.002</td><td>0.1</td><td>1</td></tr><tr><td>s10eEr10</td><td>10</td><td>0.945</td><td>η<sub>eff</sub></td><td>0.002</td><td>10</td><td>10</td></tr><tr><td>s4eEr1</td><td>4</td><td>0.887</td><td>η<sub>eff</sub></td><td>0.002</td><td>1</td><td>2</td></tr><tr><td>s1eEr1</td><td>1</td><td>0.704</td><td>η<sub>eff</sub></td><td>0.002</td><td>1</td><td>1</td></tr><tr><td>s10pic</td><td>10</td><td>0.945</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>s4pic</td><td>4</td><td>0.887</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>s1pic</td><td>1</td><td>0.704</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

While we enforced the diverge-free condition by using staggered magnetic field components via the Constrained Transport (CT) method (Londrillo & Del Zanna 2004), we used the latter equation to compute the electric charge without enforcing the electric field divergence (Bucciantini & Del Zanna 2013). All our models were performed using also WENOZ (Castro et al. 2011) reconstruction and the MHLLC Riemann solver (Mignone et al. 2019).

### 2.1.1. Effective resistivity model

In this subsection, we describe our implementation in the PLUTO code of the effective resistivity model proposed by Selvi et al. (2023), which pertains exclusively to a pair plasma without guide field. Using 2D fully kinetic relativistic PIC models of reconnecting Harris sheets, Selvi et al. (2023) measured the resistive component of the electric field produced within the current sheet and compared it to the general expression for the relativistic Ohm law Hesse & Zenitani (2007). In their model, the sheet is aligned with the x-axis, the y direction is transverse to it, and the z-axis is a symmetry axis perpendicular to the plane containing the initial magnetic field (i.e., no guide field is included). They identified the effective resistivity parameter as the ratio between the out-of-plane components of the comoving electric field and the current density and then estimated the dominant terms contributing to the electric field in the co-moving frame of the fluid (which is expected to vanish only in the limit of a perfectly conducting plasma), obtaining:

$$ \eta_{\text{eff}} \approx \frac{m}{n_t e^2} \frac{\langle u_{\alpha z} \rangle}{\langle v_{\alpha z} \rangle} \partial_y \langle v_{\alpha y} \rangle, \quad (7) $$

where m is the mass of the electron, n_t is the plasma number density, e is the electron charge, u_αz is the z-component of the four-velocity of the species α (in our case, either an electron or a positron), v_αz is the z-component of its three-velocity, and ⟨⟩ represents an average over the particle distribution. It should be noted that Eq. (7) is not covariant and its validity is limited (a priori) to the fluid comoving frame (given the expression of the Ohm law used by Selvi et al. 2023). While a covariant formulation would be much preferable, in the absence of fast bulk motions (as is the case here) there should not be any strong deviations in the effective resistivity. To express the effective resistivity

exclusively in terms of fluid quantities that are available in an RMHD simulation, we can first approximate the ratio $\langle u_{\alpha z} \rangle / \langle v_{\alpha z} \rangle$ with the electron Lorentz factor, $\Gamma_e$, and introduce the fluid rest-mass density, $\rho = mn_t$, as well as the y-component of the fluid velocity $v_y \approx \langle v_{\alpha y} \rangle$, using charge neutrality and symmetry reasons. If we go on to express $\Gamma_e$ in terms of the current density (i.e., $\rho e \langle v_{ei} \rangle / m \approx j_i$) we can use the fact that for a 2D current sheet with no guide field, the only non-vanishing component of the current density is the out-of-plane one. We thus obtain

$$ \eta_{\text{eff}} = \frac{mc}{e} \frac{\partial_y v_y}{\sqrt{(\rho\text{ec}/m)^2 - j_z^2}}. \quad (8) $$

This relation corresponds to a form of the effective resistivity for a fluid code proposed in an earlier version of Selvi et al. (2023) and can be rewritten as

$$ \eta_{\text{eff}} = \frac{E_0}{\sqrt{J_0^2 - j_z^2}}, \quad (9) $$

where we define the quantities $E_0 = (mc/e)\partial_y v_y$ and $J_0 = \rho\text{ec}/m$. By substituting $j_z = e_z/\eta_{\text{eff}}$ and solving for $\eta_{\text{eff}}$, we can finally obtain

$$ \eta_{\text{eff}} = \frac{1}{J_0} \sqrt{E_0^2 + e_z^2}, \quad (10) $$

where

$$ e = \Gamma \left[ E + \frac{1}{c} v \times B \right] \quad (11) $$

is the electric field in the fluid frame (e.g., Mignone et al. 2019).

Compared to Eq. (8), this expression for the effective resistivity has no singularities and is completely determined by fluid quantities readily available at each time step (whereas the current density would, in principle, require us to know the instantaneous value of the displacement current). More importantly, Eq. (10) has the benefit of producing a profile of resistivity that systematically peaks within the current sheet and tends to smoothly decrease once dissipation flattens the magnetic field profile (see Appendix A.2 for a 1D numerical test), thus leading to physically motivated dynamics, despite the strong approximations used in its derivation.

### 2.1.2. Numerical implementation

Equation (10) needs to be properly expressed in a dimensionless form, if we want to employ it within a resistive RMHD code. This is due to the fact that the fluid equations are scale-free in the limit of a constant scalar resistivity that does not depend on fundamental physical constants. We use $L_u$ and $\rho_u$ to indicate, respectively, the unit length and mass density that we adopted to rescale the variables appearing in Eq. (10) in terms of corresponding dimensionless quantities (indicated in the following with a bar notation), thus obtaining:

$$ \eta_{\text{eff}} = \frac{4\pi}{\omega_u \bar{\rho}} \sqrt{\langle \bar{\delta}_u \partial_{\bar{y}} \bar{v}_y \rangle^2 + \bar{e}_z^2}. \quad (12) $$

In the previous expression, we rescale the electric (and magnetic) field according to $E = \sqrt{4\pi\rho_u c^2 \bar{E}}$ and we introduce the dimensionless plasma skin-depth, $\bar{\delta}_u = c/(\omega_u L_u)$, with $\omega_u = \sqrt{\frac{4\pi e^2 \rho_u}{m^2}}$ being the plasma frequency associated with the chosen unit of
mass density. We can then obtain an expression for the dimensionless resistivity $\bar{\eta}_{\text{eff}}$ by factoring out the physical units in Ampère-Maxwell's law as

$$ \frac{\partial \bar{E}}{\partial \bar{t}} - \bar{\nabla} \times \bar{B} = -\bar{J} \equiv -\bar{\rho}_e \bar{v} - \frac{1}{\bar{\eta}_{\text{eff}}} \Gamma [\bar{E} + \bar{v} \times \bar{B} - (\bar{E} \cdot \bar{v})\bar{v}], \quad (13) $$

which leads to this final form for the effective resistivity employed in the PLUTO code

$$ \bar{\eta}_{\text{eff}} \equiv \frac{c}{4\pi L_u} \eta_{\text{eff}} = \frac{\bar{\delta}_u}{\bar{\rho}} \sqrt{(\bar{\delta}_u \partial_{\bar{y}} \bar{v})^2 + \bar{e}_z^2}. \quad (14) $$

From Eq. (14), we can see how adopting the effective resistivity prescription introduces a characteristic scale to the problem, $\bar{\delta}_u$, which is set by the specific choice of physical units for lengths and mass densities adopted by the code. This length scale represents a crucial parameter that determines the magnitude of the effective resistivity and it has a clearly defined physical meaning. For this reason, when we use Eq. (14), we set the amount of magnetic dissipation not by arbitrarily choosing a constant resistivity, $\eta_0$; rather, we properly set an equivalent plasma skin-depth length scale and a density profile for the current sheet. Since the value of $\bar{\delta}_u$ is a direct consequence of the physical units adopted by the code and the initial scale of the problem, we treat it as a fixed parameter that does not dynamically change with the local plasma skin depth (the dependence of $\eta_{\text{eff}}$ on $\rho$ already captures this effect). We compute the spatial derivative of the velocity component transverse to the current sheet via central differences and applying a minmod limiter to avoid spurious oscillations that would otherwise affect the resistivity profile. We also set a floor value of $\eta_{\text{floor}} = 10^{-6}$ in order to avoid excessively small values of resistivity that would simply hinder the stability of the code without significantly affecting the simulation's dynamics (Ripperda et al. 2019a).

## 2.2. PIC models

Our simulations of collisionless magnetic reconnection are produced with the 2D version of the relativistic electromagnetic PIC code Zeltron (Cerutti et al. 2013; Cerutti & Werner 2019), which solves the time-dependent Maxwell equations

$$ \frac{\partial \bar{E}}{\partial t} = c \nabla \times \mathbf{B} - 4\pi \mathbf{J}, \quad (15) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = -c \nabla \times \mathbf{E}, \quad (16) $$

using a second-order finite-difference time-domain method (Yee 1966). While the divergence-free condition for the magnetic field is ensured to machine precision at any given time, the Gauss equation $\nabla \cdot \mathbf{E} = 4\pi\rho_e$ is enforced by applying at each time step a correction $\delta \mathbf{E}$ to the electric field obtained by solving Poisson's equation:

$$ \nabla^2(\delta\phi) = -(4\pi\rho_e - \nabla \cdot \mathbf{E}), \quad (17) $$

where $\rho_e$ is the electric charge density and $\delta \mathbf{E} = -(\nabla \delta \phi)$. The electromagnetic field is coupled to the dynamics of the charged particles through their equation of motion

$$ m \frac{d\mathbf{u}_p}{dt} = q \left( \mathbf{E}_p + \frac{\mathbf{u}_p \times \mathbf{B}_p}{c\Gamma_p} \right), \quad (18) $$

where $m$ is the particle's mass, $\mathbf{u}_p = \Gamma_p \mathbf{v}_p$ is its 4-velocity, $\mathbf{v}_p$ is its spatial velocity, $q$ is its electric charge, $\Gamma_p = \sqrt{1 + (u_p/c)^2}$

is its Lorentz factor, and $\{E_p, B_p\}$ are the electric and magnetic fields at the particle's location. The equation of motion is integrated using the classic Boris pusher (Birdsall & Langdon 1991). Charges and currents generated by each particle are then deposited on the nodes of the grid via linear interpolation, where they are used to compute the charge and current densities for Maxwell's equations. In this work, we assume a pure electron-positron pair plasma, as in Selvi et al. (2023).

## 2.3. Initial conditions

We initialized our domain with a standard Harris current sheet of width $a$. In the following, quantities with subscript “0” refer to the upstream region of the domain (i.e., far from the current sheet). The magnetic field has non-vanishing component only along the x-axis, which is defined as

$$ B_x(y) = B_0 \tanh\left(\frac{y}{a}\right), \quad (19) $$

where $B_0$ is set by the plasma “cold” magnetization to $\sigma_0 = B_0^2/\rho_0$ and the plasma density, $\rho_0$ (here $B_0/\sqrt{4\pi} \rightarrow B$ for convenience). Once we have chosen the ratio of thermal to magnetic pressure $\beta_0 = 2p_0/B_0^2$ (where $p_0$ is the upstream plasma thermal pressure), we can compute the Alfvén speed as in Del Zanna et al. (2016)

$$ c_{A,0} = (1/\sigma_0 + 2\beta_0 + 1)^{-1/2}, \quad (20) $$

where we assume an ideal equation of state for a relativistic plasma with adiabatic index of 4/3. The thermal pressure, $p$, is then computed by requiring a constant total pressure of $p_{\text{tot}} = p + B_x^2/2$ in all the domain. This gives

$$ p(y) = \frac{1}{2} B_0^2 (\beta + 1) - \frac{1}{2} B_x^2(y), \quad (21) $$

which we then use to compute the density across the domain by assuming an initial constant temperature, $\Theta = p/\rho = 1$. This choice leads the upstream plasma density $\rho_0$ to increase by a factor of $p(0)/p_0 = (\beta + 1)/\beta \approx \beta^{-1} = 100$ within the current sheet.

The magnetic field is expressed in terms of a vector potential along the z direction

$$ A_z(x, y) = aB_0 \log(\cosh y/a) + \frac{\epsilon B_0 L_0}{2\pi} \sin\left(\frac{2\pi x}{L_0}\right) \operatorname{sech}\left(\frac{y}{a}\right), \quad (22) $$

where the second term of the z-component introduces a sinusoidal perturbation with wavelength, $L_0$, and amplitude, $\epsilon = 0.05$, that is localized in proximity of the current sheet to seed the tearing instability in the ResRMHD models and predefine the location of the main X and O points along the sheet.

Our numerical domain consists of a Cartesian box with a uniform grid spacing equal to $\Delta x \approx \delta_0/20$ (unless otherwise stated), where $\delta_0 = c/\omega_0$ is the upstream plasma skin-depth and $\omega_0$ is the corresponding plasma frequency. Our box covers the range of $[0, 4L_0]$ in the x-direction, where $L_0 = 1$ is the characteristic length of the current sheet, while we set the current sheet width $a = 0.01L_0 = \delta_0/2$. This length is 500 times the current sheet skin-depth $\delta_p$, since $\omega_p \approx 10\omega_0$. Throughout this work we consider values for the dimensionless skin-depth parameter $\bar{\delta}_u$ in Eq. (14) corresponding to either the upstream or the current sheet plasma skin-depth; namely, $\bar{\delta}_u = \{\delta_0/L_0, \delta_p/L_0\} = \{0.02, 0.002\}$. Imposing a specific value to $\bar{\delta}_u$ and assuming $L_0 = L_u$ is equivalent to setting the unit of mass density to $\rho_u = m/r_e(\bar{\delta}_u L_u)^{-2}$, with
<figure><img src="image_2.png" /><figcaption>Fig. 1. Spatial profiles of (from top to bottom) mass density, ratio of comoving electric field over current density $e_z/j_z$, and $e_z$ for the kinetic model s10pic at $t \approx 2.4$. The density is normalized by the upstream value, $\rho_0$, the electric field by the upstream magnetic field, $B_0$, and the current density by $J_0 = \frac{c}{4\pi L_0}$.</figcaption></figure>

$r_e$ being the classical electron radius. Since the rest of the MHD equations remain scale-free, this adjustment does not affect the dynamics of the system any further. The size of our numerical box in the direction transverse to the current sheet needs to be sufficiently large so as to avoid artificial high-energy cutoffs in the spectra of particles accelerated by the reconnection electric field. Werner et al. (2016) showed that when the size of the domain in the direction transverse to the current sheet is not sufficiently large, namely, when

$$L_y \lesssim 40\sigma_0 r_L, \qquad (23)$$

the spectra of particles with mean Larmor radius, $r_L$, that are accelerated by the reconnection electric field exhibit an artificial high-energy cutoff. We conservatively doubled the constraint from Eq. (23) and set $y \in [-80 \sqrt{\sigma_0}a, 80 \sqrt{\sigma_0}a]$, where we applied $\delta = \sqrt{\sigma_0}r_L$. We imposed in the x- and y-direction periodic and reflective boundary conditions, respectively.

In the PIC simulations, the plasma is generated with 32 particles per cell everywhere in the domain. The initial density gradient in the layer is achieved thanks to variable particle weights. The initial particle energy distribution follows a relativistic Maxwellian at rest in the upstream medium, while the plasma in the layer is generated according to two counter-streaming relativistic drifting Maxwellian carrying the current (one beam for each species drifting at $v_{drift} = \pm 0.6c$ along the z-direction). All kinetic models start with the same initial plasma skin depth as the ResRMHD case (i.e., $0.02L_0$ upstream, $0.002L_0$ at the center of the current sheet), the same numerical box, and a uniform grid spacing of $\Delta x = \delta_0/10$.

# 3. Effective resistivity models

## 3.1. Current sheet dynamics

We start with a brief description of the main properties of kinetic model s10pic, which represents a realization from first principles

of the reconnection of the current sheet considered in this study. In the following, we measure space and time in units of $L_0$ and $L_0/c$, respectively. The initial equilibrium is promptly destroyed by a violent fragmentation of the current sheet in multiple magnetic islands that start merging and producing ever larger plasmons, with a single large magnetic island remaining once saturation is reached. Figure 1 shows a snapshot taken at $t \approx 2.4$ to show case this well-known dynamics. A chain of four large high-density magnetic islands is interrupted by a series of small plasmons that are about to merge with the larger O-points (top panel). Regions of intense resistive electric field develop in between plasmons and at the interface of merging magnetic islands (bottom panel). The ratio between this electric field component and the current density (middle panel) represents the quantity that we refer to as the effective resistivity of the system. To avoid spurious high values in the upstream regions due to exceedingly low current densities, we do not render this ratio wherever the current density's magnitude drops below a factor $10^{-2}$ with respect to its maximum at the beginning of the simulation. We can see a clear anti-correlation with mass density in the low resistivity regions at the core of each plasmon, while high values of $e_z/j_z$ are reached in correspondence to high values of the resistive electric field. This relations are both in agreement with the dependencies included in Eq. (14) that we previously introduced.

We go on to analyze models s10eEr1 and s10eEd0r1 to provide a reference for the reconnection dynamics produced by the effective resistivity prescription (a complete list of the simulations included in this work is presented in Table 1). Here, we consider a simulation with $N_x = 4096$ points along the x-axis (which corresponds to $\Delta x = \delta_0/20$, namely, twice the resolution used in the corresponding PIC simulation). A deeper analysis of the impact of numerical dissipation is undertaken in Sect. 4.1. As a reference, for a given resolution the PIC models require about twice as much computational resources with respect to the ResRMHD models. The employment of AMR in the grid structure would have further
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Space-time diagrams of different quantities at the center of the current sheet ($y = 0$) over time for the ResRMHD model s10eEr1. From top left to bottom right: Rest-mass density $\rho$, effective resistivity $\eta_{\text{eff}}$, z-component of the electric field in the fluid frame $e_z$, and Lorentz invariant $B^2 - E^2$. The white dotted lines indicate the instant where the profiles in Fig. 3 are shown.</figcaption></figure>

reduced the computational cost of the fluid models (especially in the upstream region, where it is enough to use a coarser resolution), but we opted for a uniform grid in order to more closely match the one used by PIC simulations.

In Fig. 2, we show a space-time diagrams of several quantities at the center of the reconnecting layer (similarly to Fig. 2 of Nalewajko et al. 2015, although for ResRMHD simulations rather than PIC). At the beginning of the simulation the current sheet becomes thinner by ≈20% due to the onset of the tearing instability, with the initial perturbation inducing the formation of a small magnetic island at $x = 1$ around $t \approx 1.7$ (top-left panel of Fig. 2). During this phase, the resistivity remains relatively low (~$10^{-5}$) close to the two large X and O points forming in $x = \pm 1$ (top right panel). After this phase, the tearing instability enters its non-linear stage, with the current sheet beginning to also fragment into smaller plasmons. The resistivity increases by more than an order of magnitude in between the smaller magnetic islands (where density decreases), as does the electric field in the plasma comoving frame. At $t \approx 2.5$, the initial chain of small plasmons starts to merge towards $x = \pm 1$, and by $t \approx 4$ two large mag-

netic islands are left. Then, newly formed plasmons continue to quickly merge toward them throughout the simulation. Both the resistivity and the non-ideal electric field spike in the low-density regions of the current sheet's center at the sides of the plasmons, with a clear correlation between them (top-right and bottom-left panels of Fig. 2). These regions, which are sparsly distributed, also produce favorable conditions for the acceleration of particles, as indicated by the negative values assumed by the Lorentz invariant $B^2 - E^2$ (Sironi 2022). The O-point formed in $x = -1$ around $t \approx 4$ instead slowly approaches the lower periodic boundary of the x axis and merges across it with the original island around $t \approx 18$. This halts the further reconnection of magnetic field lines, which leads to a drop of the non-ideal electric field on the sides of the last plasmoid and the relaxation towards a global saturated state with a single large plasmoid.

Figure 3 shows the profiles of a series of quantities at $t = 6.3$ for model s10eEr1, when the two largest plasmons are steadily merging with smaller ones. The morphology of the reconnecting layer is consistent with the standard scenario produced by both kinetic and fluid models: the fragmentation of the current sheet produces a continuous coalescence of newly formed small structures into larger magnetic islands, with thin current filaments connecting them. The effective resistivity is larger around and in between the plasmons, with strong peaks emerging at the interface of merging islands and in agreement with the results obtained by model s10pic. The non-ideal electric field also reaches its maximum value in these regions, but it is otherwise stronger within the magnetic islands, rather than at the edges. These regions also host the most promising sites for the initial phases of particle acceleration, with negative values of $B^2 - E^2$ reached in small patches in the middle of highly resistive parts of the current sheet. The last panel of Fig. 3 shows how the highest temperatures are reached not only around the dense core of the two largest plasmons, but also in the rarefied regions of strong magnetic dissipation and non-ideal electric field, followed by the regions between magnetic islands where the electric field and the resistivity both peak.

Model s10eEd0r1, for which the value of $\bar{\delta}_u$ is ten times greater, shows a similar initial transient in terms of current sheet thinning and fragmentation. However, both resistivity and electric field are initially weaker than the previous case (top right and bottom left panels of Fig. 4). In the non-linear stage of the simulation, we still have two large plasmons forming at $x = \pm 1$ and merging together on the right side of the numerical box, but the process is faster and leads (by $t \approx 12L_0/c$) to a saturated state. Both resistivity and electric field still exhibit a very good correlation, but they are now much stronger and cover larger fractions of the current sheet's length. This is accompanied by a clear lack of production of short-lived small plasmons between the two main magnetic islands, where the density generally drops by more than an order of magnitude. The regions with favorable conditions for particle acceleration (bottom right panel) are now much larger and last longer in time, consistently with the aforementioned dissipation dynamics.

The spatial profiles of such quantities at $t = 4.5$ (when the activity within the current sheet is most intense) are also qualitatively different than the case using $\bar{\delta} = 0.002$ (Fig. 5). Both resistivity and non-ideal electric field reach their maximum values in the regions of low density (i.e., between the high-density merging plasmons), which is consistent with the effective resistivity being inversely proportional to $\rho$. However, rather than having thin patches in between merging magnetic islands, we now have large, rarefied cavities with turbulent edges, which disrupt the current sheet and prevent the formation of chains of small-scale
<figure><img src="image_4.png" /><figcaption>Fig. 3. Spatial profiles of (from top to bottom) the rest-mass density, $\rho$, the effective resistivity, $\eta_{eff}$, the z-component of the electric field in the fluid frame, $e_z$, the Lorentz invariant, $B^2 - E^2$, with the contours of the vector potential, $A_z$, and the temperature for the ResRMHD model s10eEr1 at $t = 6.3$.</figcaption></figure>

plasmoids. These regions also host the most promising sites for particle acceleration, with negative values of $B^2 - E^2$ reached in patches that are considerably larger than those produced in model s10eEr1. The high temperatures reached between plasmoids (last panel) confirm that such bubbles are produced by an excess of magnetic dissipation due to the much larger resistivity that comes with a larger value of $\delta$. Such an evolutionary scenario deviates significantly from the reconnection dynamics produced by the self-consistent PIC model, thereby representing a limiting case for the effective resistivity formulation. This effect suggests that the layer thickness is set by the plasma skin depth scale in the sheet rather than in the upstream medium, which makes more sense physically. Model s10eEr1 will therefore serve as our reference run for the rest of the paper.

## 3.2. Global evolution of the magnetic field

The onset of the tearing instability triggered by the magnetic seed perturbation leads to an initial fast growth of the RMS y-component of the magnetic field (top panel of Fig. 6). The initial

transient corresponding to the linear phase of the instability lasts roughly until $t \approx 2.5$, at which point the fragmentation of the current sheet into smaller plasmoids induces an even faster increase of the reconnected magnetic field. The growth rate in this stage is roughly the inverse of the ideal time scale $\tau_A = L_0/c_A$, which is even faster than the value predicted by the “ideal tearing” regime (Pucci & Velli 2014; Landi et al. 2015), namely, $\approx 0.6\tau_A^{-1}$. Such a value corresponds, in principle, to a specific choice of aspect ratio for the current sheet (i.e., $L/a = S^{1/3}$) and becomes smaller for thicker current sheets according to (Del Zanna et al. 2016; Komissarov & Phillips 2024):

$$ \gamma_{\max} \tau_A \approx 0.6S^{-1/2}(a/L)^{-3/2}. $$

However, since the system is unstable on dynamical time scales when the current sheet reaches this thickness, obtaining a smaller value of $a/L$ (hence faster reconnection) is hindered by the fragmentation of the current sheet itself. The growth rate $\gamma_{max} = 0.6\tau_A^{-1}$ can therefore be considered an upper limit to what can be produced by a reconnection event in the ResRMHD regime
<figure><img src="image_5.png" /><figcaption>Fig. 4. Same as Fig. 2, but for the ResRMHD model s10eEd0r1. The white dotted lines indicate the instant at which the profiles in Fig. 5 are shown.</figcaption></figure>

with constant resistivity. After $t \approx 7.5$, there is no significant growth of this quantity for model s10eEd0r1 as the system has reached a fully saturated state. The model with a lower value for $\bar{\delta}$ produces almost the same evolution for $\langle B_y^2 \rangle$, with only marginally slower growth after $t \approx 4$. Overall, the corresponding PIC simulation (black dotted line) produces a very similar evolution as well. The growth of $\langle B_y^2 \rangle$ in the kinetic model is fast and steady, with a rate of increase that is close to the one experienced by the ResRMHD simulations after $t \approx 2.5$; namely, this is the moment when their current sheet starts fragmenting. When the RMS field reaches the value $\sim 0.1 B_0$, the field in the PIC simulation quickly saturates to the same intensity as for the fluid models ($\approx 0.3 B_0$). The difference between the initial transients occurring in the fluid and kinetic models stems from the fact that the former needs to go through a thinning of the current via the tearing instability, before reaching the point where the first plasmoids can form. On the other hand, the current sheet in model s10pic fragments immediately, therefore producing without delays a powerful reconnection event.

As previously seen, the effective resistivity assumes a wide range of values during the simulation. In the bottom panel of Fig. 6, we show the evolution over time of the average resistivity in the current sheet (specifically, within a distance $a$ from

the center) and the maximum value it assumes within the same region for models s10eEr1 and s10eEd0r1. We report it in units of $B_0/\rho_0$, since it is a quantity appearing from Eq. (14) that provides a natural scaling unit. During the linear phase of the tearing instability the average resistivity increases immediately from the initial floor value of $\approx 3 \times 10^{-7} B_0/\rho_0$, due to the local resistive electric field building up (bottom left panels of Figs. 2 and 4). The increase ranges from a few (model s10eEr1) up to ten times (model s10eEd0r1), depending on the value of $\bar{\delta}$. Around $t \approx 2.5$, the resistivity then grows by up to three orders of magnitude due to the current sheet beginning to fragment, with a steady mean value at the center of the current sheet of roughly either $3 \times 10^{-5} B_0/\rho_0$ (s10eEr1) or $2 \times 10^{-3} B_0/\rho_0$ (s10eEd0r1). There is also a difference between the two models that is greater than the ratio between the two plasma skin-depths used by these models. This is due to the fact that the non-ideal electric field in model s10eEd0r1 is stronger than the one produced with a smaller value of $\bar{\delta}$, hence, the local dissipation rate is increased even more. After this fast increase in resistivity, both the average and maximum values do not vary much over time, with a ratio between the two rising up to $\sim 500$. However, in model s10eEr1, the maximum value of $\eta_{\text{eff}}$ right after the formation of the first plasmoid chains ($t \approx 2.5$) decreases steadily by an order of magnitude over the period required by the reconnected magnetic field to saturate, whereas such a decrease does not occur in the case of larger $\bar{\delta}$.

To quantify the efficiency of reconnection in our models, we first computed (for both the ResRMHD and PIC simulations) the reconnected magnetic flux $\Phi_{\text{rec}}$ along the current sheet ($y=0$) as

$$ \Phi_{\text{rec}} = - \int_{x_X}^{x_O} B_y(x, 0) dx = A_z|_{X} - A_z|_{O}, \quad (24) $$

where we identify the position of the X (O) point by the maximum (minimum) value of $A_z$. For both models with effective resistivity, the reconnected magnetic flux starts increasing significantly from $t \approx 2.5$ (i.e., when the current sheets start fragmenting). Model s10eEr1 produces a rather steady growth of the magnetic flux, slowly saturating towards $\Phi_{\text{rec}} \approx 0.5\Phi_0$ by $t \approx 40$, where $\Phi_0 = B_0L_y/2$ is a measure of the magnetic flux available for reconnection at the beginning of the simulation. The flux from model s10eEd0r1 increases instead at a faster pace, it slows down after $t \approx 8.5$ (due to the merger between the two last plasmoids being slower at the beginning) and then saturates at a similar value. Such deviation from the model with lower $\bar{\delta}$ reflects the faster and more violent reconnection dynamics produced by this model. Similarly to the evolution of $B_y$, the PIC simulation shows an immediate growth of the reconnected flux shortly after the beginning of the simulation, since the current sheet does not undergo a preliminary linear stage of the tearing instability. Besides this small offset, however, the evolution of the reconnected flux of the PIC model is very similar to that of the ResRMHD ones. The initial growth of $\Phi_{\text{rec}}$ falls roughly in between the two fluid models, halting briefly before the merger of the two last large magnetic islands and then saturating close to $0.5\Phi_0$ as well.

We also computed the corresponding reconnection rates, $\beta_{\text{rec}}$, as

$$ \beta_{\text{rec}} = \frac{v_{\text{rec}}}{c_{A,0}} = \frac{1}{B_0 c_{A,0}} \frac{d\Phi_{\text{rec}}}{dt}, \quad (25) $$

where we identify the inflow velocity of the magnetic field lines into the X-point with $v_{\text{rec}} = B_0^{-1} d\Phi_{\text{rec}}/dt$. Besides the delay
<figure><img src="image_6.png" /><figcaption>Fig. 5. Same as Fig. 3, but for the ResRMHD model s10eEd0r1 at t = 4.5.</figcaption></figure>

between PIC and ResRMHD models, the reconnection rates obtained with the effective resistivity formulation are remarkably similar to the kinetic result. In the s10eEd0r1 simulation, we see that $\beta_{\text{rec}}$ spikes to a maximum value of $\approx 0.2$ at the start of the current sheet fragmentation, which lasts until $t \approx 6$. Model s10eEr1 produces a qualitatively comparable scenario, with a peak around $\approx 0.1$ and a steadier mean value over time, which is consistent with the less dramatic increase in resistivity with respect to the case having $\bar{\delta} = 0.02$. Finally, PIC and ResRMHD models exhibit a very similar time variability, on both shorter and longer timescales.

### 3.3. Comparison with the uniform resistivity case

Our four models with constant and uniform resistivity profiles (s10eC2r1 to s10eC5r1) show the profound impact that the adoption of effective resistivity has on the reconnection dynamics. The magnetic dissipation $\eta_0$ ranges from $2.5 \times 10^{-2}$ to $2.5 \times 10^{-5}$, with the corresponding macroscopic Lundquist number, $S = L_c/\eta$ going from $\approx 37.8$ to $\approx 3.78 \times 10^4$. All models were produced with the same resolution as the previous

ones, namely, $N_x = 4096$. For the highest values of resistivity we tested ($\eta = 2.5 \times 10^{-2}$), the reconnected field $B_y$ experiences only a marginal initial growth, then decaying continuously after $t \sim 5$ (top panel of Fig. 6). This is due to the fact that the dissipation of the initial magnetic field across the whole domain leads to a widening of the current sheet, therefore acting against its fragmentation into plasmons. For lower values of $\eta$, the reconnected field grows faster during the transient set by the initial perturbations, followed by a sudden rise once the fastest mode selected by the tearing instability becomes sufficiently strong. This transition occurs earlier for lower resistivity, but it remains still very slow compared to the growth produced by the models s10eEr1 and s10eEd0r1. It is only the simulation with the lowest constant resistivity (i.e., s10eC5r1) that shows the start of faster growth around $t \approx 4L_0/c$ on timescales that are compatible with the onset of the ideal tearing mode (i.e., $\approx 1.66\tau_A$); therefore, they are still slower than the effective resistivity cases. This simulation unfortunately crashed due to the low value of resistivity imposed; therefore, it remains unclear to what extent such growth would continue over time.
<figure><img src="image_7.png" /><figcaption>Fig. 6. RMS value of the transverse component of the magnetic field (top) and values of the resistivity (bottom) over time for ResRMHD models s10eEr1 (blue curve), s10eEd0r1 (red), models with constant resistivity (green, yellow, light purple, and purple), and model s10pic (dotted). The maximum and average values of the resistivity are calculated over a rectangular box of sizes $L_x \times 2a$ centered on $y = 0$.</figcaption></figure>

Despite the fact that the average resistivity within the current sheet of models s10eEr1 and s10eEd0r1 falls within the range sample by our simulations with constant $\eta$ (bottom panel of Fig. 6), none of the latter models lead to reconnection dynamics that is as fast as in the PIC case. For all simulations with constant resistivity the reconnected flux $\Phi_{rec}$ grows much more slowly than the effective resistivity case (top panel of Fig. 7) and produce reconnection rates that are at least one order of magnitude smaller than PIC or effective resistivity cases (bottom panel).

The sole exception to this trend is model s10eC5r1, which starts producing a significant increase of $\Phi_{rec}$ around $t \approx 4$. This is consistent with the growth of $B_y$ for this simulation, which coincides also with a faster reconnection rate peaking at ~0.05. This is due to the fact that the resistivity is approaching the threshold for the onset of the ideal tearing instability, which in our case corresponds to $\eta_{ideal} \approx 0.945 \times 10^{-6}$ for $S = 10^6$, $\sigma_0 = 10$, $\beta_0 = 0.01$, and $a/L = 100$. However, the ideal tearing growth rate is still only 60% of the inverse Alfvén crossing time $\tau_A^{-1}$, while for models with effective resistivity the instability develops on a timescale $\tau_A$. It also remains unclear (given the limited duration of the simulation) to what extent the reconnection rate might remain significantly higher over time than in cases with higher resistivity.

Our simulations with constant resistivity allow us to estimate the numerical resistivity for a grid spacing $\Delta x \approx 0.002$, that is $\eta_{num} \lesssim 10^{-5}$. This can be seen from the fact that model s10eC5r1 shows a qualitatively different behavior than the cases with higher resistivity, namely, the onset of the ideal

<figure><img src="image_8.png" /><figcaption>Fig. 7. Reconnected magnetic flux (top) and reconnection rate (bottom) over time for the same models as those in Fig. 6.</figcaption></figure>

tearing instability. This estimate seems also to be consistent with Eq. (101) from Komissarov & Phillips (2024)

$$ \eta_{num} = A_{\eta} \frac{\Delta x}{\Delta t} \mathcal{L} \left( \frac{\Delta x}{\mathcal{L}} \right)^r \approx 10^{-5}, \qquad (26) $$

where we used a grid spacing of $\Delta x = 0.002$, characteristic length of $\mathcal{L} = 1$, time step of $\Delta t = 10^{-4}$, order $r = 2$, and normalization factor of $A_{\eta} = 0.1$ (a conservative estimate based on the normalization reported by Komissarov & Phillips 2024).

# 4. Impact of numerical and physical parameters

## 4.1. Grid resolution

We explored the impact of grid resolution on our effective resistivity models by reproducing model s10eEr1 using different numbers of grid points, $N_x = \{256, 512, 1024, 2048, 4096\}$. As we kept all the other parameters fixed, these models resolve the upstream plasma skin depth with 1.25, 2.5, 5, 10, and 20 grid points, respectively.

From the top panel of Fig. 8 we can see that already with $N_x = 2048$ we obtain fully converged results for the growth and saturation of the RMS $y$-component of the magnetic field. The development of the tearing instability before its saturation between $t \approx 5$ and $t \approx 10$ is already quantitatively well captured with $N_x = 1024$ points. Even at a very modest resolution ($N_x = 256$), the effective resistivity prescription is capable to capture the fast onset of magnetic reconnection (although the field is never able to reach the same saturation level by the end of the simulation). In contrast, models with constant resistivity approaching the fast reconnection regime require generally much higher resolutions to reach convergence (Puzzoni et al. 2021). The fact that a coarser resolution leads to lower values of
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** RMS value of the transverse component of the magnetic field (top) and values of the resistivity (bottom) over time for model s10eEr1 reproduced with different resolutions. The maximum and average values of the resistivity are calculated over a rectangular box of sizes $L_x \times 2a$ centered on $y = 0$.</figcaption></figure>

the reconnected magnetic field at saturation stems mainly from two effects. The increased numerical dissipation leads, on the one hand, to a stronger diffusion of $B_x$ at the current sheet interface with the external plasma, reducing the magnetic field reser-voir available for reconnection and opposing the thinning of the current sheet. On the other hand, a coarser resolution is also asso-ciated with a lower effective resistivity within the current sheet during its fragmentation (see bottom panel of Fig. 8 between $t \approx 2.5$ and $t \approx 5$); despite still being larger than the numerical one, it still produces a weaker reconnection event.

The reconnecting magnetic flux and the associated reconnection rate shown in Fig. 9 provide a more detailed picture. Results obtained with $N_x \gtrsim 2000$ are quantitatively in agreement in terms of the amount of magnetic field reconnected by the $t = 40L_{u}/c$ mark and the corresponding reconnection rate. However, for $N_x = 4096$, we obtained a faster and steadier reconnection, reaching the PIC saturation level by $t \approx 30$. At the lowest resolution ($N_x = 256$), the reconnection reaches a halt after the initial fragmentation and merging of the first plasmons around $t \approx 20L_{u}/c$, at which point the two last magnetic islands remain unmerged for the rest of the simulations and the reconnected magnetic flux stalls. At twice the resolution ($N_x = 512$), the increase in reconnected flux in the first half of the simulation is not significantly impacted, as shown also by the similar reconnection rate of $\approx 0.05$. However, with 512 points along the x-direction, the model is able to continue the merging of the last large magnetic islands and approach the expected saturation level. Overall, the resulting reconnection rates are qualitatively similar among models at different resolutions, with simulations

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Reconnected magnetic flux (top) and reconnection rate (bottom) over time for the same models as Fig. 8.</figcaption></figure>

running on the coarser grids showing a decrease by at most a factor $\sim 2$ in the first stages of the reconnection event.

## 4.2. Initial upstream density

The effective resistivity expressed in Eq. (14) is completely defined by the dynamic state of the plasma and the characteristic length scale, $\bar{\delta}_u$. However, the initial choice of upstream density, $\rho_0$, can, in principle, have a direct impact on the properties of magnetic dissipation, since the effective resistivity scales as $\eta_{\text{eff}} \propto \rho^{-1}$. To quantify the impact of this dependence, we reproduced model s10eEr1 with a value of $\rho_0$ decreased and increased by a factor 10 (models s10eEr01 and s10eEr10, respectively). Changing this parameter modifies not only the magnitude of the magnetic dissipation, but also the intensity of the upstream magnetic field $B_0 = \sqrt{\rho_0\sigma_0}$ (for a fixed magnetization).

The evolution of the RMS reconnected magnetic field, $B_y$ (measured in units of $B_0$), appears to be rather insensitive to the variation in $\rho_0$ during the initial exponential growth phase (top panel of Fig. 10). This is consistent with the fact that the upstream Alfvén velocity in all these models is the same, since neither $\sigma_0$ nor $\beta_0$ have changed. However, models with higher density tend to start the saturation phase earlier and lead to slightly less reconnection for the magnetic field. Similarly to model s10eC5r1, our simulation with high background density did not reach its conclusion, as a generally much lower value of the resistivity decreases the numerical stability of the code. The temporal profiles of the maximum value assumed by the effective resistivity (bottom panel of Fig. 10) show how this quantity scales indeed with $B_0/\rho_0$. The average resistivity along the current sheet shows instead some deviations from such scaling behavior, with again higher densities being associated with even
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** RMS value of the transverse component of the magnetic field (top) and values of the resistivity (bottom) over time for models s10eEr10, s10eEr1, and s10eEr01. The maximum and average values of the resistivity are calculated over a rectangular box of sizes $L_x \times 2a$ centered on $y = 0$.</figcaption></figure>

lower resistivities. A similar dependence on $\rho_0$ can be found by looking at the reconnected magnetic flux (Fig. 11), with the least and the most dense cases delivering the fastest and slowest growth of $\Phi_{\text{rec}}$, respectively. While the differences between the two less dense models (s10eEr01 and s10eEd0r1) do not exceed the level of 10% by the time the two last plasmons are fully formed ($t \approx 7$), the highest density case presents a systematically slower increase of reconnected flux, which directly impacts the corresponding reconnection rates (bottom panel), with $\beta_{\text{rec}}$ ranging from about ~0.085 (model s10eEr10) to $\approx 0.11$ (model s10eEr01).

### 4.3. Initial magnetization

We go on to analyze the role played by the upstream magnetization $\sigma_0$ in determining the reconnection dynamics by considering two lower values of initial magnetization in addition to the previously discussed benchmark, namely, $\sigma_0 \in \{1, 4, 10\}$. For these setups, we produced ResRMHD models with effective resistivity as well as the corresponding PIC simulations. Our chosen values for $\sigma_0$ imply an asymptotic magnetic field $B_0 \in \{1, 2, 3.16\}$ and an upstream Alfvén velocity of $c_{A,0} \in \{0.887, 0.704, 0.945\}$, respectively. Both the kinetic and fluid models exhibit similar trends with varying magnetization, with a slower growth of the RMS $B_y$ field for lower values of $\sigma_0$ (Fig. 12). This is consistent with the expectation that the dynamic time scale of reconnection in the linear phase of the ideal tearing regime should scale with the Alfvén crossing time, $\tau_A = L_u/c_{A,0}$ (Del Zanna et al. 2016). A lower magnetization leads also to a lower saturation of the transverse magnetic field $B_y$, an effect

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Reconnected magnetic flux (top) and reconnection rate (bottom) over time for the same models as Fig. 10.</figcaption></figure>

that is more pronounced in model s1eEr1. The evolution of the resistivity presents a similar dilation of the dynamic time scale of the reconnection event (bottom panel of Fig. 12). Both the average and maximum values of $\eta_{\text{eff}}$ are also consistent among different models once they are normalized by $B_0/\rho_0$, which is a further confirmation of the scaling behavior of $\eta_{\text{eff}}$ with the upstream magnetic field and plasma density.

The evolution over time of the reconnected magnetic flux in Fig. 13 clearly shows that our ResRMHD models have a systematic initial delay due to the initial linear phase of the tearing instability, whose duration increases for lower magnetizations. This is reflected in the current sheet fragmentation into smaller plasmons occurring at later times for smaller values of $\sigma_0$. The PIC models do not experience any such delay for decreasing magnetizations, but have a similar decrease in the growth of the reconnected flux. The corresponding reconnection rates confirm that our ResRMHD models with effective resistivity lead to reconnection events with efficiency in quantitative agreement with PIC models, corroborating the results presented in Sect. 3.1. Moreover, the reconnection rate after the initial transient does not appear to be overly dependent on the magnetization (bottom panel of Fig. 13). This is to be expected when considering that lower values of sigma lead not only to a slower increase of the reconnected flux, but also lower Alfvén speeds.

### 4.4. Reconnection rate trends

Next, we can compare the reconnection rates obtained in this work in a more quantitative way, namely, by computing a representative time average for each one of simulations we performed. To exclude the initial transient and the later stage close to the saturation, we calculated the average in the time interval where the reconnected magnetic flux is between 30% and 60% of its
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** RMS value of the transverse component of the magnetic field (top) and values of the resistivity (bottom) over time for models s10eEr1, s4eEr1, and s1eEr1. The maximum and average values of the resistivity are calculated over a rectangular box of sizes $L_x \times 2a$ centered on $y = 0$.</figcaption></figure>

maximum value. We applied this procedure to a larger set of simulation with respect to the ResRMHD and PIC models presented in the previous sections, which leads to the scatter plot of $\langle\beta_{\text{rec}}\rangle$ against the spatial resolution, shown in Fig. 14.

The dependence on the grid resolution is evident for both models s10eEr1 and s10eEd0r1, as the reconnection rate increases by up to a factor ~1.5 between the most coarse grid and the reference cases discussed in Sect. 3.1. However, the run with $N_x = 8192$ for model s10eEd0r1 suggests that there should be no further quantitative deviations for higher resolutions and that indeed our measured rates are numerically converged. All fluid models produce reconnection rates within the same order of magnitude of the kinetic simulations, with the limit case of $\bar{\delta} = 0.02$ being in quantitative agreement with the corresponding PIC realizations. We assessed the variance of the rates produced by the kinetic models by performing an additional set of simulations with the same range of values for $\sigma_0$, but with weaker initial perturbations ($\epsilon = 0.001$ rather than 0.05); these are represented by the smaller stars in Fig. 14 and deviate from the runs with stronger perturbation by 10–20% at most. The same trend with magnetization is reproduced by all PIC and ResRMHD models, namely, there is little variance seen for different values of $\sigma_0$. Simulations with $\bar{\delta} = 0.02$ and $N_x = 4096$ appear to be an exception, with the reconnection rate increasing by up to 60% with higher magnetizations. The PIC models with $\sigma_0 = 1$ also deviate from the other kinetic simulations, this time by a factor of ~2. Changing the upstream plasma density impacts the corresponding reconnection rate for both values of the plasma skin depth parameter explored. A variation over two orders of magnitude in terms of $\rho_0$ introduces a deviation in the average reconnec-

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Reconnected magnetic flux (top) and reconnection rate (bottom) over time for the same models as Fig. 12.</figcaption></figure>

tion rate of $\Delta\langle\beta_{\text{rec}}\rangle$ between ~0.02 and ~0.04. Model s10eEr10 (red square) appears to have larger $\langle\beta_{\text{rec}}\rangle$ than expected from the discussion presented in Sect. 4.2. This is likely an artifact that is due to the short duration of the simulation, which prevents a more significant decrease in the reconnection rate over time.

## 5. Conclusions

We presented a series of ResRMHD numerical simulations of magnetic reconnection that address the gap between fluid and kinetic models of astrophysical plasmas by using, for the first time, an effective resistivity formulation derived from recent first-principle PIC models. We adapted the closure obtained by Selvi et al. (2023) to express the resistivity in terms of fluid quantities and known constants, which allows us to constrain the dissipation properties of the plasma by its local state and the characteristic length scales of the reconnecting current sheet.

All our ResRMHD models with effective resistivity lead to a fast reconnection of the initial magnetic field, with a quick fragmentation of the current sheet into plasmoids and the saturation into a large magnetic island within ≈20 dynamical timescales. The evolution of the current sheet is in good agreement among the fluid and kinetic models, which confirms the capability of the effective resistivity to capture the intrinsically collisionless nature of magnetic reconnection. To capture the onset of the ideal tearing instability on dynamical timescales, constant resistivity models require very low values of $\eta$; therefore, we need sufficiently high grid resolutions to describe the smaller dissipation length scale of the simulation. On the other hand, the effective resistivity formulation allows us to retrieve fast reconnections, even at a very modest resolution. This is due to the fact that the enhanced magnetic dissipation within the current sheet induces strong reconnection without requiring the tearing instability to thin the current
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Average reconnection rate against resolution for the ResRMHD models with effective resistivity (circles and squares) and the PIC models (stars). Small and large stars refer to PIC models having initial perturbations with amplitude 0.001 and 0.05, respectively.</figcaption></figure>

sheet's width. Future works will be aimed at making a more thorough comparison between the effective resistivity models and the dynamics of relativistic reconnection at high Lundquist numbers. This will require higher resolutions and will benefit from recently developed higher order schemes for RMHD and ResRMHD models Berta et al. (2024), Mignone et al. (2024).

While both strong non-ideal electric fields and velocity gradients a priori can contribute to enhance the resistivity (see Eq. 14), the former contribution dominates in the regions with the highest dissipation. This is because it is at first order in the dimensionless plasma skin-depth parameter, $\bar{\delta}_u$, rather than second. This suggests that our effective resistivity formulation might be simplified even more to include exclusively the non-ideal electric field, especially when applied to astrophysical systems with larger size than an isolated current sheet (e.g., a relativistic jet or a black hole magnetosphere). Such formulation would be similar to the one presented in Ripperda et al. (2019b), where the resistivity was taken to be proportional to the local current density, but with the coupling constant ($\bar{\delta}$) representing the “effective” plasma skin depth that the fluid system is allowed to see.

The value of this coupling constant can have a significant impact on the resulting reconnection dynamics. Taking $\bar{\delta}$ equal to the plasma skin depth in the current sheet leads to a fragmentation into plasmons and merging phase that qualitatively resembles both PIC models and constant resistivity ResRMHD simulations. The corresponding reconnection rates are much higher than in the case with constant magnetic dissipation, but slightly lower than the kinetic counterparts. On the other hand, a plasma skin depth that is ten times greater and matches the upstream conditions of the problem reproduces the diagnostics obtained by first-principles simulations in a more quantitative

way; however, the underlying plasmoid dynamics deviates qualitatively from the standard scenario, as the increased dissipation leads to the formation of hot, rarefied cavities between large magnetic islands that prevent the production of smaller plasmons connected by thinning reconnecting regions. These differences are likely to have a significant impact on how particles are accelerated during the reconnection event, which we will assess in an upcoming work. The choice of the numerical value for the $\bar{\delta}_u$ parameter becomes less evident for problems where the current sheet is not already present from the beginning. In our reconnection models, we are, in practice, imposing a characteristic mass density unit that makes this parameter equal to the actual current sheet skin depth set by the corresponding PIC model. Since a different value of $\bar{\delta}_u$ is likely to affect the formation of current sheets themselves (e.g., during the merge of two flux tubes), setting its value a priori requires more dedicated comparative studies between PIC and ResRMHD models. Future works will have to explore the dependence of our effective resistivity model on $\bar{\delta}_u$ more thoroughly for different initial pressure-balanced equilibria (e.g., with constant density or different temperatures), as well as for merging flux tubes (Ripperda et al. 2019a).

Our simulations show that a higher background density tends to weaken the reconnection and decrease the reconnection rate, as the mean resistivity within the current sheet is lowered by a factor of a few. This is a direct consequence of the prescription we used for $\eta_{eff}$ being inversely proportional to the plasma density. We also verified that lower magnetizations produce slower reconnection events and lower reconnection rates (which is also expected in the constant resistivity case), with the mean resistivity scaling with the quantity $B_0/\rho_0$.

The fast reconnection of magnetic flux occurring in our ResRMHD simulations is reflected by high reconnection rates (between ~0.05 and ~0.2). These measurements are not only an order of magnitude larger than those produced in the constant resistivity case, but they are also in quantitative agreement with the corresponding PIC models. This result showcases the potential of the effective resistivity formalism to reproduce (within the standard ResRMHD framework) the results of first-principles kinetic models. This enables a more accurate modeling of large-scale astrophysical systems, where reconnection is expected to occur. Accretion flows on black holes would be a prime example of such an application, in particular in the so-called MAD state Tchekhovskoy et al. (2011). Since the amount of dissipation regulates the magnetic flux available at the event horizon, the inclusion of an effective resistivity might change the properties of the saturated state. Moreover, a faster reconnection dynamics could profoundly impact the formation of possible particle acceleration sites in relativistic jets launched by newly formed compact objects (Mattia et al. 2023, 2024). However, it should be noted that the prospects of resolving the actual skin depth of current sheets in large-scale astrophysical models are quite poor, given the large scale separation involved. The application of our formulation to models of accreting compact objects or relativistic jets will thus require a compromise, namely, adopting a much larger “effective” skin depth; otherwise, the low resulting resistivity would simply have no dynamical effects against the grid numerical dissipation (in addition to it leading to numerical instabilities). This approach would be keen to the one adopted successfully by large scale GR-PIC models of relativistic magnetospheres (Crinquand et al. 2021, 2022; El Mellah et al. 2022), which settle for a similar trade-off for the plasma skin depth. Given this similarity, a closer comparison between such models and GRMHD simulations with effective resistivity might provide
important insights on the dynamics of magnetic dissipation in a regime of strong gravity.

Although our results show a remarkable agreement between fluid and kinetic models, it should be noted that the derivation of Eq. (14) relies on some strong assumptions. In particular, we introduced the current density (and, therefore, the electric field) in the expression for $\eta_{\text{eff}}$ to more directly tie the enhanced dissipation to the reconnection layer; however, our derivation should be more rigorously verified to fully characterize the impact of non-collisional effects on magnetic dissipation within the ResRMHD framework. Moreover, the “kinetic” formulation of the effective resistivity presented by Selvi et al. (2023) has been obtained in the absence of a guide field and for an electron-positron plasma; meanwhile, it remains unclear how the scenario would change with different magnetic field configurations or in the presence of ions. Another limitation comes from the fact that the effective resistivity measurements were performed specifically in the X-points of the reconnecting current sheet. Although they are undoubtedly the most relevant locations to study, further investigations with fully kinetic models of relativistic reconnection are required to assess the general character of the proposed formulation.

Furthermore, our 2D results capture the main properties of how magnetic field is dissipated and particles are accelerated during a reconnection event, but they should be further tested in 3D to verify the extent of the validity of the prescription we used. Finally, it should be noted that Eqs. (7) and (14) assume the current sheet to be locally aligned with the x-axis, as the only velocity gradients appearing are those along the transverse direction, $y$. This means that our prescription is likely to underestimate the contribution to the effective resistivity of the fluid velocity gradients, since the non-linear stage of the tearing instability leads to the formation of secondary current sheets perpendicular to the x-axis between merging plasmons. This aspect might be relevant when quantifying the efficiency of particle acceleration when using the effective resistivity formulation. A more consistent approach would require the identification of the local orientation of the current sheet to correctly include the velocity gradients that are perpendicular to that direction. We leave such an improvement to a follow-up study.

**Acknowledgements.** We thank the anonymous referee for improving this work with useful comments and suggestions. M.B. would like to thank Bart Ripperda, Sasha Philippov, Niccolò Bucciantini, and Emanuele Papini for useful discussions. This project has received funding from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101064953 (GR-PLUTO). B.C. and E.F. are supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 863412). L.D.Z. acknowledges support from the ICSC-Centro Nazionale di Ricerca in High-Performance Computing; Big Data and Quantum Computing, funded by European Union - NextGenerationEU. This work has received funding from the European High Performance Computing Joint Undertaking (JU) and Belgium, Czech Republic, France, Germany, Greece, Italy, Norway, and Spain under grant agreement No. 101093441. This paper is also supported by the Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations, National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN_00000013 "Italian Research Center on High-Performance Computing, Big Data and Quantum Computing" funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di "campioni nazionali di R&S (M4C2-19)" - Next Generation EU (NGEU). Computational resources on the OCCAM cluster (https://c3s.unito.it/index.php/super-computer) were provided by the Centro di Competenza sul Calcolo Scientifico (C3S) of the University of Torino. We also acknowledge CINECA for the availability of high performance computing resources and support on the LEONARDO supercomputer (https://leonardo-supercomputer.cineca.eu/) through a CINECA-INFN agreement, providing the allocations INF23_teongrav and INF24_teongrav. Computing resources for performing the PIC simulations were provided by TGCC under the allocation A0150407669 made by GENCI. The figures presented

in this paper were generated using the PyPLUTO (https://github.com/GiMattia/PyPLUTO) package (Mattia et al., in prep.) while all the line colors have been checked through a color blindness simulator (https://www.color-blindness.com/coblis-color-blindness-simulator/).

## References

Berta, V., Mignone, A., Bugli, M., & Mattia, G. 2024, J. Comput. Phys., 499, 112701

Birdsall, C. K., & Langdon, A. B. 1991, Plasma Physics via Computer Simulation (New York: McGraw-Hill)

Bucciantini, N., & Del Zanna, L. 2013, MNRAS, 428, 71

Castro, M., Costa, B., & Don, W. S. 2011, J. Comput. Phys., 230, 1766

Cerutti, B., & Werner, G. 2019, Astrophysics Source Code Library [record ascl:1911.012]

Cerutti, B., Uzdensky, D. A., & Begelman, M. C. 2012a, ApJ, 746, 148

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

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

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

Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., & Philippov, A. 2022, Phys. Rev. Lett., 129, 205101

Del Zanna, L., & Bucciantini, N. 2018, MNRAS, 479, 657

Del Zanna, L., Papini, E., Landi, S., Bugli, M., & Bucciantini, N. 2016, MNRAS, 460, 3753

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

Giannios, D. 2013, MNRAS, 431, 355

Hesse, M., & Zenitani, S. 2007, Phys. Plasmas, 14, 112102

Komissarov, S. 2007, MNRAS, 382, 995

Komissarov, S., & Phillips, D. 2024, ArXiv e-prints [arXiv:2409.03637]

Landi, S., Del Zanna, L., Papini, E., Pucci, F., & Velli, M. 2015, ApJ, 806, 131

Londrillo, P., & del Zanna, L. 2004, J. Comput. Phys., 195, 17

Lyubarskii, Y. E. 1996, A&A, 311, 172

Mahlmann, J. F., Philippov, A. A., Levinson, A., Spitkovsky, A., & Hakobyan, H. 2022, ApJ, 932, L20

Mattia, G., Del Zanna, L., Bugli, M., et al. 2023, A&A, 679, A49

Mattia, G., Del Zanna, L., Pavan, A., & Ciolfi, R. 2024, A&A, 691, A105

McKinney, J. C., & Uzdensky, D. A. 2012, MNRAS, 419, 573

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Mignone, A., Mattia, G., Bodo, G., & Del Zanna, L. 2019, MNRAS, 486, 4252

Mignone, A., Berta, V., Rossazza, M., et al. 2024, MNRAS, 533, 1670

Most, E. R., & Philippov, A. A. 2023, Phys. Rev. Lett., 130, 245201

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

Nathanail, A., Mpisketzis, V., Porth, O., Fromm, C. M., & Rezzolla, L. 2022, MNRAS, 513, 4267

Olmi, B., Del Zanna, L., Amato, E., Bucciantini, N., & Mignone, A. 2016, J. Plasma Phys., 82, 635820601

Palenzuela, C., Lehner, L., Reula, O., & Rezzolla, L. 2009, MNRAS, 394, 1727

Pareschi, L., & Russo, G. 2005, J. Sci. Comput., 25, 129

Petropoulou, M., Giannios, D., & Sironi, L. 2016, MNRAS, 462, 3325

Pucci, F., & Velli, M. 2014, ApJ, 780, L19

Puzzoni, E., Mignone, A., & Bodo, G. 2021, MNRAS, 508, 2771

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

Ripperda, B., Bacchini, F., Porth, O., et al. 2019a, ApJS, 244, 10

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

Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJ, 924, L32

Selvi, S., Porth, O., Ripperda, B., et al. 2023, ApJ, 950, 169

Sironi, L. 2022, Phys. Rev. Lett., 128, 145102

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

Speiser, T. W. 1970, Planet. Space Sci., 18, 613

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

Tomei, N., Del Zanna, L., Bugli, M., & Bucciantini, N. 2020, MNRAS, 491, 2346

Uzdensky, D. A., & Spitkovsky, A. 2014, ApJ, 780, 3

Vos, J. T., Olivares, H., Cerutti, B., & Mościbrodzka, M. 2024, MNRAS, 531, 1554

Watanabe, N., & Yokoyama, T. 2006, ApJ, 647, L123

Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8

Yee, K. 1966, IEEE Trans. Antennas Propag., 14, 302

Zanotti, O., & Dumbser, M. 2011, MNRAS, 418, 1004

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

Zenitani, S., Hesse, M., & Klimas, A. 2010, ApJ, 716, L214

Zhang, B., & Yan, H. 2011, ApJ, 726, 90
# Appendix A: Numerical benchmarks

Here, we present a series of tests designed to validate our implementation of the effective resistivity within the ResRMHD module of the PLUTO code.

## A.1. Uniform stationary resistive plasma

Using our prescription for the effective resistivity in Eq. 10 we can express Ampère's law as

$$ \partial_t E - \nabla \times B = -\Gamma \frac{|J_0|}{\sqrt{E_0^2 + e_z^2}} [E + v \times B - (E \cdot v)v] - qv. \quad (A.1) $$

If we now assume to have a stationary, unmagnetized, and uniform plasma with

$$ v = 0; \quad B = 0; \quad P = 1; \quad \rho = 1, $$

we can reduce Eq. A.1 to

$$ \frac{\partial E}{\partial t} = -|J_0| \frac{E}{|E_z|}. \quad (A.2) $$

This equation requires only the integration of the z component of the electric field, since we can take the other two as proportional to $E_z$ (i.e., $E_{x,y} = \alpha_{x,y}E_z$) and obtain the exact same differential equation for the electric field, namely,

$$ \frac{dE_z}{dt} = -|J_0| \operatorname{sign}(E_z). \quad (A.3) $$

The proportionality factors for $E_{x,y}$ are simply retrieved from the initial condition by requiring $E_{x,y}(t)/E_{x,y}(0) = E_z(t)/E_z(0)$. By integrating Eq. A.3 in time, we can finally obtain the analytical solution

$$ E_z(t) = \begin{cases} E_z(0) - |J_0| \operatorname{sign}(E_z)t & \text{for } t < \tau, \\ 0 & \text{for } t \ge \tau, \end{cases} $$

where $\tau = E_z(0)/|J_0|$ is the timescale required for the electric field to be completely dissipated. We note that such solution differs substantially from the case with constant resistivity $\eta_0$, which would instead lead to an exponential decay of the electric field, according to $\partial_t E = -\eta_0^{-1} E$.

The left panel of Fig. A.1 shows the results of a 1d simulation of such configuration using our effective resistivity prescription and an initial uniform electric field given by $E = (3, 2, -1)$ and setting $|J_0| = 1$. The field correctly decays linearly with time, completely vanishing at $t = \tau = 1$. The effective resistivity starts with a value equal to $E_z/|J_0| = 1$ and decays proportionally with the electric field, reaching the minimum floor value $\eta_{\text{floor}} = 10^{-6}$ at $t = 1$ (right panel of Fig. A.1).

We further tested our implementation of the effective resistivity by extending the previous test to non-vanishing velocity gradients, namely, $\partial_y v_y \neq 0$. In this case, Ampère's law takes the form of

$$ \frac{\partial E}{\partial t} = -|J_0| \frac{E}{\sqrt{E_0^2 + E_z^2}}, \quad (A.4) $$

which can be reduced to a single equation for the z component of the electric field as it was for Eq. A.1. Equation A.4 can be then integrated to obtain

$$ \sqrt{E_0^2 + E_z^2}(t) + \frac{1}{2} E_0 \ln \left| \frac{\sqrt{E_0^2 + E_z^2}(t) - E_0}{\sqrt{E_0^2 + E_z^2}(t) + E_0} \right| = -|J_0|t + K, \quad (A.5) $$

<figure><img src="image_16.png" /><figcaption>Fig. A.1. Electric field decay (left panel) and effective resistivity (right) over time for the uniform stationary resistive test with vanishing velocity gradients. The analytical solution is represented with the blue dashed line.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. A.2. Same as Fig. A.1 but for the non-vanishing velocity gradients case.</figcaption></figure>

which can be solved numerically once the integration constant K is set by the initial conditions.

The left panel of Fig. A.2 shows both the field described by Eq. A.5 and the solution obtained by the PLUTO code with $J_0 = 1$, $E_0 = 0.5$, and the same initial electric field as the previous test. Also in this case the code's calculation agrees with the analytic solution, which leads to the complete dissipation of the electric field by $t \approx 4$. However, its decay is no longer linear but rather exponential-like, as the non-vanishing $E_0$ introduces a minimum constant resistivity acting during the whole evolution of the system. This can be seen in the evolution over time of the resistivity (right panel of Fig. A.2), which starts from an initial value $\eta(t=0) = \sqrt{E_0^2 + E_z^2/|J_0|} \approx 1.12$ and then decays down to the terminal value $E_0/|J_0| = 0.5$.

## A.2. 1D self-similar current sheet

As second test for our implementation of the effective resistivity, we considered the self-similar dissipation of a 1D current sheet (Komissarov 2007). While this problem has become a standard numerical test for resistive RMHD codes, we use it here as a benchmark to characterize the action of the effective resistivity in a more realistic and physically motivated setup.

In the limit of infinite resistivity and thermal pressure only the magnetic field evolves, with the induction equation being the
only one that needs to be solved. In this resistive regime, the magnetic field evolves in time according to

$$B_x(y, t) = B_{0\text{erf}}\left(\frac{y}{2\sqrt{\eta t}}\right). \qquad (A.6)$$

We start our simulations at the instant $t = 1$ imposing $\rho = 1$, $p = 50$, $B_0 = 1$, and $\mathbf{B} = (B_x(y, t = 1), 0, 0)$. Along with the test using the effective resistivity prescription, we also run a benchmark with constant resistivity $\eta_0 = 0.1$ which will serve as reference. The same value is used in Eq. A.6 to set identical initial conditions for the effective resistivity case. The 1D computational domain covers with $N_y = 4096$ grid points a transversal slice of a current sheet with $y \in [-20, 20]$. As for the previous tests, we used RK3 timestepping, MP5 reconstruction, and we set $|J_0| = 1$.

The top panel of Fig. A.3 shows the evolution of the current sheet when using a constant value for the resistivity $\eta = 0.1$, which reproduces the analytical profile described by Eq. A.6 (dotted curves). When we use the effective resistivity prescription (mid panel), there is a localized dissipation of the current sheet closer to its center, with a continuous expanding of the region affected by the magnetic dissipation. The magnetic field relaxes then towards a profile that increases linearly with the domain coordinate $x$. This behaviour is consistent with the actual evolution of $\eta$ (bottom panel), which starts from its imposed floor value $\eta_{\text{floor}} = 10^{-6}$. It then increases to values larger than 0.1 by $t = 1.5$ but only within $\Delta x \approx 1$ from the current sheet's center, whereas it remains weaker by at least two orders of magnitude for further distances. As time progresses, the resistivity profile becomes flatter over a larger region, while also decreasing in amplitude and approaching a mean value within the current sheet of $\sim 0.01$ by $t \approx 3$.

<figure><img src="image_18.png" /><figcaption>Fig. A.3. 1D self-similar current sheet. Top panel: Spatial profiles of $B_y$ for a constant resistivity at different times. The analytical solution given by Eq. A.6 is indicated with black dotted lines. Mid panel: Same as the top panel, but using an effective resistivity. Bottom panel: Spatial profiles of effective resistivity at different times.</figcaption></figure>