In situ formation of SgrA stars via disk fragmentation: parent cloud properties and thermodynamics
The formation of the massive young stars surrounding SgrA is still an open question. In this paper, we simulate the infall of a turbulent molecular cloud towards the Galactic Center (GC). We adopt two different cloud masses ( M and M). We run five simulations: the gas is assumed to be isothermal in four runs, whereas radiative cooling is included in the fifth run. In all the simulations, the molecular cloud is tidally disrupted, spirals towards the GC, and forms a small, dense and eccentric disk around SgrA. With high resolution simulations, we follow the fragmentation of the gaseous disk. Star candidates form in a ring at pc from the super-massive black hole (SMBH) and have moderately eccentric orbits (), in good agreement with the observations. The mass function of star candidates is top-heavy only if the local gas temperature is high ( K) during the star formation and if the parent cloud is sufficiently massive ( M). Thus, this study indicates that the infall of a massive molecular cloud is a viable scenario for the formation of massive stars around SgrA, provided that the gas temperature is kept sufficiently high ( K).
Subject headings:Methods: numerical - Galaxy : center - Stars: formation - ISM: clouds
The origin of young massive stars that crowd the Galactic Center (GC) has been a puzzle for a long time. Most of the massive stars observed in the central parsec reside in one or perhaps two disks: a well-defined clockwise-rotating disk, and a counter-clockwise disk, whose existence is still debated (Genzel et al. 2003, hereafter G03; Paumard et al. 2006, hereafter P06; Lu et al. 2009; Bartko et al. 2009, 2010; see Genzel et al. 2010 for a recent review). These disks have well-defined inner ( pc) and outer radii ( pc). The stars observed at distances pc from SgrA, the source identified with the super massive black hole (SMBH), have randomly oriented motions and do not belong to the disks (G03; Ghez et al. 2005; Eisenhauer et al. 2005). The massive stars inside the disks are young ( Myr, P06) and must have formed over a short period ( Myr, P06). Their estimated initial mass function (IMF) is heavier than Salpeter’s one (P06; Bartko et al. 2009). The total mass in the disks cannot exceed M, but is more likely of the order of M (P06).
Such stars cannot have formed in situ in ‘normal’ conditions, as the tidal forces exerted from the SMBH would have disrupted the parent molecular cloud (Levin & Beloborodov 2003; G03). Thus, an alternative scenario has been proposed, according to which a young cluster spiraled towards the GC and deposited its stars around SgrA (Gerhard 2001; McMillan & Portegies Zwart 2003; Portegies Zwart, McMillan & Gerhard 2003; Kim & Morris 2003; Kim, Figer & Morris 2004; Gürkan & Rasio 2005; Fujii et al. 2008, 2009). However, even the latter scenario suffers from various shortcomings, such as the premature disruption of the cluster and the excessively long dynamical friction time. These problems have only been partially solved by assuming that the original clusters host intermediate-mass black holes (Portegies Zwart et al. 2006), or by using new computational schemes (Fujii et al. 2008).
On the other hand, the problem of tidal forces exerted by the SMBH can be overcome if, at some point in the past, a dense gaseous disk existed around SgrA. Such a disk could have formed because of the infall and tidal disruption of a molecular cloud. If the density in the disk was high enough, it might have become unstable to fragmentation and formed stars (Poliachenko & Shukhman 1977; Kolykhalov & Sunyaev 1980; Shlosman & Begelman 1987; Sanders 1998; Levin & Beloborodov 2003; G03; Goodman 2003; Milosavljevic & Loeb 2004; Nayakshin & Cuadra 2005; Rice et al. 2005; Alexander et al. 2008; Collin & Zahn 2008; Wardle & Yusef-Zadeh 2008; Yusef-Zadeh & Wardle 2011). The (nearly) absence of massive stars at distances pc from SgrA supports this idea of in situ formation (Nayakshin & Sunyaev 2005; P06). This scenario is also favored by the existence of two massive molecular clouds within pc from the dynamical center of our Galaxy (Solomon et al. 1972). On-going star formation (SF) is associated with one of these two clouds, named M (Yusef-Zadeh et al. 2010). The other cloud (M) is highly elongated toward SgrA and has a ‘finger-like’ extension pointing in the direction of the circumnuclear ring (Okumura et al. 1991; Ho et al. 1991; Novak et al. 2000, and references therein). The circumnuclear ring of dense molecular gas lies at a projected distance of pc from the GC (e.g., Genzel et al. 1985; Güsten et al. 1987; Scoville et al. 2003; Christopher et al. 2005) and is on the verge of forming stars (Yusef-Zadeh et a. 2008; 2011). Furthermore, the innermost 3 pc of the GC are extremely rich of ionized gas, organized in at least three streams (e.g., Scoville et al. 2003; Zhao et al. 2009, and references therein). Recently, a small ( Earth masses) and relatively cold cloud (gas temperature K) has been observed on its way towards the SMBH, with a orbital plane coincident with that of the clockwise stellar disk (Gillessen et al. 2012; Burkert et al. 2012).
Nayakshin, Cuadra & Springel (2007, hereafter NCS07) and Alexander et al. (2008) simulated SF in a gaseous disk around SgrA, and found encouraging results for this scenario. However, NCS07 and Alexander et al. (2008) assume that the gaseous disk was already in place when it started forming stars, and do not consider the process that lead to the formation of the disk itself. Bonnell & Rice (2008, hereafter BR08), Mapelli et al. (2008, hereafter M08) and Alig et al. (2011, hereafter A11) simulate the infall of a molecular cloud toward SgrA and study the formation of a dense gaseous disk around the SMBH. Hobbs & Nayakshin (2009, hereafter HN09) simulate the collision of two molecular clouds leading to the formation of one or more disks around the SMBH. A11 study the density distribution of the gas in the circumnuclear disk, but do not trace the evolution of the system up to the fragmentation and to the formation of the proto-stars. BR08, M08 and HN09 use the sink-particle technique to model SF in the gaseous disk around the SMBH. This was required by the huge density range between the ‘background’ gas and the particles that are collapsing into stars. However, the sink-particle technique might be affected by various uncertainties, because it does not trace the fragmentation of the cloud, but it replaces this process with the formation of sink particles. In this paper, we simulate the infall of a molecular cloud toward SgrA and we study the ongoing SF around the SMBH without adopting the sink-particle method. Our results indicate that the possibility of forming massive stars around SgrA strongly depends on the initial mass of the infalling cloud and on the temperatures that are locally reached by the gas.
2. Models and simulations
We ran N-body/Smoothed Particle Hydrodynamics (SPH) simulations of a molecular cloud evolving in a potential dominated by the SMBH. We used the SPH code GASOLINE (Wadsley, Stadel & Quinn 2004), upgraded with the Read et al. (2010) OSPH modifications, to address the SPH limitations outlined, most recently, by Agertz et al. (2007).
In the simulations, the SMBH is represented by a sink particle, with initial mass M (Ghez et al. 2003), sink radius pc and softening radius pc. We also add a rigid potential, according to a density distribution , where M, pc, and at and at (G03; although recent papers indicate a more complex scenario, e.g., Yusef-Zadeh, Bushouse & Wardle 2012). Before this paper, the cuspy rigid potential was introduced only by M08 and by HN09, whereas other papers (e.g., BR08; A11) consider the SMBH and the cloud in isolation.
The clouds used in this experiment are spherical with a radius of 15 pc, are seeded with supersonic turbulent velocities and marginally self-bound. To simulate interstellar turbulence, the velocity field of the cloud was generated on a grid as a divergence-free Gaussian random field with an imposed power spectrum P(), varying as . This yields a velocity dispersion , varying as , chosen to agree with the Larson scaling relations (Larson 1981). The velocities were then interpolated from the grid to the particles.
Note. – Notes. is the number of gas particles in each run, M and T are the total mass and the initial temperature of the simulated molecular cloud, respectively. E.O.S. indicates the equation of state of gas.
The center-of-mass of the cloud is initially at 25 pc from the SMBH. The orbit of the cloud was chosen so that the impact parameter111The adopted impact parameter is very small, to guarantee a low orbital angular momentum and to increase the chance that the core of the parent cloud engulfs the SMBH (Wardle & Yusef-Zadeh 2008). This choice is biased, as it likely favors the formation a dense gas disk. The properties of the resulting gas disk are expected to depend on the impact parameter, as well as on the other (both orbital and intrinsic) properties of the cloud. On the other hand, BR08 adopt a 10 times larger value for the impact parameter (0.1 pc) than in this paper, and the differences in the properties of the gaseous disk (e.g., scale-length, thickness, outer and inner radius) are not significant. The role of the impact parameter as well as that of the other orbital properties deserve further studies. with respect to the SMBH is pc and the initial velocity is close to the escape velocity from the SMBH at the initial distance (i.e. the orbit is marginally bound and highly eccentric). As the tidal density at pc from the GC is higher than the average initial density of our cloud, our marginally bound cloud must have formed further out and then migrated closer to the GC. Various processes could have brought the cloud in such position. For example, the cloud might have achieved this orbit after a collision with another cloud (e.g., Wardle & Yusef-Zadeh 2008; HN09; A11). On the other hand, the existence of two giant molecular clouds, M and M, at and pc, respectively, from the GC (Solomon et al. 1972; Okumura et al. 1991; Ho et al. 1991; Novak et al. 2000) shows that molecular clouds, probably unbound, exist at pc from the GC.
The cloud is assumed to be isothermal in runs AD. This assumption neglects the local variations of the effective equation of state, that might occur since the balance between heating and cooling processes depends on local conditions (see Section 3.3). In run E, we include radiative cooling. The radiative cooling algorithm is the same as that described in Boley (2009) and in Boley et al. (2010). The cooling is calculated from , where and , for the local opacity , particle mass and density . D’Alessio et al. (2001) opacities are used, with a 1 m maximum grain size. The irradiation temperature is T K everywhere.
In this paper, we run five simulations, with two different masses for the cloud ( and M) and two different initial gas temperatures (100 and 500 K). The initial temperature of the gas is set to be quite high for a Galactic molecular cloud, as observational data (see e.g. Nagai et al. 2007; Oka et al. 2007) suggest that the gas temperature is K in the vicinity of Sgr A. The initial masses of the simulated clouds, the initial temperature of the gas and the adopted equation of state in the various runs are listed in Table 1. In all the runs, the mass of the gas particles is 0.04 M and the softening length pc.
3.1. Evolution of the gaseous disk
In all the runs, the cloud, orbiting around SgrA, rapidly ( yr) stretches towards the SMBH and is partially disrupted. The branch of the cloud that points toward the SMBH begins to spiral in towards the GC (see Fig. 1, in the case of run A). A dense, rotating gaseous disk forms at yr at the location of the in-spiraling branch. In runs A, B and E (C and D), the initial density222Here and in the following, we assume molecular weight , the same as in our simulations, consistently with Galactic molecular clouds. of the gaseous disk is cm ( cm) and its initial mass is (). The outer radius of the disk is pc. At this stage, the parent cloud is still feeding it through a finger-like structure (see Fig. 1).
The gaseous disk is not homogeneous or regular, but it is the assembly of many concentric annuli, that spiral around the SMBH, with slightly different inclination and thickness (see Fig. 2). However, all the structures that form inside pc lie almost on the same plane as the disk: we do not see the formation of any relevant structure with significant inclination () in the very central region.
The initial average thickness of the disk, defined as the ratio between the average scaleheight and , is .
At yr the disk has the average density of cm in runs A, B and E, and of cm in runs C and D (see Fig. 2). However, local densities cm are reached (see the density map of gas in Fig. 2). At this stage, the total mass of the gaseous disk is in runs A, B and E, and in runs C and D.
The outer radius and the thickness of the disk are still pc and : they do not change during the entire simulation. The size of the gaseous disk agrees with the estimate by Wardle & Yusef-Zadeh (2008, 2011), based on the Bondi-Hoyle-Lyttleton formalism (Hoyle & Lyttleton 1939) for the capture radius of a black hole. Applying this formalism to our simulated clouds, we predict a disk radius pc (where is the velocity of the cloud when entering the influence radius of the SMBH and is the average ratio of the specific angular momentum of a fluid element settled in the final disk to its initial angular momentum, see equation 5 of Wardle & Yusef-Zadeh 2011 for details). The disk appears distorted at the edges, where fresh gas is being fed by the parent cloud. The orbits of gas particles have eccentricity in the range .
3.2. Formation of the stars
In our analysis, we assume that a star (or, more precisely, a stellar candidate) is formed, when a spherical clump of particles (where 32 is the number of neighbors in our code) and radius pc reaches an average density cm (assuming molecular weight ), i.e. it has an average density times higher than the average density of the disk. Under such conditions, the spherical clump is gravitationally bound and completely decoupled from the surrounding particles. The threshold values and have been selected on the basis of the density maps (see Fig. 2). We notice that the threshold value cm is well above the maximum tidal density inside the simulated disk, which is cm. accounts for the behavior of the Jeans instability (Jeans 1919). There are no significant changes in the stellar mass function (MF) for , where is the softening, and for . Furthermore, our choice of the threshold values is quite conservative, as it allows to exclude marginally bound objects. For further details about the choice of and , see the discussion in Appendix A.
3.2.1 The stellar disk
The first gravitationally bound clumps, i.e. the first stellar candidates form at yr. We continue the integration up to yr, in order to account for possible mass accretion of the proto-star and to allow most of gravitational instabilities in the disk to evolve. For yr the integration slows down dramatically in run A, B and E, because of the large number of collapsed clumps.
In Fig. 3 the star candidates are marked by filled circles, superimposed to the density map. Run C is the only simulation in which no stars have formed yet at yr. In run A, B, D and E the star candidates are distributed in a ring, corotating with the gaseous disk. The stellar ring has a thickness . In run A we can almost distinguish two concentric rings. In runs A, B, D and E the ring (and therefore the orbits of star candidates) is eccentric, with minimum (maximum) eccentricity (). The average eccentricity is . Fig. 4 shows the eccentricity versus the semi-major axis for all the stars in run E. The eccentricity and the semi-major axis were derived assuming that each star was in a binary system with the SMBH and neglecting the perturbations from other stars and gas. We notice that the maximum eccentricity tends to increase with the semi-major axis: goes from for pc to for pc. The simulated eccentricity range of the disk stars () is in good agreement with the observed one for the clockwise disk surrounding SgrA (; P06; Bartko et al. 2009). However, we stress that our simulations cover only the first Myr since the formation of the gaseous disk, whereas the estimated age of the stars surrounding SgrA is Myr: the long-term dynamical evolution of the stellar disk might have affected the orbital properties of the stars.
Fig. 4 shows that, in our run E, there is a gap in the distribution of semi-major axes between pc and pc. Furthermore, only a handful of stars have pc. This depends on the structure of the gaseous disk, which is composed of many concentric annuli (as we said in Section 3.1). In particular, stars with pc in run E come essentially from the same gaseous annulus, which has an outher radius pc (see Fig. 3). Out of this annulus, there is a lower density region, surrounded by an outer high-density annulus (Fig. 3). The bunch of stars with pc formed in this outer annulus.
In runs A, B, D and E, the first stars form in the outer parts of the ring, whereas stars form in the inner part (where the shear from the SMBH is stronger) only at later times. Furthermore, the clumps of gas particles tend to collapse while at the apocenter, where the distance from the SMBH is maximum. Figs. 2 and 3 show that run B and run E are very similar.
Fig. 5 shows the radial distribution of the stars that are closest to the pericenter, i.e. in the quadrant with and (top panel), and the radial distribution of the stars that are closest to the apocenter, i.e. in the quadrant with and (bottom panel), at yr. The top panel of this Figure shows that the stars at the pericenter are at a minimum distance from the SMBH pc in run A, B and E, and pc in run D. The bottom panel of the same Figure indicates that the stars at the apocenter are at a maximum distance from the SMBH pc in all the simulations.
In the simulations, the existence of the inner radius is probably due to the fact that even the high central density of gas cannot counteract the Keplerian velocity at such small distances from the SMBH. In fact, Toomre’s parameter, defined as (where is the angular velocity, the sound speed, the gravitational constant and is the local surface density, Toomre 1964), is at radii , , and pc for run A, B, D and E, respectively (see Fig. 6). For higher values of , the growth of gravitational instabilities in the disk is unlikely. The behavior of Toomre’s parameter also explains why gas clumps collapse predominantly at the apocenter of the orbit (see Fig. 6). Furthermore, we notice that, in Fig. 6, is always larger than 1 for run C: this is consistent with the fact that no star candidates form in run C. In addition, Fig. 6 indicates that pc is a sort of cut-off for the SF in our simulations, as for pc in all the runs. On the other hand, this result strongly depends on the initial properties of the simulated cloud (density, mass, orbit, etc.). Moreover, pc should be interpreted as the tail end of a continuous distribution rather than a genuine cut-off.
The gas inside the inner radius remains in place till the end of the simulations. It would be interesting to assess whether this gas can be evaporated by the ultra-violet (UV) radiation and by the first supernovae of the young massive stars. Recently, Alexander et al. (2011) suggested that a fraction of the innermost gas (inside pc) might have been used to power Sgr A in the past.
We notice that the outer radius ( pc) and the thickness () of the stellar ring are slightly smaller than those of the gaseous disk (for which pc and pc), but the stellar ring is inside the gaseous disk, has the same orientation, and the stellar orbits have the same eccentricity as those of the gaseous particles. Similarly to the gaseous disk, the stellar disk is irregular, quite warped and tilted. To quantify the importance of the warp and/or of the tilt, we study the radial change of the angular momentum orientation of stars and gas particles. In particular, we define the angle as where is the modulus of the total angular momentum of a star or gas particle and is the component of the angular momentum along the axis which maximizes the symmetry of the gas disk inside 1 pc radius (i.e., if the disk was perfectly axisymmetric, would be along the symmetry axis of the disk). Fig. 7 shows that there is a radial change by in in both the gas and the stellar component, in all the considered runs (in the Figure, we show only runs A and E, but the other three runs have very similar behavior). This indicates that both the gaseous and the stellar disks are slightly warped/tilted. Fig. 7 cannot be directly compared with figure 13 of Bartko et al. (2009), as the angle shown therein is defined in a slightly different way (on the basis of the analysis described in Bartko et al. 2009). On the other hand, we notice that the trend of the angle in our paper and that shown in figure 13 of Bartko et al. (2009) are similar.
The ring distribution of the star candidates is in good agreement with the observations (G03; P06). The outer radius of the stellar ring is also in agreement with the observations (G03; P06), whereas the inner radius is a factor of larger in the simulations than in the observations (which suggest pc). The fact that the stellar ring is irregular and that the stellar candidates have eccentric orbits with a large spread of eccentricities are also consistent with the observations (P06; Cuadra, Armitage & Alexander 2008; Bartko et al. 2009).
3.2.2 The mass function (MF)
Fig. 8 shows the MF of the stellar candidates for the different runs at yr. Although we cannot exclude that the MF slightly evolves for yr, Fig. 8 is particularly interesting as it shows the differences among the various runs. Furthermore, the MF in Fig. 8 slightly depends on our definition of star candidates and on the adopted threshold values. Thus, the normalization of the masses in Fig. 8 might change for different definitions of the proto-stars, but the differences among various runs are independent of such normalization.
Fig. 8 indicates that the MF is heavier in run A than in run B, D and E. Fig. 8 also shows that the number of star candidates formed in runs B (1363) and E (1252) is much larger than in runs A (238) and D (49). However, the total mass of stars in run A, , is times larger than that in run D () and only a factor of lower than that in run B () and in run E (), since the MF in run A is heavier than in the other runs. We note that the total mass of stars in runs A, B and E matches the value predicted by the observations (), although we have already stressed the of the normalization of our simulated MFs.
Table 2 shows the least square fit of for runs A, B and E333In Table 2, we do not report run C, because no stars form in this run. In run D the number of stellar candidates is too low to distinguish among different MFs: we report the results of run D in Table 2, but the uncertainty on the fit parameters is huge.. The best fitting MF is heavier than the Salpeter MF () for all the considered runs. The MF of run E (with radiative cooling) is quite similar to that of run B (isothermal), whereas the MF of run A is significantly heavier than in both run B and E. The MF of run E (which differs from run B only for the inclusion of radiative cooling) deviates from that of run B mainly because of a high mass tail, that is not present in run B: 19 (10) stars form in run E (B) with mass M. The fact that runs A, B and E have a MF heavier than the Salpeter MF is important, as the observations of young stars in the GC indicate that they are relatively massive. We notice that the best-fitting slope of the MF in run A () is in fair agreement with the most recent observations of the stellar disk(s) in the GC (see, e.g., from fig. 3 of Bartko et al. 2010).
Intuitively, the reason why the MF in run A is heavier than in runs B and E (and the MF in all the performed runs is heavier than the Salpeter MF) is connected with the gas temperature. In fact, the Jeans mass (Jeans 1919) scales with . This implies that there should be a factor of difference in the typical mass of stars formed in a MC at K (which is a typical temperature for MCs, in absence of strong UV background) with respect to stars formed in a MC at K (which is a typical temperature in the neighborhoods of the GC, Nagai et al. 2007). Similarly, there should be a factor of difference between the typical mass of stars in run A and that of stars in runs B and E. Actually, the average star mass in run A is , only a factor of larger than the average star mass in run B (). The discrepancy between the difference expected from the Jeans instability and the difference obtained from the simulations is probably due to the resolution limit of the simulations. In fact, we cannot resolve stars with mass M.
Fig. 9 shows the radial dependence of the star mass () in the simulations at yr. The top panel of Fig. 9 shows the radial distribution of mass only for the stars closest to the pericenter. All these stars are at . However, this distribution is not relevant to understand their formation, as most of stars form while at the apocenter (see explanation in the previous section). Therefore, stars that are at the pericenter at yr have likely formed at larger radii. The bottom panel is more interesting, as it shows the distribution of stars that are at the apocenter at yr. In fact, such distribution approximately mirrors that of the stellar birthplaces. In all the simulations there is no significant correlation between the star mass and the radial distance from the SMBH. In run D, stars form at , whereas in runs A, B and E there is a wider spread (). In runs B and E there are also a few stars at larger radii than pc.
In Fig. 9 we also compare the radial mass distribution with the estimates from Toomre’s most unstable wavelength. Toomre’s most unstable wavelength, defined as (where is the local surface density and is the epicyclic frequency, Binney & Tremaine 1987) in the case of a thin gaseous disk, is the wavelength at which instability first appears, when drops below unity in a differentially rotating disk. The dot-dashed line in Fig. 9 shows the Toomre’s mass, defined as , in the case of run A. The distribution of stars in runs A, B and E is in partial agreement with the predictions from Toomre’s mass. The differences are likely due to the fact that Toomre’s mass does not account for the eccentricity of the disk, for further gas accretion and for changes in the local surface density. Moreover Toomre’s mass assumes that the disk is in equilibrium.
3.2.3 The SF efficiency
We can quantify the efficiency of SF. The surface density of stars (filled circles) and of gas (open squares) are both shown in the top panel of Fig. 10: the mass density of stars is comparable with that of gas only at pc and falls abruptly at smaller and larger radii. The SF efficiency (defined as , where and are the gas mass and the stellar mass, respectively) reaches a maximum of at pc in run E, and stays above 0.1 in the pc range. Out of this range, it drops very quickly. A possible issue may be the removal and/or evaporation of the remaining gas surrounding the massive stars. On the other hand, our code does not include radiative feedback (e.g., UV emission) from the newly formed stars and cannot describe the evolution of the gas after the formation of the proto-stars. For the same reason, we cannot quantify the importance of further gas accretion onto the formed stars, over longer timescales than yr. Furthermore, the estimated age for the young stars surrounding SgrA is about 6 Myr: at this epoch, the supernovae from the most massive stars are expected to have already taken place.
3.3. Properties of gas
Run E includes accurate recipes for gas cooling (see Section 2). Therefore, we focus on run E to analyze the local properties of gas in the disk. Fig. 11 shows the radial behavior of gas temperature (top panel), density (central) and Jeans mass (bottom). The temperature plot shows that most of gas is close to the temperature floor K. Gas temperature increases only very close to the SMBH, and at pc inside the collapsed clumps (i.e. the star candidates), because of the much higher density and of the consequent compressional heating.
The plot of the density (central panel of Fig. 11) is particularly interesting because it shows three (almost distinct) regimes. For pc, the local density in the disk is high ( cm) and almost constant, with a very slow rise towards the SMBH. For pc, the local density in most of the disk drops quite fast (, where the slope is ). Independently of this global behavior, a bunch of gas particles at pc reach very high densities ( cm): this represents the phase of star candidates.
Finally, Jeans mass (bottom panel of Fig. 11) confirms that we actually resolve most of the star candidates in run E, at least for pc. In fact, the expected Jeans mass is always M for pc. This means that we cannot resolve the small stars that might form at pc, as our mass resolution is 1.28 M. For larger radii the Jeans mass increases, as a consequence of the fast decrease in density, and becomes larger than our mass resolution. We stress that the subsample of gas particles with Jeans mass M in Fig. 11 is spurious: these gas particles are within the already collapsed clumps and the Jeans mass is no longer a valid estimator for them, as their density is already decoupled from the background density. We do not plot the behavior of density and Jeans mass in run B, because it is very similar to run E.
Fig. 12 shows the behavior of density and Jeans mass for run A. The behavior of density (top panel) is not very different from that of run E, apart from the fact that the number of collapsed clumps (star candidates) is lower. Instead, the Jeans masses are a factor of higher. This means that the Jeans mass is M for pc. We notice that the minimum mass of star candidates in run A from Fig. 8 is M, that is we resolve most of the star candidates in run A.
3.4. Comparison with previous papers
In this Section, we briefly compare our simulations and our main results with those reported in previous papers. In particular, we consider the four papers whose approach is closest to ours, i.e. M08, BR08, HN09 and A11. We first compare the methods adopted in these papers, and then their results.
3.4.1 Different methods
All these papers, as well as ours, simulate the infall of a molecular cloud towards the GC, its disruption from the SMBH, and the formation of a dense gaseous disk. The adopted softening and mass resolution are comparable in all these simulations. The simulations in M08, HN09, A11 and this paper were run for similar evolutionary times ( yr), whereas the simulations in BR08 were run only for yr. This difference is mainly due to the orbit of the simulated cloud(s).
M08, BR08 and this paper adopt models of gas cloud that are turbulently supported, while HN09 and A11 consider a simplified model of spherical and homogeneous cloud. The radius of the cloud is quite different in these studies: M08 and the current paper adopt a radius of 15 pc, consistent with the radii of Galactic molecular clouds; A11 consider clouds with initial radius of 3.5 pc; BR08 consider very compact clouds with a radius of 1 pc. In the case of HN09 the involved clouds are small (with a radius of pc), but this cannot be compared directly with the other papers, as HN09 simulate a collision between the two clouds, which changes deeply the properties of the involved gas.
M08, HN09 and this paper include a cuspy rigid potential to account for the other stars close to the GC (see, e.g., G03), whereas BR08 and A11 do not include any underlying rigid potentials. The inclusion of a spherical potential contributes to stabilize the disk (delaying or preventing local fragmentation), although the mass of the SMBH dominates any other gravitational contributions in the region of interest ( pc).
The simulations in M08 are isothermal, with T K, respectively. Such temperature is likely too low, if compared to the background radiation field in the GC (see e.g. Nagai et al. 2007). In BR08 the simulations include an approximate radiative transfer formalism, with compressional heating balanced by cooling rates derived from estimated optical depths. In HN09, the simulations include a very simplified model of cooling (see their equation 1), and the initial temperature of the cloud is very low (T K). HN09 do not specify whether there is a temperature floor in their simulations and do not provide a quantitative description of temperature evolution (see their appendix for a discussion about cooling). In A11 and in this paper, the simulations include different thermodynamical treatments for the gas, including both isothermal and radiative cooling cases (even an adiabatic run in A11). The floor temperature for the simulations with radiative cooling is set to be 50 K in A11 and 100 K in both BR08 and this paper. For a further analysis of the differences between our radiative transfer formalism and that used in BR08, see Appendix B.
A11 stop their simulation before fragmentation takes place in the disk, whereas the other considered papers study the formation of star candidates in the disk. M08, BR08 and HN09 adopt the sink particle technique, to model SF. Only in this paper we follow the fragmentation of the disk, although our resolution is border-line for the lower star masses ( M). The attempt of directly following the fragmentation, without using the sink particles, is the main improvement of this paper with respect to the previous ones.
3.4.2 Comparison of the results
All the considered papers indicate that a dense gaseous disk may form around the GC as a consequence of the infall and of the disruption of a massive molecular cloud. If the impact parameter is non-zero, the resulting gaseous disk will be eccentric, consistently with the observations of the stellar orbits in Sgr A. The eccentricity and the size of the final gaseous disk are similar in all the considered papers.
The main differences among various simulations refer to the formation of star candidates. M08 adopted a value of the critical density for converting gaseous particles into sink particles ( cm) that is too low if compared to the tidal density induced by the presence of the SMBH inside the disk. Such assumption generally leads to the formation of sink particles with larger masses than those we expect for cold ( K) clouds. HN09 assume a more robust critical density, weighted with the local tidal density ( cm, where is assumed and is the tidal density). This choice of the critical density for converting gaseous particles into sink particles is very close to our criterion for SF (see Section 3.2). On the other hand, the MF in HN09 is quite low-mass (especially in their runs S3 and S4), because of the approximations in the cooling recipe and because of the absence of opacity prescriptions.
BR08 adopt a very conservative value of the critical density for converting gaseous particles into sink particles ( M pc cm, assuming ), well above the critical tidal density. Therefore, their MF is consistent with that predicted by the Jeans mass for the local density and temperature of the clouds.
Similarly, the MFs derived in this paper (Fig. 8) are consistent with the predictions calculated from Jeans mass and Toomre instability. Furthermore, the fact that we do not use sink particles prevents any possible bias connected with that method. The selection of star candidates on the basis of (see Section 3.2) is shown to be robust by the analysis of Fig. 11. Confirming the results by BR08, we show that the MF in the stellar ring may be top-heavy only if the parent cloud is sufficiently massive ( M) and if sufficiently high temperatures ( K) are reached during the SF.
On the other hand, the differences between our run B (massive isothermal cloud with T K) and our run E (massive cloud with radiative cooling and T K) are quite small. Furthermore, we do not observe in our run E the formation of the very massive stars ( M) that were found in the massive cloud simulated by BR08. The MF derived from run E is consistent with a single power-law with index , whereas the MF by BR08 is clearly bimodal, showing two distinct stellar populations (see fig. 4 by BR08). The very massive stars in BR08 are all formed at pc, where high temperatures and Jeans masses are reached (as it can be seen from their figs. S1, S2 and S3). In our simulations, star candidates never form at pc, because the shear from the SMBH prevents local collapse. Instead, in the simulation by BR08, the formation of sink particles is very efficient at pc (see their fig. 2). We point out that the mass of the cloud, the mass resolution of gas particles, the softening and smoothing lengths and the temperature floor in the two simulations are very similar. The reasons of the difference are likely the following. First, the orbit of the two clouds are different, leading to different distributions of position and angular momentum around the SMBH. We start integrating when the cloud is 25 pc far from the GC, whereas BR08 assume a much smaller initial distance (3 pc). Furthermore, we adopt a cloud radius pc, whereas BR08 assume that the cloud radius is only pc (for approximately the same cloud mass). This implies that the density of gas streams in the disrupting cloud of BR08 is significantly higher than in our simulations. This is evident also from Figs. 11 and 12, where only few gas clumps reach the density cm, that is the density threshold for the conversion of gas into sinks in BR08. Furthermore, our implementation of the opacity and that used by BR08 have some differences. In particular, Stamatellos et al. (2007, from which BR08 take their formalism) have a deep opacity gap at temperature T K, when the refractory elements start to vaporize (see Appendix B). In addition, Bell & Lin (1994) Rosseland opacity (adopted by BR08) is too low at T K (see, e.g., Alexander & Ferguson 1994), that is the temperature at which most of massive stars form in BR08. This might justify a heavier MF in BR08 than in our paper. In addition, we do not use the sink particle method (apart from the SMBH), whereas BR08 do. In particular, BR08 adopt a very robust minimum density ( cm) to create new sink particles. Instead, we assume that star particles form where the local density of gas is cm. Our density threshold is three orders of magnitude smaller than the threshold by BR08, but is well motivated, because clumps of particles with density above are gravitationally bound and have a thermal behavior completely distinct from the rest of the disk (see our Fig 11 and Section 3.3). It may be that the higher density threshold required by BR08, together with the sink recipes, introduces a selection bias in the formation of stars and suppresses the formation of low-mass stars (see, e.g. the discussion in Appendix A). We stress that all very massive stars in BR08 form very close ( pc) to the GC, where massive stars have not been observed in the Milky Way (the observed ring of young stars having an inner radius of pc). From this point of view, our run E seems to match better the observational properties of the Milky Way.
4. Summary and conclusions
We simulate the infall of a molecular cloud toward SgrA. We ran five different clouds, with various masses (, M), initial temperatures (100, 500 K) and thermodynamics (isothermal or radiative cooling). In all the simulations, the cloud is disrupted by the tidal forces of the SMBH and spirals towards it. In less than yr, more than one tenth of the gas in the parent cloud ends up into a dense, eccentric and distorted disk around the SMBH, with a small outer radius ( pc). Locally, the surface density of the gaseous disk may overcome the tidal shear from the SMBH and fragmentation may take place. We define as ‘star candidates’ those gas clumps that decouple from the background gas and collapse. This is the first attempt to trace the fragmentation of the disk, without adopting the sink particle technique.
Star candidates form in four of the five simulations within the first yr. Only run C, in which the cloud is small ( M) and warm (500 K), does not host star candidates at yr, but stars might form at later times.
In the remaining four simulations, star candidates are distributed in a thin ring at a distance of pc from the SMBH. They have eccentric orbits (), with average eccentricity , in agreement with the properties of the young massive stars observed around Sgr A (, Bartko et al. 2009).
In the three runs (A, B and E) where the cloud is more massive ( M), the total mass of star candidates ( M) is consistent with the observations (P06). Instead, the total mass of stars is only M in run D, where the cloud is smaller ( M). Therefore, our simulations suggest that a rather massive molecular cloud ( M) is required, in order to form the stars observed around the GC.
Furthermore, the MF of stellar candidates is very top-heavy (fitted by a single power-law with ) only in run A, for which T K. In runs B and E (for which T K) the MF (fitted by a single power-law with ) is still heavier than the Salpeter MF. Run E, which includes radiative cooling, differs from run B (isothermal) only for the presence of a high mass tail in the MF.
This is consistent with the predictions from Jeans instability. Therefore, high gas temperatures (T K) should be reached during the SF, in order to generate a top-heavy MF, such as the one observed in the GC. Such high temperatures are not unrealistic in the neighborhood of Sgr A. We stress that, in addition to the total mass of the cloud, the density of the original cloud can strongly affect the stellar MF. This is likely the main difference between our simulations and those by BR08, together with the usage of sink particles.
A possible limit of our simulations is the mass resolution: we barely resolve M stars. Higher resolution simulations are important, to improve our understanding of SF around SMBHs. Furthermore, in our five runs the cloud has always the same orbital parameters. It would be interesting to study the influence of the cloud orbit onto the formation of the gaseous disk and of the stellar ring. This issue will be also addressed in a forthcoming paper. In this paper, we consider only a single molecular cloud. HN09 find that the collision between two gas clouds is naturally able to explain the existence of more than one stellar disk. Therefore, the scenario of the cloud collision deserves further investigations, with accurate cooling and opacity prescriptions.
A number of other open questions deserve further studies. For example, this paper focuses on the formation of the stars in the disk, but the long-term dynamical evolution of these stars cannot be studied with the same numerical approach (and we remind that the observed massive stars in the GC have an estimate age of Myr). Various papers addressed this issue before, with somewhat different results (see, e.g., Alexander et al. 2007; Cuadra et al. 2008; Löckman et al. 2010).
A possibly related issue is the origin of the so-called ‘S-stars’ (e.g., Schödel et al. 2002; Ghez et al. 2003; Bartko et al. 2010). These are B stars that do not belong to the disk(s), as have isotropic spatial distribution and higher eccentricities. They are predominantly located in the inner 0.05 pc, but recent studies (Bartko et al. 2010) show that their distribution extends further out. Our simulations cannot explain at the same time young disk stars and isotropic S-stars, because they are not consistent with the same stellar population. On the other hand, it is likely that the disruption of two different clouds at different times can account for both populations (e.g., HN09; Perets & Gualandris 2010). In this scenario, the S-stars are the survivors of a previous disk and their orbits have been randomized by the influence of a perturber (e.g., Perets et al. 2007; Yu et al. 2007; Perets et al. 2008; Gualandris & Merritt 2009; Perets & Gualandris 2010), or by a Kozai-like resonance between two disks (Löckmann et al. 2008, 2009), or by spiral density wave migration in the gaseous disk (Griv 2010). We propose that even the evaporation of the parent gaseous disk (due, e.g., to the local UV and X-ray emission) might randomize the orbits of the stars, if the SF efficiency was sufficiently low. On the other hand, we cannot exclude on the basis of our simulations a completely different origin for the young disk(s) and for the S-stars (e.g., Alexander & Livio 2004; Fujii et al. 2010; Perets & Gualandris 2010; Baruteau, Cuadra & Lin 2011; Madigan et al. 2011).
We thank the anonymous referee for providing constructive comments that improved the paper, and D. Semenov for enlightening discussions about opacity. We also thank A. Gualandris, B. Moore, E. D’Onghia, E. Ripamonti and P. Englmaier for useful discussions, and acknowledge J. Stadel and D. Potter for technical support. The simulations were run onto the CRAY XT3 cluster at the Swiss National Supercomputing Center in Lugano, onto the cluster at the Institute for Theoretical Physics of the University of Zurich and onto the ‘lagrange’ cluster at the Consorzio Interuniversitario Lombardo per L’Elaborazione Automatica (CILEA). MM, TH and LM acknowledge support from the Swiss National Science Foundation.
- () Agertz, O., et al. 2007, MNRAS, 380, 963
- () Alexander, R. D., Armitage, P. J., Cuadra, J., Begelman, M. C. 2008, ApJ, 674, 927
- () Alexander, D. R., Ferguson, J. W. 1994, ApJ, 437, 879
- () Alexander, T., Livio, M. 2004, ApJ, 606L, 21
- () Alexander, R. D., Begelman, M. C., Armitage, Ph. J., 2007, ApJ, 654, 907
- () Alexander, R. D., Smedley, S. L., Nayakshin, S., King, A. R., 2011, MNRAS, accepted, arXiv:1109.4148
- () Alig, C., Burkert, A., Johansson, P. H., Schartmann, M. 2011, MNRAS, 412, 469 (A11)
- () Bartko, H., et al. 2009, ApJ, 697, 1741
- () Bartko, H., et al. 2010, ApJ, 708, 834
- () Baruteau, C., Cuadra, J., Lin, D. N. C. 2011, ApJ, 726, 28
- () Bell, K. R., Lin, D. N. C. 1994, ApJ, 427, 987
- () Binney, J., Tremaine, S. 1987, Galactic dynamics, Princeton University Press
- () Boley, A. C., Hayfield, T., Mayer, L., Durisen, R. H. 2010, Icarus, 207, 509
- () Boley, A. C. 2009, ApJ Letters, 695, 53
- () Bonnell, I. A., Rice, W. K. M. 2008, Science, 321, 1060 (BR08)
- () Burkert, A., Schartmann, M., Alig, Ch., Gillessen, S., Genzel, R., Fritz, T., Eisenhauer, F., 2012, ApJ, submitted, arXiv:1201.1414
- () Collin, S., Zahn, J.-P. 2008, A&A, 477, 419
- () Christopher, M. H., Scoville, N. Z., Stolovy, S. R., Yun, M. S. 2005, ApJ, 622, 346
- () Cuadra, J., Armitage, P. J., Alexander, R. D. 2008, MNRAS, 388L, 64
- () D’Alessio, P., Calvet, N., Hartmann, L. 2001, ApJ, 553, 321
- () Eisenhauer, F., et al. 2005, ApJ, 628, 246
- () Fujii, M., Iwasawa, M., Funato, Y., Makino, J. 2008, ApJ, 686, 1082
- () Fujii, M., Iwasawa, M., Funato, Y., Makino, J. 2009, ApJ, 695, 1421
- () Fujii, M., Iwasawa, M., Funato, Y., Makino, J. 2010, ApJ, 716L, 80
- () Genzel, R., Crawford, M. K., Townes, C. H., Watson, D. M., 1985, ApJ, 297, 766
- () Genzel, R., et al. 2003, ApJ, 594, 812 (G03)
- () Genzel, R., Eisenhauer, F., Gillessen, S. 2010, Reviews of Modern Physics, 82, 3121
- () Gerhard, O. 2001, ApJ, 546, L39
- () Ghez, A. M., et al. 2003, ApJ, 586L, 127
- () Ghez, A. M., Salim, S., Hornstein, S. D., Tanner, A., Morris, M., Becklin, E. E., Duchene, G. 2005, ApJ, 620, 744
- () Gillessen, S., et al. 2012, Nature, in press, arXiv:1112.3264
- () Goodman, J. 2003, MNRAS, 339, 937
- () Griv, E. 2010, ApJ, 709, 597
- () Gualandris, A., Merritt, D. 2009, ApJ, 705, 361
- () Gürkan, M. A., Rasio, F. A. 2005, ApJ, 628, 236
- () Güsten, R., Genzel, R., Wright, M. C. H., Jaffe, D. T. 1987, ApJ, 318, 124
- () Ho, P. T. P., Ho, L. C., Szczepanski, J. C., Jackson, J. M., Armstrong, J. T., Barrett, A. H. 1991, Nature, 350, 309
- () Hobbs, A., Nayakshin, S. 2009, MNRAS, 394, 191 (HN09)
- () Hoyle, F., Lyttleton, R. A. 1939, Proc. Camb. Phil. Soc.,, 34, 405
- () Jeans, J. H. 1919, Problems of cosmogony and stellar dynamics, Cambridge, University press
- () Kim, S. S. Figer, D. F. Morris, M. 2004, ApJ, 607, L123
- () Kim, S. S., Morris, M 2003, ApJ, 597, 312
- () Kolykhalov, P. I., Syunyaev, R. A 1980, SOVIET ASTRON. LETT (TR. PISMA ASTR. ZH.), 6, 357
- () Larson, R. B. 1981, MNRAS, 194, 809
- () Levin, Y., Beloborodov, A. M. 2003, ApJ, 590L, 33
- () Löckmann, U., Baumgardt, H., Kroupa, P. 2008, ApJ, 683, L151
- () Löckmann, U., Baumgardt, H., Kroupa, P. 2009, MNRAS, 398, 429
- () Löckmann, U., Baumgardt, H., Kroupa, P. 2010, MNRAS, 402, 519
- () Lu, J. R., Ghez, A. M., Hornstein, S. D., Morris, M. R., Becklin, E. E., Matthews, K. 2009, ApJ, 690, 1463
- () Madigan, A.-M., Hopman, C., Levin, Y. 2011, ApJ, 738, 99
- () Mapelli, M., Hayfield, T., Mayer, L., Wadsley, J. 2008, arXiv0805.0185 (M08)
- () McMillan, S. L. W., Portegies Zwart, S. F. 2003, ApJ, 596, 314
- () Milosavljevic, M., Loeb, A. 2004, ApJ, 604L, 45
- () Nagai, M., Tanaka, K., Kamegai, K., Oka, T. 2007, PASJ, 59, 25
- () Nayakshin, S., Cuadra, J. J. 2005, A&A, 437, 437
- () Nayakshin, S., Cuadra, J., Springel, V. 2007, MNRAS, 379, 21 (NCS07)
- () Nayakshin, S., Sunyaev, R. 2005, MNRAS, 364L, 23
- () Novak, G., Dotson, J. L., Dowell, C. D., Hildebrand, R. H., Renbarger, T., Schleuning, D. A. 2000, ApJ, 529, 241
- () Oka, T., Nagai, M., Kamegai, K., Tanaka, K., Kuboi, N. 2007, PASJ, 59, 15
- () Okumura, S. K., Ishiguro, M., Fomalont, E. B., Hasegawa, T., Kasuga, T., Morita, K. I., Kawabe, R., Kobayashi, H. 1991, ApJ, 378, 127
- () Paumard, T., et al. 2006, ApJ, 643, 1011 (P06)
- () Perets, H. B., Hopman, C., Alexander, T. 2007, ApJ, 656, 709
- () Perets, H. B., Gualandris, A., Merritt, D., Alexander, T. 2008, Memorie della Società Astronomica Italiana, 79, 1100
- () Perets, H. B., Gualandris, A. 2010, ApJ, 719, 220
- () Poliachenko, V. L., Shukhman, I. G. 1977, Soviet Astronomy Letters, 3, 134
- () Portegies Zwart, S. F., Baumgardt, H., McMillan, S. L. W., Makino, J., Hut, P., Ebisuzaki, T. 2006, ApJ, 641, 319
- () Portegies Zwart, S. F., McMillan, S. L. W., Gerhard, O. 2003, ApJ, 593, 352
- () Read, J. I., Hayfield, T., Agertz, O. 2010, MNRAS, 405, 1513
- () Rice, W. K. M., Lodato, G., Armitage, P. J. 2005, MNRAS, 364, L56
- () Sanders, R. H. 1998, MNRAS, 294, 35
- () Schödel, R., et al. 2002, Nature, 419, 694
- () Scoville, N. Z., Stolovy, S. R., Rieke, M., Christopher, M., Yusef-Zadeh, F. 2003, ApJ, 594, 294
- () Shlosman, I., Begelman, M. C. 1987, Nature, 329, 810
- () Solomon, P. M., Scoville, N. Z., Jefferts, K. B., Penzias, A. A., Wilson, R. W. 1972, ApJ, 178, 125
- () Stamatellos, D., Whitworth, A. P., Bisbas, T., Goodwin, S. 2007, A&A, 475, 37
- () Toomre, A. 1964, ApJ, 139, 1217
- () Wadsley, J. W., Stadel, J., Quinn, T. 2004, New Astronomy, 9, 137
- () Wardle, M., Yusef-Zadeh, F. 2008, ApJ, 683L, 37
- () Wardle, M., Yusef-Zadeh, F. 2011, ApJ Letters, submitted, arXiv:1108.2175
- () Yu, Q., Lu, Y., Lin, D. N. C. 2007, ApJ, 666, 919
- () Yusef-Zadeh, F., Braatz, J., Wardle, M., Roberts, D. 2008, ApJ, 683, L147
- () Yusef-Zadeh, F., Lacy, J. H., Wardle, M., Whitney, B., Bushouse, H., Roberts, D. A., Arendt, R. G., 2010, ApJ, 725, 1429
- () Yusef-Zadeh F., Wardle M., 2011, to appear in ”The Central Kiloparsec in Galactic Nuclei: Astronomy at High Angular Resolution 2011”, open access Journal of Physics: Conference Series (JPCS), published by IOP Publishing, arXiv:1112.3279
- () Yusef-Zadeh, F., Bushouse, H., Wardle, M. 2012, ApJ, 744, 24
- () Zhao, J.-H., Morris, M. R., Goss, W. M., An, T. 2009, ApJ, 699, 186
- () Zhao, J.-H., Blundell, R., Moran, J. M., Downes, D., Schuster, K. F., Marrone, D. P. 2010, ApJ, 723, 1097
Appendix A How critical is the choice of the threshold density?
We adopt a threshold density cm (assuming molecular weight ) and a threshold radius pc, in order to identify star candidates. Such choices are reasonable, as the adopted is of the same order of magnitude as the softening length (therefore, we cannot consider smaller radii) and the adopted is at least times larger than the local tidal density (therefore, the gravitational torque by the SMBH cannot stop the collapse of these clumps).
Slightly different values of and of do not change significantly our results. Fig. 13 shows that the MF of star candidates in run A is affected only for stellar masses , if we change from cm to cm. In particular, lower mass stars are suppressed by larger choices of . In the case of run E, where the gas is initially at a lower temperature (T K), the MF is substantially unchanged, even when changing by one order of magnitude.
Interestingly, BR08 adopt a significantly higher threshold value ( cm), in order to convert gas clumps into sink particles. If we use the same density as in our simulations, we obtain the results shown in Fig. 14. A much lower number of star candidates is derived, by assuming this threshold, and most of the low-mass star candidates are suppressed.
Fig. 15 shows how different choices of the threshold radius can affect the MF. Smaller values of (down to about one tenth of the fiducial value) do not change significantly the MF. Similarly, the MF is not significantly affected by larger values of , up to times the fiducial value. The choice of larger values of is more critical: the MF becomes significantly heavier in all the considered runs. However, there are no physical reasons to adopt such large values of .
In conclusion, the choice of and of is crucial for the final MF. We adopt physically well motivated values for both and . Slightly different choices of and do not affect significantly the MF. Significantly different choices of and (e.g., pc and cm) change both the normalization and the slope of the MF. However, in the current study, we are mainly interested in the relative differences between various runs, that are almost not affected by the choice of and of .
Appendix B Different treatments of opacity
In run E, we use the same radiative cooling algorithm as that described in Boley (2009) and in Boley et al. (2010). D’Alessio et al. (2001) Rosseland and Planck opacities are used in our code.
BR08 use the formalism described in Stamatellos et al. (2007). They compute Rosseland opacities using Bell & Lin (1994) recipes.
The primary differences are: (i) Stamatellos et al. (2007) fold the gas potential (but not the star potential) and density into an (effective) pseudo density, and we do not. (ii) Stamatellos et al. (2007) gas energy equation includes a whole range of source terms related to dissociation and ionization of H, H and He. Our gas energy equation does not. (iii) Bell & Lin (1994) opacities are (orders of magnitude) too low at T K, because they truncate the water opacities at too short a wavelength (D. Semenov, private communication; see also Alexander & Ferguson 1994).
The fact that we do not include source terms related to dissociation of H is an issue for our model, but only for temperatures above K, that are basically never reached in our run E (see Fig. 11). Instead, the fact that Bell & Lin (1994) opacities (adopted by BR08) are too low at T K may have serious consequences on the MF, as most of massive stars in BR08 form in that range of temperatures (see their figure S1). Too low opacities might artificially boost fragmentation, inducing spurious SF at temperatures (T K) where the Jeans mass is high.
In Fig. 16, we plot our opacity versus temperature, along isopycnals from , to g cm (bottom to top). This does not totally match figure 5 of Stamatellos et. al (2007), that have a very deep opacity gap at T K, when the refractory elements start to vaporize. This might partially explain why BR08 MF is top-heavier than ours.