Issue 
A&A
Volume 691, November 2024



Article Number  A105  
Number of page(s)  17  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202451528  
Published online  31 October 2024 
Magnetic dissipation in short gammarayburst jets
I. Resistive relativistic MHD evolution in a model environment
^{1}
INFN, Sezione di Firenze, Via G. Sansone 1, I50019 Sesto Fiorentino, FI, Italy
^{2}
Dipartimento di Fisica e Astronomia, Università di Firenze, Via G. Sansone 1, 50019 Sesto Fiorentino, FI, Italy
^{3}
INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
^{4}
INAF, Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
^{5}
INFN, Sezione di Padova, Via F. Marzolo 8, 35131 Padova, Italy
^{⋆} Corresponding author; mattia@fi.infn.it
Received:
16
July
2024
Accepted:
15
September
2024
Aims. Short gammaray bursts originate when relativistic jets emerge from the remnants of binary neutron star (BNS) mergers, as observed in the first multimessenger event GW170817–GRB 170817A, which coincided with a gravitational wave signal. Both the jet and the remnant are believed to be magnetized, and the presence of magnetic fields is known to influence the jet propagation across the surrounding postmerger environment. In the magnetic interplay between the jet and the environment itself, effects due to a finite plasma conductivity may be important, especially in the first phases of jet propagation. We aim to investigate these effects, from jet launching to its final breakout from the postmerger environment.
Methods. We used the PLUTO numerical code to perform 2D axisymmetric and full 3D resistive relativistic magnetohydrodynamic (MHD) simulations, employing spherical coordinates with spatial radial stretching. We considered and compared different models for physical resistivity, which must be small but still dominating over the intrinsic numerical dissipation (which yields unwanted smearing of structures in any ideal MHD code). Stiff terms in the current density are treated with IMplicitEXplicit Runge Kutta algorithms for timestepping. A Syngelike gas (Taub equation of state) is also considered. All simulations were performed using an axisymmetric analytical model for both the jet propagation environment and the jet injection; we leave the case of jet propagation in a realistic environment (i.e., imported from an actual BNS merger simulation) to a future study.
Results. As expected, no qualitative differences are detected due to the effect of a finite conductivity, but significant quantitative differences in the jet structure and induced turbulence are clearly seen in 2D axisymmetric simulations, and we also compare different resistivity models. We see the formation of regions with a resistive electric field parallel to the magnetic field, and nonthermal particle acceleration may be enhanced there. The level of dissipated Ohmic power is also dependent on the various recipes for resistivity. Most of the differences arise before the breakout from the inner environment, whereas once the jet enters the external weakly magnetized environment region, these differences are preserved during further propagation despite the lower grid refinement. Finally, we show and discuss the 3D evolution of the jet within the same environment in order to highlight the emergence of nonaxisymmetric features.
Key words: diffusion / magnetic fields / magnetohydrodynamics (MHD) / relativistic processes / gammaray burst: general / stars: jets
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The gravitational wave (GW) detection of the binary neutron star (BNS) merger event GW170817, followed 1.74 s later by the detection of the short gammaray burst (sGRB) GRB 170817A, with subsequent followup electromagnetic observations across the whole spectrum, was the first multimessenger event involving a GW source (Abbott et al. 2017a,b,c). This breakthrough event confirmed in particular that BNS mergers can produce a powerful relativistic jet associated with a sGRB (Mooley et al. 2018; Ghirlanda et al. 2019). However, various aspects of the postmerger dynamics remain elusive, including the precise nature of the merger remnant acting as the central engine for launching the jet. In general, prompt black hole (BH) formation will occur for highmass BNS mergers, while lowermass systems may lead to hypermassive neutron stars (HMNSs) in a temporary equilibrium before collapsing to a BH, or even to stable NSs. In the GW170817 case, models tend to favor the formation of a remnant surviving the collapse for about 0.5−1 seconds (e.g., Margalit & Metzger 2017), leaving the question of whether the jet was launched before or after BH formation. In the former case, the jet could be powered by a central magnetarlike object (e.g., Ciolfi et al. 2017, 2019) in which the field has been enormously amplified by shear during the merging process itself (e.g., Palenzuela et al. 2022 and references therein) via magnetorotational instability (e.g., Mösta et al. 2020), or via internal meanfield dynamo action (e.g., Franceschetti & Del Zanna 2020; Del Zanna et al. 2022; Kiuchi et al. 2024). However, this scenario is a matter of debate, and serious difficulties have been pointed out in the last few years (Ciolfi 2020). The alternative BH central engine scenario seems to find stronger support in numerical relativity simulations, where the possibility of forming at least a magnetically dominated and mildly relativistic incipient jet has been demonstrated (e.g., Ruiz et al. 2016).
GRB 170817A was also a peculiar sGRB, the first observed 15−30° away from the jetpropagation axis and well outside the jet core opening angle (of only a few degrees; Mooley et al. 2018; Ghirlanda et al. 2019). This significantly boosted investigations into how incipient jets propagate through postmerger environments and acquire their specific angular structure (e.g., Lazzati et al. 2018, 2021; Xie et al. 2018; Geng et al. 2019; Kathirgamaraju et al. 2019; Gottlieb et al. 2020, 2021, 2022; Nathanail et al. 2020, 2021; Pavan et al. 2021, 2023; MurguiaBerthier et al. 2021; Urrutia et al. 2021, 2023; Hamidani & Ioka 2023).
A “realistic” environment (i.e., directly imported from the outcome of a BNS merger simulation) has recently been shown to be a necessary ingredient in accurately modeling the above propagation and the properties of the emerging jet (Pavan et al. 2021, 2023; Lazzati et al. 2021). Moreover, magnetic fields in the environment need to be included (along with those within the jet), as they can have a significant impact on the postbreakout angular structure and the final energetics of the jetcocoon system (Pavan et al. 2023; GarcíaGarcía et al. 2023). In this case, however, the complex interaction between the magnetic field of the jet and that of the surrounding environment may suffer from excessive smoothing due to numerical viscosity (or more properly, resistivity) in magnetohydrodynamic (MHD) simulations, especially in three dimensions, where resolution limits must be imposed.
Such “spurious” dissipation, originated by truncation errors and highly susceptible to the grid resolution and the numerical algorithms employed, can strongly impact our understanding of the evolution of astrophysical environments (see, e.g., Puzzoni et al. 2021, 2022). On the other hand, numerical simulations of relativistic plasmas encompassing a physical resistivity model that dominates over the numerical diffusion (see, e.g., Palenzuela 2013; Dionysopoulou et al. 2015; Qian et al. 2018; IndaKoide et al. 2019; Vourellis et al. 2019) demonstrated the importance of a dissipative model for the simulation outcome and how a physical resistivity is crucial not only for the reproducibility of the results (which, for instance, are not affected by the algorithms employed) but also for the consistency of the different physical processes taking place.
Therefore, the goal of the present series of two papers is to investigate, for the first time, the magnetic dissipation and in general the effect of a finite plasma conductivity in twodimensional (2D) axisymmetric and threedimensional (3D) relativistic MHD (RMHD) simulations of sGRB jets propagating through a BNS merger environment, in a similar fashion to our previous work on generic astrophysical jets in a uniform medium (Mattia et al. 2023). The main result of the latter work is the finding that the fine current sheet structures and the level of turbulence are affected by the value of the plasma conductivity, for which both uniform and variable models (based either on a tracer or on the local Alfvén speed) were tested. Regions with a nonideal electric field parallel to the local magnetic field form, and plasmoids inside current sheets are visible if the resistivity is low and the resolution is sufficiently high. We also computed the magnetic dissipated power.
In the present paper (Paper I), the BNS merger environment is provided by an analytical axisymmetric model describing a magnetized wind in radial expansion –with a homologous velocity profile– surrounded by a static external atmosphere (recipes in Pavan et al. 2021). The magnetic field of the wind, in particular, is provided by a dipole with an initial strength of ∼4 × 10^{13} G superimposed on the hydro structure. Into this wind, an incipient magnetized and uniformly rotating sGRB jet is continuously injected from the inner boundary at a given radius, covering an opening angle of 10°, with an initial Lorentz factor of 3 and a total luminosity of about 10^{52} erg/s. The evolution is followed for 0.5 s after launching, when the jet has eventually broken out from the expanding environment, up to a distance of about 10^{5} km. We present a detailed comparison between 2D axisymmetric simulations and 3D ones, assuming the same setup for both the environment and the incipient jet, while we leave the presentation and discussion of the fully 3D case – with a realistic (nonaxisymmetric) environment imported following the guidelines in Pavan et al. (2023) – to the forthcoming Paper II.
The resistive relativistic MHD (RRMHD) simulations were performed using spherical coordinates with the PLUTO code (Mignone et al. 2007). We employed stateoftheart methods for resistivity in shockcapturing evolution schemes, assuming a Syngelike equation of state, which is more appropriate than the usual ideal gas condition for a relativistic plasma (for references and details see Mattia et al. 2023).
The present Paper I is structured as follows: In Sect. 2, we describe the system of equations solved, the numerical methods employed, and the setup for the initial and boundary conditions; Sect. 3 is devoted to the presentation and discussion of the results of axisymmetric 2.5D simulations; in Sect. 4 we compare those results against full 3D runs using the same axisymmetric initial setup. Conclusions are drawn in Sect. 5. Appendix A shows a comparison with an ideal run.
2. Equations and numerical setup
In this section, we describe the equations adopted (i.e., the equations of RRMHD), as well as the numerical setup employed throughout the entire paper.
2.1. Equations of RRMHD
The system of RRMHD equations consists of a set of hyperbolic partial differential equations with nonzero source terms (Komissarov 2007; Del Zanna et al. 2016; Mignone et al. 2019)
$$\begin{array}{c}\hfill \frac{\partial U}{\partial t}=\mathrm{\nabla}\xb7\mathsf{F}(U)+{S\phantom{\rule{0.166667em}{0ex}}}_{e}(U)+{S\phantom{\rule{0.166667em}{0ex}}}_{i}(U),\end{array}$$(1)
where U is the vector of conserved variables, F represents the fluxes and S_{e} and S_{i} represent, respectively, the standard (explicit) and stiff (implicit) source terms. By splitting the equations for the different conserved variables we get
$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\partial D}{\partial t}+\mathrm{\nabla}\xb7(D\mathit{v})=0,\hfill \\ \hfill & \frac{\partial \mathit{m}}{\partial t}+\mathrm{\nabla}\xb7[\rho h{\gamma}^{2}\mathit{v}\mathit{v}\mathit{E}\mathit{E}\mathit{B}\mathit{B}+(p+\frac{{E}^{2}+{B}^{2}}{2})\phantom{\rule{0.166667em}{0ex}}\mathsf{I}]=\mathit{f},\hfill \\ \hfill & \frac{\partial \mathcal{E}}{\partial t}+\mathrm{\nabla}\xb7\mathit{m}={\mathcal{P}}_{\mathit{f}},\hfill \\ \hfill & \frac{\partial \mathit{B}}{\partial t}+\mathrm{\nabla}\times \mathit{E}=0,\hfill \\ \hfill & \frac{\partial \mathit{E}}{\partial t}\mathrm{\nabla}\times \mathit{B}=\mathit{J},\hfill \\ \hfill & \mathrm{\nabla}\xb7\mathit{B}=0,\hfill \\ \hfill & \mathrm{\nabla}\xb7\mathit{E}=q,\hfill \end{array}\end{array}$$(2)
where we have set c = 1 and we have absorbed a factor $1/\sqrt{4\pi}$ into the definitions of the electromagnetic fields.
The conserved fluid variables (D, m, ℰ) are, respectively, the laboratory frame density, momentum density, total energy density, B, E and q represent the magnetic and electric fields and the electric charge density. The primitive fluid variables (ρ, v, p) are, respectively, the restmass density, the velocity, and the pressure of the fluid, h is the relativistic enthalpy function and γ is the Lorentz factor of the flow. Here f is a generic external force per unit volume, and 𝒫_{f} the corresponding power per unit volume (see the PLUTO user guide^{1} for a general description and Sect. 2.6 for the particular implementation needed for the present problem).
The conserved fluid variables to be evolved in time are defined in terms of the primitive variables as
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D& =\rho \gamma ,\hfill \\ \hfill \mathit{m}& =\rho h{\gamma}^{2}\mathit{v}+\mathit{E}\times \mathit{B},\hfill \\ \hfill \mathcal{E}& =\rho h{\gamma}^{2}p+\frac{{E}^{2}+{B}^{2}}{2},\hfill \end{array}\end{array}$$(3)
while, due to the high nonlinearity of the RRMHD equation, the conversion from conserved to primitive variables must be performed numerically (Mignone et al. 2019).
The electric current density J appearing in the source terms of the RRMHD equations is defined in terms of the electromagnetic fields as
$$\begin{array}{c}\hfill \mathit{J}=q\mathit{v}+{\eta}^{1}\gamma [\mathit{E}+\mathit{v}\times \mathit{B}(\mathit{E}\xb7\mathit{v})\mathit{v}].\end{array}$$(4)
This expression comes from the simplest possible form of Ohm’s law written for the electric field and current density in the frame comoving with the fluid, that is
$$\begin{array}{c}\hfill {e}^{\mu}=\eta {j}^{\mu},\end{array}$$(5)
where η is the resistivity coefficient, supposed to be a scalar for the sake of simplicity, and where the comoving electromagnetic fields have been split as
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {e}^{\mu}& =\gamma (\mathit{v}\xb7\mathit{E},\phantom{\rule{0.166667em}{0ex}}\mathit{E}+\mathit{v}\times \mathit{B}),\hfill \\ \hfill {b}^{\mu}& =\gamma (\mathit{v}\xb7\mathit{B},\phantom{\rule{0.166667em}{0ex}}\mathit{B}\mathit{v}\times \mathit{E}).\hfill \end{array}\end{array}$$(6)
The above expression for J is finally obtained from the projection j^{μ} = J^{μ} + (J^{ν}u_{ν})u^{μ}, where J^{μ} = (q, J) is the fourcurrent (Bucciantini & Del Zanna 2013; Del Zanna & Bucciantini 2018; Mignone et al. 2024).
Following Mattia et al. (2023), the set of RRMHD equations is closed by assuming the Taub equation of state (Mignone & McKinney 2007; Mizuno 2013), which is the appropriate one when a relativistic fluid (the jet) interacts with a nonrelativistic medium
$$\begin{array}{c}\hfill h=\frac{5}{2}\mathrm{\Theta}+\sqrt{\frac{9}{4}{\mathrm{\Theta}}^{2}+1},\end{array}$$(7)
where Θ = p/ρ is a temperature function.
In addition, the evolution of a passive scalar tracer ξ is also considered, in order to properly disentangle the jet from the ambient medium (the tracer will be initialized only at the injection boundary). By combining the appropriate advection equation for ξ with the continuity equation for D, the code actually evolves in time the following additional conservative equation
$$\begin{array}{c}\hfill \frac{\partial (D\xi )}{\partial t}+\mathrm{\nabla}\xb7(D\xi \mathit{v})=0.\end{array}$$(8)
2.2. Numerical setup
All the simulations have been performed with the PLUTO code (Mignone et al. 2007) on a static mesh in spherical coordinates (r, θ, ϕ). We refer to the code’s user guide for the recipes employed in the code to handle nonCartesian geometries (in particular, the ϕ component of the momentum equation is multiplied by r). In some of the plots, cylindrical coordinates (R = r sin θ, z = rcosθ) are introduced when a generic meridional plane is displayed.
The numerical algorithms adopted in this paper show a few slight differences from the ones described in Mattia et al. (2023). In particular, we carry out our simulations using the IMplicitEXplicit RungeKutta SSP3 (3, 3, 2) scheme (Pareschi & Russo 2005; Palenzuela et al. 2009) coupled with a fifthorder (instead of fourthorder) Piecewise Parabolic reconstruction method (Mignone 2014). Instead of the Harten, Lax, and van Leer Riemann solver (Harten et al. 1983) adopted in Mattia et al. (2023), we chose to employ an adaptation of the FORCE Riemann solver (Mattia & Mignone 2022) where the fast waves are approximated with the speed of light (which does not only represent the maximum speed reachable but also the propagation of the fast waves in the frozen limit, as shown in Mignone et al. 2018). As in Mattia et al. (2023), we preserve the solenoidality of the magnetic field through the divergence cleaning algorithm. We note that the divergence cleaning method requires an additional evolutionary equation
$$\begin{array}{c}\hfill \frac{\partial \psi}{\partial t}+{c}_{h}^{2}\mathrm{\nabla}\xb7\mathit{B}=\frac{{c}_{h}^{2}}{{c}_{p}^{2}}\psi ,\end{array}$$(9)
which, as in Mignone et al. (2010a), is solved in two steps:

the advective step is solved within the IMEX framework,

the source term is solved using its analytical solution,
$$\begin{array}{c}\hfill {\psi}^{\mathrm{\Delta}{t}^{n}}={\psi}^{(0)}exp(\frac{{c}_{h}^{2}}{{c}_{p}^{2}}\mathrm{\Delta}t).\end{array}$$(10)
For this simulations, we set c_{h} = 1, c_{p}^{2} = Δh/α^{2}, α = 0.01, and Δh = min(Δr). We point out that the GLM method also impacts the evolutionary equation of the magnetic field as well, which becomes
$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}+\mathrm{\nabla}\times \mathit{E}+\mathrm{\nabla}\psi =0.\end{array}$$(11)
Conversely, the divergence of the electric field is not explicitly enforced, as in Del Zanna et al. (2016), Tomei et al. (2020), Mattia et al. (2023).
Due to the high Lorentz factors involved, the implicit IMEX step may become brittle due to the high nonlinearity of the RRMHD equations. In particular, we found two potential issues:

As already shown in Mattia et al. (2023), if the fourvelocity is nonnegligible in any of the components, convergence may not be reached due to an endless increase in the Lorentz factor. Therefore, if an iteration during the implicit step becomes higher than 10^{3} (since, as shown in the following subsections, we expect the maximum Lorentz factor achievable by our simulations to be < 300) we reset the fourvelocity to 0 and the other primitive variables to the ones of the previous timestep;

The initial guess provided for the implicit step may not always be physical, since it is computed from the explicit step of the IMEX scheme. Such a step, without the implicit update, is still incomplete and therefore the conserved variable resulting from it may not represent a physically valid configuration. If the velocity computed from the intermediate step exceeds the speed of light we set its initial guess to 0.9 (in units of c). Such a fix involves only the velocity since all the other fluid variables do not require a guess during the implicit step.
The code robustness is additionally enhanced by employing the shock flattening methods present in the PLUTO code (see the code’s user guide).
2.3. Numerical grid
In this paper, we present two types of simulations. The first kind is a 2D axisymmetric parameter study where half of the angular domain is explored, namely θ ∈ [0, π/2] covered by N_{θ} = 256 grid points, and therefore no counterjet is shown. The second kind is a half 3D grid (i.e., with no counterjet, as for the 2D simulations) rotated by π/2 so that the jet symmetry axis is on the equatorial plane, in order to fully explore the jet dynamics without the potential numerical issues or artifacts coming from the polar axis singularity. In 3D simulations the angular directions θ ∈ [0.1, π − 0.1] and ϕ ∈ [0, π] are covered by, respectively, N_{θ} = 480 and N_{ϕ} = 512 grid points in order to preserve the resolution adopted in the 2D axisymmetric simulations.
The radial direction is spaced in the same way for all the simulations and, as shown in Fig. 1, is designed to capture the dissipative effects closer to the excision radius and to maintain their impact up to the larger scales. To achieve this goal, we constructed a uniform grid from r_{min} = 562.6 km to r_{bound} = 2036.3 km covered by 512 grid points (yielding a spacing of Δr ≈ 3 km). On top of the uniform grid, we constructed a stretched grid up to r_{max} = 95439.3 km, reaching the coarsest grid spacing of ∼10^{3} km near the outer boundaries. The 2D and 3D simulations are run, respectively until t = 10^{5} and t = 2 × 10^{4} in code units, which corresponds to t = 0.4926 s and t = 0.0985 s.
Fig. 1. Radial spacing as a function of the distance from the center. The zoomed panel highlights the boundary between the uniform and the stretchedspaced regions, plotted, respectively, in blue and red for the sake of clarity. 
In the remainder, we describe both the initial conditions for t = 0 and the jet’s injection recipes for t > 0. For the sake of simplicity, we just describe the setup for the 2D axisymmetric case. In the 3D case, the vector fields are tilted by 90° to arrange the orbital axis orthogonally to the cylindrical zaxis, as anticipated. For instance, the rotated expression for the magnetic field becomes (the full calculation can be found in Appendix C of Pavan et al. 2023)
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathit{B}=& {B}^{r}{\mathit{e}}_{{r}^{\prime}}+({B}^{\theta}sin\psi sin\varphi +{B}^{\varphi}\frac{cos\psi}{sin\theta}){\mathit{e}}_{{\theta}^{\prime}},\hfill \\ \hfill & +({B}^{\theta}\frac{cos\psi}{sin\theta}+{B}^{\varphi}sin\psi sin\varphi ){\mathit{e}}_{{\varphi}^{\prime}},\hfill \end{array}\end{array}$$(12)
where here the primes indicate the new 3D reference system, with
$$\begin{array}{c}\hfill \psi =\mathrm{sign}(cos\theta )arccos\left(\frac{\mathrm{sign}(cos\varphi sin\theta )}{\sqrt{1+{(cos\varphi tan\theta )}^{2}}}\right),\end{array}$$(13)
and the same rotation is applied to the velocity.
2.4. Initial conditions
Due to the highly complex structure of the postmerger magnetic field in realistic BNS merger simulations, a detailed numerical study of resistive/dissipative processes in such cases is rather challenging. In this paper, for the sake of simplicity, the jet is launched into an idealized analytic environment, similar to the test case presented (along with the more realistic one) in Pavan et al. (2021). We postpone to a subsequent work (Paper II) the study of resistive effects within the realistic postmerger environment directly imported from the outcome of a BNS merger simulation.
We start by defining a radius r_{b} = 7383.2 km, which represents the (initial) boundary between the radially expanding, magnetized postmerger environment and the external, unperturbed, unmagnetized static atmosphere. The density profile is defined as
$$\begin{array}{c}\hfill \rho \propto \{\begin{array}{ccc}{\left({\displaystyle \frac{r}{{r}_{min}}}\right)}^{s}\hfill & \phantom{\rule{2em}{0ex}}& r\le {r}_{\mathrm{b}}\hfill \\ {\left({\displaystyle \frac{r}{{r}_{\mathrm{b}}}}\right)}^{6.5}\hfill & \phantom{\rule{2em}{0ex}}& r\ge {r}_{\mathrm{b}},\hfill \end{array}\end{array}$$(14)
with continuity for r_{b}, where ρ(r_{min}) = 1.056 × 10^{8} g cm and s = 3.981. A similar formula applies for the pressure, with p(r_{min}) = 6.408 × 10^{25} dyn cm and s = 3.320.
The postmerger wind is described by a radial velocity with a homologous increase up to r_{b}
$$\begin{array}{c}\hfill {v}^{r}/c=0.047\left(\frac{r}{{r}_{min}}\right)0.037,\end{array}$$(15)
while we assume v = 0 in the static atmosphere. As for the velocity, we assume B = 0 in the atmosphere, while the postmerger magnetic field is assumed to be a dipole up to r_{b}, with
$$\begin{array}{c}\hfill \begin{array}{c}\hfill {B}^{r}=2{B}_{\mathrm{d}}cos\theta {\left(\frac{{r}_{min}}{r}\right)}^{3}\\ \hfill {B}^{\theta}={B}_{\mathrm{d}}sin\theta {\left(\frac{{r}_{min}}{r}\right)}^{3},\end{array}\end{array}$$(16)
with B_{ϕ} = 0, where B_{d} = 4.013 × 10^{13} G. Finally, the electric field is initialized everywhere to its ideal value (E = −v × B) (thus null in the static atmosphere), while the passive tracer is set to 0, given that the jet will be injected from the boundary only for t > 0, as described below.
2.5. Boundary conditions and jet injection
Due to the different types of simulations presented in this paper, the boundary conditions are also depending on the configuration adopted. In the 3D simulation, we adopted, respectively, zerogradient (outflow) and periodic boundaries for the θ′ and ϕ′ coordinates. Conversely, in the 2D simulations, we impose axisymmetric and equatorial symmetric boundary conditions at θ = 0 and θ = π/2, respectively. The outer radial boundaries are set to outflow in all the simulations, while the inner radial boundaries are defined in order to inject the magnetized jet after the collapse phase (30 ms after the NS collapse) into the computational domain. The recipes for the injection are summarized in the following steps (a more complete derivation can be found in Pavan et al. 2023), in which dimensional units are introduced:

We set the jet opening angle to θ_{j} = 10°, and hence the jet is injected where θ < θ_{j}, and we also set the tracer ξ = 1 there (zero elsewhere).

We set the injection parameters as in Pavan et al. (2023). In particular, the injection and termination Lorentz factors are set, respectively to γ_{j} = 3 and γ_{∞} = 300, the ratio between the maximum toroidal and radial magnetic field strength is set to B_{rat} = 0.5, the jet angular velocity is set to $\overline{\mathrm{\Omega}}\sim 10$ rad/s, the maximum azimuthal magnetic field is set to B_{j} = 1.55 × 10^{13} G, the magnetic field halfopening angle is set to θ_{j, m} = 4°, and the jet density at the edge is set to ρ(θ_{j}) = 2.203 × 10^{3} g/cm^{3}.

We compute the injection velocity and magnetic field at the initial injection time
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{v}}_{0}& =(c\sqrt{1{\gamma}_{j}^{2}},0,\overline{\mathrm{\Omega}}rsin\theta )\hfill \\ \hfill {\mathit{B}}_{0}& =({B}_{\mathrm{rat}}{B}_{j},0,\frac{2{B}_{j}(\theta /{\theta}_{j,m})}{1+{(\theta /{\theta}_{j,m})}^{2}}),\hfill \end{array}\end{array}$$(17)
and we also compute b^{2}(θ) using the expression for b^{μ} in the case of the ideal MHD approximation, where we reintroduce the correct units,
$$\begin{array}{c}\hfill {b}^{\mu}=\frac{1}{\sqrt{4\pi}}[\frac{\gamma (\mathit{v}\xb7\mathit{B})}{c},\frac{\mathit{B}}{\gamma}+\frac{\gamma (\mathit{v}\xb7\mathit{B})\mathit{v}}{{c}^{2}}].\end{array}$$(18)

We compute the jet density by assuming transverse equilibrium between total pressure gradient, centrifugal force, and magnetic tension,
$$\begin{array}{c}\hfill \rho (\theta )=[\rho ({\theta}_{j})+{\displaystyle {\int}_{{\theta}_{j}}^{\theta}}exp(F(\theta \prime ))g(\theta \prime )d\theta \prime ]\phantom{\rule{0.166667em}{0ex}}exp(F(\theta )),\end{array}$$(19)
where
$$\begin{array}{c}\hfill F(\theta )=\frac{4{h}_{0}^{\ast}{\gamma}^{2}{\overline{\mathrm{\Omega}}}^{2}{r}_{min}^{2}}{{c}^{2}({h}_{0}^{\ast}1)}\frac{{sin}^{2}(\theta ){sin}^{2}({\theta}_{j})}{2}ln\frac{sin(\theta )}{sin({\theta}_{j})},\end{array}$$(20)
and
$$\begin{array}{c}\hfill g(\theta )=\frac{4{({b}^{\varphi})}^{2}{b}^{2}+sin\theta \frac{d{b}^{2}}{d\theta}}{{c}^{2}({h}_{0}^{\ast}1)sin\theta},\end{array}$$(21)
with h_{0}^{*} = γ_{∞}/γ_{0}. Subsequently, the jet’s initial twosided luminosity is computed as
$$\begin{array}{c}\hfill {L}_{j}=4\pi {r}_{min}^{2}{\displaystyle {\int}_{0}^{{\theta}_{j}}}[\rho {h}_{0}^{\ast}{\gamma}_{j}^{2}{c}^{2}(p+\frac{{b}^{2}}{2}){({b}^{0})}^{2}]{v}^{r}sin\theta d\theta ,\end{array}$$(22)
where the initial pressure is computed as
$$\begin{array}{c}\hfill p=[\rho {c}^{2}({h}_{0}^{\ast}1)+{b}^{2}]/4.\end{array}$$(23)
In all the simulations, the jet parameters yield a luminosity of L_{j} = 1.08 × 10^{52} erg/s, which corresponds to an initial jet local magnetization of
$$\begin{array}{c}\hfill \sigma =\frac{{b}^{2}}{4\pi \rho {c}^{2}h}\end{array}$$(24)
between σ_{min} = 0.066 and σ_{max} = 11.52, depending on the angle θ. We note that, here, we are using the fully relativistic definition of the (hot) magnetization parameter, as appropriate for highspeed flows. The corresponding relative magnetic luminosity contribution is
$$\begin{array}{c}\hfill \frac{{L}_{j}{L}_{j,\mathrm{HD}}}{{L}_{j}}\sim 30\%,\end{array}$$(25)
where L_{j, HD} is the jet luminosity without the electromagnetic field contribution. Therefore, at the injection site, the jet can be considered as strongly magnetized (Gottlieb et al. 2020). We note that this step is computed only at the beginning of the simulation and not at every evolutionary step.

We apply the decay of the luminosity by adopting
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \gamma (t)& =\sqrt{1+{({\gamma}_{j}{v}_{0}^{r}{e}^{t/(2{\tau}_{d})}/c)}^{2}}\hfill \\ \hfill {v}^{r}(t)& =\frac{{\gamma}_{j}{v}_{0}^{r}{e}^{t/(2{\tau}_{d})}}{\gamma (t)}\hfill \\ \hfill {B}^{\varphi}(t)& ={B}_{0}^{\varphi}\sqrt{\frac{{v}_{0}^{r}}{{v}^{r}(t)}}{e}^{t/(2{\tau}_{d})}\gamma (t)\hfill \\ \hfill {B}^{r}(t)& ={B}_{0}^{r}\sqrt{\frac{{v}_{0}^{r}}{{v}^{r}(t)}}{e}^{t/(2{\tau}_{d})}\gamma (t)\hfill \\ \hfill {h}^{\ast}(t)& =\frac{{h}_{0}^{\ast}{\gamma}_{j}{e}^{t/(2{\tau}_{d})}}{\gamma (t)},\hfill \end{array}\end{array}$$(26)
where the quantities with the label 0 refer to the initial time, using the values specified in steps 3 and 4, while τ_{d} = 300 ms represents the temporal decay timescale. The pressure is computed again from Eq. (23) with the timeevolved values of the magnetic field (in both the laboratory and comoving frames) and velocity.

We compute the electric field under the ideal approximation E = −v × B.
In the full 3D case, we note that all vectors must be rotated according to the recipes previously mentioned. The condition for the jet’s injection simply becomes
$$\begin{array}{c}\hfill {r}_{\mathrm{cyl}}=rcos\theta \sqrt{1+{(cos\varphi tan\theta )}^{2}}<=rsin{\theta}_{j}.\end{array}$$(27)
2.6. Gravity
As shown in the simulations of Pavan et al. (2021), in the vicinity of the denser parts of the BNS remnant, gravity plays a role and must be considered. Following Pavan et al. (2021, 2023), gravity is treated like an external force in the RMHD equations, with the corresponding acceleration g which is not simply provided by the standard Newtonian term due to the central mass M_{0}/M_{⊙} = 2.596 of the remnant, but it also takes into account the (outward) radial pressure gradients caused by hot and magnetized material during the BH formation, a term that is gradually decreasing in time. As done in the cited works, we model this scenario by providing an effective mass and then we use
$$\begin{array}{c}\hfill \mathit{g}=\frac{G[{M}_{0}{M}_{\mathrm{eff}}(r,t)]}{{r}^{2}}\widehat{\mathit{r}}.\end{array}$$(28)
The explicit value of M_{eff} is computed by first recovering its value at r = r_{min} and the collapse time t = t_{c} s by solving the angleaveraged radial momentum equilibrium, i.e., yielding
$$\begin{array}{c}\hfill {M}_{\mathrm{eff}}({r}_{\mathrm{min}},{t}_{c})={\left[{r}^{2}\frac{1}{\overline{\rho}G}(\frac{\mathrm{d}\overline{p}}{\mathrm{d}r}+\frac{1}{8\pi}\frac{\mathrm{d}\overline{{B}^{2}}}{\mathrm{d}r}\frac{\overline{(\mathit{B}\xb7\mathrm{\nabla}){B}_{r}}}{4\pi})\right]}_{{r}_{\mathrm{min}},0}.\end{array}$$(29)
Consequently, the radial and temporal dependence is provided by:
$$\begin{array}{c}\hfill {M}_{\mathrm{eff}}={M}_{0}max(0.1\frac{r{r}_{\mathrm{min}}}{{r}_{G}{r}_{\mathrm{min}}}){M}_{\mathrm{eff}}({r}_{\mathrm{min}},{t}_{c})\phantom{\rule{0.166667em}{0ex}}{\mathrm{e}}^{\frac{t+{\tau}_{j}}{\tau}},\end{array}$$(30)
where r_{G} ∼ 1033.65 km is the gravity characteristic radius, τ_{j} = −t_{c} = 30 ms represents the delay time between the neutron stars collapse and the jet launch, and
$$\begin{array}{c}\hfill \tau ={\tau}_{j}+({\tau}_{d}{\tau}_{j}){sin}^{2}\theta \end{array}$$(31)
represents the gravitational timescale^{2} with τ_{j} = 0.3 s the accretion timescale of the black holedisk system. We note that, in all the simulations, we adopt the pointmass approximation for Newtonian gravity.
3. Axisymmetric parameter study
To understand the impact of the resistivity, we performed a set of different axisymmetric 2D simulations, thus reducing the necessary computational cost. Despite the issues described in the literature (e.g., LópezCámara et al. 2013; Harrison et al. 2018) related to the propagation of the jet along the axis at θ = 0, the relatively short computational time required to evolve the jet up to r ∼ 90 000 km allowed us to perform a parameter study. To properly disentangle the impact of the magnetic dissipation on the jet propagation and dynamics, we chose 4 configurations, two with Lundquist numbers (defined as S = v_{a}L/η where v_{a} is the Alfvén speed and L = 150 km is the length scale) S = 10^{2} (acronym S2) and two with S = 10^{3} (acronym S3). For each chosen value of the Lundquist number, we performed one simulation assuming v_{a} = c (acronym C), that is, with constant resistivity in time and space, and one where we recover v_{a} in its exact form (acronym V):
$$\begin{array}{c}\hfill {v}_{A}/c=max(\frac{\mathit{B}}{\sqrt{\rho {c}^{2}h+{\mathit{B}}^{2}}},0.001),\end{array}$$(32)
where the lower limit is set to avoid potential cells where η = 0, which may not be properly solved by the IMEX scheme (Bugli et al. 2014; Ripperda et al. 2019a). We show, in the left panel of Fig. 2, the values of the resistivity at the final time of each simulation as a function of space, overlapped by the contour lines of the passive scalar tracer (dashed and dotted, respectively corresponding to the values of 0.1 and 0.01). As a preliminary consideration, we note that since within the jet, the Alfvén speed is close to the speed of light, cases with the same Lundquist number yield similar values of resistivity (only within the jet). Additionally, we report a schematic summary of the key parameters in Table 1.
Fig. 2. Resistivity η (left panel) and restmass density ρ (right panel) at the final time of the simulation (t ∼ 0.5 s) for the different resistivity cases. The dashed black line represents the tracer contour line at level 1%, while the solid lines are only designed to delimit the different cases. We note that in this and all the 2D plots, the angle θ is between 0° and 90° and the negative cylindrical coordinates are simply a postprocessing transformation to enhance the comparison quality. 
Relevant parameters for the simulations discussed in the paper.
Within the single cases with variable resistivity, the jet is more diffusive when compared to the ambient medium, due to the higher values of the jet’s local Alfvén speed. Therefore, due to the impact of the numerical resistivity (see Appendix A), the case S3V is dominated by the physical nonideal processes only within the jet (e.g., where the Alfven speed is closer to the speed of light) or at shorter distances, while the numerical one dominates only beyond the interface between the jet and the ambient medium at larger distances due to the lower resolution caused by the stretched grid (as shown next subsection and in Appendix A). Conversely, the dynamics of the case S2V is still dominated by the physical resistivity also outside the jet and through the entire radial domain extension. We point out that, due to the magnetization of the postmerger wind, resistive processes can also occur at larger distances from the injection site and in the postmerger environment.
The density at the final time of the simulations is shown in the right panel of Fig. 2. First, we note a ringlike structure, which corresponds to the “slow” waves traveling in the expanding wind (the quotes are adopted since the resistivity affects the wave propagation, as shown in Mignone et al. 2018), which is expanding into the unmagnetized atmosphere. Due to the expansion of the magnetized environment into the atmosphere, different propagating waves yield a series of shocks and discontinuities in the external regions, which propagate at different characteristic speeds. Since the propagation of the environment is, although slower than the jet, still relativistic, at the final time of the simulation the jet does not break into the unmagnetized atmosphere but travels into different environment regions characterized by the propagation of the slow waves. This propagation of the environment leads to an accumulation of material and a magnetic field, which affects the jet propagation and dynamics. Due to the different behavior between these two propagation stages, we report the time of the slow breakout in Table 1. The impact of the propagation of the jet into the strongly magnetized (inner) environment region and the weakly magnetized (outer) environment region is discussed in the following subsections of this paper.
Moreover, the magnetic resistivity affects the shape and the position (see the following subsection) of the termination jet lobes and the jet backflow, as well as the shape and number of the recollimation regions present (although we point out that these regions cannot be associated with the m = 0 plasma columns instabilities since, due to the axisymmetric assumption, the jet shapes like a funneled column subjected to the higher pressure from the ambient medium from both inside and outside the jet). In particular, we note that a constant resistivity is associated with a more uniform dissipation of the magnetic field close to the lobes, resulting in less (maximum 2 in each C case) recollimation regions. Conversely, a variable (in time and space) resistivity leads to a more fragmented outer jet (3 and 5 recollimation zones can be detected, respectively, in the cases S2V and S3V). While in these simulations the jet shape is constrained by the assumption of axisymmetry, we expect such features to impact the jet dynamics in fully 3D simulations strongly (see Sect. 4), both on the jet dynamics and the radiation signature associated (e.g., Dong et al. 2020).
Finally, in Fig. 3 we show the magnetization (defined as in Eq. 24). Due to the dissipation caused by the resistivity, all the cases show a lower magnetization at larger distances and much lower values than the magnetization injected at t = 0. Moreover, the resistivity affects the jet magnetization both at the jet axis and the termination regions, yielding higher values of σ for lower values of η. In particular, following the criteria of Bromberg & Tchekhovskoy (2016), the S3 cases show a moderate jet magnetization (σ ≲ 1) even at larger distances, while the S2 cases can be considered as weakly magnetized (i.e., σ ≪ 1, more similar to the jets investigated in Gottlieb et al. 2020). A comparison with an ideal run is done in Appendix A, while a more quantitative analysis of the magnetic energy dissipation is performed in the following section.
Fig. 3. Magnetization σ at the final time of the simulation (t ∼ 0.5 s) for the different resistivity cases. 
3.1. Jet energy and dynamics
As shown in Mattia et al. (2023), the magnetic resistivity is expected to play a role in terms of the shape and dynamics of potential current sheets, energy conversion (from magnetic to thermal), and production of unaligned electric field (from the ideal −v × B direction) and dissipated power. Due to the lower grid resolution in the outer domain regions, the development of smallscale (i.e., few km) turbulent structures is strongly suppressed for t ≳ 0.1 s in our simulations. Nevertheless, the physical diffusion of the electromagnetic field remains dominant over the numerical one (see Appendix A for a comparison with an ideal run where the electric field is imposed to be E = −v × B).
The top and bottom left panels of Fig. 4 show the evolution of the jet electromagnetic energy with respect to the case S3C at the final time and the nonelectromagnetic energy (i.e., the thermal + kinetic without the restmass contribution, see, e.g., Ricci et al. 2024):
Fig. 4. Time evolution of the jet energetics and dynamics. The electromagnetic energy contribution (compared to that of the case S3C at the final time in the top left panel and to the kinetic + internal energy in the bottom left panel) and the jet front (position in the top panel, and Lorentz factor in the bottom panel) are shown for different 2D runs. The circle dot represents the time of slow breakout. 
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {E}_{\mathrm{int}}& ={\displaystyle \int \xi}[\rho {\gamma}^{2}(h1){c}^{2}p]{r}^{2}sin\theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi ,\hfill \\ \hfill {E}_{\mathrm{mag}}& ={\displaystyle \int \xi}\left[\frac{{E}^{2}+{B}^{2}}{2}\right]{r}^{2}sin\theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi ,\hfill \\ \hfill {E}_{\mathrm{kin}}& ={\displaystyle \int \xi}\left[D(\gamma 1){c}^{2}\right]{r}^{2}sin\theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi .\hfill \end{array}\end{array}$$(33)
As a first consideration, we note that the S2 cases, those shown in the topleft panel of Fig. 4, show a quasiconstant behavior until t ∼ 0.2 s and (only the S2C case) a decrease around t ∼ 0.05 s. Such a lack of increase is the result of the strong diffusivity processes that balance the continuous energy injection at the inner boundary. However, at later times (i.e., when the jet has pierced the remnant and starts propagating through the weakly magnetized environment region), the diffusivity plays a much less crucial role, which is due to the lack of the magnetic field in that region, leading to an increase in the electromagnetic energy in every case during the later evolutionary stages (despite its decay at the injection site described by Eq. 26). On the other hand, the S3 cases show an increase in electromagnetic energy up to t ∼ 0.4 s, with the S3C case showing such an increase throughout the entire evolution. As for the previous case, two distinct phases can be detected: in the early stage of the simulation (i.e., t ≲ 0.1 s), the dissipative processes play a crucial role in the interaction between the magnetized wind and the jet; once the jet has passed the slow breakout radius, the magnetic resistivity becomes less relevant, because the external environment region is very weakly magnetized and the injection process becomes dominant. The reason behind the “plateau” visible at the end of the simulation S3V may be related to a shift in the dominant diffusive mechanism, since at larger distances the numerical resistivity may overcome the physical one in this particular case. Therefore, an uncontrolled numerical resistivity may slow even further the increase of the magnetic energy related to the jet injection.
In the bottom left panel of Fig. 4 the impact of electromagnetic energy with respect to the kinetic and thermal jet energy is shown. We note that the jet remains hydrodynamically (thermally) dominated since the electromagnetic contribution remains below 7% for all the simulations and settles around 4% and ≲2% respectively for the S3 cases and the S2 cases. In the very initial evolutionary stages, the electromagnetic contribution to the jet’s total energy increases (concerning both its initial value and the nonelectromagnetic contribution, as shown also in the top left panel); however, this increase of the electromagnetic energy fraction is reverted already at t ∼ 0.015 s, due to the combination of the interaction between the jet and the remnant medium and the magnetic dissipation caused by the magnetic resistivity.
All the cases (but S2C) show a similar peak in the electromagnetic energy impact, while the suppression of the magnetic energy contribution (which is reverted around t ∼ 0.1 s for the S3 cases and t ∼ 0.2 s for the S2V case) yields a different minimum before a subsequent increase due to the propagation of the jet in the weakly magnetized environment region. At this point, due to the lack of a magnetic field, the jet propagation proceeds unaltered (i.e., no strong turbulence is present outside the outer lobes). As a result, the only physical processes able to impact the electromagnetic field are magnetic dissipation and jet injection, leading to a slow increase in the contribution of electromagnetic energy. However, as stated previously, due to the lack of magnetization, electromagnetic field injection becomes the dominant process, especially for low Lundquist numbers.
We point out that these simulations cannot reveal a heating process by solely looking at the energy output, since energy conversion can also occur outside the resistive regime due to the simple interaction between the jet and the surrounding environment. However, this is a clear dissipation process triggered by the magnetic resistivity which affects the impact of the electromagnetic component of the jet energy up to ∼30% on the jet energy balance. Nevertheless, as pointed out at the beginning of the section, the jet propagation and dynamics are also affected by the magnetic resistivity. Due to the energy balance resulting from the different simulations, one would expect that a more dissipative jet (because of the energy conversion) would propagate on different timescales through the external medium. If on one hand, the dissipation of the magnetic field quenches the magnetic acceleration, on the other hand, the heating dissipative processes are able to increase the heat reservoir converted to kinetic energy during the jet expansion. In order to assess more quantitatively the impact of the magnetic dissipation, we computed (see the top right panel of Fig. 4) the jet front as the maximum radial distance from the core at which the Lorentz factor is above 2 (note that the remnant has an initial maximum Lorentz factor of ∼1.22, well below our constraint). The position of the jet front is also used to compute the time at which the slow breakout occurs, that is, when the jet front has pierced the discontinuity associated with the environment’s slow wave propagation. The time and position of the slow breakout are reported in Table 1.
As expected, once the jet has reached the weakly magnetized environment region, the resistivity impacts the jet propagation by ≈10 000 km (which corresponds to approximately the 11% of the distance covered by the jet at t ∼ 0.5 s). We note that the jet front spatial propagation is not linear, suggesting a nontrivial temporal evolution of its propagation. By looking at the Lorentz factor associated with the jet front propagation, shown in the bottom right panel of Fig. 4, we can divide the jet propagation into three different stages:

An initial stage (t < 0.05 s), wherein the jet is still confined within the (more dense) postmerger environment, and the jet front Lorentz factor decreases since the matter is slowed down by both the accumulation of matter at the jet axis and the turbulence caused by the interaction between the jet and the external medium in all our simulations.

A slow breakout stage (0.05 s < t < 0.2 s), where the jet breaks through the postmerger environment and starts propagating through the weakly magnetized environment region; all the simulations show a similar increase in the Lorentz factor.

A final stage (t > 0.2 s), where the jet propagates through the lowdensity and lowpressure weakly magnetized environment region, keeping the same internal structure. The resistivity strongly affects this stage because the magnetic resistivity is the only magnetic process that can impact the jet energy.
We note that the positions and times of the slow breakout show very little or no difference between the cases if compared with the evolution at the later stages. Such division differs partially from the one obtained by looking at the electromagnetic energy because changes in the electromagnetic energy do not immediately reflect on their hydrodynamical counterpart (and most of the Ohmic dissipation takes place in the inner radii, see next subsection). Nevertheless, the three evolutionary stages detected in the bottom left panel of Fig. 4 are consistent with the different trends in jet propagation.
The different trend of the Lorentz factor leads to a difference of about ∼7% between the S2 and the S3 cases at the final time of our simulation. Since this relative difference in the Lorentz factor increases with time (at t = 0.1 and t = 0.3 the differences are, respectively, around ∼1% and 3%), we expect the dynamics of jets characterized by different resistivity models, to show differences also at larger times.
3.2. Parallel field and dissipated power
The magnetic resistivity is also responsible for the generation and amplification of electric field which deviates from the −v × B assumption (Del Zanna et al. 2016), as shown in Fig. 5 as dissipated power (left panel, normalized to the maximum value obtained from the S2V case at the end of the simulation) and parallel electromagnetic field (right panel), computed, respectively, as (see also Eq. (6))
$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{P}& =\eta {j}^{2}={e}^{2}/\eta ={\gamma}^{2}\frac{{(\mathit{E}+\mathit{v}\times \mathit{B})}^{2}{(\mathit{v}\xb7\mathit{E})}^{2}}{\eta},\hfill \\ \hfill E{B}_{\Vert}& =\frac{\mathit{E}\xb7\mathit{B}}{\sqrt{{E}^{2}{B}^{2}}+\u03f5},\hfill \end{array}\end{array}$$(34)
Fig. 5. Dissipated power (left panel) and normalized parallel component of the electric magnetic field (right panel) for the different resistivity cases at t ∼ 0.5 s. The zoomed regions show the inner region up to 3000 km in both the R and z directions for both panels. 
where ϵ serves the purpose of avoiding undefined divisions.
In agreement with Mattia et al. (2023), a stronger resistivity yields a higher dissipated power in the cases where the resistivity is constant in time and space. In all cases, most of the power is dissipated within the innermost region and slowly decays as the jet propagates until the termination region (i.e., ∼3 × 10^{4} km, where the 𝒫 is several orders of magnitude below the innermost one). However, the power is not dissipated only within the jet, since for all the cases a nonnegligible power is dissipated in the remnant ambient region. In particular, a filamentary structure arising from the equator (θ = 90°) is present. The structure of such filaments is much more complex in the variable resistivity cases, where multiple gradients are present. In contrast, in the constant resistivity cases, the structure outside the jet shows less turbulence and a much smoother behavior. This is not surprising since narrow regions of high local Alfvén speed, outside the jet, may yield a higher level of dissipated power that propagates outwards, leading to weakly resolved features whose spatial size corresponds to few numerical cells and depends on the local magnetic resistivity. Moreover, the S3V case shows significant dissipated power outside the jet and also at greater distances from the inner boundary (near both the equator and the jet axis) and a higher level of dissipated power when compared to the S3C case. Such a feature is not present in any of the other cases, and its nature seems to be numerical (similar issues in the resistive electric field were found in the low resistivity cases of Mattia et al. 2023).
By looking at the right panel of Fig. 5, we note that the S2C case yields a much more diffused “parallel field region” than the other 3 cases, especially close to the midplane region and in the weakly magnetized environment region. Such a feature is a direct consequence of the strong magnetic resistivity through the entire domain, meaning not only within the jet but also in the magnetized ambient medium. As expected, the S2 cases show a stronger parallel component close to the jet axis, extending to ∼30 000 km (the maximum jet position along the axis θ = 0). Conversely, the S3 cases do not lead to a strong nonideal electric field (≳1% of its magnitude) within the jet and show deviations from the ideal regime only from the interface between the jet and the ambient medium, where the KelvinHelmholtz instability and magnetic reconnection are likely to take place and generate and amplify resistive electric field.
All the cases (but S3V, likely due to the low resistivity outside the jet) show a clear double parallel field inversion: within the jet, the resistive electric field features a parallel component, which becomes antiparallel outside the interface between the jet and the ambient. Additionally, another resistive field reversal, with different magnitudes, is present in the inner weakly magnetized environment regions. Such inversion of the resistive electric field is clearly enhanced in the presence of a stronger physical resistivity. and is strictly correlated to a velocity inversion (radial between the jet and remnant, toroidal within the remnant region, see the next section) At the early evolutionary stages, such inversion can trigger the development of the KelvinHelmholtz Instability (KHI, Chow et al. 2023) which is a favorable environment for the amplification of nonideal magnetic field and magnetic dissipation and reconnection (Sironi et al. 2021). Conversely, the S3V case shows no inversion at the jetambient interface, although a strong gradient in the parallel electromagnetic field magnitude is still present.
3.3. Initial evolution
To provide a more meaningful investigation of short GRB jets, we need to combine the scales at which the resistive processes take place with the large propagation scale reached by the jet after the surrounding environment is pierced. Although the lower resolution, at large distances, smears out much of the turbulence present within the jet, we can still investigate the impact of the resistivity more in detail when looking at the jet at early times (e.g., t ∼ 0.1 s, which corresponds to the 20% of the final time of our simulations). In particular, the electromagnetic features of the resistive simulations can provide a clear indication of magnetic reconnection and particle acceleration sites. The ratio between the electric and magnetic field is reported in the left panel of Fig. 6, superimposed by the contour lines of the passive tracer at level 0.1.
Fig. 6. Electromagnetic field properties at the early simulation stages. The ratio between the electric and magnetic energy (left panel) and the radial component of the magnetic field (right panel) are shown for the different resistivity cases, at t = 0.1 s, superimposed by the contour lines of the scalar tracer at 0.1. 
Here the jet structure shows, already at this evolutionary stage, significant differences in shape, both inside and outside the regions where the jet has interacted with the postmerger environment. The interrelation between the jet magnetic field and its dynamics has been extensively studied in a more general context (Leismann et al. 2005; MoyaTorregrosa et al. 2021) and in more specific jet types, such as AGN jets (Ricci et al. 2024) and GRBjets (Gottlieb et al. 2020; Pavan et al. 2023) The region affected by the interaction between the jet and the external medium shows a high level of turbulence and a ratio E^{2}/B^{2} between 1% and 10% (as shown in the left panel of Fig. 6), suggesting that magnetic reconnection (which is likely to occur in such location due to the KHI and the plasma column instabilities) may play a crucial role in particle acceleration (as already shown in Beniamini & Giannios 2017) in GRBs as well.
Reversals in the radial component of the magnetic field are also an indicator of potential magnetic reconnection sites. Due to the injected magnetic field structure, we expect the magnetic field to play a significant role in the jet structure and collimation. The presence and the shape of magnetic field reversal zones, where the magnetic field reaches a minimum due to a change of polarity, are strongly affected by the resistivity adopted (Ripperda et al. 2019b; Mattia et al. 2023). In particular, since the formation of thin magnetic layers with different polarity takes place where the jet and the postmerger environment have strongly interacted, the physical resistivity must overcome the numerical one so that the reconnecting region does not depend on the resolution and the numerical algorithms adopted. This is the case for the S3V model, where such regions are not resolved and their width is, for the most part, limited to very few computational cells. Conversely, the cases S3C and S2V show regions favorable for magnetic reconnection associated with wider magnetic layers. In addition, in the turbulent ambient region (where the Alfvén speed is of the order of ∼10% of the speed of light), we find that the resistivity associated with the reconnecting regions is comparable for the cases S3C and S2V (as in the left panel of Fig. 2), while within the jet (where the Alfvén speed approaches c) the cases with the same Lundquist number also show a similar resistivity Lastly, the S2C case resistivity is too strong to favor the formation of turbulent reconnecting regions, yielding a more diffused structure similar to the one found in the more diffusive cases of Mattia et al. (2023).
4. Full 3D simulation
As shown in fully 3D numerical simulations of relativistic jets (e.g., Mignone et al. 2010b; LópezCámara et al. 2013; Bromberg & Tchekhovskoy 2016), the assumption of axisymmetry can be broken already at the early stages of jet propagation. The absence of axisymmetric boundary conditions, in particular, prevent the unphysical accumulation of dense material along the jet injection axis, which can alter the propagation itself leading to unrealistic configurations of the jet angular structure (see, e.g., Fig. 1 in Mignone et al. 2010b). Moreover, the turbulent interaction (which we expect to be strongly affected by the resistive and dissipative processes considered in this paper) between the jet and its surrounding environment can fully develop only in 3D. For these reasons, we performed a fully 3D simulation with analytical initial and boundary conditions (as for the 2D axisymmetric simulations) until t = 0.1 s, which is approximately the time at which the jet started piercing the postmerger environment and reached a distance of r ∼ 10^{4} km from the central object.
In Fig. 7, we compare a 3D simulation cut (at ϕ = π/2 and θ ≥ π/2) with the corresponding 2D axisymmetric case^{3}. On the left panel, we show the density (top part) and the radial magnetic field (bottom part), while on the right panel, the parallel electromagnetic field (top part) and the dissipated power (bottom part) are shown. Here the radial and scalar components are unaffected by the setup rotation, so the comparison does not require any coordinate transformation. We note that the accumulation of material closer to the jet axis here present is not related to numerical instabilities, but to the injection of a lighter jet into a more massive environment. Due to the onset of the KHI, which takes place already in the merger (Kiuchi et al. 2015) and the earlier injection phase and represents a powerful magnetic field amplification mechanism, the 3D case shows a higher level of turbulence at the layers between the jet and its surroundings. In particular, a stronger radial magnetic field with steeper gradients is present, as shown by the multiple field reversals. A stronger magnetic field can be associated with more efficient magnetic reconnection and particle acceleration (Sironi & Spitkovsky 2014).
Fig. 7. Comparison against the axisymmetric and the fully 3D simulations. The density (top part) and radial magnetic field (bottom part) are shown in the left panel, while the normalized parallel electromagnetic field (top part) and the jetdissipated power (bottom part) are shown in the right panel. The left part of each panel represents the case S2V, while the right part represents the cut of the 3D simulation at ϕ = π/2 and θ ≥ π/2. 
Due to the enhanced turbulence, the 3D simulation shows a different resistive electromagnetic field (Fig. 7, right panel). More specifically, closer to the axis the electric field is quasiperpendicular to the magnetic field (the parallel component is below 0.01%), while outside the jet the turbulent structure is, as for the previous variables here considered, more prominent. Here charged particles can also be accelerated due to the nonperpendicularity of the electromagnetic field, provided that their gyroradius is much smaller than the variations in the electromagnetic field. Since the minimum jet magnetic field is approximately B_{min} ≈ 10^{11} G, the maximum gyroradius for a charged particle traveling at the speed of light would be 3.3 × 10^{−4} m for a TeV electron/proton, with a linear dependence of the gyroradius on the particle energy.
In Fig. 8, we report the azimuthal distribution at different radii of several fluid variables, to assess more deeply the level of axisymmetry present in this simulation (which, as described in the previous section, starts with an axisymmetric jet launched into an axisymmetric environment). All the right panels, which correspond to a radial distance of 5238.9 km from the center, show an instability with mode m = 4. Such a feature should be related to the RayleighTaylor (RT) instability, triggered by the accelerated jet expansion in the denser environment, which drives mixing between the two fluids (Mignone et al. 2010b; Gottlieb et al. 2021). For this to occur, the magnetic tension due to the toroidal component of the magnetic field should not be the dominant one (in the presence, as in our case, of a magnetized jetenvironment). Here, the ratio between the toroidal component (retrieved by applying the inverse transformation of Eq. 12) and total magnetic field, shown in the second row of Fig. 8, can help to understand the origin of the instability itself better. At lower radii (left panel), two regions where the magnetic field is toroidally dominated are detected: within the jet and the surrounding cocoon; such zones are separated by a thin poloidally dominated layer. At larger radii, the toroidal component is no longer always dominant in the second region, and the RT instability is free to develop there, leading to a strong mixing of material with the external medium.
Fig. 8. Jet angular profiles at different radii. From left to right, the distance from the center is 1357.5 km, 2715.6 km, and 5238.9 km, respectively. The restmass density, toroidal over total magnetic field, jet rotation velocity, and normalized parallel electromagnetic field are shown from top to bottom. The dashed black lines represents the cut at ϕ = π/2 and θ ≥ π/2, as in Fig. 7. 
We also computed the jet rotation, shown in the third row of Fig. 8. A double inversion in the jet rotation is visible here: the inner rotation change at the jet boundary and an additional flip in the rotation direction at larger angles. We note that such a trend is recovered also in the 2D simulations, where there is no equivalent behavior of the magnetic field. At larger distances (middle column) the jet rotation becomes more “patchy” with a general counterclockwise predominance and some localized regions with clockwise movement. Once the jet has approached its peak (right column), the m = 4 instability has fully taken place, and the backflow of material is divided into four regions, each with two opposite rotation directions. The presence of such strong gradients in shear velocities cannot exclude a contribution to the enhancement of the turbulence level due to KHI, which also in the relativistic regime, as in classical MHD, is affected by the relative importance of magnetic field components, aligned or transverse to the shear (Bucciantini & Del Zanna 2006).
Finally, the parallel electromagnetic field EB_{∥} is shown in the last row of Fig. 8. Here, the symmetry with respect to the jet injection axis is broken already at r = 1357.5 km, due to the higher level of turbulence outside the jet, as seen in Fig. 7. At larger distances, the parallel electromagnetic field becomes weaker due to the magnetic field diffusion. Nevertheless, few regions where particles can undergo nonthermal acceleration are present even at r = 5238.9 km. This indicates that resistivity becomes even more relevant in fully 3D simulations due to their enhanced turbulence (compared with 2D simulations).
In conclusion, Fig. 9 shows the 3D distributions of the Lorentz factor (left panel) and magnetic energy density (right panel), in logarithmic scale, in the lab frame at the end of the simulation. A clear anticorrelation between the two quantities can be noted. In particular, a stronger magnetic field is present closer to the jet front, while at the inner propagation region (i.e., where the jet has maximum velocity) the magnetic energy density has lower values. As shown in Pavan et al. (2023) and our axisymmetric simulations, we expect the magnetic energy to be converted into kinetic energy after the jet has pierced the environment and broken out into the weakly magnetized environment region. Moreover, where the magnetic energy has not been converted into kinetic energy, the jet shows a stronger accumulation of the former at the boundaries between the jet and the environment, while the core remains weakly magnetized. Such hollow core is also found in different simulations (e.g., Nathanail et al. 2021). Here, its origin is associated with the lateral expansion of the jet due to the interaction with the denser environment since the jet has not broken into the less dense atmosphere.
Fig. 9. Fully 3D jet structure. The logarithms of the Lorentz factor and the magnetic energy density (in G^{2} units) are displayed in the left and right panels, respectively. 
By looking at the hydrodynamic (left panel of Fig. 9) and the electromagnetic (right panel) jet components, we note a similar shape in the presence of enough jet material (which is correlated with a high Lorentz factor). However, due to the environment’s magnetic field, the turbulent interaction outside the jet cone yields a more turbulent structure and extends much further in the electromagnetic component, suggesting that reconnection and dissipation processes can more efficiently take place in such a scenario. The enhanced turbulent level, coupled with a more complex behavior of the resistive quantities, is expected to be enhanced even more by injecting the jet into a nonaxisymmetric realistic environment.
5. Conclusions
In this paper, we present the first resistive relativistic simulations of short GRB jets. By combining a suite of highorder numerical algorithms with a twocomponent radial grid, we were able to properly capture the dissipative effects occurring at the early propagation stages, that is, when the jet has not yet pierced the magnetized remnant. Our results can be summarized as follows:

We performed four different 2D axisymmetric simulations up to t ∼ 0.5 s with varying prescriptions of resistivity (i.e., = 10^{2} or S = 10^{3} and uniform or variable resistivity), allowing us to disentangle the resistive processes during the jet propagation.

A higher resistivity is related to a higher degree of electromagnetic dissipation, yielding a more thermal jet that propagates faster through the remnant and the weakly magnetized environment region. While at the earlier propagation stages, the resistivity plays a key role mostly in the electromagnetic jet components, once the jet has pierced the magnetized remnant environment, the resistivity becomes crucial not only in the diffusion of the electromagnetic field but also in the position and velocity of the jet front.

Most resistive effects, that is, the generation and amplification of parallel electromagnetic field and dissipated power, are seen in the innermost jet region. However, a filamentary structure is visible closer to the equator, whose gradients strongly depend on the resistivity profile. The parallel electromagnetic field shows a double inversion in most cases, with a tight interrelation between a lower diffusivity and steeper gradients.

Already at the early evolutionary stages, magnetic resistivity strongly affects the electromagnetic jet counterpart. In particular, a lower resistivity yields a higher level of turbulence and plasma instabilities at the interface between the jet and the ambient medium, which, coupled with the strong electric field, makes them an upandcoming site for the acceleration of nonthermal particles. Conversely, a higher resistivity is associated with a smoother magnetic field whose configuration is less favorable to the onset of magnetic reconnection.

We performed a fully 3D simulation up to t ∼ 0.1 s in order to assess the deviation from axisymmetry even in a fully axisymmetric initial configuration, as well as the features related to the absence of boundaries at the jet axis. As a consequence of the loss of axisymmetry, there is no accumulation of material at the jet axis or amplification of the magnetic turbulence farther from it, yielding a more turbulent electric field outside the jet.

At larger distances from the injection sites, the jet undergoes a plasma instability of mode 4 due to the RayleighTaylor instability, which affects its shape and generates a filamentary backflow that undergoes a series of recollimation regions. We also find that the resistive components of the electromagnetic field are more susceptible to the breaking of the axisymmetry, strongly deviating from axisymmetry already closer to the jet base, indicating that fully 3D simulations will augment the impact of the resistivity on the jet evolution.
In conclusion, the importance of physical resistivity, which is not dependent on the grid resolution or the suite of numerical algorithms adopted (provided that it is higher than the numerical dissipation), is made clear by the findings presented here, as is the importance of fully 3D simulations. Due to the complexity of the physical processes involved here, which should encompass very different lengths and timescales, the grid resolution can strongly impact the outcome in the presence of only an uncontrolled numerical resistivity. Conversely, adopting a physically motivated, consistent resistivity model makes it possible to more accurately reproduce the observational features found in nature.
In a companion paper (Paper II), we will explore the impact of the resistivity on jet propagation in a more realistic nonaxisymmetric environment, which will be imported consistently from BNS numerical relativity simulations as in Pavan et al. (2021, 2023, and in prep.).
Data availability
The PLUTO code is publicly available. Simulation data will be shared at a reasonable request by the corresponding author.
A typo is present in Pavan et al. (2021, 2023), although this affects only the manuscript and not the corresponding code.
Acknowledgments
We thank the anonymous referee for the valuable comments and suggestions that greatly improved the quality of the manuscript. GM acknowledges D. Crocco for the assistance and valuable comments on the Python routines. LDZ acknowledges support from the ICSC–Centro Nazionale di Ricerca in HighPerformance Computing; Big Data and Quantum Computing, funded by European Union – NextGenerationEU. AP and RC acknowledge support from the European Union under NextGenerationEU, via the PRIN 2022 Project “EMERGE: Neutron star mergers and the origin of short gammaray bursts”, Prot. n. 2022KX2Z3B (CUP C53D23001150006), and further support by the Italian Ministry of Foreign Affairs and International Cooperation (MAECI), grant number US23GR08. All the simulations were performed on GALILEO100 at CINECA (Bologna, Italy). In particular, we acknowledge CINECA for the availability of high performance computing resources and support through awards under the ISCRA initiative (Grant IsB27_BALJET) and through a CINECA–INFN agreement (providing the allocation INF24_teongrav). The 1D and 2D figures presented in this paper were generated using the (https://github.com/GiMattia/PyPLUTO) PyPLUTO package (Mattia et al., in prep.); volume renderings were generated using ParaView (Ahrens et al. 2005).
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJ, 848, L13 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Phys. Rev. Lett., 119, 161101 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, ApJ, 848, L12 [Google Scholar]
 Ahrens, J., Geveci, B., & Law, C. 2005, Visualization Handbook (Elsevier) [Google Scholar]
 Beniamini, P., & Giannios, D. 2017, MNRAS, 468, 3202 [Google Scholar]
 Bromberg, O., & Tchekhovskoy, A. 2016, MNRAS, 456, 1739 [Google Scholar]
 Bucciantini, N., & Del Zanna, L. 2006, A&A, 454, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bucciantini, N., & Del Zanna, L. 2013, MNRAS, 428, 71 [Google Scholar]
 Bugli, M., Del Zanna, L., & Bucciantini, N. 2014, MNRAS, 440, L41 [Google Scholar]
 Chow, A., Davelaar, J., Rowan, M. E., & Sironi, L. 2023, ApJ, 951, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Ciolfi, R. 2020, MNRAS, 495, L66 [Google Scholar]
 Ciolfi, R., Kastaun, W., Giacomazzo, B., et al. 2017, Phys. Rev. D, 95, 063016 [Google Scholar]
 Ciolfi, R., Kastaun, W., Kalinani, J. V., & Giacomazzo, B. 2019, Phys. Rev. D, 100, 023005 [Google Scholar]
 Del Zanna, L., & Bucciantini, N. 2018, MNRAS, 479, 657 [NASA ADS] [Google Scholar]
 Del Zanna, L., Papini, E., Landi, S., Bugli, M., & Bucciantini, N. 2016, MNRAS, 460, 3753 [Google Scholar]
 Del Zanna, L., Tomei, N., Franceschetti, K., Bugli, M., & Bucciantini, N. 2022, Fluidika, 7, 87 [Google Scholar]
 Dionysopoulou, K., Alic, D., & Rezzolla, L. 2015, Phys. Rev. D, 92, 084064 [Google Scholar]
 Dong, L., Zhang, H., & Giannios, D. 2020, MNRAS, 494, 1817 [Google Scholar]
 Franceschetti, K., & Del Zanna, L. 2020, Universe, 6, 83 [Google Scholar]
 GarcíaGarcía, L., LópezCámara, D., & Lazzati, D. 2023, MNRAS, 519, 4454 [Google Scholar]
 Geng, J.J., Zhang, B., Kölligan, A., Kuiper, R., & Huang, Y.F. 2019, ApJ, 877, L40 [Google Scholar]
 Ghirlanda, G., Salafia, O. S., Paragi, Z., et al. 2019, Science, 363, 968 [Google Scholar]
 Gottlieb, O., Bromberg, O., Singh, C. B., & Nakar, E. 2020, MNRAS, 498, 3320 [Google Scholar]
 Gottlieb, O., Nakar, E., & Bromberg, O. 2021, MNRAS, 500, 3511 [Google Scholar]
 Gottlieb, O., Moseley, S., RamirezAguilar, T., et al. 2022, ApJ, 933, L2 [Google Scholar]
 Hamidani, H., & Ioka, K. 2023, MNRAS, 520, 1111 [Google Scholar]
 Harrison, R., Gottlieb, O., & Nakar, E. 2018, MNRAS, 477, 2128 [Google Scholar]
 Harten, A., Lax, P., & Leer, B. 1983, SIAM Rev., 25, 35 [Google Scholar]
 IndaKoide, M., Koide, S., & Morino, R. 2019, ApJ, 883, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Kathirgamaraju, A., Tchekhovskoy, A., Giannios, D., & Barniol Duran, R. 2019, MNRAS, 484, L98 [Google Scholar]
 Kiuchi, K., CerdáDurán, P., Kyutoku, K., Sekiguchi, Y., & Shibata, M. 2015, Phys. Rev. D, 92, 124034 [Google Scholar]
 Kiuchi, K., ReboulSalze, A., Shibata, M., & Sekiguchi, Y. 2024, Nat. Astron., 8, 298 [Google Scholar]
 Komissarov, S. S. 2007, MNRAS, 382, 995 [Google Scholar]
 Lazzati, D., Perna, R., Morsony, B. J., et al. 2018, Phys. Rev. Lett., 120, 241103 [Google Scholar]
 Lazzati, D., Perna, R., Ciolfi, R., et al. 2021, ApJ, 918, L6 [Google Scholar]
 Leismann, T., Antón, L., Aloy, M. A., et al. 2005, A&A, 436, 503 [Google Scholar]
 LópezCámara, D., Morsony, B. J., Begelman, M. C., & Lazzati, D. 2013, ApJ, 767, 19 [CrossRef] [Google Scholar]
 Margalit, B., & Metzger, B. D. 2017, ApJ, 850, L19 [Google Scholar]
 Mattia, G., & Mignone, A. 2022, MNRAS, 510, 481 [Google Scholar]
 Mattia, G., Del Zanna, L., Bugli, M., et al. 2023, A&A, 679, A49 [Google Scholar]
 Mignone, A. 2014, J. Comput. Phys., 270, 784 [Google Scholar]
 Mignone, A., & McKinney, J. C. 2007, MNRAS, 378, 1118 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mignone, A., Tzeferacos, P., & Bodo, G. 2010a, J. Comput. Phys., 229, 5896 [Google Scholar]
 Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010b, MNRAS, 402, 7 [Google Scholar]
 Mignone, A., Mattia, G., & Bodo, G. 2018, Phys. Plasmas, 25, 092114 [Google Scholar]
 Mignone, A., Mattia, G., Bodo, G., & Del Zanna, L. 2019, MNRAS, 486, 4252 [Google Scholar]
 Mignone, A., Berta, V., Rossazza, M., et al. 2024, MNRAS, 533, 1670 [Google Scholar]
 Mizuno, Y. 2013, ApJS, 205, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018, Nature, 561, 355 [Google Scholar]
 Mösta, P., Radice, D., Haas, R., Schnetter, E., & Bernuzzi, S. 2020, ApJ, 901, L37 [Google Scholar]
 MoyaTorregrosa, I., Fuentes, A., Martí, J. M., Gómez, J. L., & Perucho, M. 2021, A&A, 650, A60 [Google Scholar]
 MurguiaBerthier, A., RamirezRuiz, E., De Colle, F., et al. 2021, ApJ, 908, 152 [Google Scholar]
 Nathanail, A., Gill, R., Porth, O., Fromm, C. M., & Rezzolla, L. 2020, MNRAS, 495, 3780 [Google Scholar]
 Nathanail, A., Gill, R., Porth, O., Fromm, C. M., & Rezzolla, L. 2021, MNRAS, 502, 1843 [Google Scholar]
 Palenzuela, C. 2013, MNRAS, 431, 1853 [Google Scholar]
 Palenzuela, C., Lehner, L., Reula, O., & Rezzolla, L. 2009, MNRAS, 394, 1727 [Google Scholar]
 Palenzuela, C., AguileraMiret, R., Carrasco, F., et al. 2022, Phys. Rev. D, 106, 023013 [Google Scholar]
 Pareschi, L., & Russo, G. 2005, J. Sci. Comput., 25, 129 [Google Scholar]
 Pavan, A., Ciolfi, R., Kalinani, J. V., & Mignone, A. 2021, MNRAS, 506, 3483 [Google Scholar]
 Pavan, A., Ciolfi, R., Kalinani, J. V., & Mignone, A. 2023, MNRAS, 524, 260 [Google Scholar]
 Puzzoni, E., Mignone, A., & Bodo, G. 2021, MNRAS, 508, 2771 [Google Scholar]
 Puzzoni, E., Mignone, A., & Bodo, G. 2022, MNRAS, 517, 1452 [Google Scholar]
 Qian, Q., Fendt, C., & Vourellis, C. 2018, ApJ, 859, 28 [Google Scholar]
 Ricci, L., Perucho, M., LópezMiralles, J., Martí, J. M., & Boccardi, B. 2024, A&A, 683, A235 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ripperda, B., Bacchini, F., Porth, O., et al. 2019a, ApJS, 244, 10 [Google Scholar]
 Ripperda, B., Porth, O., Sironi, L., & Keppens, R. 2019b, MNRAS, 485, 299 [Google Scholar]
 Ruiz, M., Lang, R. N., Paschalidis, V., & Shapiro, S. L. 2016, ApJ, 824, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Sironi, L., & Spitkovsky, A. 2014, ApJ, 783, L21 [Google Scholar]
 Sironi, L., Rowan, M. E., & Narayan, R. 2021, ApJ, 907, L44 [Google Scholar]
 Tomei, N., Del Zanna, L., Bugli, M., & Bucciantini, N. 2020, MNRAS, 491, 2346 [Google Scholar]
 Urrutia, G., De Colle, F., MurguiaBerthier, A., & RamirezRuiz, E. 2021, MNRAS, 503, 4363 [Google Scholar]
 Urrutia, G., De Colle, F., & LópezCámara, D. 2023, MNRAS, 518, 5145 [Google Scholar]
 Vourellis, C., Fendt, C., Qian, Q., & Noble, S. C. 2019, ApJ, 882, 2 [Google Scholar]
 Xie, X., Zrake, J., & MacFadyen, A. 2018, ApJ, 863, 58 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Comparison with the ideal run
To assess the contribution of the physical resistivity, as opposed to the numerical one, we performed two additional runs with S = 10^{4} (S4V and S4C) and a run by solving the ideal equations with the same numerical methods adopted for the resistive simulations (except for the time integration, which is now an explicit 3^{rd}order RungeKutta method). We note that, despite the same general structure, the FORCE Riemann solver differs from its resistive version since the computation of the fastest wave propagation, as well as the fluxes, is not the same.
A.1. Energy budget and magnetic turbulence
The numerical resistivity, which is intrinsically present in every simulation, may differ, especially at larger distances from the merger where the grid is less refined due to the stretching. Such difference is clearly visible in Fig. A.1, where the magnetic energy is reported as a function of time (similarly to Fig. 4). Here we compared only cases with constant resistivity (in time and space) to neglect the impact of the jet dynamics on the resulting resistivity. Although these simulations are less realistic due to the simplified magnetic resistivity model adopted, they provide a better insight into the impact of the numerical diffusion. All the cases but the S4C one show a very good agreement with a linear dependence, provided that the jet has pierced the postmerger environment (i.e., after t ∼ 0.1 s). Due to the grid stretching, we expect the numerical resistivity to increase at larger radii, eventually overcoming the physical one for very high values of the Lundquist number S. This seems to be the case of the S4C run, where a strong trend change can be detected around t ∼ 0.25. We note that the energy increase of the case S4C and S3V (top left panel of Fig. 4), despite both being dominated by the numerical diffusion at large distances, show a very different behavior, very likely because of the interaction between the numerical and physical resistivity. We note that, before t ∼ 0.25 s, the ideal and S4C cases almost overlap, indicating that the contributions of the numerical and physical (only in the S4C simulation) are comparable. Moreover, due to the unpredictable nature of the numerical resistivity, the ideal case slope may strongly change with different grid resolutions or more/less accurate numerical algorithms.
Fig. A.1. Electromagnetic energy contribution (compared to the one of the case S3C at the final time in the top panel and to the kinetic + internal energy in the bottom panel) for different values of the Lundquist number and constant resistivity (colored lines) and ideal run (black line). The dashed lines are the linear approximation of the evolution of electromagnetic energy at later times. The circle dot represents the slow breakout. 
The impact of the resistivity on the role of the jet’s electromagnetic energy (compared to the thermal and kinetic energy) is shown in the bottom panel of Fig. A.1. Here a higher Lundquist number is related to more (although the jet is dominated by the thermal energy in all simulations) magnetized jets, due to the lower amount of diffusivity. As shown in the top panel, the ideal and S4C cases show a similar behavior up to t ∼ 0.25 s, where the latter case becomes more dominated by the numerical resistivity. We note that the slow breakout occurs at slightly later times in the S4C and ideal cases (see also Table 1), although no additional statement can be made due to the intricate interrelation between the jet and the environment.
As in Mattia et al. (2023), we also checked the impact of the resistivity on the turbulence resolution by computing the normalized gradient < ∇B^{2}/B^{2}> for different runs with both constant and variable resistivity and the ideal run. As shown in Fig. A.2, the turbulence level between the ideal and the resistive cases differs by a factor ≳2. We also note that the similarities between the cases S2V and S3C are also reflected here since the two cases show a similar level of turbulence at the end of the simulations. The same statement holds for the cases S3V and S4C, which both seem to be affected by the numerical diffusion at larger distances. From the energy budget, we suggested that both S3V and S4C cases are affected by the numerical diffusion in the outer domain regions (where the resolution is the poorest); this hypothesis is strengthened by the very good agreement in terms of turbulence refinement that the two simulations reach. Conversely, all the other simulations show a level of turbulence that is dictated by a physical, and therefore controllable, resistivity.
Fig. A.2. Effective turbulence resolution of the magnetic energy for different resistivity and ideal cases. The circle dot represents the slow breakout. 
Fig. A.3. Electromagnetic field properties at the early and late simulation stages. The radial component of the magnetic field is shown, for the different constant resistivity and ideal cases, at t = 0.1 s (left panel) and t = 0.1 s (right panel). 
A.2. The electromagnetic field
The main difference between the resistive and the ideal runs is the presence/absence of physical resistivity controlling the level of dissipation throughout the entire simulation. The shape and size of current sheets or steep magnetic field gradients are strongly affected by the amount of dissipation present throughout the simulation. Therefore, in Fig. A.3, we compared the radial component of the magnetic field among different resistive simulations (with constant resistivity) and the corresponding ideal simulation. The left and right panel show, respectively, the jet propagation at t = 0.1 s and t = 0.5 s.
Both bottom panels, corresponding to the S2C and S3C simulations, show thicker (when compared with the other two cases) regions of the magnetic field with the same polarity, both at the early and later stages of the simulation. Nevertheless, the radial magnetic field shows a comparable magnitude in all the simulations, regardless of the employed resistivity. We note that, in agreement with Fig. A.2, the S2C case shows the lowest turbulence resolution due to the high diffusion in both the jet and the ambient medium.
Conversely, the S4C and IDEAL cases show a more refined turbulence due to the lack of diffusion. As expected, at lower distances from the center (left panel), the two cases yield a similar level of dissipation, confirming that the numerical diffusivity can be assessed to be around a Lundquist number of S = 10^{4}. At larger distances and times (right panel), the case S4C suffers from higher numerical diffusion (compared to the ideal case) caused by the more diffusive Riemann solver and timestepping algorithm. However, the two simulations remain very qualitatively similar and are affected by less dissipation than the cases with high resistivity, confirming that all the cases (but S3V) feature a higher physical resistivity (than the numerical one) at all the stages of the simulations.
Finally, in Fig. A.4, we report the spatial distribution of the magnetization σ, computed as in Eq. 24. In agreement with the previous subsection, we found that the highest magnetization levels are found in the S4C and ideal cases due to the weak diffusion of the magnetic field. In both cases, the jet is moderately magnetized at the final stages of the simulation, although in the termination stage (where the magnetic energy has not been fully converted to kinetic energy) the jet shows a stronger magnetization. We note that the upper cases (S4C and IDEAL) show when compared to the lower cases (S2C and S3C) a stronger magnetization in both the injection, propagation, and outer phases due to the lack of resistivity.
Fig. A.4. Magnetization σ at the final time of the simulation (t ∼ 0.5 s) for the constant resistivity and ideal cases. 
All Tables
All Figures
Fig. 1. Radial spacing as a function of the distance from the center. The zoomed panel highlights the boundary between the uniform and the stretchedspaced regions, plotted, respectively, in blue and red for the sake of clarity. 

In the text 
Fig. 2. Resistivity η (left panel) and restmass density ρ (right panel) at the final time of the simulation (t ∼ 0.5 s) for the different resistivity cases. The dashed black line represents the tracer contour line at level 1%, while the solid lines are only designed to delimit the different cases. We note that in this and all the 2D plots, the angle θ is between 0° and 90° and the negative cylindrical coordinates are simply a postprocessing transformation to enhance the comparison quality. 

In the text 
Fig. 3. Magnetization σ at the final time of the simulation (t ∼ 0.5 s) for the different resistivity cases. 

In the text 
Fig. 4. Time evolution of the jet energetics and dynamics. The electromagnetic energy contribution (compared to that of the case S3C at the final time in the top left panel and to the kinetic + internal energy in the bottom left panel) and the jet front (position in the top panel, and Lorentz factor in the bottom panel) are shown for different 2D runs. The circle dot represents the time of slow breakout. 

In the text 
Fig. 5. Dissipated power (left panel) and normalized parallel component of the electric magnetic field (right panel) for the different resistivity cases at t ∼ 0.5 s. The zoomed regions show the inner region up to 3000 km in both the R and z directions for both panels. 

In the text 
Fig. 6. Electromagnetic field properties at the early simulation stages. The ratio between the electric and magnetic energy (left panel) and the radial component of the magnetic field (right panel) are shown for the different resistivity cases, at t = 0.1 s, superimposed by the contour lines of the scalar tracer at 0.1. 

In the text 
Fig. 7. Comparison against the axisymmetric and the fully 3D simulations. The density (top part) and radial magnetic field (bottom part) are shown in the left panel, while the normalized parallel electromagnetic field (top part) and the jetdissipated power (bottom part) are shown in the right panel. The left part of each panel represents the case S2V, while the right part represents the cut of the 3D simulation at ϕ = π/2 and θ ≥ π/2. 

In the text 
Fig. 8. Jet angular profiles at different radii. From left to right, the distance from the center is 1357.5 km, 2715.6 km, and 5238.9 km, respectively. The restmass density, toroidal over total magnetic field, jet rotation velocity, and normalized parallel electromagnetic field are shown from top to bottom. The dashed black lines represents the cut at ϕ = π/2 and θ ≥ π/2, as in Fig. 7. 

In the text 
Fig. 9. Fully 3D jet structure. The logarithms of the Lorentz factor and the magnetic energy density (in G^{2} units) are displayed in the left and right panels, respectively. 

In the text 
Fig. A.1. Electromagnetic energy contribution (compared to the one of the case S3C at the final time in the top panel and to the kinetic + internal energy in the bottom panel) for different values of the Lundquist number and constant resistivity (colored lines) and ideal run (black line). The dashed lines are the linear approximation of the evolution of electromagnetic energy at later times. The circle dot represents the slow breakout. 

In the text 
Fig. A.2. Effective turbulence resolution of the magnetic energy for different resistivity and ideal cases. The circle dot represents the slow breakout. 

In the text 
Fig. A.3. Electromagnetic field properties at the early and late simulation stages. The radial component of the magnetic field is shown, for the different constant resistivity and ideal cases, at t = 0.1 s (left panel) and t = 0.1 s (right panel). 

In the text 
Fig. A.4. Magnetization σ at the final time of the simulation (t ∼ 0.5 s) for the constant resistivity and ideal cases. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.