Resistive relativistic MHD simulations of astrophysical jets

Aims. The main goal of the present paper is to provide the first systematic numerical study of the propagation of astrophysical relativistic jets, in the context of high-resolution shock-capturing resistive relativistic magnetohydrodynamics (RRMHD) simulations. We aim at investigating different values and models for the plasma resistivity coefficient, and at assessing their impact on the level of turbulence, the formation of current sheets and reconnection plasmoids, the electromagnetic energy content, and the dissipated power. Methods. We use the PLUTO code for simulations and we assume an axisymmetric setup for jets, endowed with both poloidal and toroidal magnetic fields, and propagating in a uniform magnetized medium. The gas is assumed to be characterized by a realistic Synge-like equation of state (Taub equation), appropriate for such type of astrophysical jets. The Taub equation is combined here for the first time with the Implicit-Explicit Runge-Kutta time-stepping procedure, as required in RRMHD simulations. Results. The main result is that turbulence is clearly suppressed for the highest values of resistivity (low Lundquist numbers), current sheets are broader, and plasmoids are barely present, while for low values of resistivity results are very similar to ideal runs, where dissipation is purely numerical. We find that recipes employing a variable resistivity based on the advection of a jet tracer or on the assumption of a uniform Lundquist number improve on the use of a constant coefficient and are probably more realistic, preserving the development of turbulence and of sharp current sheets, possible sites for the acceleration of the non-thermal particles producing the observed high-energy emission.


Introduction
Astrophysical jets, i.e. supersonic collimated outflows, represent a ubiquitous phenomenon in the Universe, characterizing various classes of celestial sources at very different spatial scales and evolutionary stages.In particular, the interplay of strong gravity, rotation, and magnetic fields of compact objects may cause the launching of relativistic jets, even with high Lorentz factors, as spectacularly demonstrated by their apparent superluminal motion when propagation is close to the line of sight, and more recently by the imaging of the supermassive black hole as the inner engine of the kiloparsec-scale jet of M87 (Event Horizon Telescope Collaboration 2019).Other than in Active Galactic Nuclei (AGN), relativistic outflows and jets are also present in several types of stellar-sized astrophysical bodies, such as the collapsing sources of long Gamma-Ray Bursts (GRB) (Woosley 1993), Pulsar Wind Nebulae (PWN) (Weisskopf et al. 2000), and recently we had the first multi-messenger indication that also Binary Neutron Star (BNS) merger events are able to drive relativistic outflows in the form of short GRBs (Abbott et al. 2017).The electromagnetic spectral signature of such objects often features highly non-thermal radiation, mostly synchrotron and in-verse Compton processes, due to accelerated relativistic particles (electrons).In the case of AGNs, particles may be accelerated either in the inner engine and advected, or along the jet and at the terminal radio lobes (e.g.Begelman et al. 1984;Blandford et al. 2019).More recently, even protostellar jets have shown similar radiation features (Lee et al. 2018).
Different processes have been invoked in order to explain such particles acceleration in astrophysical sources.In particular, relativistic magnetic (fast) reconnection, i.e. the impulsive topological rearrangement of field lines in a magnetically dominated or very hot plasma (Lyutikov & Uzdensky 2003;Lyubarsky 2005;Komissarov et al. 2007;Del Zanna et al. 2016;Ripperda et al. 2019b) is a very general and promising acceleration process, which has been extensively investigated in recent years, especially through particle-in-cell (PIC) simulations (Sironi & Spitkovsky 2014;Cerutti et al. 2014;Comisso & Sironi 2018;Zhang et al. 2021).The precise mechanism regulating the acceleration of particles during the reconnection process, namely the formation of X and O points and the subsequent merging of plasmoids (e.g.Loureiro & Uzdensky 2016), is still debated, as well as their impact on the large-scale scenario (e.g.Medina-Torrejón et al. 2023).
Because of the unavoidable numerical resistivity present in every MHD simulation (even within the ideal framework), the disentanglement of the effects due to the truncation error from the physical process is a key point in order to understand the role of magnetic reconnection in astrophysical jets.The most critical issue of the ideal approximation (which constraints the electric field to the standard −v × B) is the lack of control on the numerical diffusivity, which becomes strongly susceptible on the grid resolution and numerical algorithms (Núñez-de la Rosa & Munz 2016;Felker & Stone 2018).
It is therefore important to perform numerical simulations of astrophysical plasmas going beyond the usual assumption of ideal MHD, almost invariably employed especially in the relativistic regime.The simplest possible framework in which the role of a finite plasma conductivity can be investigated is that of the resistive relativistic magnetohydrodynamics (RRMHD) regime, in which the electric field in the coming frame is not vanishing as in the ideal MHD case but is simply proportional to the local current density through a (scalar for simplicity) resistivity coefficient η.By doing so, the electric field in the laboratory frame is allowed to depart from its ideal configuration (which imposes a vanishing component along the magnetic field lines) especially in current sheets, where reconnection events and particle acceleration are expected to take place.Obviously any current numerical modeling will never achieve the realistic low values of η typical of collisionless astrophysical plasmas, though it is well known that small-scale unresolved turbulent fluctuations may act as sources of additional dissipation (see the Conclusions for additional comments on sub-grid modeling).
Here we just want to stress that the RRMHD regime allows to single out the regions where resistive effects may be important, by adopting a physically motivated model.This is instead impossible in numerical ideal MHD, where magnetic dissipation is invariably present at the grid level (due to numerical round-off errors), and any non-ideal aspect cannot be controlled.Moreover, employing an explicit resistivity allows for more meaningful comparisons between models produced with different codes and numerical algorithms, as the uncertainties related to magnetic numerical dissipation are removed.Thus, this approach promotes the reproducibility of jet simulations and strengthens the physical interpretation of their results.Most of the RRMHD numerical investigations of reconnection, either due to the tearing instability in pre-existing current sheets (e.g.Del Zanna et al. 2016), or formed in the process of a turbulent cascade (e.g.Chernoglazov et al. 2021), involve localized portions of plasma and periodic boundary conditions.Other studies assume reconnection induced by current-driven kink instabilities in (periodic) columns of magnetized plasma (Medina-Torrejón et al. 2021;Bodo et al. 2021).A few investigations are concerned with the accretion disk and launching region at the base of the jet, some of them also addressing the problem of the flaring activity observed in Sgr A * (Qian et al. 2017;Vourellis et al. 2019;Ripperda et al. 2020;Nathanail et al. 2022;Ripperda et al. 2022), though the last works assume ideal conditions.Other numerical investigations combine RRMHD with mean-field dynamo action in order to explain the growth of the magnetic field in accretion disks up to the required equipartition values (Tomei et al. 2020(Tomei et al. , 2021;;Del Zanna et al. 2022).
On the other hand, to our knowledge there are no studies of propagating jets using RRMHD simulations with finite plasma conductivity.The aim of the present work is hence precisely to perform the first systematic numerical investigation of propagat-ing magnetized jets via RRMHD simulations.The combination of high-Mach number collimated ouflows, strong magnetization, and turbulence development is particularly challenging in nonideal relativistic MHD.Robust schemes and high resolution are both required, and this must be the reason why such effort has not been faced yet, to our knowledge.
In the present work we perform and discuss high-resolution axisymmetric simulations (in cylindrical coordinates and flat Minkowskian metric) of jets endowed with both poloidal and toroidal magnetic fields, propagating in a uniform magnetized medium.Simulations are obtained with the PLUTO1 shockcapturing code (Mignone et al. 2007(Mignone et al. , 2012)).The gas is assumed to be characterized by a realistic Synge-like equation of state (Taub equation), believed to be appropriate for such type astrophysical jets (Mignone & McKinney 2007).The Taub equation of state is combined here for the first time with Implicit-Explicit Runge-Kutta (Pareschi & Russo 2005;Palenzuela et al. 2009) routines for time-stepping, which allow a proper treatment of stiff terms in the evolutionary equation for the electric field in the limit of small resistivity coefficients.
Various levels of magnetization and plasma beta are investigated, and different values and models for the resistivity coefficient (assumed to be a scalar, but not necessarily a constant) are proposed, taking care that numerical diffusivity remains lower than the physical one.Results are compared, in terms of morphology and turbulence level in the jet.Regions of high values of the current density and the possible sites for particle acceleration inside current sheets are singled out, zooming promising regions with evolving plasmoids in selected cases.In order to quantify and characterize the role of resistivity, we then show important quantities like the electromagnetic energy content in the jet, the dissipated Ohmic power, and the non-ideal component of the electric field, which is no longer perpendicular to both the magnetic field and the velocity in the resistive case.We show that our combination of numerical algorithms is able to ensure robustness for the variety of the explored parameters, and highly accurate results, in terms of small-scale turbulent features and reconnection sites.All results are reported in Section 3, whereas Section 2 is devoted to the illustration of the system of equations solved and the numerical recipes employed, and conclusions and future applications will be discussed in Section 4. Numerical benchmarks, resolution tests and comparisons with other numerical schemes and with ideal simulations are reported in the appendices.

Equations and numerical methods
In this section, we briefly describe the fundamental equations of resistive relativistic MHD (Anile 2005;Komissarov 2007), first in covariant form and then specialized to Minkowskian flat space, splitting time and space components.We then shortly summarize the numerical algorithms employed for this work.

The system of Resistive Relativistic MHD (RRMHD)
The equations adopted in order to describe resistive relativistic plasmas, consist in the Maxwell equations: where F µν and F * µν are, respectively, the Faraday tensor and its dual, and J µ is the four-current density, coupled with a set of conservation laws for, respectively, the mass and (total) energymomentum: where ρ is the rest mass density, u µ is the four-velocity of the fluid and T µν is the total energy-momentum tensor, with its fluid (gas) and fields (em) components, where the latter is given by with g µν the metric tensor.The Faraday tensor (and its dual) can be conveniently written as functions of the four-velocity and the magnetic b µ and electric e µ fields measured in the fluid rest frame: where ϵ µνλκ represents the Levi-Civita pseudo-tensor.It is then convenient to express the two energy-momentum tensor contributions as2 : where h, ε and p are the gas specific enthalpy, energy density and pressure, respectively, measured in the fluid rest frame and related by ρh = ε+ p.Notice that in the non-resistive case of ideal MHD the condition of a vanishing electric field must hold in the frame comoving with the fluid, that is e µ = 0, and the expression for T µν em simplifies considerably.In numerical relativity, it is convenient to split time and space according to a so-called Eulerian observer (e.g.Del Zanna et al. 2007;Rezzolla & Zanotti 2013).Here we assume for simplicity a flat spacetime, and this duty will be easily accomplished.The fluid four-velocity is split as: where γ is the Lorentz factor and v is the three-dimensional velocity.The relation between the rest-frame (e µ , b µ ) and the standard laboratory (Eulerian) frame (E, B) electromagnetic field components is: In the resistive case e µ is not vanishing, and the simplest possible form of the (relativistic) Ohm's law assumes that it is proportional to the current density measured in the fluid rest frame j µ = J µ + (J ν u ν )u µ , that is: where η is the resistivity coefficient (which here, for the sake of simplicity, is assumed to be a scalar function, or even constant).
In the assumed split, the explicit expression of the four-current density for flat metric becomes: in which q represents the charge density measured in the laboratory frame.Notice that the spatial current density J contains a standard advection term (of q) and the resistive term, proportional to 1/η, that is expected to be large in high-Reynolds number astrophysical plasmas.Splitting the time and space components of the equations themselves, the evolutionary system of resistive relativistic MHD (RRMHD from now on) is, in vectorial form: respectively the equations for the conservation of mass, momentum, energy, and the two evolutionary Maxwell equations, with u em = (E 2 + B 2 )/2 the electromagnetic energy density.The conserved quantities are, respectively, the density in the laboratory frame D, the total momentum density m, the total energy density E, the magnetic B and electric E fields.The relation from primitive (ρ,v,p) to conserved fluid quantities is: while, contrary to non-relativistic MHD, an analytical relation from conserved to primitive variables is not possible.The set of RRMHD equations is completed by the non-evolutionary Maxwell's constraints: and by a suitable equation of state relating the thermodynamical variables, in particular providing h = h(ρ, p) (see next subsection).
Notice that in the ideal MHD case, when η → 0, the vanishing of e µ translates to the condition E = −v × B, the same as in classical MHD.Therefore, the field E is now a derived quantity, and the last evolution equation of the RRMHD set does not need to be solved.Moreover, in this case Eqs.(7) become: (13)

The Taub Equation of State (EoS)
The RRMHD system must be closed by an Equation of State (EoS), describing the thermodynamical properties of the gas through a relation between the (specific) enthalpy and the temperature function Θ = p/ρ.In numerical relativity hydro and MHD simulations the constant Γ−law for an ideal gas is typically assumed: where Γ is the specific heat ratio, whose value lies between 4/3 (ultra-relativistic plasma) and 5/3 (non-relativistic plasma).A major drawback of such approach is that the Taub's inequality (Taub 1948), which is required in order to gain consistency with the relativistic kinetic theory : may not be fulfilled.
From the theory of relativistic perfect gas (Synge 1957), a rigorous yet time-consuming relation (which includes modified Bessel functions) can be recovered.An alternative approach, which always ensures the Taub's relation (with the equal sign), is to assume (Taub EoS from now on): The above EoS is an approximation of the rigorous Synge's function h(Θ) (they differ for few percent values), it is much faster to compute numerically, and it retains the two physically important limiting values.It has been applied to the ideal relativistic MHD equations by Mignone & McKinney (2007) and later extended by Mizuno (2013) to resistive relativistic MHD applications.All the simulations presented in this paper, if not stated otherwise, involve the above Taub EoS.

The IMEX scheme and numerical methods
Let us start by discussing the time-stepping algorithm.The set of RRMHD equations can be rewritten in a more compact form: where U is the vector of the conserved variables, R contains the divergence of the fluxes F and the standard source terms S e , whereas the S includes the stiff terms, i.e. those proportional to 1/η ≫ 1 (the resistive part of the current density, in the equation for the electric field), requiring some special numerical treatment.Because of the stiffness of the RRMHD equations, purely explicit algorithms would require a drastically small timestep in order to encompass the time-scales typical of the nearly ideal regime.Several strategies have been developed in order to overcome such issue.For instance, Komissarov (2007); Takamoto & Inoue (2011); Mizuno (2013) employed a Strang splitting algorithm.However, as pointed in Ripperda et al. (2019a), due to the strong assumptions of such approach (i.e. that the magnetic field and the fluid velocity remain unaltered during the evolution of the electric field), the timestep required to achieve convergence again decreases dramatically in the regimes typical of astrophysical plasmas.
On the other hand, the Strong Stability Preserving (SPP) IMplicit-EXplicit (IMEX) Runge-Kutta (RK) schemes by Pareschi & Russo (2005) have shown great stability properties without additional assumptions, and are frequently adopted for both special and general relativistic resistive MHD simulations (Palenzuela et al. 2009;Del Zanna et al. 2016;Tomei et al. 2020).In this paper we adopt the IMEX-RK SSP3 (3,3,2) scheme to solve Eq. 17: where a = 1−1/ √ 2, which is second order accurate in time.Note that Eq. 18 holds for any general set of hyperbolic differential equations with stiff terms, with the coefficient derived from the Butcher's tableau (Pareschi & Russo 2005).The extension of the IMEX scheme to a generic EoS, and to the Taub one in particular, can be found in the next sub-section.
All the simulations have been performed by using the PLUTO code (Mignone et al. 2007(Mignone et al. , 2012)), adopting a fourthorder Piecewise Parabolic Method (PPM) (Mignone et al. 2005;Mignone 2014) for spatial reconstruction of variables at intercell faces, to compute numerical fluxes.For the sake of simplicity and ease of reproducibility of results, we choose a Harten, Lax, van Leer (HLL) (Harten et al. 1983) approximate Riemann solver, with local maximum characteristic velocities equal to the speed of light (recall that we are dealing with a hot, magnetized plasma with high flow velocities), since it proves to be very robust and no qualitative differences have been noticed by employing the more accurate estimates for the maximum eigenvalues (Leismann et al. 2005;Del Zanna et al. 2007).For similar reasons, we choose not to adopt a specific Maxwell Riemann solver for the induction equation (Mignone et al. 2019).The use of more diffusive Riemann solvers is expected to be compensated here by performing high-resolution runs and using high-order spatial reconstruction methods (Del Zanna et al. 2003, 2007).
As far as the methods for preserving the constraints in Maxwell equations, we adopt the divergence cleaning (GLM) method (Dedner et al. 2002;Mignone & Tzeferacos 2010) and we do not evolve the charge density, which is directly replaced by the divergence of the electric field (Bucciantini & Del Zanna 2013;Bugli et al. 2014).A comparison between the divergence cleaning method adopted in this paper and the Upwind Constrained Transport (UCT) for both the magnetic (Londrillo & Del Zanna 2004;Del Zanna et al. 2007) and electric field (Mignone et al. 2019) is shown in Appendix A, through a Kelvin-Helmholtz instability benchmark.These choices are again motivated by simplicity, ease of implementation and reproducibility (in different codes where not all methods may be available), given that from our benchmarks accuracy and robustness are preserved anyway.
In order to ensure stability in the more troublesome zones, we switch on the shock flattening procedures of the PLUTO code (Mignone & Bodo 2006).In presence of zones with negative pressure, we set the latter to a floor value of p f = 10 −5 while redefining conservative variables accordingly.

The IMEX scheme for a generic EoS
Here we derive, for the first time to our knowledge, the extension of the IMEX scheme for any EoS, in particular for the Taub equations of state.The implicit step of the IMEX scheme is solved by finding the roots of the equation: through a Newton-Broyden algorithm (see eq. 59 of Mignone et al. 2019), where D, m, E, and B are the conserved variables, known in the conversion step to the primitive variables.The explicit form of the Jacobian is: where ϵ ilm is the 3D Levi-Civita symbol.The electric field components and their derivatives can be easily calculated by using the prescription of Tomei et al. (2020) (the electric field does not depend on the EoS).
Once the electric field is computed (as function of the fourvelocity u), the fluid pressure can also be derived from the energy equation.If we define E g = E − u em , which in turn depends on u, we can obtain p = p(γ), hence p = p(u).In the ideal case the expression is: while, if the Taub EoS is employed, it can be derived by solving the quadratic equation: The pressure is needed to compute the temperature function Θ = p/ρ = pγ/D, hence the enthalpy h = h(Θ) is recovered, and it is possible to compute f(u).The derivatives of the specific enthalpy, to be inserted in the Jacobian J i j , can be computed as: independently on the EoS, where (24) The implementation of the coupling between the IMEX scheme and the EoS proposed has been tested through a set of numerical benchmarks which are shown in Appendix A.
From the computational point of view, a crucial point is the initial guess for the 4-velocity in the iteration scheme, which may not be trivial for more complex configurations of the RRMHD variables.Towards the ideal regime (e.g.∆t/η ≳ 10), the ideal guess allows for convergence in less than 5 iteration.On the other hand, in the more diffusive regimes, the ideal solution becomes brittle and therefore the solution at the previous timestep is taken.However, when the velocity is non-negligible in any of the three components, the convergence may not be reached with one of the two guesses mentioned previously.In such case, the Lorentz factor grows exponentially until it reaches unphysical values.To circumvent this scenario, we put a flag on the Lorentz factor.In case the upper value of γ = 1000, which in the simulations presented in the paper is a clear sign of lack of convergence of the implicit step, we reset the primitive variables to their value at the previous timestep, while the velocity is set to zero in all the components.By adopting this fix, the stability of the simulations has strongly improved.As a drawback, the number of iterations in the challenging cells has increased from ∼ 5 to ∼ 15.However, we point out that this increase in the number of iterations occurs only in few cells and at very few steps, producing negligible additional computational time.

Jet propagation and resistive models
Here we discuss the propagation of relativistic jets with finite conductivity.At first we briefly describe the numerical setup, then we investigate the effects of the resistivity.Consequently we introduce more consistent diffusivity models and we discuss the results obtained by employing non-constant physical resistivity.

Numerical setup
We consider a 2D axisymmetric setup in cylindrical coordinates (r, ϕ, z) with r ∈ [0, 25 r j ], z ∈ [0, 50 r j ] (invariance is assumed in the toroidal direction ϕ), and uniform grid resolution of 1200 × 2400 computational zones.The magnetized relativistic jet is injected into the domain from the lower boundary (z = 0) in the region r < r j , where r j = 1 is the characteristic length used for normalization, which corresponds to 48 grid cells in both r and z directions.
The jet is initialized with uniform density ρ j = 1 and vertical velocity v j = (0, 0, v j ), specified by the Lorentz factor γ j = 10 (i.e.v j ≃ 0.99).Its magnetic field structure is the same as in Mignone et al. (2009).It consists of a variable toroidal component, defined as: where r m = r j /2 is the radius where the toroidal field attains its maximum value γ j b m .This is determined by: where we set σ ϕ = 0.3.The value of the pressure p j is retrieved by imposing the Mach number: and here we choose M = 6 and Γ = 5/3, independently on the EoS employed.The thermal pressure distribution inside the jet can then be recovered assuming radial momentum equilibrium, and we find: hence the pressure has a maximum at the axis for r = 0, reaches its minimum value p j for r = r m , and remains constant beyond it.The poloidal magnetic field is instead purely vertical and constant everywhere: with σ z = 0.7.The two free magnetization parameters σ ϕ and σ z are actually defined through radial averages of quantities like 2b 2 /p, see Mignone et al. (2009) for their actual definitions.The electric field, unless stated otherwise, is initialized inside the jet to its ideal value E = −v×B, providing a purely radial component The external ambient medium is static and uniform, with density ρ a = 10 3 ρ j , gas pressure p j , purely vertical magnetic field B z (as in Leismann et al. 2005;Beckwith & Stone 2011), vanishing electric field.Hence, at the jet boundary r = r j , the quantities ρ, v z , B ϕ , and E r are discontinuous.
All the simulations are run until t = 300 (time will be expressed from now on in units r j /c).In order to easily distinguish the jet plasma from that coming from the ambient medium, we follow the evolution of an additional equation describing the advection of a tracer f , which vanishes in the ambient and is injected with the value f = 1 in the jet region.

Constant resistivity models
In order to investigate the impact of the resistivity on the jet propagation process, we performed three simulations with constant (in space and time) resistivity (low, η = 10 −6 , medium, η = 10 −3 , and high, η = 10 −2 ).We also ran a simulation (not shown here for the sake of simplicity) with the same numerical setup, numerical schemes and infinite conductivity as an additional validation.Differences between the ideal simulation and the simulation with the lowest resistivity value η = 10 −6 showed no qualitative difference.Therefore we expect the numerical diffusivity to lie between 10 −6 and 10 −3 , although we chose a lower value of eta as proof of code robustness.
In Fig. 1 we show maps, at the final time and for the three reference resistive simulations, of the rest mass density and of the inverse of the plasma−β, that is β −1 = B 2 /2p, the ratio of magnetic to thermal pressure.Together with the proper relativistic magnetization (defined as σ = B 2 /ρh, or σ = b 2 /ρh), the above classical definition of the inverse plasma-β is very often used as a proxy for magnetization in studies of the reconnection process in the non-relativistic and in the relativistic regime (Del Zanna et al. 2016;Ripperda et al. 2019b;Puzzoni et al. 2021).As a preliminary consideration we notice that the level of refinement (determined by the turbulent structures shown in the density and the plasma-β maps) of the low resistivity simulation, is comparable with the one obtained by (Mignone et al. 2009, see also Appendix B), who employed the less diffusive HLLD Riemann solver, linear reconstruction and a 3 rd -order time integration using a Runge-Kutta algorithm (Gottlieb et al. 2001).More accurate and less diffusive Riemann solvers (Miranda-Aranguren et al. 2018;Mignone et al. 2019) would be likely to suppress the numerical diffusivity and enhance even more the turbulent jet structure (when looking at the fluid variables); nevertheless, in presence of physical resistivity, the spatial scale of the magnetic turbulence should not be dictated by the numerical methods (see the resolution study of Appendix B).However, we confirm here that for relativistic MHD the HLL approximate Riemann solver combined to high-order spatial reconstruction (Del Zanna et al. 2003, 2007) is a robust and efficient compromise also for resistive simulations.
The most prominent differences are related to the electromagnetic field, which is strongly affected by the physical resistivity present in our simulations.We notice (see the zoomed panels of Fig. 1) a much steeper gradient in the magnetization at the interface (contact discontinuity) between the jet plasma and the shocked ambient medium (from 2 orders of magnitude in the low-η case to a smooth transition in the high-η case).Such a feature is not visible while looking at the density profile, and it is related to the suppression, due to the smearing of magnetic gradients, of the magnetic field turbulence caused by the interaction between the jet gas and the ambient medium.As a result, the jet magnetization is suppressed for the largest value of η.
For low and moderate values of the resistivity (i.e.η ≲ 10 −3 ), the jets show multiple zones with an average magnetization between 0.1 and 10, which is comparable to previous RMHD reconnection studies (Del Zanna et al. 2016;Ripperda et al. 2019b), suggesting that large-scale magnetic reconnection and plasmoids formation may occur in such locations (see the next sub-section).Note that in the cited works the plasma magnetization was set as an initial parameter, which is independent on the resistivity; here, despite the same initial physical conditions both at injection and in the ambient plasma for all our simulations, the final jet magnetization is strictly interrelated with the physical resistivity.
The impact of the resistivity reflects on the large-scale evolution of the field, as it can be shown by computing the injected electromagnetic energy, defined as (see Perucho et al. 2019) where f is the jet tracer.This is shown in Fig 2 as a function of time.The solid line represent the values of the injected energy normalized to the final value of the run with η = 10 −6 .In the very early evolutionary phases, the middle and low resistivity cases are almost indistinguishable.However, around t ∼ 35, the two cases start showing some differences since the jet head has already started interacting with the surrounding medium, slowing down the propagation and triggering the turbulence.Such difference is much more prominent while looking at the high diffusivity cases, in terms of both time scale and suppression of the electromagnetic energy.At time t ∼ 300 the ratio between the injected energy is ∼ 0.64 between the middle and the low resistivity and ∼ 0.15 between the high and the low resistivity cases.
The dashed lines of Fig. 2 represent the fraction of injected energy normalized respect to the final value of each simulation (the cyan solid and dashed lines overlap since their normalization is equivalent).We notice that the cases with higher resistivity reach a higher fraction of their final injected energy at much earlier times.Such differences can be also related to the contribution provided by the resistive field, which becomes less negligible for higher values of η (see a wider discussion in 3.2.2).This, in addition to the higher field diffusion (as shown by the maps of β −1 ), may contribute to explain this behavior.
On the other hand, the impact of the lower jet magnetization and higher resistive electric field for increasing values of η on the jet shape is much less evident (see Fig. 3).This is shown by the low and medium resistivity cases whose jet width (defined as the largest r−coordinate of the interface between the jet and the cocoon) differs from each other by ≲ 10%.On the other hand, the high resistivity case yields a slower (for most of the evolution) and wider jet.The impact of the resistivity on the jet shape becomes more relevant after the turbulence has fully developed within the jet.In particular, differences in the jet width or head (defined as the largest z−coordinate of the interface between the jet and the cocoon) are not seen (≲ 1%) until t ≳ 70.

The formation of current sheets
As we have seen, the most prominent differences among the propagating jets at various values of (constant) η are related to the electromagnetic field, which is strongly affected by its diffusion caused by the non-ideal processes occurring in our simulations.More specifically, the higher diffusion of the magnetic field leads to a smearing of the poloidal component, resulting in the suppression of the flow and magnetic vortexes typical of MHD turbulence (see the zoomed panels of Fig. 4).
Turbulence in plasmas is known to lead to the formation of current sheets and to reconnection events, and this is particularly important in relativistic magnetically-dominated plasmas where the release of energy could be huge.Relativistic magnetic  In order to more sistematically detect possible reconnection sites, we focus on the structure of the poloidal magnetic field and its curl, that is the non-relativistic current J ϕ = (∇ × B) ϕ , as often done in the literature, even by devising specific detectors based on this quantity (e.g.Zhdankin et al. 2013; Nurisso et al. 2023).Therefore, in Fig. 4 we show maps of the non-relativistic toroidal current, and of the poloidal field from which the former quantity has been computed.As a first consideration, we point out that current sheets and magnetic islands can also appear in simulations without an explicit physical resistivity (ideal MHD), due to the intrinsic and unavoidable diffusivity of the numerical schemes.Properties such as the width of the current sheets arising from the ideal simulations are strictly related to the grid resolution and the numerical schemes.Here we have verified that our ideal runs are very close to the case with η = 10 −6 , meaning that high-order reconstruction is able to maintain at the minimum level the numerical diffusivity (comparable with the ideal simulations), in spite of the use of the HLL approximate Riemann solver.
Because of the lack of turbulence caused by the physical resistivity in the most dissipative case with η = 10 −2 , the formation of current sheets is largely reduced compared with the other two cases.In particular, thicker current sheets are formed only in selected regions, e.g.near the axis, the jet head and between the jet  , where the poloidal magnetic field (color) and magnetic field lines (white lines) are plotted in order to highlight the reversal of the magnetic field and the formation of magnetic islands due to magnetic reconnection.material and the bow shock.Conversely, a lower resistivity leads to more ubiquitous and thinner current sheet regions.The strong dependence of the formation of current sheets with the higher magnetization, maintained only in the lower resistivity runs, is promising, given that astrophysical plasmas are expected to be quasi-ideal and favourable for the turbulent reconnection scenario.
In Figure 5 we have isolated four different regions from the case η = 10 −3 where the magnetic field shows an inversion of its polarity and becomes prone to the formation of relativistic current sheet (which are a primary ingredient for magnetic reconnection).Each region shows one or multiple inversions of the magnetic field lines due to the turbulence within the jet, and the consequent tearing instability leading to the plasmoids.Despite the lack of resolution required in order to properly resolve the complex internal structure of the magnetic islands, we are able to select several portions of each region where the magnetic field polarity is suitable for magnetic reconnection and plasmoid instability, allowing us to properly detect the X-and O-points in current sheets (as shown in the zoomed panels).All the panels selected (but the top right, number 2) show an average magnetization between 0.1 and 10 (see Fig. 1), since a lower magnetization may not lead to the formation of plasmoids (Ripperda et al. 2019b).Moreover, as shown in the right panel of Fig. 5, notice that the electric and magnetic contributions to the electromagnetic density are comparable in the jet region, mostly due to the high speed flows.Higher resolution simulations are expected to resolve the X-and O-points to a better accuracy, showing the regions inside current sheets where the electric field becomes stronger than the magnetic field.
Despite the strong jet magnetization in terms of magnetic energy, comparable to the thermal one (at least in the lower resistivity cases), when the interaction with the ambient medium starts, the resulting relativistic magnetization is strongly suppressed by the mixing, given that the ambient medium is a thousand times denser than the jet.The interface at the contact discontinuity is prone to the Kelvin-Helmholtz instability, due to the strong shear velocities, and this is the main trigger for this mixing.However, the same instability enhances the turbulence level and favours in turn magnetic reconnection, as shown by the presence of the Xpoints and magnetic islands in the zoomed panels also present in these regions.

The resistive electric field
In addition to the formation of current sheets and plasmoids, the resistive component of the electric field is also an indicator for sites where non-ideal processes may take place.Since this term can only appear if a physical explicit resistivity coefficient is present, magnetic reconnection in relativistic plasmas should be addressed through consistent numerical recipes, i.e. employing resistive RMHD schemes with explicit resistivity, as properly achieved here.
The presence of such non-ideal component of the electric field can be highlighted by displaying quantities like E • v and E • B, since they both vanish in ideal MHD (and RMHD).These are shown in Fig. 6, the parallel electric field component with respect to the velocity in the left panel block and the parallel electric field component with respect to the magnetic field in the right panel block, for the intermediate (left panels) and high (right panels) diffusivity models.Notice that those components vanish at the injection region by construction (we recall that an ideal electric field is assumed there, in the radial direction).In presence of a high resistivity, the resistive electric field starts to arise closer to the jet base (see Fig. 6).
Since the numerical resistivity is dominant in the lowest resistivity model, we do not show here the corresponding resistive electric field.Instead, we focus on the models whose physical diffusivity overcomes the numerical dissipation, since the resistive properties can be derived more consistently.
Because of the very turbulent structure of the magnetic field, the intermediate diffusivity case yields small layers with different polarity (see the zoomed panels of Figure 6).Due to the "poor" resolution, such layers are separated by a very limited number of computational cells.The most evident difference between the medium and the high resistivity cases is given by the the electric field orientation.As shown by the left and middle panel blocks of Fig. 6, the parallel component in the η = 10 −3 case shows, throughout the whole jet, an angle (given by the normalization of the scalar product between electric field and magnetic/velocity field) closer to 90 • than the high resistivity case.
We also notice that in both resistive regimes the resulting spatial structure of E•B is quite smoother with respect to E•v, especially for the η = 10 −2 case where larger scales are clearly produced.As discussed previously, this is not surprising, since the dissipation scale of the magnetic field is determined by the values of the physical resistivity, while the turbulent velocity structure is more affected by the intrinsic diffusion of the numerical algorithms.Note that a stronger relevance of the parallel component does not necessarily mean a stronger resistive electric field, since the electromagnetic energy depends on the diffusivity as well.
In addition to the different relevance of the parallel component, the case with high resistivity shows a smoother transition of the parallel component from the jet to the outer ambient medium (where it is absent by construction).In particular, the parallel electric field, in presence of higher resistivity, is diffused outwards more rapidly also in the region between the turbulent material and the bow shock (especially when looking at the parallel component with respect to the magnetic field).
A final consideration can be done on the impact of the finite conductivity on the jet dissipated power, whose definition in the fluid comoving frame is: Non-ideal MHD processes (e.g.resistivity) usually involve an energy conversion; in this case the magnetic energy is converted into thermal energy.As shown in Fig. 7, a higher resistivity leads to a more efficient energy conversion rate and heat generation through all the jet up to the bow shock.On the other hand, the case η = 10 −3 shows significant power dissipated near the jet head and, more generally, in the jet regions where magnetic reconnection is more likely to take place.Such behavior is consistent with the fact that the reconnection process involves a conversion from magnetic energy to kinetic, thermal and particles (on much smaller scales than the one typical of these simulations) energy.

Variable resistivity models
A constant resistivity profile, albeit a very simple and effective strategy, may not consistently model the differences between the jet and the ambient medium.In this section we compare two more consistent, non-constant in both time and space, resistivity models.The assumption of a resistivity confined in the jet is justified, as in Fendt & Čemeljić (2002), by the turbulent (and thus resistive) nature of the accretion disks.2020) in the context of mean-field dynamo in thick accretion tori.Moreover, turbulence is also supposed to be generated in the vicinity of black holes (Nathanail et al. 2022), leading to the formation of plasmoids for up to ∼ 40 gravitational radii.It is therefore natural to expect that the turbulent/resistive nature of the accretion flow is present also in the outflow at greater distance from the accreted matter.
A simple yet effective way to model a more consistent diffusivity profile is to relate it to the passive scalar tracer: Since the passive tracer has values between 0 (ambient) and 1 (jet), the resistivity will always be bound between 10 −6 and 10 −3 .Here the resistivity is mostly determined by the mixing of jet and ambient medium, while the electromagnetic components are only included from the full evolution of the RRMHD equations.
A more consistent resistivity profile is computed, similar to Fendt & Čemeljić (2002), by fixing the Lundquist number S , yielding: where v A = |B|/ ρh + |B| 2 is the Alfvén speed and L, a typical length scale of the system, is set to be equal to the injection radius r j = 1.Here, in contrast to Fendt & Čemeljić (2002), we consider the whole magnetic field (i.e.toroidal + poloidal), instead of only the poloidal field.We choose a value of S = 10 3 in order to reproduce a proper physical resistivity (i.e. higher than the numerical dissipation) in a suitable range to model astrophysical plasmas.

Magnetic field and magnetization
In Fig. 8, we show the toroidal current and the poloidal magnetic field of the constant low resistivity model (here used as a proxy for the numerical dissipation) compared with the two more consistent resistivity profiles (henceforth Tr-case and S-case).As a first consideration, we notice how the turbulent pattern of the poloidal magnetic field features the same level of turbulence (see Appendix B).Another similarity between all these runs is the formation of current sheets, highlighted by the toroidal component of the (classical) current (∇ × B) ϕ , through the entire jet domain.The jet head and axis, as well as the region between the jet and the bow shock, remain the locations where current sheets are more likely to be located, but the turbulent region (especially in the S-case) is equally promising as source of non-thermal particles due to magnetic reconnection.However, the injected energy (shown in Fig. 9), in both cases with non-uniform resistivity, saturates at lower values (∼ 75% of the low constant resistivity case).This different saturation is we notice a different saturation of such value.More specifically, already at t ∼ 30, the average jet resistivity has reached a value of 2.5 × 10 −4 for the S-case, which remains quite constant through the rest of the simulation.On the other hand, the average jet resistivity in the Tr-case starts with a value of ∼ 10 −3 as the jet is injected (t ≳ 0), since the jet and ambient matter have not started to interact yet; around t ∼ 100, the jet and the ambient medium have already interacted and the average jet resistivity slowly converges to 4.5 × 10 −4 until the end of the simulation.
This confirms the trend seen in Fig. 2, where a higher resistivity lead to a weaker injected energy.Moreover, it explains why the case S = 10 3 reaches the saturation of the injected energy at earlier time; since the resistivity decreases because of the mixing, the saturation value increases as the turbulence develops.
We also notice that the spatial resistivity profile shows clear differences between the different resistivity models.More specifically, the Tr-case shows a resistivity dominated by the jet at earlier time (see the left panel of Fig. 10), and even at later times a difference of almost 2 orders of magnitude are present between the resistivity at the lower regions and the jet head.On the other hand, the S-case shows a more uniform diffusivity distribution, since a stronger magnetic field triggers a higher resistivity which diffuses the magnetic field.This "cycle" of magnetic field and diffusivity has already converged at t ≲ 50, as shown in the left panels of Fig. 10.

The resistive electric field
As for the constant diffusivity models, the non-ideal electric field is a key component in order to investigate the role of the resistive features.In particular, in Fig. 11 we show the parallel component of the electric field with respect to the the magnetic field.
In both cases, the polarity of the parallel electric field near the boundary between jet and environment is the opposite with respect to the cases shown in Section 3.2.2(i.e.η = 10 −3 and η = 10 −2 ).This feature is not surprising, since in while investigating constant diffusivity models (see Section 3.2.2) we noticed such inversion in presence of a numerical resistivity (which is the case at the interface between jet and ambient medium for both variable diffusivity models).
However, the parallel electric field arises near the jet base, because of the higher physical resistivity at the jet base.The development of a strong parallel component in such region has been seen also in the η = 10 −3 and η = 10 −2 cases As a result, the jet location where the parallel field becomes more relevant is near the axes, in agreement with the η = 10 −3 case.
We notice that the spatial structure of E • B is significant less smooth than the cases with constant resistivity because of the lower values of the resistivity outside the very inner propagation  region.Moreover, the two models show a consistent difference in terms of spatial distribution of the parallel field.In particular, both cases show a relevant parallel component in the turbulent region, which is suppressed in the η = 10 −3 case.
The sharp gradient of the parallel field at the interface between jet and ambient is likely due to the numerical diffusivity, as in both models the value of η drops abruptly in the jet's environment.For the same reason, there is also a lack of diffused parallel field beyond that interface, as opposed to the models with constant resistivity (see the first panel of Fig. 6).

Discussion and conclusions
In this paper we presented the first systematic numerical study, in the regime of resistive relativistic magnetohydrodynamics (RRMHD), of astrophysical relativistic jets endowed with both poloidal and toroidal magnetic field components propagating through a uniform magnetized medium.We combined a high resolution grid in cylindrical coordinates with high-order numerical algorithms and realistic physical assumptions in order to achieve the best level of refinement, while preserving the robustness required to overcome the numerical issues often arising when simulating relativistic flows in a strongly magnetized regime.
The following summarizes our approaches and results: 1. we extended the IMEX scheme to the use of the Taub equation of state, in order to couple an accurate and robust evolution scheme for the electric field in non-ideal conditions and the proper treatment of the gas varying conditions, from the jet to the ambient.We have also augmented the stability of the IMEX scheme by improving the choice of the initial primitive variables during the implicit step; 2. we have initially investigated the role of a constant (in space and time) resistivity coefficient.More specifically, we found that a high resistivity coefficient, say η ∼ 10 −2 in typical code units, corresponding to roughly a Lundquist number of S ∼ 100 (normalized to the injection length and for an Alfvén speed close to the speed of light), suppresses almost completely the turbulence in the propagating jet and therefore the formation of relativistic current sheets.Conversely, a lower resistivity (at least η ≲ 10 −3 , or S ≳ 10 3 ) triggers the formation of multiple zones which could be prone (in terms of magnetization and field topology) to the plasmoid instability; 3. the electric field, in presence of a physical resistivity, develops a component parallel to the magnetic field which can be used as an indicator of non-ideal processes occurring at the turbuent scale, even in absence of magnetic reconnection events.We have found that different values of the resistivity strongly impact the arising of the resistive/parallel electric field, leading to a different spatial structure and dissipated power; 4. the combination of high-order numerical algorithms and high resolution enables us to find with increased accuracy the formation of X-and O-points within the jets, despite the lack of resolution required to resolve the particles acceleration scales.However, selected zones of the large scale jets can provide more consistent recipes in order to investigate the reconnection process at smaller scales; 5. we have also adopted and compared more consistent diffusivity models with a variable coefficient η, one based on a passively evolved scalar jet tracer, and another one fixing a constant value for the Lundquist number.Such models feature a more refined level of turbulence with respect to the case of intermediate constant diffusivity, consistent with the fact that the resistivity in the jet varies in the range 10 −5 < η < 10 −3 .However the injected energy saturates at different levels and at different times, because of the spatial and temporal variations of the resistivity.In our opinion the recipe based on a constant Lundquist number, that is η = v A L/S (L is a typical scale, here the injection width, and v A is the Alfvén speed), is to be preferred for both computational reasons (it is not subject to the numerical diffusion of an advected passive scalar) and physical ones (the resistivity is basically proportional to the local magnetization), so we recommend its use.
In conclusion, let us briefly discuss some important physical implications and possible lines of future development of the present work.
As illustrated by our simulations, the importance of addressing the large scale-separation between jet and fundamental plasma dynamics cannot be overstated.The use of high-accuracy numerical schemes is thus crucial to quantitatively assess the impact of turbulent diffusion in relativistic jets.Future work will therefore need to improve upon this aspect, e.g. by considering global high-order schemes (Berta et al. 2023).
Another debated issue is the possible role played by the resistive electric field, that is the non-ideal component proportional to the current density in the simplest Ohmic-type modelization, especially for the acceleration of the highest energy particles (see, e.g., Kowal et al. 2009;Medina-Torrejón et al. 2021;Sironi 2022;Puzzoni et al. 2022 and references therein).A possible scenario arising from PIC simulations is that the nonideal field (at X-points) plays an important role in non-thermal particle acceleration in microscopic current sheets as it starts the acceleration process.After this, the Fermi-type acceleration via the ideal electric field in the evolving plasmoids takes over, and most of the particle energy is gained (e.g.Guo et al. 2019).At the much larger macroscopic scale of astrophysical jets, further acceleration is likely to take place via the ideal electric field of MHD turbulence (e.g.Medina-Torrejón et al. 2023).This largescale Fermi acceleration can be studied using test particles or PIC-MHD methods (Mignone et al. 2018a;Vaidya et al. 2018;Medina-Torrejón et al. 2023), even in the context of propagating relativistic jets.
Due to the large separation in scales, the precise interplay between the microscopic ones (i.e. the current sheets typically studied via PIC simulations) and the large scales typical of MHD simulations of astrophysical sources is still very unclear.However, we would like to stress here that since the real resistivity is orders of magnitude lower than the values adopted in this paper (and the numerical diffusivity typical of any ideal RMHD simulation), the latter should not be strictly interpreted as the microscopic resistivity coming from the PIC simulations, but rather as an effective "sub-grid" resistivity, attempting to model the im-pact of unresolved and non-ideal kinetic processes on the larger scales.The effectiveness of such sub-grid models (which in several cases combine turbulent diffusion and mean-field dynamo action) has been positively tested in different astrophysical contexts by means of high-resolution studies (e.g.Liska et al. 2020;Tomei et al. 2020;Guilet et al. 2022;Reboul-Salze et al. 2022;Kiuchi et al. 2023).More refined resistivity models (see, e.g.Selvi et al. 2023) will be able to mimic with even better consistency the impact of microscopic kinetic process on the large scale scenario.
Other than AGN and blazar jets, the type of resistive simulations described here could also apply to different type of sources.A noticeable example is provided by the short-GRB jets arising from BNS merger multimessenger events (Pavan et al. 2021(Pavan et al. , 2023)), where the magnetic field is known to play a crucial role (Ciolfi 2020).For relativistic jets propagating in such a complex and structured magnetized environment, a finite conductivity is expected to be important and to lead to very different observational signatures, so this will certainly be a direction for future investigation.the one used for our models).All the quantities are normalized to the ones at the final time of the ideal simulation.As expected, all the runs show excellent agreement (i.e.discrepancies are ≲ 15%) in density (left panel) regardless of the specific value of η, while the magnetic field turbulence is affected by a sufficiently strong physical resistivity.In particular, we notice that the low-resistivity case, as well of the non-constant diffusivity runs show a turbulence which is comparable with the ideal simulation, while the runs with higher resistivity show a lower level of turbulence (respectively, 65% and 40% for the cases η = 10 −3 and η = 10 −2 ).
On the other hand, the turbulent scale required a sufficiently high resolution in order to be properly resolved.A "poor" resolution would lead to a higher numerical diffusivity which would strongly affect the turbulent scale and therefore the dynamical evolution of the jet.In order to prove that our resolution is enough to capture the magnetic dissipation scales dictated by the physical resistivity, we compared the normalized gradients for the case η = 10 −3 and η = 10 −2 by adopting different grid resolutions.As shown in Fig. B.2, the magnetic turbulence (right panel) shows very good agreement even at half (i.e.24 cells per r j ) of the reference resolution, while the poorest resolution model (i.e. 12 cells per r j ) shows significant qualitative differences.Concerning the density turbulence, despite slight differences between 36 and 48 cells per r j (≲ 10%), we find that our reference resolution is high enough to reach convergence.Lower grid resolutions, on the other hand, lead to a more suppressed fluid turbulence.

Fig. 2 .
Fig. 2. Jet injected electromagnetic energy.The solid lines represent the energy normalized to the final value of the run with η = 10 −6 , while the dashed lines represent the same quantities normalized to the final values of each corresponding simulation.

Fig. 3 .
Fig. 3. Jet head (circles) and width (crosses) positions for different values of the resistivity computed as functions of time.

Fig. 4 .
Fig. 4. Toroidal component of the magnetic field's curl (left panel blocks, representing a proxy for the toroidal current) and poloidal magnetic field (right panel blocks) for different values of the resistivity, in logarithmic scale at t = 300.The left, middle and right panel blocks show results computed with resistivity values of η = 10 −6 , η = 10 −3 and η = 10 −2 respectively.

Fig. 5 .
Fig.5.Current sheet forming in the turbulent jet region for the case η = 10 −3 at t = 300.The right panel shows the ratio between the electric and magnetic components of the electromagnetic energy.The 4 left panels show zoomed portions of the jet domain (see the rectangles in the right plot), where the poloidal magnetic field (color) and magnetic field lines (white lines) are plotted in order to highlight the reversal of the magnetic field and the formation of magnetic islands due to magnetic reconnection.
For instance, Qian et al. (2018); Vourellis et al. (2019) assumed a resistivity caused by the thin accretion disk rotation which triggers the magnetorotational instability.A similar assumption is made by Bugli et al. (2014); Tomei et al. (

Fig. 6 .
Fig. 6.Normalized parallel component of the electric field.For each panel group, E • B (left panel blocks) and E • v (right panel blocks) are shown for different values of the diffusivity at t = 300 in symmetric logarithmic scale.The parallel components are normalised respect to the corresponding quantities.The left and right panels of each block show results computed with resistivity values of η = 10 −3 and η = 10 −2 respectively.

Fig. 8 .
Fig. 8. Toroidal current (left panel blocks, computed from the curl of the magnetic field) and poloidal magnetic field (right panel blocks) for different values of the resistivity, in logarithmic scale.The left, middle and right panel blocks show results computed respectively with a constant resistivity η = 10 −6 , the Tr-case and the S-case.

Fig. 9 .
Fig. 9. Jet injected electromagnetic energy for different resistivity models at t = 300 in logarithmic scale.The solid lines represent the energy normalized to the final value of the run with constant (in space and time) resistivity η = 10 −6 , while the dashed lines represent the same quantities normalized to the final values of each corresponding simulation.

Fig. 10 .
Fig. 10.Resistivity computed at different times (t = 100, 200, 300 respectively in the left, middle and right panel blocks) with the two different variable resistivity model (respectively, Tr-case on the left panels and S-case on the right panels).

Fig. 11 .
Fig. 11.Parallel (normalized, as in Fig. 6 component of the electric field respect to the magnetic field for the Tr-case (left panel) and the Scase (right panel) resistivity models at t = 300 in symmetric logarithmic scale.

Fig. B. 1 .
Fig. B.1.Effective turbulence resolution (computed as ⟨|∇Q|/Q⟩) of density (left panel) and magnetic energy (right panel) computed for different resistivity models and the ideal simulation using the same setup as Mignone et al. (2009).