Amplification and Dissipation of Magnetic Fields in Accreting Compact Objects

Magnetic fields play a crucial role in shaping the dynamics of accreting compact objects. Whether we consider the formation of a proto-neutron star during the gravitational collapse of a massive star or the accretion disk around a black hole after a compact binary merger, a key process that remains challenging to include in large-scale simulations is the amplification and dissipation of magnetic fields driven by turbulent fluid motions. Despite the enormous increase in computational power currently available, the large separation between all the relevant spatial and temporal scales still poses severe limits to what can be achieved with ideal fluid simulations. One way to tackle such issue is to rely on sub-grid models, which however need to be appropriately tuned in light of models probing the small-scale dynamics. In this work we present the current state-of-the-art of dynamo models in proto-neutron stars, which aim at describing the amplification of magnetar-like magnetic fields during the gravitational collapse of a massive star. We also review some of the works from the past few years that included turbulent dynamos in accretion disks around a black holes, relying on a mean-field formalism. Finally, we will present a recent study on polar jets with explicit turbulent resistivity which showcases the importance of employing highly accurate numerical schemes.


Introduction
Accretion onto compact objects is one of the most efficient mechanisms in the Universe to convert gravitational energy into other forms of energy (either kinetic, magnetic or thermal), energetic particles, and gravitational waves.The efficiency for converting accretion power to luminous output can be as high as ∼ 42% in the case of a maximally rotating black hole (as opposed to the ∼ 0.7% efficiency of mass to energy conversion during hydrogen fusion in stellar cores) which motivates the central role of accreting black holes and neutron stars in the modeling of astrophysical high-energy sources such as X-ray binaries, active galactic nuclei, and gamma-ray bursts.Observational evidences for non-thermal spectral components and significant degrees of polarization are found in virtually all high-energy sources, and are usually interpreted as smoking guns for particle acceleration processes and large-scale magnetic fields, respectively.The role of magnetic fields is particularly crucial for accreting compact objects.They enable the accretion of matter via magnetohydrodynamic (MHD) turbulence produced by the onset of the Magneto-Rotational Instability (MRI, [1]) and when distributed over large-scales they can launch relativistic outflows.Moreover, magnetic fields can directly energize the outflow through reconnection events, where particles are accelerated at the expenses of magnetic energy due to the change in magnetic field topology [2].
The numerical modeling of accretion and ejection flows around compact objects remains a challenging endeavour to this day, given the large scale separation (both in space and time) that divides the local plasma dynamics, the central engine, and the emission sites of these sources A self-consistent description of the non-collisional nature of astrophysical plasmas (such as in PIC simulations) allows to capture fundamental kinetic and dissipative effects.However, the elevated computational costs that come with it limits this approach to relatively short temporal and spatial scales compared to those of typical astrophysical compact objects.At the same time, fluid approximations such as MHD can probe the accretion flow dynamics and the launch of relativistic jets, but they generally ignore the microphysics behind the dissipative and kinetic properties of the plasma.
There are several approaches that can be adopted to face these challenges and try to bridge the gaps that separate all the relevant scales of the problem.One might attempt to focus on different parts of the accreting systems, resolving smaller and larger scales with dedicated numerical simulations.The connection between the different scales can then be established either with suitable boundary or initial conditions (such as in the case of jet launching models not including the central engine ( [3], [4]) or with a well-calibrated sub-grid model (e.g. to include unresolved turbulent dynamics of the plasma [5]- [7]).Hybrid MHD-PIC methods are another useful tool to describe the acceleration of non-thermal particles embedded in a thermal fluid (such as in the case of cosmic rays [8] or electrons accelerated at collisionless shocks [9]).Finally, the advent of GPU-accelerated astrophysical codes [10]- [12] has paved the way to numerical studies with unprecedented resolutions [13]- [15] that probed dynamical regimes so far unattainable.
Whether we consider a still forming proto-neutron star (PNS) during the onset of a corecollapse supernova (CCSN) or the accretion disk around a supermassive black hole, the fundamental ingredients describing the amplification and dissipation of magnetic fields and their dynamical impact on the accretion/ejection dynamics are remarkably similar.In this work we review some of the most recent developments in the numerical modeling of magnetized accreting compact objects, which involve some of the aforementioned strategies to tackle the scale-separation problem.In §2 we will present results from studies focusing on PNS dynamo studies and CCSN large-scale simulations, i.e. the early stages of the formation of stellar compact objects.We will then show in §3 current state-of-the-art simulations of accreting black holes and relativistic jets that include non-ideal effects to model dynamo and diffusive processes in the plasma.Finally, in §4 we will present a new 4 th -order accurate finite volume method for classical and relativistic MHD that can be implemented to boost the effective resolution of numerical models and increase their efficiency.
2. Core-collapse supernovae CCSN are the result of the gravitational collapse of massive stars at the end of their life, and represent one the main channels for the formation of stellar compact objects.Once nuclear densities are reached at the core of the star, the equation of state (EoS) suddenly stiffens and a shock wave starts propagating from a central newly formed PNS.Due to the ram pressure of the infalling stellar material and the loss of energy via photo-dissociation of Fegroup elements [16], the shock stalls and cannot further propagate unless some other process Figure 1.Left panel: 3D rendering of the magnetic field produced by the convective dynamo in the strong regime [22].Blue (red) isosurfaces represent the radial velocity, i.e. downflows (outflows) of the order of 10 8 cm/s.Right panel: magnetic field lines within the MRI unstable zone of the PNS from a numerical simulation in the incompressible MHD regime [23].Blue and red lines indicate, respectively, low and high intensity magnetic fields.provides a suitable amount of energy.In the vast majority of CCSN events the shock is revived by the so-called neutrino-heating mechanism: a small fraction of the flux of neutrinos emitted from the PNS surface (which overall amount to ∼ 99% of the total gravitational binding energy lost by the explosion) can be absorbed behind the shock, stirring convection and aiding the outward propagation of the shock wave [16].However, outstanding stellar explosions such as superluminous supernovae [17] and Type Ic-BL [18] cannot be explained by neutrinos alone, and need to consider an additional source of energy.A common suggestion is that stellar progenitor with fast rotation can produce, in presence of strong large-scale magnetic fields, a magneto-rotational explosion: the rotational energy within the stellar core is extracted either via magnetic braking (for a proto-neutron star) or a Blandord-Znajek mechanism (for a black hole with and accretion disk), launching energetic bipolar outflows [19]- [21].

PNS dynamos
The exact origin of the magnetic fields required to obtain magneto-rotational CCSN remains unclear, as they could be produced within the stellar progenitor in the pre-collapse phase ( [24], [25]), amplified during stellar mergers [26] or directly in-situ within the PNS during the explosion via dynamo processes [27]- [29].In a series of recent studies we explored the latter scenario using the pseudo-spectral code MagIC ( [30], [31]) to focus on specific layers of the PNS where magnetic field amplification can occur through the action of mainly two physical processes.The first one is the convection triggered in the inner regions of the PNS (10 ≲ r ≲ 25 km) predominantly by negative lepton gradients due to persistent losses of electron neutrinos ν e [32] about ∼ 50 ms after shock formation.When rotation in the PNS is sufficiently fast, the Coriolis and Lorentz forces balance each other, thus producing a magnetostrophic balance that leads to the saturation of a dynamo magnetic field close to magnetar-like values (up to ∼ 10 16 G) [22].The resulting field is dominated by the axisymmetric toroidal component (left panel of Fig. 1) and non-axisymmetric structures, and can lead to a significant emission of GW over the spans of several seconds [33].
Dynamo action can also take place in the stably stratified layers of the PNS that surround its convective region (25 ≲ r ≲ 50 km) in presence of differential rotation, as the MRI can develop and further amplify the magnetic fields [34].Starting from small-scale magnetic seeds (but with intensity sufficiently high so as to properly resolve the MRI) the onset of the instability leads to the amplification of large-scale magnetic fields [35], which are also dominated by the toroidal component (right panel of Fig. 1).Moreover, the low multipolar (i.e.large-scale) components of the field are generally misaligned with respect to the PNS rotational axis.When the density stratification of the PNS is taken into account, the amplification process assumes characteristics typical of an αΩ-dynamo, with magnetic field reversals and periodic patterns that give rise to butterfly-diagrams [23].

Magneto-rotational explosions and multi-messenger signals
Regardless of the specific process at the core of PNS dynamos, the large-scale magnetic fields that they amplify tend to be dominated by non-axisymmetric structures and multipolar components beyond a simple magnetic dipole.To self-consistently evaluate the impact of such magnetic fields on the explosion mechanism and the ejecta launch would require to cover a range of spatial and temporal scales spanning at least 5 orders of magnitude (from a few tens of meters to thousands of km).A more approximated (but also affordable) approach is to perform CCSN simulations with initial magnetic fields having a topology more complex than a standard aligned dipolar field (as usually assumed in magneto-rotational explosion models) to qualitatively capture the possible effects of PNS dynamos on the explosion mechanism [37].We considered the stellar progenitor 35OC [38] with its original fast rotation rate and superimposed different magnetic field configurations: aligned dipole (model L1-0), aligned quadrupole (L2-0A and L2-0B), equatorial dipole (L1-90), and without any magnetic field (H).The strength of the magnetic field in the iron core is set to ∼ 10 12 G, which is a typical value employed in magneto-rotational explosion models [21], [39].We performed our 3D relativistic MHD simulations with the AENUS-ALCAR code [40], which includes nuclear EoS and a multi-dimensional M1-neutrino transport scheme to include all the relevant interactions between matter and neutrinos.The two-moment closure approximation we use to describe the neutrino dynamics is adopted by several other research groups [41]- [43], as it allows to capture the multidimensional effects of neutrino-matter interactions at a significantly lower computational cost w.r.t. the full-Boltzmann approach.However, it should be noted that this approximation can become inaccurate when transitioning between the optically thick and thin regimes if the structure of the flow is spatially complex (as in the PSN and the post-shock region; [44], [45]) and lead to quantitatively different results [46], [47].We include gravity via a simple Newtonian treatment with general relativistic corrections to the gravitational potential [48] (which captures the main impact that general relativity has on the compression of the PNS) and compute GW strains by means of the quadrupole formula described in [49].As shown in Fig. 2, configurations with higher multipoles (L2-0) or inclined magnetic axis (L1-90) produce less energetic explosions and less collimated outflows.Moreover, magnetic fields play a crucial role in quenching the onset of the low-T /|W | instability, i.e. a hydrodynamic co-rotational instability that can affect a PNS with very fast rotation [50].While in the hydrodynamic case we see around 200 ms post-bounce the emergence of large-scale spiral structures that protrude from the inner regions of the PNS through its surface (left panel of Fig. 3), in none of the magnetized models we observe such dynamics (right panel).This is a direct consequence of the magnetic transport of angular momentum, which is quite efficient at flattening the rotation profile within the PNS, thus stabilizing the PNS against the onset of the low-T /|W | instability.
This qualitative deviation between magnetized and unmagnetized models reflects also on the properties of the GW and neutrino emission during the explosion [51].The onset of the low-T /|W | instability corresponds to a burst of GW that lasts for several hundreds of seconds, which has a distinctive spectral signature in the frequency domain that appears around ∼ 200 ms p.b. and increases frequency with time due to the compression of the PNS (left column of Fig. 4).On the other hand, no such strong emission is observed in the magnetized models, which produce a GW signal with a rather broad spectrum that goes beyond the rotational frequencies within the PNS (white curves), as they have been narrowed by the outer transport of angular momentum.The neutrino emission is also deeply affected by the action of magnetic fields.All magnetized models have neutrino lightcurves along the equatorial plane which are systematically lower than the hydrodynamic case.This is due to the more oblate shape of the PNS, which is induced once again by the Maxwell stresses acting in the PNS and leads to the emission along the equatorial plane of neutrinos with a lower mean energy (since they decouple from the matter at larger distances).Moreover, the time-frequency diagrams of such lightcurves show features that, in  the hydrodynamic case, correlate quite well with the onset of the low-T /|W | instability and the its associated GW emission (see Fig. 12 of [51]), showcasing the important synergy between GW and neutrino detectors.

Accretion disks and jets
We now move to a later stage in the activity of compact objects, namely the case of a black hole surrounded by an accretion disk.Since at least two decades general relativistic magnetohydrodynamics (GRMHD) provide the fundamental framework for the study of such systems.The scientific community has independently developed several ideal GRMHD codes that have proved to be able to consistently produce models in quantitative agreement which each other [52].Nonetheless, in recent years there have been several efforts to go beyond the ideal GRMHD formulation by introducing kinetic dissipative effects [53], [54], electron thermodynamics [55], viscosity/anisotropic pressure [56], collisionless moment closures [57], magnetic resistivity [58]- [60], and mean-field dynamo processes [61], [62].In the following we focus on the last two effects, which are examples of sub-grid models aim at including physical processes that occur at unresolved spatial scales.

Resistive (G)RMHD
The introduction of an explicit resistivity η allows to model the physical dissipation of magnetic fields and avoid to be dominated by the numerical diffusivity (which is always present in any grid-based code).This is particularly relevant for highly magnetized environments such as the corona of an accretion disk or a relativistic jet close to its central engine, since the formation of current sheets and reconnection events play a crucial role in particle acceleration and magnetic dissipation [2], [64].However, the relativistic regime poses an additional challenge with respect to the classical case, since the displacement current cannot, in principle, be neglected in Ampère's law.Moreover, the electric field E cannot simply be expressed in terms of the magnetic field B and the fluid velocity v as in the ideal GRMHD case (where E = −v × B), but it is instead coupled to the current density J via the resistive Ohm's law [58] where Γ is the fluid's Lorentz factor and q is the electric charge density The combination of these two aspects promotes Ampère's law to a new evolution equation for the electric field  [58]: where γ is the determinant of the 3D spatial metric tensor, α is the lapse function and β is the shift vector.Since Eq. 2 can become stiff for low values of η, a common strategy to ensure stability to the integration of the resistive (G)RMHD equations is to employ IMEX schemes [65], which are designed to treat implicitly the stiff equations and explicitly the rest of the system.We recently used this framework (specialized to a flat space-time) to perform, for the first time, axisymmetric simulations of relativistic magnetized jets with a finite physical resistivity propagating through a magnetized ambient medium [63] with the resistive RMHD module of the PLUTO code.The resistivity we considered in this study does not represent the Ohmic dissipation within the plasma (which is generally extremely small and occurring at microscopic scales) but rather a turbulent resistivity, which tracks the dissipation of magnetic fields due to unresolved turbulent motions.While such effects are generally dependent on the dynamic and thermodynamic state of the plasma, for simplicity we first treated η as a completely free parameter.We tested the impact of a constant η = 10 −6 , 10 −3 , 10 −2 across the domain (in code units, i.e. in units of 4πLc with L = 10 −6 pc) and verified that for a higher resistivity the jet injected less electromagnetic energy and the magnetic field was distributed over larger spatial scales.In all cases, the bulk of the dissipation occurred in the jet's spine, close to the symmetry axis (left panel of Fig. 5).We also tested two different profiles for a non-constant resistivity, by either: 1) modulating its value between 10 −6 and 10 −3 using a passive tracer for the jet; 2) imposing a constant Lundquist number S = v A L/η = 10 3 with v A the Alfvèn speed and L the characteristic length scale of the system which we identified with the injection radius r j = 1.The resultant spatial distributions for η (right panel of Fig. 5) differ significantly from one another, with the tracer prescription leading to a slightly stronger mean value for η over the course of the simulation.

Mean-field dynamos in accretion disks
A further generalization of Eq. 1 can include the effects of a mean-field dynamo action, which models unresolved small-scale fluctuations in velocity and magnetic field as an additional largescale electromotive force in the fluid comoving frame that is proportional to the magnetic field [61].Such an assumption leads to a slightly modified Ohm's law [61] where the parameter ξ is the dynamo parameter defined as −α dyn , i.e. the opposite of the parameter after which the α-effect of αΩ-dynamos is named.On the other hand, all the rest of the relevant equations coincide with the purely resistive case.With tested this prescription by studying both kinematic [67] and fully dynamic [5] phases of an αΩ-dynamo in thick accretion disks around black holes.An arbitrarily weak magnetic field seed can be efficiently amplified by such mechanism (even in axisymmetry) and the process self-regulates itself through the modulation of the mass accretion rate onto the black hole, without any explicit quenching mechanism.The resulting magnetic and mass fluxes and the black hole event horizon are consistent with those retrieved by 3D ideal GRMHD models with higher resolution, validating its implementation [68].Moreover, when applied to 3D models the relativistic αΩ-dynamo induces large-scale non-axisymmetric structures [66] that can sometimes be found in ideal 3D models (Fig. 6).The emergence of αΩ-dynamos in relativistic disks has been recently obtained by very high-resolution 3D models [69], which furthermore demonstrates the potential of a relativistic dynamo closure to study more efficiently accretion flows around black holes.

A new 4 th -order FV scheme
We finally present another recent development that aims at increasing the efficiency of numerical simulations of astrophysical plasmas and therefore mitigating the difficulties of studying systems with a large scale-separation.We developed a new genuinely 4 th -order accurate FV method based on pointwise reconstructions that can significantly increase the effective resolution of a typical astrophysical simulation and thus reduce its computational cost [70].Consider the following FV discretization of a generic conservation law over a uniform Cartesian grid with spacing (∆x, ∆y, ∆z) where ⟨U ⟩ c is the cell volume average of the conservative quantity U , Fx f is the surface average of the x-component of the flux over a face perpendicular to the x-direction, Fx f −êx is the surface average of the same component over the opposite side of the cell, and similarly for the other flux terms along the y-and z-directions.For schemes that are less than 4 th -order accurate, one can identify the volume and surface averages in Eq. 4 with the point values of U and F evaluated at the cell center and at its sides centers, respectively, which effectively makes low-order FV and finite differences (FD) formulations equivalent to each other.However, this is no longer true for schemes that aim at achieving higher orders of accuracy, which then need to calculate suitable high-order approximations of both sides of Eq. 4. Starting from a high-order cell volume average of the conserved quantity U , we adopt the following general strategy (see [70] for a complete description): expression where ∆() is the Laplacian operator introduced by [71].• From the ensemble of cell centered conservative variables {U c } we retrieve the primitive variables {V c }. • We apply a pointwise reconstruction to obtain the primitive variables at the cell's interfaces and compute the upwind fluxes.• We convert the point-valued fluxes into high-order surface averages using where ∆ x ⊥ is a transverse Laplacian operator along the x-direction (and similarly for the other components).
• We integrate Eq. 4 in time using a 4 th -order accurate Runge-Kutta scheme.
As shown in Fig. 7, this scheme exhibits a clear 4 th -order scaling and greatly improves the solution accuracy even at modest resolutions.Its numerical dissipation rate (right panel) is more than 10 times lower with respect to a 2 nd -order scheme even at very low resolutions (∼ 8 points), which demonstrate its potential at allowing a more realistic description of magnetic diffusion when applied to a resistive framework.

Conclusion
We presented a series of recent results from numerical simulations aimed at studying the dynamics of magnetic fields in different kinds of accreting compact objects, such as protoneutron stars in CCSN and black holes surrounded by an accretion disk.Despite the different scientific context, microphysics, and numerical tools required, all these studies demonstrate how an accurate description of the amplification and dissipation of magnetic fields can significantly impact the accretion process, the launch of polar outflows, and the acceleration of particles.Several strategies can be implemented to overcome the challenges posed by the large scaleseparation between all the relevant physical processes involved, from multi-scale studies to subgrid models, from the use of GPU-accelerated HPC code to the adoption of more accurate and higher order numerical schemes.

15thFigure 2 .
Figure 2. Volume rendering of the specific entropy for CCSN models with different initial magnetic field configurations: no magnetic field (top left panel), aligned dipole (top right), equatorial dipole (lower left panel), and aligned quadrupole (lower right panel) at t = 410 ms p.b.[36]

Figure 4 .
Figure 4. Top row: Gravitational wave strains (+ and × polarizations) for the hydrodynamic CCSN model (left) and the one with an aligned quadrupolar magnetic field (right)[51].Bottom row: time-frequency spectrograms of the characteristic strain h char relative to the the strains above.The white lines track twice the rotational frequency of the PNS surface and center[51].

Figure 6 .
Figure 6.Toroidal magnetic field distribution in the equatorial (left panel) and meridional (right) plane for a GRMHD model of accretion disk around a Kerr black hole with spin parameter a = 0.9375 and mean-field dynamo action [66].The xy−axes and the magnetic field are express in code units, i.e. in units of the black hole gravitational radius R G = GM BH /c 2 and 4π/G 3 c 4 /M BH , respectively.

Figure 7 .
Figure 7. Left panel: convergence scaling computed for 3D Alfvén waves propagating along a diagonal for a 2 nd -order (triangles) and 4 th -order (squares and circles) FV schemes[70].Right panel: scaling with resolution of the estimated numerical dissipation from the aforementioned test problem.