© ESO 2018

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⁻¹³ g cm⁻³, 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₂ 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²¹ cm² s⁻¹) is three orders of magnitude higher than those of young stars (~10¹⁸ cm² s⁻¹; 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 rotationnally 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 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 RAMSESand 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:

where ρ is the density, u the fluid velocity, P the thermal pressure, B the magnetic field, and Φ the gravitational potential. E_tot and P_tot are the total energy and total pressure:

where ϵ is the specific internal energy and γ the adiabatic index. E_NIMHD represents the electric field created by the three non-ideal MHD effects, and reads (Nakano et al. 2002; Masson et al. 2012) (7)

with ηΩ, η_H, and η_AD the Ohmic, Hall, and ambipolar resistivities, and J = ∇×B the electric current. The MHD equations have been written in the Gaussian cgs-units convention, such that the permeability constant μ₀ = 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 B = Be_z. We introduced small perturbations δu_x, δu_y, δB_x, δB_y of the velocity and magnetic fields in the x- and y-directions, propagating along the z-direction at the frequency ω. Each perturbation is therefore proportional to exp(iωt − ikz). Considering only the Hall effect, the induction Eq. (3) reads (8)

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

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

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

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

with the Alfvén speed. This equation can also be written as

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

Contrarily to the ambipolar and the Ohmic diffusions, the Hall effect is purely dispersive. Any perturbation creates two waves whose frequencies ω are linked to their wavelength λ = 2π∕k. They are called whistler waves. They propagate the magnetic field perturbations along the field lines at the whistler speed: (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 , then the Hall effect is weak and c_w ≈ c_A. In this case, the whistler waves behave as Alfvén waves. Otherwise, the Hall effect plays a significant role in the magnetic field evolution. This defines the Hall length below whichthe 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 , and , 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 (20)

In RAMSES, since Δx = Δy = Δz, we use only Δ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 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.

Fig. 1 Visualization of the constrained transport (CT) method for the calculation of the magnetic field Bx at cell interface .

3.2 Principle of the implementation

The Hall induction equation can be rewritten (21)

where 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 u by the “effective” velocity u + 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 λ = 2Δx, where Δ x is the size of the local cell. Following Eq. (19), the associated whistler speed is (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∕Δ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 u_H also needs to be computed on cell edges. We take the example of the edge located at . 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:

where Eⁿ = uⁿ ×Bⁿ is the electric field calculated on cell edges using simple interpolation of uⁿ and Bⁿ at these locations, for example,

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 yields

where 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,

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 = ∇×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 B^pred the averaged magnetic field, therefore calculated as (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 B^pred. On the edge located at , this yields

Figure 2 depicts the magnetic fields used in the calculation of J on edge . The faces represented are situated in (see left picture).

We thencomputed J as ∇×B

The Hall resistivity η_H depends on the density ρ, the temperature T, the cosmic-ray ionization rate ζ, 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 (39)

Finally, the Hall speed is (40)

with every quantity calculated at .

Fig. 2 Magnetic fields used to calculate J on edge . The right picture represents the cell faces at , as highlighted on the 12 cells block (i − 1⋯i, j − 1⋯j, k − 1⋯k + 1) on the left (the front right middle cell is cell (i, j, 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 , 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:

for any flow variable U, and the magnetic field is averaged on the interfaces: (45) (46) (47) (48)

We write these quantities , , , and . 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 E = (u + u_H) ×B. On the edge, we needed the z-component:

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, λ = u ± c_A.

• Two fast magnetosonic waves, representing the correlation between B and ρ by the coupling between Lorentz force and thermal pressure, (with the component of the Alfvèn speed in the direction of propagation and the sound speed).

• Two slowmagneto sonic waves, representing the anti-correlation between B and ρ by the coupling between Lorentz force and thermal pressure, .

• One contact discontinuity wave λ = 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~Δx⁻¹. 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) (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).

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 .

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 , and the corresponding time-step (54)

The scaling of the time-step in Δx² 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_H 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 unconditionnally stables, 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 (∇×B) ×B is evolved att^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 u × B into a modified constrained transport method. The tests in Sano & Stone (2002) show that the dispersion relation of the Hall effect is verified, butwe 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 performthe 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

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 η_H, typically between 5 and 20 for k and 5 × 10⁻⁴ to 0.1 cm² s⁻¹ for η_H. We used a uniform resolution ofΔ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 η_H = 5 × 10⁻³ cm² s⁻¹, 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 ω ≈ 2 s⁻¹ and ω ≈ 0.3 s⁻¹. By varying the Hall resistivity η_H and the wave-number k, we obtain a set of frequencies. Figure 6 shows the frequencies, normalized by the Hall frequency , 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.

Fig. 4 Typical evolution of uy (in blue) andBy (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 cm−1, ηH = 5 × 10−3 cm2 s−1.

Fig. 5 Fourier transform of By(t) for a random cell in the box (same as Fig. 4).

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.

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 the signal clearly appears in Fig. 4 with the decrease of the wave amplitudes.

The transverse magnetic field can then be estimated as (55)

where f(t) is the periodic contributionand e^−at the damping, with a the damping coefficient in s⁻¹. 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⁻¹, η_H = 0.005 cm² s⁻¹. We also used the slope limitersminmod, moncen, and generalized moncen (β = 1.5), for resolutions of 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∕Δ 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 . It means that at a given number of points per l_H, the damping rate increases as k⁴. We evaluated the damping rate for η_H = 0.005 cm² s⁻¹ at a fixed resolution of cm and k = 5,  7,  8,  10,  20 cm⁻¹. 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⁴, which indicates that smaller wavelengths are substantially more affected by the numerical dissipation than the larger ones.

Fig. 7 Numerical damping rates as a function of resolution for different slope limiters, and Δ x and Δ x2 scalings for comparison.

Fig. 8 Numerical damping rates versus wave-number k. Solid line: numerical results, dashed line: k4 scaling.

4.1.3 Grid refinement

Here we investigate the behavior of whistler waves on non-regular grids, using the same plane-parallel setup of theprevious tests, while refining the grid on specific locations. We once more took k = 5 cm⁻¹ and η_H = 0.005 cm² s⁻¹ 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 ℓ corresponds to a resolution Δ x = 1∕2^ℓ cm. For the two-levels tests, the levels of resolution are ℓ = ℓ_min + 1 for − 0.25 < x < 0.25, and ℓ = ℓ_min everywhere else. For three-levels tests, we impose ℓ = ℓ_min + 2 between − 0.125 < x < 0.125, ℓ = ℓ_min + 1 in − 0.25 < x < 0.25, and ℓ = ℓ_min in the rest of the box. We chose ℓ_min = 5, 6 for the two-levels simulations and ℓ_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 Δ x² scaling as inthe 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.

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 ℓmin = 4, blue for ℓmin = 5, green for ℓmin = 6, and red for ℓmin = 7. The circlepoints represent the regular grids, and the dashed line the least-squares fit of these points (∝ Δ x2.11).

4.2 Shock test

We then tested the Hall effect in a more physical situation. We considered a unidimensional isothermal shock in a 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 cm s⁻¹, B_x = 1 G, and we used the resistivities ηΩ = 10⁻⁹, η_H = −2 × 10⁻²B, and cm² s⁻¹. 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 ∈ [0;0.3] while the region x ∈ [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 hasa 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 ℓ_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.

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: ρ and ux, uy and uz, By and Bz.

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 (56)

This expression shows that Hall effect generates a toroidal magnetic field B_θ, 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⁻¹⁵ g cm⁻³. 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. u_H being directly proportional to 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 − η_HJ_θ. 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 −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 ρ₀ = 4.15 × 10⁻¹⁸ g cm⁻³ for a total mass of 1.5 M_⊙. The thermal over gravitational energy ratio is α = 0.25 with a uniform temperatureT = 10 K. The external medium is 100 times less dense. The magnetic field is uniform at B = 90.3 μG in the whole box along the z-direction, which is equivalent to a mass-to-flux ratio of μ = 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 η_H = 10²⁰ cm² s⁻¹. 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: (57)

The initial grid is refined at AMR level 5 (32³) 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 Δ 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⁻¹³ g cm⁻³. Finally, we refer to a structure rotationally supported against gravity as “disk”. We used the criteria of Joos et al. (2012):

• density criterion: ρ > 3.8 × 10⁻¹⁵ g cm⁻³,

• rotationally supported: u_θ > 2u_r, with u_θ and u_r the azimuthal and radial velocities, respectively,

• hydrostatic equilibrium within the disk: u_θ > 2u_z,

• not within the first core: ρu²∕2 > 2P_therm, with P_therm the thermal pressure.

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_θ ∕B and of the azimutal velocity u_θ, 0.230 kyr after the formation of the first core. The twisting of the field 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_θ of the fluid to the azimuthal Hall speed and to the Keplerian velocity. The Hall speed is defined accordingly in Sect. 3.2: (58)

Using the second Newton law, the Keplerian velocity was calculated as (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 ≈ 160 au, down to a density of ρ ≈ 10⁻¹⁴ g cm⁻³. 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 . 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, is lower (in negative value) while 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 ρ = 5 × 10⁻¹² g cm⁻³. The magnetic field in ideal MHD is roughly 0.5 G at the center, while it is bounded to 3 × 10⁻² G with the Hall effect. On the contrary, the magnetic field amplitudes are very similar at densities below 10⁻¹⁵–10⁻¹⁴ g cm⁻³, 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.

Fig. 11 Toroidal magnetic field Bθ (top panel) and azimuthal velocity uθ (bottom panel) on the plane x = 0, 0.230 kyr after the formation of the first core.

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 = tc + 0.500 and t =tc + 0.800 kyr (top and bottom panels). The black dotted line represents the density profile (right axis).

Fig. 13 As Fig. 12, at larger radius and without the Keplerian velocity.

Fig. 14 Density and velocity map of the reference case at t = tC + 0.800 kyr.

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 × 10−12 g cm−3 (at t = tC + 0.220 kyr for the Hall simulation).

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 L the angular momentum vector (60)

where m_i, r_i, and u_i represent the mass, the position vector, and the velocity vector of cell i. The systembeing nearly axisymmetric, the main component of this vector is L_z (i.e., the initial magnetic field direction), that can also be calculated as (61)

with r_c the cylindrical radius and u_θ the azimuthal velocity. We define the “positive” and “negative” momenta and as follows:

We therefore have .

In our setup, there is initially no angular momentum in the system. thus represents the “positive” rotation created by the Hall effect, balanced by the “negative” rotation , and the numerical solution should verify . 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 and ).

In Fig. 16, we have plotted the evolution of , , 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 . The angular momentum is therefore not conserved in the simulation box.

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, which we define as the gas with density ρ > 10⁻¹³ g cm⁻³. The disk+core angular momentum shows the same time evolution as the excess of angular momentum created in the simulation box, with L_z ≈ 0.85L_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).

Fig. 16 Positive and negative 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 Lz, and the blue lines, the angular momentum of the core+disk region.

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(λ_J∕Δx > 8). We ran another simulation with a doubled Jeans length resolution criterion (λ_J ∕Δ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 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 λ_J ∕Δ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 (Δ x= 16 au) to 15 (Δ 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 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 whistlerwaves 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¹. 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 ofthis issue.

Fig. 17 Top panel: angular momentum evolution of the two simulations with different refinement criteria (red: λJ ∕Δ x > 8, blue: λJ ∕Δx > 16). The solid curves represent and dashed curves represent . Bottom panel: temporal variation of the positive angular momentum in solid lines, and evolution of the maximum density in dotted lines.

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.

Fig. 19 As top panel of Fig. 17 but comparing the reference AMR case (red) and the uniform 2563 grid (blue).

5.4.3 Effect of the Hall resistivity magnitude

We repeated the experiment with various Hall resistivities, from η_H = 10¹⁸ cm² s⁻¹ to η_H = 10²¹ cm² s⁻¹, 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¹⁹ cm² s⁻¹, 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.

Fig. 20 Evolution of the angular momentum from the first Larson core formation for various Hall resistivities.

5.4.4 Effect of the Ohmic and ambipolar diffusions

We now aimto include ambipolar and Ohmic diffusions with the Hall term to see whether the other two non-ideal MHD effects, 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 ηΩ = η_H = η_AD = 10²⁰ cm² s⁻¹. 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 ∥ (J ×B) ×B and E Ω ∥J, are both perpendicular to the Hall electric field E_H ∥J ×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 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.

Fig. 21 Angular momentum evolution in the reference case (red curves) and with the three non-ideal MHD effects (blue curves).

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 η_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. Thereis however no mention of this specific problem in previous studies, most of them using different numerical methods. Inorder 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 PhysiqueStellaire” (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).

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 in Eqs. (1)–(4). This yields

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

with (A.9)

and (A.10)

Eqs. (A.7) and (A.8) are integrated using a Runge–Kutta scheme. We then used the constant quantities Q = ρu_x (Eq. A.1) and (Eq. A.3) to compute the density and velocity profiles.

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.

Fig. B.1 Positive (solid lines) and negative (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.

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 and 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.

References

All Tables

All Figures

	Fig. 1 Visualization of the constrained transport (CT) method for the calculation of the magnetic field Bx at cell interface .
In the text

	Fig. 2 Magnetic fields used to calculate J on edge . The right picture represents the cell faces at , as highlighted on the 12 cells block (i − 1⋯i, j − 1⋯j, k − 1⋯k + 1) on the left (the front right middle cell is cell (i, j, k)).
In the text

	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 .
In the text

	Fig. 4 Typical evolution of uy (in blue) andBy (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 cm−1, ηH = 5 × 10−3 cm2 s−1.
In the text

	Fig. 5 Fourier transform of By(t) for a random cell in the box (same as Fig. 4).
In the text

	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.
In the text

	Fig. 7 Numerical damping rates as a function of resolution for different slope limiters, and Δ x and Δ x2 scalings for comparison.
In the text

	Fig. 8 Numerical damping rates versus wave-number k. Solid line: numerical results, dashed line: k4 scaling.
In the text

	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 ℓmin = 4, blue for ℓmin = 5, green for ℓmin = 6, and red for ℓmin = 7. The circlepoints represent the regular grids, and the dashed line the least-squares fit of these points (∝ Δ x2.11).
In the text

	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: ρ and ux, uy and uz, By and Bz.
In the text

	Fig. 11 Toroidal magnetic field Bθ (top panel) and azimuthal velocity uθ (bottom panel) on the plane x = 0, 0.230 kyr after the formation of the first core.
In the text

	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 = tc + 0.500 and t =tc + 0.800 kyr (top and bottom panels). The black dotted line represents the density profile (right axis).
In the text

	Fig. 13 As Fig. 12, at larger radius and without the Keplerian velocity.
In the text

	Fig. 14 Density and velocity map of the reference case at t = tC + 0.800 kyr.
In the text

	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 × 10−12 g cm−3 (at t = tC + 0.220 kyr for the Hall simulation).
In the text

	Fig. 16 Positive and negative 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 Lz, and the blue lines, the angular momentum of the core+disk region.
In the text

	Fig. 17 Top panel: angular momentum evolution of the two simulations with different refinement criteria (red: λJ ∕Δ x > 8, blue: λJ ∕Δx > 16). The solid curves represent and dashed curves represent . Bottom panel: temporal variation of the positive angular momentum in solid lines, and evolution of the maximum density in dotted lines.
In the text

	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.
In the text

	Fig. 19 As top panel of Fig. 17 but comparing the reference AMR case (red) and the uniform 2563 grid (blue).
In the text

	Fig. 20 Evolution of the angular momentum from the first Larson core formation for various Hall resistivities.
In the text

	Fig. 21 Angular momentum evolution in the reference case (red curves) and with the three non-ideal MHD effects (blue curves).
In the text

Fig. B.1 Positive (solid lines) and negative (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.

In the text