Issue 
A&A
Volume 619, November 2018



Article Number  A37  
Number of page(s)  16  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201832907  
Published online  30 October 2018 
Impact of the Hall effect in star formation and the issue of angular momentum conservation
^{1}
Department of Earth and Space Science, Osaka University,
Toyonaka,
Osaka
5600043,
Japan
email: marchand@astroosaka.jp
^{2}
CRAL, Ecole normale supérieure de Lyon, UMR CNRS 5574, Université de Lyon,
Lyon,
France
^{3}
School of Physics, University of Exeter,
Exeter,
EX4 4QL,
UK
Received:
26
February
2018
Accepted:
29
July
2018
We present an implementation of the Hall term in the nonideal magnetohydrodynamics equations into the adaptivemeshrefinement 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 twodimensional 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 midplane that depends on the sign of the Hall resistivity, while counterrotating envelopes develop on each side of the midplane. However, we find that the predictability of our numerical results is severely limited. The angular momentum is not conserved in any of our dense corecollapse 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
© 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 overdense 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 gasdust 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 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 fluxfreezing approximation of ideal MHD causes a magnetic pileup 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; JohnsKrull 2008).
Among the possible solutions to solve both the angular momentum and the magnetic flux problem, nonideal MHD is considered to be the most efficient and robust. Nonideal 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 nonideal 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 starformation 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 gridbased 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 adaptivemeshrefinement (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 nonideal MHD
We write the complete set of nonideal MHD equations in conservative form with selfgravity 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 nonideal 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 cgsunits convention, such that the permeability constant μ_{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 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 ydirections, propagating along the zdirection 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 nonconservative form: (9)
We introduced the perturbations in Eqs. (8) and (9), and cancel out the zeroth and secondorder terms:
The system of equations is linear and can be written in matrix form: (14)
A nontrivial 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 xcomponent 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 timestep 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 B_{x} at cell interface . 

Open with DEXTER 
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 secondorder 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 wellknown 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^{n} = u^{n} ×B^{n} is the electric field calculated on cell edges using simple interpolation of u^{n} and B^{n} 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 xdirection 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).
The Hall resistivity η_{H} depends on the density ρ, the temperature T, the cosmicray ionization rate ζ, and the magnetic field amplitude B. These quantities are interpolated on cell edges from the cellcentered 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)). 

Open with DEXTER 
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 zcomponent:
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 anticorrelation 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^{−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) (53)
where S_{L}, S_{R}, S_{B}, and S_{T} are the lowest and highest wave speeds in the x and ydirections. 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 zdirection. The blue area represents the states at the corner for each cell at . 

Open with DEXTER 
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 timestep, otherwise it would be lost since information exchange happens between adjacent cells. At high resolution, we have roughly , and the corresponding timestep (54)
The scaling of the timestep in Δx^{2} induces a fast decrease of the timestep as resolution increases, which is a limitation of current numerical simulation codes. Techniques that allow a faster integration such as Super TimeStepping cannot be employed for the Hall effect because a larger timestep would cut off the highfrequency 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 secondorder 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 timestep. This emancipation results in a faster code because larger timesteps 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 HallEMF 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 nonideal MHD terms were implemented as source terms. Testing them with the ambipolar diffusion shows secondorder 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 shocktest 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 wavenumber k, and study its propagation in the xdirection 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^{−4} to 0.1 cm^{2} s^{−1} 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^{−3} cm^{2} 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 ydirection 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^{−1} and ω ≈ 0.3 s^{−1}. By varying the Hall resistivity η_{H} and the wavenumber 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 u_{y} (in blue) andB_{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 cm^{−1}, η_{H} = 5 × 10^{−3} cm^{2} s^{−1}. 

Open with DEXTER 
Fig. 5 Fourier transform of B_{y}(t) for a random cell in the box (same as Fig. 4). 

Open with DEXTER 
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. 

Open with DEXTER 
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 welldesigned 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^{−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}, η_{H} = 0.005 cm^{2} s^{−1}. 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 secondorder accuracy.
A secondorder 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 secondorder 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^{4}. We evaluated the damping rate for η_{H} = 0.005 cm^{2} s^{−1} at a fixed resolution of 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 wavenumber 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.
Fig. 7 Numerical damping rates as a function of resolution for different slope limiters, and Δ x and Δ x^{2} scalings for comparison. 

Open with DEXTER 
Fig. 8 Numerical damping rates versus wavenumber k. Solid line: numerical results, dashed line: k^{4} scaling. 

Open with DEXTER 
4.1.3 Grid refinement
Here we investigate the behavior of whistler waves on nonregular grids, using the same planeparallel setup of theprevious tests, while refining the grid on specific locations. We once more took k = 5 cm^{−1} and η_{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 ℓ corresponds to a resolution Δ x = 1∕2^{ℓ} cm. For the twolevels tests, the levels of resolution are ℓ = ℓ_{min} + 1 for − 0.25 < x < 0.25, and ℓ = ℓ_{min} everywhere else. For threelevels 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 twolevels simulations and ℓ_{min} = 4, 5, 6 for the threelevels. The whistler wave propagates along the xaxis.
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^{2} scaling as inthe previous test. However, nonregular 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 smallscale 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 leastsquares fit of these points (∝ Δ x^{2.11}). 

Open with DEXTER 
4.2 Shock test
We then tested the Hall effect in a more physical situation. We considered a unidimensional isothermal shock in a Halldominated 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 xmagnetic field are uniform, with cm s^{−1}, B_{x} = 1 G, and we used the resistivities ηΩ = 10^{−9}, η_{H} = −2 × 10^{−2}B, and cm^{2} s^{−1}. We considered a cubic box with a discontinuity in the xdirection, 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 zdirections, and are fixed with the left and right states in the xdirection. 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 steadystate profile of the shock. The analytical calculation of the theoretical profile is detailed in Appendix A.
Figure 10 displays both the numerical steadystate 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.
Initial and boundary flow variables for the shock test.
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 u_{x}, u_{y} and u_{z}, B_{y} and B_{z}. 

Open with DEXTER 
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 “pseudodisk” (Galli & Shu 1993), with a density typically higher than 10^{−15} g cm^{−3}. It is perpendicular to the largescale 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 midplane of the pseudodisk. u_{H} being directly proportional to J, the field lines are then twisted in the plane of the pseudodisk. 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 − η_{H}J_{θ}. 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 corecollapse simulation, and study the role of the Hall effect. We used the standard Boss & Bodenheimer (1979) initial conditions with a nonrotating uniform sphere of gas. The sphere radius is 3712 au and its density ρ_{0} = 4.15 × 10^{−18} g cm^{−3} 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 zdirection, which is equivalent to a masstoflux 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^{20} cm^{2} 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: (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 Δ 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: ρ > 3.8 × 10^{−15} g cm^{−3},
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}∕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 midplane. A large quantity of angular momentum is created in the disk while the whistler waves transport its negative counterpart 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 counterrotation of an hourglassshaped 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 counterrotating 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 zoomout 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 ringshaped 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^{−14} g cm^{−3}. As explained in Sect. 5.1, the region we could define as the pseudodisk 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 counterpart without the Hall effect (i.e., ideal MHD), when both simulations reach a maximum density of ρ = 5 × 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 × 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 nonideal 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. 

Open with DEXTER 
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). 

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

Open with DEXTER 
Fig. 14 Density and velocity map of the reference case at t = t_{C} + 0.800 kyr. 

Open with DEXTER 
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 = t_{C} + 0.220 kyr for the Hall simulation). 

Open with DEXTER 
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 nonconservation 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^{−13} 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} ≈ 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 L_{z}, and the blue lines, the angular momentum of the core+disk region. 

Open with DEXTER 
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^{3} (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^{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 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. 

Open with DEXTER 
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. 

Open with DEXTER 
Fig. 19 As top panel of Fig. 17 but comparing the reference AMR case (red) and the uniform 256^{3} grid (blue). 

Open with DEXTER 
5.4.3 Effect of the Hall resistivity magnitude
We repeated the experiment with various Hall resistivities, from η_{H} = 10^{18} cm^{2} s^{−1} to η_{H} = 10^{21} cm^{2} 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} cm^{2} 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.
Larson core formation time for different Hall resistivities.
Fig. 20 Evolution of the angular momentum from the first Larson core formation for various Hall resistivities. 

Open with DEXTER 
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 nonideal 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^{20} cm^{2} 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} ∥ (J ×B) ×B and E Ω ∥J, are both perpendicular to the Hall electric field E_{H} ∥J ×B, so they do not counteract 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 nonideal MHD effects (blue curves). 

Open with DEXTER 
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 secondorder 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 midplane is set into rotation (in a direction that depends on the sign of η_{H} and the direction of the magnetic field), counterrotating envelopes develop on each side of the midplane 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 clearcut 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 (ANR10LABX0066).
Appendix A Analytical solution of the steadystate 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 divergencefree condition. Integrating Eq. (A.2) gives two coupled differential equations in B_{y} and B_{z}. Isolating the derivatives yields
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. 

Open with DEXTER 
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
 Allen, A., Shu, F. H., & Li, Z.Y. 2003, ApJ, 599, 351 [NASA ADS] [CrossRef] [Google Scholar]
 ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Bodenheimer, P. 1995, ARA&A, 33, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P., & Bodenheimer, P. 1979, ApJ, 234, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [NASA ADS] [CrossRef] [Google Scholar]
 Falle, S. A. E. G. 2003, MNRAS, 344, 1210 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuente, A., Gerin, M., Pety, J., et al. 2017, A&A, 606, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galli, D., & Shu, F. H. 1993, ApJ, 417, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Galli, D., Lizano, S., Shu, F. H., & Allen, A. 2006, ApJ, 647, 374 [NASA ADS] [CrossRef] [Google Scholar]
 Gerin, M., Pety, J., Fuente, A., et al. 2015, A&A, 577, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gerin, M., Pety, J., Commeron, B., et al. 2017, A&A, 606, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harten, A. 1983, J. Comput. Phys., 49, 357 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Harten, A., Lax, P., & Van Leer B. 1983, SIAM Rev., 25 [Google Scholar]
 Hennebelle, P., & Fromang, S. 2008a, A&A, 477, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., & Teyssier, R. 2008b, A&A, 477, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hopkins, P. F. 2015, MNRAS, 450, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Hopkins, P. F. 2017, MNRAS, 466, 3387 [NASA ADS] [CrossRef] [Google Scholar]
 JohnsKrull, C. M. 2008, in 14th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. G. van Belle, ASP Conf. Ser., 384, 145 [NASA ADS] [Google Scholar]
 Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krasnopolsky, R., Li, Z.Y., & Shang, H. 2010, ApJ, 716, 1541 [NASA ADS] [CrossRef] [Google Scholar]
 Krasnopolsky, R., Li, Z.Y., & Shang, H. 2011, ApJ, 733, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. W., & Mouschovias, T. C. 2010, MNRAS, 408, 322 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, Z.Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Londrillo, P., & del Zanna L. 2004, J. Comput. Phys., 195, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Machida, M. N., Inutsuka, S., & Matsumoto, T. 2006, ApJ, 647, L151 [NASA ADS] [CrossRef] [Google Scholar]
 Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masson, J., Teyssier, R., MuletMarquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commercon, B. 2016, A&A, 587, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsumoto, T., & Tomisaka, K. 2004, ApJ, 616, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Mellon, R. R., & Li, Z. 2009, ApJ, 698, 922 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Murillo, N. M., Lai, S.P., Bruderer, S., Harsono, D., & van Dishoeck, E. F. 2013, A&A, 560, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nakano, T. 1984, Fund. Cosmic Phys., 9, 139 [NASA ADS] [Google Scholar]
 Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Pezzuto, S., Elia, D., Schisano, E., et al. 2012, A&A, 547, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Price, D. J., & Bate, M. R. 2007, Astrophys. Space Sci., 311, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., Wurster, J., Nixon, C., et al. 2017, PHANTOM: Smoothed Particle Hydrodynamics and Magnetohydrodynamics Code, Astrophysics Source Code Library [Google Scholar]
 Sano, T., & Stone, J. M. 2002, ApJ, 570, 314 [NASA ADS] [CrossRef] [Google Scholar]
 Shannon, C. 1949, Proc. IRE, 37, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Srinivasan, B.,& Shumlak, U. 2011, Phys. Plasmas, 18, 092113 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Tóth, G., De Zeeuw, D. L., Gombosi, T. I., & Powell, K. G. 2006, J. Comput. Phys., 217, 722 [NASA ADS] [CrossRef] [Google Scholar]
 Tóth, G., Ma, Y., & Gombosi, T. I. 2008, J. Comput. Phys., 227, 6967 [NASA ADS] [CrossRef] [Google Scholar]
 Troland, T. H., & Crutcher, R. M. 2008, ApJ, 680, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179 [NASA ADS] [CrossRef] [Google Scholar]
 Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, ApJ, 810, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Tsukamoto, Y., Okuzumi, S., Iwasaki, K., Machida, M. N., & Inutsuka, S.i. 2017, PASJ, 69, 95 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer B. 1976, in Computing in Plasma Physics and Astrophysics, ed. D. Biskamp (Amsterdam: Elsevier) [Google Scholar]
 Wurster, J. 2016, PASA, 33, e041 [NASA ADS] [CrossRef] [Google Scholar]
 Wurster, J., Price, D., & Ayliffe, B. 2014, MNRAS, 444, 1104 [NASA ADS] [CrossRef] [Google Scholar]
 Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037 [NASA ADS] [CrossRef] [Google Scholar]
 Wurster, J., Price, D. J., & Bate, M. R. 2017, MNRAS, 466, 1788 [NASA ADS] [CrossRef] [Google Scholar]
 Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 475, 1859 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Visualization of the constrained transport (CT) method for the calculation of the magnetic field B_{x} at cell interface . 

Open with DEXTER  
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)). 

Open with DEXTER  
In the text 
Fig. 3 Scheme of a 2D Riemann problem in the zdirection. The blue area represents the states at the corner for each cell at . 

Open with DEXTER  
In the text 
Fig. 4 Typical evolution of u_{y} (in blue) andB_{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 cm^{−1}, η_{H} = 5 × 10^{−3} cm^{2} s^{−1}. 

Open with DEXTER  
In the text 
Fig. 5 Fourier transform of B_{y}(t) for a random cell in the box (same as Fig. 4). 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 7 Numerical damping rates as a function of resolution for different slope limiters, and Δ x and Δ x^{2} scalings for comparison. 

Open with DEXTER  
In the text 
Fig. 8 Numerical damping rates versus wavenumber k. Solid line: numerical results, dashed line: k^{4} scaling. 

Open with DEXTER  
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 leastsquares fit of these points (∝ Δ x^{2.11}). 

Open with DEXTER  
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 u_{x}, u_{y} and u_{z}, B_{y} and B_{z}. 

Open with DEXTER  
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. 

Open with DEXTER  
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 = 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). 

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

Open with DEXTER  
In the text 
Fig. 14 Density and velocity map of the reference case at t = t_{C} + 0.800 kyr. 

Open with DEXTER  
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 = t_{C} + 0.220 kyr for the Hall simulation). 

Open with DEXTER  
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 L_{z}, and the blue lines, the angular momentum of the core+disk region. 

Open with DEXTER  
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. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 19 As top panel of Fig. 17 but comparing the reference AMR case (red) and the uniform 256^{3} grid (blue). 

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

Open with DEXTER  
In the text 
Fig. 21 Angular momentum evolution in the reference case (red curves) and with the three nonideal MHD effects (blue curves). 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.