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

HAL
open science

# Impact of the Hall effect in star formation and the issue of angular momentum conservation

P. Marchand, B. Commerçon, G. Chabrier

► To cite this version:

P. Marchand, B. Commerçon, G. Chabrier. Impact of the Hall effect in star formation and the issue of angular momentum conservation. Astronomy & Astrophysics - A&A, 2018, 619, pp.A37. 10.1051/0004-6361/201832907. [hal-02371328](https://hal.science/hal-02371328v1)

HAL Id: hal-02371328

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

Submitted on 19 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Impact of the Hall effect in star formation and the issue of angular momentum conservation

P. Marchand<fnref n="1" />,<fnref n="2" />, B. Commerçon<fnref n="2" />, and G. Chabrier<fnref n="2" />,<fnref n="3" />

Received 26 February 2018 / Accepted 29 July 2018

## ABSTRACT

We present an implementation of the Hall term in the non-ideal magnetohydrodynamics equations into the adaptive-mesh-refinement code RAMSES to study its impact on star formation. Recent works show that the Hall effect heavily influences the regulation of the angular momentum in collapsing dense cores, strengthening or weakening the magnetic braking. Our method consists of a modification of the two-dimensional constrained transport scheme. Our scheme shows convergence of second order in space and the frequency of the propagation of whistler waves is accurate. We confirm previous results, namely that during the collapse, the Hall effect generates a rotation of the fluid with a direction in the mid-plane that depends on the sign of the Hall resistivity, while counter-rotating envelopes develop on each side of the mid-plane. However, we find that the predictability of our numerical results is severely limited. The angular momentum is not conserved in any of our dense core-collapse simulations with the Hall effect: a large amount of angular momentum is generated within the first Larson core, a few hundred years after its formation, without compensation by the surrounding gas. This issue is not mentioned in previous studies and may be correlated to the formation of the accretion shock on the Larson core. We expect that this numerical effect could be a serious issue in star formation simulations.

**Key words.** stars: formation – magnetohydrodynamics (MHD) – ISM: magnetic fields – methods: numerical

## 1. Introduction

Star formation originates in the collapse and fragmentation of molecular clouds. These giant structures of the interstellar medium harbor thousands of solar masses of gas, mostly hydrogen (~74% by mass) and helium (~25%). This fragmentation leads to the creation of over-dense regions, the dense cores, that dynamically decouple from the rest of the cloud. When a dense core gravitationally collapses onto itself, the accumulation of matter and the rise of the temperature form a star embryo, a protostar.

Larson (1969) describes the early stages of the collapse of dense cores. Initially, the collapse is isothermal as the gas-dust mixture is transparent to its own blackbody radiation. When the density reaches $10^{-13}$ g cm$^{-3}$, this radiation is trapped and heats up the gas, leading to the formation of an hydrostatic core called the first Larson core. The core contracts adiabatically until the temperature reaches 2000 K. At this moment, $H_2$ molecules dissociate in an endothermic reaction, thus creating a gap of energy that allows a second collapse. Once all hydrogen is atomic, the adiabatic contraction resumes in what is called the second Larson core. While the second core is considered to be a protostar, the first Larson core remains a theoretical object not confirmed by observations because of its short lifetime (few thousands years at maximum). Nevertheless, serious candidates have been revealed by recent observations (e.g., Pezzuto et al. 2012; Gerin et al. 2015).

The dynamics of the collapse is affected by magnetic fields threading the dense cores. It has been well established that they play a major role during the contraction of dense cores,

and they are considered as the best way to solve the angular momentum problem. Observations show that the specific angular momentum of dense cores (~$10^{21}$ cm$^2$ s$^{-1}$) is three orders of magnitude higher than those of young stars (~$10^{18}$ cm$^2$ s$^{-1}$; Bodenheimer 1995). The magnetic field is able to transport the angular momentum from the inner parts of the dense cores to their envelop through the magnetic braking. In the last 20 yr, there have been several studies investigating the impact of the magnetic braking during protostellar collapse through numerical simulations using the approximation of ideal magnetohydrodynamics (MHD), in which the charged particles are perfectly coupled with the neutral gas. However, the braking is so efficient in removing the angular momentum that, at best, only small rotationally supported disks could be formed (Allen et al. 2003; Matsumoto & Tomisaka 2004; Galli et al. 2006; Price & Bate 2007; Hennebelle & Teyssier 2008; Hennebelle & Fromang 2008; Commerçon et al. 2010; Masson et al. 2016). While the conditions of their formation is still matter to debate, such disks evolve in protoplanetary disks, environments of planet formation that have been widely observed (Murillo et al. 2013; ALMA Partnership et al. 2015; Gerin et al. 2017; Fuente et al. 2017). Additionally, the flux-freezing approximation of ideal MHD causes a magnetic pile-up in the protostar, leading to values of the magnetic flux up to three orders of magnitude higher than those observed (Nakano 1984; Troland & Crutcher 2008; Johns-Krull 2008).

Among the possible solutions to solve both the angular momentum and the magnetic flux problem, non-ideal MHD is considered to be the most efficient and robust. Non-ideal MHD describes the resistance encountered by the charged particles

<footnote n="1">¹ Department of Earth and Space Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
e-mail: marchand@astro-osaka.jp</footnote>

<footnote n="2">² CRAL, Ecole normale supérieure de Lyon, UMR CNRS 5574, Université de Lyon, Lyon, France</footnote>

<footnote n="3">³ School of Physics, University of Exeter, Exeter, EX4 4QL, UK</footnote>
into a not fully ionized fluid. The first studies accounted for the Ohmic and ambipolar diffusion during the early stages of the collapse (i.e., Machida et al. 2006; Duffin & Pudritz 2008; Mellon & Li 2009; Kunz & Mouschovias 2010; Tomida et al. 2015; Masson et al. 2016). They showed that these processes could indeed temper the magnetic braking and regulate the magnetic flux. The third non-ideal MHD process, the Hall effect, has been included only recently in numerical simulations. Its resistivity value indicates that its influence on the system can be similar in strength to the ambipolar and Ohmic diffusion at densities ranging from the isothermal stage to the end of the first Larson core (Marchand et al. 2016; Wurster 2016). It is therefore necessary to quantify its impact through numerical simulations. Recent studies have shown that the Hall effect does indeed play a major role in the regulation of the angular momentum (Krasnopolsky et al. 2011; Tsukamoto et al. 2015, 2017; Wurster et al. 2016) and in the formation of binary stars (Wurster et al. 2017).

Several astrophysical codes, both Eulerian and Lagrangean, already include the Hall effect, through various numerical methods. The first star-formation simulations, including the Hall term (Tsukamoto et al. 2015, 2017) as well as the most recent studies (Wurster et al. 2016, 2017, 2018) have been performed with SPH codes. The protoplanetary disks shearing box simulations by Lesur et al. (2014) uses the grid-based code PLUTO (Mignone et al. 2007). The Hall effect is also implemented in the AMR codes BASTRUS (Tóth et al. 2008) and ZEUS (Stone & Norman 1992; Krasnopolsky et al. 2010), as well as in the multimethods code GIZMO (Hopkins 2015, 2017).

The aim of this paper is first to present an original implementation of the Hall effect in the adaptive-mesh-refinement (AMR) code RAMSES and to apply it in protostellar collapse studies. The paper is organized as follows. First, we briefly recall the basics of the theory about the Hall effect in Sect. 2. Second, we detail our numerical implementation in Sect. 3 and present validation tests of our scheme in Sect. 4. Section 5 is devoted to the protostellar collapse models in which we incorporate the Hall effect. Finally, in Sect. 6 we give our conclusions.

## 2. Theory

### 2.1. Equations of non-ideal MHD

We write the complete set of non-ideal MHD equations in conservative form with self-gravity as a source term:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot [\rho \mathbf{u}] = 0, \qquad (1) $$

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot \left[ \rho \mathbf{u}\mathbf{u} + \left(P + \frac{B^2}{2}\right)\mathbf{I} - \mathbf{B}\mathbf{B} \right] = -\rho \nabla \Phi, \qquad (2) $$

$$ \frac{\partial \mathbf{B}}{\partial t} - \nabla \times [\mathbf{u} \times \mathbf{B} + \mathbf{E}_{\text{NIMHD}}] = 0, \qquad (3) $$

$$ \frac{\partial E_{\text{tot}}}{\partial t} + \nabla \cdot [(E_{\text{tot}} + P_{\text{tot}})\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\mathbf{B} - \mathbf{E}_{\text{NIMHD}} \times \mathbf{B}] = -\rho \mathbf{u} \cdot \nabla \Phi, \qquad (4) $$

where $\rho$ is the density, $\mathbf{u}$ the fluid velocity, $P$ the thermal pressure, $\mathbf{B}$ the magnetic field, and $\Phi$ the gravitational potential. $E_{\text{tot}}$ and $P_{\text{tot}}$ are the total energy and total pressure:

$$ E_{\text{tot}} = \rho\epsilon + \frac{1}{2}\rho u^2 + \frac{1}{2}B^2, \qquad (5) $$

$$ P_{\text{tot}} = (\gamma - 1)\rho\epsilon + \frac{1}{2}B^2, \qquad (6) $$

where $\epsilon$ is the specific internal energy and $\gamma$ the adiabatic index. $\mathbf{E}_{\text{NIMHD}}$ represents the electric field created by the three non-ideal MHD effects, and reads (Nakano et al. 2002; Masson et al. 2012)

$$ \mathbf{E}_{\text{NIMHD}} = -\eta_{\Omega}\mathbf{J} - \frac{\eta_{\text{H}}}{\|\mathbf{B}\|} (\mathbf{J} \times \mathbf{B}) + \frac{\eta_{\text{AD}}}{\|\mathbf{B}\|^2} ((\mathbf{J} \times \mathbf{B}) \times \mathbf{B}), \qquad (7) $$

with $\eta_{\Omega}$, $\eta_{\text{H}}$, and $\eta_{\text{AD}}$ the Ohmic, Hall, and ambipolar resistivities, and $\mathbf{J} = \nabla \times \mathbf{B}$ the electric current. The MHD equations have been written in the Gaussian cgs-units convention, such that the permeability constant $\mu_0 = 1$.

### 2.2. The Hall effect

We first performed a perturbation analysis to characterize the Hall effect. We ignored all forces but the Lorentz force, as well as the Ohmic and ambipolar diffusion terms. Let us consider a fluid at rest in an equilibrium state, threaded by a uniform magnetic field $\mathbf{B} = B\mathbf{e}_z$. We introduced small perturbations $\delta u_x$, $\delta u_y$, $\delta B_x$, $\delta B_y$ of the velocity and magnetic fields in the x- and y-directions, propagating along the z-direction at the frequency $\omega$. Each perturbation is therefore proportional to $\exp(i\omega t - ikz)$. Considering only the Hall effect, the induction Eq. (3) reads

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left[ \mathbf{u} \times \mathbf{B} - \frac{\eta_{\text{H}}}{\|\mathbf{B}\|} \mathbf{J} \times \mathbf{B} \right]. \qquad (8) $$

It is coupled with the momentum Eq. (2), with only the Lorentz force that we write in its non-conservative form:

$$ \rho \frac{\partial \mathbf{u}}{\partial t} = \mathbf{J} \times \mathbf{B}. \qquad (9) $$

We introduced the perturbations in Eqs. (8) and (9), and cancel out the zeroth- and second-order terms:

$$ \rho i \omega \delta u_x = -ikB \delta B_x, \qquad (10) $$

$$ \rho i \omega \delta u_y = -ikB \delta B_y, \qquad (11) $$

$$ i \omega \delta B_x = ikB \delta u_x + \eta_{\text{H}} k^2 \delta B_y, \qquad (12) $$

$$ i \omega \delta B_y = ikB \delta u_y - \eta_{\text{H}} k^2 \delta B_x. \qquad (13) $$

The system of equations is linear and can be written in matrix form:

$$ \begin{pmatrix} \rho i \omega & 0 & ikB & 0 \\ 0 & \rho i \omega & 0 & ikB \\ -ikB & 0 & i\omega & -\eta_{\text{H}}k^2 \\ 0 & -ikB & \eta_{\text{H}}k^2 & i\omega \end{pmatrix} \begin{pmatrix} \delta u_x \\ \delta u_y \\ \delta B_x \\ \delta B_y \end{pmatrix} = 0. \qquad (14) $$

A non-trivial solution is possible if the matrix determinant is zero. We obtain

$$ \omega^4 - \omega^2(2k^2c_A^2 + \eta_{\text{H}}^2k^4) + k^4c_A^4 = 0, \qquad (15) $$

with $c_A = \frac{B}{\sqrt{\rho}}$ the Alfvén speed. This equation can also be written as

$$ (\omega^2 - k^2 c_A^2)^2 - (\eta_{\text{H}} k^2 \omega)^2 = 0, \quad \text{then} \qquad (16) $$

$$ (\omega^2 - \eta_{\text{H}} k^2 \omega - k^2 c_A^2)(\omega^2 + \eta_{\text{H}} k^2 \omega - k^2 c_A^2) = 0. \qquad (17) $$

We only keep the positive solutions for $\omega$, which yields the following dispersion relation (Balbus & Terquem 2001):

$$ \omega = \pm \frac{\eta_{\text{H}} k^2}{2} + \sqrt{\left(\frac{\eta_{\text{H}} k^2}{2}\right)^2 + k^2 c_A^2}. \qquad (18) $$
Contrarily to the ambipolar and the Ohmic diffusions, the Hall effect is purely dispersive. Any perturbation creates two waves whose frequencies $\omega$ are linked to their wavelength $\lambda = 2\pi/k$. They are called whistler waves. They propagate the magnetic field perturbations along the field lines at the whistler speed:

$$c_w = \frac{\omega}{k} = \pm \frac{\eta_H k}{2} + \sqrt{\left(\frac{\eta_H k}{2}\right)^2 + c_A^2}. \qquad (19)$$

One of these two waves is always faster than the Afvén speed, and can theoretically travel at an unbound speed for very high frequencies and very small wavelengths. In practice, it is physically limited by the ion cyclotron frequency (Srinivasan & Shumlak 2011). From Eq. (19), we can define two regimes. If $\frac{\eta_H k}{2} \ll c_A$, then the Hall effect is weak and $c_w \approx c_A$. In this case, the whistler waves behave as Afvén waves. Otherwise, the Hall effect plays a significant role in the magnetic field evolution. This defines the Hall length $l_H = \frac{\eta_H k}{c_A}$ below which the Hall effect is dominant. Note that by being purely dispersive, the Hall effect does not directly impact the fluid energy or momentum. It only modifies the topology of the magnetic field.

### 3. Numerical methods – implementation of the Hall effect

We implemented the Hall effect in the Eulerian code RAMSES (Teyssier 2002). Since the Hall effect is dispersive, the introduction of new waves in the system modifies the Riemann problem that needs to be solved on each cell interface. For this reason we have adapted the method of Lesur et al. (2014) implemented in the PLUTO code (Mignone et al. 2007). This method is similar to Tóth et al. (2008) and is based on a direct modification of the Riemann solver. Therefore, we do not consider the Hall effect as a source term, as it is the case for the ambipolar and Ohmic diffusions in RAMSES (Masson et al. 2012).

#### 3.1. The constrained transport in RAMSES

In RAMSES, the magnetic field is naturally defined at the center of each orthogonal face, meaning that for the cell centered on indexes $(i, j, k)$ (for the $x$, $y$, and $z$ directions), $B_x$, $B_y$, and $B_z$ are defined at $(i - \frac{1}{2}, j, k)$, $(i, j - \frac{1}{2}, k)$, and $(i, j, k - \frac{1}{2})$, respectively. The code uses the constrained transport (CT; Evans & Hawley 1988) to solve both the induction equation of ideal MHD (Teyssier et al. 2006; Fromang et al. 2006) and the ambipolar and Ohmic diffusions (Masson et al. 2012). This method exploits the Stokes theorem by using the value of the electromotive forces (EMF) on cell edges to compute the magnetic flux through cell interfaces. For instance, Fig. 1 shows which EMFs are used to evolve the $x$-component of the magnetic field at the cell interface $(i - 1/2, j, k)$. The constrained transport update reads

$$\frac{B^{n+1}_{x,i-\frac{1}{2},j,k} - B^{n}_{x,i-\frac{1}{2},j,k}}{\Delta t} = \frac{E^{n+\frac{1}{2}}_{z,i-\frac{1}{2},j+\frac{1}{2},k} - E^{n+\frac{1}{2}}_{z,i-\frac{1}{2},j-\frac{1}{2},k}}{\Delta y} - \frac{E^{n+\frac{1}{2}}_{y,i-\frac{1}{2},j,k+\frac{1}{2}} - E^{n+\frac{1}{2}}_{y,i-\frac{1}{2},j,k-\frac{1}{2}}}{\Delta z}. \qquad (20)$$

In RAMSES, since $\Delta x = \Delta y = \Delta z$, we use only $\Delta x$ in the remainder of the manuscript. As RAMSES uses the MUSCL scheme (van Leer 1976), a predictive step is first performed to compute the flow variables from time-step $n$ to $n + 1/2$. These predictions

<figure><img src="image_3.png" /><figcaption>Fig. 1. Visualization of the constrained transport (CT) method for the calculation of the magnetic field $B_x$ at cell interface $(i - \frac{1}{2}, j, k)$.</figcaption></figure>

are then used to obtain the fluxes (here the EMFs). In the case of the CT, the EMFs are located at cell edges, at the intersection of four different states. In RAMSES, a 2D Riemann solver is used to deal with the discontinuity and compute the EMFs from these states. We explicitly detail these steps in the description of our implementation.

#### 3.2. Principle of the implementation

The Hall induction equation can be rewritten

$$\frac{\partial \mathbf{B}}{\partial t} = -\nabla \times [(\mathbf{u} + \mathbf{u}_H) \times \mathbf{B}], \qquad (21)$$

where $\mathbf{u}_H = -\eta_H \frac{\nabla \times \mathbf{B}}{\|\mathbf{B}\|}$ is the Hall speed (not to be confused with the whistler speed !), at which the Hall effect moves the field lines. Our goal is to replace the fluid velocity $\mathbf{u}$ by the “effective” velocity $\mathbf{u} + \mathbf{u}_H$ in the 2D Riemann solver that computes the EMFs on cell edges for the CT scheme.

The Hall effect is a dispersive term which induces the propagation of waves at every wavelength. The Shannon–Nyquist theorem (Shannon 1949) states that whistler waves with a wavelength smaller than twice the local resolution cannot be described and propagated. This defines the minimum resolved whistler wavelength as $\lambda = 2\Delta x$, where $\Delta x$ is the size of the local cell. Following Eq. (19), the associated whistler speed is

$$c_w = \frac{\eta_H \pi}{2 \Delta x} + \sqrt{\left(\frac{\eta_H \pi}{2 \Delta x}\right)^2 + c_A^2}. \qquad (22)$$

Since RAMSES is an explicit code in time, this wave speed needs to be accounted for in the computation of the characteristic speeds. $c_w$ however scales as $1/\Delta x$, which consequently decreases the spatial order of the scheme by one. To get a consistent implementation, an at least second-order scheme is hence necessary, which is the case here with the MUSCL scheme.

#### 3.3. The Hall speed

Since the CT scheme uses the EMFs on cell edges, the Hall speed $\mathbf{u}_H$ also needs to be computed on cell edges. We take the example of the edge located at $(i - \frac{1}{2}, j - \frac{1}{2}, k)$. The first step is to determine
the predicted states of the flow variables on cell faces and edges using a TVD slope. As part of the MUSCL scheme, the principle is to use a linear function whose slope is determined by considering neighboring cells. The TVD method (Harten 1983) limits this slope to prevent the generation of new, unwanted, extrema in the system. For example, a well-known slope limiter is minmod. For the variable $U$ in the cell with index $i$, the minmod slope is 0 if $U_i$ is a local extremum, and the lowest magnitude between $(U_i - U_{i-1})$ and $(U_{i+1} - U_i)$ otherwise. For the magnetic field, the predicted states are computed using the induction equation without the Hall term, coupled with the Stokes theorem:

$$ \frac{B_{x,i-\frac{1}{2},j,k}^{n+\frac{1}{2}} - B_{x,i-\frac{1}{2},j,k}^{n}}{\Delta t} = \frac{E_{z,i-\frac{1}{2},j+\frac{1}{2},k}^{n} - E_{z,i-\frac{1}{2},j-\frac{1}{2},k}^{n}}{\Delta x} - \frac{E_{y,i-\frac{1}{2},j,k+\frac{1}{2}}^{n} - E_{y,i-\frac{1}{2},j,k-\frac{1}{2}}^{n}}{\Delta x}, \quad (23) $$

$$ \frac{B_{y,i,j-\frac{1}{2},k}^{n+\frac{1}{2}} - B_{y,i,j-\frac{1}{2},k}^{n}}{\Delta t} = \frac{E_{x,i,j-\frac{1}{2},k+\frac{1}{2}}^{n} - E_{x,i,j-\frac{1}{2},k-\frac{1}{2}}^{n}}{\Delta x} - \frac{E_{z,i+\frac{1}{2},j-\frac{1}{2},k}^{n} - E_{z,i-\frac{1}{2},j-\frac{1}{2},k}^{n}}{\Delta x}, \quad (24) $$

$$ \frac{B_{z,i,j,k-\frac{1}{2}}^{n+\frac{1}{2}} - B_{z,i,j,k-\frac{1}{2}}^{n}}{\Delta t} = \frac{E_{y,i+\frac{1}{2},j,k-\frac{1}{2}}^{n} - E_{y,i-\frac{1}{2},j,k-\frac{1}{2}}^{n}}{\Delta x} - \frac{E_{x,i,j+\frac{1}{2},k-\frac{1}{2}}^{n} - E_{x,i,j-\frac{1}{2},k-\frac{1}{2}}^{n}}{\Delta x}, \quad (25) $$

where $E^n = u^n \times B^n$ is the electric field calculated on cell edges using simple interpolation of $u^n$ and $B^n$ at these locations, for example,

$$ \mathbf{u}_{i-\frac{1}{2},j-\frac{1}{2},k}^{n} = \frac{1}{4} \left( \mathbf{u}_{i,j,k}^{n} + \mathbf{u}_{i-1,j,k}^{n} + \mathbf{u}_{i-1,j-1,k}^{n} + \mathbf{u}_{i-1,j-1,k}^{n} \right), \quad (26) $$

$$ B_{x,i-\frac{1}{2},j-\frac{1}{2},k}^{n} = \frac{1}{2} \left( B_{x,i-\frac{1}{2},j,k}^{n} + B_{x,i-\frac{1}{2},j-1,k}^{n} \right), \quad (27) $$

$$ B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n} = \frac{1}{2} \left( B_{y,i,j-\frac{1}{2},k}^{n} + B_{x,i-1,j-\frac{1}{2},k}^{n} \right). \quad (28) $$

Only the perpendicular component of the magnetic field is calculated on each face. We reconstruct the other two by interpolating the state at the cell center using the TVD slope. For instance, obtaining $B_y$ on the face located at $(i - \frac{1}{2}, j, k)$ yields

$$ B_{y,i,j,k}^{n+\frac{1}{2}} = \frac{1}{2} \left( B_{y,i,j-\frac{1}{2},k}^{n+\frac{1}{2}} + B_{y,i,j+\frac{1}{2},k}^{n+\frac{1}{2}} \right), \quad (29) $$

$$ B_{y,i-\frac{1}{2},j,k}^{n+\frac{1}{2}} = B_{y,i,j,k}^{n+\frac{1}{2}} - \frac{\Delta_x(B_y)_{i,j,k}}{2}, \quad (30) $$

where $\Delta_x(B_y)_{i,j,k}$ is the slope of $B_y$ in the x-direction for cell $(i, j, k)$. As a result, two different magnetic fields are defined on each face, computed independently using the flow variables of the cells on each side of the interface. From here on, we have used the LI and RI exponents (left and right interfaces) to indicate the component coming from the lower index and higher index cells respectively, for example,

$$ \mathbf{B}_{i-\frac{1}{2},j,k}^{n+\frac{1}{2},\mathrm{LI}} = \mathbf{B}_{i-1,j,k}^{n+\frac{1}{2}} + \frac{\Delta_x(\mathbf{B})_{i-1,j,k}}{2}, \quad (31) $$

$$ \mathbf{B}_{i-\frac{1}{2},j,k}^{n+\frac{1}{2},\mathrm{RI}} = \mathbf{B}_{i,j,k}^{n+\frac{1}{2}} - \frac{\Delta_x(\mathbf{B})_{i,j,k}}{2}, \quad (32) $$

<figure><img src="image_4.png" /><figcaption>Fig. 2. Magnetic fields used to calculate J on edge $(i - \frac{1}{2}, j - \frac{1}{2}, k)$. The right picture represents the cell faces at $(i - \frac{1}{2})$, as highlighted on the 12 cells block $(i - 1 \cdots i, j - 1 \cdots j, k - 1 \cdots k + 1)$ on the left (the front right middle cell is cell $(i, j, k)$).</figcaption></figure>

with the same logic for y- and z- directions. We note that in each case, the component perpendicular to the face is naturally defined and does not need to be interpolated with the slope. This discontinuity of B is a problem here since the Hall speed is proportional to the electric current $J = \nabla \times B$. The jumps at cell interfaces may be assimilated to perturbations with an infinitely small wavelength that would propagate at an infinite speed, which is of course unphysical. To tackle this problem, we averaged B and consider that J is constant on cell faces and edges (Tóth et al. 2008; Lesur et al. 2014), at the price of a greater numerical diffusivity since we smoothed out discontinuities. We note Bpred the averaged magnetic field, therefore calculated as

$$ \mathbf{B}_{i-\frac{1}{2},j,k}^{\mathrm{pred}} = \frac{\mathbf{B}_{i-\frac{1}{2},j,k}^{n+\frac{1}{2},\mathrm{LI}} + \mathbf{B}_{i-\frac{1}{2},j,k}^{n+\frac{1}{2},\mathrm{RI}}}{2}. \quad (33) $$

Next we interpolated J on cell edges. For this, we also need the magnetic field at the two ends of the edges that we interpolated from Bpred. On the edge located at $(i - \frac{1}{2}, j - \frac{1}{2}, k)$, this yields

$$ \mathbf{B}_{i-\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}}^{\mathrm{int}} = \frac{\mathbf{B}_{i-\frac{1}{2},j-1,k}^{\mathrm{pred}} + \mathbf{B}_{i-\frac{1}{2},j-1,k+1}^{\mathrm{pred}} + \mathbf{B}_{i-\frac{1}{2},j,k+1}^{\mathrm{pred}} + \mathbf{B}_{i-\frac{1}{2},j,k}^{\mathrm{pred}}}{4}, \quad (34) $$

$$ \mathbf{B}_{i-\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}^{\mathrm{int}} = \frac{\mathbf{B}_{i-\frac{1}{2},j-1,k-1}^{\mathrm{pred}} + \mathbf{B}_{i-\frac{1}{2},j-1,k+1}^{\mathrm{pred}} + \mathbf{B}_{i-\frac{1}{2},j,k-1}^{\mathrm{pred}} + \mathbf{B}_{i-\frac{1}{2},j,k+1}^{\mathrm{pred}}}{4}, \quad (35) $$

Figure 2 depicts the magnetic fields used in the calculation of J on edge $(i - \frac{1}{2}, j - \frac{1}{2}, k)$. The faces represented are situated in $x = (i - \frac{1}{2})$ (see left picture).

We then computed J as $\nabla \times B$

$$ J_{x,i-\frac{1}{2},j-\frac{1}{2},k} = \frac{B_{z,i-\frac{1}{2},j,k}^{\mathrm{pred}} - B_{z,i-\frac{1}{2},j-1,k}^{\mathrm{pred}}}{\Delta x} - \frac{B_{y,i-\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}}^{\mathrm{int}} - B_{y,i-\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}^{\mathrm{int}}}{\Delta x}, \quad (36) $$

$$ J_{y,i-\frac{1}{2},j-\frac{1}{2},k} = \frac{B_{x,i-\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}}^{\mathrm{int}} - B_{x,i-\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}}^{\mathrm{int}}}{\Delta x}, $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Scheme of a 2D Riemann problem in the z-direction. The blue area represents the states at the corner for each cell at $i - \frac{1}{2}$, $j - \frac{1}{2}$, $k$.</figcaption></figure>

$$-\frac{B^{\text{pred}}_{z,i,j-\frac{1}{2},k} - B^{\text{pred}}_{z,i-1,j-\frac{1}{2},k}}{\Delta x}, \qquad (37)$$

$$J_{z,i-\frac{1}{2},j-\frac{1}{2},k} = \frac{B^{\text{pred}}_{y,i,j-\frac{1}{2},k} - B^{\text{pred}}_{y,i-1,j-\frac{1}{2},k}}{\Delta x} - \frac{B^{\text{pred}}_{x,i-\frac{1}{2},j,k} - B^{\text{pred}}_{x,i-\frac{1}{2},j-1,k}}{\Delta x}. \qquad (38)$$

The Hall resistivity $\eta_H$ depends on the density $\rho$, the temperature $T$, the cosmic-ray ionization rate $\zeta$, and the magnetic field amplitude $B$. These quantities are interpolated on cell edges from the cell-centered predicted states. For each of these variables $A$, we computed

$$\bar{A}_{i-\frac{1}{2},j-\frac{1}{2},k} = \frac{1}{4} (A_{i,j,k}^{n+\frac{1}{2}} + A_{i-1,j,k}^{n+\frac{1}{2}} + A_{i-1,j-1,k}^{n+\frac{1}{2}} + A_{i,j-1,k}^{n+\frac{1}{2}}). \qquad (39)$$

Finally, the Hall speed is

$$\mathbf{u}_{H,i-\frac{1}{2},j-\frac{1}{2},k} = -\frac{\eta_H(\bar{\rho}, \bar{T}, \bar{\zeta}, \bar{B})}{\bar{B}} \mathbf{J}, \qquad (40)$$

with every quantity calculated at $(i - \frac{1}{2}, j - \frac{1}{2}, k)$.

### 3.4. Solving the Riemann problem

We then solved the 2D Riemann problem to obtain the value of the total electromotive forces on cell edges. The notations are defined in Fig. 3. The blue area represents the corner states of the four cells at $(i - \frac{1}{2}, j - \frac{1}{2}, k)$, as calculated in Eqs. (41)–(52), and the red label is the emf computed in Eq. (53).

The flow variables were extrapolated on edges from the predicted states at cell centers using Eqs. (1)–(4) and the TVD slopes. Similar to the magnetic field on cell faces, this process yields four different states at each corner:

$$U_{i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},LB} = U_{i-1,j-1,k}^{n+\frac{1}{2}} + \frac{\Delta_x(U)_{i-1,j-1,k}}{2} + \frac{\Delta_y(U)_{i-1,j-1,k}}{2}, \qquad (41)$$

$$U_{i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RB} = U_{i-1,j-1,k}^{n+\frac{1}{2}} - \frac{\Delta_x(U)_{i,j-1,k}}{2} + \frac{\Delta_y(U)_{i,j-1,k}}{2}, \qquad (42)$$

$$U_{i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},LT} = U_{i-1,j,k}^{n+\frac{1}{2}} + \frac{\Delta_x(U)_{i-1,j,k}}{2} - \frac{\Delta_y(U)_{i-1,j,k}}{2}, \qquad (43)$$

$$U_{i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RT} = U_{i,j,k}^{n+\frac{1}{2}} - \frac{\Delta_x(U)_{i,j,k}}{2} - \frac{\Delta_y(U)_{i,j,k}}{2}, \qquad (44)$$

for any flow variable $U$, and the magnetic field is averaged on the interfaces:

$$B_{x,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},LB} = \frac{B_{x,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RB} + B_{x,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RB}}{2} \\ = \frac{\left[ B_{x,i-\frac{1}{2},j-1,k}^{n+\frac{1}{2},LI} + \frac{\Delta_y(B_x)_{i-1,j-1,k}}{2} \right] + \left[ B_{x,i-\frac{1}{2},j-1,k}^{n+\frac{1}{2},RI} + \frac{\Delta_y(B_x)_{i-1,j-1,k}}{2} \right]}{2}, \qquad (45)$$

$$B_{x,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},T} = \frac{B_{x,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},LT} + B_{x,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RT}}{2} \\ = \frac{\left[ B_{x,i-\frac{1}{2},j,k}^{n+\frac{1}{2},LI} - \frac{\Delta_y(B_x)_{i-1,j,k}}{2} \right] + \left[ B_{x,i-\frac{1}{2},j,k}^{n+\frac{1}{2},RI} - \frac{\Delta_y(B_x)_{i-1,j,k}}{2} \right]}{2}, \qquad (46)$$

$$B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},L} = \frac{B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},LB} + B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},LT}}{2} \\ = \frac{\left[ B_{y,i-1,j-\frac{1}{2},k}^{n+\frac{1}{2},LI} + \frac{\Delta_x(B_y)_{i-1,j-1,k}}{2} \right] + \left[ B_{y,i-1,j-\frac{1}{2},k}^{n+\frac{1}{2},RI} + \frac{\Delta_x(B_y)_{i-1,j-1,k}}{2} \right]}{2}, \qquad (47)$$

$$B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},R} = \frac{B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RB} + B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RT}}{2} \\ = \frac{\left[ B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},LI} - \frac{\Delta_x(B_y)_{i-1,j-1,k}}{2} \right] + \left[ B_{y,i-\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2},RI} - \frac{\Delta_x(B_y)_{i-1,j-1,k}}{2} \right]}{2}. \qquad (48)$$

We write these quantities $B_x^B$, $B_x^T$, $B_y^L$, and $B_y^R$. The EMFs were then computed at each corner state, with the addition of the Hall speed (40) (that is the same for the four states), so that $\mathbf{E} = (\mathbf{u} + \mathbf{u}_H) \times \mathbf{B}$. On the $(i - \frac{1}{2}, j - \frac{1}{2}, k)$ edge, we needed the z-component:

$$E_z^{n+\frac{1}{2},LB} = (u_x^{n+\frac{1}{2},LB} + u_{H,x})B_y^L - (u_y^{n+\frac{1}{2},LB} + u_{H,y})B_x^B, \qquad (49)$$

$$E_z^{n+\frac{1}{2},LT} = (u_x^{n+\frac{1}{2},LT} + u_{H,x})B_y^L - (u_y^{n+\frac{1}{2},LT} + u_{H,y})B_x^T, \qquad (50)$$

$$E_z^{n+\frac{1}{2},RB} = (u_x^{n+\frac{1}{2},RB} + u_{H,x})B_y^R - (u_y^{n+\frac{1}{2},RB} + u_{H,y})B_x^B, \qquad (51)$$

$$E_z^{n+\frac{1}{2},RT} = (u_x^{n+\frac{1}{2},RT} + u_{H,x})B_y^R - (u_y^{n+\frac{1}{2},RT} + u_{H,y})B_x^T, \qquad (52)$$

The resulting EMF on the edge was obtained by solving the Riemann problem at this discontinuity of four different states. The Hall effect is dispersive and introduces two waves into the system, for a total of 10:

- Two Alfvén waves that propagate the perturbations of the magnetic field, $\lambda = u \pm c_A$.
- Two fast magnetosonic waves, representing the correlation between **B** and $\rho$ by the coupling between Lorentz force and thermal pressure,

$$ \lambda = u \pm \sqrt{\frac{1}{2}(c_s^2 + c_A^2) + \frac{1}{2}\sqrt{(c_s^2 + c_A^2) - 4c_s^2 c_{A,x}^2}} $$

(with $c_{A,x} = B_x / \sqrt{\rho}$ the component of the Alfvèn speed in the direction of propagation and $c_s = \sqrt{\gamma P / \rho}$ the sound speed).

- Two slow magneto sonic waves, representing the anti-correlation between **B** and $\rho$ by the coupling between Lorentz force and thermal pressure,

$$ \lambda = u \pm \sqrt{\frac{1}{2}(c_s^2 + c_A^2) - \frac{1}{2}\sqrt{(c_s^2 + c_A^2) - 4c_s^2 c_{A,x}^2}} $$

- One contact discontinuity wave $\lambda = u$.

- Two dispersive whistler waves.

Given the high number of waves, we used approximate Riemann solvers to compute the flux at cell faces and the EMF on cell edges. Tóth et al. (2008) demonstrated that the Lax-Friedrich Riemann solver is inconsistent with the Hall effect, because the whistler speed scales as $c_w \sim \Delta x^{-1}$. We chose to use the HLL solver (Harten et al. 1983), which considers only the lowest and highest (algebraic) wave speeds in each direction, and a constant state is assumed between these waves. It is a satisfying compromise between accuracy and numerical complexity, compared, for instance, to the HLLD solver (Miyoshi & Kusano 2005) which takes five wave speeds into account. In the case of a 2D Riemann problem, the HLL solution for the EMF is (Londrillo & del Zanna 2004)

$$ E_z^{n+\frac{1}{2}} = \frac{S_L S_B E_z^{n+\frac{1}{2},LB} - S_L S_T E_z^{n+\frac{1}{2},LT} - S_R S_B E_z^{n+\frac{1}{2},RB} + S_R S_T E_z^{n+\frac{1}{2},RT}}{(S_R - S_L)(S_T - S_B)} - \frac{S_T S_B}{S_T - S_B}(B_x^T - B_x^B) + \frac{S_R S_L}{S_R - S_L}(B_y^R - B_y^L), \quad (53) $$

where $S_L, S_R, S_B$, and $S_T$ are the lowest and highest wave speeds in the x- and y-directions. Once this work has been done for every edge, the magnetic field is updated with the constrained transport method, using Eq. (20).

### 3.5. The speed of whistler waves

In order to correctly propagate the perturbations caused by the Hall effect, the whistler wave speed (22) needs to be taken into account in the Courant Friedrichs Lewy (CFL) condition (Courant et al. 1928). This condition ensures that information does not travel further than one cell in one time-step, otherwise it would be lost since information exchange happens between adjacent cells. At high resolution, we have roughly $c_w \approx \frac{\eta_{hi} \pi}{\Delta x}$, and the corresponding time-step

$$ \Delta t_{\text{Hall}} \approx \frac{\Delta x}{c_w} \propto \Delta x^2. \quad (54) $$

The scaling of the time-step in $\Delta x^2$ induces a fast decrease of the time-step as resolution increases, which is a limitation of current numerical simulation codes. Techniques that allow a faster integration such as Super Time-Stepping cannot be employed for the Hall effect because a larger time-step would cut off the high-frequency whistler waves.

At high resolutions, the fast whistler waves have the fastest characteristic speeds of the Riemann problem. We then have

$S_L = u_x = c_w, S_R = u_x + c_w, S_B = u_y - c_w$ and $S_T = u_y + c_w$. This impacts the 2D Riemann problem (as described in the implementation (Sect. 3.4)), as well as the 1D Riemann problems that are solved at cell interfaces during the prediction step of the MUSCL scheme.

### 3.6. Summary

The different steps of the Hall resolution are in the following order:

- Compute the predicted states of the flow variables (Eqs. (23)–(28)).

- Interpolate **B** on cell faces and corners with Eqs. (29)–(35).

- Calculate the electric current **J** on edges with Eqs. (36)–(38).

- Calculate the Hall speed on edges with the averaged flow variables and the interpolated **J** (Eqs. (39)–(40)).

- Extrapolate the flow variables on cell corners (Eqs. (41)–(48)).

- Use the effective speed **u** + **u**<sub>H</sub> instead of **u** to compute the EMFs at cell corners (Eqs. (49)–(52)).

- Solve the 2D Riemann problem by considering the whistler waves in the calculation of $S_L, S_R, S_B$, and $S_T$, to get the electric fields $E^{n+1/2}$ with Eq. (53).

- Update the flow variables with the constrained transport scheme (20).

### 3.7. Comparison to other implementation

Our implementation of the Hall effect in RAMSES followed the work of Lesur et al. (2014), who included it in the Eulerian code PLUTO (Mignone et al. 2007). While both codes use the CT algorithm, Lesur et al. (2014) use a 1D Riemann solver to compute the EMF while we have included the Hall effect into 2D Riemann solver, which, to the best of our knowledge, has never been done before. Both codes use a second-order scheme (Runge-Kutta for PLUTO), and share the same convergence properties for the Hall effect. The approach is also similar to the implementation of Tóth et al. (2008) for the BASTRUS code (Tóth et al. 2006), although they implemented the Hall effect using an explicit time integration within an implicit code. Since implicit schemes are unconditionally stable, they are not strictly bound to the whistler speed to choose their time-step. This emancipation results in a faster code because larger time-steps are allowed. Furthermore, without the constraint of the whistler speed, the scheme order is increased by one. The drawback of the BASTRUS code is a less precise scheme at small scales because the shortest whistler waves cannot be correctly propagated.

The Hall effect has been implemented in the Eulerian code ZeusTW (Stone & Norman 1992; Krasnopolsky et al. 2010) according to Sano & Stone (2002). The term $(\nabla \times \mathbf{B}) \times \mathbf{B}$ is evolved at $t^{n+1/2}$ using a finite difference scheme to evaluate the curls in an explicit way. The resulting Hall-EMF is then added to the ideal EMF $\mathbf{u} \times \mathbf{B}$ into a modified constrained transport method. The tests in Sano & Stone (2002) show that the dispersion relation of the Hall effect is verified, but we have not found any information about the accuracy of the scheme.

Due to their nature, Lagrangean codes cannot use the constrained transport method. In the SPH code PHANTOM (Wurster et al. 2014; Price et al. 2017) as well as in Tsukamoto et al. (2015, 2017), the three non-ideal MHD terms were implemented as source terms. Testing them with the ambipolar diffusion shows second-order convergence, but no information is given for the Hall effect, nor the propagation of whistler waves
at various resolution. They do, however, recover the theoretical dispersion relation, and perform the shock-test described in Sect. 4.2 with results fitting the expected profile within 3% (Wurster et al. 2016).

## 4. Tests

In the following section, we describe how we tested our implementation of the Hall effect, especially the propagation of the whistler waves. Except in the shock test, $E_{NIMHD}$ is always reduced to the Hall effect.

### 4.1. Propagating whistler waves

#### 4.1.1. Dispersion relation test

We tested the propagation of whistler waves with our scheme by following the setup described in Kunz & Lesur (2013). First, we aimed to retrieve the dispersion relation (18). In a cubic box of length 1 cm with uniform density and pressure, we set a whistler wave with a wave-number $k$, and study its propagation in the $x$-direction and its oscillations. The initial MHD quantities are

$$ 
\begin{aligned}
\rho &= 1 \text{ g cm}^{-3}, & B_x &= 0.1 \text{ G}, \\
P &= 1.5 \times 10^{-5} \text{ dyne cm}^{-2}, & B_y &= 10^{-3} \cos(2\pi kx) \text{ G}, \\
\gamma &= \frac{5}{3}, & B_z &= 0 \text{ G}, \\
\mathbf{u} &= 0 \text{ cm s}^{-1}.
\end{aligned}
$$

We imposed periodic boundary conditions, which requires $k$ to be an integer to ensure continuity. To retrieve the dispersion relation (18), we vary $k$ and $\eta_H$, typically between 5 and 20 for $k$ and $5 \times 10^{-4}$ to $0.1 \text{ cm}^2 \text{ s}^{-1}$ for $\eta_H$. We used a uniform resolution of $\Delta x = 1/128$ cm and the moncen slope limiter. The question of resolution and slope limiters are discussed in the next section.

Figure 4 shows the time evolution of $u_y$ and $B_y$ in the case $k=7$ and $\eta_H = 5 \times 10^{-3} \text{ cm}^2 \text{ s}^{-1}$, for a random cell in the computational domain (here $x=y=32.5/128$, $z=63.5/128$). As expected, two whistler waves propagate with different frequencies. The fluid is set into motion at the same frequencies as the magnetic field. Figure 5 shows the Fourier transform of the temporal evolution of the magnetic field in the $y$-direction for the same cell as Fig. 4. It contains two peaks corresponding to the frequencies at which the waves propagate. These peaks are visible around $\omega \approx 2 \text{ s}^{-1}$ and $\omega \approx 0.3 \text{ s}^{-1}$. By varying the Hall resistivity $\eta_H$ and the wave-number $k$, we obtain a set of frequencies. Figure 6 shows the frequencies, normalized by the Hall frequency $\omega_H = \frac{c^2}{\eta_H}$, plotted as a function of $kl_H$. The agreement between our results and the theoretical dispersion relation is very good (within 5% for most points, see blue points in Fig. 6). However, we could not extract good approximations of the low frequency waves at high $kl_H$ because of the fast dissipation of the signal (see Sect. 4.1.2). For a similar reason, the eigenfrequencies of $kl_H < 1$ waves are less precise because of their slow oscillations.

#### 4.1.2. Truncation errors and order of the scheme

Numerical schemes are based on the approximation of Taylor series, which creates truncation errors. In the context of wave propagations, they take the form of an artificial dissipation of the signal. For a well-designed numerical scheme, this error should decrease when the resolution increase. The damping of

<figure><img src="image_6.png" /><figcaption>Fig. 4. Typical evolution of $u_y$ (in blue) and $B_y$ (in red) through time for a random cell in the box (here $x = 48.5/128$, $y = 0.5/128$, $z = 63.5/128$). $k = 7 \text{ cm}^{-1}$, $\eta_H = 5 \times 10^{-3} \text{ cm}^2 \text{ s}^{-1}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Fourier transform of $B_y(t)$ for a random cell in the box (same as Fig. 4).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Correspondence between the theoretical dispersion relation of the whistler waves (solid line), and the main oscillation frequencies of the magnetic field (circles). The blue square points represent the relative error.</figcaption></figure>

the signal clearly appears in Fig. 4 with the decrease of the wave amplitudes.

The transverse magnetic field can then be estimated as

$$ B_y(t) = f(t)e^{-at}, \qquad (55) $$
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Numerical damping rates as a function of resolution for different slope limiters, and $\Delta x$ and $\Delta x^2$ scalings for comparison.</figcaption></figure>

where $f(t)$ is the periodic contribution and $e^{-at}$ the damping, with $a$ the damping coefficient in s$^{-1}$. Ideally, this coefficient needs to be as small as possible. To achieve this, we wanted to make sure that our scheme is at the highest possible order in space, meaning, at least second order in our implementation. We were able to control the convergence of the calculation by changing the grid resolution, as well as the slope limiter as suggested in Tóth et al. (2008) and Lesur et al. (2014). We use the same setup as previously, with $k = 5$ cm$^{-1}$, $\eta_H = 0.005$ cm$^2$ s$^{-1}$. We also used the slope limiters minmod, moncen, and generalized moncen ($\beta = 1.5$), for resolutions of $\Delta x = \frac{1}{32}$, $\frac{1}{64}$, and $\frac{1}{128}$ cm. The resulting damping rates $a$ are plotted in Fig. 7. Similar to Tóth et al. (2008) and Lesur et al. (2014), we find that the minmod limiter is only first order, while the moncen and generalized moncen limiters show second-order accuracy.

A second-order scheme like MUSCL, combined with the appropriate slope limiter, appears to be necessary to treat the Hall effect. As already mentioned in Sect. 3.2, the characteristic whistler speed scales as $1/\Delta x$, which lowers the order of the scheme by one. The test results show that without a TVD slope, the numerical dissipation does not decrease with the resolution, and the scheme is inconsistent.

Tóth et al. (2008) show that for a second-order slope limiter, the numerical dissipation scales as $\frac{\partial^4 U}{\partial x^4}$. It means that at a given number of points per $l_H$, the damping rate increases as $k^4$. We evaluated the damping rate for $\eta_H = 0.005$ cm$^2$ s$^{-1}$ at a fixed resolution of $\Delta x = \frac{1}{128}$ cm and $k = 5, 7, 8, 10, 20$ cm$^{-1}$. We used the moncen slope limiter. Figure 8 shows the evolution of the numerical damping rate as a function of the wave-number $k$. We recovered the theoretical scaling in $k^4$, which indicates that smaller wavelengths are substantially more affected by the numerical dissipation than the larger ones.

### 4.1.3. Grid refinement

Here we investigate the behavior of whistler waves on non-regular grids, using the same plane-parallel setup of the previous tests, while refining the grid on specific locations. We once more took $k = 5$ cm$^{-1}$ and $\eta_H = 0.005$ cm$^2$ s$^{-1}$ as the reference case, with the moncen slope limiter. We typically used two or three different levels of refinement, with the maximum resolution centered on $x = 0$ and gradually decreasing on each side. The simulation box extends between $x = -0.5$ cm and $x = +0.5$ cm. A refinement level $\ell$ corresponds to a resolution $\Delta x = 1/2^\ell$ cm. For the two-levels tests, the levels of resolution are $\ell = \ell_{min} + 1$

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Numerical damping rates versus wave-number $k$. Solid line: numerical results, dashed line: $k^4$ scaling.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Damping coefficient of the whistler wave as a function of the maximum (triangles) and minimum (squares) resolution of the grid. Colors indicate the minimum resolution: purple for $\ell_{min} = 4$, blue for $\ell_{min} = 5$, green for $\ell_{min} = 6$, and red for $\ell_{min} = 7$. The circle points represent the regular grids, and the dashed line the least-squares fit of these points ($\propto \Delta x^{2.11}$).</figcaption></figure>

for $-0.25 < x < 0.25$, and $\ell = \ell_{min}$ everywhere else. For three-levels tests, we impose $\ell = \ell_{min} + 2$ between $-0.125 < x < 0.125$, $\ell = \ell_{min} + 1$ in $-0.25 < x < 0.25$, and $\ell = \ell_{min}$ in the rest of the box. We chose $\ell_{min} = 5, 6$ for the two-levels simulations and $\ell_{min} = 4, 5, 6$ for the three-levels. The whistler wave propagates along the x-axis.

Figure 9 shows the damping rates of the whistler waves as a function of the maximum (triangles) and minimum (squares) resolution for the five runs described above, and for regular grids (circles) with minimum refinement levels from four to seven (colors). The damping in regular grids follows the same $\Delta x^2$ scaling as in the previous test. However, non-regular grids produce a higher damping rate that does not always scale down with the maximum resolution. Instead, the damping rate seems to be only affected by the minimum resolution. While the AMR is an efficient tool to capture small-scale phenomena at a relatively low cost, it does not seem to be adapted for the propagation of whistler waves, whose dissipation depends only on the coarser resolution they cross.

### 4.2. Shock test

We then tested the Hall effect in a more physical situation. We considered a unidimensional isothermal shock in a
<table><caption>**Table 1.** Initial and boundary flow variables for the shock test.</caption><thead><tr><th>Variable</th><th>Left state</th><th>Right state</th></tr></thead><tbody><tr><td>ρ (g cm<sup>-3</sup>)</td><td>1.7942</td><td>1.0</td></tr><tr><td>P (dyn cm<sup>-3</sup>)</td><td>0.017942</td><td>0.01</td></tr><tr><td>u (cm s<sup>-1</sup>)</td><td>-0.9759</td><td>-1.751</td></tr><tr><td></td><td>-0.6561</td><td>0.0</td></tr><tr><td></td><td>0.0</td><td>0.0</td></tr><tr><td></td><td>1.0</td><td>-1.0</td></tr><tr><td></td><td>1.74885</td><td>0.6</td></tr><tr><td></td><td>0.0</td><td>0.0</td></tr></tbody></table>

Hall-dominated regime in presence of the Ohmic and the ambipolar diffusions, following the setup of Falle (2003) and Wurster et al. (2016). The sound speed and x-magnetic field are uniform, with $c_s = \sqrt{\rho} = 0.1$ cm s$^{-1}$, $B_x = 1$ G, and we used the resistivities $\eta_\Omega = 10^{-9}$, $\eta_H = -2 \times 10^{-2}$ B, and $\eta_{AD} = 3.5 \times 10^{-3} \frac{B^2}{\rho}$ cm$^2$ s$^{-1}$. We considered a cubic box with a discontinuity in the x-direction, and characterized by the left and right states presented in Table 1.

Initially, the left state spans $x \in [0; 0.3]$ while the region $x \in [0.3; 1]$ is in the right state. We put the discontinuity at $x = 0.3$ to let the fluid oscillate in the low density region. The boundaries are periodic in the y- and z-directions, and are fixed with the left and right states in the x-direction. The box has a 1 cm size. We used the AMR in order to gain computation time while keeping a good precision in the sensitive zones. The minimum resolution is 1/64 cm (level $l_{min} = 6$ of AMR), and we refine the grid to keep the magnetic field variation lower than 5% in one cell, up to a resolution of 1/256 cm (level eight of AMR). We let the system evolve until the initial discontinuity converges toward a steady-state profile of the shock. The analytical calculation of the theoretical profile is detailed in Appendix A.

Figure 10 displays both the numerical steady-state profiles and the analytical solution. The agreement is very good with a difference lower than 1% at shock interface. The larger errors in outer regions are oscillations due to the low resolution and the fixed boundary conditions, which happen even without the Hall effect.

## 5. Collapse of a dense core

In this section, we describe our simulations of dense core collapse with the Hall effect as a practical test of our implementation. Simulations are presented in Sects. 5.2 and 5.3, following elements of Hall effect theory in this context in Sect. 5.1.

### 5.1. Hall effect in star formation: theory

Assuming axisymmetry and accounting only for the Hall effect, the azimuthal component of Eq. (8) reads

$$ \begin{aligned} \frac{\partial B_{\theta}}{\partial t} &= -\left[\nabla \times \left(\frac{\eta_{\mathrm{H}}}{\|\mathbf{B}\|} (\nabla \times \mathbf{B}) \times \mathbf{B}\right)\right]_{\theta} \\ &= -\frac{\eta_{\mathrm{H}}}{B} \left[ \frac{\partial}{\partial z} \left[ B_z \left( \frac{\partial B_r}{\partial z} - \frac{\partial B_z}{\partial r} \right) - \frac{B_{\theta}}{r} \frac{\partial(r B_{\theta})}{\partial r} \right] \right. \\ & \qquad \left. + \frac{\partial}{\partial r} \left[ B_{\theta} \frac{\partial B_{\theta}}{\partial z} + B_r \frac{\partial B_r}{\partial z} + B_r \frac{\partial B_z}{\partial r} \right] \right]. \end{aligned} \quad (56) $$

This expression shows that Hall effect generates a toroidal magnetic field $B_\theta$, provided that its radial and vertical components $B_r$ and $B_z$ are not uniform.

During the collapse of a magnetized dense core, the pinching of the field lines creates an accumulation of matter known as the “pseudo-disk” (Galli & Shu 1993), with a density typically higher than $10^{-15}$ g cm$^{-3}$. It is perpendicular to the large-scale magnetic fields and is located where the pinching is the strongest, at the neck of the “hourglass” formed by the field lines. The large spatial variations of the magnetic field at this location generates a strong toroidal electric current in the mid-plane of the pseudo-disk. $\mathbf{u}_H$ being directly proportional to $\mathbf{J}$, the field lines are then twisted in the plane of the pseudo-disk. If the core is initially rotating, the fluid also twists the lines. The Hall effect therefore strengthens or weakens the toroidal magnetic field generated by the rotation, depending on the sign of $-\eta_H J_\theta$. The consequence is a stronger or weaker magnetic braking, resulting in a slower or faster rotation speed of the fluid (see, e.g., Tsukamoto et al. 2015). If the core is not rotating initially, the Hall effect still generates a toroidal magnetic field, which induces the rotation of the fluid, in the same way that the magnetic braking slows the rotation via the magnetic tension. In either case, the fluid is accelerated in the direction of $-\mathbf{u}_H$.

### 5.2. Setup

The objective is to test our implementation with a simple core-collapse simulation, and study the role of the Hall effect. We used the standard Boss & Bodenheimer (1979) initial conditions with a non-rotating uniform sphere of gas. The sphere radius is 3712 au and its density $\rho_0 = 4.15 \times 10^{-18}$ g cm$^{-3}$ for a total mass of 1.5 $M_\odot$. The thermal over gravitational energy ratio is $\alpha = 0.25$ with a uniform temperature $T = 10$ K. The external medium is 100 times less dense. The magnetic field is uniform at $B = 90.3 \ \mu\mathrm{G}$ in the whole box along the z-direction, which is equivalent to a mass-to-flux ratio of $\mu = 7$. We solved the complete set of Eqs. (1)–(4) without considering the ambipolar and the Ohmic diffusions. The Hall resistivity is set at a fixed value of $\eta_H = 10^{20} \mathrm{cm}^2 \mathrm{s}^{-1}$. The length of the simulation box is four times the sphere radius and we use the generalized moncen slope limiter with a 1.5 coefficient. This setup is now refereed as our “reference case”. To mimic the evolution of temperature during the collapse up to the first Larson core, we used the following barotropic equation of state:

$$ T = 10 \left[ 1 + \left( \frac{\rho}{10^{-13} \mathrm{g} \mathrm{cm}^{-3}} \right)^{\frac{5}{3}} \right]. \quad (57) $$

The initial grid is refined at AMR level 5 ($32^3$) and we imposed at least eight points per Jeans length to satisfy the Truelove et al. (1997) condition of at least four points per Jeans length. The maximum resolution is $\Delta x = 2$ au at the end of the simulation. The “formation of the first core” is the time $t_C$ at which the maximum density reaches $10^{-13}$ g cm$^{-3}$. Finally, we refer to a structure rotationally supported against gravity as “disk”. We used the criteria of Joos et al. (2012):

- density criterion: $\rho > 3.8 \times 10^{-15}$ g cm$^{-3}$,

- rotationally supported: $u_\theta > 2u_r$, with $u_\theta$ and $u_r$ the azimuthal and radial velocities, respectively,

- hydrostatic equilibrium within the disk: $u_\theta > 2u_z$,

- not within the first core: $\rho u^2/2 > 2P_{\mathrm{therm}}$, with $P_{\mathrm{therm}}$ the thermal pressure.
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Comparison between the analytical solution of the standing shock (solid black line) and the numerical result (red squares). The AMR level is represented as the dashed black line (right axis). *From top to bottom panels:* $\rho$ and $u_x$, $u_y$ and $u_z$, $B_y$ and $B_z$.</figcaption></figure>

5.3. Results

A differential rotation of the cloud starts about an axis in the
same direction as the initial magnetic field, with an increased

amplitude over time. Figure 11 shows maps on the plane $x = 0$
of the ratio of the toroidal magnetic field over the total magnetic
field amplitude $B_{\theta}/B$ and of the azimuthal velocity $u_{\theta}$, 0.230 kyr
after the formation of the first core. The twisting of the field
<figure><img src="image_13.png" /><figcaption>Fig. 11. Toroidal magnetic field $B_\theta$ (top panel) and azimuthal velocity $u_\theta$ (bottom panel) on the plane $x = 0$, 0.230 kyr after the formation of the first core.</figcaption></figure>

lines is the strongest 200 au above and below the mid-plane. A large quantity of angular momentum is created in the disk while the whistler waves transport its negative counter-part along the field lines. This negative angular momentum is transferred to the regions above and below the disk, thus setting a backward rotation that is visible in blue colors in the bottom panel. The momentum of the core and disk are therefore compensated by the counter-rotation of an hourglass-shaped envelope of 1000 au in height and diameter. With a similar setup, Krasnopolsky et al. (2011) and Li et al. (2011) obtain the same qualitative results, with counter-rotating regions in the envelope.

We then compared the azimuthal velocity $u_\theta$ of the fluid to the azimuthal Hall speed and to the Keplerian velocity. The Hall speed is defined accordingly in Sect. 3.2:

$$ \mathbf{u}_H = \eta_H \frac{\mathbf{J}}{\|\mathbf{B}\|}. \qquad (58) $$

Using the second Newton law, the Keplerian velocity was calculated as

$$ u_{\text{Kepler}}(\mathbf{r}) = \left[ r \left( \sum_{i=1, \text{ncells}} \frac{G m_i}{\|\mathbf{r}_i - \mathbf{r}\|^3} (\mathbf{r}_i - \mathbf{r}) \right) \right]^{\frac{1}{3}} \qquad (59) $$

for each cell in our discretized system (assuming the center of mass is at the origin).

Figure 12 shows the radial profile of the three velocities 0.500 and 0.800 kyr after the first core formation. Figure 13 is a zoom-out of Fig. 12 without the Keplerian velocity. All radial profiles are azimuthally averaged. In the inner regions, the fluid rotation speed is equal to the Keplerian velocity and larger than the Hall speed. Eventually, the centrifugal force overcomes gravity, which results in a positive radial velocity of the fluid and the formation of a ring-shaped core, as shown in Fig. 14. All the velocities drop outside the Larson core, beyond 30 au, and the fluid is roughly rotating at the Hall velocity, as seen in Fig. 13. The two values do not perfectly match but differ by less than 50%. In the absence of every other forces but the Lorentz force, the two velocities should correspond, but in our case, the presence of gravity and thermal pressure prevent a perfect match. The coupling is effective up to $r \approx 160$ au, down to a density of $\rho \approx 10^{-14}$ g cm$^{-3}$. As explained in Sect. 5.1, the region we could define as the pseudo-disk is rotating at the Hall velocity. Around 30 au at $t = t_c + 0.500$ kyr and 40 au at $t = t_c + 0.800$ kyr, the azimuthal Hall velocity changes sign in a small interval. This component of $u_H$ is proportional to $\mathbf{J}_\theta = \frac{\partial B_z}{\partial z} - \frac{\partial B_r}{\partial r}$. While $B_z$ decreases with $r$, there is a small interval at the outer edge of the core, before the accretion shock, in which its decrease is much lower. As a result, $\frac{\partial B_z}{\partial r}$ is lower (in negative value) while $\frac{\partial B_r}{\partial z}$ keeps an approximately constant value, causing a temporary change of sign.

The Hall effect redistributes the angular momentum within the envelope, and we know to address the case of magnetic flux. Figure 15 shows the magnetic field amplitude as a function of the density for each cell in the box for this simulation and its counter-part without the Hall effect (i.e., ideal MHD), when both simulations reach a maximum density of $\rho = 5 \times 10^{-12}$ g cm$^{-3}$. The magnetic field in ideal MHD is roughly 0.5 G at the center, while it is bounded to $3 \times 10^{-2}$ G with the Hall effect. On the contrary, the magnetic field amplitudes are very similar at densities below $10^{-15}-10^{-14}$ g cm$^{-3}$, confirming that the Hall effect becomes relevant above these densities. Masson et al. (2016) found a comparable “plateau” using the ambipolar diffusion, which also redistributes the magnetic flux and yields a similar profile for the magnetic field. Both non-ideal MHD effects then act to regulate the magnetic flux.

## 5.4. The issue of angular momentum conservation

### 5.4.1. Presentation of the problem

Here we introduce some notations that are useful to the remainder of the discussion. We note $\mathbf{L}$ the angular momentum vector

$$ \mathbf{L} = \sum_{i=1, \text{ncell}} m_i \mathbf{r}_i \times \mathbf{u}_i, \qquad (60) $$

where $m_i$, $\mathbf{r}_i$, and $\mathbf{u}_i$ represent the mass, the position vector, and the velocity vector of cell $i$. The system being nearly axisymmetric, the main component of this vector is $L_z$ (i.e., the initial magnetic field direction), that can also be calculated as

$$ L_z = \sum_{i=1, \text{ncell}} m_i r_{c,i} u_{\theta,i}, \qquad (61) $$

with $r_c$ the cylindrical radius and $u_\theta$ the azimuthal velocity. We define the “positive” and “negative” momenta $L_z^+$ and $L_z^-$.
<figure><img src="image_14.png" /><figcaption>Fig. 12. Radial profiles ($z = 0$, azimuthal averaging) of the azimuthal velocity (red curve), the azimuthal Hall speed (blue curve), and the Keplerian velocity (green curve) at $t = t_c + 0.500$ and $t = t_c + 0.800$ kyr (top and bottom panels). The black dotted line represents the density profile (right axis).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. As Fig. 12, at larger radius and without the Keplerian velocity.</figcaption></figure>

as follows:

$$L_z^+ = \sum_{i, v_{\theta,i} > 0} m_i r_{c,i} u_{\theta,i}, \qquad (62)$$

$$L_z^- = \sum_{i, v_{\theta,i} < 0} m_i r_{c,i} u_{\theta,i}. \qquad (63)$$

We therefore have $L_z^+ + L_z^- = L_z \approx \|\mathbf{L}\|.$

In our setup, there is initially no angular momentum in the system. $L_z^+$ thus represents the “positive” rotation created by the Hall effect, balanced by the “negative” rotation $L_z^-$, and the numerical solution should verify $\|L_z^+\| \approx \|L_z^-\|$. Usually, $L_z$ does not exactly equal zero because of truncation errors, but this is not an issue if $L_z$ remains small (typically, compared with $L_z^+$ and $L_z^-$).

In Fig. 16, we have plotted the evolution of $L_z^+$, $L_z^-$, and $L_z$ during the simulation after the formation of the first Larson core. In the remainder of the manuscript, this reference case is plotted in red. From $t = t_c + 0.230$ kyr, the positive momentum skyrockets and almost doubles its value by $t = t_c + 1.300$ kyr, while the negative momentum increases by less than 10%. The total angular momentum $L_z$ then increases as well and reaches 2/3 of the value of $L_z^-$. The angular momentum is therefore not conserved in the simulation box.

<figure><img src="image_16.png" /><figcaption>Fig. 14. Density and velocity map of the reference case at $t = t_c + 0.800$ kyr.</figcaption></figure>

Krasnopolsky et al. (2011) also encounter a non-conservation of angular momentum in their simulation with a similar setup. They attribute this phenomenon to torsional Alfvén waves leaving the simulation box. It is however not the case here since the boundary conditions are periodic. The blue curve of Fig. 16 represents the angular momentum of the disk and the core,
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Magnetic field amplitude as a function of the density for two simulations with and without Hall (red and blue domains, respectively). The snapshot is taken when the maximum density reaches $5 \times 10^{-12} \text{ g cm}^{-3}$ (at $t = t_C + 0.220$ kyr for the Hall simulation).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Positive $L_z^+$ and negative $L_z^-$ angular momentum evolution after the formation of the first Larson core, in red dashed and dotted lines, respectively. The red solid line shows the total angular momentum $L_z$, and the blue lines, the angular momentum of the core+disk region.</figcaption></figure>

which we define as the gas with density $\rho > 10^{-13} \text{ g cm}^{-3}$. The disk+core angular momentum shows the same time evolution as the excess of angular momentum created in the simulation box, with $L_z \approx 0.85 L_{\text{core+disk}}$. This strongly suggests that the problem takes its origin in the central region, and we have calculated that whistler waves do not have time to travel from the center to the outskirt of the box. The same issue arises for different boundary conditions, box sizes, and slope limiters (see Appendix B).

### 5.4.2. Whistler waves on discretized space

We investigated the role of the spatial resolution in this problem. In the previous run, the refinement criterion ensures at least 8 cells per local Jeans length ($\lambda_J/\Delta x > 8$). We ran another simulation with a doubled Jeans length resolution criterion ($\lambda_J/\Delta x > 16$). Figure 17 shows the comparison between the angular momenta in absolute value for the two runs (top panel) as well as their temporal variation (bottom panel). While the most resolved simulation shows momenta larger by 50% (probably because the Hall effect can act at smaller scales during the isothermal stage of the collapse), the increase of angular momentum is lower by roughly 50% in the most resolved case. The local

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

<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Top panel: angular momentum evolution of the two simulations with different refinement criteria (red: $\lambda_J/\Delta x > 8$, blue: $\lambda_J/\Delta x > 16$). The solid curves represent $L_z^+$ and dashed curves represent $L_z^-$. Bottom panel: temporal variation of the positive angular momentum in solid lines, and evolution of the maximum density in dotted lines.</figcaption></figure>

maximum of the angular momentum variation, for the lower resolution case at $t = t_C + 0.400$ kyr, corresponds to the time at which the core expands due to the centrifugal force, as seen in the evolution of maximum density in dotted line. This “bouncing” also happens at a smaller scale for the $\lambda_J/\Delta x > 16$ simulation at $t = t_C + 0.300$ kyr, but is tempered by the lower excess of angular momentum and the additional gravitational force (due to the higher density).

We observed the same behavior when we limit the maximum resolution. In the following set of simulations, we use the 16 points per Jeans length criterion, but we did not allow the code to refine above a given level. Figure 18 represents the angular momentum evolution for maximum levels from 10 ($\Delta x = 16$ au) to 15 ($\Delta x = 0.5$ au). The increase of the angular momentum starts at the same time after the formation of the first Larson core for every simulation. The local maximum is lower and happens sooner as the resolution increases. This trend might be explained by the increase of the angular momentum starting in the most refined cells. With a lower resolution, cells are larger, and therefore hold more mass for the same density. Thus, the gas they harbor takes more time to be set in rotation due to inertia, and once it reaches the critical rotation speed at which the centrifugal force dominates gravity, there is more mass in rotation at a larger effective distance from the center of the simulation (because flow variables are computed at the cell centers), meaning a larger
<figure><img src="image_21.png" /><figcaption>**Fig. 18.** As top panel of Fig. 17, but comparing different maximum resolutions. The vertical bars show when the Truelove et al. (1997) condition is no longer verified for maximum resolution levels of 10 and 11. Please note that the level 15 line in black is the same as the blue curve in Fig. 17.</figcaption></figure>

angular momentum. Nonetheless, a higher resolution seems to temper the increase.

The dissipation of whistler waves on a discretized space could be at stake. As a consequence of the AMR method, short and fast whistler waves generated in high resolution regions could simply be dissipated at refinement level interfaces due to a lack of resolution (the Shannon theorem). To test this hypothesis, we ran the same simulation with a uniform grid of 256³ (level eight of refinement, or a resoluton of 64 au). The aim is purely the numerical study since such a low resolution cannot provide meaningful physical results, and a better uniform resolution would be prohibitively expensive. Figure 19 displays the angular momentum of this simulation compared to the reference case. The increase of angular momentum occurs at a greater intensity, and at an earlier time after the formation of the first core<fnref n="1" />. We thus conclude that the “Shannon dissipation” does not play a significant role in the problem. As for now, we yet have to determine the origin of this issue.

### 5.4.3. Effect of the Hall resistivity magnitude

We repeated the experiment with various Hall resistivities, from $\eta_H = 10^{18} \text{ cm}^2 \text{ s}^{-1}$ to $\eta_H = 10^{21} \text{ cm}^2 \text{ s}^{-1}$, which are relevant values in this context (Marchand et al. 2016). Table 2 shows the formation time of the first Larson core for the four cases, and the ideal MHD case for comparison. Unsurprisingly, a higher Hall resistivity induces a faster rotation speed, which delays the collapse because of the increasing centrifugal force. Below $10^{19} \text{ cm}^2 \text{ s}^{-1}$, the induced rotation is weak and barely contributes to the centrifugal support in the isothermal stage. Figure 20 displays the total angular momentum evolution for all the cases. A higher resistivity leads to an earlier and a steeper increase of the angular momentum, and the ideal case shows perfect conservation.

### 5.4.4. Effect of the Ohmic and ambipolar diffusions

We now aim to include ambipolar and Ohmic diffusions with the Hall term to see whether the other two non-ideal MHD effects,

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** As top panel of Fig. 17 but comparing the reference AMR case (red) and the uniform 256³ grid (blue).</figcaption></figure>

<table><caption>**Table 2.** Larson core formation time for different Hall resistivities.</caption><thead><tr><td>$\eta_H$ (cm$^2$ s$^{-1}$)</td><td>$t_c$ (kyr)</td></tr></thead><tbody><tr><td>0 (ideal case)</td><td>33.550</td></tr><tr><td>$10^{18}$</td><td>33.130</td></tr><tr><td>$10^{19}$</td><td>33.146</td></tr><tr><td>$10^{20}$</td><td>33.287</td></tr><tr><td>$10^{21}$</td><td>34.352</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Evolution of the angular momentum from the first Larson core formation for various Hall resistivities.</figcaption></figure>

can change the angular momentum increase. For the implementation of the ambipolar and Ohmic diffusions in RAMSES, we refer the reader to Masson et al. (2012). Figure 21 shows the evolution of the angular momentum for $\eta_Q = \eta_H = \eta_{AD} = 10^{20} \text{ cm}^2 \text{ s}^{-1}$. The divergence occurs ~400 yr after the first Larson core formation, which is a short 100 yr delay compared with the reference case, and with a sharper initial increase. While the ambipolar and Ohmic diffusions influence the dynamics of the collapse, their respective electric fields, $E_{AD} \parallel (\mathbf{J} \times \mathbf{B}) \times \mathbf{B}$ and $E_Q \parallel \mathbf{J}$, are both perpendicular to the Hall electric field $E_H \parallel \mathbf{J} \times \mathbf{B}$, so they do not counter-act the twisting of the field lines induced by the Hall effect. However, the ambipolar and Ohmic diffusions also dissipate and redistribute the magnetic

<footnote n="1">¹ This is a later absolute time, but the first core formation time is difficult to evaluate at this low resolution.</footnote>
<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Angular momentum evolution in the reference case (red curves) and with the three non-ideal MHD effects (blue curves).</figcaption></figure>

field, effectively weakening the Hall effect in the core, which could be at the origin of the short delay. Nonetheless, our tests show that unless the Hall resistivity is several orders of magnitudes lower than the ambipolar and Ohmic resistivities, the two later effects do not naturally damp the short whistler waves faster than the numerical scheme.

## 6. Conclusions

We have implemented the Hall effect in the AMR code RAMSES to study its impact on star formation and circumstellar disks. The method relies on the modification of the 2D HLL Riemann solver, at the heart of the constrained transport scheme. Our scheme shows a second-order convergence in space and the whistler waves are propagated at the correct frequencies. However, the truncation errors lead to a strong numerical dissipation of short whistler waves, as the fourth power of the wave number, and the AMR method seems inefficient to reduce this dissipation.

During the collapse of a dense core, the Hall effect generates the rotation of the fluid, and thus plays a significant role in the regulation of the angular momentum. While the gas in the mid-plane is set into rotation (in a direction that depends on the sign of $\eta_H$ and the direction of the magnetic field), counter-rotating envelopes develop on each side of the mid-plane by conservation of the angular momentum. This behavior has been observed in previous studies and was therefore expected.

Roughly 300 yr after the first Larson core formation, the angular momentum increases rapidly in the simulation box, meaning that it is no longer conserved. We show that contrarily to the interpretation of Krasnopolsky et al. (2011) of a similar phenomenon, this problem is not related to our boundary conditions. It emerges for every one of our simulations including the Hall effect. Since this divergence begins shortly after the formation of the first Larson core, it may be correlated to the development of the accretion shock, with only few cells describing the sharp transition. However, to date, we have no clear-cut answer concerning the origin of this problem and various leads will be explored in a further work.

Nonetheless, angular momentum conservation is a key issue in numerical simulations exploring the formation of protostellar disks. There is however no mention of this specific problem in previous studies, most of them using different numerical methods. In order to unambiguously assess the validity of disk

formation in such simulations, this issue must at the very least be discussed in detail. A future work will feature additional tests to continue investigating this problem.

**Acknowledgements.** We gratefully thank Geoffroy Lesur, Patrick Hennebelle, Yusuke Tsukamoto, and Kengo Tomida for the interesting ideas and discussions. We also thank the anonymous referee for his thorough report and comments which helped to improve the clarity of the article. We acknowledge financial support from “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, CEA, and CNES, France, and from an International Research Fellowship of the Japan Society for the Promotion of Science. This work was granted access to the HPC resources of CINES (Occigen) under the allocation DARI A0020407247 made by GENCI. Computations were also performed at the Common Computing Facility (CCF) of the LABEX Lyon Institute of Origins (ANR-10-LABX-0066).

## References

* Allen, A., Shu, F. H., & Li, Z.-Y. 2003, ApJ, 599, 351

* ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3

* Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235

* Bodenheimer, P. 1995, ARA&A, 33, 199

* Boss, A. P., & Bodenheimer, P. 1979, ApJ, 234, 289

* Commercon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3

* Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32

* Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659

* Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659

* Falle, S. A. E. G. 2003, MNRAS, 344, 1210

* Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

* Fuente, A., Gerin, M., Pety, J., et al. 2017, A&A, 606, L3

* Galli, D., & Shu, F. H. 1993, ApJ, 417, 243

* Galli, D., Lizano, S., Shu, F. H., & Allen, A. 2006, ApJ, 647, 374

* Gerin, M., Pety, J., Fuente, A., et al. 2015, A&A, 577, L2

* Gerin, M., Pety, J., Commercon, B., et al. 2017, A&A, 606, A35

* Harten, A. 1983, J. Comput. Phys., 49, 357

* Harten, A., Lax, P., & Van Leer B. 1983, SIAM Rev., 25

* Hennebelle, P., & Fromang, S. 2008a, A&A, 477, 9

* Hennebelle, P., & Teyssier, R. 2008b, A&A, 477, 25

* Hopkins, P. F. 2015, MNRAS, 450, 53

* Hopkins, P. F. 2017, MNRAS, 466, 3387

* Johns-Krull, C. M. 2008, in 14th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. G. van Belle, ASP Conf. Ser., 384, 145

* Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128

* Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2010, ApJ, 716, 1541

* Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2011, ApJ, 733, 54

* Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295

* Kunz, M. W., & Mouschovias, T. C. 2010, MNRAS, 408, 322

* Larson, R. B. 1969, MNRAS, 145, 271

* Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

* Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180

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

* Machida, M. N., Inutsuka, S., & Matsumoto, T. 2006, ApJ, 647, L151

* Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18

* Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24

* Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commercon, B. 2016, A&A, 587, A32

* Matsumoto, T., & Tomisaka, K. 2004, ApJ, 616, 266

* Mellon, R. R., & Li, Z. 2009, ApJ, 698, 922

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

* Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315

* Murillo, N. M., Lai, S.-P., Bruderer, S., Harsono, D., & van Dishoeck, E. F. 2013, A&A, 560, A103

* Nakano, T. 1984, Fund. Cosmic Phys., 9, 139

* Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199

* Pezzuto, S., Elia, D., Schisano, E., et al. 2012, A&A, 547, A54

* Price, D. J., & Bate, M. R. 2007, Astrophys. Space Sci., 311, 75

* Price, D. J., Wurster, J., Nixon, C., et al. 2017, PHANTOM: Smoothed Particle Hydrodynamics and Magnetohydrodynamics Code, Astrophysics Source Code Library

* Sano, T., & Stone, J. M. 2002, ApJ, 570, 314

* Shannon, C. 1949, Proc. IRE, 37, 10

* Srinivasan, B., & Shumlak, U. 2011, Phys. Plasmas, 18, 092113

* Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753

* Teyssier, R. 2002, A&A, 385, 337
Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44

Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117

Tóth, G., De Zeeuw, D. L., Gombosi, T. I., & Powell, K. G. 2006, J. Comput. Phys., 217, 722

Tóth, G., Ma, Y., & Gombosi, T. I. 2008, J. Comput. Phys., 227, 6967

Troland, T. H., & Crutcher, R. M. 2008, ApJ, 680, 457

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179

Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, ApJ, 810, L26

Tsukamoto, Y., Okuzumi, S., Iwasaki, K., Machida, M. N., & Inutsuka, S.-i. 2017, PASJ, 69, 95

van Leer B. 1976, in *Computing in Plasma Physics and Astrophysics*, ed. D. Biskamp (Amsterdam: Elsevier)

Wurster, J. 2016, PASA, 33, e041

Wurster, J., Price, D., & Ayliffe, B. 2014, MNRAS, 444, 1104

Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037

Wurster, J., Price, D. J., & Bate, M. R. 2017, MNRAS, 466, 1788

Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 475, 1859

## Appendix A: Analytical solution of the steady-state shock

In this section, we calculate the analytical profile of the shock presented in Sect. 4.2. We first cancel all derivatives except $\frac{\partial}{\partial x}$ in Eqs. (1)–(4). This yields

$$ \frac{d\rho u_x}{dx} = 0, \qquad (A.1) $$

$$ \frac{d}{dx} \begin{pmatrix} u_x B_y - u_y B_x \\ u_x B_z - u_z B_x \end{pmatrix} = \frac{d}{dx} \begin{pmatrix} \left[ \eta_Q + \eta_{AD} \left(1 - \frac{B_z^2}{B^2}\right) \right] \frac{dB_y}{dx} + \left[ \eta_H \frac{B_x}{B} + \eta_{AD} \frac{B_y B_z}{B^2} \right] \frac{dB_x}{dx} \\ \left[ -\eta_H \frac{B_x}{B} + \eta_{AD} \frac{B_y B_z}{B^2} \right] \frac{dB_y}{dx} + \left[ \eta_Q + \eta_{AD} \left(1 - \frac{B_y^2}{B^2}\right) \right] \frac{dB_x}{dx} \end{pmatrix}, \qquad (A.2) $$

$$ \frac{d((\rho u_x)u_x + \rho c_s^2 + B^2/2)}{dx} = 0, \qquad (A.3) $$

$$ \frac{d((\rho u_x)u_y - B_x B_y)}{dx} = 0, \qquad (A.4) $$

$$ \frac{d((\rho u_x)u_z - B_x B_z)}{dx} = 0, \qquad (A.5) $$

$$ \frac{dB_x}{dx} = 0. \qquad (A.6) $$

Equation (A.6) is needed to satisfy the divergence-free condition. Integrating Eq. (A.2) gives two coupled differential equations in $B_y$ and $B_z$. Isolating the derivatives yields

$$ \frac{\partial B_y}{\partial x} = \frac{M_1 R_{22} - M_2 R_{12}}{R_{11} R_{22} - R_{21} R_{12}}, \qquad (A.7) $$

$$ \frac{\partial B_z}{\partial x} = \frac{M_2 R_{11} - M_1 R_{21}}{R_{11} R_{22} - R_{21} R_{12}}, \qquad (A.8) $$

with

$$ M = \begin{pmatrix} u_x B_y - u_y B_x - u_{y,0} B_x + u_{x,0} B_{y,0} \\ u_x B_z - u_z B_x - u_{z,0} B_x + u_{x,0} B_{z,0} \end{pmatrix} \qquad (A.9) $$

and

$$ R = \begin{pmatrix} \eta_Q + \eta_{AD} \left(1 - \frac{B_z^2}{B^2}\right) & \eta_H \frac{B_x}{B} + \eta_{AD} B_y B_z \\ -\eta_H \frac{B_x}{B} + \eta_{AD} B_y B_z & \eta_Q + \eta_{AD} \left(1 - \frac{B_y^2}{B^2}\right) \end{pmatrix}. \qquad (A.10) $$

Eqs. (A.7) and (A.8) are integrated using a Runge–Kutta scheme. We then used the constant quantities $Q = \rho u_x$ (Eq. A.1) and $K = Qu_x + \frac{Q}{u_x}c_s^2 + \frac{B^2}{2}$ (Eq. A.3) to compute the density and velocity profiles.

$$ u_x = \frac{1}{2Q} \left( K - \frac{B^2}{2} - \sqrt{\left(K - \frac{B^2}{2}\right)^2 - 4Q^2 c_s^2} \right), \qquad (A.11) $$

$$ u_y = u_{y,0} + \frac{B_x(B_y - B_{y,0})}{Q}, \qquad (A.12) $$

$$ u_z = u_{z,0} + \frac{B_x(B_z - B_{z,0})}{Q}, \qquad (A.13) $$

$$ \rho = \frac{Q}{u_x}. \qquad (A.14) $$

## Appendix B: Angular momentum issue for different numerical parameters

In this section, we present the evolution of angular momentum for three variations of the numerical parameters compared to the reference case (moncen and periodic). The angular momentum is displayed in Fig. B.1, with the reference case in red. In all cases, the divergence stars 300 yr after the formation of the first Larson core and is roughly the same.

With the minmod slope limiter (blue curves), the angular momentum evolution is roughly the same. In this case, the order of the scheme is first order instead of second, which indicates that the truncation errors do not play a significant role in this problem.

Outflow boundary conditions (green curves) show a larger initial difference between $L_z^+$ and $L_z^-$ because a fraction of the gas leaves the simulation box. The angular momentum divergence still presents the same evolution as in other cases.

The purple curves represent a larger box (with a length of 8 times the dense core radius instead of 4). Both angular momentum components have a higher value due to the larger amount of available mass in the box. The divergence starts again 300 yr after $t_C$ with a comparable increase. These last two cases show that the divergence does not seem related to boundary conditions.

<figure><img src="image_25.png" /><figcaption>Fig. B.1. Positive $L_z^+$ (solid lines) and negative $L_z^-$ (dashed lines) angular momentum evolution after the formation of the first Larson core for several variations of the reference case. Red curves: reference case (generalized moncen slope limiter and periodic boundary conditions), green curves: outflow boundary conditions, blue curves: minmod slope limiter, purple curves: simulation box twice as large. The value of the angular momentum has been divided by two for the latter to improve the readability of the graph.</figcaption></figure>