Magnetic fields during the early stages of massive star formation – I. Accretion and disc evolution
We present simulations of collapsing 100 M mass cores in the context of massive star formation. The effect of variable initial rotational and magnetic energies on the formation of massive stars is studied in detail. We focus on accretion rates and on the question under which conditions massive Keplerian discs can form in the very early evolutionary stage of massive protostars. For this purpose, we perform 12 simulations with different initial conditions extending over a wide range in parameter space. The equations of magnetohydrodynamics (MHD) are solved under the assumption of ideal MHD. We find that the formation of Keplerian discs in the very early stages is suppressed for a mass-to-flux ratio normalised to the critical value below 10, in agreement with a series of low-mass star formation simulations. This is caused by very efficient magnetic braking resulting in a nearly instantaneous removal of angular momentum from the disc. For weak magnetic fields, corresponding to , large-scale, centrifugally supported discs build up with radii exceeding 100 AU. A stability analysis reveals that the discs are supported against gravitationally induced perturbations by the magnetic field and tend to form single stars rather than multiple objects. We find protostellar accretion rates of the order of a few 10 M yr which, considering the large range covered by the initial conditions, vary only by a factor of 3 between the different simulations. We attribute this fact to two competing effects of magnetic fields. On the one hand, magnetic braking enhances accretion by removing angular momentum from the disc thus lowering the centrifugal support against gravity. On the other hand, the combined effect of magnetic pressure and magnetic tension counteracts gravity by exerting an outward directed force on the gas in the disc thus reducing the accretion onto the protostars.
keywords:hydrodynamics – magnetic fields – methods: numerical – stars: formation – stars: massive.
The question of how massive stars form is still a highly debated field of research (e.g. Zinnecker & Yorke, 2007). It is believed that massive star formation takes place in high mass molecular cloud cores with masses ranging from roughly 100 M up to a few 1000 M. Characteristic for such cores are sizes of few 0.1 pc and peak densities up to 10 cm (e.g. Beuther et al., 2007). Furthermore, from observations it is known that the interstellar medium as a whole is magnetised (see Beck, 2009, for a recent overview). Also the star forming cloud cores partly reveal a significant magnetisation. The importance of the magnetic field can be estimated by the mass-to-flux ratio normalised to the critical mass-to-flux ratio (Mouschovias & Spitzer, 1976):
Observed mass-to-flux ratios in high-mass star forming cores are typically only slightly supercritical with (Falgarone et al., 2008; Girart et al., 2009; Beuther et al., 2010) indicating a significant influence of magnetic fields on the star formation process. In magneto-hydrodynamical simulation, however, also higher values of up to 20, i.e. weaker magnetic fields, are found (e.g. Padoan et al., 2001; Tilley & Pudritz, 2007). Another common feature of star forming cores are their slow rotation velocities. Observed cores have rotational energies normalised to the gravitational energy which scatter around a mean of about 0.01 (Goodman et al., 1993; Pirogov et al., 2003).
In the field of low-mass star formation there is a great number of observations of discs and large-scale outflows which are the basic keystones of the widely accepted disc accretion scenario for low-mass star formation (see, e.g. the reviews by Mac Low & Klessen, 2004; McKee & Ostriker, 2007). A similar formation scenario for massive stars is supported by a growing number of discs and bipolar outflows observed around high-mass protostellar objects (see Beuther & Shepherd, 2005; Cesaroni et al., 2007, for recent reviews). Therefore, it is worthwhile to study the influence of magnetic fields and rotation on the formation of discs and outflows in the context of massive star formation with numerical simulations.
For low-mass star forming regions (M 1 M), the influence of magnetic fields on the collapse of rotating cloud cores, the subsequent formation and evolution of discs and the launching of outflows has received extensive attention (e.g. Allen et al., 2003; Banerjee & Pudritz, 2006; Price & Bate, 2007; Mellon & Li, 2008; Hennebelle & Fromang, 2008; Hennebelle & Teyssier, 2008; Hennebelle & Ciardi, 2009; Duffin & Pudritz, 2009; Machida et al., 2011). All authors find a more or less significant influence of magnetic fields on the evolution of discs surrounding the protostars. The perhaps most important result of these studies is that for a mass-to-flux ratio of the formation of Keplerian discs is largely suppressed. This so-called magnetic braking catastrophe (Allen et al., 2003; Mellon & Li, 2008) turned the traditional angular momentum problem upside down: In highly magnetised cores magnetic braking seems to be so efficient that large-scale Keplerian discs, commonly observed around low-mass protostars, cannot form. Low-mass star formation simulations also reveal a strong impact of magnetic fields on the fragmentation properties of discs. In particular, strong magnetic fields tend to suppress disc fragmentation even in the presence of initial density perturbations (e.g. Hosking & Whitworth, 2004; Machida et al., 2005; Hennebelle & Teyssier, 2008; Duffin & Pudritz, 2009). These results are in contrast with the observational fact that a large fraction of low-mass stars are binaries (e.g. Duquennoy & Mayor, 1991).
Numerical studies of the influence of magnetic fields on massive star formation, however, have received attention only recently (Banerjee & Pudritz, 2007; Peters et al., 2011; Hennebelle et al., 2011). Banerjee & Pudritz (2007) study the very early evolution of a protostar, its surrounding disc and the outflow performing a simulation with extremely high spatial resolution. The protostellar evolution over a timescale of some 10 yr is examined by Peters et al. (2011) and Hennebelle et al. (2011), the latter authors focussing on the effect of magnetic fields and turbulence while Peters et al. (2011) study the interplay of magnetic fields and radiation. Here, we systematically study the influence of rotation and magnetic fields on the formation and accretion history of massive stars and examine the question under which conditions massive Keplerian discs can form in an already very early stage of protostellar evolution. We perform a series of collapse simulations with different initial rotational and magnetic energies following the protostellar evolution over 4000 yr. Rotational and magnetic energies are selected in a way to cover a large range in parameter space in accordance with observations and numerical simulations. Thus, we are able to detect systematic dependencies on the initial conditions. Furthermore, the simulations serve as a useful guide to select representative parameter sets for subsequent and more detailed studies.
The paper is organised as follows. In Section 2 the numerical techniques and the initial conditions of the runs performed are described briefly. The results of the simulations are presented in Section 3. First, the time evolution of the discs is presented for two representative simulations. Next, we analyse the velocity structure and the magnetic properties of the gas in the midplane. Afterwards, the accretion histories of the protostars as well as the effects of disc fragmentation are examined. In Section 4 the results are discussed in a broader context and are compared to other numerical and observational studies before we summarise our results in Section 5.
2 Numerical techniques and initial conditions
2.1 Numerical methods
We present 3-dimensional, magnetohydrodynamical (MHD) simulations of collapsing molecular cloud cores using the AMR code FLASH (Fryxell et al., 2000). We solve the equations of ideal MHD including self-gravity. The MHD-solver used preserves positive states and applies well for highly supersonic, astrophysical problems (Bouchut et al., 2007; Waagan, 2009; Bouchut et al., 2010; Waagan et al., 2011). Using 13 levels of refinement results in a maximum spatial resolution of 4.7 AU. We apply a refinement criterion that guarantees that the Jeans length
is always resolved with at least 8 grid cells. Here denotes the sound speed and the gravitational constant. A resolution of 4.7 AU guarantees that we can resolve the Jeans length at g cm corresponding to the density limit where molecular gas starts to get optically thick (see also the Appendix). In addition, a second refinement criterion implemented in FLASH is applied making use of the second derivative of the density field. This criterion guarantees that protostellar discs and outflows are resolved on the highest level of refinement allowing us to draw detailed conclusions about their properties.
We follow the collapse of the molecular cloud core with increasingly higher resolution until the maximum refinement level is reached. As soon as the density exceeds the critical value of
we create a sink particle (for details see Federrath et al., 2010). Subsequently, all gas with within a radius of AU around the sink particle is accreted onto the sink provided that it is gravitationally bound. The density cut-off at guarantees that the Jeans length is resolved everywhere by about 8 grid cells. The magnetic field is not changed during the accretion process to avoid the violation of the divergence-free condition. Keeping the magnetic field unchanged is motivated also physically by the onset of ambipolar diffusion at densities of 10 g cm (e.g. Nakano et al., 2002). Ambipolar diffusion would allow the gas to slip against the field lines, hence leaving them in the environment of the protostar. We briefly mention that due to numerical diffusion the increase of the magnetic field strength in the centre ceases at some point during the simulations therefore avoiding an unlimited growth of the Alfvénic velocity
In order to follow the thermodynamical evolution of the gas properly, we apply a cooling routine which includes the effects of molecular cooling and dust cooling as well as a treatment for optically thick gas (see Banerjee et al., 2006, for details). In this routine the cooling rates are computed after each hydrodynamical timestep and the state variables are updated accordingly for all grid cells. Using a subcycling scheme ensures that during one subcycle timestep the thermal energy of each cell does not change by more than 20%. We mention that the cooling is not applied to the ambient medium enclosing the molecular cloud core to prevent it from collapsing as well.
As soon as an outflow is launched from the protostellar disc, we introduce a minimum gas density threshold within a radius of 67 AU around the simulation centre. Within this radius the gas density is kept artificially above a minimum value of g cm. Whenever a cell density falls below this threshold during one timestep, mass is added to this particular cell until the limit of g cm is reached. Hence, we avoid the hydrodynamical timestep falling to prohibitively small values which otherwise would result from high Alfvénic velocities caused by low densities (see Eq. 4). The effects of this artificial density threshold will be discussed in detail in Section 4.1.
2.2 Initial conditions
With the numerical methods described above we are able to follow the collapse of a magnetised molecular cloud core down to a resolution of 4.7 AU. The core has a mass of 100 M and a diameter of 0.25 pc surrounded by gas with a density 100 times lower than the density at the core edge. The cubic simulation box has a length of 0.75 pc with the core located in its centre. Initially the core is resolved by a grid with a spacing of 302 AU fulfilling the applied refinement criteria. To assure pressure equilibrium at the edge of the core, the ambient gas has a temperature of 2000 K which is 100 times higher than the initially constant core temperature of 20 K. Calculating the Jeans mass using a temperature of 20 K and the average density of the core shows that the core contains about 56 Jeans masses and is therefore highly gravitationally unstable. The core itself has a density profile
Hence, the density in the centre
Beside the mass distribution, which is not changed in this work, there are two more crucial parameters affecting the evolution of the cloud core, namely the magnetic field strength and the rotational energy. Initially, all simulations have a magnetic field pointing into the z-direction where declines outwards with the cylindrical radius R as
This corresponds to a constant in the equatorial plane where is the thermal pressure. To guarantee , in the initial configuration is constant along the z-direction. We note that the chosen field configuration is compromise between observations and a numerically feasible setup. However, by the time the first sink particles form, the field has reached a self-consistent, hourglass-like configuration as observed in massive star forming regions (e.g. Girart et al., 2009; Alves et al., 2011). Therefore, we expect the chosen initial configuration of the magnetic field not to change the results significantly compared to an initially hourglass-shaped configuration. Initially, the cores are also rotating rigidly with the rotation axis parallel to the magnetic field. Neither the density field nor the velocity field have any perturbations at the beginning of the simulation so that a uniform and inward directed collapse can be expected.
We performed 12 simulations changing the initial magnetic field strength and the initial rotation frequency . For comparative purposes we also performed a simulation with zero magnetic field and zero rotation. The initial values of all simulations presented here are shown in Fig. 1. As the magnetic field is changing with radius, Fig. 1 shows the initial ratio of magnetic energy to gravitational energy and the corresponding normalised mass-to-flux ratio (Eq. 1) averaged over the entire core. The initial magnetic field strength at the centre of the cores ranges from G up to 1 mG. For comparison with , Fig. 1 also shows the ratio of rotational energy to gravitational energy beside the angular frequency. Both and are calculated numerically and are always well below 1 as the gravitational energy is significantly larger than the magnetic and rotational energies. In addition, it can be seen from Fig. 1 that magnetic and rotational energies extend over more than two orders of magnitude and bracket the line where equals . The initial values and the names of all runs are listed in Table 1. In the nomenclature of each run the first number give the mass-to-flux ratio and the second number multiplied by a factor of 100.
The considered values of rotational energies coincide reasonably well with values from observations, ranging from up to 1.4 with a mean around 0.01 (Goodman et al., 1993; Pirogov et al., 2003). Furthermore, observations of high-mass star forming regions typically reveal mass-to-flux ratios 5 and magnetic field strengths between a few 100 G and a few mG (Lai et al., 2001; Curran & Chrysostomou, 2007; Falgarone et al., 2008; Girart et al., 2009; Beuther et al., 2010). Here, we explore this range with a number of simulations () but also consider initial configurations with significantly weaker magnetic fields (). We recognize that turbulent motions are significant for massive clumps (e.g. Caselli & Myers, 1995). However, since the focus of this work is to investigate the formation of discs in magnetised environments, we will delay the inclusion of turbulence to a subsequent paper.
In this section, we present the results of our collapse simulations where we in particular focus on the evolution of the discs and the accretion onto the protostars. The evolution and properties of the outflows launched from the protostellar discs will be discussed in a subsequent paper. In the following, we mainly limit our consideration to the phase after the first sink particle has been formed and only shortly describe the initial collapse phase. After the formation of the first sink particle the simulations run for further 4000 yr (except the runs 2.6-20 and 2.6-4 which are run for 3000 yr only) to study the time evolution of the outflow, the formation of a protostellar disc and of potential secondary sink particles as well as the accretion history of the sink particles. The computational cost of all simulations presented here sum up to about 1 000 000 CPU-h.
The typical timescale for the collapse of the considered cloud core is its free-fall time
where we have used the central density of g cm. In Table 2 we list the formation time t of the first sink particle for all runs performed. As can be seen, t is longer than by a factor of 1.1 to 1.4. As run inf-0 has a prolonged collapse time compared to as well, we performed a careful numerical test showing that t can be reduced by a few 100 yr using a higher initial resolution. However, the main reason for the delay of the collapse is probably not a numerical issue but rather the build-up of a strong pressure gradient in the centre of the core counteracting gravity. In contrast, the unexpected fact that t for run inf-0 is somewhat longer than for the runs 26-0.4 and 26-0.04 is most like a numerical issue. Turning on weak magnetic field in a test run with no rotation results in a slightly shortened collapse time compared to inf-0. This suggest that the somewhat longer t in run inf-0 is an intrinsic effect of the numerical scheme we will not follow up further here. However, in each simulation subset (equal or ) there are physically well motivated trends in t recognisable. As can be inferred from Table 2, the collapse gets slowed down with increasing magnetic field strength and increasing amount of rotational energy, both counteracting gravity. Interestingly, the variation of the rotational energy by two orders of magnitude at a fixed magnetic field strength changes t by no more than roughly 0.8 kyr, i.e. by 5%. In contrast, increasing the magnetic energy by a factor of 100 for fixed rotational energies prolongates the collapse time by roughly 5 kyr. However, as all cores are supercritical () and have rotational energies well below the gravitational energy (see Figure 1), one cannot expect the collapse time t to be significantly longer than the free-fall time. In the following, all times refer to the time elapsed since the formation of the first sink particle at t.
3.1 Global disc properties
As a first step, we analyse the time evolution of two representative simulations after the formation of the first sink particle. For this purpose, we consider run 26-4 and run 5.2-4 with equal initial rotational energies but different magnetic field strengths. To begin with, we compare the properties of the gas in a well defined region of interest. This region is described by a cylinder around the simulation centre with a height of 47 AU above and below the midplane. We mention that, although this height does not reflect the real disc scale height, we use it for three reasons. Firstly, this choice provides a numerically simple calculation of the desired quantities. Furthermore, it guarantees that for different simulations and times identical areas are compared. Finally, it is used as we cannot determine a well defined disc scale height as even in a single run the disc height varies in time and radial position. For example, the frequently used approximation of the scale height (e.g. Cesaroni et al., 2007)
varies between a few AU and 100 AU. This is in accordance with other estimates of the disc height made the authors. Hence we argue that our choice of 47 AU, located in the middle of this range, is reasonable.
In the following we will denote this region simply as the disc. The quantities in the disc are averaged vertically and azimuthally before consideration.
From the results shown in Fig. 2, it can be inferred that in run 26-4 an accretion shock occurs where both, column density and temperature experience a sudden increase. This is also true for run 5.2-4, although here the jumps in density and temperature is somewhat smoother. In both cases, however, the shock front moves outwards as time evolves. In run 26-4 the density profile inside the accretion shock is nearly flat with values around a few 10 cm corresponding to densities of 10 g cm. For run 5.2-4 the density profile is declining outwards and seems to decrease slightly over time. In both cases the mass-weighted temperature in the midplane increases by about one order of magnitude up to a few 100 K.
The temperature increase in the inner region is caused by the gas getting optically thick at density around 10 g cm (see Banerjee et al., 2006, for details of the applied cooling function). Therefore, the gas in the disc looses its ability to cool efficiently and cannot fast enough radiate away the thermal energy transfered to it by compression work. This results in temperatures of up to a few 100 K in the inner disc region consistent with observational results (see e.g. the review of Cesaroni et al., 2007, and references therein). Further out in the disc the gas experiences a strong temperature increase due to shock heating. At the accretion shock kinetic energy is converted quite fast into thermal energy. Assuming that gas gets decelerated by 1 km s (see Fig. 2) and the corresponding kinetic energy is immediately transformed into thermal energy, this would result in a temperature increase of roughly 90 K. This is in accordance with the observed increase of about 30 K regarding to the fact that a good fraction of the energy is radiated away and that the conversion into heat does not happen all at once.
In Fig. 3 we show the magnetic properties of the discs. For the calculation of and the absolute values are used for averaging as both components have opposite signs above and below the disc and hence would cancel themselves out. For comparative purposes we plot (initially the only component) and either or in each panel. For the comparison we only consider the radial range above 10 AU as the region further inside is subject to resolution effects. As can be seen, the different components in the midplane reach values up to a about 1 G. is somewhat larger in run 5.2-4 than in run 26-4 due to the five times higher, initial field strength in run 5.2-4. In both runs, however, (bottom panel of Fig. 3) is of the order of . is created by the inwards drag of the magnetic field during the collapse and later during the accretion process. In contrast to , the toroidal components of both runs differ significantly (see top panel of Fig. 3). In run 26-4 is the dominant component in the region within the accretion shock being larger than and by up to one order of magnitude. This is due to the fast rotation of the disc winding up the poloidal field (see Fig. 4 and text below). In contrast, in the strongly magnetised run 5.2-4, is mostly smaller than which is attributed to the lower rotation velocity in this case (see Fig. 4). We also mention that all components show signs of an accretion shock moving outwards with time in accordance with the density and temperature field shown in Fig. 2. This is caused by the tight coupling of magnetic fields and matter due to ideal MHD. Furthermore, a comparison with the thermal pressure shows that in both runs the magnetic pressure in the region within the accretion shock is larger than or at least equal to the thermal pressure. This implies that the magnetic field must play a significant role in the evolution of the circumstellar disc and in the accretion history of the protostars.
In Fig. 4 the time evolution of the rotation (), radial () and Keplerian velocity () is shown for both runs.
For run 26-4 (left panel) a rotationally supported disc builds up with rotation velocities close to the Keplerian velocity. The maximum radius, where equals increases over time consistent with the build-up of a Keplerian disc. At the same time the radial velocity in the inner region nearly drops to zero. We emphasise that considering the behaviour within the central 10 AU is not conclusive as we reach the resolution limit at such small radii. The overall situation changes dramatically when considering run 5.2-4 (right panel). Here, no rotationally supported disc builds up with staying significantly below all the time. In fact, the absolute value of is nearly as high as denoting gas which is almost in free-fall. Nevertheless, a flat disc-like density structure builds up.
The dramatically different evolution of the gas in the midplane of both runs is also indicated in Fig. 5. Here we show the integrated column density, the magnetic field strength and the velocity field along a slice in the midplane.
For run 26-4 (upper panel of Fig. 5) a well defined Keplerian disc with more or less sharp boundaries develops. The velocity field within the disc shows a significant rotational component as already seen in Fig. 4. Despite the strong fragmentation occurring in the inner 200 AU, the overall disc-like structure is maintained. In contrast, in run 5.2-4 (lower panel of Fig. 5) there is no evidence for the development of a Keplerian disc and the column density increases more or less smoothly towards the centre with nearly radial infall. The right panel of Fig. 5 shows the z-component of the magnetic field in the midplane. As can be seen, a field strength of up to 1 G is reached with slightly higher values for run 5.2-4. The close coupling of magnetic fields and matter due to ideal MHD is especially pronounced in the top panel of Fig. 5 where the strong increase in the column density at r 350 AU coincides quite well with a jump in . In Fig. 6 the coupling of gas and magnetic field over the whole density range is demonstrated even more clearly. Apart from local variations, the average of and , calculated in density bins of equal size in log-space, scales roughly as or slightly weaker over more than 6 order of magnitude in density. This is in accordance with the scaling of a magnetic field in case of a spherical collapse under the condition of ideal MHD.
The bubble like features seen in run 5.2-4 can be explained as follows: during the accretion process the magnetic field is dragged inwards with the gas. As we consider ideal MHD, the field lines cannot diffuse outwards leading to a constant increase in the field strength in the centre. Beyond a certain point in time, the magnetic pressure is strong enough to effectively counteract gravity and starts to push material outwards. Indeed, analysing the component of the magnetic field in the bottom right panel of Fig. 5 shows that the outwards moving gas is associated with a strong . This causes to decrease slightly in the very centre whereas at larger radii the magnetic field strength increases over time (see right panel of Fig. 3).
The two simulations presented here show significant differences although the initial magnetic field strength is varied by a factor of 5 only. These systematic differences also show up in a more or less pronounced way when comparing all simulations with each other. We will discuss these differences in detail in the next section where the complete set of simulations with varying initial conditions is considered.
3.2 Velocity structure
Next, we focus on the effect of the initial conditions on the velocity field in the midplane around the central sink particle. We omit the time evolution as the behaviour of each individual run is qualitatively similar to one of the two runs shown before and concentrate on the situation after 4000 yr. As already shown in Section 3.1, even small changes in the initial configuration of the core cause characteristical differences. The effects of rotational and magnetic energy on the velocity structure at the end of the simulations can be seen in Fig. 7.
Beside the radial velocity, the ratio / is shown as well. Decreasing the initial amount of rotational energy for runs with fixed magnetic field strength (see individual panels of Fig. 7) reduces the centrifugal support against gravity resulting in lower values of / and consistently in higher infall velocities. Additionally, the differences for runs with varying but fixed seem to decrease when the initial field strength is increased. Furthermore, comparing runs with fixed but varying field strength shows that the centrifugal support, i.e. / decreases with increasing field strength. This effect is pronounced strongest for high initial rotational energies.
In principle, the simulations can be divided into two sets depending on the initial mass-to-flux-ratio . In general, for runs with centrifugally supported discs develop whereas for the formation of Keplerian discs is largely suppressed and only sub-Keplerian discs form at this early stage (see Section 4.2 for a detailed comparison with other numerical work). An overview over the observed dependency of disc formation and fragmentation on the initial conditions is given in Fig. 8.
We note that in the runs 26-0.04 and 10-0.4 the structure formed in the midplane strongly resembles that of centrifugally supported discs but with rotational velocities well below the Keplerian velocity. Thus, we do not denote them as Keplerian discs.
In the following, we qualitatively describe the track of a fluid particle moving along the midplane towards the centre. For all simulations considered in this work, the fluid particle first gets accelerated inwards until a radius of some 100 AU (depending on the specific simulation considered) is reached. At this radius the gas experiences a deceleration, meaning its infall motion slows down. This region can be identified as a so-called magnetic barrier (Mellon & Li, 2008) and is found in all simulations. This kind of barrier is different from a centrifugal barrier as here the rotation velocity is well below the Keplerian velocity which would be necessary to balance gravity. The reason for deceleration at such a magnetic barrier is twofold. Firstly, an outward directed thermal pressure gradient accounts for a part of the deceleration. This can be inferred from the significant increase in density and temperature at the corresponding radius (compare Fig. 2). Furthermore, the magnetic field itself slows down the infall motion via the combination of an outward directed magnetic pressure gradient and the effect of magnetic tension which can be inferred from Fig. 3. All components of the magnetic field experience a sudden increase in the region of deceleration, thus resulting in outward directed magnetic forces. Both magnetic pressure and magnetic tension increase in strength for smaller with the latter starting to dominate for low . For cases with not considered here, we expect the collapse of the core perpendicular to the field lines to be prevented completely by the Lorentz force (Mouschovias & Spitzer, 1976).
For radii within the magnetic barrier the velocity profiles start to differ significantly from each other. For weak magnetic fields (top panel of Fig. 7) the gas infall speed stays roughly constant. Due to angular momentum conservation this results in an increased rotation velocity which in consequence leads to another slow down of the infall motion. For the highest rotational energies a centrifugal barrier is encountered bringing the infall completely to halt. As seen in the top panel of Fig. 7, the occurrence and extension of the centrifugal barrier depend strongly on the initial setup. On the other hand, an increase in centrifugal support and hence a slowdown of infall gives magnetic braking more time to operate (Mouschovias & Paleologou, 1980). This is seen in the sudden drop of / for runs with causing a loss of centrifugal support and thus again speed-up of the infall motion. This is also the case for runs with low (bottom panel of Fig. 8) where magnetic braking starts to act efficiently directly after passing the magnetic barrier so that no centrifugal barrier occurs any more. The innermost drop of (r 10 AU) without a corresponding increase in / is certainly caused by the limiting effect of numerical resolution as this region is only marginally resolved by about 3 grid cells.
The reason for the effective magnetic braking can be seen in Fig. 9 showing the edge-on view of the runs 26-4 and 5.2-4 at two different times. As can be seen, the collapse of the gas has dragged in the magnetic field in the midplane producing a long magnetic lever arm (Allen et al., 2003) with a large radial component (see also Fig. 3). This lever arm connects the outer, slowly rotating region with the inner, fast rotating region, hence significantly enhancing the magnetic braking efficiency.
As shown before, in the case of strong magnetic fields () magnetic braking is so efficient that large centrifugally supported discs do not form at this early stage. This is known as the magnetic braking catastrophe and was reported previously by several authors for ideal MHD simulations (Allen et al., 2003; Mellon & Li, 2008; Hennebelle & Fromang, 2008; Hennebelle & Ciardi, 2009). For a better qualitative and quantitative understanding of the magnetic braking effect we calculate the z-component of the two main torques acting on the disc, i.e. the torque exerted by the infalling gas
and the torque exerted by magnetic fields
For the calculations we integrate over a disc with a height of 47 AU above and below the midplane. We omit the gravitational torque as it is by two to four orders of magnitude smaller than and due to the nearly axisymmetric gravitational potential.
In Fig. 10 the time evolution of and is shown for the runs 26-4 and 5.2-4. The torques are averaged azimuthally and plotted against the radius. To allow for a better comparison with we plot . A positive denotes a flux of angular momentum into the disc while for a negative angular momentum is removed from it. Hence, in both runs the magnetic field is removing angular momentum from the disc, i.e. slowing down its rotation () whereas the gas exerts a positive torque on the disc. It can be seen that the torques increase steadily with time (except for run 26-4 at 3000 yr). For run 26-4 is nearly always larger than the magnetic torque whereas in run 5.2-4 is roughly balanced by from the very beginning. These differences show up even more clearly in Fig. 11 where the torques at the end of the simulations are shown.
For runs with a weak magnetic field (top left panel of Fig. 11) the torque exerted by the gas is nearly everywhere above the magnetic torque. Thus, there is a net flux of angular momentum into the disc leading to the observed build-up of centrifugally supported discs. Only for run 26-20 equals or even exceeds in the inner region although only on a low level. This is due to the large Keplerian disc which has already built up in this run having only very small infall velocities (see top left panel of Fig. 7). The sharp jump of around r = 200 AU is caused by the accretion shock at the edge of the disc where drops to zero (see top left panel of Fig. 7).
Analysing the different panels in Fig. 11, it can be seen that at large radii the gas torques are always larger than . This implies that the magnetic field has only a small effect on the collapse in the outer parts. In contrast, at smaller radii the effect of magnetic fields gets more and more pronounced. With decreasing the magnetic torque approaches and in particular for the runs with and 2.6 is very close to . This is attributed to the fact that an equilibrium between and is reached where as much angular momentum is removed by magnetic braking as it is added due to the gas infall. Hence, in the case of strong magnetic fields the net angular momentum flux is roughly zero preventing a Keplerian disc from forming.
To summarise, centrifugally supported discs at very early times only form in simulations with . In these cases magnetic braking is not strong enough to remove angular momentum at high enough rates and hence is not able to prevent Keplerian disc formation (see Fig. 8 for an overview).
3.4 Accretion rates
Closely related to the velocity structure of the matter around the sink particles is the accretion onto the particles themselves. The general behaviour of the accretion rate with varying initial conditions can be inferred from Table 2 where we list the totally accreted masses and the corresponding time averaged accretion rates. Additionally, in Fig. 12 the accretion rates of all runs are plotted against the initial mass-to-flux ratio.
As listed in Table 2, in neither of the runs more than 4 M are accreted during the simulation resulting in time averaged accretion rates around a few 10 M yr. Interestingly, the accretion rates do not vary by more than a factor of about 3 between the different simulations. This is remarkable, considering the large range in parameter space covered by the initial conditions (Fig. 1). For each set of simulations with equal there is also a rough correspondence between the accretion rate and the infall velocity shown in Fig. 7. As expected, higher infall motions result in higher accretion rates.
For a more detailed analysis of the accretion rates, we consider their time evolution in Fig. 13. We mention that for the runs 26-20 and 26-4 where more than one sink particle is formed only the accretion onto the first particle created is shown. A more detailed analysis of the accretion history in the case of fragmentation will be carried out in Section 3.5.
As can be seen in Fig. 13, there seems to be a slight decrease in the accretion rates by a few 10% over time. Numerical studies (e.g. Klessen, 2001; Schmeja & Klessen, 2004) as well as observational results (see Andre et al., 2000, for an overview) indicate that there is indeed a decline in the accretion rate by orders of magnitude, although timescales for this decline are typically of the order of several 10 yr and are therefore much longer than in our study. Additionally, there occur fast variations within a factor of about 2 around the mean value. This fast variations are caused by moderate density perturbations developing in the midplane. Each time a perturbation moves through the centre, it causes the accretion to vary around its mean. The variations are generally rather small (with one exception in run 26-0.4) and would probably be smoothed out in time by viscous effects in the inner disc not resolved here. However, we cannot exclude that the varying accretion rates would influence the protostellar evolution as proposed by stellar evolution calculations (e.g. Wuchterl & Klessen, 2001; Baraffe & Chabrier, 2010).
The only exception where the accretion rates decrease significantly over time are the runs with showing a sharp drop in the mass accretion after roughly 2500 yr. This is caused by the occurrence of magnetically driven bubbles in the midplane as shown exemplarily in the bottom panel of Fig. 5. However, in the other runs the accretion rates seem to decrease only slightly showing no sign that the magnetically driven outflow, which is launched from the protostellar disc shortly after the formation of the protostar, can significantly reduce mass accretion over time. Similar results in related work on magnetic fields in massive star formation are found by Peters et al. (2011) and Hennebelle et al. (2011) as well. This behaviour is due to the ongoing accretion through the disc which is nearly unaffected by the shut-off of gas infall from below or above the disc due to the outflow. This can be seen in Fig. 9 where we plot the velocity structure and magnetic field lines in a slice along the z-axis.
As can be inferred from Table 2 and from Fig. 12, there are some systematical trends in the accretion rates with changing initial conditions, although the overall variation is not very large. Increasing the overall rotational support against gravity, i.e. , for fixed results in lower accretion rates. This is in agreement with the increase in / and the decrease in in the surrounding disc as shown in Fig. 7. As already observed for the velocity structure in the midplane, the differences in the accretion rates for runs with different but fixed decrease with increasing magnetic field strength as can be seen in Fig. 12. This is due to the efficient magnetic braking removing angular momentum at roughly the same rate as it is transported inwards (see Fig. 11). As a consequence, the accretion rates are roughly independent of the initial amount of angular momentum. An even further increase in the field strength would probably result in even lower accretion rates than the 4 - 5 M yr observed for = 2.6 due to stronger magnetic forces counteracting gravity.
The accretion rates for varying but fixed show an interesting and less clear behaviour. For high initial rotational energies, i.e. and 0.04, respectively, the accretion rates first increase with decreasing and drop again for the case with . For the accretion rate increases from to before declining again for lower . We attribute this behaviour to two competing effects of the magnetic field. On the one hand, magnetic fields act to enhance accretion onto the protostar by magnetic braking reducing the centrifugal support against gravity. Hence, the effect of magnetic braking alone would cause increasing accretion rates with decreasing as it is indeed observed for low field strengths. The second effect influencing the accretion rates is the Lorentz force, i.e. the combination of magnetic pressure and magnetic tension, induced by strongly bent field lines (see Fig. 9). Magnetic pressure and magnetic tension counteract gravity by exerting an outward force on the gas resulting in reduced accretion rates. The strength of this effect increases with the field strength thus tending to lower accretion rates with lower .
The combination of both effects – magnetic braking enhancing accretion and the Lorentz force counteracting accretion – results in the observed behaviour of the accretion rates: By increasing the magnetic field strength at a given up to a certain critical value an equilibrium between the torques acting on the disc is reached where the removal of angular momentum by magnetic braking balances its inwards transport (see Fig. 11). Up to this value, an increase in the field strength is associated with increasing accretion rates as observed in our simulations. Further increase in the field strength beyond this point cannot enhance the magnetic braking efficiency any more. In fact, now the increase (decrease, if is considered) leads to declining accretion rates due to the growing strength of the Lorentz force counteracting gravity, in accordance with our findings of declining accretion rates for strong fields with . The exact value of , where this turnover occurs, depends on the initial amount of rotational energy and decreases with increasing .
3.5 Disc fragmentation
As mentioned earlier, the runs 26-20 and 26-4 show rapid fragmentation of the protostellar disc after the first protostar has formed (see Fig. 8). Due to a high amount of rotational energy and a weak magnetic field, magnetic braking can only remove a small amount of angular momentum leading to the formation of Keplerian discs with considerable extensions of a few 100 AU (see top left panel of Fig. 7 for comparison). As the mass load onto these discs exceeds their capability to transport material inwards by gravitational or viscous torques, the discs become unstable and fragment (e.g. Kratter et al., 2010). At the end of the simulations, i.e. after 4000 yr there are 10 sink particles in run 26-4 and 13 in run 26-20. All other simulations show no fragmentation so far although some of them form a Keplerian disc (compare Fig. 12).
The accretion histories for run 26-20 and run 26-4 are shown in Fig. 14. Run 26-20 exhibits a very symmetric disc fragmentation forming pairs of protostars at roughly the same time and opposite positions (as seen from the centre). Even in their further evolution each pair develops very similar as can be seen from the left panel of Fig. 14 where the lines of each pair are nearly indistinguishable. For run 26-4 only the evolution of the second and third sink particle is symmetric, while at later times there is no pairwise formation of sink particles any more due to a nonsymmetric evolution of the disc.
We note that in both runs only the first sink particle created has reached a mass above 1 M so far whereas all other particles have masses well below 0.1 M. For comparative purposes we also plot the totally accreted mass of all sinks formed in Fig. 14 (see also Table 2). As can be seen, in both runs more than 80% of the total mass is accreted onto the first sink particle.
After about 2500 yr, i.e. after the creation of further sink particles, there occurs a slight but nevertheless noticeable decrease in the accretion onto the first sink particle (see the red and black lines in the top left panel of Fig. 13). This fragmentation-induced starvation (Peters et al., 2010b) is caused by the surrounding sink particles soaking up the infalling material. This behaviour was also observed recently in related work on massive star formation (Peters et al., 2010a, 2011; Girichidis et al., 2011).
Beside the sink masses, Fig. 14 shows the calculated disc masses as well. For the calculation we define the disc as follows: Firstly, we determine the maximum cylindrical radius R where the gas falls below a density of g cm. The disc mass is then defined as the total mass of gas within a cylinder with a height of 47 AU above and below the midplane and a radius of R around the centre. As can be seen, the disc masses lie around 1 M and are therefore somewhat below the totally accreted sink masses.
Next, we study the stability of the discs against gravitationally induced perturbations. Disc stability is described by the Toomre parameter (Toomre, 1964) defined as
with the epicyclic frequency , sound speed , surface density , and gravitational constant . Gravitational instability sets in when . As magnetic fields are present in the discs, one can define the magnetic Toomre parameter (Kim & Ostriker, 2001)
where is the Alfvénic velocity taking into account all components of the magnetic field. In the following, we analyse the stability of the midplane in the runs 26-20, 10-20, 5.2-20 and 26-0.4 as these simulations cover a wide range of initial conditions. In particular, we concentrate on the stabilising effect of the magnetic field by comparing and .
In Fig. 15 the face-on view of the discs after 4000 yr as well as the radial dependence of the Toomre parameters are plotted. and are calculated numerically by the azimuthally averaged values of , , , and and are shown for several times. We neglect the contribution of secondary sink particles to the column density in run 26-0.4 as their only effect would be to lower and in regions which are already unstable and fragmenting and therefore they would not qualitatively change the situation. The plot shows that within a radius of r 500 AU the Toomre parameter is almost everywhere below 1 for all runs and all times considered. Hence, in all runs fragmentation would be expected to occur. As this is not the case except for run 26-20, the Q-parameter seems not able to properly describe the stability properties of the discs. Taking into account the effect of magnetic fields, the situation displayed in the column density plots is described much more adequately. The value of is higher than by roughly a factor of 10 in all cases thus pointing to much more stable configurations. Nevertheless, for run 26-20 even is below 1 with one spatially confined exception at 2000 yr. Hence, in this case even the combination of thermal and magnetic pressure cannot stabilise the disc against fragmentation as shown in the upper left panel of Fig. 15. Fragmentation occurs in the innermost part (r 50 AU) close to the first sink particle and in a filamentary-like ring at r 200 AU. We note that the situation is somewhat different to run 26-4 where fragmentation seems to occur mainly in spiral like density perturbations closer to the centre (see top panel of Fig. 5). However, the stability analysis of run 26-4 shows that and are below 1 as well.
All other runs considered in Fig. 15 show no fragmentation so far. This agrees quite well with the findings of . However, there are some situations when is marginally smaller than 1 with the disc still showing no signs of fragmentation. Kratter et al. (2010), analysing the stability of discs under purely hydrodynamical conditions, find that discs can be indeed stable if Q is locally smaller than 1. The authors attribute this to the fact that = 1 indicates instability of axisymmetric perturbations in infinitely thin discs (Toomre, 1964). For thick discs, as it is the case in our simulations, they argue that the instability criterion is expected to be somewhat relaxed (Goldreich & Lynden-Bell, 1965). Hence, we expect our discs to be stable for Toomre values even slightly smaller than 1.
To summarise, the comparison of and shows that the stabilising effect of magnetic fields is needed to effectively prevent the discs from fragmenting. For weak magnetic fields () and high amounts of angular momentum, however, even magnetic fields cannot stabilise the discs against fragmentation. Furthermore, we expect that fragmentation would occur for other runs with low magnetic field strengths () as well if the evolution of the discs would be followed beyond 4000 yr. This assumption is confirmed by the results that the discs show a continuing growth due to ongoing infall and are only marginally stable ( 1).
4.1 Numerical caveats
The simulations considered so far were run with a maximum resolution of 4.7 AU. To test the resolution dependency of our results we performed two more simulations with the same initial setup as in run 26-4 but with a resolution varied by a factor of 4 in either direction. A detailed comparison of the results is shown in the appendix. The results of the resolution study suggest that the spatial resolution of 4.7 AU used in the simulations presented in this work is sufficient to properly follow the dynamical evolution of the protostellar discs and of the accretion rates.
The bubble like features seen in the lower panel of Fig. 5 also occur in other runs with strong initial magnetic fields and are artifacts of our assumption of ideal MHD. In particular, for runs with the bubbles show a significant influence on the accretion rates (bottom right panel of Fig. 13) reducing the accretion after 2500 yr. However, for runs with the dynamical influence of such bubbles seems to be rather limited as we cannot detect any significant changes in the accretion rates (see bottom left panel of Fig. 13). We therefore conclude that the accretion rates are reliable for runs with down to 5 and up to 2500 yr even for the both runs with . A possible solution to avoid the formation of bubbles like features could be the inclusion of the effects of non-ideal MHD such as ambipolar diffusion or ohmic dissipation as done in recent work (Duffin & Pudritz, 2009; Mellon & Li, 2009; Dapp & Basu, 2010). However, Nakano et al. (2002) argue that ohmic dissipation starts to act efficiently at particle densities above 10 cm which is well above the typical densities found in our discs. Hence, ohmic dissipation is not expected to be capable of significantly reducing the magnetic flux in the centre. Ambipolar diffusion, however, starts to act at lower densities (Duffin & Pudritz, 2009) and thus might be able to reduce the field strength in the centre before magnetically driven bubbles can occur.
As mentioned in Section 2.1, a minimum density threshold of g cm is applied after the outflow is launched. The exact amount of mass added artificially during the 4000 yr depends on the individual simulation but never exceeds a few 10 M corresponding to a mass rate of M yr which is significantly below the observed accretion rates of a few 10 M yr. Hence, we suppose that the dynamical influence of the density threshold should be negligible and its application is acceptable in order to avoid very small timesteps and in consequence significantly higher computational costs.
It is not clear to what extent our results are affected by the chosen density profile . Girichidis et al. (2011) find that their results strongly depend on the initial profile which they attribute to the varying importance of turbulence compared to gravity. As we do not include turbulent motions in our simulations, our results might depend somewhat less on the density profile, in particular with regard to the accretion rates which also vary rather moderately in the work of Girichidis et al. (2011).
4.2 Disc evolution
As demonstrated in Fig. 7, a transition from early-type, large-scale Keplerian discs to sub-Keplerian discs occurs around a normalised initial mass-to-flux ratios of 10 independent of the initial amount of rotational energy. This is in good accordance with a series of papers studying the evolution of low-mass magnetised disc. While Hennebelle & Fromang (2008) and Hennebelle & Ciardi (2009) find a value for between 5 - 10, below which Keplerian disc formation is suppressed, Allen et al. (2003) and Mellon & Li (2008) find no Keplerian discs for up to at least 10. Duffin & Pudritz (2009) studying the possible fragmentation of magnetised discs with an initial = 3.5 find sub-Keplerian rotation profiles as well in agreement with the results mentioned before. Although these simulations apply to low-mass star formation with core masses around 1 M, the observed maximum value of , for which the formation of large-scale Keplerian discs is suppressed, agrees remarkable well with our finding of .
For the sub-Keplerian discs observed in our simulations the infall velocities are generally significantly larger than the rotation speed, i.e. /, and close to free-fall (see right panel of Fig. 4 and Fig. 7). This is attributed to the highly gravitationally unstable configuration of the cores containing about 56 Jeans masses. Hence, they cannot be stabilised against gravitational collapse by thermal pressure alone in contrast to low-mass cores containing only 1 Jeans mass. Consistently, for highly magnetised low-mass cores Hennebelle & Fromang (2008) and Duffin & Pudritz (2009) find smaller infall velocities as well as / contrary to our results.
For , we find large-scale Keplerian discs in our simulations. This is probably a consequence of not having turbulence in our initial setup. Turbulence would most likely hamper the early formation of large rotating structures and delay it to later stages. However, in recent years there is an increasing number of observations of rotationally supported discs in massive star forming regions (e.g. Cesaroni et al., 2007, and references therein). As an example, discs with sizes of a few 100 AU and masses between 0.15 M (Fuller et al., 2001) and about 10 M (Shepherd et al., 2001) have been observed. These discs are similar in size to the discs obtained in our runs with weak magnetic fields and also show Keplerian-like rotation profiles. In addition, our disc masses of 1 M calculated in run 26-20 and run 26-4 (see Fig. 14) fit perfectly in the mass interval spanned by these observations. However, the protostars observed have masses around 5 - 10 M and are therefore a factor of a few more massive than our most massive sink particles. We attributed this to the fact that the observed disc/star-systems are in much later evolutionary stages than our systems.
In contrast, for sub-Keplerian discs are observed in our simulations. Again, we mention that our simulations end 4000 yr after sink particle formation, thus in a very early phase. An observer looking at such a system from edge-on will observe a flattened structure similar to a Keplerian disc but without the typical signatures of rotation (see right panel of Fig. 9). Indeed, there is a growing number of observations of such flattened structures reported in literature (see Cesaroni et al., 2007, and references therein), although these observations often refer to more massive structures and more evolved protostellar objects than present in our simulations. Due to the insufficient centrifugal support, these structures are not in equilibrium and show considerable radial infall motions, in accordance with our findings (compare bottom panel of Fig. 7).
The reason for the lack of centrifugal support in such sub-Keplerian discs is the strong magnetic torque acting on the midplane (see Fig. 11). Angular momentum is removed by magnetic braking at roughly the same rate as it is added by the infalling gas. The angular momentum removed from the midplane is partly deposited in the outflow and partly in regions further out which are connected to the inner parts by the magnetic lever arm (Allen et al., 2003) created through the equatorial pinching of the magnetic field (see Fig. 9). A possible way trying to reduce the magnetic braking efficiency would be to include non-ideal MHD effects like ohmic dissipation or ambipolar diffusion in the simulations. However, recent numerical work including ambipolar diffusion (Mellon & Li, 2009; Duffin & Pudritz, 2009), ohmic dissipation (Dapp & Basu, 2010) and also both effects (Li et al., 2011) show that even in the case of non-ideal MHD it is not possible to form Keplerian discs in such early stages. In fact, these authors find that the formation of rotationally supported discs in the case of strong magnetic fields is suppressed down to scales well below our resolution limit of roughly 10 AU. Hence, at the evolutionary stage considered in this work, we cannot expect to resolve a proper Keplerian disc even by including the effects of ambipolar diffusion or ohmic dissipation in our calculations.
The question now arises how Keplerian discs on 100 AU scale, frequently observed around massive protostars, can form if the cores have mass-to-flux ratios . The situation even intensifies as observed star forming regions usually have mass-to-flux ratios which are only slightly supercritical with a mean (e.g. Falgarone et al., 2008; Girart et al., 2009; Beuther et al., 2010). Following the long term evolution of highly magnetised low-mass cores, Machida et al. (2011) find large-scale Keplerian discs occurring after roughly 10 yr. The authors argue that in their case magnetic braking only redistributes the angular momentum within the collapsing cores but does not remove it. Additionally, they find that the outflows are too weak to ultimately remove a significant fraction of mass and angular momentum from the cores. Thus, most of the cores mass and angular momentum finally fall onto the midplane resulting in large-scale Keplerian discs even for high magnetic field strengths. However, we note that the authors might underestimate the total amount of mass and angular momentum transfered into the interstellar medium as only low velocity outflows are modeled in their simulations. Considering high velocity jets and radiation-driven outflows, the amount of material ejected into the interstellar medium could be considerably larger, thus impeding the formation of large Keplerian discs. For massive stars we expect this to be even more severe as here radiation-driven outflows are even more powerful and will be able to eject angular momentum – deposited in the envelope by magnetic effects – at even higher rates than low-mass protostellar outflows (e.g. Arce et al., 2007).
At the same time outflows also provide a possible solution of the disc formation problem by diluting the envelope in which the magnetic field lines are anchored (see also 6.2.2 in Mellon & Li, 2008). As the magnetic braking time (Mouschovias & Paleologou, 1980)
increases with decreasing envelope density, this allows large-scale, centrifugally supported discs to form in particular at later stages not considered in this work when most of the envelope has been blown away by powerful outflows. This would agree with the fact that most of the observed disc/star-systems (see Cesaroni et al., 2007, and references therein) are in a later evolutionary stage than our systems although the observations could also be biased by the fact that massive Class 0 objects are more difficult to detect. The picture of a successive growth of centrifugally supported discs during the evolution into Class I / II is also supported from the theoretical side by Dapp & Basu (2010). Therefore, we expect large-scale Keplerian discs to form in our runs as well if the evolution would be followed over a much longer time, thus relaxing the problem of catastrophic magnetic braking.
In summary, our simulations suggest that Keplerian discs do not form around massive protostars in the very early stage except for unusually weak magnetic fields. Nevertheless, we expect that centrifugally supported discs will build up during the subsequent evolution of the collapsing cores.
Another question not addressed so far is the problem of massive binary formation. In neither of our simulations do we see indications for massive binaries so far which of course could be a consequence of the early stage the simulations end. Nevertheless, our results indicate that massive binaries possibly form in later evolutionary stages and that the initial mass ratio should be far from unity with the more massive star sitting in the centre. Observations, however, reveal that a significant fraction of massive binaries consists of components with nearly equal masses (e.g Mason et al., 1998; Pinsonneault & Stanek, 2006). This apparent contradiction could be resolved by the work of Artymowicz & Lubow (1996) and Bate & Bonnell (1997). These authors simulating circumbinary discs propose that the masses of binary components, even when unequal in the beginning, tend to equal as accretion within the disc occurs preferentially onto the lower-mass companion. Initial density perturbations not consider here could also be a natural way to from massive binaries. However, there are indications that fragmentation induced by initial density perturbations might be suppressed by magnetic fields (Hennebelle & Teyssier, 2008). Thus, different binary formation scenarios might have to be at work.
4.3 Accretion rates
As shown in Fig. 12, the accretion rates of the different runs vary only by a factor of 3 which is attributed to the two competing effects of the magnetic field namely magnetic braking and the Lorentz force. Adopting accretion rates of a few 10 M yr as observed in our simulations, stars of about 30 M would form within some 10 up to a few 10 yr roughly independent of the initial conditions. However, this only holds if the accretion rates stay roughly constant over the whole formation process which might be an oversimplification as indicated by the slight decrease seen in Fig. 13 (see also Klessen, 2001; Schmeja & Klessen, 2004). A possible way to significantly change accretion rates would be to vary the initial density profile and mass of the molecular cloud core, parameters which are not explored in this work in order to limit the computational costs (see e.g. Girichidis et al., 2011).
However, our observed accretion rates agree well with accretion rates from a number of massive star formation simulations. To begin with, our accretion rates match those necessary to overcome radiation pressure deduced from 1-dimensional calculations (Kahn, 1974; Wolfire & Cassinelli, 1987). Radiation-hydrodynamical collapse simulations in 2 dimensions (Yorke & Sonnhalter, 2002; Kuiper et al., 2010) and 3 dimensions (Krumholz et al., 2007, 2009; Kuiper et al., 2011) with similar core masses reveal accretion rates of a few 10 M yr up to 10 M yr very similar to ours. Peters et al. (2010a, b, 2011) simulating the long term evolution of H ii-regions around massive protostars find similar accretion rates as well. This suggests that in our simulations accretion would continue even if the effect of radiation would be included. Girichidis et al. (2011) studying the effect of different initial conditions on massive star formation find accretion rates around 10 M yr, somewhat higher than ours possibly caused by the omission of rotation or magnetic fields counteracting gravity and accretion. Similar accretion rates were also found by Banerjee & Pudritz (2007) studying the very early evolution of a collapsed cloud core with initial magnetisation and rotation using a significantly higher resolution than in our work. Hennebelle et al. (2011) observe accretion rates of the order of 10 - 10 M yr somewhat smaller than ours which we attribute to the larger core size of 1 pc used in their setup.
Our accretion rates also agree quite well with theoretical estimates and observations. Calculating the accretion rates with the formula given in the theoretical work of McKee & Tan (2003) adapted to our setup we find a value of about 3 10 M yr very similar to the actually observed accretion rates. Observational results for accretion rates in massive star forming regions are hard to obtain and are often calculated only indirectly via observed outflow mass rates. Nevertheless, results from several high-mass star forming regions indicate accretion rates of the order of 10 - 10 M yr (e.g. Beuther et al., 2002a, c; Beltrán et al., 2006) again in good accordance with our results.
As shown in Fig. 13, the accretion rates do not drop significantly over time (except for runs with , see Section 4.1). Hence, the outflows launched seem to be not capable of significantly reducing mass accretion onto the sinks over time. In order to analyse to what extent the combined effect of magnetic fields and rotation influences mass accretion, we performed the reference calculation inf-0 with no magnetic field and zero rotation. In this run the sink particle accretes 3.55 M within the first 4000 yr corresponding to a time averaged accretion rate of M yr. Hence, the accretion rates in runs with magnetic fields and rotation are reduced to a level of about 35% - 95% of this value (see Table 2). Our simulations show that feedback due to outflows is not able to stop accretion onto the protostars and hence outflows are not able to determine the low star formation efficiency observed in molecular clouds.
We have studied the collapse of massive molecular cloud cores with varying initial rotational and magnetic energies. The cores are supercritical with mass-to-flux ratios between 2.6 and 26 and have rotational energies well below the gravitational energy. Containing about 56 Jeans masses the cores are highly gravitationally unstable and hence might be sites of protostellar cluster formation. We focussed our discussion on the formation of protostellar discs and on protostellar accretion as measured by sink particles. We find that disc properties are highly sensitive to the initial magnetic field strength whereas protostellar accretion rates are only marginally influenced by varying initial conditions. In the following we summarise our main findings.
For normalised mass-to-flux ratios below 10 the formation of centrifugally supported discs is completely suppressed. Instead, for sub-Keplerian discs are formed showing nearly no rotational support but radial infall velocities close to free-fall. For weak magnetic fields (), however, well defined Keplerian discs with sizes of a few 100 AU build up over time. This agrees well with several studies of collapsing low-mass cores.
Sub-Keplerian rotation for strong magnetic fields () is caused by magnetic braking. Analysing the torques acting on the midplane we find that angular momentum is removed due to magnetic braking at roughly the same rate as it is added due to the gas infall thereby preventing Keplerian discs from forming.
Observed accretion rates are of the order of a few 10 M yr varying only within a factor of 3 between the individual runs. This variation is remarkably small considering the large differences in the initial conditions varying over more than two orders of magnitude. We attributed this to two competing effects of the magnetic field. Increasing the magnetic field strength results in an increased accretion rate due to an enhanced magnetic braking efficiency lowering centrifugal support. Above a certain field strength, however, a further increase leads to declining accretion rates due to the effect of magnetic pressure and tension. This results in rather moderate changes in the accretion rates for different initial conditions. Furthermore, accretion rates for different amounts of angular momentum converge with increasing field strength due to the effect of magnetic braking.
For the majority of the simulations disc fragmentation does not occur. Analysing the Toomre parameter in these discs reveals that fragmentation is mainly suppressed due to the magnetic pressure. Thermal pressure alone cannot account for the stabilisation as is well below 1 whereas the magnetic Toomre parameter settles around 1 indicating stability. In two simulations, which are subject to fragmentation, neither thermal nor magnetic pressure can stabilise the discs. Accordingly, both and are well below 1.
In the two runs with disc fragmentation more than 10 sink particles are formed during the first 4000 yr. Among these sinks only the first one created reaches a mass above 1 M thus containing more than 80% of the totally accreted mass. All other particles have masses well below 0.1 M. The discs formed in both cases reach masses around 1 M somewhat below the totally accreted sink particle masses.
The outflows launched from the protostellar discs are not capable of significantly reducing mass accretion over time. Compared to a run with zero magnetic field and zero rotation, mass accretion is reduced to a level of 35% - 95%.
Radial profiles of column density and temperature exhibit accretion shock features moving outwards as time evolves. The shocks occur when the infalling gas hits either centrifugal or magnetic barriers thus strongly decreasing the infall speed.
A growing number of observations of discs and bipolar outflows around high-mass protostars (see Beuther & Shepherd, 2005; Cesaroni et al., 2007, for recent reviews) support a high-mass star formation scenario via disc accretion. On the other hand, observations also reveal that prestellar cores with masses ranging from 2 - 2000 M are usually only slightly supercritical with (Falgarone et al., 2008; Girart et al., 2009; Beuther et al., 2010). Together with our numerical results this suggest that there should be no Keplerian discs in the very early stages of typical high-mass star forming regions but rather flattened, strongly sub-Keplerian structures. This apparent dichotomy has an important impact on the formation of discs around massive stars. To enable the observed presence of centrifugally supported discs in later stages, effects in the later evolution of the system are required reducing the efficiency of magnetic braking. We discussed possible effects in the context of massive star formation like outflows and non-ideal MHD leading to a successive growth of discs in the later evolution.
So far, the influence of turbulence is completely neglected in our setup. Hence, we plan to include initial velocity perturbations in our setup to obtain more realistic initial conditions as indicated by observations (e.g. Caselli & Myers, 1995). For subsequent simulation it would be also interesting to study the effect of ambipolar diffusion which is expected to act efficiently in the high density regime of protostellar discs. For this purpose, the existing set of simulations serves as a useful guide to select representative simulations to be repeated with increased resolution and additional physics.
The authors like to thank the anonymous referee for comments which helped to significantly improve the paper. D.S. likes to thank Philipp Girichidis and Christoph Federrath for many useful discussions and suggestions. The simulations presented here were performed on HLRB2 at the Leibniz Supercomputing Centre in Garching and on JUROPA at the Supercomputing Centre in Jülich. The FLASH code was developed partly by the DOE-supported Alliances Center for Astrophysical Thermonuclear Flashes (ASC) at the University of Chicago. D.S. and R.B. acknowledge funding of Emmy-Noether grant BA3706 by the DFG. R.S.K. acknowledges subsidies from the Baden-Württemberg-Stiftung (grant P-LS-SPII/18) and from the German Bundesministerium für Bildung und Forschung via the ASTRONET project STAR FORMAT (grant 05A09VHA). R.S.K. furthermore gives thanks for subsidies from the Deutsche Forschungsgemeinschaft (DFG) under grants no. KL 1358/1, KL 1358/4, KL 1359/5, KL 1358/10, and KL 1358/11, as well as from a Frontier grant of Heidelberg University sponsored by the German Excellence Initiative. D.D. is supported by McMaster University and R.E.P by a Discovery grant from NSERC of Canada.
Here we compare simulation 26-4 to two more runs with identical initial conditions but a maximum spatial resolution varied by a factor 4 in either direction. In particular, the initial resolution is identical to that in run 26-4, i.e. in the beginning the mesh in the core has a spacing of 302 AU. We list the runs and their corresponding parameters in Table 3. The critical value of the density above which sink particles are created is adapted in accordance with the resolution. For all runs performed the refinement criterion applied guaranties that the midplane region is resolved on the highest level used. In particular, we focus in this analysis on accretion properties and radial profiles of different quantities in the midplane. Due to computational cost reasons run 26-4-H is followed for 2000 yr only. Hence, we compare the results at this time.
|(AU)||(10 g cm)||(kyr)||(M)|
First, we consider radial profiles of the density, temperature and velocity in the midplane.
The quantities in each run are taken in the midplane. In the density profile (left panel of Fig. 16) the accretion shock occurring at 150 AU is clearly resolved in the runs 26-4-H and 26-4. In run 26-4-L, however, the shock is somewhat smoothed out due to the limited resolution. Hence, a resolution of 4.7 AU seems required to properly resolve the accretion shock. Within the accretion shock, however, the density increases with resolution. We attribute this to the fact that the vertical structure of the disc is not fully resolved at least in the runs 26-4 and 26-4-L. Here in large parts the vertical disc height is represented by a few grid points only (compare Fig. 9). Therefore, to fully resolve the vertical disc structure, a higher resolution, probably even above that in run 26-4-H, would be needed which is not feasible for computational costs reasons.
A similar result holds for the temperature profiles as well (middle panel of Fig. 16). The temperature jump at the accretion shock seems be to reasonably well resolved in run 26-4 whereas in run 26-4-L it is smoothed out markedly. Within the accretion shock the temperature reaches as higher values as higher the final resolution is. This is due to the strong coupling of temperature and density above 10 g cm where the gas gets optically thick resulting in higher temperatures at higher gas densities.
Next, we analyse the velocity structure in the midplane (see right panel of Fig. 16). In run 26-4-L no Keplerian disc has built up yet. However, we mention that in its further evolution the rotation reaches Keplerian velocities as well. In contrast, in run 26-4 the rotation is already Keplerian up to a radius of 60 AU in good agreement with run 26-4-H. Furthermore, the comparison shows that in run 26-4 the decline in / at 20 AU is most likely a resolution effect as in run 26-4-H this decline occurs at a roughly three times smaller radius of 8 AU. This supports the statement made in Section 3.1 that the inner 10 AU in runs with a resolution of 4.7 AU are strongly affected by numerical resolution (see also Fig. 7). The radial velocities at radii larger 10 AU show a reasonably well agreement between all runs. Hence, we conclude that regarding the velocity structure in the midplane run 26-4 is reasonably well converged.
In Fig. 17 we show the time evolution of the total accretion rate for the three runs considered. In run 26-4-H two more sink particles are created after roughly 1500 yr which cause the large variations in the accretion rate. In general, however, the accretion rates of the three runs are of the same order of magnitude. From the formation time t listed in Table 3 it can be inferred that t increases with spatial resolution. This behaviour is expected as for higher resolution sink particles are created at higher densities and thus later times during the collapse. It can also be inferred from Table 3 that the mean accretion rate, i.e. the total accreted mass, decreases with increasing spatial resolution. The accretion rate of run 26-4-L is higher than that of run 26-4 by about 35%. We therefore conclude that regarding accretion properties run 26-4-L is not yet converged. The difference in accreted mass between run 26-4-H and run 26-4 is of the order of 2% thus significantly lower than the difference between run 26-4-L and run 26-4. Hence, there seems to be a convergence of the accretion rates with increasing resolution and, although run 26-4 seems not fully converged, we conclude that a resolution of 4.7 AU is sufficiently high to properly describe accretion properties of the protostars.
In summary, we can say that a resolution of 4.7 AU seems sufficiently high to properly follow protostellar accretion rates, resolve the accretion shock at the edge of the disc, and display the velocity structure in the disc down to about 10 - 15 AU. In contrast, run 26-4-L with a resolution of 18.9 AU reveals significant differences from run 26-4-H showing that is not yet converged. However, also run 26-4 seems to be not fully converged regarding to the density and temperature structure in the disc. We attribute this partly to the poor resolution of the vertical disc structure. Hence, a higher resolution, probably even above the one used in run 26-4-H, would be necessary to reach convergence regarding this point. However, due to significantly higher computational costs this has not been feasible wherefore a spatial resolution of 4.7 AU is used throughout the paper. This particular choice is also motivated physically as we want to resolve the first core. Hence, the Jeans length at densities above 10 g cm has to be resolved requiring a resolution of a few AU.
- pagerange: Magnetic fields during the early stages of massive star formation – I. Accretion and disc evolution–Appendix
- pubyear: 2011
- To avoid unphysically high densities in the interior of the core, we cut off the -profile at a radius of 0.0125 pc. Within this radius the density distribution follows a parabola .
- Allen, A., Li, Z., & Shu, F. H. 2003, ApJ, 599, 363
- Alves, F. O., Girart, J. M., Lai, S.-P., Rao, R., & Zhang, Q. 2011, ApJ, 726, 63
- Andre, P., Ward-Thompson, D., & Barsony, M. 2000, Protostars and Planets IV, 59
- Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, Protostars and Planets V, 245
- Artymowicz, P. & Lubow, S. H. 1996, ApJ, 467, L77+
- Banerjee, R. & Pudritz, R. E. 2006, ApJ, 641, 949
- Banerjee, R. & Pudritz, R. E. 2007, ApJ, 660, 479
- Banerjee, R., Pudritz, R. E., & Anderson, D. W. 2006, MNRAS, 373, 1091
- Baraffe, I. & Chabrier, G. 2010, A&A, 521, A44+
- Bate, M. R. & Bonnell, I. A. 1997, MNRAS, 285, 33
- Beck, R. 2009, Astrophysics and Space Sciences Transactions, 5, 43
- Beltrán, M. T., Cesaroni, R., Codella, C., et al. 2006, Nature, 443, 427
- Beuther, H., Churchwell, E. B., McKee, C. F., & Tan, J. C. 2007, Protostars and Planets V, 165
- Beuther, H., Schilke, P., Gueth, F., et al. 2002a, A&A, 387, 931
- Beuther, H., Schilke, P., Menten, K. M., et al. 2002b, ApJ, 566, 945
- Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002c, A&A, 383, 892
- Beuther, H. & Shepherd, D. 2005, in Cores to Clusters: Star Formation with Next Generation Telescopes, ed. M. S. N. Kumar, M. Tafalla, & P. Caselli, 105–119
- Beuther, H., Vlemmings, W. H. T., Rao, R., & van der Tak, F. F. S. 2010, ApJ, 724, L113
- Bouchut, F., Klingenberg, C., & Waagan, K. 2007, Numerische Mathematik, 108, 7, 10.1007/s00211-007-0108-8
- Bouchut, F., Klingenberg, C., & Waagan, K. 2010, Numerische Mathematik, 115, 647, 10.1007/s00211-010-0289-4
- Caselli, P. & Myers, P. C. 1995, ApJ, 446, 665
- Cesaroni, R., Galli, D., Lodato, G., Walmsley, C. M., & Zhang, Q. 2007, Protostars and Planets V, 197
- Curran, R. L. & Chrysostomou, A. 2007, MNRAS, 382, 699
- Dapp, W. B. & Basu, S. 2010, A&A, 521, L56+
- Duffin, D. F. & Pudritz, R. E. 2009, ApJ, 706, L46
- Duquennoy, A. & Mayor, M. 1991, A&A, 248, 485
- Falgarone, E., Troland, T. H., Crutcher, R. M., & Paubert, G. 2008, A&A, 487, 247
- Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010, ApJ, 713, 269
- Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
- Fuller, G. A., Zijlstra, A. A., & Williams, S. J. 2001, ApJ, 555, L125
- Girart, J. M., Beltrán, M. T., Zhang, Q., Rao, R., & Estalella, R. 2009, Science, 324, 1408
- Girichidis, P., Federrath, C., Banerjee, R., & Klessen, R. S. 2011, MNRAS, 413, 2741
- Goldreich, P. & Lynden-Bell, D. 1965, MNRAS, 130, 125
- Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528
- Hennebelle, P. & Ciardi, A. 2009, A&A, 506, L29
- Hennebelle, P., Commerçon, B., Joos, M., et al. 2011, A&A, 528, A72+
- Hennebelle, P. & Fromang, S. 2008, A&A, 477, 9
- Hennebelle, P. & Teyssier, R. 2008, A&A, 477, 25
- Hosking, J. G. & Whitworth, A. P. 2004, MNRAS, 347, 1001
- Kahn, F. D. 1974, A&A, 37, 149
- Kim, W. & Ostriker, E. C. 2001, ApJ, 559, 70
- Klessen, R. S. 2001, ApJ, 550, L77
- Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585
- Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959
- Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754
- Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010, ApJ, 722, 1556
- Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20
- Lai, S., Crutcher, R. M., Girart, J. M., & Rao, R. 2001, ApJ, 561, 864
- Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ArXiv e-prints
- Mac Low, M.-M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
- Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2011, PASJ, 63, 555
- Machida, M. N., Matsumoto, T., Hanawa, T., & Tomisaka, K. 2005, MNRAS, 362, 382
- Mason, B. D., Gies, D. R., Hartkopf, W. I., et al. 1998, AJ, 115, 821
- McKee, C. F. & Ostriker, E. C. 2007, ARA&A, 45, 565
- McKee, C. F. & Tan, J. C. 2003, ApJ, 585, 850
- Mellon, R. R. & Li, Z. 2008, ApJ, 681, 1356
- Mellon, R. R. & Li, Z. 2009, ApJ, 698, 922
- Mouschovias, T. C. & Paleologou, E. V. 1980, ApJ, 237, 877
- Mouschovias, T. C. & Spitzer, Jr., L. 1976, ApJ, 210, 326
- Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199
- Padoan, P., Nordlund, Å., Rögnvaldsson, Ö. E., & Goodman, A. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 243, From Darkness to Light: Origin and Evolution of Young Stellar Clusters, ed. T. Montmerle & P. André, 279–+
- Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M. 2011, ApJ, 729, 72
- Peters, T., Banerjee, R., Klessen, R. S., et al. 2010a, ApJ, 711, 1017
- Peters, T., Klessen, R. S., Mac Low, M., & Banerjee, R. 2010b, ApJ, 725, 134
- Pinsonneault, M. H. & Stanek, K. Z. 2006, ApJ, 639, L67
- Pirogov, L., Zinchenko, I., Caselli, P., Johansson, L. E. B., & Myers, P. C. 2003, A&A, 405, 639
- Pirogov, L. E. 2009, Astronomy Reports, 53, 1127
- Price, D. J. & Bate, M. R. 2007, MNRAS, 377, 77
- Schmeja, S. & Klessen, R. S. 2004, A&A, 419, 405
- Shepherd, D. S., Claussen, M. J., & Kurtz, S. E. 2001, Science, 292, 1513
- Tilley, D. A. & Pudritz, R. E. 2007, MNRAS, 382, 73
- Toomre, A. 1964, ApJ, 139, 1217
- Waagan, K. 2009, Journal of Computational Physics, 228, 8609
- Waagan, K., Federrath, C., & Klingenberg, C. 2011, Journal of Computational Physics, 230, 3331
- Wolfire, M. G. & Cassinelli, J. P. 1987, ApJ, 319, 850
- Wuchterl, G. & Klessen, R. S. 2001, ApJ, 560, L185
- Yorke, H. W. & Sonnhalter, C. 2002, ApJ, 569, 846
- Zinnecker, H. & Yorke, H. W. 2007, ARA&A, 45, 481