Non-ideal magnetohydrodynamics on a moving mesh II: Hall effect

In this work we extend the non-ideal magnetohydrodynamics (MHD) solver in the moving mesh code AREPO to include the Hall effect. The core of our algorithm is based on an estimation of the magnetic field gradients by a least-square reconstruction on the unstructured mesh, which we also used in the companion paper for Ohmic and ambipolar diffusion. In an extensive study of simulations of a magnetic shock, we show that without additional magnetic diffusion our algorithm for the Hall effect becomes unstable at high resolution. We can however stabilise it by artificially increasing the Ohmic resistivity, $\eta_{\rm OR}$, so that it satisfies the condition $\eta_{\rm OR} \geq \eta_{\rm H} /5$, where $\eta_{\rm H}$ is the Hall diffusion coefficient. Adopting this solution we find second order convergence for the C-shock and are also able to accurately reproduce the dispersion relation of the whistler waves. As a first application of the new scheme, we simulate the collapse of a magnetised cloud with constant Hall parameter $\eta_{\rm H}$ and show that, depending on the sign of $\eta_{\rm H}$, the magnetic braking can either be weakened or strengthened by the Hall effect. The quasi-Lagrangian nature of the moving mesh method used here automatically increases the resolution in the forming core, making it well suited for more realistic studies with non-constant magnetic diffusivities in the future.


INTRODUCTION
Magnetic fields play an important role at different scales in the Universe.They can interact with a conducting fluid, decisively influence the evolution of the fluid (Ferriere 2001;Cox 2005) and at the same time amplify the field through dynamo effects.
A very simple approach to combine magnetic field with the theory of fluid dynamics is to use the ideal MHD equations.The latter extend the hydrodynamic Euler equations by adding the Lorentz force to the momentum equation and by introducing the induction equation, which describes the temporal evolution of the magnetic field.Because of its simplicity ideal MHD is often used in numerical studies, for example in simulations of large cosmological boxes (e.g.Marinacci et al. 2015;Marinacci & Vogelsberger 2016;Dolag et al. 2016;Marinacci et al. 2018b), galaxy clusters (e.g.Dolag et al. 1999Dolag et al. , 2002)), the magnetic field within individual galaxies (e. e.g.Pakmor et al. 2014Pakmor et al. , 2017)), the interstellar medium (e.g.Kim & Ostriker 2018;Simpson et al. 2023), molecular clouds (e.g.Grudić et al. 2021) and the collapse of protostellar cores (e.g.Mellon & Li 2008).
The ideal MHD equations contain however several assumptions that can be violated, in which case modifications to the underlying equations are required.This is the case, for example, in weakly ionised regions where collisions between neutral and charged particles become rare enough to allow drift velocities between different species.In the low frequency limit one can still derive single fluid equations describing the evolution of the bulk fluid (for a detailed discussion see Pandey & Wardle 2008).If we neglect the pressure terms in the induction equation, we find three new terms in the induction ★ E-mail: ozier@mpa-garching.mpg.deequation: The resistive Ohmic and ambipolar diffusion and the Hall effect.While the first two terms arise from collisions of charged particles that allow the magnetic field to diffuse, the Hall term becomes important in situations where ions are decoupled from the magnetic field while electrons are still coupled to it.In weakly ionised systems, decoupling can occur because collisions between neutrals decouple ions from the magnetic field more easily than electrons.Decoupling can also occur in the high frequency limit, when ions cannot follow magnetic fluctuations due to their lower cyclotron frequency compared to electrons.For a full discussion of both regimes see Pandey & Wardle (2008).Unlike the other two effects, which only damp MHD waves, the Hall effect introduces two new waves, the whistler (electron-cyclotron) and the ion-cyclotron wave.
A prototypical example of a system in which all three non-ideal MHD effects play a significant role is the collapse of a protostellar core and the subsequent evolution of the weakly ionised protostellar disk.In the ideal MHD approximation, magnetic fields can efficiently transport angular momentum outwards during collapse and prevent disk formation in this so-called magnetic braking catastrophe (see e.g.Allen et al. 2003).Non-ideal effects can weaken the magnetic field and reduce the efficiency of angular momentum transport, which supports the formation of rotationally supported discs (Wurster et al. 2016).For a recent review of the influence of magnetic fields on the formation of protostellar disks, see Wurster & Li (2018).
In the ideal MHD approximation, the newly formed disk can become unstable to the magnetorotational instability (MRI, Velikhov 1959;Chandrasekhar 1960;Fricke 1969;Balbus & Hawley 1991), which creates turbulence that can act as an effective viscosity, allowing the outward transport of angular momentum and thus accretion onto the central object (Shakura & Sunyaev 1973;Lynden-Bell & Pringle 1974).The MRI can be suppressed by the diffusive non-ideal MHD effects, but the Hall effect can also influence its linear properties, leading to the Hall-shear instability (HSI, Kunz 2008).Outflows from protostellar disks, which can remove angular momentum, can also be affected by non-ideal MHD effects (e.g.Bai 2014).Unlike the other two non-ideal effects, the sign of the Hall effect can change, resulting in a bimodal behaviour depending on whether the axes of rotation and the vertical magnetic field are (anti-)aligned.Depending on this condition, magnetic braking can be enhanced or weakened (Wurster et al. 2021).
Star formation is highly non-linear and arises from the combination of several different physical processes, such as gravity, fluid dynamics or magnetic fields, that all play an important role.Computer simulations are therefore a powerful tool for analysing such systems.Because of their importance, the three non-ideal MHD effects discussed above have been implemented in Eulerian codes such as PLUTO (Lesur et al. 2014), ATHENA (Bai 2014), ATHENA++ (Stone et al. 2020), RAMSES (Masson et al. 2012;Marchand et al. 2018), ZeusTW (Li et al. 2011) and MANCHA3D (González-Morales et al. 2018), as well as Lagrangian particle methods such as the SPH code of Tsukamoto et al. (2017), PHANTOM (Price et al. 2018) and GIZMO (Hopkins 2017).Eulerian codes typically suffer from larger advection in regions with high bulk velocities whereas the popular Lagrangian method called smoothed particle hydrodynamics (SPH) produces noisier results with lower accuracy but offers a densityadaptive resolution.
A relatively new approach that attempts to combine the advantages of a Galilei-invariant Lagrangian method with the high accuracy of the finite volume methods typically used in Eulerian codes is the moving mesh method as implemented in the AREPO code (Springel 2010;Weinberger et al. 2020).It solves the fluid equations on an unstructured, evolving grid, which makes it particularly interesting for the analysis of systems with large density gradients, due to its adaptive cell sizes, and systems with large bulk velocities, such as rotationally supported disks.As we have shown in Zier & Springel (2022b) and Zier & Springel (2023) the moving-mesh approach is indeed able to resolve the MRI as well as the gravitational instability (GI) with high accuracy in the shearing box approximation (Hill 1878;Goldreich & Lynden-Bell 1965;Zier & Springel 2022a).
In Zier et al. (2023) (Paper I from now on) we extended the AREPO code to include resistive Ohmic and ambipolar diffusion.The module is based on a least-squares fit to obtain accurate estimates of the magnetic field gradients and shows comparable accuracy to the stateof-the-art code ATHENA++, which uses a Eulerian method.The aim of this paper is to add the Hall effect to our implementation, which is numerically quite challenging as it is known to easily destabilise numerical algorithms that use second-order accurate time integration as in AREPO.However, we show that the algorithm can be stabilised by locally increasing the magnetic diffusion in a judicious fashion.
The paper is structured as follows.In Section 2 we introduce the non-ideal MHD equations and derive the linear dispersion relation of the new MHD waves introduced by the Hall effect.In Section 3 we introduce the moving mesh method and discuss the implementation of the term describing the Hall effect.For discretization we use the same method that we have already used for Ohmic and ambipolar diffusion.In Section 4 we show by simulating a magnetic C-shock that our scheme becomes unstable when the Hall effect is strongly dominant over the magnetic diffusion originating from explicit diffusion or from numerical resistivity.We also show that by locally increasing the Ohmic diffusion we can stabilize our scheme and find convergence to the semi-analytic solution in this case.In Section 5, we demonstrate that our new scheme with the additional Ohmic diffusion is able to reproduce the analytical dispersion relation for the whistler and ion-cylcotron waves.We also compare the diffusivity of our scheme with an implementation using the more diffusive HLL Riemann solver instead of an Ohmic diffusion.As a first application, in Section 6 we present simulations of the collapse of a magnetised cloud with constant non-ideal MHD parameters.We show that, depending on the sign of the Hall effect, the magnetic braking can be enhanced or weakened.Finally, we summarise our results in Section 7.

The non-ideal MHD equations
The non-ideal MHD equations can be written as a conservation law: Here, we introduced the state vector , the ideal MHD flux function  ideal , the flux function  res for the non-ideal MHD effects of Ohmic and ambipolar resistivity, and the flux function  H describing the Hall effect.They are given in Heaviside-Lorentz units by: where , , , ,  = /||, , and  = ∇ ×  are the density, velocity, total energy per unit mass, magnetic field strength, direction of the magnetic field, pressure, and electric current, respectively.The specific energy  =  + 1 2  2 + 1 2  2 consists of the thermal energy per mass , the kinetic energy density 1 2  2 , and the magnetic field energy density 1 2  2 .The pressure  =  gas + 1 2  2 includes a thermal and a magnetic component.The system of equations is closed by the equation of state (EOS), which expresses  gas as a function of the other thermodynamical quantities.In this paper, we mostly use an isothermal EOS, with constant isothermal sound speed   .In some cases we also use an adiabatic equation of state where we introduced the adiabatic coefficient  = 5/3.The diffusivities  OR ,  AD and  H describe the strength of the Ohmic and ambipolar diffusion and the Hall effect, respectively, which are in general a function of the magnetic field strength and the local chemical composition.In Zier et al. (2023) (Paper I) we discussed the resistive term  res and in this paper we will concentrate on the Hall term  H .

Linear MHD waves in Hall MHD
As one can see in equation ( 4), the Hall effect does not affect the energy equation but only leads to a change in the topology of the magnetic field.In contrast to the other two non-ideal MHD effects it is not diffusive but dispersive, which has not only an impact on its physical nature but also on the stability of numerical schemes.This can be seen by assuming a system with a purely vertical magnetic field  =  ê , and by adding perturbations   ,   ,   , and   to the velocity and magnetic field, which are travelling in the -direction, i.e. they are proportional to   (  −) .If one only takes into account the induction equation with the Hall effect, and the momentum equation with the Lorentz force, one ends up with the linearized equations (Marchand et al. 2018): This system only has a non-trivial solution if the dispersion relation with the Alfvén velocity   = / √  is fulfilled.It has two solutions for a positive frequency  (Balbus & Terquem 2001): which simplify to the standard Alfvén waves for  H = 0.These new waves are the whistler wave (+ sign) and the ion cyclotron wave (− sign), and their velocities are given by: which shows that the Hall effect is dispersive.One of the waves always travels faster than the Alfvén speed, and for small perturbations (large ) its velocity is even unbound.This unphysical behaviour is a direct consequence of neglecting the electron inertia and assuming effectively an infinite speed of light (charge neutrality).By using a two-fluid model one can show that the frequency of whistler waves is actually bounded by the electron cyclotron frequency (Srinivasan & Shumlak 2011).From a numerical point of view the unbounded velocity is also problematic since this means that perturbations on the grid scale (including unphysical noise) travel the fastest, which can make numerical schemes unstable.In pseudo-spectral codes that solve the system in Fourier space shortrange wavelengths can be easily removed, and they therefore show better stability behaviour.In fact, Kunz & Lesur (2013) demonstrate that their pseudo-spectral algorithm is stable if one uses a third-order time integration scheme, while Falle (2003) showed that first and second order explicit schemes are unconditionally unstable.Methods that integrate in real space and only use a second-order accurate time integration scheme therefore require some kind of dissipation which efficiently damps small scale whistler waves.This can be a physical (Ohmic or ambipolar diffusion), an artificial but explicit (e.g.artificial resistivity in SPH) or an implicit numerical diffusion (e.g. through the Riemann solver).

Ideal MHD in AREPO
In this paper, we use the moving mesh code ARPEO (Springel 2010;Weinberger et al. 2020) to solve equation (1) on a moving, unstructured Voronoi mesh with the finite volume method using the Riemann HLLD solver (Miyoshi & Kusano 2005).Springel (2010) introduced the code and described the general algorithms for mesh construction and hydrodynamic simulations.Pakmor et al. (2011) and Pakmor & Springel (2013) extended it to ideal MHD, and Pakmor et al. (2016) introduced a new gradient estimate and time integration to achieve in most test problems second-order convergence in time and space.
Zier & Springel (2022a) further improved the method by introducing a higher-order flux integration that reduces the production of noise on the grid scale, which is especially useful in shear flows.
The continuum equations (1) conserve the condition ∇ •  = 0, but this property can get lost once the equations are discretized.This can in turn lead to purely numerical instabilities and therefore to large deviations from the ∇ •  = 0 condition.Pakmor & Springel (2013) introduced the Powell scheme (Powell et al. 1999) in ARPEO for convergence control that diffuses the magnetic monopoles whereas Pakmor et al. (2011) implemented the Dedner scheme (Dedner et al. 2002) that additionally damps it.In this paper, we will use the former approach in all simulations.
We note that static grid codes often use the constrained transport method that is able to conserve the condition ∇ •  = 0 up to machine precision even for the discretized equations (Evans & Hawley 1988).Implementing this method for dynamic meshes becomes very complicated, however, and also tends to be more diffusive (Mocz et al. 2014).

Magnetic diffusion in ARPEO
In Paper I we introduced a new diffusion solver to integrate the resistive flux function  res in equation (1) using operator splitting and the finite volume method.We use the same basic approach also for the Hall effect.By integrating equation (1) over the volume   of a Voronoi cell  we find for the non-ideal part: Here we introduced the vector of conserved quantities   = ∫    and applied Gauss' theorem.For the time integral we use the standard second-order accurate Runge-Kutta scheme while the surface integral can be rewritten as a sum over all interfaces of one cell.The integral over one interface can be approximated by evaluating the Flux function at the center of the interface and by multiplying it with its area.We note that in Zier & Springel (2022a) we showed that this approximation can become inaccurate for the ideal MHD flux in shear flows, where a higher order approximation is required.For the non-ideal MHD terms we did not find a significant improvement by using those more expensive integration rules and therefore we will not use them here.For the integration of the non-ideal terms we use again the Strang-splitting scheme.
To apply the Hall term we are left with evaluating  H at the interface, which requires according to equation (4) accurate estimates for the magnetic field and its gradients.In Paper I we introduced a method based on performing a least square fit at the centre of the interface, taking into account the values of the magnetic field at all neighbouring cells.We will use the same method in this paper and therefore refer to Paper I for full details.

Time step constraints
The Hall effect leads to an additional time step constraint which can be written as (Bai 2014): where Δ  = [3  /(4)] 1/3 is the effective radius of the Voronoi cell  with volume   in three dimensions ( √︁   / in two dimensions) and  is the number of dimensions. H is the Courant number for the Hall effect which we typically set to 0.5.Like the other non-ideal MHD terms, the Hall effect leads for high resolutions to a very restrictive time step constraint.In contrast to the other two effects the Hall effect cannot be accelerated by using super-time stepping (see Berlok et al. (2020) for an implementation in AREPO for Braginskii viscosity), so the only option would be subcycling.In this paper we will not use such methods but rather observe the time step constraint, which is not yet a problem at the only moderately high resolution we employ.

The strength of the Hall effect
The strength of the Hall effect can be parameterised by the diffusion coefficient  H .It strongly depends on the local chemical composition and magnetic field strength.As a consequence factors such as the cosmic ray ionization rate and the dust size distribution can have a large impact on  H .In contrast to the coefficients for Ohmic and ambipolar diffusion, the Hall coefficient can be both positive and negative, and will generally have different signs throughout the simulation domain.Although there exist library solutions (e.g.Wurster 2016) that allow the calculation of the non-ideal MHD diffusion coefficients in equilibrium the values of  H can vary over several orders of magnitudes.It is therefore difficult to anticipate when and where the Hall effect is dominating (Wurster & Li 2018).
The goal of this paper is to demonstrate the accuracy of our nonideal MHD solver, and therefore we will use only constant or parameterized coefficients: where  is a constant, for simplicity and definiteness.The Hall effect dominates over the ideal MHD induction term if  H /2 ≫   holds, which can be translated to a critical length scale  H =  H /  .This length scale allows us to define the Lundquist number for perturbations with wave vector .

STABILIZING THE HALL TERM IN THE HALL-DOMINATED REGIME
As we have discussed in Section 2.2, the Hall-MHD equations without resistivity allow for whistler waves whose velocity becomes arbitrarily large for small wavelengths.This can destabilize our scheme since numerical noise that is inherent in simulations on a moving mesh (see e.g.Zier & Springel 2022a) travels the fastest.In fact, we found that our scheme becomes unstable in simulations of a gravitational collapse of a protostellar core if we only take into account the Hall effect and no diffusion.The instability can be observed first in the densest region, which is equivalent to the region with the highest spatial resolution.This behaviour can be understood because the numerical diffusion, which always exists in our scheme, decreases with better spatial resolution, which in turn means that its stabilizing influence on the Hall effect diminishes in the core.We thus expect that it first falls there below the stability threshold for the Hall effect.Kunz & Lesur (2013) show that for a third-order time integration scheme no additional diffusion would be required, but this is complicated to implement in a multi-physics code such as AREPO.In the literature different methods were therefore introduced to stabilize the Hall effect if there is not enough physical resistivity but a secondorder accurate time integration scheme is still kept.O' Sullivan & Downes (2006) proposed a semi-implicit integration method, which is marginally stable.In this scheme, one Cartesian component of the magnetic field is updated with a fully explicit timestep.The second component uses this updated field for a semi-implicit step.Finally, the last step for the third component is fully implicit, using the two already updated components.The generalization to 3D was performed in O' Sullivan & Downes (2007) and implemented into ATHENA by Bai (2014).Bai & Stone (2017) found that in spherical coordinates this scheme becomes unstable and therefore used instead a modified HLL Riemann solver that takes into account the whistler waves.The HLL Riemann solver is significantly more diffusive but its increased numerical resistivity helps to stabilize the Hall effect.Marchand et al. (2019) showed that in simulations of the collapse of protostellar cores, the modified HLL solver can lead to the artificial creation of angular momentum if the whistler wave is accounted for in all variables.We will show in Section 5.2 and Appendix B that the HLL Riemann solver is indeed more diffusive, and that this can reduce the numerical accuracy even in regions where the Hall effect is subdominant.
Another option to stabilize the Hall effect is to use an artificial resistivity, which is anyway required in SPH simulations.One can also introduce a hyper-resistivity, which is very efficient in damping noise on small scales but does not significantly affect larger scales (Srinivasan & Shumlak 2011).This would require the calculation of higher derivatives of the magnetic field, which is however very expensive and not particularly accurate on an unstructured mesh.As we will discuss in the following we employ instead a standard artificial resistivity.

The stability of a C-shock
In ideal MHD, shocks contain discontinuous jumps in the parallel component of the magnetic field and in the hydrodynamic quantities.The non-ideal MHD terms smooth these discontinuities, making them continuous ("C-type" shock, Draine 1980).As we have shown in Paper I, Ohmic and ambipolar diffusion lead to a pure smoothing of the jump, and all quantities stay monotonous functions in space in the direction normal to the shock front.This changes for the Hall effect, which introduces oscillations close to the shock front.These oscillations have to be damped by Ohmic or ambipolar diffusion, otherwise they would continue in the steady state throughout the whole system.This can be seen by solving for the static profile of a C-shock, which we further discuss in Appendix A1.The smaller the explicit diffusion is compared to the Hall effect, the more oscillations can be observed.In simulations lacking any explicit diffusion the oscillations will nevertheless be damped by numerical diffusion  N , making the C-shock an ideal system to measure  N but also for investigating the numerical stability of the Hall term.
We use a similar setup as in Paper I, which means we will simulate a box of size   ×   ×   with   =   .The shock is travelling in the -direction and we use periodic boundary conditions in the -and -directions.In the -direction we extend the box by 0.5 on both sides and set for −0.5 <  < 0 all primitive variables to the pre-shock values, and for   <  <   +0.5 to the post-shock values.We use for the extended box periodic boundary conditions, which ensures that the total number of cells stays constant.With a fixed box size it becomes expensive to perform resolution studies, since to double the spatial resolution in the -direction we would have to use 8 times the amount of cells.To circumvent this problem we use a box size   = 10/, where  3 is the desired number of cells within a unit box.For strong Ohmic diffusion we can use   = 1 but for simulations with weaker diffusion the shock becomes more extended and we have to increase   .We present the box sizes for simulations with explicit diffusion in Table 1.In Table 2 we give the pre-and post-shock values for our simulations.We use  H = 0.02 || =  ||, and as we show in Appendix A1 the structure is fully determined by the ratio  OR /.The setup is inspired by Marchand et al. (2018), and we use an isothermal equation of state with sound speed   = 0.1.
We run all simulations until  = 5 and afterwards bin the Voronoi cells in the -direction.We also calculate for each bin the standard deviation as a measure of the noise in our simulations, which is only significant in the density profile.We use 2 bins per unit length which means that on average 50 cells can be found in each bin.Afterwards, we calculate profiles for different  OR and try to fit the profile that minimizes the error for   .We introduced here the average magnetic field B,, of our simulations in the -direction in each bin, and the value at the center of each bin for the semi-analytic profile  ,, .
For the initial conditions we use two different setups: In the first one we set the primitive variables to the pre-shock value for  <  0 and to the post-shock value for  0 < .In the second one we smooth the initial conditions by using the equation for the primitive variable  , where the subscript  refers to the preshock value and  the post-shock value.We typically use  0 =   −0.5 as the initial position of the shock.As in Paper I we construct initially an unstructed mesh by replicating a box with 10 3 cells.We perform simulations with explicit diffusion for which we expect convergence to the semi-analytical profile and simulations without explicit diffusion, which allows us to measure   and the maximum ratio  H / OR for which our scheme is stable.

With explicit diffusion
As one can see in Fig. 1 our method is able to accurately resolve the structure of a C-shock even in the strongly Hall-dominated regime.The largest deviations are found in the density profile, which is also the only quantity showing significant noise.As we have already discussed in Paper I this noise can also be observed without non-ideal MHD and therefore is not caused by inaccuracies in our non-ideal MHD solver.
We performed a parameter study in which we varied  OR and the resolution, for which we show the most important parameters in Table 3.The numerical resistivity can be defined as the difference between the measured one and the explicit one we set at the beginning of the simulation.As expected it strongly decreases with increasing resolution, and we find convergence to the physical resistivity.This can also be seen in Fig. 2, where we find a similar scaling of   with the resolution .But we also observe that the trend cannot be described well by a  −2 power law.For  OR / ≥ 0.1, all simulations stay stable whereas for smaller explicit diffusivity some simulations with high resolution become unstable.We observe that simulations with discontinuous initial conditions become immediately unstable at the beginning of the calculation while simulations with continuous initial conditions first form a C-shock profile before they become unstable. OR .We use on average  cells per unit length in the -direction and measure the total resistivity  tot by fitting the semi-analytic profile to the simulation data.We performed simulations with discontinuous (D) and continuous (C) initial conditions, and also give the final  1 error for   , which we try to minimize.

Without explicit diffusion
To ensure that the full C-shock can be resolved in simulations without explicit diffusion, we use in this section a fixed box size   = 10.As one can see in Fig. 3, the numerical resistivity acts indeed similarly to a physical one which allows us to measure   by close inspection of the shock profile.In Table 4 we give an overview of all simulations we performed without physical resistivity.  is more or less independent of the initial conditions and decreases with increasing resolution.For  = 90 the simulations are still stable but we can observe larger noise especially in the density, which is also reflected in the strongly increased  1 error compared to  = 80.For even higher resolution the simulations become unstable, so we conclude that  = 90 is close to the instability.In Fig. 4 we show the dependence of   on the resolution.We find that   does not follow a simple power law, which prevents us from defining a simple lower boundary for the numerical resistivity.

Convergence
To verify the convergence of our method we define the  1 error norm for each quantity at time  = 5 following equation ( 18), and we compare our profiles directly with the semi-analytic result for the same physical resistivity.As one can see in Fig. 5 we find for strong Ohmic resistivity at least second-order convergence.For  OR ≤ 0.001, the errors do not decrease for low resolution which can be explained by the fact that the numerical resistivity is dominating and therefore the profile in the simulation corresponds to one with larger  OR .By further increasing the resolution we can find for  OR = 0.001 decreasing errors, although we cannot further increase the resolution since otherwise the simulations become unstable.4. Overview of simulations of a C-shock without explicit Ohmic resistivity.We use on average  cells per unit length in the -direction and measure the numerical resistivity  N by fitting the semi-analytic profile to the simulation data.The simulation can be stable or unstable for discontinous (D) or continous (C) initial conditions.We use a box size   = 10 to make sure that the boundary does not affect the results and give the average  1 error for   at  = 5.

Artificial resistivity
Based on a better understanding of the stability of our implementation for the Hall effect we opt for an approach where we deliberately increase the magnetic diffusivity by the small amount required for stabilization.This prevents us from using the HLL-modified Riemann solver in the full domain.Instead we locally increase the Ohmic resistivity.The condition for stability of the Hall effect can be formulated as where we introduced a coefficient  1 of order unity.The contribution of the spatial discretization to the numerical resistivity can in principle be parameterized by (Rembiasz et al. 2017): for a second order accurate scheme, where  2 is another parameter of order unity, Δ  = [3  /(4)] 1/3 is the approximate size of a cell ,  denotes the characteristic length scale of the problem, and   gives the local wave speed, which we set equal to the velocity of the fast magneto-sonic wave,   = √︃  2  +  2 /.But as we have seen in the last section, the factor  2 can also depend on Δ in our scheme, and in general it is not always obvious what the characteristic length scale of a problem is (Rembiasz et al. 2017).
This makes its hard for us to find the required lower boundary for  2 .We will therefore neglect the numerical diffusivity in equation (20) and instead use the condition: Following the results presented in the last section we use which is also found to be stable for all simulations of the collapse of a protostellar core we present later in the paper.We note that by neglecting the numerical diffusivity we also apply some artificial resistivity in regions where   is already sufficient to stabilize the Hall term, but in this case we also expect  OR <   , and therefore this should not influence our results.

PROPAGATION OF A WHISTLER WAVE
As we have seen in Section 2.2, the Hall effect introduces two additional MHD waves.We perform here a similar study as in Marchand The frequency of the oscillation is modified to: Both frequencies are linear combinations of the two frequencies defined in ( 12), and therefore a Fourier decomposition of the temporal signal allows use to measure the dispersion relation for the whistler and the ion cyclotron waves.We simulate a cubic box of size   ×   ×   = 1 × 10/ × 10 with uniform density and pressure. 3 is here again the number of cells we expect in a cubic box with size unity.The initial conditions are given by: and we use periodic boundary conditions, which restricts  to be an integer times 2.We use an adiabatic equation of state with  = 5/3, and the full analytical solution including  is given in Appendix A2.

Linear dispersion relation
As a first test we try to measure the dispersion relation ( 12) using a resolution of 160×10×10 cells initially ( = 160).We introduce the Hall frequency   =  2  / H , which allows us to rewrite equation ( 12) as: The normalized frequency is now only a function of the product   , with   =  H /  .This product is the inverse of the Lundquist number we introduced in equation ( 17).In Table 5 we give an overview of the parameters we used in the simulations to measure the dispersion relation and Fig. 6 shows our results.We find perfect agreement with the theoretical dispersion relation.For large   we are only able to measure the high-frequency branch since the lower frequencies require a much longer measurement.

Diffusivity of our scheme
Both explicit and numerical diffusion lead to a damping of the amplitude of the wave.The damping rate is in general a function of the resolution and the wave vector , which we will further analyze in this section.For comparison, we also implemented the Hall modified HLL solver from Lesur et al. (2014), which we also shortly discuss in Appendix B1.As an example, we show in Fig. 7 the temporal evolution of the magnetic and kinetic energy for a setup with  = 160,  = 10 (2) and  H = 0.005.At the beginning, we can observe for the HLLD Riemann solver a strong decay of the magnetic field, and at a later time a decay of the magnetic field with the same decay rate as the kinetic energy.For the HLL Riemann solver both types of energy decay at a similar rate but in general significantly faster than for the HLLD Riemann solver.The faster decay of the kinetic energy in this case can also be explained by the more diffusive behaviour of the HLL Riemann solver since it also increases the numerical viscosity.
The evolution of the energy can be approximated by an exponential decay with a damping rate , which we define as half of the decay rate of the magnetic energy.By fitting a function proportional to  −2 to the evolution of the magnetic field energy, we can measure this decay rate as a function of  and resolution.We perform simulations with  = 160 and  = 320, different  and  H = 0.005, and show the results in Fig. 8.As expected for the HLL Riemann solver the damping rate is proportional to  4 (Marchand et al. 2018) and decreases by a factor of 4 if we double the resolution.For the HLLD Riemann solver, the damping rate is a significantly weaker function of .For small  it is approximately independent of the resolution since in this case the damping is dominated by the Ohmic diffusion.For larger  the numerical diffusion plays a role and the damping rate can be reduced by increasing the resolution.

COLLAPSE OF A MAGNETIZED CLOUD
To test our new module for non-ideal MHD in an environment closer to typical science applications, we repeat in this section the simulation of the collapse of a rotating, magnetized cloud under its own gravity performed in Paper I but now with the Hall effect.This setup acts as a simplified model for star formation, allowing us to test the stability of our methods in a situation with a large dynamic range of scales.
The initial conditions are again taken from Pakmor et al. (2011) and Marinacci et al. (2018a), which were in turn adapted from Hennebelle & Teyssier (2008).They consist of a rigidly rotating, homogeneous sphere with radius  0 = 0.015 pc and mass 1 M ⊙ embedded via a small transition region in a low density environment with a density contrast of 100.The total size of the simulation box is 0.06 pc.The initial density of the cloud is 4.8 × 10 −18 g cm −3 , which translates to a free fall time of 3 × 10 4 years.We choose a ratio of 0.045 between rotational and gravitational energy, which is equivalent to a rotational period of 4.7×10 5 years for the rigid body rotation around the -axis.The cloud is initially penetrated by a homogeneous, purely vertical magnetic field of magnitude 107 G, equivalent to a mass-to-flux over critical mass-to-flux ratio of 5.6.We use a barotropic equation of state (Hennebelle & Teyssier 2008) with sound speed  ,0 = 0.2 km s −1 and critical density   = 10 −13 g cm −3 , which represents a transition from an initially isother-  .Formation of a dense core in the collapse of a magnetized cloud under its own self-gravity with different models for the magnetic field.The left column shows a simulation with ideal MHD, while the middle one adds a constant positive, and the right one a constant negative Hall coefficient as well as Ohmic diffusion for stability.All the panels show projected slices in the -plane with the slice having a thickness of 0.1 times the side length of the plot.In the top row, we display a volume-weighted density projection of the gas in the central region (0.03  0 in diameter).The middle and bottom row show density-weighted projections of the -and the azimuthal component of the magnetic field, respectively, in a region 10 times as large.Note that the -component of the magnetic field is negative (inverted with respect to the initial condition) due to the outflow in some regions.Results from the three simulations are shown after an identical amount of elapsed time, approximately 3.4 × 10 4 years.mal EOS to an adiabatic one at higher densities.We start with an initial resolution of 128 3 cells on a uniform Cartesian mesh and apply periodic boundary conditions at the edges of the simulation box for all quantities except gravity.Cells are refined if their free fall time is smaller than ten times their sound crossing time.We also introduce a minimum volume of 5 × 10 −17 pc which is equivalent to an effective resolution of 16384 3 cells.The latter condition reduces the computational cost, especially considering the quadratic dependence of the diffusive timestep on cell radius.Here, we set a constant coefficient of  H = ±10 18 cm 2 s −1 to demonstrate the aforementioned bimodality of the Hall effect in either promoting or suppressing disk formation.To stabilize our simulation we also use  OR = 2 • 10 17 cm 2 s −1 .
The drift of magnetic field lines in comparison to ideal MHD due to the Hall effect is In our setup, the initial rotation and magnetic field axes are aligned.
In this case, the azimuthal current resulting from the radial bending of the magnetic field is parallel to the rotational velocity of the gas.Depending on the sign of  H , the magnetic braking caused by the azimuthal bending of the magnetic field is either weakened ( H > 0) or further enhanced ( H < 0) by the presence of the Hall effect.The effect of this bimodality is shown in Figure 9, where positive and negative Hall diffusivity show diverging trends with respect to ideal MHD.In the case with positive coefficient, the drift of magnetic field lines caused by the Hall effect leads to a weaker azimuthal field and thereby to a broader and slower outflow compared to ideal MHD.On the other hand, the azimuthal field is enhanced and the outflow is both faster and more strongly collimated if the coefficient is negative.
The effect on the azimuthal field can in Fig. 9 most clearly be seen in the fact that in the bottom panels a much larger section around the midplane has a very small value (black in the image) in the simulation with the positive coefficient than with the nagative one.
Note that when the Hall coefficient is (locally) given by a chemical model, the bimodality is instead usually characterized by an (anti-)alignment of the axes of rotation and magnetic field.For in-stance, Wurster et al. (2021) show that since in their model the Hallcoefficient is negative throughout the early evolution of envelope and disk, an aligned configuration of magnetic field and rotation axes enhances the azimuthal field and removes angular momentum from the disk -this corresponds to the situation in the right panel of Figure 9.They even find a counter-rotating inner disk at the end of the simulation.The middle panels can on the other hand be seen as a simplified model of the anti-aligned configuration (in our setup the signs of the Hall coefficient and the current are both reversed with respect to theirs, cancelling out).Similarly, Zhao et al. (2020) observe either enhanced or suppressed disk formation, depending on alignment.Bai & Stone (2017), in contrast to the two previously mentioned works, assume a positive Hall coefficient throughout their disks and consequently find instead enhanced outward angular momentum transport in the anti-aligned configuration.

SUMMARY AND CONCLUSIONS
In this paper we have extended the non-ideal MHD module in the moving mesh code AREPO to include the Hall effect.The Hall effect is not diffusive but dispersive, leading to the introduction of two new MHD waves, the whistler and the ion-cylcotron waves.It is known to destabilise second order accurate time integration schemes without sufficient magnetic diffusion, leading to various approaches in the past to reach stability by increasing the diffusion.
Our implementation of the Hall effect follows the method we introduced in Zier et al. (2023) for modelling Ohmic and ambipolar diffusion, which relies on a least-squares fit to estimate the magnetic field gradient at the interfaces between neighbouring cells.By simulating a magnetic C-shock, we show in Section 4 that our scheme suffers from the typical numerical instability at high resolution.By adding sufficient Ohmic resistivity we can stabilize our method and obtain second-order accuracy.The critical ratio below which we find the instability for our scheme is  OR / H < 0.1.To be on the safe side, we therefore conservatively propose to enforce the condition  OR ≥  H /5 in all cells.For this choice and for these parameters, we did not find any signs of instability.
Using this new condition, we have shown in Section 5 that we are able to accurately resolve the dispersion relation of the new MHD waves that appear for the Hall effect.We compared the damping rate for the magnetic field with a similar implementation using a modified HLL Riemann solver instead of an artificially increased Ohmic diffusion.We showed that the increased numerical diffusion of the HLL Riemann solver is particularly efficient on small scales, but also increases the numerical viscosity.As a first example of a more complicated setup, we simulated in Section 6 the collapse of a magnetised gas cloud with a constant Hall diffusion coefficient.We demonstrated that, depending on the coefficient of the Hall effect, the magnetic braking can be weakened or strengthened.
We developed our module with the collapse of protostellar cores in mind, since in this case we can directly profit from the high accuracy and adaptive resolution of the moving mesh approach.In these simulations, the resolution within a potentially forming disk is rather low, and we need a low numerical viscosity to accurately follow their evolution.The Hall effect may only dominate within small regions of the system (Wurster & Li 2018), and therefore we feel that using the more diffusive HLL Riemann solver globally is excessive, rather we consider a localized increase of the Ohmic diffusion preferable.We are currently also working on more realistic simulations with nonconstant diffusivities (Mayer et al., in preparation).For the future we also plan to implement a third-order accurate time integration scheme for the Hall effect similar to the one we already implemented in Zier & Springel (2022a) for ideal MHD.This scheme might be stable without requiring any artificially increased diffusion.To overcome the inherent stability problems of the operator-splitting approach, we have alternatively implemented a modified HLL Riemann solver to incorporate the Hall effect directly into the ideal MHD flux calculation.Tóth et al. (2008) introduced this method which was later implemented in Lesur et al. (2014) in the PLUTO code.We follow their implementation and solve the equation: which is problematic because the flux depends on the magnetic field and its derivative via the current  = ∇ × .This makes the Riemann problem ill-defined, since  is not defined if there is a discontinuity in the magnetic field.To get around this problem, we treat  as an external parameter to the Riemann problem and estimate it as the average of the currents of the two cells of the interface, using the standard gradient estimates in AREPO.
As other studies in the past (Lesur et al. 2014;Bai & Stone 2017;Marchand et al. 2018Marchand et al. , 2019) ) we have not been able to modify the accurate HLLD or Roe Riemann solvers to deal with this ill-defined Riemann problem and therefore we use the diffusive HLL Riemann solver as a basis instead.The HLL Riemann solver requires only the smallest algebraic signal velocity of the left state and the largest algebraic signal velocity of the right state.We define them as: where  is the fluid velocity,  f is the fast magnetosonic speed and  w is the whistler wave speed.Since the Hall effect is dispersive, we use the cell radius  to calculate the fastest whistler wave.
For the HLL Riemann solver the flux to update the conservative variables is given by: where  L,R are the left and right fluxes.In contrast to our standard scheme with the HLLD Riemann solver and Ohmic diffusion, this method reduces the accuracy even in the absence of the Hall effect.

B2 Comparison of the HLL and HLLD Riemann solver
To better understand the influence of the Riemann solver on our results we repeated the simulations presented in Section 6 with the Hall modified HLL Riemann solver instead of an artificial Ohmic diffusion.As we can see in Fig. A1 the HLL Riemann solver not only is more diffusive for the magnetic field but also directly influences hydrodynamic quantities such as the density.We also repeated the simulations with the Hall effect and show in Fig. A2 the azimuthal magnetic field.In comparison to the simulations with HLLD Riemann solver, the influence of the Hall effect is also reduced, showing only very small differences in the outflow structure.While the general trends expected from the inclusion of the Hall effect are still captured by the modified HLL solver, it appears that the increased diffusivity relative to the HLLD solver that we showed in Fig. 7 for the wave test also has a significant effect in the context of more complicated applications.

B3 Effect of the artificially increased Ohmic diffusion
As we have discussed in Section 4, simulations with the Hall effect become unstable if there is too little diffusion (numerical or physical) to stabilize it.But by adding an unphysical Ohmic diffusion we also might influence the outcome of our simulations.We therefore performed simulations with  H = 10 18 cm 2 s −1 for which the numerical resistivity of the HLLD Riemann solver alone is high enough to stabilize the code for our resolution.Nevertheless, as one can see in Fig. A3 there is some noise observable in the magnetic field.We repeated the simulation with an additional Ohmic diffusion  OR = 2 × 10 17 cm 2 s −1 and show in Fig. A3 that this additional diffusion removes the noise and leads to a smoother magnetic field.But  the large scale structure is still very similar.The additional Ohmic diffusion reduces the efficiency of magnetic braking which leads to a stronger disk like structure observable in the density field.We note that in simulations with realistic values for the diffusivity coefficients, the Hall effect is dominated almost everywhere by either ambipolar diffusion (low density) or Ohmic diffusion (high density).Therefore, instabilities resulting from the Hall effect are suppressed in most situations by the presence of strong physical diffusion.Otherwise, they can be mitigated by locally increasing the diffusion above its physical value, ensuring that the threshold for stability ( OR,AD > 0.2 | H | based on our C-shock tests) is guaranteed everywhere.

Figure 1 .
Figure 1.The structure of the C-shock at  = 5 for a simulation with physical resistivity  OR = 0.001 and  = 160.The black line corresponds to the best fit semi-analytically expected structure of the shock with  OR = 0.0012, while the red line shows the average values in a simulation with AREPO at time  = 5.The blue lines show the ±1 standard deviation in each bin.

Figure 2 .
Figure 2. The numerical resistivity  N as a function of resolution for different simulations of the C-shock with physical resistivity  OR .We show the averaged values for continuous and discontinuous initial conditions.

Figure 3 .
Figure 3.The structure of the C-shock at  = 5 for a simulation without physical resistivity,  = 80 and continuous initial conditions.The black line corresponds to the best fit semi-analytically expected structure of the shock with  OR = 0.0009, while the red line shows the average values in a simulation with AREPO at time  = 5.The blue lines show the ±1 standard deviation in each bin

Figure 4 .Figure 5 .
Figure 4.The numerical resistivity   as a function of resolution for simulations of the C-shock without physical resistivity.We show the results for continuous and discointous inital conditions.

Figure 6 .
Figure 6.The dispersion relation of the whistler and ion cyclotron waves (black line, see equation 26) and the measured relation in our simulations (red symbols).More information can be found in Section 5.1.

Figure 7 .Figure 8 .
Figure7.The temporal evolution of the magnetic and kinetic energy for simulations of an MHD wave with  = 10 (2  ) and resolution  = 160.We perform the simulation once with the HLLD Riemann solver together with Ohmic diffusion and once with the HLL Riemann solver.

Figure 9
Figure9.Formation of a dense core in the collapse of a magnetized cloud under its own self-gravity with different models for the magnetic field.The left column shows a simulation with ideal MHD, while the middle one adds a constant positive, and the right one a constant negative Hall coefficient as well as Ohmic diffusion for stability.All the panels show projected slices in the -plane with the slice having a thickness of 0.1 times the side length of the plot.In the top row, we display a volume-weighted density projection of the gas in the central region (0.03  0 in diameter).The middle and bottom row show density-weighted projections of the -and the azimuthal component of the magnetic field, respectively, in a region 10 times as large.Note that the -component of the magnetic field is negative (inverted with respect to the initial condition) due to the outflow in some regions.Results from the three simulations are shown after an identical amount of elapsed time, approximately 3.4 × 10 4 years.

Figure A1 .Figure A2 .
Figure A1.A comparison of the density (top) and azmithuthal magnetic field (bottom) structure in ideal MHD simulations of the collapse of a magnetized cloud core with the HLLD Riemann solver (left side) and the HLL Riemann solver (right side).We use the same parameters as we use in Fig. 9 to produce the plot.

Figure A3 .
Figure A3.Density and magnetic field slices in simulations of the collapse of a magnetized cloud using different resistivity coefficients.The left panels have solely  H = 10 18 cm 2 s −1 , while the right panels show simulations that additionally include  OR = 2 × 10 17 cm 2 s −1 .

Table 1 .
The box size for simulations of a C-shock with explicit diffusion

Table 2 .
Initial conditions for the C-shock tests with sound speed   = 0.1.

Table 3 .
Overview of simulations of a C-shock with explicit Ohmic resistivity

Table 5 .
Overview of simulations we performed to measure the dispersion relation of the Hall MHD waves.We give the wave vector , the Hall diffusivitiy  H , the product   , the total simulation time  max and the minimum resolved time  to measure the frequencies.etal.(2018), i.e. we try to measure the numerical dispersion relation and the damping rates for those waves.We will simulate a polarized, travelling wave in the -direction with a guide field   and wave number .By assuming  , ≪   the Hall coeffcient  H becomes constant and the Hall effect leads to a modulation of the amplitude of the magnetic field with frequency: