EDP Sciences
Free Access
Issue
A&A
Volume 587, March 2016
Article Number A32
Number of page(s) 20
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201526371
Published online 12 February 2016

© ESO, 2016

1. Introduction

In the past few decades, star formation studies have been struggling to properly describe the mechanism of angular momentum transport, regulated by angular momentum conservation and magnetic braking. The formation of rotationally supported structures, that is, the protoplanetary disks that are expected to give birth to planets and/or binaries, is directly related to the amount of angular momentum in protostellar systems. While there is plenty of observational evidence for large Keplerian disks1 around Class-II and Class-I objects, their presence around younger Class-0 objects is still subject to debate (as discussed in Maury et al. 2010 or Tobin et al. 2012, 2013; see the review by Li et al. 2014a).

The simplified framework of the ideal magnetohydrodynamics (MHD) used in the first studies led to the disappearance of the large Keplerian disks that are easily formed in hydrodynamical simulations as a result of the very effective magnetic braking created by a strong pile-up of magnetic flux towards the centre of the collapsing system (Galli et al. 2006; Allen et al. 2003, Price & Bate 2007; Hennebelle & Teyssier 2008; Matsumoto & Tomisaka 2004; Hennebelle & Fromang 2008; Commerçon et al. 2010). In these simulations, disks were found to form only for unrealistically weak magnetic field intensities (corresponding to a mass-to-flux ratio more than 10 times the critical value derived by Mouschovias & Spitzer 1976). Other consequences of the magnetic flux freezing assumption inherent to the ideal MHD approximation include the strong resulting magnetisation of protostars compared to the low observed values in stars (Crutcher 2012), the generation of violent interchange instabilities at the protostar-disk interface (Li et al. 2014b), the growth of the pseudo-disk (Hennebelle & Fromang 2008), and the distortion of the upper layer of the pseudo-disk due to reconnection and the split monopole (Galli & Shu 1993; Li et al. 2014b).

In the framework of ideal MHD, there is no possibility to regulate the magnetic flux pile-up and its consequences, except for the intrinsic numerical resistivity due to both the numerical method used to solve the induction equation and the numerical grid resolution. To circumvent this problem, the magnetic field redistribution needs to be correctly addressed in the framework of complete non-ideal MHD. Following the pioneering work of Mestel & Spitzer (1956), the study of non-ideal MHD effects has been the focus of intense research in the recent years, as illustrated by the studies of Duffin & Pudritz (2008) and Mellon & Li (2009) for ambipolar diffusion, or Machida et al. (2006) and Machida & Matsumoto (2011) for a generic resistivity. While it is still unclear if non-ideal MHD can, by itself, solve the problems regarding the formation of disks (Krasnopolsky et al. 2010), a physical dissipative scale for the magnetic flux definitely improves the regulation of magnetic flux pile-up in star formation studies.

In this work, we study the effects of ambipolar diffusion in the context of star formation, more specifically, of low-mass star formation, and highlight the differences compared to ideal MHD simulations. The influence of turbulent initial conditions will be studied in a forthcoming paper (Paper II). The article is organised as follows. In Sect. 2 we discuss the framework and numerical setup of the study. In Sect. 3 we focus on the general description and properties of the collapsing core until formation and early evolution of the first Larson core, for various cases. In Sect. 5 we discuss the long-term evolution of the structures and highlight the limits of the ideal MHD framework that are due to numerical issues. In Sect. 4, we examine the formation of rotationally supported structures in non-ideal MHD. Last, we summarise our findings in Sect. 7.

2. General context

2.1. Physical framework

In star formation, the ionisation fraction is low, and quasi-equilibrium holds between the Lorentz force and the plasma-neutrals friction force as a result of collisions. Therefore, we can drop the pressure and gravitational forces for charged particles when writing the equations of motion for each fluid particle. In the case of positively charged species of number density ni, velocity vi, atomic number Zi , and collision rate νij with the species j, and using the subscript e for the negatively charged species, the equations of motion read: (1)while the generalized Ohm’s law is written (see e.g. Balbus & Terquem 2001) (2)with The symbols vn, ρ, and ρi denote the velocity of neutral particles, the total density, and the ion density, respectively. In this framework, the drag force per unit volume exerted by the neutrals on the ions compensates for the Lorentz force and reads Fdrag = γADρiρ(vivn). When considering multiple species (see Kunz & Mouschovias 2009; Nakano et al. 2002), Ohm’s law can be written as (5)with |||| standing for the L2 norm, and the Ohmic, Hall, and ambipolar diffusivities are defined as The parallel, perpendicular, and Hall conductivity components (σ, σ, and σH, respectively) compose the conductivity tensor (9)from Faraday’s law j = σ(E + v × B). These resistivities reduce to the values displayed in Eq. (2) in a three-fluid description. In this case, The present work is essentially devoted to the first collapse and the formation of the first core. At this stage and on this scale, only ambipolar diffusion plays a role. In the following, we therefore only focus on this term. Rewriting the above equations with ηΩ = ηH = 0 and (13)Eq. (5) reduces to (14)In contrast to a Laplace operator, as in Ohmic dissipation, Eq. (14) does not allow magnetic reconnection, since it corresponds to another flux-freezing condition at a different speed . We stress that this is true only under the one-fluid approximation, as pointed out by Tsap et al. (2012).

2.2. Magnetic braking

The main effect of magnetic braking in the context of prestellar core collapse is to slow down the rotationally supported structures that develop during the gravitational collapse of matter with non-zero net angular momentum; this occurs because of the conservation of angular momentum. One of the consequences is that the formation of large and rapidly rotating disks is hampered, as found in purely hydrodynamical simulations. The braking arises essentially from the magnetic tension that prevents the distortion of the field lines, which are coupled to the flow through the ionised species. The angular momentum is transported along the field lines at a speed close to the Alfvén speed, depending on the topology of the magnetic field. The braking efficiency can be characterised by deriving an Alembert equation of propagation for the angular momentum (e.g. Gillis et al. 1974, 1979; Mouschovias 1978, 1979; Mouschovias & Paleologou 1980), but can be estimated by order-of-magnitude calculations (Joos et al. 2012). In cylindrical coordinates, with the z-axis aligned with the mean angular momentum direction vector, the angular momentum flux due to the magnetic field reads (Joos et al. 2012) (15)It is proportional to the radial and poloidal components of the field Br and Bz, while scaling with the square of the toroidal component Bθ, making this latter the main focus of magnetic braking studies.

2.3. Numerical setup

2.3.1. Methods

To carry out our study of collapsing magnetised molecular cloud cores, we used the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002; Fromang et al. 2006) in its non-ideal MHD extension (Masson et al. 2012). RAMSES solves the complete set of MHD equations (self-gravity, Euler’s fluid flow equations, and the induction equation with non-ideal terms) using the constrained transport method, which preserves the divergence-free condition for the magnetic field to machine precision. The adaptive mesh is extremely well suited to protostellar collapse calculations, where many levels of refinement are needed to efficiently describe spatial scales spanning 104−105 orders of magnitude in a single simulation. AMR is also a powerful tool for fragmentation studies and disk formation in a turbulent medium where nested grids are difficult to use.

The grid refinement criterion is based on the Jeans mass, ensuring the Jeans length is always sampled by at least eight cells. The coarse grid has a resolution of 323, and 11 levels of AMR were used, resulting in a maximum resolution of 0.15 au at the finest level. Our general set of equations includes the conservation of mass, the mean neutral gas dynamics (Euler equation), the induction equation with the ambipolar diffusion electromotive force, self-gravity through the Poisson equation, and the divergence-free constraint: Here ρ is the mean fluid density, v the mean fluid velocity, P the thermal pressure of the gas, and is the electromagnetic force due to the ambipolar diffusion, as derived in Eq. (5). The energy equation is approximated by a barotropic equation of state (see below).

Magnetic resistivities are calculated using a reduced chemical network including neutral and charged species, as well as dust grains. We followed Kunz & Mouschovias (2009) to compute the relevant charged species abundances including grains sizes in a classical MRN distribution (Mathis et al. 1977), which were sampled using 50 bins. For an exhaustive description of the chemical model used and its application in the context of star formation, we refer to Marchand et al. (2015). We computed a three-dimensional table of density, temperature, and magnetic field dependent resistivities covering the ranges 10-24 < ρ < 10-10 g cm-3, 5 < T < 2000 K, and 10-6 < B < 102 G, respectively. During the simulations, the resistivities in each grid cell are interpolated on-the-fly according to the local state variables, which greatly reduces computational cost but implies thermodynamical equilibrium.

Table 1

Summary of the initial conditions for the different simulations.

Table 2

Summary of the main properties of the first core for fiducial cases with solid-body rotation.

2.3.2. Initial conditions

We adopted initial conditions similar to those in Commerçon et al. (2010), who followed Boss & Bodenheimer (1979). A magnetised uniform-density sphere of molecular gas, rotating about the z-axis with solid body rotation, is placed in a surrounding medium a hundred times less dense with equal pressure. The prestellar core mass has a mass of 1 M, a radius R0 = 2500 au2 and a ratio of rotational over gravitational energy of βrot = 0.02. The magnetic field is initially parallel to, and invariant along, the rotation z-axis. The field strength is stronger in a cylinder of radius R0 (with the dense core at its centre) than in the surrounding medium, with Bz(r>R0) = Bz(R0) ∗ 1002/3, where the factor of 100 comes from the difference in density between the core and the surroundings3. We define a mass-to-flux ratio parameter similar to the one defined by Mouschovias & Spitzer (1976) to measure the importance of the magnetisation in the core: (20)with the critical value (Mouschovias & Spitzer 1976). We note that μ(r = R0) is strictly equal to the theoretical value for a homogeneous cloud permeated by vertical field lines.

Even though a Bonnor-Ebert (BE) density profile may better fit observations of dense cores (see Andre et al. 2000; Belloche et al. 2002) and has an analytical foundation (Li & Shu 1996; Hunter 1977), we assumed a magnetised medium permeated by straight parallel field lines, and it is hardly possible to end up with a BE density profile without bending the field lines during the formation of the density enhancement (the dense core). Some authors (e.g., van Loo et al. 2008) found non-linear density enhancements in simulated turbulent molecular sheets via slow-mode magnetic waves, which left the magnetic field unchanged in the cores, but it remains unclear whether this process does occur in molecular clouds. Our own tests tend to show that the initial density profile is not a critical matter because the infalling material adopts a BE-like density profile very shortly after the beginning of the dense core collapse.

As the aim of the present paper is to focus on the effect of non-ideal MHD effects, notably ambipolar diffusion, on prestellar core collapse and disk formation, we used a barotropic equation of state (EOS) to mimic the effect of radiative transfer instead of solving the full set of radiation magnetohydrodynamics equations. This was done for the sake of simplicity and to save computational cost. The gas pressure is thus related to the density as (21)where cs is the gas sound speed and nH the number density of atoms. After an initial phase of isothermal collapse up to a density nH ≃ 10-13 g cm-3, the first hydrostatic Larson core (Larson 1969) forms when the gas becomes optically thick enough to stop radiative cooling, which was until then counter-balancing the compressional heating. This adiabatic phase is reproduced with the barotropic EOS by using a polytropic index n for a gas (which is equal to the ratio of specific heats γ) , higher than the critical value for stability against gravitational collapse . The temperature and density of the gas begin to rise inside the core as it accretes material from the surrounding envelope. The energy sink provided by the dissociation of H2 molecules, which occurs around T ~ 2000 K, triggers a second phase of collapse (Larson 1969). We here focus on the properties of the first Larson core. Our general approach allows us to accurately describe the first Larson core and its surroundings and the magnetic flux transport without needing to introduce a sink particle that would limit the resolution. It is beyond the scope of the present paper to include long-term evolution, which leads to the formation of the second core. This would imply the accurate treatment of jets, high-energy radiation, and remaining non-ideal MHD terms and will be addressed in a forthcoming study.

We have performed eight simulations for which we varied the mass-to-flux ratio μ = 2 and 5, the angle between the rotation axis and the magnetic field initial direction (0 and 40 degree) and used ideal (iMHD) or non-ideal (niMHD) magnetohydrodynamics (accounting only for ambipolar diffusion). We also performed additional simulations with increased and decreased resolutions to study the convergence of our results for our fiducial case (, aligned case). The various run parameters are given in Table 1.

3. Early evolution

3.1. Aligned case

In this section, we study the general properties of the collapsing system by comparing the iMHD case and a case with ambipolar diffusion in the fiducial aligned case. We first focus on the case, which we describe in detail. We then report the differences or similarities with the case.

3.1.1. First Larson core

We determined the onset of the first core formation as the point when the density in the computational domain reaches 10-12 g cm-3. We subsequently defined the first core itself by the cells fulfilling the criterion ρ> 10-12 g cm-3. The core radius, mass, mass-to-flux ratio, and peak magnetic field value were measured 200 years after its formation. The properties of the first core are summarised in Table 2.

The first core is oblate with a radius ranging from 8 to 9 au, independently of the initial magnetisation ( or ). While pure hydrodynamical simulations (e.g. Bate 2011; Tomida et al. 2010a) and some radiation iMHD (RMHD) calculations (Tomida et al. 2010b) found a similar result with a flattened first core, other RMHD studies (Commerçon et al. 2010) did not find flattened cores in the case μ = 5. As a result of the additional magnetic support, the free-fall (thus core formation) timescale is twice longer for than for . The additional support also hinders accretion, and the first core is significantly less massive for than for . Similarly, for the stronger magnetisation (), the core mass is lower for niMHD than for iMHD. Indeed, the strongest magnetic field inside the core remains of the order of 10-1 G in niMHD, while it is ten times stronger in iMHD. This effect of ambipolar diffusion on the magnetisation during the core formation is illustrated by the differences in the mass-to-flux ratio μ at 10 and 100 au between the ideal and non-ideal MHD cases. While μ is higher than 10 at 10 au and at ≲ 4 at 100 au in niMHD, it is at most ≲ 5 in iMHD.

3.1.2. Magnetic field repartition

Ambipolar diffusion enables neutral particles to overcome magnetic field lines and redistribute the magnetic field. We first focus on the magnetic field intensity and topology to understand its consequences for the dynamics of the collapsing core. Figure 1 shows the magnetic field repartition as a function of gas density for the iMHD and niMHD runs. The flux-freezing behaviour, Bρ1 /ξ with ξ ~ 3/2, is obvious for iMHD (red) with a highest magnetic field value of ~1 G. The niMHD simulation (blue) shows an identical initial evolution for densities below ~ 10-14 g cm-3. After this point, ambipolar diffusion starts to dominate the later evolution of the core, and a very well defined diffusion plateau forms that is usually referred to as the decoupling stage in star formation (Desch & Mouschovias 2001). In both the iMHD and niMHD calculations, the generation of outflows can be identified by the high values of the magnetic field at high density, compared with the expected perfect flux-freezing result.

As seen in Fig. 1, a scaling relation Bρ2/3 seems to represent the evolution of the magnetic field during the collapse better than the traditional relation (see the differences in the scaling for the different components of the field in Hennebelle & Fromang 2008). Further details about this scaling relation and the analytical derivation of an estimate of the field saturation value are given in Appendix A.

Figure 2 illustrates the comparison between μ = 2 and μ = 5 niMHD. Both simulations show similar initial evolutions (Fig. 2a) when flux-freezing is still dominant (μ = 2 naturally shows a higher magnetic field intensity at a given density), but have slightly different saturation values. At later stages (Fig. 2b) a plateau forms in both low magnetisation and high magnetisation runs with Bplateau ≃ 0.1 G.

thumbnail Fig. 1

Magnetic field distribution as a function of gas density for aligned, using the iMHD (red) and the niMHD (blue) formalisms, 200 years after first core formation. The solid black line shows the scaling of the magnetic field Bρ2/3 and the dashed black line the scaling Bρ1/2.

Open with DEXTER

thumbnail Fig. 2

a) Early time, when flux freezing is still relevant. b) After the decoupling stage. B(ρ) for (red) and (blue).

Open with DEXTER

This plateau arises when ambipolar diffusion becomes efficient enough to overcome the dynamical effects. The non-linearity of the effective diffusion coefficient, , explains the self-regulation and the formation of the diffusion plateau. Indeed, for any non-potential ( × B ≠ 0) field configuration, an increase in the magnetic field leads to a rise of the local diffusion coefficient and thus in turn to a decrease in field intensity. The plateau first forms when the dynamical evolution of the magnetic field, v × B, no longer dominates the diffusion term, (with |||| standing for the L2 norm in this paper). Assuming flux-freezing holds during the first isothermal phase of the collapse (i.e. Bρ1 /ξ), we can estimate the saturation value for the magnetic field (22)where γAD is the drag coefficient (Eq. (3)), C the ionisation fraction such that , G the gravitational constant, B0 and ρ0 the typical initial conditions for the magnetic field and gas density, and ξ the power index for the proportionality law B(ρ) ∝ ρ1 /ξ. The initial thermal support and mass-to-flux ratio are linked to the values of B0 and ρ0. Further details on the derivation of Bsat can be found in Appendix A. The values we obtain for Bsat with for various values of thermal support α, that is, the ratio of thermal over gravitational energy, are shown Fig. 3 (thick black solid line). For , which corresponds to a homogeneous sphere permeated by a uniform magnetic field, we find that the saturation value does not depend on α. In this case, the saturation value estimate for is Bsat = 0.13 G, very close to the numerical value B ≲ 0.1 G (see Fig. 1). According to this simple order-of-magnitude calculation, the saturation value should even be lower in more magnetised cores, with μ< 5. For example, the case yields Bsat = 8 × 10-3 G. In our numerical simulation, we do not observe such a strong diminution, but the saturation value is lower for than for (see Fig. 2). The discrepancy between the analytical prediction and the numerical value for Bsat in the strongly magnetised case originates from the fact that in the latter case the magnetic field distribution departs appreciably from a simple power-law parametrisation Bρ2/3 (see Fig. 2). Using a slightly different power index that better fits the effective B(ρ) repartition in the case, for instance ξ = 1.73, as shown on the inset in Fig. 3, we obtain Bsat = 0.06 (coloured lines in Fig. 3), which is closer to the value obtained in the simulation.

thumbnail Fig. 3

Saturation value as a function of μ for various values of thermal support α. The points corresponding to μ = 2 (lower values) and μ = 5 (upper values) are marked in the figures.

Open with DEXTER

Figure 4 shows snapshots of the mass-to-flux ratio , as defined Eq. (20), as a function of radius at different times in the simulations. At large radii (r ≳ 1000 au), the magnetisation is similar for iMHD (red) and niMHD (black), confirming the fact that at these scales the AD timescale is orders of magnitude longer than the dynamical timescale. There is a kink in the AD case at a few tens of au (highlighted by the green square in the figure) that slowly propagates outwards (compare the blue dashed and solid curves); this corresponds to the efficient diffusion of the field at these densities.

thumbnail Fig. 4

, aligned case. Mass-to-flux ratio () after first core formation (at ~ 24 kyr) for the niMHD (black) and IMHD (red) cases (solid lines). We also display the initial condition (dotted blue line) and two later outputs for the niMHD case (dashed and dotted solid black lines). The decoupling stage region due to magnetic field pile-up is emphasised with a green square.

Open with DEXTER

3.1.3. Outflows

The piling-up of the toroidal component of the magnetic field during the collapse of the dense core ultimately leads to the formation of a growing vertical structure, called magnetic tower. We note that we use the generic term magnetic tower to describe any outflowing structure. In reality, there are several ways to launch outflows in a magnetised environment, as studied initially by Lynden-Bell & Boily (1994) and in contemporary studies by Hennebelle & Fromang (2008), Ciardi & Hennebelle (2010). This outflow is launched shortly (≲ 1 kyr) after the formation of the first core. Figure 5 (top row) displays the azimuthally averaged density and velocity fields (panels a and c) and Alfvén speed with the magnetic field direction (panels b and d) for the iMHD and niMHD μ = 5 runs. The simulations are compared 500 years after the formation of the first Larson core, but the same conclusions hold during all the time evolution. In all panels, the z-axis is aligned with the direction of the average angular momentum vector computed in a sphere of radius 100 au, centred around the densest cell in the simulation. The velocity vectors are overlaid on the density maps, while the direction of the magnetic field in the (r,z) plane is superimposed on the Alfvén speed colour maps. The colour coding for the magnetic field illustrates the intensity of the toroidal component Bφ.

The magnetically driven outflow is reinforced in the iMHD case by the more important magnetic field pile-up. This accumulation stems both from the radial bending of the field lines, as clearly seen in Fig. 5b, and from the increasing toroidal field (of prime importance for magnetic braking, as shown by Eq. (15)) that is generated by the rotation, as clearly seen in Fig. 5b and d, especially close to the first Larson core where the toroidal support is significantly stronger in (b). While the gas inside the magnetic tower is slightly less dense in the AD case, the lower magnetic field intensity (compared to the iMHD case) yields an overall Alfvén velocity that is lower by one order of magnitude. This in turn explains the weaker magnetic braking, since the angular momentum is carried away by slower Alfvén waves4.

thumbnail Fig. 5

Top row: , aligned case, snapshots at a time t = 1.02 tfree - fall (500 years after first core formation). Panel a) shows a map of the density in the ideal MHD simulation. The black solid lines are logarithmically spaced density contours. The arrows represent the velocity field where blue to green corresponds to the increasing magnitude of the velocity in the (r,z) plane. Panel b) is a map of the Alfvén speed (blue), with the magnetic field direction in the (r,z) plane indicated by the segments. Yellow to red in these segments corresponds to increasing relative magnitude of the toroidal component of the field compared to the total field. Panels c) and d) are the same as a) and b), respectively, but for the niMHD case. Bottom row: same as the top row for μ = 2, and at the same time t = 1.01 tfree - fall (500 years after first core formation). Every quantity is azimuthally averaged.

Open with DEXTER

thumbnail Fig. 6

Four main panels show the aligned (top left, a), misaligned (top right, b), aligned (bottom left, c), and misaligned (bottom right, d) cases 500 years after first core formation. In each main panel, subpanel a) shows a colour map of the azimuthal average of . Green to blue means that dynamical processes (collapse, Keplerian disk) dominate, while red to yellow means that ambipolar diffusion dominates. The grey segments show the direction of the magnetic field, and black solid lines are isodensity contours. Panels b) and c) display the magnetic field magnitude B and the value of Am as a function of density. Red cells in the Am scatter plot are cells where B(ρ) > 0.08 G, while grey cells have B(ρ) < 0.08 G.

Open with DEXTER

In the iMHD simulation, the velocity field in the tower (in the rz plane) is vertical and the gas mainly flows out, the launching mechanism taking root in the strong toroidal field and differential rotation close to the Larson core. The field is dominantly vertical, and neutral matter follows the field lines. When AD is included, the growth of the magnetic tower is still magnetically regulated, but the growth is slower and the velocity field in the tower is almost null, or slightly directed towards the protostar. A detailed analysis shows that there is no real outflowing gas motion of gas per se, but that the magnetic tower structure or growth is supported by magnetic pressure. We also note that close to the core (r< 50 au), while the field lines are pinched in the iMHD case (split monopole topology), field re-distribution has operated in the niMHD run because of the ambipolar diffusion (the magnetic field vectors are much more vertical).

μ = 2 :

in a more magnetised case, the picture is very different. The bottom row of Fig. 5 shows the structure of the outflow for μ = 2. Magnetic braking is much more effective (the magnetic field is overall stronger, yielding a stronger magnetic braking, as seen Sect. 2.2), causing the field topology to come closer to a split monopole configuration at large scale, with strong pinching of the lines in the equatorial plane. The enhanced braking weakens the outflow-launching mechanism because the pile-up of the toroidal field is less effective. For iMHD, there is still an outflowing feature, but its velocity field is either null in the rz plane or even falls back onto the core. The boundary of the outflow region is less well defined, but is characterised by a discontinuity in the velocity field. In the AD case, the braking is stronger for μ = 2 than for μ = 5 at the early stages of the collapse (see Sect. 3.1.4), which results in a weaker magnetic field (the diffusion plateau is lower) and reduces the pile-up of toroidal field close to the core. This weakens the launching process to the point that no outflow is produced in this case. The Alfvén speed close to the core is very similar to the weak field μ = 5 run (Fig. 5d), as expected from the field magnitude distribution as a function of density (Fig. 2), where B(ρ> 10-12 g cm-3) distributions are almost identical for the strongly and weakly magnetised cases.

3.1.4. Regions of active ambipolar action

Figure 5 also shows that for R> 100 au the iMHD and AD runs look similar. This resemblance arises from the strong dependence of ambipolar diffusion efficiency upon density and the sharp transition between the flux-freezing and AD dominated regimes. The initial isothermal phase of the collapse thus remains very similar in both cases. To examine the effect of flux dissipation in detail, we now focus on the regions where ambipolar diffusion dominates.

These regions are highlighted in Fig. 6a, which shows maps of Am, an adimensional number that characterises the efficiency of the diffusion process compared to the dynamical ones, defined as (23)The typical length-scale L is taken as the distance to the protostar and the velocity V as the local velocity along the field lines, V = (v·B)/ ||B||. Ambipolar diffusion essentially plays no role in regions where Am > 1, which is the case for r > 100 au, where niMHD and iMHD calculations produce identical structures. The region of the outflow, in contrast, is a region of very active ambipolar diffusion (Am < 1) because of a lower density that corresponds to a higher resistivity ηAD (see Eq. (2)). Dissipative effects are also very strong around the mid-plane for r < 50 au. This reduces magnetic braking by relaxing the split monopole configuration close to the dense core. The right panels in Fig. 6 show the magnetic field distribution as a function of density, along with the Am number (right axis). All cells with B > 0.08 G (horizontal dot-dashed line) have Am values coloured in red (compared to non-flagged cells, which are grey). We note that almost every cell belonging to the magnetic field plateau is indeed dominated by ambipolar diffusion effects rather than dynamical effects (Am ≲ 1).

μ = 2:

the same diagram for the more magnetised case is shown in Fig. 6b. Dynamical motions are less dominant because of the additional magnetic support. As a consequence, the active ambipolar diffusion region is stretched into an hourglass shape with a pinched equatorial waist in the mid-plane for r < 50 au, of about the same size as for μ = 5 . The pinching of the field lines is also weaker in this central region than for iMHD. However, compared to , the split monopole configuration remains, with a more radial field in the midplane, especially in the domain 50 < r < 150 au. This eventually leads to increased magnetic braking, which significantly hampers the formation of large Keplerian disks. Cells belonging to the diffusion plateau are, as previously, dominated by ambipolar diffusion effects and characterised by Am ≲ 1.

3.2. Misaligned case

The effects of tilting the rotation axis of the molecular cloud dense core with respect to the orientation of the magnetic field has consequences on the properties of the first Larson core and the accretion disks that can form around it (see Joos et al. 2012). We have repeated this analysis in MHD calculations including ambipolar diffusion.

3.2.1. First Larson core

The properties of the first Larson core in the misaligned case are summarised in Table 2. As in Sect. 3.1.1, the listed quantities were measured 200 years after peak density reached 10-12 g cm-3.

In the misaligned case, the regulation of the magnetic flux occurs as in the aligned case, yielding similar properties for the mass-to-flux ratio repartition and the value of the strongest magnetic field. For the collapse phase, the additional magnetic tension leads to a slightly longer free-fall time. For the weakly magnetised , the mass of the core is increased in the misaligned case. Indeed, the mass-to-flux ratio reaches the same value, but the magnetic field is twice stronger than in the aligned case, yielding a core mass almost twice higher.

Therefore, misalignment of the magnetic field and rotation axis does not change the properties of the first core significantly. It affects the formation and properties of rotationally supported structures, however, because of the weaker magnetic breaking, as we examine in Sect. 4.2.

3.2.2. Magnetic field repartition

thumbnail Fig. 7

Same as Fig. 1 for the misaligned case.

Open with DEXTER

Magnetic field repartition as a function of density is qualitatively unchanged compared to the aligned case, as seen in Fig. 7. The same non-linear self regulation process yields a diffusion plateau. The region above and below the core is strongly depleted, leading to a highly magnetised and low-density region (ρ< 10-16 g cm-3, B> 10-2 G) that corresponds to the polar cavity mentioned in Tomida et al. (2015). It is more extended in the more magnetised case (see Fig. 7a compared to Fig. 7b), as in the aligned configuration (see Figs. 6c and d). In iMHD, flux freezing still yields high peak magnetic field values, but the mixing due to the misalignment increases the numerical reconnection. As a result, a plateau somewhat similar to the niMHD case starts to develop (especially in the μ = 2 case). We will discuss this in further detail in a forthcoming paper (Paper II), which will include the effect of turbulence.

3.2.3. Outflows

We consider only niMHD in this paragraph. Because of the misalignment angle, outflows, when present, exhibit a wider opening angle in the misaligned case than in the aligned one, and their growth is hindered. Misalignment yields less pile-up of the toroidal field, which weakens the outflow-launching mechanism. As for the aligned case, we do not see any outflowing motions in the μ = 2 calculations.

We conclude that the presence or absence of outflowing gas is closely linked to conservation of angular momentum and to the regulation of magnetic field accumulation (and thus magnetic braking). Low magnetisation (μ = 5) allows for enough toroidal field accumulation to grow a magnetically and density-enhanced tower, while higher magnetisation (μ = 2) produces too effective braking, preventing the formation of magnetically launched and supported outflows.

3.2.4. Region of active ambipolar action

Figure 6c shows the efficiency of the ambipolar diffusion along with density contours and magnetic field orientation for . The results are similar to the aligned case (cf. Fig. 6a), with a larger pseudo-disk5. At these mid- to large scales (>50 au), ambipolar diffusion has little influence (Am ≫ 1). The polar cavities above and below the core are still present and well defined, and the region of ambipolar action in the mid-plane leads to the formation of a rotationally supported structure. Again, the magnetic plateau that forms at B ≲ 0.1 G corresponds to regions where Am ≲ 1, dominated by ambipolar diffusion.

The misaligned case is almost indistinguishable from the aligned one, as can be seen when comparing Figs. 6b and d. This suggests that in more magnetised cases the final configuration of the protostellar system results essentially from a saturation process (self-regulation from the ambipolar diffusion) and is independent of the alignment or misalignment. Another way to formulate this result is that all the relevant regions, including the diffusion plateau and most of the simulated box, are dominated by ambipolar diffusion (Am ≲ 1).

4. Long-term evolution in non-ideal MHD: disk formation

In this section, we examine the rotationally supported structures that can form around the first Larson core. We focus on niMHD calculations; iMHD calculations are discussed in Sect 5.

thumbnail Fig. 8

Visualisation of the disk structures. Top row: panels a) and b) show side and top views, respectively, of the disk density (orange), with the velocity vectors superimposed for μ = 5. Panels c) and d) are the same side and top views, but showing the plasma β parameter (colour map), over which we plot the magnetic field direction for the side view c) and the velocity vectors for the top view d). Bottom row: same as for the top row, but for the μ = 2 simulation. Each row has a different spatial scale.

Open with DEXTER

4.1. Aligned case

Conservation of angular momentum during the collapse yields a very high specific angular momentum close to the protostar, which eventually leads to the formation of rotationally dominated structures with a strong toroidal velocity component and in some cases to genuine disks with Keplerian velocity profiles. To define a “disk” in the simulations, we used the criteria defined in Joos et al. (2013). Structures in which the magnetic pressure either dominates or is at least not negligible, however, may satisfy these criteria while being very different from the flat rotationally supported structures characteristic of Keplerian disks. It is important to clearly identify and distinguish these two types of structures, since the existence of flat disks around Class-0 objects remains a subject of heated debate, and this issue is addressed below. Following Joos et al. (2013), a piece of fluid belongs to the disk if it fulfils all the following constraints (using f = 2):

  • it is close to hydrostatic equilibrium:vθ>fvz;

  • it has strong rotational support: vθ>fvr;

  • its rotational support is stronger than the thermal support: ;

  • it has a high density: ρ> 3.8 × 10-15 g cm-3.

In our fiducial aligned case with , we do observe the formation of a disk according to the above criteria that strongly resembles a Keplerian disk. Its characteristic properties are given in Table 3 and the results of the simulations are shown in Fig. 8 (top row). The disk is well defined in density, beginning just at the edge of the core and exhibiting sharp edges at its periphery. The breaking of symmetry observed in the top view (Fig. 8b) is due to the evolution over several orbital periods at the core radius, with small numerical errors leading to a bar-like instability that evolves into two spiral arms. The disk then grows to a radius of r ~ 80 au and is flat. The velocity profile in the disk is displayed Fig. 9 and is very close to a Keplerian one, with an estimated mass for the central object of 0.17 M compared to 0.23 M (Table 3) (see Appendix C for an important remark about assuming a Keplerian profile to estimate the mass of the central object).

Table 3

Properties of the disks that are formed in the non-ideal MHD simulations.

The evolution of the disk for our two values of magnetisation is displayed in Fig. 10 (solid lines for μ = 5, dashed lines for μ = 2). In niMHD, the pinching of the field lines the midplane is greatly reduced and the radial component of the field is almost null. In this case, there is almost no magnetic flux accreted onto the central object in the radial direction.

Figures 8c and d show the plasma β, defined as the ratio of thermal over magnetic pressure, β = P/Pmag, in the disk and its surroundings. We note that everywhere inside the disk and the core, β ≫ 1, meaning that the thermal pressure dominates magnetic pressure in these regions.

These results are of prime importance both from a physical and numerical point of view. Physically, they bear major consequences on our understanding of fragmentation and angular momentum transport in disks. In ideal MHD calculations, magnetic fields have been shown to inhibit fragmentation (Hennebelle & Teyssier 2008 and Commerçon et al. 2010). The fact that in more realistic non-ideal MHD calculations β ≫ 1 suggests that magnetic fields have less impact than expected according to iMHD results and that fragmentation may be facilitated. A second important consequence of the present calculations, of the characterisation of plasma β, and of the detailed topology of the field is the major role of these quantities in determining viscous transport in disks, in particular for the magneto-rotational instability (Balbus & Terquem 2001). A better knowledge of these quantities will enable us to define more accurate initial conditions in such studies (see e.g. Lesur et al. 2014).

thumbnail Fig. 9

, aligned case. Toroidal (vθ) velocity as a function of radius for disk cells. The pink dashed line is the fit (in the log/log space) of the blue points by a power law axb. The grey dashed line is the same fit but by enforcing b = −1/2. The vertical dotted line at 7 au is the inner radius for the fit.

Open with DEXTER

thumbnail Fig. 10

Disk mass evolution. Solid lines: μ = 5; dashed lines: μ = 2. Black: aligned case. Red: misaligned case. t = 0 corresponds to the first core formation to allow direct comparison between cases.

Open with DEXTER

Form the numerical point of view, the present studies are very important for the use of sink particles in simulations. Current sink particle models do not treat magnetic flux transport at the sink-gas interface in a consistent way. The flux is not accreted, which creates an effective barrier of infinite diffusion and spurious numerically driven interchange instabilities. As mentioned above, in the AD case, the radial component of the magnetic field is almost non-existent, yielding a configuration of zero net radial flux, allowing the accurate characterisation of the accreted (or in this case not accreted) magnetic flux onto the sink particle.

μ = 2:

in the more magnetised case, a disk-like structure around the first core is still observed, according to the above criteria, but the disk mass is an order of magnitude lower than in the less magnetised case. Figures 8e and f represent the disk cells in an edge-on and in a top view. As for the less magnetised case, the disk inner radius coincides with the outer layers of the core, and the disk slowly grows as matter is accreted. The final radius is about 20 au with a mass of 0.03 M. The aspect ratio remains close to unity, and the structure exhibits the same characteristic properties as in the previous case: it has a close-to-Keplerian velocity field with a high plasma β in the rotationally dominated structure, while in contrast, the regions outside the disk are magnetically dominated. No spiral arms have developed in this case, even at the end of the simulation, but steady state has not been reached yet as the disk is still slowly accreting.

thumbnail Fig. 11

Same as Fig. 8, but for the misaligned configuration. Each row has a different spatial scale.

Open with DEXTER

As an important global remark, we note that since values of β ≫ 1 in the disk are found in all our simulations in niMHD for any initial condition, a β-based criterion could thus be used as a robust disk criterion in niMHD simulations6.

4.2. Misaligned case

In the misaligned case, a disk can form more easily and is more massive, as already noted in previous studies (Joos et al. 2012; Li et al. 2013). Indeed, misalignment yields a thicker pseudo-disk, which reduces the magnetic braking. Figure 11 (top row) shows the disk structure at the end of the simulation. It exhibits a flat inner morphology, with a Keplerian velocity profile and high plasma β values. In this case we note large accretion spiral arms (with significant infalling motions) above and below the disk plane, with β ≪ 1. These spiral arms, however, do not represent a significant part of the disk mass. We note that using a β-based disk selection criterion would exclude these arms from the disk (see Fig. 11d compared to 11b).

For (Figs. 11eh), the early evolution is similar to the aligned case, the disk is at first hardly distinguishable from the core (outer radius is about 10 to 15 au) and then grows by accretion. Significant differences compared to the aligned case, however, occur at later stages, as seen in the figures. In this case, the disk is massive enough to trigger the formation of spiral arms, with a mean radius of 20 to 50 au. As for the previous cases, the disk is very well defined by the region of high β plasma.

4.3. Emerging consistent picture

Figure 12 shows a direct visual comparison of disks in each simulation. The disk formation and evolution during the collapse in the aligned and misaligned cases is presented in Fig. 10. In the misaligned cases, the disk forms at the same time as or even before the first core forms. This confirms, as discussed above, that the disk emerges from the outer layers of a distorted first core embryo. After about a hundred years, the disk is well defined, with a mass ≳ 10-3M in all cases. At this stage, the disk growth histories for the aligned and misaligned cases are essentially indistinguishable.

Figure 13 displays the angular momentum evolution in spatially restricted regions, separating the disk and outflow components in the system. The disk plane includes every cell fulfilling the disk criteria inside a sphere of radius 100 au. Conversely, the outflow component corresponds to every other cell in a sphere of radius 100 au. The μ = 5 and μ = 2 cases exhibit the same pattern. The mass-averaged angular momentum is slightly larger in the outflow region in the aligned case (dashed blue line) than in the misaligned one (dashed red line) during the first thousand years after the core formation. In the disk plane, the definition of the disk and the presence of spiral arms introduce variability, but the trend remains similar for the aligned and misaligned configurations. The misaligned case in iMHD releases the strong gradients in the pseudo-disk (strongly pinched magnetic field lines) found in the aligned case and thus enables larger disks to form. Non-ideal MHD produces the same smoothing of gradients, yielding eventually similar structures in the aligned and misaligned cases.

This similarity between aligned and misaligned simulations in niMHD shows that the formation and evolution of a disk is independent of the initial misalignment. Ambipolar diffusion, when it is efficient, regulates the angular momentum transport and the magnetisation during the collapse. A stronger magnetic braking will operate in the aligned case, leading to more angular momentum being evacuated in the vertical direction and ultimately yielding the same disk mass and final angular momentum once a physical equilibrium has been reached. The larger the ambipolar diffusion, the more efficient the regulation mechanism, leading eventually to similar structures. This is well illustrated by the even more pronounced similarity between the aligned and misaligned cases for μ = 2.

thumbnail Fig. 12

Representation of disk sizes and shapes for every niMHD simulation. The density colour map and contours are the same as in Figs. 8 and 11.

Open with DEXTER

thumbnail Fig. 13

Angular momentum evolution. Blue and red: μ = 5; grey and black: μ = 2. Solid lines: angular momentum in the disk plane; dashed lines: angular momentum in the outflow (see text). t = 0 corresponds to the formation of the first core.

Open with DEXTER

5. Limits of ideal MHD

In this section, we highlight several limits of the ideal MHD approximation in the context of prestellar magnetised collapse and disk formation, justifying in passing the fact that we did not discuss iMHD simulations in Sect. 4.

5.1. Counter-rotation

In ideal MHD, the increase in the toroidal field component is due to the rotation of the gas around the forming protostar and is only hindered by numerical reconnection. The resulting magnetic braking can be efficient enough to completely stall the rotation of the core. During the time the information propagates outwards, the core can start to counter-rotate, creating a new braking, until eventually co-rotation with its surrounding is reached again. This is illustrated in Fig. 14, where we slightly increased the initial rotation in our fiducial case to better illustrate our purpose: μ = 5, α = 0.25, βrot. = 0.03. Time is evolving from top to bottom. To trace the rotation, we calculated the angular momentum for each cell with respect to the mean direction of rotation in a sphere of radius 200 au. The first and second columns show the negative and positive toroidal velocity (vθ) maps, respectively, in a iMHD simulation. The third column displays the positive toroidal velocity for niMHD7. Strong counter-rotation occurs in the iMHD simulation, starting near the core (second row), and propagates into the outflow (third row). There is no sign of counter-rotation in the niMHD run. The mechanism continues and generates zones that extend to very large radii with alternate positive and negative toroidal velocities, producing a butterfly-shaped outflow (bottom row). This has important consequences on the expansion of the outflow, which is narrower for iMHD (this is most visible in the third row) because the counter-rotation hampers the pile-up of the toroidal component of the magnetic field.

We carried out similar studies in the misaligned case and found out that counter-rotation still develops after formation of the first core, although with a more limited spatial extent than in the aligned case. We also examined lower initial magnetisations and found out that counter-rotation is greatly reduced when μ> 5.

thumbnail Fig. 14

μ = 5, aligned case. Four snapshots (time increases from top to bottom) from an ideal MHD simulation (two left columns) and the ambipolar diffusion case (right column) with the magnetic field initially aligned with the rotation axis. For the iMHD simulation, the angular momentum is plotted in the left or right panel, depending on the sign of the azimuthal speed vθ in cylindrical coordinates.

Open with DEXTER

5.2. Flux redistribution

thumbnail Fig. 15

Evolution of the disk mass (solid line) and the peak magnetic field value (dashed line) for the , aligned case, in iMHD. The four red dots marked 1 to 4 indicate the times at which the snapshots in Fig. 16 were taken. We recall that the formation of the first Larson core occurs at t = 24300 years (see Table 2). The number 1 indicates the beginning of the release of magnetic flux from inside the core.

Open with DEXTER

thumbnail Fig. 16

a) Evolution of the density/velocity map for the four times labelled 14 in Fig. 15. b) The three first columns represent the relative importance of each component of the magnetic field. The colour scale shows the value of each component as follows: black for high values (≳ 1), white for low values (≲ 10-1.5); blue and red (in-between ≳ 10-1). The rightmost column displays the plasma parameter β = P/Pmag. Snapshots for different times of interest in the iMHD μ = 5 aligned case, as labelled Fig. 15.

Open with DEXTER

Figure 15 illustrates the evolution of the disk mass and central magnetic field for iMHD. At t ≳ 26.5 kyr, the field intensity decreases by almost an order of magnitude, producing a drastic increase of the disk mass. This stems from the flux accumulation at the centre of the collapsing system because the magnetic flux is frozen with the flow, which produces a major interchange instability. After this event, the core environment is weakly magnetised and quite disorganised (looking similar to turbulent runs), enabling the formation of a massive rotationally supported disk. The evolutionary sequence is portrayed in Fig. 16 (top), where the times corresponding to snapshots 1 to 4 in Fig. 15 are indicated in Fig. 15. The thin disk-like structure in 1 and 2, with a mass Mdisk ≲ 5 × 10-2M, evolves into a thick structure (3 and 4) with strong magnetic support, although it does not fulfil the disk criteria we defined above. The violent release of magnetic flux produced by the interchange instability between epochs 2 and 3 breaks the top-down symmetry of the system and displaces the core by several tens of au (this is not visible in Fig. 15 since we have integrated the maps by always using the densest cells as the origin). Detailed maps of the different magnetic components and the plasma β are displayed in Fig. 16 (second and third rows) for snapshots 1 and 4. The symmetry breaking, as well as the toroidal support, are apparent in the figure. Throughout most of the structure, β is lower than unity, showing again that magnetic support is ubiquitous in the system. As examined in the previous section, for AD the pile-up of the field occurs outside the core and is much less extended (see Fig. 4, at r ~ 30 au), and is controlled by the physical, rather than numerical, resistivity. Furthermore, it does not lead to the growth of the pseudo-disk8. Fig. 17 compares the density/velocity and plasma β maps, viewed from above, in the iMHD (top) and AD (bottom) cases. We clearly see the contrast between the highly magnetised (β ≪ 1) disk-like structure that arises in iMHD and the rotationally supported disk that develops spiral arms with β ≫ 1 when ambipolar diffusion is taken into account.

thumbnail Fig. 17

a) iMHD. b) niMHD. , aligned case. View of the disk plane in a cylindrical volume of height h = 20 au. Left and right: density/velocity field and plasma β, respectively, with the velocity field.

Open with DEXTER

thumbnail Fig. 18

iMHD, , aligned case. View of the disk plane in a cylindrical control volume of height h = 20 au. First row: t ~ 25.6 kyr. Second row: t ~ 27 k years. From left to right, the same simulation with increased resolution: 6, 8, and 16 cells per Jeans length.

Open with DEXTER

5.3. Numerical convergence

Numerical convergence is always a problem in numerical simulations. We carried out a resolution study in the aligned case and, while good agreement was found between the various AD runs, convergence in iMHD was not satisfactory. Numerical diffusion effects, as well as the choice of Riemann solver to compute fluxes at the cell interfaces, prevent a clear identification of the source of this problem. This issue has been highlighted in Commerçon et al. (2010), Hennebelle & Fromang (2008) for the influence of different solvers, and explored in detail by Li et al. (2014b), to which we refer to for further information and discussions. As part of our resolution study, we performed new runs of our benchmark aligned μ = 5 case (hereafter denominated case-0), with a mesh resolution increased by a factor 2 (16 cells per Jeans length) and a resolution decreased to six cells per Jeans length (this latter is only applied for densities above 10-12 g cm-3, while keeping the original Jeans length sampling elsewhere).

By comparing the high-resolution iMHD calculation to case-0, we found a good agreement until t ~ 25 kyr (the disk mass is slightly lower by a factor 1.5 at this stage, but this probably stems from the uncertainties of the criteria used to define the disk when no clearly identifiable rotationally supported structures are present; see Sect. 4), at which point magnetic field oscillations at the centre of the core start to appear due to the flux freezing condition (earlier in the high-resolution run). This is reassuring in the sense that it implies that the resolution we have used in case-0 is high enough for numerical resistivity effects to be negligible, which agrees with the absence of a strong field accumulation at the boundary of the core. The higher resolution, however, further diminishes the numerical resistivity, and thus flux freezing and flux accumulation are increased, facilitating the development of instabilities, which now appear earlier in the simulation.

In the low-resolution run, the numerical resistivity is increased, and we observe the appearance of a diffusion plateau, similar to the one in AD runs, but developing at a slightly higher value (B ≃ 1 G). In this case, the disk mass is overestimated but the symmetry is preserved longer than in the high-resolution runs. We observe transient ejections similar to those produced by the interchange instability, however, that are due to a strong pile-up of the magnetic field close to the core. These ejections ultimately lead to a similar quasi-steady final state.

A visual comparison of the different resolutions (6, 8, and 16 cells per Jeans length) at two different times is shown in Fig. 18. As seen, the resolution strongly affects the simulation output when numerical resistivity is not negligible. This is of prime importance in turbulent simulations, where it is never clear that the turbulent reconnection accurately describes the sub-grid unresolved physics. Another parameter that leads to numerical diffusivity is the numerical method chosen to solve the induction equation. More details on this issue are given in Appendix B.

The numerical convergence is much better for niMHD. Figure 19 shows the disk mass evolution using different mesh resolutions and Riemann solvers. The long-term evolutions are similar, with both the same global variations and final disk mass. The differences also remain small when changing the solver from HLLD to HLL, as long as the resolution is ≳ 8 cells per Jeans length. We used the HLLD solver for all our simulations to ensure that ambipolar diffusion dominates the numerical one.

thumbnail Fig. 19

niMHD runs. Comparison of the disk mass in case-0 for different numerical resolutions using the HLLD solver (black) and HLL solver (red).

Open with DEXTER

The conclusion of our convergence study is as follows. While in ideal MHD the isothermal first phase of the collapse is accurately described, with the formation of the first core and the early growth of the magnetic bubble/tower, the later evolution strongly depends on the numerical resolution. The key issue is the sudden release of the magnetic energy accumulated in the core that leads to oscillations in the peak magnetic field value (see Fig. 15 dashed blue line). The release of this energy, depending on the numerical tools, ultimately leads to either a quasi-steady state in which a massive disk-like rotationally and magnetically supported structure forms (see Fig. 17 in the top panels) or, for HLLD, to the appearance of a high-velocity jet and high Alfvén speeds that are due to depleted cavities above and below the core (see Appendix B for a discussion of the Riemann solver). The long-term evolution of these structures, which directly inherit their properties from the angular momentum and magnetic flux repartition close to the core, is thus of highly questionable validity. In misaligned configurations, the pseudo-disk is thicker and the spurious effects due to very high gradients are less acute (see e.g. Joos et al. 2012). The interchange instability is less visible than in Fig. 18, but the flux release and the growth of a rotationally supported structure (as in Fig. 15) is still present in the simulation.

In niMHD, in contrast, neither the resolution (as long as it is good enough) nor the choice of the solver have significant effects on the final outputs. This brings confidence to our global study and reinforces our general conclusion that ambipolar diffusion in the protostellar first collapse plays a major role in regulating angular momentum transport and magnetic flux diffusion, not to mention preventing spurious instabilities found in ideal MHD.

6. Comparison with previous work

In this section, we compare our results with those of four previous studies that are representative of the exploration of star formation in non-ideal magnetohydrodynamics.

Krasnopolsky et al. (2010) first explored the impact of enhanced resistivities on disk formation. Their study was motivated by the fact that in ideal MHD, magnetic fields prevent the formation of large Keplerian disks, which were ubiquitous in hydrodynamical simulations. Resistive MHD was identified as a physical process able to limit this effect. Their setup uses an isothermal equation of state, does not account for self-gravity, and assumes axisymmetry. Their main results, compared to ours, are the following:

  • in their Figs. 58, the pinching of the field lines in the equator and therefore the split monopole configuration is released as the resistivity is increased.

  • they found that a resistivity ≳ 1019 cm2 s-1 is needed to form a substantial disk.

While they did not directly link the release of the split-monopole configuration to the formation of rotationally supported disks, they concluded that increasing the resistivity helps both to form disks and to reduce numerical artefacts through physical rather than numerical dissipation, and it also changes the topology of the field that is prone to reconnect close to the star where the pinching was strong. We agree with these conclusions and note that the resistivity above which they formed disks corresponds to or is slightly stronger than the typical values of ambipolar resistivity we used (see Marchand et al. 2015) that lead to the formation of disks. However, these authors assumed a constant enhanced resistivity in the induction equation in the Laplace operator. Consequently, they could not grasp any non-linear effect, which, as we showed, can lead to a saturation of the magnetic field and further facilitate the formation of rotationally supported structures.

Li et al. (2014b) continued the work started by Krasnopolsky et al. (2010) and studied the mechanisms of disk formation. While they focused on iMHD, they raised many of the questions we developed in Sect. 5. In particular, they insisted on the role of the warping and the reduced reconnection close to the protostar, and they emphasized that a correct treatment of magnetic flux accretion on the protostar (especially if a sink particle is used) and reconnection at the boundary of the disk/pseudo-disk are mandatory and are lacking in many previous studies. The problem of turbulent reconnection in iMHD was also discussed and placed in perspective. We agree with their findings and insist on the fact that numerical questions regarding flux conservation and reconnection are crucial for a proper study of disk formation.

Tsukamoto et al. (2015) presented the first second core calculations with both Ohmic dissipation and ambipolar diffusion. They performed three-dimensional simulations using an SPH code with a Godunov module and divergence cleaning methods for the induction equation. Their main conclusions are the following:

  • 1)

    a value βplasma ≫ 104 in a flattened first core;

  • 2)

    formation of a circumstellar disk (with a radius of 1 au) at the formation of the protostar;

  • 3)

    occurrence of a saturation plateau at B ~ a few10-3 G in the evolution of the magnetic field at densities 10-15<ρ< 10-14 g cm-3 followed by a sheet-like collapsing phase (in which Bρ1/2) until the formation of the first core;

  • 4)

    they did not observe disks around the first Larson core.

Their values for βplasma (point 1) are similar to ours (see Table 2). So far, we cannot comment on point 2 since we did not address the second core formation in this paper. When studying the evolution of the central magnetic field as a function of density, they only plotted the data at the centre of the system (while we showed every cell in the simulation), and they thus failed to report a saturation of the magnetic field that is seen in our Fig. 1. Through private communication, we have obtained the confirmation that they indeed observed a diffusion plateau that was due to ambipolar diffusion, but its origin or consequences are not discussed in their work. We note that they assumed a sheet-like accretion through the disk, whereas in our simulations accretion onto the first core mostly occurs in channels driven by the hourglass configuration. There is, however, no real evidence for their assumption. Finally, the lifetime of the first core in their simulations is comparable to the time of integration after the first core formation in our runs (of the order of a thousand years for μ = 4). However, they did not form any mid-scale (a dozen to a few dozen au) rotationally supported structure around the first core. Instead, in their Ohmic and ambipolar+Ohmic calculations, the flattened first core itself evolves into a rotationally dominated structure that they called a disk. These differences in the dynamics of the collapse highlight the difficulty of conducting reliable numerical calculations during prestellar core collapse, which involves a wide range of temporal and spatial scales, and the necessity of including all relevant physical processes and of using reliable numerical schemes. These studies, and the present ones, open the route to such explorations.

Tomida et al. (2015) performed a study similar to ours, focusing on the first core formation and the surrounding structures, with the addition of radiation transfer instead of a barotropic equation of state. They used only one value of magnetisation that is close to our low-magnetisation case. Their model with ambipolar diffusion also included Ohmic dissipation, and they did not study the effects of ambipolar diffusion alone. We agree with the following of their results:

  • a flattening of the first core, especially in their OA(Ohmic+ambipolar diffusion) run.

  • depleted polar cavities where the ambipolar action is dominant.

  • in the simulation labelled I (ideal MHD), they found counter-rotation of the core (as visible in their Fig. 4), as in the present Fig. 14, but explained it by the interchange instability mixing the physical quantities. While we also found counter-rotation in iMHD, it occurs before the interchange instability develops. Following Galli et al. (2006), we explain this counter-rotation by a strong pile-up of the magnetic field yielding enough magnetic braking (or torque) to completely stop the rotation of the core and spin it backwards. The general conclusion, however, remains similar: iMHD models are unable to accurately describe angular momentum transport after the instability has developed.

  • an increase of the mass-to-flux ratio up to or above a factor 10 due to Ohmic and Ohmic+ambipolar diffusion action, especially at the first core scales.

  • no splitting of the field lines in the midplane.

  • the use of a genuine physical dissipation scale helps numerical convergence and thus assesses the reliability of the results.

Some of their results are more difficult to discuss. For instance, their explanation for a higher first core mass in the resistive cases where the magnetic pressure and the thermal and rotational supports are higher. In our experience, the interplay between magnetic braking (yielding stronger or weaker rotational support) and removal of magnetic flux (yielding stronger or weaker magnetic support) strongly affects the mass of the core. Although it is possible that in their set of simulation, the iMHD models lead to less massive first cores, in our runs we found that the first core mass is higher in the iMHD framework for μ = 2, but somewhat similar to or lower than the niMHD value for μ = 5 (see Table 2). The physical explanation remains the same: less support in total, but it is not straightforward to predict the results over a broad range of initial conditions. Although their discussion of the disk is not expanded enough to allow detailed comparison with our results, the main result is similar: non-ideal MHD enables disk formation with no need to artificially relax the flux-freezing constraint or to invoke enhanced resistivities.

7. Conclusions

We have thoroughly explored the role of ambipolar diffusion in magneto-hydrodynamic collapse calculations of a dense molecular cloud core in the context of star formation. Our setup involved a magnetised core of uniform density in solid-body rotation, which collapsed under its own gravity. The gas had a barotropic equation of state to mimic radiative transfer, and the magnetic resistivities entering the calculations of the non-ideal MHD terms were calculated using a reduced chemical network. We performed eight simulations, with two different magnetisations (μ = 5 and 2) and two tilting values of the rotation axis ( and 40°) with respect to the magnetic field direction, both in ideal MHD and in non-ideal MHD with ambipolar diffusion. We paid particular attention to the problems of magnetic flux conservation, magnetic braking, and flux release due to diffusion processes. Our main findings are summarised as follows.

  • Ambipolar diffusion creates a magnetic diffusion barrier at about the time the first Larson core forms, preventing the magnetic field from being amplified above 0.1 G. Flux freezing, however, still holds during the initial stages of the collapse, when resistivities are low. The mass and radius of the first Larson core remain rather similar between iMHD and niMHD models.

  • The magnitude of B at the diffusion plateau can be estimated by simple order-of-magnitude arguments. This saturation value appears to depend very weakly on the initial cloud magnetisation, suggesting a convergence of the final state once the initial conditions have been “forgotten” by the system.

  • The occurrence of a diffusion plateau has crucial consequences on magnetic braking processes. Not only does it prevent a catastrophic amplification of the field B, which controls the braking efficiency, but it also reorganises the field topology, reducing the pinching of field lines close to the protostellar object, which also hinders the braking mechanism.

  • Magnetic flux freezing and magnetic braking play a central role in the formation and development of the structures during the collapse. While in iMHD strongly amplified fields launch powerful magnetically supported outflows, these latter are much weaker or even disappear in niMHD.

  • Misalignment between the initial rotation axis and the magnetic field direction does not appear to affect the niMHD results, showing that the physical flux dissipation due to ambipolar diffusion dominates the effects of initial configuration or numerical diffusion. For iMHD models, the additional mixing due to a tilted configuration also produces a diffusion plateau, similar to the ambipolar calculations.

  • The disks that form in iMHD, characterised by strong magnetic support and an inflated shape (due to toroidal field lines loaded by infalling matter) strongly differ from the flat disks that form in niMHD. Formation and long-term evolution of disks in ideal MHD, however, are of dubious validity. Furthermore, the excessive magnetic braking generates unphysical counter-rotation inside the outflow or the magnetic tower. Interchange instabilities develop at the interface between the core and the disk, producing a redistribution of the flux and displacement of the core and disrupting the top-bottom symmetry in the system. Mesh resolution is found to strongly affect the simulation results in iMHD. None of these effects is observed in the niMHD simulations if the resolution is high enough.

  • Disks with Keplerian velocity profiles around the protostar form in all our niMHD simulations for all different magnetisations and inclination angles. Their size and mass, however, is significantly reduced (by a factor ~ 10) in the more magnetised case (μ = 2) because of the increased braking. Such a magnetisation value seems to be typical of most molecular clouds. Aligned and misaligned initial configurations have no consequence on the disk properties in niMHD and yield disks with very similar properties.

Magnetic flux diffusivity due to ambipolar diffusion thus appears to play a dominant role during the first protostellar collapse and the formation of the first Larson core and its surrounding disk, and it appears to dominate processes such as magnetic field and rotation axis orientations. Flux diffusion during the collapse allows the formation of quasi-Keplerian disks, solving the “disk formation crisis” found in ideal MHD. The mass, size, and magnetic properties (plasma β) of these disks, however, strongly depend on the initial magnetisation (μ parameter) of the cloud, and in all cases the disks are significantly smaller and less massive than those found in pure hydrodynamics calculations. Characterisation of this plasma β is of major importance to characterise viscous transport in disks, notably in the MRI. To complete our study, we will explore the effect of turbulence on the first collapse of a magnetised cloud core in a forthcoming paper.


1

See Appendix C for details on our definition of a Keplerian disk.

2

The initial condition corresponds to a ratio of the thermal over gravitational energies α = 0.25 and a density of 9.4 × 10-18 g cm-3.

3

This is chosen to try to reproduce the dragging-in of field lines that would occur in the formation of the dense core (see Gillis et al. 1974, for example), while also retaining in the simplest manner the divergence-free condition for the MHD.

4

In a simple representation of bent field lines, it is possible to derive the angular momentum transport equation. Angular momentum propagates along the field lines and follows a wave equation with a velocity defined by the local Alfvén speed and the topology of the field. For more details see Masson (2013) or the original article by Gillis et al. (1979).

5

The pseudo-disk is the overdensity region resulting from the loading of pinched field lines in a collapsing core that resembles a disk. See Galli & Shu (1993) for the first study of the pseudo-disk.

6

This criterion encompasses the core itself, which needs to be removed from the disk by considering the thermal to kinetic energy ratio, as done in Joos et al. (2013).

7

We do not display the counter-rotating vθ< 0 for AD since there are almost no counter-rotating cells except very close to the rotation axis at the very end of the simulation.

8

Growth of the pseudo-disk due to pile-up of the field has been studied in Hennebelle & Fromang (2008) in the appendix with a model based on quasi-equilibrium growth of a magnetised self-gravitating rotationally dominated structure in iMHD.

9

Discontinuities either in density, or when the wave speed for the intermediate states for HLLD are ill defined due to a very small denominator.

Acknowledgments

We thank the anonymous referee for the suggestions and remarks that contributed to improve the quality of this manuscript. The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 247060). B.C. gratefully acknowledges support from the French ANR Retour Postdoc program (ANR-11-PDOC-0031). We finally acknowledge financial support from the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France.

References

Appendix A: Calculation of the saturation value

In this section, we provide an analytical estimate of the value of B at which the magnetic diffusion starts to overcome the dynamical MHD effects. We rewrite the initial values for B0 and ρ0 in terms of the ratio of thermal to gravitational energy α and mass-to-flux ratio μ, according to The relationship between the magnetic field and density is assumed to follow ρBξ. We also assume that in the early phases of collapse, the typical length scale is Lcstfree - fall and the velocity is Vcs. We then seek the value of the magnetic field Bsat for which (A.3)with . Using the fact that (from Shu 1992), we can then write (A.4)which yields the saturation value for the magnetic field (A.5)

Appendix B: Influence of the solver

thumbnail Fig. B.1

Density maps of the region around the first core using four different Riemann solvers: a) Lax-Friedrich; b) HLL; c) HLLD+Switch; d) HLLD. The velocity vectors are overlaid on density maps with green and blue arrows.

Open with DEXTER

The choice of the Riemann solver to compute the conservative fluxes at the interfaces between domain cells can greatly affect the results of a simulation, especially when studying diffusionprocesses. We compare here four test simulations that use different strategies to compute the fluxes. The resolution for each run is 8 cells per Jeans length. In Fig. B.1 from left to right, the solver is increasingly less diffusive.

In the HLLD+switch simulation, we used the HLLD solver over most of the computational domain, but switched to a Lax-Friedrich solver whenever large discontinuities9 are present in the flow variables between neighbouring cells. This allowed us to evolve the simulations for long periods of time, up to ~ 1.5 free-fall time. However, it raises crucial problems on both the numerical reconnection in the vicinity of the core and the consequent interchange instability, as well as a violent core displacement that accompanies an unphysical release of magnetic energy at the centre. The formation of large inflated disks is also questionable.

In a second calculation, we removed the switch to a more diffusive solver, keeping HLLD in the entire domain. In this case, the discontinuities (mostly magnetic and density related) continue to sharpen, until the time integration in the simulation is frozen due to very high Alfvén speeds in the outflows or jets (it can reach speeds of the order of ≳ 104 km s-1). Whether this is a correct behaviour is open to debate (see also Li et al. 2014b). We reach convergence of results quicker than in the HLLD+switch case, since we basically prevent numerical diffusion, provided we have enough resolution at a given scale and density. On the other hand, there are physical diffusive processes at play in star formation, either from microphysics (e.g. ambipolar diffusion) or from turbulent reconnection. Therefore, a description that avoids any diffusion lacks physical mechanisms. Finally, we performed simulations using only HLL that takes into account only three waves instead of five with HLLD. This case fits in between the two previous ones; convergence is reached with 16 cells per Jeans length, but we do witness diffusion in the outflow and the vicinity or interior of the Larson core. In this case, Alfvén speeds remain below 103 km s-1. For simulations of the second collapse and formation of the second core, for which time integration is critical, our former tests suggest that using HLL with high enough resolution, at least in ideal MHD, is better than using HLLD (with or without the switch) because it enables both numerical convergence and the probing of longer timescales. For the sake of comparison, we performed one test case using the Lax-Friedrich solver alone. The results are similar to the HLL case.

Additionally, Fig. B.2 shows the first core mass evolution for each simulation, highlighting significant differences in an outcome as well defined as the first core mass. Until the time when the least diffusive solver (HLLD) does not allow continuing the simulation anymore because the time step is too small, the results are similar. However, after this point, depending on the choice of solver, results can significantly differ both in the qualitative picture (as seen in Fig. B.1) or for a precise outcome (Fig. B.2). Therefore, depending on the timescale of interest and the precise point of study, a careful choice of the numerical method has to be made.

thumbnail Fig. B.2

First core mass evolution for various solvers. Standard resolution of 8 cells per Jeans length.

Open with DEXTER

Appendix C: Disk velocity profile using ideal MHD

In this appendix, we caution about the definition of a Keplerian disk and the estimate of the central core mass obtained when fitting a priori the disk radial velocity distribution with a Keplerian profile (r− 1/2), in particular in structures whose growth is

magnetically driven, as is the case for pseudo-disks in ideal MHD. The disk criteria we used, based on Joos et al. (2012), define a rotationally dominated structure. We then proceeded by studying the radial dependence of the toroidal velocity, and compared it to the classical Keplerian velocity profiles. The velocity profile, the plasma β, and the aspect ratio allowed us to estimate how similar the disks we found are to a Keplerian disk. As emphasised below, we find that these disks can depart significantly from the classical Keplerian picture. The relatively high mass of the surroundings of the first core and the magnetic fields cause the usual Keplerian model to fail. This issue is of prime importance when trying to characterise the properties of the disk structures observed around Class-0 objects (see e.g. Tobin et al. 2013, 2015).

An example is shown in Fig. C.1. Whereas a Keplerian profile (grey dashed line) can be roughly fitted to the toroidal velocity field (blue dots), the best-fitting estimate (dashed red line) gives a significantly different value, r-0.32, yielding a central mass different from the one inferred from the model of a Keplerian velocity profile around a central point mass. This result stresses the fact that fitting a scattered velocity profile, which can depart significantly from a Keplerian one, by such a value can be very uncertain, casting doubts on estimates of central masses obtained with Keplerian formulae.

thumbnail Fig. C.1

Same as Fig. 9 for iMHD.

Open with DEXTER

All Tables

Table 1

Summary of the initial conditions for the different simulations.

Table 2

Summary of the main properties of the first core for fiducial cases with solid-body rotation.

Table 3

Properties of the disks that are formed in the non-ideal MHD simulations.

All Figures

thumbnail Fig. 1

Magnetic field distribution as a function of gas density for aligned, using the iMHD (red) and the niMHD (blue) formalisms, 200 years after first core formation. The solid black line shows the scaling of the magnetic field Bρ2/3 and the dashed black line the scaling Bρ1/2.

Open with DEXTER
In the text
thumbnail Fig. 2

a) Early time, when flux freezing is still relevant. b) After the decoupling stage. B(ρ) for (red) and (blue).

Open with DEXTER
In the text
thumbnail Fig. 3

Saturation value as a function of μ for various values of thermal support α. The points corresponding to μ = 2 (lower values) and μ = 5 (upper values) are marked in the figures.

Open with DEXTER
In the text
thumbnail Fig. 4

, aligned case. Mass-to-flux ratio () after first core formation (at ~ 24 kyr) for the niMHD (black) and IMHD (red) cases (solid lines). We also display the initial condition (dotted blue line) and two later outputs for the niMHD case (dashed and dotted solid black lines). The decoupling stage region due to magnetic field pile-up is emphasised with a green square.

Open with DEXTER
In the text
thumbnail Fig. 5

Top row: , aligned case, snapshots at a time t = 1.02 tfree - fall (500 years after first core formation). Panel a) shows a map of the density in the ideal MHD simulation. The black solid lines are logarithmically spaced density contours. The arrows represent the velocity field where blue to green corresponds to the increasing magnitude of the velocity in the (r,z) plane. Panel b) is a map of the Alfvén speed (blue), with the magnetic field direction in the (r,z) plane indicated by the segments. Yellow to red in these segments corresponds to increasing relative magnitude of the toroidal component of the field compared to the total field. Panels c) and d) are the same as a) and b), respectively, but for the niMHD case. Bottom row: same as the top row for μ = 2, and at the same time t = 1.01 tfree - fall (500 years after first core formation). Every quantity is azimuthally averaged.

Open with DEXTER
In the text
thumbnail Fig. 6

Four main panels show the aligned (top left, a), misaligned (top right, b), aligned (bottom left, c), and misaligned (bottom right, d) cases 500 years after first core formation. In each main panel, subpanel a) shows a colour map of the azimuthal average of . Green to blue means that dynamical processes (collapse, Keplerian disk) dominate, while red to yellow means that ambipolar diffusion dominates. The grey segments show the direction of the magnetic field, and black solid lines are isodensity contours. Panels b) and c) display the magnetic field magnitude B and the value of Am as a function of density. Red cells in the Am scatter plot are cells where B(ρ) > 0.08 G, while grey cells have B(ρ) < 0.08 G.

Open with DEXTER
In the text
thumbnail Fig. 7

Same as Fig. 1 for the misaligned case.

Open with DEXTER
In the text
thumbnail Fig. 8

Visualisation of the disk structures. Top row: panels a) and b) show side and top views, respectively, of the disk density (orange), with the velocity vectors superimposed for μ = 5. Panels c) and d) are the same side and top views, but showing the plasma β parameter (colour map), over which we plot the magnetic field direction for the side view c) and the velocity vectors for the top view d). Bottom row: same as for the top row, but for the μ = 2 simulation. Each row has a different spatial scale.

Open with DEXTER
In the text
thumbnail Fig. 9

, aligned case. Toroidal (vθ) velocity as a function of radius for disk cells. The pink dashed line is the fit (in the log/log space) of the blue points by a power law axb. The grey dashed line is the same fit but by enforcing b = −1/2. The vertical dotted line at 7 au is the inner radius for the fit.

Open with DEXTER
In the text
thumbnail Fig. 10

Disk mass evolution. Solid lines: μ = 5; dashed lines: μ = 2. Black: aligned case. Red: misaligned case. t = 0 corresponds to the first core formation to allow direct comparison between cases.

Open with DEXTER
In the text
thumbnail Fig. 11

Same as Fig. 8, but for the misaligned configuration. Each row has a different spatial scale.

Open with DEXTER
In the text
thumbnail Fig. 12

Representation of disk sizes and shapes for every niMHD simulation. The density colour map and contours are the same as in Figs. 8 and 11.

Open with DEXTER
In the text
thumbnail Fig. 13

Angular momentum evolution. Blue and red: μ = 5; grey and black: μ = 2. Solid lines: angular momentum in the disk plane; dashed lines: angular momentum in the outflow (see text). t = 0 corresponds to the formation of the first core.

Open with DEXTER
In the text
thumbnail Fig. 14

μ = 5, aligned case. Four snapshots (time increases from top to bottom) from an ideal MHD simulation (two left columns) and the ambipolar diffusion case (right column) with the magnetic field initially aligned with the rotation axis. For the iMHD simulation, the angular momentum is plotted in the left or right panel, depending on the sign of the azimuthal speed vθ in cylindrical coordinates.

Open with DEXTER
In the text
thumbnail Fig. 15

Evolution of the disk mass (solid line) and the peak magnetic field value (dashed line) for the , aligned case, in iMHD. The four red dots marked 1 to 4 indicate the times at which the snapshots in Fig. 16 were taken. We recall that the formation of the first Larson core occurs at t = 24300 years (see Table 2). The number 1 indicates the beginning of the release of magnetic flux from inside the core.

Open with DEXTER
In the text
thumbnail Fig. 16

a) Evolution of the density/velocity map for the four times labelled 14 in Fig. 15. b) The three first columns represent the relative importance of each component of the magnetic field. The colour scale shows the value of each component as follows: black for high values (≳ 1), white for low values (≲ 10-1.5); blue and red (in-between ≳ 10-1). The rightmost column displays the plasma parameter β = P/Pmag. Snapshots for different times of interest in the iMHD μ = 5 aligned case, as labelled Fig. 15.

Open with DEXTER
In the text
thumbnail Fig. 17

a) iMHD. b) niMHD. , aligned case. View of the disk plane in a cylindrical volume of height h = 20 au. Left and right: density/velocity field and plasma β, respectively, with the velocity field.

Open with DEXTER
In the text
thumbnail Fig. 18

iMHD, , aligned case. View of the disk plane in a cylindrical control volume of height h = 20 au. First row: t ~ 25.6 kyr. Second row: t ~ 27 k years. From left to right, the same simulation with increased resolution: 6, 8, and 16 cells per Jeans length.

Open with DEXTER
In the text
thumbnail Fig. 19

niMHD runs. Comparison of the disk mass in case-0 for different numerical resolutions using the HLLD solver (black) and HLL solver (red).

Open with DEXTER
In the text
thumbnail Fig. B.1

Density maps of the region around the first core using four different Riemann solvers: a) Lax-Friedrich; b) HLL; c) HLLD+Switch; d) HLLD. The velocity vectors are overlaid on density maps with green and blue arrows.

Open with DEXTER
In the text
thumbnail Fig. B.2

First core mass evolution for various solvers. Standard resolution of 8 cells per Jeans length.

Open with DEXTER
In the text
thumbnail Fig. C.1

Same as Fig. 9 for iMHD.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.