Free Access
Issue
A&A
Volume 615, July 2018
Article Number A5
Number of page(s) 18
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201732075
Published online 03 July 2018

© ESO 2018

1. Introduction

Angular momentum transport, and its regulation through magnetic braking, is one of the most important, yet poorly understood, physical mechanisms in star formation (e.g., Hennebelle & Charbonnel 2013). Under the ideal magnetohydrodynamic (hereafter MHD) approximation, magnetic fields typically observed in molecular clouds (Crutcher 2012) are powerful enough to remove all angular momentum from collapsing dense stellar progenitors: a problem known as the “magnetic braking catastrophe” (Matsumoto & Tomisaka 2004; Hennebelle & Fromang 2008; Hennebelle & Teyssier 2008; Mellon & Li 2008; Commerçon et al. 2010). Angular momentum is needed to form protoplanetary disks around young stars, and three possible solutions are currently being investigated by theoretical studies to try and solve the magnetic braking puzzle.

The first invokes the omnipresent turbulence in the molecular clouds, which, through turbulent reconnection, is thought to effectively regulate the concentration of magnetic flux and lead to the formation of protoplanetary disks (Santos-Lima et al. 2012, 2013; Leão et al. 2013; Lazarian 2013; Joos et al. 2013). Indeed, the first numerical studies of low-mass star formation were carried out in a rather simplified setup where the collapsing cloud was in solid body rotation, permeated by a uniform magnetic field. It has also been proposed that a disorganised field is simply less efficient at removing angular momentum from the system (Seifried et al. 2013, 2015). The second solution is once again related to the simulation setup; it is argued that the situation where the magnetic field direction is aligned with the parent body’s rotation axis is a very special case, with its own peculiarities, and unlikely to happen in nature. While the alignment between magnetic field and large density structures in molecular clouds has been studied with recent observations (Planck Collaboration Int. XXXV 2016; Hull et al. 2017), the spatial resolution does not allow to perform the same quantitative analysis at the cloud dense core level. It is however perfectly possible that rotation axis and magnetic field are misaligned, especially if the magnetization is weak (Mocz et al. 2017; Hull et al. 2017). Hull et al. (2013) present dust-polarization observations towards 16 nearby low-mass protostars and conclude that their data are consistent with disks that are not aligned with the magnetic fields in the cores from which they formed. This scenario was investigated by several authors (Hennebelle & Ciardi 2009; Joos et al. 2012; Li et al. 2013; Krumholz et al. 2013; Masson et al. 2016) and was found to also be conducive to disk formation. Nevertheless, we note that as the magnetic dissipation relies on numerical diffusion, these studies do not always yield resolution converged results in the ideal MHD framework.

Finally, resistive effects in the induction equation were suggested as a means to reduce the pile-up of magnetic field around the central object (Duffin & Pudritz 2008; Mellon & Li 2009; Krasnopolsky et al. 2010; Li et al. 2011; Machida & Matsumoto 2011; Dapp & Basu 2012). The gas inside protostellar envelopes and protoplanetary disks is poorly ionised, and ion-neutral collisions, which act as a diffusive process in the MHD equations, are omnipresent. While in early 2D studies, neither ohmic nor ambipolar diffusion were able to circumvent the magnetic braking catastrophe without requiring abnormally large resistivities (Krasnopolsky et al. 2010; Li et al. 2011), more recent 3D calculations have shown that magnetic diffusion with realistic resistivities can facilitate the formation of flat rotationally dominated structures, with radii of about 50–60 astronomical units (AU; Tomida et al. 2015; Tsukamoto et al. 2015a; Masson et al. 2016; Hennebelle et al. 2016)1. This third pathway provides a physical diffusion mechanism which does not depend on the numerical resolution or the orientation of the magnetic field, it is simply governed by the microphysics of molecular cloud.

The vast majority of the works listed above have studied the first hydrostatic core stage of star formation (scales of ~10 AU), and very few have considered the scales typical of the protostellar seed; the second Larson core (<0.1 AU; Larson 1969; Masunaga & Inutsuka 2000; Vaytet et al. 2013). The first full 3D hydrodynamical simulations of the formation of the second Larson core were carried out by Bate (1998). Since then, only a limited number of studies have reached the second core stage, with different numerical methods (nested grid codes, smoothed particle hydrodynamics), incorporating increasingly complex microphysics including magnetic fields, radiative transfer, magnetic diffusion. We summarise the list of these papers in Table 1. The recent works by Tomida et al. (2015), using a nested-grid code, and Tsukamoto et al. (2015a), using smoothed particle hydrodynamics, were the first ones to include radiative transfer coupled to MHD with both ambipolar and ohmic diffusion2. Even more recently, Wurster et al. (2018) went a step further by adding the Hall effect in their calculations of the second core formation. To help establish theoretical results, it is crucial to verify computational results across different codes and numerical methods. This paper aims to do precisely this, expanding on the latest Japanese and British studies to strengthen the validity of the star formation process. We follow the gravitational collapse of a dense sphere of magnetised gas, from molecular cloud densities to the formation of the protostar, including ambipolar and ohmic diffusion. We compare the results to the classical ideal MHD (IMHD) framework, and illustrate why magnetic diffusion is of paramount importance in low-mass star formation.

Table 1.

3D numerical studies of the formation of the second Larson core.

2. Numerical method and initial conditions

2.1. RAMSES with non-ideal MHD and flux-limited diffusion

The simulations were carried out using a modified version of the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002; Fromang et al. 2006) which incorporates the effects of ambipolar and ohmic diffusion (Masson et al. 2012), and radiative transfer via a time-implicit flux-limited diffusion (FLD) approximation (Commerçon et al. 2011b, 2014). The governing equations are(1) (2) (3) (4) (5) (6) (7)

The quantities are (in order of appearance): the gas density ρ, time t, the gas velocity v, the gas pressure p, the magnetic field B, the identity matrix 𝕀, the gravitational potential Φ, the radiative flux limiter λ, the radiative energy Er. The total gas energy is defined as E = ∈ + ρv  · v/2 + B · B/2 where is the internal gas energy. ηO and ηA are the ohmic and ambipolar magnetic resistivities, κP is the Planck mean opacity, c is the speed of light, ar is the radiation constant, while T represents the gas temperature, G is the gravitational constant, ℙr is the radiation pressure, and κR is the Rosseland mean opacity.

Equations (1)(3) describe the conservation of mass, momentum, and energy, respectively. Equation (4) is the induction equation, Eq. (5) is the divergence-free condition, Eq. (6) is the Poisson equation for self-gravity, and Eq. (7) is the conservation of radiative energy density. In this work, we used the HLL Riemann solver for the MHD, and the Minerbo flux limiter (Minerbo 1978) for the FLD which is defined as (8)

where R = |Er|/(ρκREr). The radiation pressure is given by ℙr = 𝔻Er, and the Eddington tensor is(9)

with χ = λ + λ2R2 and n = Er/|Er| (Levermore 1984). The code incorporates the gas equation of state of Saumon et al. (1995), and its extension to low densities (see Vaytet et al. 2013), for a mixture of hydrogen (73%) and helium (27%, in mass). The interstellar dust and gas opacities were taken from Vaytet et al. (2013). These comprise the dust opacities of Semenov et al. (2003; assuming a 1% dust content, by mass) at low temperatures (below 1500 K), the molecular gas opacities of Ferguson et al. (2005) for temperatures between 1500 and 3200 K, and the atomic gas opacities from the OP project (Badnell et al. 2005) above 3200 K. To aid the convergence of the implicit radiative transfer solver, we artificially limited the optical depth per cell to a minimum value of 10−4. When the gas is optically thin, it is not crucially important for the heating and cooling mechanisms whether the opticaly depth is 10−8 or 10−4, but we observed that choosing the latter can typically cut the number of iterations in the conjugate gradient solver by a factor of 4 or more. We show a validation of this acceleration scheme in Appendix A.

The magnetic resistivities were computed from a reduced chemical network including neutral and charged species, as well as dust grains, using an earlier version of the Marchand et al. (2016) model. It is in fact identical to the fiducial model of (Marchand et al. 2016; with a cosmic ray ionisation rate of 10−17 s−1) for densities below 10−8 g cm−3, but features a smooth decay in both η A and η O beyond this point, following Machida et al. (2007) who use this to represent the thermal ionization of alkali metals, instead of taking into account the effects of grain evaporation, thermal ionisation of potassium, sodium, and hydrogen, and grain thermionic emission. Using this tool, a three-dimensional table of density, temperature, and magnetic field dependent resistivities was computed. During the simulations, the resistivities in each grid cell were interpolated on-the-fly according to the local state variables, greatly reducing computational cost but implying thermodynamical equilibrium. The resistivities severely limit the integration timestep, and a stable super-time stepping method for ambipolar diffusion on an AMR grid with level-by-level sub-cycling is still lacking from the literature. To speed up the calculations, the timestep was prevented from going below a fraction of the ideal MHD timestep. It is taken to be the minimum of the three timescales:(10)

where Δx is the cell size, ξ = 0.1, and(11)

is the fast magnetosonic speed in direction i, where v A =  is the Alfvén speed, and the sound speed (12)

includes the contribution from the radiation pressure (see Commerçon et al. 2011b). The idea is that the exact amount of magnetic diffusion included is not crucially important, as long as some diffusion is operating (see Appendix B for more details). It is, however, necessary to compute the resistivity coefficients accurately with a chemical network, as in Marchand et al. (2016), as the densities and temperatures at which they either rise or fall are important. The mesh refinement criterion was defined so that the local Jeans length was always sampled with a minimum of 32 cells everywhere in the computational domain. Initial tests with lower resolutions yielded spurious heating between the first and second core stages, due to inefficient cooling (see Appendix C and Vaytet & Haugbølle 2017).

2.2. Simulation set-up

We adopt initial conditions similar to those in Commerçon et al. (2010). A magnetised isothermal sphere of molecular gas with quasi uniform density, rotating about the z-axis with solid body rotation, is placed in a surrounding medium a hundred times less dense with equal temperature. The sphere has a mass M0 = 1 M, a radius R0 = 2753 AU, and a temperature T0 = 10 K, for an initial ratio of thermal to gravitational energies of(13)

where kB is Boltzmann’s constant, µ is the mean molecular weight (=2.31 initially for the H2 + He mixture), and mH is the hydrogen atomic mass. The density in the domain is defined by(14)

where ρ0 = 6.76 × 10−18 g cm−3 and includes an m = 2 perturbation of amplitude Δρ = 0.1, which has been used in many of our previous works to favour fragmentation in the collapsing system (see Commerçon et al. 2008; Commerçon et al. 2010). The amount of rotation given to the cloud is parametrised according to the ratio of rotational to gravitational energies, which was chosen to be (15)

where Ω0 is the angular velocity. The strength of the magnetic field is defined in terms of the mass-to-flux ratio normalised by the critical value of stability for a uniform sphere(16)

where and (Mouschovias & Spitzer 1976) and is the cylindrical radius. The magnetic field is initially parallel to, and invariant along, the axis of rotation z. The field is stronger in a cylinder of radius R 0 (with the dense core at its centre) than in the surrounding medium, with Bz (r cyl < R 0) =  B 0 = 1002/3 Bz (r cyl > R 0), where the factor of 100 comes from the difference in density between the core and the surroundings (see Masson et al. 2016). The base grid at the coarsest level counted 643 cells, and an additional 21 AMR levels yielded a final effective resolution of 8×10−5 AU.

3. Results

We performed two simulations: the first using the ideal MHD approximation (runID), and the second including ambipolar and ohmic diffusion (runAO), requiring 40 000 and 180 000 CPU hours, respectively3. In the remainder of this paper, we focus on describing the differences between the two models.

3.1. Early evolution

The evolution of a gravitationally collapsing dense molecular cloud core has been described in detail in past works (see Masunaga & Inutsuka 2000; Vaytet et al. 2013, for instance), and is displayed in Fig. 1 for our two runs. It begins with an isothermal phase of contraction, clearly visible in the lower left corner of panel b, where the compressive heating is lost via radiative cooling. As the density rises, the system’s optical thickness increases and the radiative cooling becomes less and less efficient, until it can no longer counter-balance the compressive heating. The system enters its first adiabatic phase when densities exceed ~10−13 g cm−3, where the first hydrostatic Larson core is formed. The first core continues to accrete material from its envelope, and the sustained increase in mass forces the temperature to rise in the centre. When the gas reaches 2000 K, H2 molecules begin to dissociate. The effective adiabatic index drops below the critical value of 4/3 for support against gravitational contraction, and a second, very rapid, phase of collapse takes place, at the end of which the second hydrostatic Larson core is formed. The moment where the curves in all three panels exhibit a very sharp rise marks the onset of second collapse.

thumbnail Fig. 1.

Density (panel a), temperature (panel b) and magnetic field strength (panel c) as a function of time, for the densest cell in the system. The red lines represent runID, while the blue lines are for runAO. In the top panel, the two insets show maps of the logarithm of density in runID just before (panel d) and after (panel e) the development of the interchange instability (see text).

Open with DEXTER

In the early stages (t < 28 230 yr) runID and runAO have very similar central density and temperature evolutions. Only the strength of the magnetic field differs significantly already after 28 000 yr, because the ambipolar and ohmic diffusion strongly hinders the condensation of magnetic flux. Just before the second collapse in runAO (t ≃ 28 230 yr), the discrepancy in B has grown to almost 3 orders of magnitude. The effects of a strong field amplification are visible in the subsequent evolution of runID. All three displayed quantities show a plateau after 28 250 yr, where contraction and heating is halted, delaying the second collapse. As illustrated by maps of the gas density in insets (d) and (e), this is caused by interchange instabilities that develop in the presence of extreme gradients in the magnetic field (Spruit et al. 1995). This effect was already observed in other works (e.g. Zhao et al. 2011; Tomida et al. 2015; Masson et al. 2016), and is discussed further below.

3.2. Physical picture at the time of second core formation

We now turn to describing in more detail the properties of the first and second Larson cores, at a time right after the formation of the second core. Finding a moment in both simulations where all aspects and structures of the collapsing systems can be directly compared is not trivial. The two runs reach the second core stage at slightly different times, and with different densities and temperatures in their centres. We defined the formation of the second core as the moment when a fully formed accretion shock is present, with a sharp density and velocity gradient at the core border. The justification for this somewhat arbitrary criterion will become clear in the following paragraphs. In addition, in the remainder of this work, a density threshold criterion – favoured for its simplicity and robustness – will be used to define the first and second Larson cores (see Appendix E). All the cells with a density higher than 10−10 g cm−3 make up the first core, while the threshold is 10−5 g cm−3 for the second core.

We first look at the evolution of the gas temperature at the centre of the system as a function of density, represented by the dashed lines in Fig. 2a. The quasi isothermal contraction at low densities (<10−13 g cm−3) is clearly visible in the lower left corner. The curves then follow an isentrope with an almost constant adiabatic index γeff ≃ 7/54 until temperatures reach 2000 K and γeff falls to ~1.1, initiating the second collapse. The value of 7/5 is recovered towards the end of the tracks, once temperatures exceed ~104 K. The evolutions in runID and runAO are very similar, following tracks which strongly resemble the results of past 1–3D studies (Masunaga & Inutsuka 2000; Vaytet et al. 2013; Tomida et al. 2013; Bate et al. 2014, to only name a few). The colour maps in Fig. 2a show a single snapshot in time of the distributions in the (ρ, T) plane of all the cells in the simulation domain, just after the formation of the second Larson core. Red colours are for runID while blue is for runAO. The cells have been divided into two regions; the equatorial region (light colours) where the polar coordinate θ = cos−1(z/r) is in the range π/4 < θ < 3π/4, and the polar region above and below the central protostellar object where θ < π/4 or θ > 3π/4. The centre of the polar coordinate system is the centre of the second Larson core, found by calculating the mean coordinate of all cells with ρ > 10−5 g cm−3. The results from the two different calculations are overall qualitatively similar. The most noticeable difference is the density at which the shock heating occurs when the gas enters the second core. The shock heating happens at densities two orders of magnitude higher in runID than in runAO, suggesting that the protostellar core is more compact in the IMHD run. We also note that the gas in the polar regions (darker colours) undergoes shock heating earlier (i.e. at lower densities) than around the equator (lighter colours), suggesting that the gas reaching the second Larson core is more diffuse close to the poles. This is actually visible below, in the density map around the second core in Fig. 3r.

thumbnail Fig. 2.

Left column: temperature (panel a), magnetic field (panel b), plasma β(panel c), and ratio of thermal to radiative pressure (panel d) as a function of density, for every cell in the computational domain at the epoch of second core formation. The IMHD simulation is represented by the red colours, while the blue shades are for the NIMHD run. The green colours correspond to areas where IMHD and NIMHD results agree within 10%. Each data set is delineated by a solid contour line which outlines the data distributions. The dark and light colours give an indication of the positions of the cells in the simulation box according to the θ = cos−1(z/r) angle: the light colours denote cells close to the equatorial region (π/4 < θ < 3π/4) while dark colours show cells in the polar regions (θ < π/4 or θ > 3π/4). The dashed lines in panels a and b represent the time evolution of the central (densest) cell inside the mesh. The thin black line in panel b is the power law predicted from magnetic flux conservation in a contracting gas sphere. Center and right columns: radial distributions of various quantities for every cell in the computational domain. As in the left column, red colours are for runID while blue colours are for runAO. In panels g and h additional lines show the integrated enclosed mass and angular momentum, respectively, in successive spherical shells going outward from the centre of the system.

Open with DEXTER
thumbnail Fig. 3.

Slices through the centre of the domain, comparing the morphologies of the protostellar system in runID (columns 1 and 3) and runAO (columns 2 and 4) on three different spatial scales at the epoch of second core formation. Panels a–h: display a wide region around the first Larson core, the typical scale of a protoplanetary disk. Panels i–p: immediate vicinity of the first Larson core. Panels q–x: present the second Larson core and its close surroundings. Two left columns: side x–z views of the system, while the two right columns display the top x–y perspective. The coloured maps in each row alternate between representing the gas density and temperature. The arrows on the density maps depict the gas velocity field. Overlayed onto the temperature maps are magnetic field lines (left column) and AMR level contours (right column).

Open with DEXTER

Figure 2b shows the distributions of the magnitude of the magnetic field vector B = |B| as a function of gas density. At low densities (<10−15 g cm−3), runID and runAO yield identical results. Above this point, we observe the same behaviour as in Masson et al. (2016). While the magnetic field follows a B ∝ ρ2/3 power law in runID (consistent with magnetic flux conservation for a contracting gas sphere), a clear magnetic diffusion plateau appears in runAO around 0.1 G. This diffusion barrier strongly limits the amplification of the magnetic field, reduces magnetic braking, and prevents several IMHD peculiarities such as counter-rotation of gas inside the envelope surrounding the first core, or the development of interchange instabilities (see Masson et al. 2016). As the resistivities begin to drop above densities of ~10−8 g cm−3 (see Sect. 2.1), B rises once again, but will remain between one and two orders of magnitude below the IMHD values. This has very important consequences for the properties of the second Larson core.

The ratio of thermal to magnetic pressure, otherwise known as the plasma β = 2p/B 2 is displayed in panel c as a function of density. The effects of magnetic diffusion are once again unequivocal. At low densities, outside of the first core, the magnetic pressure dominates everywhere in both runID and runAO. It also mostly dominates (or is comparable to the thermal pressure) inside the first and second cores in runID. However, the thermal pressure is orders of magnitude higher than the magnetic pressure when magnetic diffusion is included, as was reported in Masson et al. (2016). The first and second hydrostatic cores are genuinely supported by thermal pressure, and the two simulations are forming two completely different protostars.

Panel d displays the ratio of thermal to isotropic radiative pressure P rad = E r/3, as a function of density. The two runs yield similar results. At low densities, radiative and thermal pressures are comparable, but as the gas contracts isothermally, P rad remains constant while p scales linearly with density. As a result, the thermal pressure vastly dominates virtually everywhere in the collapsing system.

We now turn to studying in panels e to l the distributions of the fluid variables as a function of radius. Panel f shows the gas density as a function of radius, and the distributions are relatively similar between IMHD and non-ideal MHD (NIMHD) models. The densities are in general lower along the polar directions than in the equatorial plane, which is expected for a disk forming in the plane of rotation. The second core in runID appears to be more compact than its runAO counterpart, and seems to also have a different structure; its density is relatively uniform, suggesting a more spherical morphology, while the runAO core is elongated in the equatorial plane and has density peaks away from the centre. The temperature distribution in panel e shows again the more compact nature of the runID second core. It also reveals that in runID, temperatures are higher in most of the computational domain. This includes the regions inside the second core (r < 0.003 AU), around the first core border (1 < r < 10 AU) and also at larger radii (r ~ 100 AU).

Panels i and j show the radial (vr ) and azimuthal (vϕ ) components of the gas velocity, as a function of radius. Two (negative) spikes in vr around 1 and 0.01 AU in runAO mark the first and second core borders, respectively. In runID, the first core border is less well defined and has a radius 3 times larger, while the second core is clearly visible around 3 × 10−3 AU. As expected, the highest velocities are found in the polar regions, where the gas is free-falling along the magnetic field lines, meeting no resistance along its path. The IMHD model has positive vr between 2 and 100 AU, representative of an outflow; a feature absent from runAO. The positive radial velocities inside the second core in runAO are a sign that the core is expanding because of strong rotation. Indeed, panel j shows a colossal amount of rotation in and around the runAO second core, while it is effectively zero in runID. The magnetic braking is so efficient in the latter that it has removed all angular momentum from the second core (this confirms the results of Tomida et al. 2013).

Panels k and l display the vertical (Bz) and toroidal (Bϕ) components of the magnetic field, divided by the magnitude of the B field vector. This reveals that around the first core region (0.5 < r < 50 AU), the field is much more vertical in runAO (Bϕ falls to zero), while the opposite happens in runID. The magnetic diffusion allows the field lines to remain vertical without being drawn in by the fluid, unlike the IMHD model where perfect coupling between fluid and magnetic field means that the field lines are dragged into a pinched hourglass shape (see Krasnopolsky et al. 2010, for example), changing the orientation of the field and strongly reducing Bz. The picture is almost reversed for the second Larson core, but for a different reason. The field is almost entirely toroidal in runAO (Bz/B → 0 and Bϕ/B → 1), because of the strong rotation of the gas which drags the field lines along (at these densities and temperatures, the gas is almost fully ionised and the field is once again perfectly coupled to the gas). On the other hand, Bϕ remains rather small in runID because of the lack of rotation at the second core level. We also note that throughout the domain, the field remains mostly vertical in the polar regions, in both simulations, which is fully expected in a set-up where the rotation axis is initially aligned with the magnetic field.

Finally, in panels g and h we show the distribution of the mass and angular momentum, respectively, contained in the grid cells. The mass contained inside a cell may not provide much valuable information, as it is governed by our mesh refinement strategy, and the fact that it varies only lightly across the entire radial extent is simply a result of choosing to refine the grid according to the Jeans criterion. More interestingly, if we integrate the mass inside successive spherical shells around the protostar, we obtain the enclosed mass which we represent by the two solid lines in the upper half of panel g. The two systems have similar mass profiles, apart from inside the second core which is more compact in runID. In the case of the angular momentum, the main difference between the two runs is a collection of cells in runAO with much higher angular momentum than in runID, in the range −3 < log(r) < −1.5. This corresponds to the cells with high azimuthal velocities found in panel j. As a consequence, the integrated angular momentum for radii below 1 AU is orders of magnitude higher in runAO than in runID. In fact, the exceedingly strong magnetic braking in runID even forced a sign reversal of the angular momentum inside the second Larson core (dashed red line). However, the amount of rotation is so small (see also Sect. 3.3) that it is difficult to see as a bulk counter-rotating motion; the main component of the gas velocity is radially infalling at these radii.

This result is of crucial importance. It shows that magnetic diffusion (both ambipolar and ohmic) starts to become effective for radii below 10 AU, and even more so below 1 AU, indicating that a spatial resolution of at least ~1 AU is necessary to correctly study angular momentum transfer and the formation of rotationally supported disks around protostars5.

3.3. Morphologies

Figure 3 contains multiple slices through the data, comparing the morphologies of the protostellar system in runID (columns 1 and 3) and runAO (columns 2 and 4) on three different scales. The top two rows display a wide region around the first Larson core, the typical scale of a protoplanetary disk. The two middle rows show the immediate vicinity of the first Larson core, while the bottom two rows present the second Larson core and its close surroundings. The two left columns show side x–z views of the system, while the two right columns display the top x–y perspective. The simulation times are the same as in Fig. 2.

3.3.1. The first Larson core and its surroundings

Panels a–d show gas density maps with velocity vectors. An equatorial density enhancement, typical of an accretion disk, is clearly visible in the side view of both simulations. In the top view, a filamentary structure extending from the north-west to the south-east of the protostar has formed from the initial density m = 2 perturbation6. A magnetic tower with outflowing velocity arrows (corresponding to the positive radial velocities in Fig. 2e) is observed in runID (a), while it is absent from runAO (b), as was the case in the strongly magnetised simulations of Masson et al. (2016). Another large difference between the two runs, and another sign of strong amplification of the magnetic field, is the presence of “bubbles” in the x–y view (c) of runID which are caused by interchange instabilities (see Zhao et al. 2011; Krasnopolsky et al. 2012, for a detailed study of these structures). While it has been argued that misalignement between the initial B field and the rotation axis and turbulence are both able to prevent the formation of such structures (Li et al. 2013, 2014), ambipolar and ohmic diffusion provide a physical rather than numerical diffusion that dominates the dissipation processes, with no dependence on the initial direction of the B field nor the numerical resolution. The aligned case is no longer a special set-up with its strange behaviours and artefacts (see Masson et al. 2016). Further evidence of the rearrangement of magnetic field lines provided by resistive effects is seen in the second row (panels e and f), where the magnetic field lines are very pinched in runID, while they are much more vertical in runAO. This corroborates our findings above; the field lines are no longer perfectly coupled to the gas and get less dragged in by the collapsing fluid. The modification of the magnetic field topology is provoked by the ambipolar diffusion, the dominant mechanism in this region (r < 30 AU; see Appendix D)7. The temperature maps are also markedly different, with runID showing higher temperatures everywhere around the central protostar, up to a radius of ~100 AU.

Taking a closer look at the first Larson core in panels i to p, we notice that the disk is “puffed” up in runID (i) compared to runAO. The top view (k) also clearly show gas ejections from the interchange instabilities with outflowing velocity vectors. When looking at the time evolution of the gas temperature, we found that a sudden heating of the gas around the first core coincides with the development of the interchange instabilities, although we have not been able to establish if the instability is directly responsible for the heating. Other possible explanations include shock heating from waves launched by the instabilities, or irradiation from the protostar which is enhanced because the density – and hence optical thickness – of the gas around the first core drops as it gets ejected. One could even envisage a combination of the two, where shock heating raises the temperature around the core above ~1000 K where dust grains start to sublimate, abruptly lowering the opacities, which in turn intensifies the irradiation.

In runAO, all the gas is moving towards the core, and the accretion is highly anisotropic, occuring primarily along the two high-density streams seeded by the perturbation in the initial conditions. In panels m and n, the contrast in magnetic field orientation is glaring; the field in runID is pinched to the extreme, while it has become almost vertical in runAO due to the resistive effects. Panel p shows the high-density accretion streams hindering the propagation of heat from the central source, which progresses instead along the perpendicular direction. In runID, the more homogeneous density structure leads to a more homogeneous temperature distribution. The magnetic reconnection that occurs when interchange instabilities develop may also provide additional heating. However, this is not reconnection enabled by ohmic diffusion (that generates Joule heating) since it appears in the IMHD simulation; it is known as numerical reconnection. We have not been able to determine whether numerical reconnection heating is significant (or even happening at all) when compared to the irradiation from the central object, but the gas heating does appear to coincide with the development of the bubble-like ejections.

3.3.2. The second Larson core

Panels q to t show once again density maps with velocity vectors for runID and runAO, but this time in the vicinity of the second Larson core. The morphologies are here also very different. The second core border is not very well defined in runID, where the gas density shows a rather smooth transition from 10−7 to 10−3 g cm−3, as was already found in Fig. 2f. The protostellar seed also displays a loss of top-down symmetry (q), most probably due to magnetic flux redistribution during the development of the interchange driven magnetic “bubbles”. We also note the absence of any rotation in panel s, as already mentioned in Sect. 3.2. On the other hand, the runAO second core has a sharp border, strong rotation and a preserved top-down symmetry. It is flatter around the poles, due to both the rotation and the high infall speeds in the polar direction. The top view (t) also reveals the early development of a spiral structure inside the core. The second core masses for runID and runAO are 3.8 × 10−3 M and 7.4 × 10−3 M , respectively.

The temperature maps with overlayed magnetic field lines in panels u and v expose the compact nature of the second core in runID. Temperatures at the very centre are higher than in runAO, and the core surroundings are also slightly warmer. The field lines in the side view from both simulations have a very similar pinched shape, which is expected because the field is coupled to the gas in both runs as it is fully ionised at these scales. It is always a challenge to view magnetic field lines in a 2D plane, and Fig. 4 shows a 3D rendering of the magnetic field lines for both simulations, along with density isosurfaces. This view reveals the true topology of the field; a near perfect hourglass in runID, and strong winding inside the second core in runAO. The generation of toroidal field in runAO is expected to eventually lead to the launching of a fast outflow (Machida et al. 2006; Tomida et al. 2013).

thumbnail Fig. 4.

3D visualization of logarithmically spaced density isosurfaces in the inner-most region of the computational domain showing the structure of the second Larson core, in the case of ideal (top) and non-ideal (bottom) MHD. The isosurfaces have been cut half-way in the x-direction. The magnetic field lines are overlayed and have been coloured according to the magnitude of the magnetic field vector. The insets in the lower left corner of each panel show (with the same spatial scale) the central region of the system without the B field for a better view of the morphology. The density and magnetic field colour scales apply to both panels.

Open with DEXTER

3.4. Late evolution

In this section, we look at the subsequent evolution of the IMHD and NIMHD systems. Figure 5 shows density and temperature slices in the two simulations, approximately one month (24 days) after the formation of the second core. The second core in runID is still compact, has reached even higher densities and temperatures in its centre (0.1 g cm−3; 105 K), and appears to have filamentary accretion streams that are associated to the magnetic field topology (see panel b). Its mass is now 9.5 × 10−3 M , with an effective mass accretion rate of ~7 × 10−2 M yr−1.

thumbnail Fig. 5.

Slices of the gas density with velocity vectors in runID (panel a) and runAO panel c, about one month after the formation of the second Larson core. The area shown is the same as in Figs. 3s and 3t. In panel c, the yellow contour marks the disk limit, taken as ρ > 10−6.7 g cm−3, while the black dashed contour delineates the second hydrostatic core with ρ > 10−5 g cm−3. Panels b and d: slices of the gas temperature with magnetic field streamlines overlayed. Panel e: logarithmic map of the magnetic Toomre stability criterion Q mag inside the disk that forms around the second core in runAO. The grey-shaded areas indicate regions in the disk where the epicyclic frequency ω is imaginary and no Q could be computed. The yellow and dashed black contours are the same as in panel c. Panel f : radial profile of the azimuthal velocity for all the cells inside the runAO second core disk. The colours code for the mass contained in a particular region of the plot. A Keplerian velocity profile is overlayed (black solid line).

Open with DEXTER

The small spiral instability in runAO detected in Fig. 3t has developed into a small disk around the second core with two spiral arms. At this point, the second core mass has grown to 7.7 × 10−3 M , for an effective mass accretion rate of ~4 × 10−3 M  yr−1 (the core is delineated by the black dashed contour in Fig. 5c). It has a rotation period of ~22 days. The disk mass is 1.8 × 10−4 M (the disk was defined as the gas with densities in the range 10−6.7 g cm−3 < ρ < 10−5 g cm−3; this is marked by the yellow and dashed black contours). We computed the magnetic Toomre stability criterion Q mag (Kim & Ostriker 2001) for this disk according to (17)

where c s is the gas sound speed, Σ is the disk surface density, and (18)

is the epicyclic frequency of the gas with angular velocity Ω. The surface density was integrated over the height of the disk, while ω, c s, and v A in Eq. (17) actually represent the mass-weighted average values inside a given vertical column (we note that v A ≪ c s because β ≫ 1 at the densities considered). A map of Q mag is displayed in Fig. 5e, revealing that the disk is stable against gravitational contraction. This is suggesting that forming tight binaries from fragmentation inside the second core disk may be difficult, but this is at such an early stage in the protostar’s life that we cannot rule it out with the present result. Indeed, the disk is still rapidly growing in mass (see below), and may become unstable at a later stage. In addition, Fig. 5f shows the distribution of the azimuthal velocity as a function of radius of all the cells inside the disk. Even though the shape of the rotation profile is Keplerian-like, the disk is mostly sub-Keplerian, which is in agreement with the fact that the core is still accreting mass. Finally, it should also be noted that our resolution is insufficient to correctly characterise the viscous dissipation inside the disk and adequately treat the protostellar core accretion shock cooling through the disk. Moreover, following the disk evolution for many orbital periods is computationally prohibitive (see below), and by limiting ourselves to such early epochs, we are not capturing the global disk cooling. These two mechanisms can affect the disk temperature and hence its dynamics and gravitational stability.

The very stringent Courant-Friedrichs-Lewy (CFL; Courant et al. 1967) condition inside the second core (because of the high sound speed) makes it very difficult to integrate for long periods of time after the second core formation. The simulation essentially “freezes” in time, as the timestep in a central region about 0.05 AU in diameter plunges to 10–20 s, which is not tracktable on astrophysical timescales. In addition, the 27 levels of refinement needed to resolve the second core imply that the vast majority of cells lie in a tiny region in the centre of the simulation box, a situation where the CPU domain decomposition along a Hilbert space-filling curve performs poorly. Many processors end up holding no cells in the top AMR levels and spend much of their time waiting for the finer timesteps to complete on the other CPUs. Increasing the number of CPUs beyond 48 did not show convincing boosts in execution speeds, as any gain in processing power gets almost entirely counter-balanced by a heightened communications load.

4. The first and second core accretion shocks

In this final section, we investigate in more detail the accretion flows onto the first and second Larson cores, and more particularly the radiative efficiency of the accretion shocks. Over the years, this subject has been of paramount importance to early evolutionary models of low-mass stars (e.g. Baraffe et al. 2012) as well as planets forming via the core accretion scenario (e.g., Mordasini et al. 2012). Small changes in the fraction of the infalling gas energy that is either absorbed by the core, or radiated away at the accretion shock can yield significant differences in stellar and planetary luminosities and temperatures. However, the lack of accurate models of the accretion shocks which can predict the exact fraction of energy that is accreted or radiatied away in the literature have forced authors to bracket their results using two limiting cases known as “cold” (all energy is radiated away) and “hot” (all energy is absorbed) accretion. Recent numerical studies have suggested that the first Larson core accretion shock tends to be in the super-critical regime, radiating most of the infalling energy away (Commerçon et al. 2011b; Vaytet et al. 2012), while the shock at the second Larson core border is sub-critical, transfering all the energy to the core (Vaytet et al. 2013; Tomida et al. 2013).

These predictions were mainly obtained with onedimensional models of protostellar formation, and we now have the possibility to examine the 3D structure of the accretion flow and the resulting shock efficiency. Because it boasts the more complete microphysics, we consider only the runAO results in this section. In Fig. 6, panels a, b, and c show Hammer projections of the mass accretion rate per unit area , the radiative flux, and the ratio of outgoing radiative flux to incoming gas energy flux F rad/F acc just upstream of the first core accretion shock. Because the hydrostatic core is not spherical, we computed the maps by extracting density, velocity and radiative flux profiles along 64 × 128 different directions, starting from the centre of the second Larson core. The location of the accretion shock in each direction was chosen where the density and velocity gradients are at their maximum. Equations (3) and (7) give us the conservation of total and radiative energy, respectively. Figure 2d revealed that at densities of the first and second Larson cores, the radiative energy is negligible compared to the gas internal energy, and we can thus drop the λv · Er term in (3). In a similar manner, we drop all the terms involving the magnetic field because the plasma β is above 100 for all densities above 10−10 g cm−3 (see Fig. 2c). In a purely conservative form, the gravity term in the right-hand side of Eq. (3) should be included inside the left-hand side divergence. We rewrite it as (19)

thumbnail Fig. 6.

Hammer projections (runAO only) of the mass accretion rate (left column), the radiative flux (middle column), and the ratio of radiative to accretion flux (right column). The first row is for the first Larson core, while the second row shows the second Larson core just after its formation. The third and fourth rows show the second Larson core 24 days after formation and the accretion flow at the edge of the disk around the second hydrostatic core, respectively. The green colours in panels g and i indicate negative values.

Open with DEXTER

Then, because we wish to look at a snapshot of the energy balance at the shock and not an evolution in time, we can assume a stationary state at the core accretion shocks, which means that (19) reduces to  · (ρv Φ) by virtue of (1), and can be inserted directly into the left-hand side divergence. We are now able to write the energy fluxes as(20) (21)

where M enc is the mass enclosed inside the sphere of radius r and ds = r 2 sin θdθdϕ is the line of sight area element. Since we are computing an angular-dependent shock efficiency, we must measure it locally, rather than use a more global definition such as the energy balance scheme recently suggested by Marleau et al. (2017).

Panel a reveals that mass accretion onto the first core is funneled along the dense filaments that we observed in Fig. 3l. These appear as two large, almost circular, hot-spots in the panel a map, centred at longitudes of 75° and −105°. These regions dominate the total mass accretion rate, and illustrate once again that mass accretion is highly anisotropic. We also note that there are no negative values for , meaning that radial velocities are negative everywhere; there are no outflows. In contrast, the radiative flux appears strong in regions of low accretion rate, although this is not a strict correlation. The radiation appears to propagate in directions where it meets low density gas which has a low optical depth. The resulting ratio of radiative to accretion flux in panel c is fascinating. Going against the commonly accepted paradigm that the core endures either cold or hot accretion, the map shows that it can be both at the same time. The accretion flux vastly dominates over its radiative counterpart (by 3 orders of magnitude) in the accretion hot-spots, while the two become comparable elsewhere. Going back to Fig. 2e, we note that the temperature profile of runAO seems to show a temperature discontinuity for some of the gas at a radius of 1 AU, where the temperature jumps from 100 to almost 1000 K. Such a discontinuity is indicative of a radiatively inefficient accretion shock, and the fact that this gas belongs to the equatorial regions (light blue colour) is consistent with the accretion hot-spots we report here. By contrast, there is a small dark blue (polar) region in the temperature profile of Fig. 2e around ~3 AU that exhibits a less pronounced discontinuity, which corresponds to the radiatively efficient polar regions in Fig. 6c.

In the case of the second Larson core (panel d), the mass accretion rate is highest all around the equator, with no predominant hot-spots. This is a result of a higher gas density in the equatorial region just ahead of the shock (see Fig. 3r). On the other hand, the radiative flux is higher in the polar regions where the lower density gas it has to travel through allows it to escape more freely. However, when we compare the accretion and radiative fluxes, even though we see structure dividing equatorial and polar regions, the accretion flux still dominates everywhere, by at least 4 orders of magnitude. This result is thus in agreement with past 1–3D studies (Vaytet et al. 2013; Tomida et al. 2013). The surface integrated mass accretion rate is colossal at 0.28 M yr−1 (also in agreement with Vaytet et al. 2013), and it is difficult to imagine that this will be sustained for very long, as the protostar would finish accreting its entire 1 M envelope in under 4 years. Even though we have only run the simulation for ~1 month after the formation of the second core, we already observe a dramatic drop in mass accretion rate in our final snapshot. Figure 6g displays the structure of the accretion flow onto the second core once the disk seen in Fig. 5 has formed; the strong equatorial accretion has disappeared and some regions of negative accretion (corresponding to positive values of vr , shown in green) have even emerged. The disk acts as a buffer between the infalling material and the protostar; the gas is rotating in almost Keplerian fashion (see Fig. 5) inside the disk, and radial inward motion is governed primarily by viscous transport. The radial velocity – and hence the mass accretion rate – at the protostellar surface is thus considerably reduced. For this final snapshot, we measure a surface integrated mass accretion rate of 0.074 M  yr−1 onto the protostar, but neither this nor the initial mass accretion rate of 0.28 M  yr−1 are a good indication of how fast the core is growing. Indeed, the accretion flow is unsteady and the average mass accretion rate during the first 24 days is only 4 × 10−3 M  yr−1 (as mentioned in the previous section). Conversely, the mass accretion onto the disk is much more stable, with an average value of 2×10−2 M  yr−1. It should, however, be noted that we probably do not have sufficient resolution to adequately resolve instabilities such as the magnetorotational instability (Balbus & Hawley 1991), which generate turbulence and regulate material and angular momentum transport inside the disk. Nevertheless, even if the mass accretion flow is unsteady, the ratio of infalling (kinetic and gravitational) to outgoing (radiative) energy is actually very stable; the accretion shock is radiatively inefficient throughout the early evolution of the protostar (panels f and i). The accretion energy flux also dominates over the radiation flux at the edge of the disk (panel 1).

We emphasise here that these results only apply to the very early stages of the protostar’s evolution, and cannot be assumed to hold for the remainder of the main accretion phase. They merely suggest that the second core accretion shock is initially radiatively inefficient, and reveal that it is possible to have both hot and cold accretion at the same time over the surface of the first core. We are reporting on the structure of the accretion flow at the birth of the protostar, and we do not know if this accretion arrangement can be applied to protostellar evolution models. We simply hint that the picture may not be either fully hot or cold; both regimes could be operating at the same time over the surface of the hydrostatic cores.

5. Comparison with previous works

In this section, we compare the present study with previous articles that report on simulations of protostellar formation. For the sake of brevity, we limit ourselves to 3D non-ideal MHD simulations that have reached the second Larson core stage.

The first 3D models including ohmic diffusion were performed by Machida et al. (2006) using a nested-grid MHD code. The main difference between their models and our runs is that they use a barotropic equation of state, while we include radiative transfer via the FLD. They also lack ambipolar diffusion. Nevertheless, they already report a strong increase in plasma β and angular momentum when number densities exceed 1014 cm−3 in the resistive run compared to using ideal MHD. In the past five years, Tomida et al. (2013, 2015) and Tsukamoto et al. (2015a) performed simulations including radiative transfer via the FLD, as well as non-ideal MHD with ohmic diffusion and ambipolar diffusion. The most recent work by Wurster et al. (2018) includes radiative transfer and the three non-ideal MHD effects.

Table 2 shows the properties of the first and second cores formed in our simulations. Overall, our results are qualitatively similar to those reported in the recent literature within a factor of a few (since we do not use the same definition criteria for the first and second cores, we expect to have small differences). For instance, Tomida et al. (2013) reported second core mass of 2 × 10−2 M one year after its formation. Assuming the system settles on timescales much shorter than a year after formation (i.e. about a month, as observed in our simulation), this yields an average mass accretion rate of 2 × 10−2 M yr−1, which is five times our measured rate of 4 × 10−3 M yr−1. However, Tomida et al. (2013) define their protostellar core as a pressure-supported body that would also include the small disk in our simulation (see Fig. E.1). Considering the disk as part of the second core means the second core mass accretion is now the flux at the disk border, which stands at 2 × 10−2 M yr−1 (cf. Sect. 4) and is now entirely consistent with Tomida et al. (2013). The second core mass and size we derive are also roughly consistent with the results of Wurster et al. (2018) six months after the stellar core formation, who find masses of 1.5 × 10−2 M in IMHD and 3.4 × 10−3 M in NIMHD, as well as a radius of 0.013 × 10−2 AU in both cases. We note that they use a similar criterion as ours for the second core definition, but with a density threshold a factor of ten higher. In addition, Tsukamoto et al. (2015a) found plasma beta within the first cores β > 104 in NIMHD and β ~ 10 in IMHD, which is fully consistant with Fig. 2c.

Besides this qualitative agreement, there are some discrepancies in the structure of the collapsing core, as well as in the first core lifetime. First, Tomida et al. (2013, 2015) found that outflows and disks form early, even prior to the second collapse (with ohmic and ambipolar diffusion). The outflows reported in Tomida et al. have a relatively small extent 170 AU maximum at the end of the first core phase. Second they observed longer first core lifetimes and the latter increases when non-ideal MHD effects are included, whereas we find the opposite. In our models, we attribute this increase in the first core lifetime with IMHD to the development of interchange instabilities which heat up and bloat the first core (see Sect. 3.1). Interchange instabilities are reported in Tomida et al. (2015) but do not affect the first core in IMHD as in ours. We think that these differences originate from the initial conditions. While we use uniform initial density profile, Tomida et al. used Bonnor-Ebert profile which is close to equilibrium. The time spent to form the first core is much longer when the initial core mass is close to the Bonnor-Ebert mass (see Vaytet & Haugbølle 2017, Fig. 7 therein). As previously mentioned, the accretion rate is a factor ~5 higher in our models than in Tomida’s, so that the first core evolves much quicker and the dynamic is more violent, leading to powerful magnetic interchange instability. The absence of outflows and large disks in our results is also consistent with the differences excepted between models using either a uniform or a Bonnor–Ebert density profile (Machida et al. 2014). In addition, Tsukamoto et al. (2015a) used uniform initial density and found that the protostellar disk forms after the second core in their NIMHD models. Wurster et al. (2018) also report outflows at first core scales in NIMHD models using similar initial conditions as ours. However, they observe that outflows become broader and slower as the cosmic ray ionisation rate is reduced. The minimum ionisation rate they explore is 10−16 s−1 while we use 10−17 s−1. Whether outflows launching at the first core scale depends on the cosmic ray ionsitation rate remains to be studied in detail. Clearly, the effect of the initial conditions, as well as the effect of the chemical set up used to estimate the MHD resistivities, has to be investigated in the near future to truly compare results.

Table 2.

Properties of the first and second Larson cores extracted about one month after the birth of the second core.

6. Conclusions

We have performed two 3D simulations of the gravitational collapse of a dense sphere of molecular cloud gas. Both runs include the following physics: hydrodynamics, radiative transfer, self-gravity, a non-ideal gas equation of state, and magnetic fields. In the second run, the effects of ambipolar and ohmic diffusion were included in the MHD equations, and their impact on the simulation results were assessed through comparisons with the ideal MHD model. The magnetic diffusion creates a barrier which prevents amplification of the magnetic field beyond 0.1 G in the first Larson core, with many consequences for the structure and evolution of the system. In the IMHD simulation, the magnetic field dominates the energy budget everywhere inside and around the first core, spawning interchange instabilities that create bubble-like ejections, as well as driving a low-velocity outflow above and below the equatorial plane of the system. A strong magnetic field also implies a heightened magnetic braking, removing essentially all angular momentum from the second Larson core.

When ambipolar and ohmic diffusion are present, the first and second cores become genuinely thermally supported and have a large amount of rotation. This leads to the formation of a small Keplerian-like gravitationally stable disk around the second core, and rolls the magnetic field lines into a toroidal topology which is expected to propel an outflow at the second core level. Due to stringent CFL limitations, it was, however, not possible for us to follow the evolution of the system long enough to observe the launch. We were also neither able to study the formation of a protoplanetary disk and a low-velocity outflow (Gerin et al. 2017) around the first Larson core because the simulation essentially “froze” in time when the second core was formed. Future plans involve replacing the second core with a sink particle, allowing for much longer time integrations. The stark contrast between the ideal and NIMHD simulations proves that magnetic diffusion is of crucial importance to star-formation; not only does it enable the formation of disks in which planets will eventually form (Masson et al. 2016), it also shapes the protostar itself by preventing angular momentum loss and restoring thermal pressure support.

The use of idealised isolated initial conditions has been challenged by recent studies which claim that accretion processes in star formation are vastly influenced by the environment around the protostellar system (Kuffmeier et al. 2017). And while this may indeed be relevant at the first Larson core scale, we postulate that the dynamics at the second Larson core level are so disconnected, both in terms of spatial scales and evolutionary timescales, from the material 100 AU away, that the impact of large-scale turbulence would be negligible. Nevertheless, we are currently investigating the robustness of our results across different initial conditions, varying the parent cloud mass, changing the magnetic field strength and orientation, and introducing turbulence in the initial velocity field. Another shortcoming of the model presented in this paper is the lack of Hall effect in the MHD solver. Believed to be prominent in protoplanetary disks, the Hall effect has attracted much attention of late (e.g. Lesur et al. 2014; Tsukamoto et al. 2015b, 2017; Wurster et al. 2016), and is considered to play a major role in angular momentum transport both inside the disk and in the protostellar envelope. We are in the process of implementing the Hall effect in our version of RAMSES. Last but not least, large uncertainties remain in the models used to estimate the resistivity coefficients because of poor constraints on the dust size properties (charge, size distribution) and on the chemistry at play in the high density and temperature regions of protostellar collapse. As a result, it is currently not clear which non-ideal effects dominate in the different parts of the collapsing cloud, particularly for the Hall and ambipolar resistivities that strongly depend on the local physical and chemical conditions. Further work is required to better estimates of the non-ideal resitivities, which would in turn allow a more robust assessment of their impact on the star, disk, and planet formation process.

Acknowledgments

We are indebted to the anonymous referee for his/her insightful comments that have vastly improved the solidity of our study, with no stones left unturned. We also thank Troels Haugbølle for very useful discussions during the writing of this paper. NV gratefully acknowledges support from the European Commission through the Horizon 2020 Marie Skłodowska-Curie Actions Individual Fellowship 2014 programme (Grant Agreement no. 659706). The research leading to these results has also received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 247060). We acknowledge financial support from “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, CEA and CNES, France. This work was granted access to the HPC resources of CINES (Occigen) under the allocation 2016-047247 made by GENCI. We also made use of the astrophysics HPC facility at the University of Copenhagen, which is supported by a research grant (VKR023406) from Villum Fonden. In addition, we thank the Service d’Astrophysique, IRFU, CEA Saclay, and the Laboratoire Astrophysique Instrumentation Modélisation, France, for granting us access to the supercomputer IRFUCOAST where the groundwork with many test calculations were performed. All the figures were created using the OSIRIS8 visualization package for RAMSES, except Fig. 4 which was rendered with the PARAVIEW9 software.

References

  1. André, P., Di Francesco, J., Ward-Thompson, D., et al., 2014, in Protostars and Planets VI, eds. H.Beuther, R. F.Klessen, C. P.Dullemond, & T.Henning (Tucson: The University of Arizona Press), 27 [Google Scholar]
  2. Badnell, N. R., Bautista, M. A., Butler, K., et al., 2005, MNRAS, 360, 458 [NASA ADS] [CrossRef] [Google Scholar]
  3. Balbus, S. A., & Hawley, J. F.1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
  4. Baraffe, I., Vorobyov, E., & Chabrier, G.2012, ApJ, 756, 118 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bate, M. R.1998, ApJ, 508, L95 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bate, M. R.2010, MNRAS, 404, L79 [NASA ADS] [Google Scholar]
  7. Bate, M. R.2011, MNRAS, 417, 2036 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bate, M. R., Tricco, T. S., & Price, D. J.2014, MNRAS, 437, 77 [NASA ADS] [CrossRef] [Google Scholar]
  9. Béthune, W., Lesur, G., & Ferreira, J.2016, A&A, 589, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R.2008, A&A, 482, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R.2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Commerçon, B., Hennebelle, P., & Henning, T.2011a, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
  13. Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G.2011b, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Commerçon, B., Debout, V., & Teyssier, R.2014, A&A, 563, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Courant, R., Friedrichs, K., & Lewy, H.1967, IBM J. Res. Dev., 11, 215 [NASA ADS] [CrossRef] [Google Scholar]
  16. Crutcher, R. M.2012, ARA&A, 50, 29 [NASA ADS] [CrossRef] [Google Scholar]
  17. Dapp, W. B., & Basu, S.2012, A&A, 541, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Duffin, D. F., & Pudritz, R. E.2008, MNRAS, 391, 1659 [NASA ADS] [CrossRef] [Google Scholar]
  19. Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S.2011, ApJ, 731, 62 [NASA ADS] [CrossRef] [Google Scholar]
  20. Ferguson, J. W., Alexander, D. R., Allard, F., et al., 2005, A&AS, 623, 585 [Google Scholar]
  21. Fromang, S., Hennebelle, P., & Teyssier, R.2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Gerin, M., Pety, J., Commerçon, B., et al. 2017, A&A, 606, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. González, M., Vaytet, N., Commerçon, B., & Masson, J., 2015, A&A, 578, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Hennebelle, P., & Charbonnel, C.2013, EAS Pub. Ser., 62 [Google Scholar]
  25. Hennebelle, P., & Ciardi, A.2009, A&A, 506, L29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Hennebelle, P., & Fromang, S.2008, A&A, 477, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hennebelle, P., & Teyssier, R.2008, A&A, 477, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P.2016, ApJ, 830, L8 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hull, C. L. H., Plambeck, R. L., Bolatto, A. D., et al. 2013, ApJ, 768, 159 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hull, C. L. H., Mocz, P., Burkhart, B., et al. 2017, ApJ, 842, L9 [NASA ADS] [CrossRef] [Google Scholar]
  31. Joos, M., Hennebelle, P., & Ciardi, A.2012, A&A, 543, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S.2013, A&A, 554, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Kim, W.-T., & Ostriker, E. C.2001, ApJ, 559, 70 [NASA ADS] [CrossRef] [Google Scholar]
  34. Krasnopolsky, R., Li, Z.-Y., & Shang, H.2010, ApJ, 716, 1541 [NASA ADS] [CrossRef] [Google Scholar]
  35. Krasnopolsky, R., Li, Z.-Y., Shang, H., & Zhao, B.2012, ApJ, 757, 77 [NASA ADS] [CrossRef] [Google Scholar]
  36. Krumholz, M. R., Klein, R. I., & McKee, C. F.2012, ApJ, 754, 71 [NASA ADS] [CrossRef] [Google Scholar]
  37. Krumholz, M. R., Crutcher, R. M., & Hull, C. L. H.2013, ApJ, 767, L11 [NASA ADS] [CrossRef] [Google Scholar]
  38. Kuffmeier, M., Haugbølle, T., & Nordlund, Å.2017, ApJ, 846, 7 [NASA ADS] [CrossRef] [Google Scholar]
  39. Larson, R. B.1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
  40. Lazarian, A.2013, in Numerical Modeling of Space Plasma Flows, 7th Int. Conf. (ASTRONUM2012), ASP Conf. Ser., 474, 15 [NASA ADS] [Google Scholar]
  41. Leão, M. R. M., de Gouveia Dal Pino, E. M., Santos-Lima, R., & Lazarian, A.2013, ApJ, 777, 46 [NASA ADS] [CrossRef] [Google Scholar]
  42. Lesur, G., Kunz, M. W., & Fromang, S.2014, A&A, 566, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. LevermoreC. D.1984, JQSRT, 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
  44. Li, Z.-Y., Krasnopolsky, R., & Shang, H.2011, ApJ, 738, 180 [NASA ADS] [CrossRef] [Google Scholar]
  45. Li, Z.-Y., Krasnopolsky, R., & Shang, H.2013, ApJ, 774, 82 [NASA ADS] [CrossRef] [Google Scholar]
  46. Li, Z.-Y., Krasnopolsky, R., Shang, H., & Zhao, B.2014, ApJ, 793, 130 [NASA ADS] [CrossRef] [Google Scholar]
  47. Machida, M. N., & Matsumoto, T.2011, MNRAS, 413, 2767 [NASA ADS] [CrossRef] [Google Scholar]
  48. Machida, M. N., Inutsuka, S.-i, & Matsumoto, T.2006, ApJ, 647, L151 [NASA ADS] [CrossRef] [Google Scholar]
  49. Machida, M. N., Inutsuka, S.-i, & Matsumoto, T.2007, ApJ, 670, 1198 [NASA ADS] [CrossRef] [Google Scholar]
  50. Machida, M. N., Tomisaka, K., Matsumoto, T., & Inutsuka, S.-i.2008, ApJ, 677, 327 [NASA ADS] [CrossRef] [Google Scholar]
  51. Machida, M. N., Inutsuka, S., & Matsumoto, T.2014, MNRAS, 438, 227 [NASA ADS] [CrossRef] [Google Scholar]
  52. Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Marleau, G.-D., Klahr, H., Kuiper, R., & Mordasini, C.2017, ApJ, 836, 221 [NASA ADS] [CrossRef] [Google Scholar]
  54. Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G.2012, ApJS, 201, 24 [NASA ADS] [CrossRef] [Google Scholar]
  55. Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B., 2016, A&A, 587, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Masunaga, H., & Inutsuka, S.-i2000, ApJ, 531, 350 [NASA ADS] [CrossRef] [Google Scholar]
  57. Matsumoto, T., & Tomisaka, K.2004, ApJ, 616, 266 [NASA ADS] [CrossRef] [Google Scholar]
  58. McKee, C. F., & Ostriker, E. C.2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
  59. Mellon, R. R., & Li, Z.-Y.2008, ApJ, 681, 1356 [NASA ADS] [CrossRef] [Google Scholar]
  60. Mellon, R. R., & Li, Z.-Y.2009, ApJ, 698, 922 [NASA ADS] [CrossRef] [Google Scholar]
  61. Minerbo, G. N.1978, J. Quant. Spec. Rad. Transf., 20, 541 [NASA ADS] [CrossRef] [Google Scholar]
  62. Mocz, P., Burkhart, B., Hernquist, L., McKee, C. F., & Springel, V.2017, ApJ, 838, 40 [NASA ADS] [CrossRef] [Google Scholar]
  63. Mordasini, C., Alibert, Y., Klahr, H., & Henning, T.2012, A&A, 547, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Mouschovias, T. C., & Spitzer, Jr.L.1976, ApJ, 210, 326 [NASA ADS] [CrossRef] [Google Scholar]
  65. Planck Collaboration Int. XXXV. 2016, A&A, 586, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Saigo, K., & Tomisaka, K.2006, ApJ, 645, 381 [NASA ADS] [CrossRef] [Google Scholar]
  67. Saigo, K., Tomisaka, K., & Matsumoto, T.2008, ApJ, 674, 997 [NASA ADS] [CrossRef] [Google Scholar]
  68. Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A.2012, ApJ, 747, 21 [NASA ADS] [CrossRef] [Google Scholar]
  69. Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A.2013, MNRAS, 429, 3371 [NASA ADS] [CrossRef] [Google Scholar]
  70. Saumon, D., Chabrier, G., & van HornH. M.1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  71. Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S.2013, MNRAS, 432, 3320 [NASA ADS] [CrossRef] [Google Scholar]
  72. Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S.2015, MNRAS, 446, 2776 [NASA ADS] [CrossRef] [Google Scholar]
  73. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E., 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Spruit, H. C., Stehle, R., & Papaloizou, J. C. B.1995, MNRAS, 275, 1223 [NASA ADS] [CrossRef] [Google Scholar]
  75. Stamatellos, D., Whitworth, A. P., Bisbas, T., & Goodwin, S.2007, A&A, 475, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Teyssier, R.2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Tomida, K., Machida, M. N., Saigo, K., Tomisaka, K., & Matsumoto, T.2010, ApJ, 725, L239 [NASA ADS] [CrossRef] [Google Scholar]
  78. Tomida, K., Tomisaka, K., Matsumoto, et al. 2013, ApJ, 763, 6 [NASA ADS] [CrossRef] [Google Scholar]
  79. Tomida, K., Okuzumi, S., & Machida, M. N.2015, ApJ, 801, 117 [NASA ADS] [CrossRef] [Google Scholar]
  80. Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S.-i.2015a, MNRAS, 452, 278 [NASA ADS] [CrossRef] [Google Scholar]
  81. Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S.-i.2015b, ApJ, 810, L26 [NASA ADS] [CrossRef] [Google Scholar]
  82. Tsukamoto, Y., Okuzumi, S., Iwasaki, K., Machida, M. N., & Inutsuka, S.-i.2017, PASJ, 69, 95 [NASA ADS] [CrossRef] [Google Scholar]
  83. Vaytet, N., Audit, E., Chabrier, G., Commerçon, B., & Masson, J.2012, A&A, 543, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Vaytet, N., Tomida, K., & Chabrier, G.2014, A&A, 563, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Vaytet, N., & Haugbølle, T.2017, A&A, 598, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Whitehouse, S. C., & Bate, M. R.2006, MNRAS, 367, 32 [NASA ADS] [CrossRef] [Google Scholar]
  88. Wurster, J., Price, D. J., & Bate, M. R.2016, MNRAS, 457, 1037 [NASA ADS] [CrossRef] [Google Scholar]
  89. Wurster, J., Bate, M. R., & Price, D. J.2018, MNRAS, 475, 1859 [NASA ADS] [CrossRef] [Google Scholar]
  90. Zhao, B., Li, Z.-Y., Nakamura, F., Krasnopolsky, R., & Shang, H.2011, ApJ, 742, 10 [NASA ADS] [CrossRef] [Google Scholar]

1

It is not clear why Krasnopolsky et al. (2010) and Li et al. (2011) were not able to form rotationally supported disks in their calculations. Possible reasons include that their models were only 2D, and did not incorporate self-gravity, although this has never been confirmed.

2

We note that Tomida et al. (2015) did not quite follow the evolution of the collapsing system all the way up to the formation of the second core.

3

The high cost for the non-ideal MHD simulation does not originate from a computationally expensive magnetic diffusion module, but comes primarily from a highly reduced integration timestep between the first and second collapse stages, as ambipolar and ohmic resistivities increase inside the first hydrostatic core (see Eq. 10).

4

It is actually closer to 5/3 for 10−13 < ρ < 10−12 g cm−3 (see Vaytet et al. 2014).

5

The maximum resolution of 0.15 AU in Masson et al. (2016) verifies this condition.

6

This may seem a little artificial but it in fact reproduces very well the density structures seen in simulations with more realistic turbulent initial conditions (Commercon et al., in prep.).

7

This was once again already observed in the simulations of Masson et al. (2016, see their Fig. 6).

10

In RAMSES, outputs are only written when a coarse step has been completed, and it is often not trivial to write snapshots at exactly the same simulation time in two different simulations.

11

Many of these errors are also due to the fact that the snapshots from the two simulations are not written at exactly the same simulation time.

Appendix A: Minimum optical depth per cell

In this section, we describe a scheme we devised to aid the convergence of the implicit radiative transfer solver. When the gas is optically thin, it is not crucially important for the heating and cooling mechanisms whether the optically depth inside a given cell is 10−8 or 10−4, as long as it is much less than unity. However, very low optical depths typically require many iterations for a time-implicit radiation solver to converge. We artificially limited the optical depth per cell to a minimum value of 10−4, by setting the mean Rosseland opacity to (A.1)

The flooring occurs in the large (low AMR level) low density cells, in the outer regions of the protostellar envelope. Figure A.1a shows the fraction of cells where the optical depth is being limited, with respect to the total number of cells in the simulation, as a function of time (red solid line). The black dashed line shows the evolution of the density at the centre of the collapsing cloud (i.e. inside the densest cell) with time. We see that while the fraction of cells with limited κRρΔx is large (~80%) at early times, it drops below 0.1 when the first Larson core forms (t ~ 28 kyr and ρ ~ 10−10 g cm−3). In panel b of Fig. A.1, we show the total number of cells per AMR level (grey histogram), for a snapshot at a time of 28.180 kyr. The red histogram shows the cells where the optical depth is being limited. We can see that the floor is operating only in the outer layers of the collapsing system, from AMR level 6 to 11, and will not impact the properties of the first and second Larson cores. In the ideal MHD simulation presented in the main part of this paper (up until a simulation time of 28.180 kyr), the total number of iterations is reduced by 25%, and the computational time reduced by 20%.

thumbnail Fig. A.1.

Panel a: fraction of cells inside the mesh where the optical depth is being limited as a function of time (red solid line). The dashed black line shows the density at the centre of the system as a function of time. Panel b: number of cells in each level (grey) and the number of cells where the optical depth floor is operating (red), at a time of 28.180 kyr, when the first Larson core is formed. Panel c: relative difference in 2D histograms of gas temperature as a function of density for all the cells in a simulation with optical depth limitation and a second simulation without, at t = 28.180 kyr. The colour scale gives a measure of ℛ = Nlimited/Nnot limited − 1, where N is the number of cells binned inside a (ρ, T) pixel for the two different simulations. Panel d: same as panel c but in the case of the Rosseland mean opacity as a function of density.

Open with DEXTER

To validate the optical depth flooring scheme, we show in panel c the temperature/density distribution of all the cells in the mesh for two simulations. The first has the optical depth limitation switched on, while it is turned off in the second. The coloured contours show the relative difference ℛ between the two simulations, for each (ρ, T) pixel in the plot. It is defined as ℛ = Nlimited/Nnot limited − 1, where N is the number of cells binned inside a (ρ, T) pixel. A red area indicates that there are more cells from the simulation with the limitation scheme than from the run without the κRρΔx floor in that particular region of the plot, and vice versa for blue areas. The differences are expected to be the largest at low densities. However, in this isothermal phase of collapse, all the gas has a constant temperature of 10 K and the optical depth limiting scheme has basically no impact on the results. Small differences, of the order of 1%, are visible at higher densities, but these mostly originate from the fact that the two simulation outputs have been written at slightly different times10. Finally, in panel d we show the Rosseland mean opacity as a function of density, using the same convention as in panel c. It is once again obvious that the limiter is only active in the outer layers of the infalling envelope, where the flow is still isothermal. The limited opacities show a stripy pattern which is due to the refinement of cells. We conclude that the optical depth limitation scheme does not appear to affect the thermodynamics of the system as it operates only in the isothermal stage of the collapse.

Appendix B: The timestep limitation scheme

One of the difficulties when working with diffusion processes on a mesh based framework is that the timestep criterion for numerical stability usually scales with the square of the mesh size ∆x. This is indeed the case for ambipolar and ohmic diffusion, and is made worse by the fact that as densities increase, not only does ∆x decrease but the resistivities can also increase by several orders of magnitude (see Fig. 5 in Marchand et al. 2016). This double effect (see Eq. 10) causes the timestep Δt to fall abruptly after the first Larson core is formed, and would require millions of timesteps to reach the second Larson core formation, making the problem non-tracktable. In the same spirit as limiting the optical depth per cell in the previous section, where we found that as long as the optical depth in a cell is much less than unity its exact value does not matter for our purposes, we postulate that as long as a strong magnetic diffusion is operating, the precise amount will not affect our results in a crucial way.

As mentioned in Sect. 2, the method we have chosen to try and prevent the MHD timestep from reaching prohibitively low values is to artificially limit the value of Δt to a fraction ξ of the ideal MHD timestep Δt ID. In practice, we found that setting the lower limit to ξ = 0.1 was a good compromise between speedup and accuracy of results. We emphasise that we have no physical justification for the value of 0.1, it was simply chosen after months of testing. To ensure consistency between the imposed value of Δt and the magnetic diffusion, one has to artificially lower the resistivities in the cells which would have Δt O,A < ξΔt ID. The resistivities are thus overwritten with . Note here that the factor of 0.1 in the numerator of the fraction on the right-hand side is different from the ξ = 0.1; it corresponds to the CFL-like factor that is used to compute the diffusion timestep, taken as a tenth of the time it would take for all the magnetic field inside the cell to diffuse.

Validation of this acceleration scheme is explicited in Fig. B.1. Panel a shows the fraction of cells inside the computational domain where the resistivities are being modified, as a function of time. The black dashed line represents the evolution of the central density, and we can see that as it reaches values characteristic of the first Larson core (~10−12 to 10−10 g cm−3), the numbers of cells where ΔtO,A is floored begin to increase. However, these fractions remain small throughout the simulation, peaking at 25% for the ambipolar diffusion (red) and 10% for the ohmic diffusion (blue). In addition, the flooring is only important during a transition phase between the formation of the first and second Larson cores, since after having increased with density, the resistivities begin to fall again once temperatures increase beyond ~1500 K where the dust grains evaporate (see Fig. 2b and Marchand et al. 2016). This is indeed reflected by the sharp fall in fractions (blue and red lines) as the density abruptly increases past 10−8 g cm−3. A histogram showing the number of cells affected by the Δt flooring for each AMR level, taken at a time of 28.2 kr where the fractions in panel a reach their maxima, is displayed in panel b. As the flooring operates only in the densest parts of the system, only the highest AMR levels are affected.

thumbnail Fig. B.1.

Panel a: fraction of cells inside the mesh where ηA (red) and ηO(blue) are being modified to prevent the MHD timestep from becoming too small, as a function of time. The dashed black line shows the evolution of the central density. Panel b: number of cells per AMR level (grey) and the number of cells where the ambipolar (red) and ohmic (blue hatched) diffusion timestep floor is operating, at a time of 28.2 kyr. Panel c: relative difference in 2D histograms of magnetic field strength as a function of density for all the cells in a simulation with Δt flooring and a second simulation without, at t = 28.2 kyr. The colour scale is analogous to that of Fig. A.1. Panel d: same as panel c but in the case of the ambipolar (red-blue) and ohmic (green-brown) resistivities.

Open with DEXTER

The resistivities affect primarily the magnetic field, and we show in panel c a distribution of the magnetic field as a function of density in every cell in two different simulations. The first has the acceleration scheme switched on, while the other is without. Because of the prohibitively small values of ∆t in the simulation without timestep acceleration, we ran both calculation with a resolution of only 12 points per Jeans length. As in the previous section, the coloured contours show the relative difference ℛ between the two simulations, for each (ρ, B) pixel in the plot. It is defined as ℛ = Naccel/Nno accel − 1, where N is the number of cells binned inside a (ρ, B) pixel. A red area indicates that there are more cells from the simulation with the acceleration scheme than from the run without the Δt floor in that particular region of the plot, and vice-versa for blue areas. As expected, the timestep limitation scheme changes the magnetic diffusion plateau at high densities (ρ > 10−13 g cm−3), but only in a very minor way. The accelerated simulation still displays a strong magnetic diffusion barrier around 0.1 G, and the values of B differ by 5% or less in the rest of the computation box, compared to the run with the correct Δt11. This, we argue, is the justification for using the acceleration scheme; the magnetic diffusion is still operating, and still dominates over any numerical diffusion. The diffusion is crucial to limiting the magnetic braking and the accumulation of magnetic flux, and this is still achieved in the accelerated run. In the last panel d, we show for informative purposes the values of the resistivities as a function of density, using the same colour convention as in panel c. The differences below ρ ~ 10−13 g cm−3 are once again due to a different simulation time output, and the resistivities are only modified by the acceleration scheme at high densities. Even though the resistivities can be modified by more than an order of magnitude, as long as they are high enough, the exact values of ηA,O do not seem to be important in the scope of our simulations.

It is of course difficult to predict the impact of such an acceleration scheme on simulation results without running the full (non Δt-limited) simulation first, as it is potentially highly problem-dependent. Even though we tested the method across a range of initial conditions (different parent cloud masses, initial magnetization, temperature, rotation) and it always gave excellent results, we limited ourselves to the problem of a gravitationally collapsing magnetised body, and we must advise caution when using it for a different kind of set-up.

Appendix C: Resolution study

In star formation studies, the refinement criterion when using an AMR mesh is usually based on the Jeans length. In other words, the Jeans length needs to be adequately sampled to properly resolve the system dynamics. There has been some debate as to how many cells per Jeans length are actually necessary, and authors commonly use 10–16 cells per Jeans length (e.g. Commerçon et al. 2011a; Krumholz et al. 2012). Vaytet & Haugbølle (2017) recently showed, using 1D simulations, that resolution can affect the thermodynamics of collapsing dense clouds, because of poor sampling of the optical depth which limits radiation cooling and causes spurious heating inside the first Larson core. If the optical depth within a cell is too large (typically >100), Vaytet & Haugbølle (2017) found that the radiative flux points inward the first core, which creates a spurious bump in the temperature profile. This numerical effect happens when the numerical resolution is too low, and Vaytet & Haugbølle (2017) showed empirically that limiting the optical depth within a cell to a few tens is enough to prevent it. We performed a resolution study to show that this effect can be also prevented in 3D simulations and to ensure it was not affecting the evolution of the protostellar system.

To determine the resolution requirements of our set-up, we ran a simulation with a lower resolution of 16 cells per Jeans length and compare it to our fiducial resolution of 32 cells per Jeans length. The results are shown in Fig C.1. The red contours are for the low-resolution run, while the blue contours are for the calculation with 32 cells per Jeans length. The dashed lines show the evolution of the densest cell in the system, and can be compared to the 1D results of Vaytet & Haugbølle (2017). In the low-resolution run, we actually observe a “turn off” in the first adiabatic phase, at densities ~10−9 g cm−3, while the high-resolution path continues along the same adiabatic track. This departure from adiabaticity actually looks identical to the phenomenon observed by Vaytet & Haugbølle (2017). We also note that the gas is hotter in the low-resolution simulation. It is obvious here that 16 cells per Jeans length is not enough to properly describe the physical processes at work. In fact, we can also see just at the top right end of the high-resolution track a small “kink” in the curve, suggesting that even 32 cells might not be enough for fully converged results. However, the simulation with ambipolar and ohmic diffusion would have been too expensive to run with anything more than NJeans = 32, and we determined that the consequences of such a small kink would only be minimal.

thumbnail Fig. C.1.

Temperature as a function of density, for every cell in the computational domain (ideal MHD case). The simulation using NJeans = 32 cells per Jeans length is represented by the blue area, while the red region is for the run with only 16 cells per Jeans length. Each data set is delineated by a solid contour line which outlines the data distributions. The two snapshots were taken at similar evolution times, chosen to be just after the NJeans = 16 run has departed from its initial adiabatic track. The dashed lines represent the time evolution of the central (densest) cell inside the mesh (these tracks continue beyond the time of the snapshots to provide a wider context). The black arrow indicates the place where the low-resolution track departs from its original adiabat.

Open with DEXTER

From this short resolution study, we see that a refinement criterion solely based on the local Jeans length is not adapted to describe the adiabatic evolution of a hydrostatic core in collapse calculations. A dedicated study of the necessary numerical resolution within the different components of a collapsing core (envelope, disk, hydrostatic cores) is clearly needed and should be the focus of future work.

Appendix D: Regions of active ambipolar and ohmic diffusion

We compute here dimensionless numbers which reveal the regions on active ambipolar and ohmic diffusion in our system. Following Tomida et al. (2015) and Masson et al. (2016), we define the ambipolar and ohmic Reynolds (or sometimes called Elsasser) numbers as (D.1)

where V is the magnitude of the gas velocity vector, and L represents the typical scale of the system, which we take as the distance from the current cell to the centre of the protostar. Figure D.1 shows a map of the logarithm of ε A (coloured contours) in the vicinity of the first Larson core (side view) with the magnetic field lines overlayed (light grey). The regions where ε A ≲ 1 (white and red) have strong ambipolar diffusive effects that modify the magnetic field topology. Indeed, the equatorial pinching of field lines, which is evident in the IMHD run (see Fig. 3m), is reduced when εA < 5 (inside 30 AU), and eventually disappears when εA < 1 (inside 10 AU).

thumbnail Fig. D.1.

Map of the ambipolar (filled blue/red contours) and ohmic (green) Reynolds numbers close to the first Larson core. The light grey lines represent the magnetic field.

Open with DEXTER

In contrast, the green region in Fig. D.1 represents areas where ohmic diffusion is active (εO < 5); it is much smaller because the ohmic resistivities peak at higher densities than their ambipolar counterpart (see Fig. B.1d). This reveals that the straightening of the field lines observed in Sect. 3.3.1 and Fig. 3f, m is due to the effects of ambipolar diffusion.

Appendix E: Definitions of the first and second proto-stellar cores

In this section, we take a look at two different definitions of the first and second Larson cores and how they may affect core morphologies, masses and radii. The cores are often referred to as “hydrostatic cores” in the literature, as they are supposedly (for the most part) in hydrostatic equilibrium. Computing the condition for hydrostatic equilibrium is often expensive in a 3D system, as pressure gradients have to be calculated in all directions, and authors have often favoured simpler criteria such as vanishing radial velocities or thermal-to-kinetic pressure equilibrium. Choosing one definition over the other can sometimes result in large differences in the extent of the core, and consequently the mass that is attributed to it. In Fig. E.1, we compare two different definitions for the proto stellar cores. These are:

  1. Thermal pressure exceeds ram pressure:

  2. Density exceeds a chosen threshold: ρ > ρ core

thumbnail Fig. E.1.

Maps and contours showing the morphologies of the cores using two different definitions. The coloured maps show the ratio of thermal to infalling ram (kinetic) pressure, while the black solid contour defines the region where the gas density exceeds density thresholds of ρcore = 10−10 g cm−3 for the first Larson core and ρcore = 10−5 g cm−3 for the second Larson core. The panels are: (a) runID first core, (b) runAO first core, (c) runID second core, (d) runID second core. We note the difference in spatial scales between panels a and b.

Open with DEXTER

The first condition characterises a thermally supported body, and is equivalent (within a factor of γ) to the definition in Tomida et al. (2010). The second definition is the one we have used throughout this paper. We chose ρcore = 10−10 g cm−3 for the first Larson core and ρcore = 10−5 g cm−3 for the second Larson core.

The left column of Fig. E.1 shows the first and second cores in runID, while the right one is for runAO. For the first core in runID (panel a), it is clear that definitions 1 is affected by the interchange instability which creates a large region of thermally supported gas. The resulting morphology is not what is usually associated with a hydrostatic core, with loops presumably connected to the magnetic field. On the other hand, definition 2 yields a close-to-spherical body. In contrast, both definitions produce similar results for the runAO first core (panel b), where the core is an unbroken/continuous body, flattened on its north and south faces by the heavy accretion streams that slam onto its surface. In the case of the second core, the situation is reversed. Both definitions agree for runID (panel c) but large discrepancies emerge for runAO (panel d). Indeed, the small disk around the second core is also pressure-supported (see Sect. 3.4) and definition 1 considers it to be part of the proto-stellar core, while definition 2 selects only a small spheroidal core, excluding the disk around it.

All Tables

Table 1.

3D numerical studies of the formation of the second Larson core.

Table 2.

Properties of the first and second Larson cores extracted about one month after the birth of the second core.

All Figures

thumbnail Fig. 1.

Density (panel a), temperature (panel b) and magnetic field strength (panel c) as a function of time, for the densest cell in the system. The red lines represent runID, while the blue lines are for runAO. In the top panel, the two insets show maps of the logarithm of density in runID just before (panel d) and after (panel e) the development of the interchange instability (see text).

Open with DEXTER
In the text
thumbnail Fig. 2.

Left column: temperature (panel a), magnetic field (panel b), plasma β(panel c), and ratio of thermal to radiative pressure (panel d) as a function of density, for every cell in the computational domain at the epoch of second core formation. The IMHD simulation is represented by the red colours, while the blue shades are for the NIMHD run. The green colours correspond to areas where IMHD and NIMHD results agree within 10%. Each data set is delineated by a solid contour line which outlines the data distributions. The dark and light colours give an indication of the positions of the cells in the simulation box according to the θ = cos−1(z/r) angle: the light colours denote cells close to the equatorial region (π/4 < θ < 3π/4) while dark colours show cells in the polar regions (θ < π/4 or θ > 3π/4). The dashed lines in panels a and b represent the time evolution of the central (densest) cell inside the mesh. The thin black line in panel b is the power law predicted from magnetic flux conservation in a contracting gas sphere. Center and right columns: radial distributions of various quantities for every cell in the computational domain. As in the left column, red colours are for runID while blue colours are for runAO. In panels g and h additional lines show the integrated enclosed mass and angular momentum, respectively, in successive spherical shells going outward from the centre of the system.

Open with DEXTER
In the text
thumbnail Fig. 3.

Slices through the centre of the domain, comparing the morphologies of the protostellar system in runID (columns 1 and 3) and runAO (columns 2 and 4) on three different spatial scales at the epoch of second core formation. Panels a–h: display a wide region around the first Larson core, the typical scale of a protoplanetary disk. Panels i–p: immediate vicinity of the first Larson core. Panels q–x: present the second Larson core and its close surroundings. Two left columns: side x–z views of the system, while the two right columns display the top x–y perspective. The coloured maps in each row alternate between representing the gas density and temperature. The arrows on the density maps depict the gas velocity field. Overlayed onto the temperature maps are magnetic field lines (left column) and AMR level contours (right column).

Open with DEXTER
In the text
thumbnail Fig. 4.

3D visualization of logarithmically spaced density isosurfaces in the inner-most region of the computational domain showing the structure of the second Larson core, in the case of ideal (top) and non-ideal (bottom) MHD. The isosurfaces have been cut half-way in the x-direction. The magnetic field lines are overlayed and have been coloured according to the magnitude of the magnetic field vector. The insets in the lower left corner of each panel show (with the same spatial scale) the central region of the system without the B field for a better view of the morphology. The density and magnetic field colour scales apply to both panels.

Open with DEXTER
In the text
thumbnail Fig. 5.

Slices of the gas density with velocity vectors in runID (panel a) and runAO panel c, about one month after the formation of the second Larson core. The area shown is the same as in Figs. 3s and 3t. In panel c, the yellow contour marks the disk limit, taken as ρ > 10−6.7 g cm−3, while the black dashed contour delineates the second hydrostatic core with ρ > 10−5 g cm−3. Panels b and d: slices of the gas temperature with magnetic field streamlines overlayed. Panel e: logarithmic map of the magnetic Toomre stability criterion Q mag inside the disk that forms around the second core in runAO. The grey-shaded areas indicate regions in the disk where the epicyclic frequency ω is imaginary and no Q could be computed. The yellow and dashed black contours are the same as in panel c. Panel f : radial profile of the azimuthal velocity for all the cells inside the runAO second core disk. The colours code for the mass contained in a particular region of the plot. A Keplerian velocity profile is overlayed (black solid line).

Open with DEXTER
In the text
thumbnail Fig. 6.

Hammer projections (runAO only) of the mass accretion rate (left column), the radiative flux (middle column), and the ratio of radiative to accretion flux (right column). The first row is for the first Larson core, while the second row shows the second Larson core just after its formation. The third and fourth rows show the second Larson core 24 days after formation and the accretion flow at the edge of the disk around the second hydrostatic core, respectively. The green colours in panels g and i indicate negative values.

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

Panel a: fraction of cells inside the mesh where the optical depth is being limited as a function of time (red solid line). The dashed black line shows the density at the centre of the system as a function of time. Panel b: number of cells in each level (grey) and the number of cells where the optical depth floor is operating (red), at a time of 28.180 kyr, when the first Larson core is formed. Panel c: relative difference in 2D histograms of gas temperature as a function of density for all the cells in a simulation with optical depth limitation and a second simulation without, at t = 28.180 kyr. The colour scale gives a measure of ℛ = Nlimited/Nnot limited − 1, where N is the number of cells binned inside a (ρ, T) pixel for the two different simulations. Panel d: same as panel c but in the case of the Rosseland mean opacity as a function of density.

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

Panel a: fraction of cells inside the mesh where ηA (red) and ηO(blue) are being modified to prevent the MHD timestep from becoming too small, as a function of time. The dashed black line shows the evolution of the central density. Panel b: number of cells per AMR level (grey) and the number of cells where the ambipolar (red) and ohmic (blue hatched) diffusion timestep floor is operating, at a time of 28.2 kyr. Panel c: relative difference in 2D histograms of magnetic field strength as a function of density for all the cells in a simulation with Δt flooring and a second simulation without, at t = 28.2 kyr. The colour scale is analogous to that of Fig. A.1. Panel d: same as panel c but in the case of the ambipolar (red-blue) and ohmic (green-brown) resistivities.

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

Temperature as a function of density, for every cell in the computational domain (ideal MHD case). The simulation using NJeans = 32 cells per Jeans length is represented by the blue area, while the red region is for the run with only 16 cells per Jeans length. Each data set is delineated by a solid contour line which outlines the data distributions. The two snapshots were taken at similar evolution times, chosen to be just after the NJeans = 16 run has departed from its initial adiabatic track. The dashed lines represent the time evolution of the central (densest) cell inside the mesh (these tracks continue beyond the time of the snapshots to provide a wider context). The black arrow indicates the place where the low-resolution track departs from its original adiabat.

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

Map of the ambipolar (filled blue/red contours) and ohmic (green) Reynolds numbers close to the first Larson core. The light grey lines represent the magnetic field.

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

Maps and contours showing the morphologies of the cores using two different definitions. The coloured maps show the ratio of thermal to infalling ram (kinetic) pressure, while the black solid contour defines the region where the gas density exceeds density thresholds of ρcore = 10−10 g cm−3 for the first Larson core and ρcore = 10−5 g cm−3 for the second Larson core. The panels are: (a) runID first core, (b) runAO first core, (c) runID second core, (d) runID second core. We note the difference in spatial scales between panels a and b.

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.